FSD.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2019 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "FSD.H"
31#include "LESModel.H"
32#include "fvcGrad.H"
33#include "fvcDiv.H"
34
35namespace Foam
36{
37namespace combustionModels
38{
39
40// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
41
42template<class ReactionThermo, class ThermoType>
44(
45 const word& modelType,
46 ReactionThermo& thermo,
48 const word& combustionProperties
49)
50:
51 singleStepCombustion<ReactionThermo, ThermoType>
52 (
53 modelType,
54 thermo,
55 turb,
56 combustionProperties
57 ),
58 reactionRateFlameArea_
59 (
61 (
62 this->coeffs(),
63 this->mesh(),
64 *this
65 )
66 ),
67 ft_
68 (
70 (
71 this->thermo().phasePropertyName("ft"),
72 this->mesh().time().timeName(),
73 this->mesh(),
74 IOobject::NO_READ,
75 IOobject::AUTO_WRITE
76 ),
77 this->mesh(),
79 ),
80 YFuelFuelStream_(dimensionedScalar("YFuelStream", dimless, 1.0)),
81 YO2OxiStream_(dimensionedScalar("YOxiStream", dimless, 0.23)),
82 Cv_(this->coeffs().getScalar("Cv")),
83 C_(5.0),
84 ftMin_(0.0),
85 ftMax_(1.0),
86 ftDim_(300),
87 ftVarMin_(this->coeffs().getScalar("ftVarMin"))
88{}
89
90
91// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
92
93template<class ReactionThermo, class ThermoType>
95{}
96
97
98// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
99
100template<class ReactionThermo, class ThermoType>
102{
103 this->singleMixturePtr_->fresCorrect();
104
105 const label fuelI = this->singleMixturePtr_->fuelIndex();
106
107 const volScalarField& YFuel = this->thermo().composition().Y()[fuelI];
108
109 const volScalarField& YO2 = this->thermo().composition().Y("O2");
110
111 const dimensionedScalar s = this->singleMixturePtr_->s();
112
113 ft_ =
114 (s*YFuel - (YO2 - YO2OxiStream_))/(s*YFuelFuelStream_ + YO2OxiStream_);
115
116
117 volVectorField nft(fvc::grad(ft_));
118
119 volScalarField mgft(mag(nft));
120
121 surfaceVectorField SfHat(this->mesh().Sf()/this->mesh().magSf());
122
123 volScalarField cAux(scalar(1) - ft_);
124
125 dimensionedScalar dMgft = 1.0e-3*
126 (ft_*cAux*mgft)().weightedAverage(this->mesh().V())
127 /((ft_*cAux)().weightedAverage(this->mesh().V()) + SMALL)
128 + dimensionedScalar("ddMgft", mgft.dimensions(), SMALL);
129
130 mgft += dMgft;
131
132 nft /= mgft;
133
134 const volVectorField& U = YO2.db().lookupObject<volVectorField>("U");
135
136 const volScalarField sigma
137 (
138 (nft & nft)*fvc::div(U) - (nft & fvc::grad(U) & nft)
139 );
140
141 reactionRateFlameArea_->correct(sigma);
142
143 const volScalarField& omegaFuel = reactionRateFlameArea_->omega();
144
145
146 const scalar ftStoich =
147 YO2OxiStream_.value()
148 /(
149 s.value()*YFuelFuelStream_.value() + YO2OxiStream_.value()
150 );
151
153 (
155 (
157 (
158 this->thermo().phasePropertyName("Pc"),
159 U.time().timeName(),
160 U.db(),
163 ),
164 U.mesh(),
166 )
167 );
168
169 volScalarField& pc = tPc.ref();
170
171 tmp<volScalarField> tomegaFuel
172 (
174 (
176 (
177 this->thermo().phasePropertyName("omegaFuelBar"),
178 U.time().timeName(),
179 U.db(),
182 ),
183 U.mesh(),
184 dimensionedScalar(omegaFuel.dimensions(), Zero)
185 )
186 );
187
188 volScalarField& omegaFuelBar = tomegaFuel.ref();
189
190 // Calculation of the mixture fraction variance (ftVar)
191 const compressible::LESModel& lesModel =
193 (
195 );
196
197 const volScalarField& delta = lesModel.delta();
198 const volScalarField ftVar(Cv_*sqr(delta)*sqr(mgft));
199
200 // Thickened flame (average flame thickness for counterflow configuration
201 // is 1.5 mm)
202
203 volScalarField deltaF
204 (
205 delta/dimensionedScalar("flame", dimLength, 1.5e-3)
206 );
207
208 // Linear correlation between delta and flame thickness
209 volScalarField omegaF(max(deltaF*(4.0/3.0) + (2.0/3.0), scalar(1)));
210
211 scalar deltaFt = 1.0/ftDim_;
212
213 forAll(ft_, celli)
214 {
215 if (ft_[celli] > ftMin_ && ft_[celli] < ftMax_)
216 {
217 scalar ftCell = ft_[celli];
218
219 if (ftVar[celli] > ftVarMin_) //sub-grid beta pdf of ft_
220 {
221 scalar ftVarc = ftVar[celli];
222 scalar a =
223 max(ftCell*(ftCell*(1.0 - ftCell)/ftVarc - 1.0), 0.0);
224 scalar b = max(a/ftCell - a, 0.0);
225
226 for (int i=1; i<ftDim_; i++)
227 {
228 scalar ft = i*deltaFt;
229 pc[celli] += pow(ft, a-1.0)*pow(1.0 - ft, b - 1.0)*deltaFt;
230 }
231
232 for (int i=1; i<ftDim_; i++)
233 {
234 scalar ft = i*deltaFt;
235 omegaFuelBar[celli] +=
236 omegaFuel[celli]/omegaF[celli]
237 *exp
238 (
239 -sqr(ft - ftStoich)
240 /(2.0*sqr(0.01*omegaF[celli]))
241 )
242 *pow(ft, a - 1.0)
243 *pow(1.0 - ft, b - 1.0)
244 *deltaFt;
245 }
246 omegaFuelBar[celli] /= max(pc[celli], 1e-4);
247 }
248 else
249 {
250 omegaFuelBar[celli] =
251 omegaFuel[celli]/omegaF[celli]
252 *exp(-sqr(ftCell - ftStoich)/(2.0*sqr(0.01*omegaF[celli])));
253 }
254 }
255 else
256 {
257 omegaFuelBar[celli] = 0.0;
258 }
259 }
260
261
262 // Combustion progress variable, c
263
264 List<label> productsIndex(2, label(-1));
265 {
266 label i = 0;
267 forAll(this->singleMixturePtr_->specieProd(), specieI)
268 {
269 if (this->singleMixturePtr_->specieProd()[specieI] < 0)
270 {
271 productsIndex[i] = specieI;
272 i++;
273 }
274 }
275 }
276
277
278 // Flamelet probability of the progress c based on IFC (reuse pc)
279 scalar YprodTotal = 0;
280 forAll(productsIndex, j)
281 {
282 YprodTotal += this->singleMixturePtr_->Yprod0()[productsIndex[j]];
283 }
284
285 forAll(ft_, celli)
286 {
287 if (ft_[celli] < ftStoich)
288 {
289 pc[celli] = ft_[celli]*(YprodTotal/ftStoich);
290 }
291 else
292 {
293 pc[celli] = (1.0 - ft_[celli])*(YprodTotal/(1.0 - ftStoich));
294 }
295 }
296
297 tmp<volScalarField> tproducts
298 (
300 (
301 IOobject
302 (
303 this->thermo().phasePropertyName("products"),
304 U.time().timeName(),
305 U.db(),
308 ),
309 U.mesh(),
311 )
312 );
313
314 volScalarField& products = tproducts.ref();
315
316 forAll(productsIndex, j)
317 {
318 label specieI = productsIndex[j];
319 const volScalarField& Yp = this->thermo().composition().Y()[specieI];
320 products += Yp;
321 }
322
324 (
325 max(scalar(1) - products/max(pc, scalar(1e-5)), scalar(0))
326 );
327
328 pc = min(C_*c, scalar(1));
329
330 const volScalarField fres(this->singleMixturePtr_->fres(fuelI));
331
332 this->wFuel_ == mgft*pc*omegaFuelBar;
333}
334
335
336template<class ReactionThermo, class ThermoType>
338{
340
341 if (this->active())
342 {
343 calculateSourceNorm();
344 }
345}
346
347
348template<class ReactionThermo, class ThermoType>
350{
352 {
353 this->coeffs().readEntry("Cv", Cv_);
354 this->coeffs().readEntry("ftVarMin", ftVarMin_);
355 reactionRateFlameArea_->read(this->coeffs());
356 return true;
357 }
358
359 return false;
360}
361
362// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
363
364} // End namespace combustionModels
365} // End namespace Foam
366
367// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scalar delta
Macros for easy insertion into run-time selection tables.
compressible::turbulenceModel & turb
const dimensionSet & dimensions() const
Return dimensions.
Internal & ref(const bool updateAccessTime=true)
Return a reference to the dimensioned internal field.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const objectRegistry & db() const noexcept
Return the local objectRegistry.
Definition: IOobject.C:500
Flame Surface Dennsity (FDS) combustion model.
Definition: FSD.H:88
virtual void correct()
Correct combustion rate.
Definition: FSD.C:337
virtual ~FSD()
Destructor.
Definition: FSD.C:94
virtual bool read()
Update properties.
Definition: FSD.C:349
Base class for combustion models using singleStepReactingMixture.
Abstract base class for turbulence models (RAS, LES and laminar).
const Type & lookupObject(const word &name, const bool recursive=false) const
Abstract class for reaction rate per flame area unit.
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
A class for managing temporary objects.
Definition: tmp.H:65
static const word propertiesName
Default name of the turbulence properties dictionary.
A class for handling words, derived from Foam::string.
Definition: word.H:68
U
Definition: pEqn.H:72
dynamicFvMesh & mesh
Calculate the divergence of the given field.
Calculate the gradient of the given field.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
word timeName
Definition: getTimeIndex.H:3
LESModel< EddyDiffusivity< turbulenceModel > > LESModel
Typedefs for turbulence, RAS and LES models for compressible flow based on the standard laminar trans...
const dimensionedScalar c
Speed of light in a vacuum.
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:49
Namespace for OpenFOAM.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedScalar exp(const dimensionedScalar &ds)
const dimensionSet dimless
Dimensionless.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:60
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:51
volScalarField & b
Definition: createFields.H:27
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333