applications/solvers/solidMechanics/elasticPlasticNonLinTLSolidFoam/elasticPlasticNonLinTLSolidFoam.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 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<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 202             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 203             << endl;
 204     }
 205
 206     Info<< "End\n" << endl;
 207
 208     return(0);
 209 }
 210
 211
 212 // ************************************************************************* //