1 /*---------------------------------------------------------------------------*\
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 -------------------------------------------------------------------------------
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;
67 mesh.solutionDict().subDict("PIMPLE").lookup("HbyAblend")
70 # include "readTimeControls.H"
74 mesh.surfaceInterpolation::deltaCoeffs()
75 *mag(phiv)/mesh.magSf()
76 ).value()*runTime.deltaT().value();
78 Info<< "Max Courant Number = " << CoNum << endl;
80 # include "setDeltaT.H"
87 fv::multivariateGaussConvectionScheme<scalar> mvConvection
92 mesh.schemesDict().divScheme("div(phiv,rhoUH)")
98 + mvConvection.fvmDiv(phiv, rho)
101 surfaceScalarField rhoUWeights =
102 mvConvection.interpolationScheme()()(magRhoU)()
105 weighted<vector> rhoUScheme(rhoUWeights);
107 fvVectorMatrix rhoUEqn
110 + fv::gaussConvectionScheme<vector>(mesh, phiv, rhoUScheme)
114 solve(rhoUEqn == -fvc::grad(p));
119 + mvConvection.fvmDiv(phiv, rhoE)
121 - mvConvection.fvcDiv(phiv, p)
124 T = (rhoE - 0.5*rho*magSqr(rhoU/rho))/Cv/rho;
128 while (pimple.correct())
130 volScalarField rrhoUA = 1.0/rhoUEqn.A();
131 surfaceScalarField rrhoUAf("rrhoUAf", fvc::interpolate(rrhoUA));
132 volVectorField HbyA = rrhoUA*rhoUEqn.H();
134 surfaceScalarField HbyAWeights =
135 HbyAblend*mesh.weights()
138 <vector, MUSCLLimiter<NVDTVD>, limitFuncs::magSqr>
139 (mesh, phi, IStringStream("HbyA")()).weights(HbyA);
143 surfaceInterpolationScheme<vector>::interpolate
144 (HbyA, HbyAWeights) & mesh.Sf()
146 + HbyAblend*fvc::ddtPhiCorr(rrhoUA, rho, rhoU, phi);
148 p.boundaryField().updateCoeffs();
150 surfaceScalarField phiGradp =
151 rrhoUAf*mesh.magSf()*fvc::snGrad(p);
155 # include "resetPhiPatches.H"
157 surfaceScalarField rhof =
158 mvConvection.interpolationScheme()()(rho)()
166 + mvConvection.fvcDiv(phiv, rho)
168 - fvm::laplacian(rrhoUAf, p)
173 phi += phiGradp + pEqn.flux();
176 mvConvection.interpolationScheme()()(rho)()
180 rhoU = HbyA - rrhoUA*fvc::grad(p);
181 rhoU.correctBoundaryConditions();
189 Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
190 << " ClockTime = " << runTime.elapsedClockTime() << " s"
194 Info<< "End\n" << endl;
200 // ************************************************************************* //