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