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 | Copyright (C) 2011-2016 OpenFOAM Foundation
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 Application
25  rhoCentralFoam
26 
27 Description
28  Density-based compressible flow solver based on central-upwind schemes of
29  Kurganov and Tadmor.
30 
31 \*---------------------------------------------------------------------------*/
32 
33 #include "fvCFD.H"
34 #include "psiThermo.H"
37 #include "directionInterpolate.H"
38 #include "localEulerDdtScheme.H"
39 #include "fvcSmooth.H"
40 
41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 
43 int main(int argc, char *argv[])
44 {
45  #define NO_CONTROL
46  #include "postProcess.H"
47 
48  #include "setRootCase.H"
49  #include "createTime.H"
50  #include "createMesh.H"
51  #include "createFields.H"
52  #include "createFieldRefs.H"
53  #include "createTimeControls.H"
54  #include "createRDeltaT.H"
55 
56  turbulence->validate();
57 
58  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
59 
60  #include "readFluxScheme.H"
61 
62  dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
63 
64  // Courant numbers used to adjust the time-step
65  scalar CoNum = 0.0;
66  scalar meanCoNum = 0.0;
67 
68  Info<< "\nStarting time loop\n" << endl;
69 
70  while (runTime.run())
71  {
72  // --- Directed interpolation of primitive fields onto faces
73 
76 
77  surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
78  surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
79 
80  volScalarField rPsi("rPsi", 1.0/psi);
81  surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
82  surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
83 
84  surfaceScalarField e_pos(interpolate(e, pos, T.name()));
85  surfaceScalarField e_neg(interpolate(e, neg, T.name()));
86 
87  surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
88  surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
89 
90  surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
91  surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
92 
93  surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
94  surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
95 
96  volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
97  surfaceScalarField cSf_pos
98  (
99  "cSf_pos",
100  interpolate(c, pos, T.name())*mesh.magSf()
101  );
102  surfaceScalarField cSf_neg
103  (
104  "cSf_neg",
105  interpolate(c, neg, T.name())*mesh.magSf()
106  );
107 
109  (
110  "ap",
111  max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
112  );
114  (
115  "am",
116  min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
117  );
118 
119  surfaceScalarField a_pos("a_pos", ap/(ap - am));
120 
121  surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
122 
123  surfaceScalarField aSf("aSf", am*a_pos);
124 
125  if (fluxScheme == "Tadmor")
126  {
127  aSf = -0.5*amaxSf;
128  a_pos = 0.5;
129  }
130 
131  surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
132 
133  phiv_pos *= a_pos;
134  phiv_neg *= a_neg;
135 
136  surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
137  surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
138 
139  // Reuse amaxSf for the maximum positive and negative fluxes
140  // estimated by the central scheme
141  amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
142 
143  #include "centralCourantNo.H"
144  #include "readTimeControls.H"
145 
146  if (LTS)
147  {
148  #include "setRDeltaT.H"
149  }
150  else
151  {
152  #include "setDeltaT.H"
153  }
154 
155  runTime++;
156 
157  Info<< "Time = " << runTime.timeName() << nl << endl;
158 
159  phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
160 
161  surfaceVectorField phiUp
162  (
163  (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
164  + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
165  );
166 
167  surfaceScalarField phiEp
168  (
169  "phiEp",
170  aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
171  + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
172  + aSf*p_pos - aSf*p_neg
173  );
174 
175  volScalarField muEff("muEff", turbulence->muEff());
176  volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
177 
178  // --- Solve density
180 
181  // --- Solve momentum
182  solve(fvm::ddt(rhoU) + fvc::div(phiUp));
183 
184  U.ref() =
185  rhoU()
186  /rho();
187  U.correctBoundaryConditions();
188  rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();
189 
190  if (!inviscid)
191  {
192  solve
193  (
194  fvm::ddt(rho, U) - fvc::ddt(rho, U)
195  - fvm::laplacian(muEff, U)
196  - fvc::div(tauMC)
197  );
198  rhoU = rho*U;
199  }
200 
201  // --- Solve energy
202  surfaceScalarField sigmaDotU
203  (
204  "sigmaDotU",
205  (
206  fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
207  + fvc::dotInterpolate(mesh.Sf(), tauMC)
208  )
209  & (a_pos*U_pos + a_neg*U_neg)
210  );
211 
212  solve
213  (
214  fvm::ddt(rhoE)
215  + fvc::div(phiEp)
216  - fvc::div(sigmaDotU)
217  );
218 
219  e = rhoE/rho - 0.5*magSqr(U);
220  e.correctBoundaryConditions();
221  thermo.correct();
222  rhoE.boundaryFieldRef() ==
223  rho.boundaryField()*
224  (
225  e.boundaryField() + 0.5*magSqr(U.boundaryField())
226  );
227 
228  if (!inviscid)
229  {
230  solve
231  (
232  fvm::ddt(rho, e) - fvc::ddt(rho, e)
233  - fvm::laplacian(turbulence->alphaEff(), e)
234  );
235  thermo.correct();
236  rhoE = rho*(e + 0.5*magSqr(U));
237  }
238 
239  p.ref() =
240  rho()
241  /psi();
242  p.correctBoundaryConditions();
243  rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();
244 
245  turbulence->correct();
246 
247  runTime.write();
248 
249  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
250  << " ClockTime = " << runTime.elapsedClockTime() << " s"
251  << nl << endl;
252  }
253 
254  Info<< "End\n" << endl;
255 
256  return 0;
257 }
258 
259 // ************************************************************************* //
autoPtr< compressible::turbulenceModel > turbulence
Definition: createFields.H:23
surfaceScalarField & phi
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
U
Definition: pEqn.H:83
volScalarField & rhoE
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
GeometricField< tensor, fvPatchField, volMesh > volTensorField
Definition: volFieldsFwd.H:59
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
static tmp< GeometricField< typename innerProduct< vector, Type >::type, fvsPatchField, surfaceMesh > > dotInterpolate(const surfaceVectorField &Sf, const GeometricField< Type, fvPatchField, volMesh > &tvf)
Interpolate field onto faces.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Read the control parameters used by setDeltaT.
tmp< surfaceScalarField > interpolate(const RhoType &rho)
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
dimensionedScalar neg(const dimensionedScalar &ds)
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
psiReactionThermo & thermo
Definition: createFields.H:31
dimensionedScalar pos(const dimensionedScalar &ds)
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
dynamicFvMesh & mesh
rhoEqn solve()
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static const char nl
Definition: Ostream.H:262
bool LTS
Definition: createRDeltaT.H:1
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const volScalarField & T
dimensionedSymmTensor dev2(const dimensionedSymmTensor &dt)
bool inviscid(true)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const dimensionedScalar c
Speed of light in a vacuum.
Calculates the mean and maximum wave speed based Courant Numbers.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const volScalarField & psi
volScalarField & p
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Execute application functionObjects to post-process existing results.
word fluxScheme("Kurganov")
Provides functions smooth spread and sweep which use the FaceCellWave algorithm to smooth and redistr...
Read the control parameters used by setDeltaT.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45