All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rhoCentralFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 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 with support for mesh-motion and topology changes.
30 
31 \*---------------------------------------------------------------------------*/
32 
33 #include "fvCFD.H"
34 #include "dynamicFvMesh.H"
35 #include "psiThermo.H"
39 #include "directionInterpolate.H"
40 #include "localEulerDdtScheme.H"
41 #include "fvcSmooth.H"
42 
43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
44 
45 int main(int argc, char *argv[])
46 {
47  #define NO_CONTROL
48  #include "postProcess.H"
49 
50  #include "setRootCaseLists.H"
51  #include "createTime.H"
52  #include "createDynamicFvMesh.H"
53  #include "createFields.H"
54  #include "createFieldRefs.H"
55  #include "createTimeControls.H"
56 
57  turbulence->validate();
58 
59  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
60 
61  #include "readFluxScheme.H"
62 
63  dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
64 
65  // Courant numbers used to adjust the time-step
66  scalar CoNum = 0.0;
67  scalar meanCoNum = 0.0;
68 
69  Info<< "\nStarting time loop\n" << endl;
70 
71  while (runTime.run())
72  {
73  #include "readTimeControls.H"
74 
75  if (!LTS)
76  {
77  #include "setDeltaT.H"
78  runTime++;
79 
80  // Do any mesh changes
81  mesh.update();
82  }
83 
84  // --- Directed interpolation of primitive fields onto faces
85 
88 
89  surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
90  surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
91 
92  volScalarField rPsi("rPsi", 1.0/psi);
93  surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
94  surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
95 
96  surfaceScalarField e_pos(interpolate(e, pos, T.name()));
97  surfaceScalarField e_neg(interpolate(e, neg, T.name()));
98 
99  surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
100  surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
101 
102  surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
103  surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
104 
105  surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
106  surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
107 
108  // Make fluxes relative to mesh-motion
109  if (mesh.moving())
110  {
111  phiv_pos -= mesh.phi();
112  phiv_neg -= mesh.phi();
113  }
114 
115  volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
116  surfaceScalarField cSf_pos
117  (
118  "cSf_pos",
119  interpolate(c, pos, T.name())*mesh.magSf()
120  );
121  surfaceScalarField cSf_neg
122  (
123  "cSf_neg",
124  interpolate(c, neg, T.name())*mesh.magSf()
125  );
126 
128  (
129  "ap",
130  max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
131  );
133  (
134  "am",
135  min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
136  );
137 
138  surfaceScalarField a_pos("a_pos", ap/(ap - am));
139 
140  surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
141 
142  surfaceScalarField aSf("aSf", am*a_pos);
143 
144  if (fluxScheme == "Tadmor")
145  {
146  aSf = -0.5*amaxSf;
147  a_pos = 0.5;
148  }
149 
150  surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
151 
152  phiv_pos *= a_pos;
153  phiv_neg *= a_neg;
154 
155  surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
156  surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
157 
158  // Reuse amaxSf for the maximum positive and negative fluxes
159  // estimated by the central scheme
160  amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
161 
162  #include "centralCourantNo.H"
163 
164  if (LTS)
165  {
166  #include "setRDeltaT.H"
167  runTime++;
168  }
169 
170  Info<< "Time = " << runTime.timeName() << nl << endl;
171 
172  phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
173 
174  surfaceVectorField phiUp
175  (
176  (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
177  + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
178  );
179 
180  surfaceScalarField phiEp
181  (
182  "phiEp",
183  aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
184  + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
185  + aSf*p_pos - aSf*p_neg
186  );
187 
188  // Make flux for pressure-work absolute
189  if (mesh.moving())
190  {
191  phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg);
192  }
193 
194  volScalarField muEff("muEff", turbulence->muEff());
195  volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
196 
197  // --- Solve density
199 
200  // --- Solve momentum
201  solve(fvm::ddt(rhoU) + fvc::div(phiUp));
202 
203  U.ref() =
204  rhoU()
205  /rho();
206  U.correctBoundaryConditions();
207  rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();
208 
209  if (!inviscid)
210  {
211  solve
212  (
213  fvm::ddt(rho, U) - fvc::ddt(rho, U)
214  - fvm::laplacian(muEff, U)
215  - fvc::div(tauMC)
216  );
217  rhoU = rho*U;
218  }
219 
220  // --- Solve energy
221  surfaceScalarField sigmaDotU
222  (
223  "sigmaDotU",
224  (
225  fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
226  + fvc::dotInterpolate(mesh.Sf(), tauMC)
227  )
228  & (a_pos*U_pos + a_neg*U_neg)
229  );
230 
231  solve
232  (
233  fvm::ddt(rhoE)
234  + fvc::div(phiEp)
235  - fvc::div(sigmaDotU)
236  );
237 
238  e = rhoE/rho - 0.5*magSqr(U);
239  e.correctBoundaryConditions();
240  thermo.correct();
241  rhoE.boundaryFieldRef() ==
242  rho.boundaryField()*
243  (
244  e.boundaryField() + 0.5*magSqr(U.boundaryField())
245  );
246 
247  if (!inviscid)
248  {
249  solve
250  (
251  fvm::ddt(rho, e) - fvc::ddt(rho, e)
252  + thermophysicalTransport->divq(e)
253  );
254  thermo.correct();
255  rhoE = rho*(e + 0.5*magSqr(U));
256  }
257 
258  p.ref() =
259  rho()
260  /psi();
261  p.correctBoundaryConditions();
262  rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();
263 
264  turbulence->correct();
265  thermophysicalTransport->correct();
266 
267  runTime.write();
268 
269  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
270  << " ClockTime = " << runTime.elapsedClockTime() << " s"
271  << nl << endl;
272  }
273 
274  Info<< "End\n" << endl;
275 
276  return 0;
277 }
278 
279 // ************************************************************************* //
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
GeometricField< tensor, fvPatchField, volMesh > volTensorField
Definition: volFieldsFwd.H:62
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
engineTime & runTime
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Read the control parameters used by setDeltaT.
rhoEqn solve()
rhoReactionThermo & thermo
Definition: createFields.H:28
const dimensionedScalar & c
Speed of light in a vacuum.
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
dimensionedScalar neg(const dimensionedScalar &ds)
phi
Definition: pEqn.H:104
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:57
dimensionedScalar pos(const dimensionedScalar &ds)
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
dynamicFvMesh & mesh
rhoReactionThermophysicalTransportModel & thermophysicalTransport
const dimensionedScalar & e
Elementary charge.
Definition: doubleScalar.H:105
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Info<< "Reading field U\"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);volScalarField rho(IOobject("rho", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), thermo.rho());volVectorField rhoU(IOobject("rhoU", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *U);volScalarField rhoE(IOobject("rhoE", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *(e+0.5 *magSqr(U)));surfaceScalarField pos(IOobject("pos", runTime.timeName(), mesh), mesh, dimensionedScalar(dimless, 1.0));surfaceScalarField neg(IOobject("neg", runTime.timeName(), mesh), mesh, dimensionedScalar(dimless, -1.0));surfaceScalarField phi("phi", fvc::flux(rhoU));Info<< "Creating turbulence model\"<< endl;autoPtr< compressible::momentumTransportModel > turbulence(compressible::momentumTransportModel::New(rho, U, phi, thermo))
Definition: createFields.H:94
static const char nl
Definition: Ostream.H:260
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
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
dimensionedSymmTensor dev2(const dimensionedSymmTensor &dt)
U
Definition: pEqn.H:72
bool inviscid(true)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
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