1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | foam-extend: Open Source CFD
4 \\ / O peration | Version: 4.0
5 \\ / A nd | Web: http://www.foam-extend.org
6 \\/ M anipulation | For copyright notice see file Copyright
7 -------------------------------------------------------------------------------
9 This file is part of foam-extend.
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.
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.
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/>.
28 Pressure-density-based compressible flow solver.
30 \*---------------------------------------------------------------------------*/
34 #include "gaussConvectionScheme.H"
35 #include "multivariateGaussConvectionScheme.H"
37 #include "LimitedScheme.H"
38 #include "boundaryTypes.H"
39 #include "pimpleControl.H"
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43 int main(int argc, char *argv[])
46 # include "setRootCase.H"
47 # include "createTime.H"
48 # include "createMesh.H"
50 pimpleControl pimple(mesh);
52 # include "readThermodynamicProperties.H"
53 # include "createFields.H"
54 # include "createTimeControls.H"
56 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
58 Info<< "\nStarting time loop\n" << endl;
60 while (runTime.loop())
62 Info<< "Time = " << runTime.value() << nl << endl;
64 scalar HbyAblend = readScalar(pimple.dict().lookup("HbyAblend"));
66 # include "readTimeControls.H"
70 mesh.surfaceInterpolation::deltaCoeffs()
71 *mag(phiv)/mesh.magSf()
72 ).value()*runTime.deltaT().value();
74 Info<< "Max Courant Number = " << CoNum << endl;
76 # include "setDeltaT.H"
83 fv::multivariateGaussConvectionScheme<scalar> mvConvection
88 mesh.schemesDict().divScheme("div(phiv,rhoUH)")
94 + mvConvection.fvmDiv(phiv, rho)
97 surfaceScalarField rhoUWeights =
98 mvConvection.interpolationScheme()()(magRhoU)()
101 weighted<vector> rhoUScheme(rhoUWeights);
103 fvVectorMatrix rhoUEqn
106 + fv::gaussConvectionScheme<vector>(mesh, phiv, rhoUScheme)
110 solve(rhoUEqn == -fvc::grad(p));
115 + mvConvection.fvmDiv(phiv, rhoE)
117 - mvConvection.fvcDiv(phiv, p)
120 T = (rhoE - 0.5*rho*magSqr(rhoU/rho))/Cv/rho;
124 while (pimple.correct())
126 volScalarField rrhoUA = 1.0/rhoUEqn.A();
127 surfaceScalarField rrhoUAf("rrhoUAf", fvc::interpolate(rrhoUA));
128 volVectorField HbyA = rrhoUA*rhoUEqn.H();
130 surfaceScalarField HbyAWeights =
131 HbyAblend*mesh.weights()
134 <vector, MUSCLLimiter<NVDTVD>, limitFuncs::magSqr>
135 (mesh, phi, IStringStream("HbyA")()).weights(HbyA);
139 surfaceInterpolationScheme<vector>::interpolate
140 (HbyA, HbyAWeights) & mesh.Sf()
142 + HbyAblend*fvc::ddtPhiCorr(rrhoUA, rho, rhoU, phi);
144 surfaceScalarField phiGradp =
145 rrhoUAf*mesh.magSf()*fvc::snGrad(p);
149 # include "resetPhiPatches.H"
151 surfaceScalarField rhof =
152 mvConvection.interpolationScheme()()(rho)()
160 + mvConvection.fvcDiv(phiv, rho)
162 - fvm::laplacian(rrhoUAf, p)
167 phi += phiGradp + pEqn.flux();
170 mvConvection.interpolationScheme()()(rho)()
174 rhoU = HbyA - rrhoUA*fvc::grad(p);
175 rhoU.correctBoundaryConditions();
183 Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
184 << " ClockTime = " << runTime.elapsedClockTime() << " s"
188 Info<< "End\n" << endl;
194 // ************************************************************************* //