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