rhoCentralFoam.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 rhoCentralFoam
29
30Group
31 grpCompressibleSolvers
32
33Description
34 Density-based compressible flow solver based on
35 central-upwind schemes of Kurganov and Tadmor with
36 support for mesh-motion and topology changes.
37
38\*---------------------------------------------------------------------------*/
39
40#include "fvCFD.H"
41#include "dynamicFvMesh.H"
42#include "psiThermo.H"
46#include "localEulerDdtScheme.H"
47#include "fvcSmooth.H"
48
49// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50
51int main(int argc, char *argv[])
52{
53 argList::addNote
54 (
55 "Density-based compressible flow solver based on"
56 " central-upwind schemes of Kurganov and Tadmor with"
57 " support for mesh-motion and topology changes."
58 );
59
60 #define NO_CONTROL
61 #include "postProcess.H"
62
63 #include "addCheckCaseOptions.H"
64 #include "setRootCaseLists.H"
65 #include "createTime.H"
66 #include "createDynamicFvMesh.H"
67 #include "createFields.H"
68 #include "createFieldRefs.H"
69 #include "createTimeControls.H"
70
71 turbulence->validate();
72
73 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
74
75 #include "readFluxScheme.H"
76
77 const dimensionedScalar v_zero(dimVolume/dimTime, Zero);
78
79 // Courant numbers used to adjust the time-step
80 scalar CoNum = 0.0;
81 scalar meanCoNum = 0.0;
82
83 Info<< "\nStarting time loop\n" << endl;
84
85 while (runTime.run())
86 {
87 #include "readTimeControls.H"
88
89 if (!LTS)
90 {
91 #include "setDeltaT.H"
92
93 ++runTime;
94
95 // Do any mesh changes
96 mesh.update();
97 }
98
99 // --- Directed interpolation of primitive fields onto faces
100
101 surfaceScalarField rho_pos(interpolate(rho, pos));
102 surfaceScalarField rho_neg(interpolate(rho, neg));
103
104 surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
105 surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
106
107 volScalarField rPsi("rPsi", 1.0/psi);
108 surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
109 surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
110
111 surfaceScalarField e_pos(interpolate(e, pos, T.name()));
112 surfaceScalarField e_neg(interpolate(e, neg, T.name()));
113
114 surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
115 surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
116
117 surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
118 surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
119
120 surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
121 // Note: extracted out the orientation so becomes unoriented
122 phiv_pos.setOriented(false);
123 surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
124 phiv_neg.setOriented(false);
125
126 // Make fluxes relative to mesh-motion
127 if (mesh.moving())
128 {
130 meshPhi.setOriented(false);
131 phiv_pos -= meshPhi;
132 phiv_neg -= meshPhi;
133 }
134
135 volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
136 surfaceScalarField cSf_pos
137 (
138 "cSf_pos",
139 interpolate(c, pos, T.name())*mesh.magSf()
140 );
141
142 surfaceScalarField cSf_neg
143 (
144 "cSf_neg",
145 interpolate(c, neg, T.name())*mesh.magSf()
146 );
147
149 (
150 "ap",
151 max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
152 );
153
155 (
156 "am",
157 min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
158 );
159
160 surfaceScalarField a_pos("a_pos", ap/(ap - am));
161
162 surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
163
164 surfaceScalarField aSf("aSf", am*a_pos);
165
166 if (fluxScheme == "Tadmor")
167 {
168 aSf = -0.5*amaxSf;
169 a_pos = 0.5;
170 }
171
172 surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
173
174 phiv_pos *= a_pos;
175 phiv_neg *= a_neg;
176
177 surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
178 surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
179
180 // Reuse amaxSf for the maximum positive and negative fluxes
181 // estimated by the central scheme
182 amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
183
184 #include "centralCourantNo.H"
185
186 if (LTS)
187 {
188 #include "setRDeltaT.H"
189
190 ++runTime;
191 }
192
193 Info<< "Time = " << runTime.timeName() << nl << endl;
194
195 phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
196
197 surfaceVectorField phiU(aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg);
198 // Note: reassembled orientation from the pos and neg parts so becomes
199 // oriented
200 phiU.setOriented(true);
201
202 surfaceVectorField phiUp(phiU + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf());
203
205 (
206 "phiEp",
207 aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
208 + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
209 + aSf*p_pos - aSf*p_neg
210 );
211
212 // Make flux for pressure-work absolute
213 if (mesh.moving())
214 {
216 meshPhi.setOriented(false);
217 phiEp += meshPhi*(a_pos*p_pos + a_neg*p_neg);
218 }
219
220 volScalarField muEff("muEff", turbulence->muEff());
221 volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
222
223 // --- Solve density
224 solve(fvm::ddt(rho) + fvc::div(phi));
225
226 // --- Solve momentum
227 solve(fvm::ddt(rhoU) + fvc::div(phiUp));
228
229 U.ref() =
230 rhoU()
231 /rho();
232 U.correctBoundaryConditions();
233 rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();
234
235 if (!inviscid)
236 {
237 solve
238 (
239 fvm::ddt(rho, U) - fvc::ddt(rho, U)
240 - fvm::laplacian(muEff, U)
241 - fvc::div(tauMC)
242 );
243 rhoU = rho*U;
244 }
245
246 // --- Solve energy
247 surfaceScalarField sigmaDotU
248 (
249 "sigmaDotU",
250 (
251 fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
252 + fvc::dotInterpolate(mesh.Sf(), tauMC)
253 )
254 & (a_pos*U_pos + a_neg*U_neg)
255 );
256
257 solve
258 (
259 fvm::ddt(rhoE)
260 + fvc::div(phiEp)
261 - fvc::div(sigmaDotU)
262 );
263
264 e = rhoE/rho - 0.5*magSqr(U);
265 e.correctBoundaryConditions();
266 thermo.correct();
267 rhoE.boundaryFieldRef() ==
268 rho.boundaryField()*
269 (
270 e.boundaryField() + 0.5*magSqr(U.boundaryField())
271 );
272
273 if (!inviscid)
274 {
275 solve
276 (
277 fvm::ddt(rho, e) - fvc::ddt(rho, e)
278 - fvm::laplacian(turbulence->alphaEff(), e)
279 );
280 thermo.correct();
281 rhoE = rho*(e + 0.5*magSqr(U));
282 }
283
284 p.ref() =
285 rho()
286 /psi();
287 p.correctBoundaryConditions();
288 rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();
289
290 turbulence->correct();
291
292 runTime.write();
293
294 runTime.printExecutionTime(Info);
295 }
296
297 Info<< "End\n" << endl;
298
299 return 0;
300}
301
302// ************************************************************************* //
Y[inertIndex] max(0.0)
Required Classes.
Calculates the mean and maximum wave speed based Courant Numbers.
surfaceScalarField & phi
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
U
Definition: pEqn.H:72
volScalarField & p
const volScalarField & psi
const volScalarField & T
bool inviscid(true)
mesh interpolate(rAU)
dynamicFvMesh & mesh
engineTime & runTime
bool LTS
Definition: createRDeltaT.H:1
Read the control parameters used by setDeltaT.
compressible::turbulenceModel & turbulence
Provides functions smooth spread and sweep which use the FaceCellWave algorithm to smooth and redistr...
const dimensionedScalar c
Speed of light in a vacuum.
tmp< surfaceScalarField > meshPhi(const volVectorField &U)
Definition: fvcMeshPhi.C:36
dimensionedSymmTensor dev2(const dimensionedSymmTensor &dt)
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
messageStream Info
Information stream (stdout output on master, null elsewhere)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
GeometricField< tensor, fvPatchField, volMesh > volTensorField
Definition: volFieldsFwd.H:87
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Execute application functionObjects to post-process existing results.
word fluxScheme("Kurganov")
volScalarField & rhoE
Read the control parameters used by setDeltaT.
CEqn solve()
volScalarField & e
Definition: createFields.H:11
scalar CoNum
scalar meanCoNum