src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Class
  25     tractionBoundaryGradient
  26
  27 Author
  28     Philip Cardiff UCD
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "tractionBoundaryGradient.H"
  33 #include "fvPatch.H"
  34 #include "Switch.H"
  35 //#include "plasticityModel.H"
  36 #include "constitutiveModel.H"
  37 #include "thermalModel.H"
  38
  39
  40 // * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * * //
  41
  42 Foam::tmp<Foam::vectorField> Foam::tractionBoundaryGradient::traction
  43 (
  44     const tensorField& gradField,
  45     const word& workingFieldName,
  46     const word& integralFieldName,
  47     const fvPatch& patch,
  48     const bool orthotropic,
  49     const nonLinearGeometry::nonLinearType& nonLinear,
  50     const bool incremental
  51 )
  52 {
  53     // Create result
  54     tmp<vectorField> ttraction
  55     (
  56         new vectorField(gradField.size(), vector::zero)
  57     );
  58     vectorField& traction = ttraction();
  59
  60     // Orthotropic material
  61     if (orthotropic)
  62     {
  63         // For now, turn off orthotropic
  64         FatalErrorIn
  65         (
  66             "tmp<vectorField> tractionBoundaryGradient::traction\n"
  67             "(\n"
  68             "    const tensorField& gradField,\n"
  69             "    const word& workingFieldName,\n"
  70             "    const word& integralFieldName,\n"
  71             "    const fvPatch& patch,\n"
  72             "    const bool orthotropic,\n"
  73             "    const nonLinearGeometry::nonLinearType& nonLinear\n"
  74             ") const"
  75         )   << "tractionBoundaryGradient::traction is not written for"
  76             << " orthotropic materials yet" << nl
  77             << "you are probably trying to use solidContact boundaries "
  78             << "with the orthotropic solver" << nl
  79             << exit(FatalError);
  80     }
  81     else
  82     {
  83         // Lookup material properties from the solver
  84         const fvPatchScalarField& mu =
  85             patch.lookupPatchField<volScalarField, scalar>("mu");
  86
  87         const fvPatchScalarField& lambda =
  88             patch.lookupPatchField<volScalarField, scalar>("lambda");
  89
  90         // only for nonlinear elastic properties
  91         // if (rheology.type() == plasticityModel::typeName)
  92         // {
  93         //     const plasticityModel& plasticity =
  94         //       refCast<const plasticityModel>(rheology);
  95
  96         //     mu = plasticity.newMu().boundaryField()[patch.index()];
  97         //     lambda = plasticity.newLambda().boundaryField()[patch.index()];
  98         // }
  99
 100         // Get patch normal
 101         vectorField n = patch.nf();
 102
 103         // Calculate traction
 104         traction = 2*mu*(n & symm(gradField)) + lambda*tr(gradField)*n;
 105
 106         // Plasticity effects
 107         const constitutiveModel& rheology =
 108             patch.boundaryMesh().mesh().objectRegistry::
 109             lookupObject<constitutiveModel>("rheologyProperties");
 110
 111         if (rheology.plasticityActive())
 112         {
 113             traction -=
 114                 2*mu*(n & rheology.DEpsilonP().boundaryField()[patch.index()]);
 115         }
 116
 117         // Thermal effects
 118         if
 119         (
 120             patch.boundaryMesh().mesh().objectRegistry::
 121             foundObject<thermalModel>("thermalProperties")
 122         )
 123         {
 124             const thermalModel& thermo =
 125                 patch.boundaryMesh().mesh().objectRegistry::
 126                 lookupObject<thermalModel>("thermalProperties");
 127
 128             const fvPatchScalarField& threeKalpha =
 129                 patch.lookupPatchField<volScalarField, scalar>
 130                 ("((threeK*rho)*alpha)");
 131
 132             // Incremental thermal contribution
 133             if (incremental)
 134             {
 135                 const fvPatchScalarField& DT =
 136                     patch.lookupPatchField<volScalarField, scalar>("DT");
 137
 138                 traction -= (n*threeKalpha*DT);
 139             }
 140             else
 141             {
 142                 const fvPatchScalarField& T =
 143                     patch.lookupPatchField<volScalarField, scalar>("T");
 144
 145                 const scalarField T0 =
 146                     thermo.T0()().boundaryField()[patch.index()];
 147
 148                 traction -= (n*threeKalpha*(T - T0));
 149             }
 150         }
 151
 152         // Non-linear terms
 153         if
 154         (
 155             nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
 156          || nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
 157         )
 158         {
 159             traction +=
 160                 (n & (mu*(gradField & gradField.T())))
 161               + 0.5*n*lambda*(gradField && gradField);
 162
 163             if
 164             (
 165                 incremental
 166              && nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
 167             )
 168             {
 169                 // Incremental total Lagrangian
 170
 171                 const fvPatchTensorField& gradU =
 172                     patch.lookupPatchField<volTensorField, tensor>
 173                     (
 174                         "grad(" + integralFieldName + ")"
 175                     );
 176
 177                 traction +=
 178                     (
 179                         n
 180                       & (
 181                           mu*
 182                           (
 183                               (gradU & gradField.T())
 184                             + (gradField & gradU.T())
 185                             )
 186                         )
 187                     )
 188                   + 0.5*n*lambda*tr
 189                     (
 190                         (gradU & gradField.T())
 191                       + (gradField & gradU.T())
 192                     );
 193             }
 194         }
 195
 196
 197         // Add old stress for incremental approaches
 198         if (incremental)
 199         {
 200             // add on old traction
 201             const fvPatchSymmTensorField& sigma =
 202                 patch.lookupPatchField<volSymmTensorField, symmTensor>
 203                 (
 204                     "sigma"
 205                 );
 206
 207             traction += (n & sigma);
 208         }
 209
 210         // Visco-elastic effects
 211         if (rheology.viscoActive())
 212         {
 213             const fvPatchSymmTensorField& DSigmaCorr =
 214                 patch.lookupPatchField<volSymmTensorField, symmTensor>
 215                 (
 216                     "DSigmaCorr"
 217                 );
 218
 219             traction += (n & DSigmaCorr);
 220         }
 221
 222         // Updated Lagrangian or total Lagrangian large strain
 223         if
 224         (
 225             nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
 226          || nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
 227         )
 228         {
 229             tensorField F = I + gradField;
 230
 231             if
 232             (
 233                 incremental
 234              && nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
 235             )
 236             {
 237                 // Incremental total Lagrangian
 238
 239                 const fvPatchTensorField& gradU =
 240                     patch.lookupPatchField<volTensorField, tensor>
 241                     (
 242                         "grad(" + integralFieldName + ")"
 243                     );
 244
 245                 F += gradU;
 246             }
 247
 248             tensorField Finv = inv(F);
 249             scalarField J = det(F);
 250
 251             // Rotate second Piola Kirchhoff traction to be Cauchy traction
 252             traction = (traction & F)/(mag(J * Finv & n));
 253         }
 254     }
 255
 256     return ttraction;
 257 }
 258
 259
 260 // * * * * * * * * * * * * * * * * Operators  * * * * * * * * * * * * * * //
 261
 262 Foam::tmp<Foam::vectorField> Foam::tractionBoundaryGradient::snGrad
 263 (
 264     const vectorField& traction,
 265     const scalarField& pressure,
 266     const word& workingFieldName,
 267     const word& integralFieldName,
 268     const fvPatch& patch,
 269     const bool orthotropic,
 270     const nonLinearGeometry::nonLinearType& nonLinear,
 271     const bool incremental
 272 )
 273 {
 274     // Create result
 275     tmp<vectorField> tgradient(new vectorField(traction.size(), vector::zero));
 276     vectorField& gradient = tgradient();
 277
 278     // Orthotropic material
 279     if (orthotropic)
 280     {
 281         // Get mechanical properties
 282         const constitutiveModel& rheology =
 283             patch.boundaryMesh().mesh().objectRegistry::
 284             lookupObject<constitutiveModel>("rheologyProperties");
 285
 286         const diagTensorField K =
 287             rheology.K()().boundaryField()[patch.index()];
 288
 289         const symmTensor4thOrderField C =
 290             rheology.C()().boundaryField()[patch.index()];
 291
 292         // Required fields
 293         vectorField n = patch.nf();
 294         const diagTensorField Kinv = inv(K);
 295         const fvPatchTensorField& gradField =
 296             patch.lookupPatchField<volTensorField, tensor>
 297             (
 298                 "grad(" + workingFieldName + ")"
 299             );
 300
 301         // Calculate the traction to apply
 302         vectorField Traction(n.size(), vector::zero);
 303         tensorField sigmaExp(n.size(), tensor::zero);
 304
 305         // Total Lagrangian, small strain
 306         if
 307         (
 308             !incremental
 309          && nonLinear == nonLinearGeometry::OFF
 310         )
 311         {
 312             // Use total traction
 313             Traction = (traction - n*pressure);
 314
 315             sigmaExp = (n*(n&(C && symm(gradField)))) - (K & gradField);
 316         }
 317         // Incremental total Lagrangian small strain
 318         else if
 319         (
 320             incremental
 321          && nonLinear == nonLinearGeometry::OFF
 322         )
 323         {
 324             const fvPatchSymmTensorField& sigma =
 325                 patch.lookupPatchField<volSymmTensorField, symmTensor>
 326                 (
 327                     "sigma"
 328                 );
 329
 330             // Increment of traction
 331             Traction = (traction - n*pressure) - (n & sigma);
 332
 333             sigmaExp = (n*(n&(C && symm(gradField)))) - (K & gradField);
 334         }
 335         // Updated Lagrangian large strain
 336         else if
 337         (
 338             nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
 339         )
 340         {
 341             const fvPatchSymmTensorField& sigma =
 342                 patch.lookupPatchField<volSymmTensorField, symmTensor>
 343                 (
 344                     "sigma"
 345                 );
 346
 347             tensorField F = I + gradField;
 348             tensorField Finv = inv(F);
 349             scalarField J = det(F);
 350             vectorField nCurrent = Finv & n;
 351             nCurrent /= mag(nCurrent) + SMALL;
 352             vectorField tractionCauchy = traction - nCurrent*pressure;
 353
 354             // Increment of 2nd Piola-Kirchhoff traction
 355             Traction = mag(J*(Finv & n))*(tractionCauchy & Finv) - (n & sigma);
 356
 357             sigmaExp = n*(n &(C && symm(gradField))) - (K & gradField);
 358         }
 359         else
 360         {
 361             FatalErrorIn
 362             (
 363                 "tmp<vectorField> tractionBoundaryGradient::snGrad\n"
 364                 "(\n"
 365                 "    const vectorField& traction,\n"
 366                 "    const scalarField& pressure,\n"
 367                 "    const word& workingFieldName,\n"
 368                 "    const word& integralFieldName,\n"
 369                 "    const fvPatch& patch,\n"
 370                 "    const bool orthotropic,\n"
 371                 "    const nonLinearGeometry::nonLinearType& nonLinear,\n"
 372                 "    const bool incremental\n"
 373                 ") const"
 374             )   << "solidTractionOrtho boundary condition not suitable for "
 375                 << " field " << workingFieldName << " and "
 376                 << nonLinearGeometry::nonLinearNames_[nonLinear]
 377                 << abort(FatalError);
 378         }
 379
 380         gradient = n & (Kinv & ( n*(Traction) - sigmaExp ));
 381     }
 382     else
 383     {
 384         // Standard isotropic solvers
 385
 386         // Lookup material properties from the solver
 387         const fvPatchScalarField& mu =
 388             patch.lookupPatchField<volScalarField, scalar>("mu");
 389
 390         const fvPatchScalarField& lambda =
 391             patch.lookupPatchField<volScalarField, scalar>("lambda");
 392
 393         // only for nonlinear elastic properties
 394         // if (rheology.type() == plasticityModel::typeName)
 395         // {
 396         //     const plasticityModel& plasticity =
 397         //       refCast<const plasticityModel>(rheology);
 398
 399         //     mu = plasticity.newMu().boundaryField()[patch.index()];
 400         //     lambda = plasticity.newLambda().boundaryField()[patch.index()];
 401         // }
 402
 403         vectorField n = patch.nf();
 404
 405         // gradient of the field
 406         const fvPatchTensorField& gradField =
 407             patch.lookupPatchField<volTensorField, tensor>
 408             (
 409                 "grad(" + workingFieldName + ")"
 410             );
 411
 412
 413         vectorField Traction(n.size(), vector::zero);
 414
 415         // Total Lagrangian, small strain
 416         if
 417         (
 418             !incremental
 419          && nonLinear == nonLinearGeometry::OFF
 420         )
 421         {
 422             // Use total traction
 423             Traction = (traction - n*pressure);
 424         }
 425         // Incremental total Lagrangian small strain
 426         else if
 427         (
 428             incremental
 429          && nonLinear == nonLinearGeometry::OFF
 430         )
 431         {
 432             const fvPatchSymmTensorField& sigma =
 433                 patch.lookupPatchField<volSymmTensorField, symmTensor>
 434                 (
 435                     "sigma"
 436                 );
 437
 438             // Increment of traction
 439             Traction = (traction - n*pressure) - (n & sigma);
 440         }
 441         else if
 442         (
 443             nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
 444          || nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
 445         )
 446         {
 447             const fvPatchSymmTensorField& sigma =
 448                 patch.lookupPatchField<volSymmTensorField, symmTensor>
 449                 (
 450                     "sigma"
 451                 );
 452
 453             tensorField F = I + gradField;
 454
 455             if
 456             (
 457                 incremental
 458              && nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
 459             )
 460             {
 461                 // Incremental total Lagrangian
 462
 463                 const fvPatchTensorField& gradU =
 464                     patch.lookupPatchField<volTensorField, tensor>
 465                     (
 466                         "grad(" + integralFieldName + ")"
 467                     );
 468
 469                 F += gradU;
 470             }
 471
 472             tensorField Finv = hinv(F);
 473             scalarField J = det(F);
 474             vectorField nCurrent = Finv & n;
 475             nCurrent /= mag(nCurrent) + SMALL;
 476             vectorField tractionCauchy = traction - nCurrent*pressure;
 477
 478             Traction = mag(J*(Finv & n))*(tractionCauchy & Finv);
 479
 480             if (incremental)
 481             {
 482                 // Increment of 2nd Piola-Kirchhoff traction
 483                 Traction -= (n & sigma);
 484             }
 485         }
 486         else
 487         {
 488             FatalErrorIn
 489             (
 490                 "tmp<vectorField> tractionBoundaryGradient::snGrad\n"
 491                 "(\n"
 492                 "    const vectorField& traction,\n"
 493                 "    const scalarField& pressure,\n"
 494                 "    const word& workingFieldName,\n"
 495                 "    const word& integralFieldName,\n"
 496                 "    const fvPatch& patch,\n"
 497                 "    const bool orthotropic,\n"
 498                 "    const nonLinearGeometry::nonLinearType& nonLinear,\n"
 499                 "    const bool incremental\n"
 500                 ") const"
 501             )   << " field " << workingFieldName << " and "
 502                 << nonLinearGeometry::nonLinearNames_[nonLinear]
 503                 << " nonLinear are not compatible!"
 504                 << abort(FatalError);
 505       }
 506
 507         //- visco-elastic
 508         const constitutiveModel& rheology =
 509             patch.boundaryMesh().mesh().objectRegistry::
 510             lookupObject<constitutiveModel>("rheologyProperties");
 511
 512         if (rheology.viscoActive())
 513         {
 514             const fvPatchSymmTensorField& DSigmaCorr =
 515                 patch.lookupPatchField<volSymmTensorField, symmTensor>
 516                 (
 517                     "DSigmaCorr"
 518                 );
 519
 520             Traction -= (n & DSigmaCorr);
 521         }
 522
 523         // Calculate the normal gradient based on the traction
 524
 525         gradient =
 526             Traction
 527             - (n & (mu*gradField.T() - (mu + lambda)*gradField))
 528             - n*lambda*tr(gradField);
 529
 530         //- Plasticity contribution
 531         if (rheology.plasticityActive())
 532         {
 533             gradient +=
 534                 2*mu*(n & rheology.DEpsilonP().boundaryField()[patch.index()]);
 535         }
 536
 537         // Thermal effects
 538         if
 539         (
 540             patch.boundaryMesh().mesh().objectRegistry::
 541             foundObject<thermalModel>("thermalProperties")
 542         )
 543         {
 544             const thermalModel& thermo =
 545                 patch.boundaryMesh().mesh().objectRegistry::
 546                 lookupObject<thermalModel>("thermalProperties");
 547
 548             const fvPatchScalarField& threeKalpha =
 549                 patch.lookupPatchField<volScalarField, scalar>
 550                 (
 551                     "((threeK*rho)*alpha)"
 552                 );
 553
 554                 if (!incremental)
 555                 {
 556                     const fvPatchScalarField& DT =
 557                         patch.lookupPatchField<volScalarField, scalar>("DT");
 558
 559                     gradient += n*threeKalpha*DT;
 560                 }
 561                 else
 562                 {
 563                     const fvPatchScalarField& T =
 564                         patch.lookupPatchField<volScalarField, scalar>("T");
 565
 566                     const scalarField T0 =
 567                         thermo.T0()().boundaryField()[patch.index()];
 568
 569                     gradient += n*threeKalpha*(T - T0);
 570                 }
 571         }
 572
 573         // Higher order non-linear terms
 574         if
 575         (
 576             nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
 577          || nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
 578         )
 579         {
 580             // no extra relaxation
 581             gradient -=
 582                 (n & (mu*(gradField & gradField.T())))
 583                 // + 0.5*n*lambda*(gradField && gradField);
 584               + 0.5*n*lambda*tr(gradField & gradField.T());
 585
 586             if
 587             (
 588                 incremental
 589              && nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
 590             )
 591             {
 592                 // gradU is const in a time step
 593                 const fvPatchTensorField& gradU =
 594                     patch.lookupPatchField<volTensorField, tensor>("grad(U)");
 595
 596                 gradient -=
 597                     (
 598                         n &
 599                         (
 600                             mu*
 601                             (
 602                                 (gradField & gradU.T())
 603                               + (gradU & gradField.T())
 604                             )
 605                         )
 606                     )
 607                   + 0.5*n*lambda*tr
 608                     (
 609                         (gradField & gradU.T())
 610                       + (gradU & gradField.T())
 611                     );
 612             }
 613         }
 614
 615         gradient /= (2.0*mu + lambda);
 616     }
 617
 618     return tgradient;
 619 }
 620
 621
 622 // ************************************************************************* //