applications/solvers/compressible/rhopSonicFoam/rhopSonicFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Application
  26     rhopSonicFoam
  27
  28 Description
  29     Pressure-density-based compressible flow solver.
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "fvCFD.H"
  34 #include "weighted.H"
  35 #include "gaussConvectionScheme.H"
  36 #include "multivariateGaussConvectionScheme.H"
  37 #include "MUSCL.H"
  38 #include "LimitedScheme.H"
  39 #include "boundaryTypes.H"
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 int main(int argc, char *argv[])
  44 {
  45     #include "setRootCase.H"
  46     #include "createTime.H"
  47     #include "createMesh.H"
  48     #include "readThermodynamicProperties.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.value() << nl << endl;
  58
  59         #include "readPISOControls.H"
  60         scalar HbyAblend = readScalar(piso.lookup("HbyAblend"));
  61
  62         #include "readTimeControls.H"
  63
  64         scalar CoNum = max
  65         (
  66             mesh.surfaceInterpolation::deltaCoeffs()
  67            *mag(phiv)/mesh.magSf()
  68         ).value()*runTime.deltaT().value();
  69
  70         Info<< "Max Courant Number = " << CoNum << endl;
  71
  72         #include "setDeltaT.H"
  73
  74         for (int outerCorr=0; outerCorr<nOuterCorr; outerCorr++)
  75         {
  76             magRhoU = mag(rhoU);
  77             H = (rhoE + p)/rho;
  78
  79             fv::multivariateGaussConvectionScheme<scalar> mvConvection
  80             (
  81                 mesh,
  82                 fields,
  83                 phiv,
  84                 mesh.divScheme("div(phiv,rhoUH)")
  85             );
  86
  87             solve
  88             (
  89                 fvm::ddt(rho)
  90               + mvConvection.fvmDiv(phiv, rho)
  91             );
  92
  93             surfaceScalarField rhoUWeights =
  94                 mvConvection.interpolationScheme()()(magRhoU)()
  95                .weights(magRhoU);
  96
  97             weighted<vector> rhoUScheme(rhoUWeights);
  98
  99             fvVectorMatrix rhoUEqn
 100             (
 101                 fvm::ddt(rhoU)
 102               + fv::gaussConvectionScheme<vector>(mesh, phiv, rhoUScheme)
 103                    .fvmDiv(phiv, rhoU)
 104             );
 105
 106             solve(rhoUEqn == -fvc::grad(p));
 107
 108             solve
 109             (
 110                 fvm::ddt(rhoE)
 111               + mvConvection.fvmDiv(phiv, rhoE)
 112              ==
 113               - mvConvection.fvcDiv(phiv, p)
 114             );
 115
 116             T = (rhoE - 0.5*rho*magSqr(rhoU/rho))/Cv/rho;
 117             psi = 1.0/(R*T);
 118             p = rho/psi;
 119
 120             for (int corr=0; corr<nCorr; corr++)
 121             {
 122                 volScalarField rrhoUA = 1.0/rhoUEqn.A();
 123                 surfaceScalarField rrhoUAf("rrhoUAf", fvc::interpolate(rrhoUA));
 124                 volVectorField HbyA = rrhoUA*rhoUEqn.H();
 125
 126                 surfaceScalarField HbyAWeights =
 127                     HbyAblend*mesh.weights()
 128                   + (1.0 - HbyAblend)*
 129                     LimitedScheme
 130                         <vector, MUSCLLimiter<NVDTVD>, limitFuncs::magSqr>
 131                         (mesh, phi, IStringStream("HbyA")()).weights(HbyA);
 132
 133                 phi =
 134                     (
 135                         surfaceInterpolationScheme<vector>::interpolate
 136                         (HbyA, HbyAWeights) & mesh.Sf()
 137                     )
 138                   + HbyAblend*fvc::ddtPhiCorr(rrhoUA, rho, rhoU, phi);
 139
 140                 p.boundaryField().updateCoeffs();
 141
 142                 surfaceScalarField phiGradp =
 143                     rrhoUAf*mesh.magSf()*fvc::snGrad(p);
 144
 145                 phi -= phiGradp;
 146
 147                 #include "resetPhiPatches.H"
 148
 149                 surfaceScalarField rhof =
 150                     mvConvection.interpolationScheme()()(rho)()
 151                    .interpolate(rho);
 152
 153                 phiv = phi/rhof;
 154
 155                 fvScalarMatrix pEqn
 156                 (
 157                     fvm::ddt(psi, p)
 158                   + mvConvection.fvcDiv(phiv, rho)
 159                   + fvc::div(phiGradp)
 160                   - fvm::laplacian(rrhoUAf, p)
 161                 );
 162
 163                 pEqn.solve();
 164
 165                 phi += phiGradp + pEqn.flux();
 166                 rho = psi*p;
 167                 rhof =
 168                     mvConvection.interpolationScheme()()(rho)()
 169                    .interpolate(rho);
 170                 phiv = phi/rhof;
 171
 172                 rhoU = HbyA - rrhoUA*fvc::grad(p);
 173                 rhoU.correctBoundaryConditions();
 174             }
 175         }
 176
 177         U = rhoU/rho;
 178
 179         runTime.write();
 180
 181         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 182             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 183             << nl << endl;
 184     }
 185
 186     Info<< "End\n" << endl;
 187
 188     return 0;
 189 }
 190
 191
 192 // ************************************************************************* //