src/solidModels/constitutiveModel/rheologyLaws/multiMaterial/multiMaterial.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 Description
  25     Zoned multi-material rheology controlled by a material indicator field.
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "multiMaterial.H"
  30 #include "addToRunTimeSelectionTable.H"
  31 #include "zeroGradientFvPatchFields.H"
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37     defineTypeNameAndDebug(multiMaterial, 0);
  38     addToRunTimeSelectionTable(rheologyLaw, multiMaterial, dictionary);
  39 }
  40
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44 Foam::tmp<Foam::scalarField> Foam::multiMaterial::indicator
  45 (
  46     const label i
  47 ) const
  48 {
  49     const scalarField& mat = materials_.internalField();
  50
  51     tmp<scalarField> tresult(new scalarField(mat.size(), 0.0));
  52     scalarField& result = tresult();
  53
  54     forAll (mat, matI)
  55     {
  56         if (mat[matI] > i - SMALL && mat[matI] < i + 1 - SMALL)
  57         {
  58             result[matI] = 1.0;
  59         }
  60     }
  61
  62     return tresult;
  63 }
  64
  65
  66 Foam::scalar
  67 Foam::multiMaterial::indicator(const label index, const label cellID) const
  68 {
  69     const scalar mat = materials_.internalField()[cellID];
  70     scalar result = 0.0;
  71
  72     if (mat > index - SMALL && mat < index + 1 - SMALL)
  73       {
  74     result = 1.0;
  75       }
  76
  77     return result;
  78 }
  79
  80 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  81
  82 // Construct from dictionary
  83 Foam::multiMaterial::multiMaterial
  84 (
  85     const word& name,
  86     const volSymmTensorField& sigma,
  87     const dictionary& dict
  88 )
  89 :
  90     rheologyLaw(name, sigma, dict),
  91     PtrList<rheologyLaw>(),
  92     materials_
  93     (
  94         IOobject
  95         (
  96             "materials",
  97             mesh().time().timeName(),
  98             mesh(),
  99             IOobject::MUST_READ,
 100             IOobject::AUTO_WRITE
 101         ),
 102         mesh()
 103     )
 104 {
 105     PtrList<rheologyLaw>& laws = *this;
 106
 107     PtrList<entry> lawEntries(dict.lookup("laws"));
 108     laws.setSize(lawEntries.size());
 109
 110     forAll (laws, lawI)
 111     {
 112         laws.set
 113         (
 114             lawI,
 115             rheologyLaw::New
 116             (
 117                 lawEntries[lawI].keyword(),
 118                 sigma,
 119                 lawEntries[lawI].dict()
 120             )
 121         );
 122     }
 123
 124     if
 125     (
 126         min(materials_).value() < 0
 127      || max(materials_).value() > laws.size() + SMALL
 128     )
 129     {
 130         FatalErrorIn
 131         (
 132             "multiMaterial::multiMaterial\n"
 133             "(\n"
 134             "    const word& name,\n"
 135             "    const volSymmTensorField& sigma,\n"
 136             "    const dictionary& dict\n"
 137             ")"
 138         )   << "Invalid definition of material indicator field.  "
 139             << "Number of materials: " << laws.size()
 140             << " max index: " << max(materials_)
 141             << ".  Should be " << laws.size() - 1
 142             << abort(FatalError);
 143     }
 144 }
 145
 146
 147 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 148
 149 Foam::multiMaterial::~multiMaterial()
 150 {}
 151
 152
 153 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 154
 155 Foam::tmp<Foam::volScalarField> Foam::multiMaterial::rho() const
 156 {
 157     tmp<volScalarField> tresult
 158     (
 159         new volScalarField
 160         (
 161             IOobject
 162             (
 163                 "rho",
 164                 mesh().time().timeName(),
 165                 mesh(),
 166                 IOobject::NO_READ,
 167                 IOobject::NO_WRITE
 168             ),
 169             mesh(),
 170             dimensionedScalar("zeroRho", dimDensity, 0),
 171             zeroGradientFvPatchScalarField::typeName
 172         )
 173     );
 174     volScalarField& result = tresult();
 175
 176     // Accumulate data for all fields
 177     const PtrList<rheologyLaw>& laws = *this;
 178
 179     forAll (laws, lawI)
 180     {
 181         result.internalField() +=
 182             indicator(lawI)*laws[lawI].rho()().internalField();
 183     }
 184
 185     result.correctBoundaryConditions();
 186
 187     return tresult;
 188 }
 189
 190
 191 Foam::tmp<Foam::volScalarField> Foam::multiMaterial::E() const
 192 {
 193     tmp<volScalarField> tresult
 194     (
 195         new volScalarField
 196         (
 197             IOobject
 198             (
 199                 "E",
 200                 mesh().time().timeName(),
 201                 mesh(),
 202                 IOobject::NO_READ,
 203                 IOobject::NO_WRITE
 204             ),
 205             mesh(),
 206             dimensionedScalar("zeroE", dimForce/dimArea, 0),
 207             zeroGradientFvPatchScalarField::typeName
 208         )
 209     );
 210     volScalarField& result = tresult();
 211
 212     // Accumulate data for all fields
 213     const PtrList<rheologyLaw>& laws = *this;
 214
 215     forAll (laws, lawI)
 216     {
 217         result.internalField() +=
 218             indicator(lawI)*laws[lawI].E()().internalField();
 219     }
 220
 221     result.correctBoundaryConditions();
 222
 223     return tresult;
 224 }
 225
 226
 227 Foam::tmp<Foam::volScalarField>
 228 Foam::multiMaterial::E(const volScalarField& epsEq) const
 229 {
 230     tmp<volScalarField> tresult
 231     (
 232         new volScalarField
 233         (
 234             IOobject
 235             (
 236                 "E",
 237                 mesh().time().timeName(),
 238                 mesh(),
 239                 IOobject::NO_READ,
 240                 IOobject::AUTO_WRITE
 241             ),
 242             mesh(),
 243             dimensionedScalar("zeroE", dimForce/dimArea, 0),
 244             zeroGradientFvPatchScalarField::typeName
 245         )
 246     );
 247     volScalarField& result = tresult();
 248
 249     // Accumulate data for all fields
 250     const PtrList<rheologyLaw>& laws = *this;
 251
 252     forAll (laws, lawI)
 253     {
 254         result.internalField() +=
 255             indicator(lawI)*laws[lawI].E(epsEq)().internalField();
 256     }
 257
 258 //     forAll(result.boundaryField(),patchI)
 259 //     {
 260 //         forAll (laws, lawI)
 261 //         {
 262 //             result.boundaryField()[patchI] +=
 263 //                 indicator(lawI)().boundaryField()[patchI]
 264 //                 *laws[lawI].E(t)().boundaryField()[patchI];
 265 //         }
 266 //     }
 267
 268     result.correctBoundaryConditions();
 269
 270     return tresult;
 271 }
 272
 273
 274 Foam::tmp<Foam::volScalarField> Foam::multiMaterial::nu() const
 275 {
 276     tmp<volScalarField> tresult
 277     (
 278         new volScalarField
 279         (
 280             IOobject
 281             (
 282                 "nu",
 283                 mesh().time().timeName(),
 284                 mesh(),
 285                 IOobject::NO_READ,
 286                 IOobject::NO_WRITE
 287             ),
 288             mesh(),
 289             dimensionedScalar("zeroE", dimless, 0),
 290             zeroGradientFvPatchScalarField::typeName
 291         )
 292     );
 293     volScalarField& result = tresult();
 294
 295     // Accumulate data for all fields
 296     const PtrList<rheologyLaw>& laws = *this;
 297
 298     forAll (laws, lawI)
 299     {
 300         result.internalField() +=
 301             indicator(lawI)*laws[lawI].nu()().internalField();
 302     }
 303
 304     result.correctBoundaryConditions();
 305
 306     return tresult;
 307 }
 308
 309
 310 Foam::tmp<Foam::volScalarField> Foam::multiMaterial::Ep() const
 311 {
 312     tmp<volScalarField> tresult
 313     (
 314         new volScalarField
 315         (
 316             IOobject
 317             (
 318                 "Ep",
 319                 mesh().time().timeName(),
 320                 mesh(),
 321                 IOobject::NO_READ,
 322                 IOobject::NO_WRITE
 323             ),
 324             mesh(),
 325             dimensionedScalar("zeroEp", dimForce/dimArea, 0.0),
 326             zeroGradientFvPatchScalarField::typeName
 327         )
 328     );
 329     volScalarField& result = tresult();
 330
 331     // Accumulate data for all fields
 332     const PtrList<rheologyLaw>& laws = *this;
 333
 334     forAll (laws, lawI)
 335     {
 336         result.internalField() +=
 337             indicator(lawI)*laws[lawI].Ep()().internalField();
 338     }
 339
 340     result.correctBoundaryConditions();
 341
 342     return tresult;
 343 }
 344
 345
 346 Foam::tmp<Foam::volScalarField>
 347 Foam::multiMaterial::Ep(const volScalarField& epsEq) const
 348 {
 349     tmp<volScalarField> tresult
 350     (
 351         new volScalarField
 352         (
 353             IOobject
 354             (
 355                 "Ep",
 356                 mesh().time().timeName(),
 357                 mesh(),
 358                 IOobject::NO_READ,
 359                 IOobject::AUTO_WRITE
 360             ),
 361             mesh(),
 362             dimensionedScalar("zeroEp", dimForce/dimArea, 0),
 363             zeroGradientFvPatchScalarField::typeName
 364         )
 365     );
 366     volScalarField& result = tresult();
 367
 368     // Accumulate data for all fields
 369     const PtrList<rheologyLaw>& laws = *this;
 370
 371     forAll (laws, lawI)
 372     {
 373         result.internalField() +=
 374             indicator(lawI)*laws[lawI].Ep(epsEq)().internalField();
 375     }
 376
 377 //     forAll(result.boundaryField(),patchI)
 378 //     {
 379 //         forAll (laws, lawI)
 380 //         {
 381 //             result.boundaryField()[patchI] +=
 382 //                 indicator(lawI)().boundaryField()[patchI]
 383 //                 *laws[lawI].E(t)().boundaryField()[patchI];
 384 //         }
 385 //     }
 386
 387     result.correctBoundaryConditions();
 388
 389     return tresult;
 390 }
 391
 392
 393 Foam::tmp<Foam::volScalarField> Foam::multiMaterial::sigmaY() const
 394 {
 395     tmp<volScalarField> tresult
 396     (
 397         new volScalarField
 398         (
 399             IOobject
 400             (
 401                 "sigmaY",
 402                 mesh().time().timeName(),
 403                 mesh(),
 404                 IOobject::NO_READ,
 405                 IOobject::NO_WRITE
 406             ),
 407             mesh(),
 408             dimensionedScalar("zeroSigmaY", dimForce/dimArea, 0.0),
 409             zeroGradientFvPatchScalarField::typeName
 410         )
 411     );
 412     volScalarField& result = tresult();
 413
 414     // Accumulate data for all fields
 415     const PtrList<rheologyLaw>& laws = *this;
 416
 417     forAll (laws, lawI)
 418     {
 419         result.internalField() +=
 420             indicator(lawI)*laws[lawI].sigmaY()().internalField();
 421     }
 422
 423     result.correctBoundaryConditions();
 424
 425     return tresult;
 426 }
 427
 428 Foam::scalar
 429 Foam::multiMaterial::sigmaY(const scalar epsilonPEq, const label cellID) const
 430 {
 431   // Accumulate data for all fields
 432   const PtrList<rheologyLaw>& laws = *this;
 433
 434   scalar result = 0.0;
 435   forAll (laws, lawI)
 436     {
 437       result +=
 438             indicator(lawI, cellID)*laws[lawI].sigmaY(epsilonPEq, cellID);
 439     }
 440
 441     return result;
 442 }
 443
 444 bool Foam::multiMaterial::plasticityModelNeeded() const
 445 {
 446   // Accumulate data for all fields
 447   const PtrList<rheologyLaw>& laws = *this;
 448
 449   bool active = false;
 450   forAll(laws, lawI)
 451     {
 452       if (laws[lawI].plasticityModelNeeded())
 453       {
 454           active = true;
 455           break;
 456       }
 457     }
 458
 459     return active;
 460 }
 461
 462 Foam::tmp<Foam::volDiagTensorField> Foam::multiMaterial::K() const
 463 {
 464     tmp<volDiagTensorField> tresult
 465     (
 466         new volDiagTensorField
 467         (
 468             IOobject
 469             (
 470                 "K",
 471                 mesh().time().timeName(),
 472                 mesh(),
 473                 IOobject::NO_READ,
 474                 IOobject::NO_WRITE
 475             ),
 476             mesh(),
 477             dimensionedDiagTensor("zeroK", dimForce/dimArea, diagTensor::zero),
 478             zeroGradientFvPatchScalarField::typeName
 479         )
 480     );
 481     volDiagTensorField& result = tresult();
 482
 483     // Accumulate data for all fields
 484     const PtrList<rheologyLaw>& laws = *this;
 485
 486     forAll (laws, lawI)
 487       {
 488     result.internalField() +=
 489       indicator(lawI)*laws[lawI].K()().internalField();
 490       }
 491
 492     result.correctBoundaryConditions();
 493
 494     return tresult;
 495 }
 496
 497 Foam::tmp<Foam::volSymmTensor4thOrderField> Foam::multiMaterial::C() const
 498 {
 499     tmp<volSymmTensor4thOrderField> tresult
 500     (
 501         new volSymmTensor4thOrderField
 502         (
 503             IOobject
 504             (
 505                 "C",
 506                 mesh().time().timeName(),
 507                 mesh(),
 508                 IOobject::NO_READ,
 509                 IOobject::NO_WRITE
 510             ),
 511             mesh(),
 512             dimensionedSymmTensor4thOrder
 513             ("zeroC", dimForce/dimArea, symmTensor4thOrder::zero),
 514             zeroGradientFvPatchScalarField::typeName
 515         )
 516     );
 517     volSymmTensor4thOrderField& result = tresult();
 518
 519     // Accumulate data for all fields
 520     const PtrList<rheologyLaw>& laws = *this;
 521
 522     forAll (laws, lawI)
 523     {
 524         result.internalField() +=
 525       indicator(lawI)*laws[lawI].C()().internalField();
 526     }
 527
 528     result.correctBoundaryConditions();
 529
 530     return tresult;
 531 }
 532
 533
 534 void Foam::multiMaterial::correct()
 535 {
 536     PtrList<rheologyLaw>& laws = *this;
 537
 538     forAll (laws, lawI)
 539     {
 540         laws[lawI].correct();
 541     }
 542 }
 543
 544
 545 // ************************************************************************* //