rhoCentralDyMFoam.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  rhoCentralDyMFoam
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 "motionSolver.H"
41 
42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43 
44 int main(int argc, char *argv[])
45 {
46  #include "setRootCase.H"
47  #include "createTime.H"
48  #include "createDynamicFvMesh.H"
49  #include "createFields.H"
50  #include "createTimeControls.H"
51 
52  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
53 
54  #include "readFluxScheme.H"
55 
56  dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
57 
58  // Courant numbers used to adjust the time-step
59  scalar CoNum = 0.0;
60  scalar meanCoNum = 0.0;
61 
62  Info<< "\nStarting time loop\n" << endl;
63 
64  while (runTime.run())
65  {
66  #include "readTimeControls.H"
67  #include "setDeltaT.H"
68 
69  runTime++;
70 
71  Info<< "Time = " << runTime.timeName() << nl << endl;
72 
73  // Do any mesh changes
74  mesh.update();
75 
76  // --- Directed interpolation of primitive fields onto faces
77 
80 
81  surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
82  surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
83 
84  volScalarField rPsi("rPsi", 1.0/psi);
85  surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
86  surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
87 
88  surfaceScalarField e_pos(interpolate(e, pos, T.name()));
89  surfaceScalarField e_neg(interpolate(e, neg, T.name()));
90 
91  surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
92  surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
93 
94  surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
95  surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
96 
97  surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
98  surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
99 
100  // Make fluxes relative to mesh-motion
101  if (mesh.moving())
102  {
103  phiv_pos -= mesh.phi();
104  phiv_neg -= mesh.phi();
105  }
106 
107  volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
108  surfaceScalarField cSf_pos
109  (
110  "cSf_pos",
111  interpolate(c, pos, T.name())*mesh.magSf()
112  );
113  surfaceScalarField cSf_neg
114  (
115  "cSf_neg",
116  interpolate(c, neg, T.name())*mesh.magSf()
117  );
118 
120  (
121  "ap",
122  max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
123  );
125  (
126  "am",
127  min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
128  );
129 
130  surfaceScalarField a_pos("a_pos", ap/(ap - am));
131 
132  surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
133 
134  surfaceScalarField aSf("aSf", am*a_pos);
135 
136  if (fluxScheme == "Tadmor")
137  {
138  aSf = -0.5*amaxSf;
139  a_pos = 0.5;
140  }
141 
142  surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
143 
144  phiv_pos *= a_pos;
145  phiv_neg *= a_neg;
146 
147  surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
148  surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
149 
150  // Reuse amaxSf for the maximum positive and negative fluxes
151  // estimated by the central scheme
152  amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
153 
154  #include "centralCourantNo.H"
155 
156  phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
157 
158  surfaceVectorField phiUp
159  (
160  (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
161  + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
162  );
163 
164  surfaceScalarField phiEp
165  (
166  "phiEp",
167  aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
168  + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
169  + aSf*p_pos - aSf*p_neg
170  );
171 
172  // Make flux for pressure-work absolute
173  if (mesh.moving())
174  {
175  phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg);
176  }
177 
178  volScalarField muEff("muEff", turbulence->muEff());
179  volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
180 
181  // --- Solve density
183 
184  // --- Solve momentum
185  solve(fvm::ddt(rhoU) + fvc::div(phiUp));
186 
187  U.dimensionedInternalField() =
188  rhoU.dimensionedInternalField()
189  /rho.dimensionedInternalField();
190  U.correctBoundaryConditions();
191  rhoU.boundaryField() == rho.boundaryField()*U.boundaryField();
192 
193  if (!inviscid)
194  {
195  solve
196  (
197  fvm::ddt(rho, U) - fvc::ddt(rho, U)
198  - fvm::laplacian(muEff, U)
199  - fvc::div(tauMC)
200  );
201  rhoU = rho*U;
202  }
203 
204  // --- Solve energy
205  surfaceScalarField sigmaDotU
206  (
207  "sigmaDotU",
208  (
209  fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
210  + (mesh.Sf() & fvc::interpolate(tauMC))
211  )
212  & (a_pos*U_pos + a_neg*U_neg)
213  );
214 
215  solve
216  (
217  fvm::ddt(rhoE)
218  + fvc::div(phiEp)
219  - fvc::div(sigmaDotU)
220  );
221 
222  e = rhoE/rho - 0.5*magSqr(U);
223  e.correctBoundaryConditions();
224  thermo.correct();
225  rhoE.boundaryField() ==
226  rho.boundaryField()*
227  (
228  e.boundaryField() + 0.5*magSqr(U.boundaryField())
229  );
230 
231  if (!inviscid)
232  {
233  solve
234  (
235  fvm::ddt(rho, e) - fvc::ddt(rho, e)
236  - fvm::laplacian(turbulence->alphaEff(), e)
237  );
238  thermo.correct();
239  rhoE = rho*(e + 0.5*magSqr(U));
240  }
241 
242  p.dimensionedInternalField() =
243  rho.dimensionedInternalField()
244  /psi.dimensionedInternalField();
245  p.correctBoundaryConditions();
246  rho.boundaryField() == psi.boundaryField()*p.boundaryField();
247 
248  turbulence->correct();
249 
250  runTime.write();
251 
252  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
253  << " ClockTime = " << runTime.elapsedClockTime() << " s"
254  << nl << endl;
255  }
256 
257  Info<< "End\n" << endl;
258 
259  return 0;
260 }
261 
262 // ************************************************************************* //
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
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.
Read the control parameters used by setDeltaT.
U
Definition: pEqn.H:82