PDRFoamAutoRefine.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-2018 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  PDRFoam
26 
27 Description
28  Solver for compressible premixed/partially-premixed combustion with
29  turbulence modelling.
30 
31  Combusting RANS code using the b-Xi two-equation model.
32  Xi may be obtained by either the solution of the Xi transport
33  equation or from an algebraic exression. Both approaches are
34  based on Gulder's flame speed correlation which has been shown
35  to be appropriate by comparison with the results from the
36  spectral model.
37 
38  Strain effects are incorporated directly into the Xi equation
39  but not in the algebraic approximation. Further work need to be
40  done on this issue, particularly regarding the enhanced removal rate
41  caused by flame compression. Analysis using results of the spectral
42  model will be required.
43 
44  For cases involving very lean Propane flames or other flames which are
45  very strain-sensitive, a transport equation for the laminar flame
46  speed is present. This equation is derived using heuristic arguments
47  involving the strain time scale and the strain-rate at extinction.
48  the transport velocity is the same as that for the Xi equation.
49 
50  For large flames e.g. explosions additional modelling for the flame
51  wrinkling due to surface instabilities may be applied.
52 
53  PDR (porosity/distributed resistance) modelling is included to handle
54  regions containing blockages which cannot be resolved by the mesh.
55 
56 \*---------------------------------------------------------------------------*/
57 
58 #include "fvCFD.H"
59 #include "dynamicFvMesh.H"
60 #include "psiuReactionThermo.H"
62 #include "laminarFlameSpeed.H"
63 #include "XiModel.H"
64 #include "PDRDragModel.H"
65 #include "ignition.H"
66 #include "Switch.H"
67 #include "bound.H"
68 #include "dynamicRefineFvMesh.H"
69 #include "pimpleControl.H"
70 
71 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
72 
73 int main(int argc, char *argv[])
74 {
75  #include "setRootCaseLists.H"
76 
77  #include "createTime.H"
78  #include "createDynamicFvMesh.H"
79 
80  pimpleControl pimple(mesh);
81 
82  #include "readCombustionProperties.H"
83  #include "readGravitationalAcceleration.H"
84  #include "createFields.H"
85  #include "initContinuityErrs.H"
86  #include "createTimeControls.H"
87  #include "compressibleCourantNo.H"
88  #include "setInitialDeltaT.H"
89 
90  turbulence->validate();
91  scalar StCoNum = 0.0;
92 
93  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
94 
95  Info<< "\nStarting time loop\n" << endl;
96 
97  bool hasChanged = false;
98 
99  while (runTime.run())
100  {
101  #include "readTimeControls.H"
102  #include "compressibleCourantNo.H"
103  #include "setDeltaT.H"
104 
105  // Indicators for refinement. Note: before runTime++
106  // only for postprocessing reasons.
107  tmp<volScalarField> tmagGradP = mag(fvc::grad(p));
108  volScalarField normalisedGradP
109  (
110  "normalisedGradP",
111  tmagGradP()/max(tmagGradP())
112  );
113  normalisedGradP.writeOpt() = IOobject::AUTO_WRITE;
114  tmagGradP.clear();
115 
116  runTime++;
117 
118  Info<< "\n\nTime = " << runTime.timeName() << endl;
119 
120  {
121  // Make the fluxes absolute
123 
124  // Test : disable refinement for some cells
125  PackedBoolList& protectedCell =
126  refCast<dynamicRefineFvMesh>(mesh).protectedCell();
127 
128  if (protectedCell.empty())
129  {
130  protectedCell.setSize(mesh.nCells());
131  protectedCell = 0;
132  }
133 
134  forAll(betav, celli)
135  {
136  if (betav[celli] < 0.99)
137  {
138  protectedCell[celli] = 1;
139  }
140  }
141 
142  // Flux estimate for introduced faces.
143  volVectorField rhoU("rhoU", rho*U);
144 
145  // Do any mesh changes
146  bool meshChanged = mesh.update();
147 
148 
149  if (meshChanged)
150  {
151  hasChanged = true;
152  }
153 
154  if (runTime.write() && hasChanged)
155  {
156  betav.write();
157  Lobs.write();
158  CT.write();
159  drag->writeFields();
160  flameWrinkling->writeFields();
161  hasChanged = false;
162  }
163 
164  // Make the fluxes relative to the mesh motion
166  }
167 
168 
169  #include "rhoEqn.H"
170 
171  // --- Pressure-velocity PIMPLE corrector loop
172  while (pimple.loop())
173  {
174  #include "UEqn.H"
175 
176 
177  // --- Pressure corrector loop
178  while (pimple.correct())
179  {
180  #include "bEqn.H"
181  #include "ftEqn.H"
182  #include "huEqn.H"
183  #include "hEqn.H"
184 
185  if (!ign.ignited())
186  {
187  hu == h;
188  }
189 
190  #include "pEqn.H"
191  }
192 
193  if (pimple.turbCorr())
194  {
195  turbulence->correct();
196  }
197  }
198 
199  runTime.write();
200 
201  Info<< "\nExecutionTime = "
202  << runTime.elapsedCpuTime()
203  << " s\n" << endl;
204  }
205 
206  Info<< "\n end\n";
207 
208  return 0;
209 }
210 
211 
212 // ************************************************************************* //
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
pimpleNoLoopControl & pimple
surfaceScalarField & phi
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
engineTime & runTime
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
Read the control parameters used by setDeltaT.
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
Info<< "Reading strained laminar flame speed field Su\"<< endl;volScalarField Su(IOobject("Su", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Reading field betav\"<< endl;volScalarField betav(IOobject("betav", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);Info<< "Reading field Lobs\"<< endl;volScalarField Lobs(IOobject("Lobs", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);Info<< "Reading field CT\"<< endl;volSymmTensorField CT(IOobject("CT", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);Info<< "Reading field Nv\"<< endl;volScalarField Nv(IOobject("Nv", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);Info<< "Reading field nsv\"<< endl;volSymmTensorField nsv(IOobject("nsv", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);IOdictionary PDRProperties(IOobject("PDRProperties", runTime.constant(), mesh, IOobject::MUST_READ_IF_MODIFIED, IOobject::NO_WRITE));autoPtr< PDRDragModel > drag
Definition: createFields.H:182
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
dynamicFvMesh & mesh
const volScalarField & betav
Bound the given scalar field if it has gone unbounded.
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("pos", dimless, 1.0));surfaceScalarField neg(IOobject("neg", runTime.timeName(), mesh), mesh, dimensionedScalar("neg", dimless, -1.0));surfaceScalarField phi("phi", fvc::flux(rhoU));Info<< "Creating turbulence model\"<< endl;autoPtr< compressible::turbulenceModel > turbulence(compressible::turbulenceModel::New(rho, U, phi, thermo))
Definition: createFields.H:94
U
Definition: pEqn.H:72
volScalarField & h
Planck constant.
void makeAbsolute(surfaceScalarField &phi, const volVectorField &U)
Make the given flux absolute.
Definition: fvcMeshPhi.C:114
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
volScalarField & p
void makeRelative(surfaceScalarField &phi, const volVectorField &U)
Make the given flux relative.
Definition: fvcMeshPhi.C:75
Read the control parameters used by setDeltaT.
autoPtr< XiModel > flameWrinkling
Create the flame-wrinkling model.
Definition: createFields.H:192