applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.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     sonicDyMFoam
  26
  27 Description
  28     Transient solver for trans-sonic/supersonic for laminar or turbulent
  29     flow of a compressible gas with support for mesh motion and
  30     topological changes.
  31
  32     Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
  33     pseudo-transient simulations.  The pressure-energy coupling is done
  34     using the Rusche manoeuvre (isentropic compression/expansion).
  35
  36     Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
  37
  38     Updated from sonicFoamAutoMotion by Hrvoje Jasak
  39
  40 Author
  41     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  42
  43 \*---------------------------------------------------------------------------*/
  44
  45 #include "fvCFD.H"
  46 #include "dynamicFvMesh.H"
  47 #include "specie.H"
  48 #include "basicPsiThermo.H"
  49 #include "turbulenceModel.H"
  50
  51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  52
  53 int main(int argc, char *argv[])
  54 {
  55 #   include "setRootCase.H"
  56
  57 #   include "createTime.H"
  58 #   include "createDynamicFvMesh.H"
  59 #   include "createFields.H"
  60 #   include "initContinuityErrs.H"
  61
  62 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  63
  64     Info<< "\nStarting time loop\n" << endl;
  65
  66     while (runTime.run())
  67     {
  68 #       include "readControls.H"
  69 #       include "readFieldBounds.H"
  70 #       include "compressibleCourantNo.H"
  71 #       include "setDeltaT.H"
  72
  73         runTime++;
  74         Info<< "deltaT = " << runTime.deltaT().value() << nl << endl;
  75         Info<< "Time = " << runTime.timeName() << nl << endl;
  76
  77         bool meshChanged = mesh.update();
  78
  79 #       include "volContinuity.H"
  80
  81         if (checkMeshCourantNo)
  82         {
  83 #           include "meshCourantNo.H"
  84         }
  85
  86         // Mesh motion update
  87         if (meshChanged)
  88         {
  89             T.max(TMin);
  90             p.max(pMin);
  91             e == max(e, thermo.Cv()*TMin);
  92
  93             thermo.correct();
  94             rho = thermo.rho();
  95
  96             if (correctPhi)
  97             {
  98 // #           include "correctPhi.H"
  99             }
 100         }
 101
 102         if (meshChanged)
 103         {
 104 #           include "compressibleCourantNo.H"
 105         }
 106
 107         // --- PIMPLE loop
 108         label oCorr = 0;
 109         do
 110         {
 111 #           include "rhoEqn.H"
 112 #           include "eEqn.H"
 113 #           include "UEqn.H"
 114
 115             // --- PISO loop
 116             volScalarField rUA = 1.0/UEqn.A();
 117
 118             surfaceScalarField psisf = fvc::interpolate(psis);
 119             surfaceScalarField rhof = fvc::interpolate(rho);
 120
 121             // Needs to be outside of loop since p is changing,
 122             // but psi and rho are not
 123             surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
 124
 125             for (int corr = 0; corr < nCorr; corr++)
 126             {
 127 #               include "pEqn.H"
 128             }
 129
 130             turbulence->correct();
 131         } while (++oCorr < nOuterCorr);
 132
 133         runTime.write();
 134
 135         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 136             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 137             << nl << endl;
 138     }
 139
 140     Info<< "End\n" << endl;
 141
 142     return(0);
 143 }
 144
 145
 146 // ************************************************************************* //