applications/solvers/solidMechanics/elasticNonLinIncrTLSolidFoam/elasticNonLinIncrTLSolidFoam.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     elasticNonLinIncrTLSolidFoam
  26
  27 Description
  28     Finite volume structural solver employing an incremental strain total
  29     Lagrangian approach.
  30     For elastic solids.
  31     Valid for finite strains, finite displacements and finite rotations.
  32
  33 Author
  34     Philip Cardiff UCD
  35     Micheal Leonard UCD
  36
  37 \*---------------------------------------------------------------------------*/
  38
  39 #include "fvCFD.H"
  40 #include "constitutiveModel.H"
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 int main(int argc, char *argv[])
  45 {
  46 #   include "setRootCase.H"
  47 #   include "createTime.H"
  48 #   include "createMesh.H"
  49 #   include "createFields.H"
  50
  51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  52
  53     Info<< "\nStarting time loop\n" << endl;
  54
  55     while(runTime.loop())
  56     {
  57         Info<< "Time = " << runTime.timeName() << nl << endl;
  58
  59 #       include "readSolidMechanicsControls.H"
  60
  61         int iCorr = 0;
  62         scalar initialResidual = 0;
  63         lduMatrix::solverPerformance solverPerf;
  64         scalar relativeResidual = 1.0;
  65         lduMatrix::debug=0;
  66
  67         do
  68         {
  69             DU.storePrevIter();
  70
  71             fvVectorMatrix DUEqn
  72             (
  73                 fvm::d2dt2(rho, DU)
  74              ==
  75                 fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
  76               + fvc::div
  77                 (
  78                   - ( (mu + lambda) * gradDU )
  79                   + (
  80                         mu *
  81                         (
  82                             gradDU.T()
  83                           + (gradDU & gradU.T())
  84                           + (gradU & gradDU.T())
  85                           + (gradDU & gradDU.T())
  86                         )
  87                     )
  88                   + ( lambda * tr(DEpsilon) * I )
  89                   + ( DSigma & gradU )
  90                   + ( (sigma + DSigma) & gradDU ),
  91                     "div(sigma)"
  92                 )
  93             );
  94
  95             solverPerf = DUEqn.solve();
  96
  97             if (iCorr == 0)
  98             {
  99                 initialResidual = solverPerf.initialResidual();
 100             }
 101
 102             DU.relax();
 103
 104             gradDU = fvc::grad(DU);
 105
 106 #           include "calculateDEpsilonDSigma.H"
 107 #           include "calculateRelativeResidual.H"
 108
 109             Info<< "\tTime " << runTime.value()
 110                 << ", Corrector " << iCorr
 111                 << ", Solving for " << DU.name()
 112                 << " using " << solverPerf.solverName()
 113                 << ", residual = " << solverPerf.initialResidual()
 114                 << ", relative residual = " << relativeResidual
 115                 << ", inner iterations = " << solverPerf.nIterations() << endl;
 116         }
 117         while
 118         (
 119             solverPerf.initialResidual() > convergenceTolerance
 120             //relativeResidual > convergenceTolerance
 121          && ++iCorr < nCorr
 122         );
 123
 124         Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()
 125             << ", Initial residual = " << initialResidual
 126             << ", Final residual = " << solverPerf.initialResidual()
 127             << ", Relative residual = " << relativeResidual
 128             << ", No outer iterations " << iCorr
 129             << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 130             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 131             << endl;
 132
 133         U += DU;
 134         gradU += gradDU;
 135         epsilon += DEpsilon;
 136         sigma += DSigma;
 137
 138 #       include "writeFields.H"
 139
 140         Info<< "ExecutionTime = "
 141             << runTime.elapsedCpuTime()
 142             << " s\n\n" << endl;
 143     }
 144
 145     Info<< "End\n" << endl;
 146
 147     return(0);
 148 }
 149
 150
 151 // ************************************************************************* //