applications/solvers/solidMechanics/elasticPlasticNonLinTLSolidFoam/elasticPlasticNonLinTLSolidFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   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 Application
  25     elasticPlasticNonLinTLSolidFoam
  26
  27 Description
  28     Finite volume structural solver employing an incremental strain total
  29     Lagrangian approach, with Mises plasticity.
  30
  31     Valid for finite strains, finite displacements and finite rotations.
  32
  33 Author
  34     Philip Cardiff UCD
  35     Aitken relaxation by T. Tang DTU
  36
  37 \*---------------------------------------------------------------------------*/
  38
  39 #include "fvCFD.H"
  40 #include "constitutiveModel.H"
  41 #include "solidContactFvPatchVectorField.H"
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 int main(int argc, char *argv[])
  46 {
  47 # include "setRootCase.H"
  48 # include "createTime.H"
  49 # include "createMesh.H"
  50 # include "createFields.H"
  51 # include "createHistory.H"
  52 # include "readDivDSigmaExpMethod.H"
  53 # include "readDivDSigmaNonLinExpMethod.H"
  54
  55 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  56
  57   Info<< "\nStarting time loop\n" << endl;
  58
  59   while(runTime.loop())
  60     {
  61       Info<< "Time: " << runTime.timeName() << nl << endl;
  62
  63 #     include "readSolidMechanicsControls.H"
  64
  65       int iCorr = 0;
  66       scalar initialResidual = 0;
  67       lduMatrix::solverPerformance solverPerf;
  68       scalar relativeResidual = GREAT;
  69       lduMatrix::debug = 0;
  70
  71       do
  72       {
  73           DU.storePrevIter();
  74
  75 #         include "calculateDivDSigmaExp.H"
  76 #         include "calculateDivDSigmaNonLinExp.H"
  77
  78           // Incremental form of the
  79           // linear momentum conservation
  80           // ensuring conservation of total momentum
  81           fvVectorMatrix DUEqn
  82           (
  83               fvm::d2dt2(rho, DU)
  84            ==
  85               fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
  86             + divDSigmaExp
  87             + divDSigmaNonLinExp
  88           //- fvc::div(2*mu*DEpsilonP, "div(sigma)")
  89             - fvc::div
  90               (
  91                   2*muf*(mesh.Sf() & fvc::interpolate(DEpsilonP))
  92               )
  93           );
  94
  95           if (largeStrainOverRelax)
  96           {
  97               // the terms (gradDU & gradU.T()) and (gradU & gradDU.T())
  98               // are linearly dependent of DU and represent initial
  99               // displacement effect
 100               // which can cause convergence difficulties when treated
 101               // explicitly
 102               // so we implicitly over-relax with gradU & gradDU here
 103               // which tends to help convergence
 104               // this should improve convergence when gradU is large
 105               // but maybe not execution time
 106               DUEqn -=
 107                   fvm::laplacian
 108                   (
 109                       (2*mu + lambda)*gradU, DU, "laplacian(DDU,DU)"
 110                   )
 111                 - fvc::div((2*mu + lambda)*(gradU & gradDU), "div(sigma)");
 112           }
 113
 114           if (nonLinearSemiImplicit)
 115           {
 116               // experimental
 117               // we can treat the nonlinear term (gradDU & gradDU.T()) in a
 118               // semi-implicit over-relaxed manner
 119               // this should improve convergence when gradDU is large
 120               // but maybe not execution time
 121               DUEqn -=
 122                   fvm::laplacian
 123                   (
 124                       (2*mu + lambda)*gradDU, DU, "laplacian(DDU,DU)"
 125                   )
 126                 - fvc::div((2*mu + lambda)*(gradDU & gradDU), "div(sigma)");
 127           }
 128
 129           solverPerf = DUEqn.solve();
 130
 131           if (iCorr == 0)
 132           {
 133               initialResidual = solverPerf.initialResidual();
 134           }
 135
 136           if (aitkenRelax)
 137           {
 138 #             include "aitkenRelaxation.H"
 139           }
 140           else
 141           {
 142               DU.relax();
 143           }
 144
 145           gradDU = fvc::grad(DU);
 146
 147           // correct plasticty term
 148           rheology.correct();
 149
 150 #         include "calculateDEpsilonDSigma.H"
 151 #         include "calculateRelativeResidual.H"
 152
 153           if (iCorr % infoFrequency == 0)
 154           {
 155               Info<< "\tTime " << runTime.value()
 156                   << ", Corrector " << iCorr
 157                   << ", Solving for " << DU.name()
 158                   << " using " << solverPerf.solverName()
 159                   << ", res = " << solverPerf.initialResidual()
 160                   << ", rel res = " << relativeResidual;
 161
 162               if (aitkenRelax)
 163               {
 164                   Info << ", aitken = " << aitkenTheta;
 165               }
 166               Info << ", iters = " << solverPerf.nIterations() << endl;
 167           }
 168       }
 169       while
 170       (
 171           iCorr++ == 0
 172        ||
 173           (
 174               //solverPerf.initialResidual() > convergenceTolerance
 175               relativeResidual > convergenceTolerance
 176            && iCorr < nCorr
 177           )
 178       );
 179
 180       Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()
 181           << ", Initial residual = " << initialResidual
 182           << ", Final residual = " << solverPerf.initialResidual()
 183           << ", Relative residual = " << relativeResidual
 184           << ", No outer iterations " << iCorr
 185           << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 186           << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 187           << endl;
 188
 189       // Update total quantities
 190       U += DU;
 191       gradU += gradDU;
 192       epsilon += DEpsilon;
 193       epsilonP += rheology.DEpsilonP();
 194       sigma += DSigma;
 195       rheology.updateYieldStress();
 196       rho = rho/det(I+gradU);
 197
 198 #     include "writeFields.H"
 199 #     include "writeHistory.H"
 200
 201       Info<< "ExecutionTime = "
 202           << runTime.elapsedCpuTime()
 203           << " s\n\n" << endl;
 204     }
 205
 206   Info<< "End\n" << endl;
 207
 208   return(0);
 209 }
 210
 211
 212 // ************************************************************************* //