applications/solvers/solidMechanics/elasticThermalSolidFoam/elasticThermalSolidFoam.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     elasticThermalSolidFoam
  26
  27 Description
  28     Transient/steady-state segregated finite-volume solver for small strain
  29     elastic thermal solid bodies. Temperature is solved and then coupled
  30     displacement is solved.
  31
  32     Displacement field U is solved for using a total Lagrangian approach,
  33     also generating the strain tensor field epsilon and stress tensor
  34     field sigma and temperature field T.
  35
  36 Author
  37    Philip Cardiff UCD
  38
  39 \*---------------------------------------------------------------------------*/
  40
  41 #include "fvCFD.H"
  42 #include "constitutiveModel.H"
  43 #include "thermalModel.H"
  44
  45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  46
  47 int main(int argc, char *argv[])
  48 {
  49 #   include "setRootCase.H"
  50 #   include "createTime.H"
  51 #   include "createMesh.H"
  52 #   include "createFields.H"
  53 #   include "readDivSigmaExpMethod.H"
  54
  55 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  56
  57     Info<< "\nStarting time loop\n" << endl;
  58
  59     while(runTime.loop())
  60     {
  61         Info<< "Time = " << runTime.timeName() << nl << endl;
  62
  63 #       include "readSolidMechanicsControls.H"
  64
  65         int iCorr = 0;
  66         scalar initialResidual = 1.0;
  67         scalar relResT = 1.0;
  68         scalar relResU = 1.0;
  69         lduMatrix::solverPerformance solverPerfU;
  70         lduMatrix::solverPerformance solverPerfT;
  71         lduMatrix::debug = 0;
  72
  73         // solve energy equation for temperature
  74         // the loop is for non-orthogonal corrections
  75         Info<< "Solving for " << T.name() << nl;
  76         do
  77         {
  78             T.storePrevIter();
  79
  80             fvScalarMatrix TEqn
  81             (
  82                 rhoC*fvm::ddt(T) == fvm::laplacian(k, T, "laplacian(k,T)")
  83             );
  84
  85             solverPerfT = TEqn.solve();
  86
  87             T.relax();
  88
  89 #           include "calculateRelResT.H"
  90
  91             if (iCorr % infoFrequency == 0)
  92             {
  93                 Info<< "\tCorrector " << iCorr
  94                     << ", residual = " << solverPerfT.initialResidual()
  95                     << ", relative res = " << relResT
  96                     << ", inner iters = " << solverPerfT.nIterations() << endl;
  97             }
  98         }
  99         while
 100         (
 101             relResT > convergenceToleranceT
 102          && ++iCorr < nCorr
 103         );
 104
 105         Info<< "Solved for " << T.name()
 106             << " using " << solverPerfT.solverName()
 107             << " in " << iCorr << " iterations"
 108             << ", residual = " << solverPerfT.initialResidual()
 109             << ", relative res = " << relResT << nl
 110             << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 111             << ", ClockTime = " << runTime.elapsedClockTime() << " s"
 112             << endl;
 113
 114         // Solve momentum equation for displacement
 115         iCorr = 0;
 116         volVectorField gradThreeKalphaDeltaT =
 117             fvc::grad(threeKalpha*(T-T0), "grad(threeKalphaDeltaT)");
 118         surfaceVectorField threeKalphaDeltaTf =
 119             mesh.Sf()*threeKalphaf*fvc::interpolate(T-T0, "deltaT");
 120
 121         Info<< "Solving for " << U.name() << nl;
 122         do
 123         {
 124             U.storePrevIter();
 125
 126 #           include "calculateDivSigmaExp.H"
 127
 128             // Linear momentum equaiton
 129             fvVectorMatrix UEqn
 130             (
 131                 rho*fvm::d2dt2(U)
 132              ==
 133                 fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
 134               + divSigmaExp
 135             );
 136
 137             solverPerfU = UEqn.solve();
 138
 139             if (aitkenRelax)
 140             {
 141 #               include "aitkenRelaxation.H"
 142             }
 143             else
 144             {
 145                 U.relax();
 146             }
 147
 148             gradU = fvc::grad(U);
 149
 150 #           include "calculateRelResU.H"
 151
 152             if (iCorr == 0)
 153             {
 154                 initialResidual = solverPerfU.initialResidual();
 155             }
 156
 157             if (iCorr % infoFrequency == 0)
 158             {
 159                 Info<< "\tCorrector " << iCorr
 160                     << ", residual = " << solverPerfU.initialResidual()
 161                     << ", relative res = " << relResU;
 162
 163                 if (aitkenRelax)
 164                 {
 165                     Info << ", aitken = " << aitkenTheta;
 166                 }
 167                 Info<< ", inner iters = " << solverPerfU.nIterations() << endl;
 168             }
 169         }
 170         while
 171         (
 172             iCorr++ == 0
 173          || (
 174                 relResU > convergenceToleranceU
 175              && iCorr < nCorr
 176             )
 177         );
 178
 179         Info<< "Solved for " << U.name()
 180             << " using " << solverPerfU.solverName()
 181             << " in " << iCorr << " iterations"
 182             << ", initial res = " << initialResidual
 183             << ", final res = " << solverPerfU.initialResidual()
 184             << ", final rel res = " << relResU << nl
 185             << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 186             << ", ClockTime = " << runTime.elapsedClockTime() << " s"
 187             << endl;
 188
 189 #       include "calculateEpsilonSigma.H"
 190 #       include "writeFields.H"
 191
 192         Info<< "ExecutionTime = "
 193             << runTime.elapsedCpuTime()
 194             << " s\n\n" << endl;
 195     }
 196
 197     Info<< "End\n" << endl;
 198
 199     return(0);
 200 }
 201
 202
 203 // ************************************************************************* //