FSD.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 \*---------------------------------------------------------------------------*/
25 
26 #include "FSD.H"
28 #include "LESModel.H"
29 #include "fvcGrad.H"
30 #include "fvcDiv.H"
32 
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 namespace combustionModels
39 {
40  defineTypeNameAndDebug(FSD, 0);
41  addToRunTimeSelectionTable(combustionModel, FSD, dictionary);
42 }
43 }
44 
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
49 (
50  const word& modelType,
51  const fluidReactionThermo& thermo,
53  const word& combustionProperties
54 )
55 :
57  (
58  modelType,
59  thermo,
60  turb,
61  combustionProperties
62  ),
63  reactionRateFlameArea_
64  (
66  (
67  this->coeffs(),
68  this->mesh(),
69  *this
70  )
71  ),
72  ft_
73  (
74  IOobject
75  (
76  this->thermo().phasePropertyName("ft"),
77  this->mesh().time().timeName(),
78  this->mesh(),
81  ),
82  this->mesh(),
84  ),
85  YFuelFuelStream_(dimensionedScalar(dimless, 1.0)),
86  YO2OxiStream_(dimensionedScalar(dimless, 0.23)),
87  Cv_(this->coeffs().template lookup<scalar>("Cv")),
88  C_(5.0),
89  ftMin_(0.0),
90  ftMax_(1.0),
91  ftDim_(300),
92  ftVarMin_(this->coeffs().template lookup<scalar>("ftVarMin"))
93 {}
94 
95 
96 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
97 
99 {}
100 
101 
102 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
103 
104 void Foam::combustionModels::FSD::calculateSourceNorm()
105 {
106  this->fresCorrect();
107 
108  const label fuelI = this->fuelIndex();
109 
110  const volScalarField& YFuel = this->thermo().composition().Y()[fuelI];
111 
112  const volScalarField& YO2 = this->thermo().composition().Y("O2");
113 
114  const dimensionedScalar s = this->s();
115 
116  ft_ =
117  (s*YFuel - (YO2 - YO2OxiStream_))/(s*YFuelFuelStream_ + YO2OxiStream_);
118 
119 
120  volVectorField nft(fvc::grad(ft_));
121 
122  volScalarField mgft(mag(nft));
123 
124  surfaceVectorField SfHat(this->mesh().Sf()/this->mesh().magSf());
125 
126  volScalarField cAux(scalar(1) - ft_);
127 
128  dimensionedScalar dMgft = 1.0e-3*
129  (ft_*cAux*mgft)().weightedAverage(this->mesh().V())
130  /((ft_*cAux)().weightedAverage(this->mesh().V()) + small)
131  + dimensionedScalar(mgft.dimensions(), small);
132 
133  mgft += dMgft;
134 
135  nft /= mgft;
136 
137  const volVectorField& U = YO2.db().lookupObject<volVectorField>("U");
138 
139  const volScalarField sigma
140  (
141  (nft & nft)*fvc::div(U) - (nft & fvc::grad(U) & nft)
142  );
143 
144  reactionRateFlameArea_->correct(sigma);
145 
146  const volScalarField& omegaFuel = reactionRateFlameArea_->omega();
147 
148 
149  const scalar ftStoich =
150  YO2OxiStream_.value()
151  /(
152  s.value()*YFuelFuelStream_.value() + YO2OxiStream_.value()
153  );
154 
156  (
158  (
159  this->thermo().phasePropertyName("Pc"),
160  U.mesh(),
162  )
163  );
164 
165  volScalarField& pc = tPc.ref();
166 
167  tmp<volScalarField> tomegaFuel
168  (
170  (
171  this->thermo().phasePropertyName("omegaFuelBar"),
172  U.mesh(),
173  dimensionedScalar(omegaFuel.dimensions(), 0)
174  )
175  );
176 
177  volScalarField& omegaFuelBar = tomegaFuel.ref();
178 
179  // Calculation of the mixture fraction variance (ftVar)
180  const compressible::LESModel& lesModel =
182  (
183  momentumTransportModel::typeName
184  );
185 
186  const volScalarField& delta = lesModel.delta();
187  const volScalarField ftVar(Cv_*sqr(delta)*sqr(mgft));
188 
189  // Thickened flame (average flame thickness for counterflow configuration
190  // is 1.5 mm)
191 
192  volScalarField deltaF
193  (
194  lesModel.delta()/dimensionedScalar(dimLength, 1.5e-3)
195  );
196 
197  // Linear correlation between delta and flame thickness
198  volScalarField omegaF(max(deltaF*(4.0/3.0) + (2.0/3.0), scalar(1)));
199 
200  scalar deltaFt = 1.0/ftDim_;
201 
202  forAll(ft_, celli)
203  {
204  if (ft_[celli] > ftMin_ && ft_[celli] < ftMax_)
205  {
206  scalar ftCell = ft_[celli];
207 
208  if (ftVar[celli] > ftVarMin_) // sub-grid beta pdf of ft_
209  {
210  scalar ftVarc = ftVar[celli];
211  scalar a =
212  max(ftCell*(ftCell*(1.0 - ftCell)/ftVarc - 1.0), 0.0);
213  scalar b = max(a/ftCell - a, 0.0);
214 
215  for (int i=1; i<ftDim_; i++)
216  {
217  scalar ft = i*deltaFt;
218  pc[celli] += pow(ft, a-1.0)*pow(1.0 - ft, b - 1.0)*deltaFt;
219  }
220 
221  for (int i=1; i<ftDim_; i++)
222  {
223  scalar ft = i*deltaFt;
224  omegaFuelBar[celli] +=
225  omegaFuel[celli]/omegaF[celli]
226  *exp
227  (
228  -sqr(ft - ftStoich)
229  /(2.0*sqr(0.01*omegaF[celli]))
230  )
231  *pow(ft, a - 1.0)
232  *pow(1.0 - ft, b - 1.0)
233  *deltaFt;
234  }
235  omegaFuelBar[celli] /= max(pc[celli], 1e-4);
236  }
237  else
238  {
239  omegaFuelBar[celli] =
240  omegaFuel[celli]/omegaF[celli]
241  *exp(-sqr(ftCell - ftStoich)/(2.0*sqr(0.01*omegaF[celli])));
242  }
243  }
244  else
245  {
246  omegaFuelBar[celli] = 0.0;
247  }
248  }
249 
250 
251  // Combustion progress variable, c
252 
253  List<label> productsIndex(2, label(-1));
254  {
255  label i = 0;
256  forAll(this->specieProd(), specieI)
257  {
258  if (this->specieProd()[specieI] < 0)
259  {
260  productsIndex[i] = specieI;
261  i++;
262  }
263  }
264  }
265 
266 
267  // Flamelet probability of the progress c based on IFC (reuse pc)
268  scalar YprodTotal = 0;
269  forAll(productsIndex, j)
270  {
271  YprodTotal += this->Yprod0()[productsIndex[j]];
272  }
273 
274  forAll(ft_, celli)
275  {
276  if (ft_[celli] < ftStoich)
277  {
278  pc[celli] = ft_[celli]*(YprodTotal/ftStoich);
279  }
280  else
281  {
282  pc[celli] = (1.0 - ft_[celli])*(YprodTotal/(1.0 - ftStoich));
283  }
284  }
285 
286  tmp<volScalarField> tproducts
287  (
289  (
290  this->thermo().phasePropertyName("products"),
291  U.mesh(),
293  )
294  );
295 
296  volScalarField& products = tproducts.ref();
297 
298  forAll(productsIndex, j)
299  {
300  label specieI = productsIndex[j];
301  const volScalarField& Yp = this->thermo().composition().Y()[specieI];
302  products += Yp;
303  }
304 
306  (
307  max(scalar(1) - products/max(pc, scalar(1e-5)), scalar(0))
308  );
309 
310  pc = min(C_*c, scalar(1));
311 
312  const volScalarField fres(this->fres(fuelI));
313 
314  this->wFuel_ == mgft*pc*omegaFuelBar;
315 }
316 
317 
319 {
320  this->wFuel_ ==
322 
323  calculateSourceNorm();
324 }
325 
326 
328 {
330  {
331  this->coeffs().lookup("Cv") >> Cv_ ;
332  this->coeffs().lookup("ftVarMin") >> ftVarMin_;
333  reactionRateFlameArea_->read(this->coeffs());
334  return true;
335  }
336  else
337  {
338  return false;
339  }
340 }
341 
342 
343 // ************************************************************************* //
scalar delta
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:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
fluidReactionThermo & thermo
Definition: createFields.H:28
virtual bool read()
Update properties.
Definition: FSD.C:327
U
Definition: pEqn.H:72
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedSymmTensor sqr(const dimensionedVector &dv)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
Base-class for multi-component fluid thermodynamic properties.
const dimensionSet dimless
LESModel< momentumTransportModel > LESModel
fvMesh & mesh
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
const dimensionedScalar c
Speed of light in a vacuum.
Macros for easy insertion into run-time selection tables.
const dimensionSet dimLength
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
const dimensionSet dimTime
const dimensionSet & dimensions() const
Return dimensions.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
FSD(const word &modelType, const fluidReactionThermo &thermo, const compressibleMomentumTransportModel &turb, const word &combustionProperties)
Construct from components.
Definition: FSD.C:49
dimensionedScalar exp(const dimensionedScalar &ds)
Calculate the gradient of the given field.
A class for handling words, derived from string.
Definition: word.H:59
const Type & value() const
Return const reference to value.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
word timeName
Definition: getTimeIndex.H:3
Calculate the divergence of the given field.
const Mesh & mesh() const
Return mesh.
const dimensionSet dimMass
Internal & ref()
Return a reference to the dimensioned internal field.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
virtual bool read()
Update properties from given dictionary.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
virtual void correct()
Correct combustion rate.
Definition: FSD.C:318
dimensioned< scalar > mag(const dimensioned< Type > &)
Base class for combustion models using basicSpecieMixture.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
A class for managing temporary objects.
Definition: PtrList.H:53
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:312
defineTypeNameAndDebug(diffusion, 0)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:98
static autoPtr< reactionRateFlameArea > New(const dictionary &dict, const fvMesh &mesh, const combustionModel &combModel)
Base class for single-phase compressible turbulence models.
addToRunTimeSelectionTable(combustionModel, diffusion, dictionary)
virtual ~FSD()
Destructor.
Definition: FSD.C:98
Namespace for OpenFOAM.