applications/solvers/solidMechanics/elasticPlasticNonLinULSolidFoam/elasticPlasticNonLinULSolidFoam.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     elasticPlasticNonLinULSolidFoam
  26
  27 Description
  28     Finite volume structural solver employing a incremental strain updated
  29     Lagrangian approach.
  30
  31     Valid for small strains, finite displacements and finite rotations.
  32
  33     Note: the reason the solver is not strictly valid for large strains is
  34     because the constitutive stiffness tensor is not rotated.
  35     For an updated Lagrangian solver which does rotate the stiffness tensor, and
  36     hence is strictly Valid for large strains, use
  37     elasticOrthoNonLinULSolidFoam.
  38
  39 Author
  40     Philip Cardiff UCD
  41     Aitken relaxation by T. Tang DTU
  42
  43 \*---------------------------------------------------------------------------*/
  44
  45 #include "fvCFD.H"
  46 #include "constitutiveModel.H"
  47 #include "volPointInterpolation.H"
  48 #include "pointPatchInterpolation.H"
  49 #include "primitivePatchInterpolation.H"
  50 #include "pointFields.H"
  51 #include "twoDPointCorrector.H"
  52 #include "leastSquaresVolPointInterpolation.H"
  53 #include "transformGeometricField.H"
  54 #include "solidContactFvPatchVectorField.H"
  55 #include "pointMesh.H"
  56 #include "symmetryPolyPatch.H"
  57
  58 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  59
  60 int main(int argc, char *argv[])
  61 {
  62 # include "setRootCase.H"
  63 # include "createTime.H"
  64 # include "createMesh.H"
  65 # include "createFields.H"
  66 # include "createHistory.H"
  67 # include "readDivDSigmaExpMethod.H"
  68 # include "readDivDSigmaNonLinExpMethod.H"
  69 # include "readMoveMeshMethod.H"
  70 # include "findGlobalFaceZones.H"
  71
  72 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  73
  74     Info<< "\nStarting time loop\n" << endl;
  75
  76     while(runTime.loop())
  77     {
  78         Info<< "Time = " << runTime.timeName() << nl << endl;
  79
  80 #       include "readSolidMechanicsControls.H"
  81
  82         int iCorr = 0;
  83         lduMatrix::solverPerformance solverPerf;
  84         scalar initialResidual = 0;
  85         scalar relativeResidual = 1.0;
  86         lduMatrix::debug = 0;
  87
  88
  89         do
  90         {
  91             DU.storePrevIter();
  92
  93 #           include "calculateDivDSigmaExp.H"
  94 #           include "calculateDivDSigmaNonLinExp.H"
  95
  96             // Updated lagrangian large strain momentum equation
  97             fvVectorMatrix DUEqn
  98             (
  99                 fvm::d2dt2(rho,DU)
 100              ==
 101                 fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
 102               + divDSigmaExp
 103               + divDSigmaNonLinExp
 104               //- fvc::div(2*mu*DEpsilonP, "div(sigma)")
 105               - fvc::div(2*muf*( mesh.Sf() & fvc::interpolate(DEpsilonP)) )
 106             );
 107
 108             if(nonLinearSemiImplicit)
 109             {
 110                 // experimental
 111                 // we can treat the nonlinear term (gradDU & gradDU.T()) in a
 112                 // semi-implicit over-relaxed manner
 113                 // this should improve convergence when gradDU is large
 114                 // but maybe not execution time
 115                 DUEqn -=
 116                     fvm::laplacian
 117                     (
 118                         (2*mu + lambda)*gradDU, DU, "laplacian(DDU,DU)"
 119                     )
 120                   - fvc::div( (2*mu + lambda)*(gradDU&gradDU), "div(sigma)");
 121             }
 122
 123             solverPerf = DUEqn.solve();
 124
 125             if(iCorr == 0)
 126             {
 127                 initialResidual = solverPerf.initialResidual();
 128             }
 129
 130             if(aitkenRelax)
 131             {
 132 #               include "aitkenRelaxation.H"
 133             }
 134             else
 135             {
 136                 DU.relax();
 137             }
 138
 139             gradDU = fvc::grad(DU);
 140
 141             // correct plasticty term
 142             rheology.correct();
 143
 144             // correct elastic properties
 145             // for nonlinear elastic materials
 146             //mu = rheology.newMu();
 147             //lambda = rheology.newLambda();
 148             //muf = fvc::interpolate(mu);
 149             //lambdaf = fvc::interpolate(lambda);
 150
 151 #           include "calculateDEpsilonDSigma.H"
 152 #           include "calculateRelativeResidual.H"
 153
 154             if(iCorr % infoFrequency == 0)
 155             {
 156                 Info<< "\tTime " << runTime.value()
 157                     << ", Corrector " << iCorr
 158                     << ", Solving for " << DU.name()
 159                     << " using " << solverPerf.solverName()
 160                     << ", res = " << solverPerf.initialResidual()
 161                     << ", rel res = " << relativeResidual;
 162                 if(aitkenRelax)
 163                 {
 164                     Info<< ", aitken = " << aitkenTheta;
 165                 }
 166                 Info<< ", iters = " << solverPerf.nIterations() << endl;
 167             }
 168         }
 169         while
 170         (
 171             iCorr++ < 2
 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             << ", Final rel residual = " << relativeResidual
 184             << ", No outer iterations " << iCorr << endl;
 185
 186         rheology.updateYieldStress();
 187
 188         U += DU;
 189         epsilon += DEpsilon;
 190         epsilonP += DEpsilonP;
 191         sigma += DSigma;
 192
 193 #       include "moveMesh.H"
 194 #       include "rotateFields.H"
 195 #       include "writeFields.H"
 196 #       include "writeHistory.H"
 197
 198         Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 199             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 200             << endl;
 201     }
 202
 203     Info<< "End\n" << endl;
 204
 205     return(0);
 206 }
 207
 208
 209 // ************************************************************************* //