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