applications/solvers/solidMechanics/elasticPlasticNonLinTLSolidFoam/elasticPlasticNonLinTLSolidFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2007 Hrvoje Jasak
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM 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 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  24
  25 Application
  26     elasticPlasticNonLinTLSolidFoam
  27
  28 Description
  29     Finite volume structural solver employing an incremental strain total
  30     Lagrangian approach, with Mises plasticity.
  31
  32     Valid for finite strains, finite displacements and finite rotations.
  33
  34 Author
  35     Philip Cardiff UCD
  36     Aitken relaxation by T. Tang DTU
  37
  38 \*---------------------------------------------------------------------------*/
  39
  40 #include "fvCFD.H"
  41 #include "constitutiveModel.H"
  42 #include "solidContactFvPatchVectorField.H"
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 int main(int argc, char *argv[])
  47 {
  48 # include "setRootCase.H"
  49 # include "createTime.H"
  50 # include "createMesh.H"
  51 # include "createFields.H"
  52 # include "createHistory.H"
  53 # include "readDivDSigmaExpMethod.H"
  54 # include "readDivDSigmaNonLinExpMethod.H"
  55
  56 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  57
  58   Info<< "\nStarting time loop\n" << endl;
  59
  60   while(runTime.loop())
  61     {
  62       Info<< "Time: " << runTime.timeName() << nl << endl;
  63
  64 #     include "readSolidMechanicsControls.H"
  65
  66       int iCorr = 0;
  67       scalar initialResidual = 0;
  68       lduMatrix::solverPerformance solverPerf;
  69       scalar relativeResidual = GREAT;
  70       lduMatrix::debug=0;
  71
  72       do
  73         {
  74             DU.storePrevIter();
  75
  76 #         include "calculateDivDSigmaExp.H"
  77 #         include "calculateDivDSigmaNonLinExp.H"
  78
  79             // Incremental form of the
  80             // linear momentum conservation
  81             // ensuring conservation of total momentum
  82             fvVectorMatrix DUEqn
  83                 (
  84                     fvm::d2dt2(rho, DU)
  85                     ==
  86                     fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
  87                     + divDSigmaExp
  88                     + divDSigmaNonLinExp
  89                     //- fvc::div(2*mu*DEpsilonP, "div(sigma)")
  90                     - fvc::div
  91                     (
  92                         2*muf*( mesh.Sf() & fvc::interpolate(DEpsilonP))
  93                         )
  94                     );
  95
  96             if (largeStrainOverRelax)
  97             {
  98                 // the terms (gradDU & gradU.T()) and (gradU & gradDU.T())
  99                 // are linearly dependent of DU and represent initial
 100                 // displacement effect
 101                 // which can cause convergence difficulties when treated
 102                 // explicitly
 103                 // so we implicitly over-relax with gradU & gradDU here
 104                 // which tends to help convergence
 105                 // this should improve convergence when gradU is large
 106                 // but maybe not execution time
 107                 DUEqn -=
 108                     fvm::laplacian
 109                     (
 110                         (2*mu + lambda)*gradU, DU, "laplacian(DDU,DU)"
 111                         )
 112                     - fvc::div( (2*mu + lambda)*(gradU&gradDU), "div(sigma)");
 113         }
 114
 115             if (nonLinearSemiImplicit)
 116             {
 117                 // experimental
 118                 // we can treat the nonlinear term (gradDU & gradDU.T()) in a
 119                 // semi-implicit over-relaxed manner
 120                 // this should improve convergence when gradDU is large
 121                 // but maybe not execution time
 122                 DUEqn -=
 123                     fvm::laplacian
 124                     (
 125                         (2*mu + lambda)*gradDU, DU, "laplacian(DDU,DU)"
 126                         )
 127                     - fvc::div( (2*mu + lambda)*(gradDU&gradDU), "div(sigma)");
 128         }
 129
 130             solverPerf = DUEqn.solve();
 131
 132             if (iCorr == 0)
 133             {
 134                 initialResidual = solverPerf.initialResidual();
 135             }
 136
 137             if (aitkenRelax)
 138             {
 139 #             include "aitkenRelaxation.H"
 140             }
 141             else
 142             {
 143                 DU.relax();
 144             }
 145
 146             gradDU = fvc::grad(DU);
 147
 148             // correct plasticty term
 149             rheology.correct();
 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++ == 0
 172               ||
 173               (//solverPerf.initialResidual() > convergenceTolerance
 174                   relativeResidual > convergenceTolerance
 175                   &&
 176                   iCorr < nCorr)
 177               );
 178
 179       Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()
 180           << ", Initial residual = " << initialResidual
 181           << ", Final residual = " << solverPerf.initialResidual()
 182           << ", Relative residual = " << relativeResidual
 183           << ", No outer iterations " << iCorr
 184           << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 185           << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 186           << endl;
 187
 188       // Update total quantities
 189       U += DU;
 190       gradU += gradDU;
 191       epsilon += DEpsilon;
 192       epsilonP += rheology.DEpsilonP();
 193       sigma += DSigma;
 194       rheology.updateYieldStress();
 195       rho = rho/det(I+gradU);
 196
 197 #     include "writeFields.H"
 198 #     include "writeHistory.H"
 199
 200       Info<< "ExecutionTime = "
 201           << runTime.elapsedCpuTime()
 202           << " s\n\n" << endl;
 203     }
 204
 205   Info<< "End\n" << endl;
 206
 207   return(0);
 208 }
 209
 210
 211 // ************************************************************************* //