applications/solvers/solidMechanics/elasticOrthoNonLinULSolidFoam/elasticOrthoNonLinULSolidFoam.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     elasticOrthoGenDirULSolidFoam
  26
  27 Description
  28     Transient/steady-state segregated finite-volume solver for large strain
  29     elastic orthotropic solid bodies allowing for general principal material
  30     directions.
  31
  32     Displacement increment field DU is solved for using an updated Lagrangian
  33     approach,
  34     also generating the Almansi strain tensor field epsilon and Cauchy stress
  35     tensor field sigma.
  36
  37     At the end of each time-step, the mesh is moved and sigma, epsilon and C
  38     are rotated to the new configuration.
  39
  40     Please cite:
  41     Cardiff P, Karac A & Ivankovic A, A Large Strain Finite Volume Method for
  42     Orthotropic Bodies with General Material Orientations, Computer Methods
  43     in Applied Mechanics & Engineering, Sep 2013,
  44     http://dx.doi.org/10.1016/j.cma.2013.09.008
  45
  46 Author
  47     Philip Cardiff UCD
  48
  49 \*---------------------------------------------------------------------------*/
  50
  51 #include "fvCFD.H"
  52 #include "constitutiveModel.H"
  53 #include "transformField.H"
  54 #include "transformGeometricField.H"
  55 #include "pointPatchInterpolation.H"
  56 #include "primitivePatchInterpolation.H"
  57 #include "pointFields.H"
  58 #include "twoDPointCorrector.H"
  59 #include "leastSquaresVolPointInterpolation.H"
  60
  61 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  62
  63 int main(int argc, char *argv[])
  64 {
  65 # include "setRootCase.H"
  66 # include "createTime.H"
  67 # include "createMesh.H"
  68 # include "createFields.H"
  69
  70 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  71
  72   Info<< "\nStarting time loop\n" << endl;
  73
  74   for (runTime++; !runTime.end(); runTime++)
  75     {
  76       Info<< "Time: " << runTime.timeName() << nl << endl;
  77
  78 #     include "readSolidMechanicsControls.H"
  79
  80       int iCorr = 0;
  81       lduMatrix::solverPerformance solverPerf;
  82       scalar initialResidual = 1.0;
  83       lduMatrix::debug = 0;
  84
  85       //- div(sigmaOld) should be zero but I will include
  86       //- it to make sure errors don't accumulate
  87       volVectorField* oldErrorPtr = NULL;
  88       if (ensureTotalEquilibrium)
  89       {
  90           oldErrorPtr = new volVectorField
  91               (
  92                   fvc::d2dt2(rho.oldTime(), U.oldTime())
  93                   - fvc::div(sigma)
  94                   );
  95       }
  96
  97       do
  98       {
  99           DU.storePrevIter();
 100
 101           //- Updated lagrangian momentum equation
 102           fvVectorMatrix DUEqn
 103               (
 104                   fvm::d2dt2(rho, DU)
 105                   + fvc::d2dt2(rho, U)
 106                   ==
 107                   fvm::laplacian(K, DU, "laplacian(K,DU)")
 108                   + fvc::div(
 109                       DSigma
 110                       - (K & gradDU)
 111                       + ( (sigma + DSigma) & gradDU ),
 112                       "div(sigma)"
 113                       )
 114                   //- fvc::laplacian(K, DU)
 115                   );
 116
 117           if (ensureTotalEquilibrium)
 118           {
 119               //- to stop accumulation of errors
 120               DUEqn += *oldErrorPtr;
 121           }
 122
 123           solverPerf = DUEqn.solve();
 124
 125           if (iCorr == 0)
 126           {
 127               initialResidual = solverPerf.initialResidual();
 128           }
 129
 130           DU.relax();
 131
 132           gradDU = fvc::grad(DU);
 133
 134           //- for 2-D plane stress simulations, the zz component of gradDU
 135           //- ensures sigma.zz() is zero
 136           //- it is assumed that z is the empty direction
 137           //#         include "checkPlaneStress.H"
 138
 139           //- sigma needs to be calculated inside the momentum loop as
 140           //- it is used in the momentum equation
 141           DEpsilon = symm(gradDU) + 0.5*symm(gradDU & gradDU.T());
 142           DSigma = C && DEpsilon;
 143
 144           if (iCorr % infoFrequency == 0)
 145           {
 146               Info << "\tTime " << runTime.value()
 147                    << ", Corr " << iCorr
 148                    << ", Solving for " << DU.name()
 149                    << " using " << solverPerf.solverName()
 150                    << ", res = " << solverPerf.initialResidual()
 151                   //<< ", rel res = " << relativeResidual
 152                    << ", inner iters " << solverPerf.nIterations() << endl;
 153           }
 154       }
 155       while
 156           (
 157               solverPerf.initialResidual() > convergenceTolerance
 158               &&
 159               ++iCorr < nCorr
 160               );
 161
 162       Info << nl << "Time " << runTime.value() << ", Solving for " << DU.name()
 163            << ", Initial residual = " << initialResidual
 164            << ", Final residual = " << solverPerf.initialResidual()
 165            << ", No outer iterations " << iCorr
 166            << nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 167            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 168            << endl;
 169
 170 #     include "moveMeshLeastSquares.H"
 171 #     include "rotateFields.H"
 172 #     include "writeFields.H"
 173
 174       Info<< "ExecutionTime = "
 175           << runTime.elapsedCpuTime()
 176           << " s\n\n" << endl;
 177     }
 178
 179   Info<< "End\n" << endl;
 180
 181   return(0);
 182 }
 183
 184
 185 // ************************************************************************* //