diffusionMulticomponent.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) 2015-2018 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
29#include "fvcGrad.H"
30#include "reactingMixture.H"
31#include "fvCFD.H"
33
34// * * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * //
35
36template<class ReactionThermo, class ThermoType>
37void Foam::combustionModels::
38diffusionMulticomponent<ReactionThermo, ThermoType>::init()
39{
40 // Load default values
41 this->coeffs().readIfPresent("Ci", Ci_);
42 this->coeffs().readIfPresent("YoxStream", YoxStream_);
43 this->coeffs().readIfPresent("YfStream", YfStream_);
44 this->coeffs().readIfPresent("sigma", sigma_);
45 this->coeffs().readIfPresent("ftCorr", ftCorr_);
46 this->coeffs().readIfPresent("alpha", alpha_);
47 this->coeffs().readIfPresent("laminarIgn", laminarIgn_);
48
49 typedef typename Reaction<ThermoType>::specieCoeffs specieCoeffs;
50
51 const speciesTable& species = this->thermo().composition().species();
52
53 scalarList specieStoichCoeffs(species.size());
54 const label nReactions = reactions_.size();
55
56 for (label k=0; k < nReactions; k++)
57 {
58 RijPtr_.set
59 (
60 k,
62 (
64 (
65 "Rijk" + Foam::name(k),
66 this->mesh_.time().timeName(),
67 this->mesh_,
70 false
71 ),
72 this->mesh_,
75 )
76 );
77
78 RijPtr_[k].storePrevIter();
79
80 const List<specieCoeffs>& lhs = reactions_[k].lhs();
81 const List<specieCoeffs>& rhs = reactions_[k].rhs();
82
83 const label fuelIndex = species.find(fuelNames_[k]);
84 const label oxidantIndex = species.find(oxidantNames_[k]);
85
86 const scalar Wu = specieThermo_[fuelIndex].W();
87 const scalar Wox = specieThermo_[oxidantIndex].W();
88
89 forAll(lhs, i)
90 {
91 const label specieI = lhs[i].index;
92 specieStoichCoeffs[specieI] = -lhs[i].stoichCoeff;
93 qFuel_[k] +=
94 specieThermo_[specieI].hc()*lhs[i].stoichCoeff/Wu;
95 }
96
97 forAll(rhs, i)
98 {
99 const label specieI = rhs[i].index;
100 specieStoichCoeffs[specieI] = rhs[i].stoichCoeff;
101 qFuel_[k] -=
102 specieThermo_[specieI].hc()*rhs[i].stoichCoeff/Wu;
103 }
104
105 Info << "Fuel heat of combustion : " << qFuel_[k] << endl;
106
107 s_[k] =
108 (Wox*mag(specieStoichCoeffs[oxidantIndex]))
109 / (Wu*mag(specieStoichCoeffs[fuelIndex]));
110
111 Info << "stoichiometric oxygen-fuel ratio : " << s_[k] << endl;
112
113 stoicRatio_[k] = s_[k]*YfStream_[k]/YoxStream_[k];
114
115 Info << "stoichiometric air-fuel ratio : " << stoicRatio_[k] << endl;
116
117 const scalar fStoich = 1.0/(1.0 + stoicRatio_[k]);
118
119 Info << "stoichiometric mixture fraction : " << fStoich << endl;
120 }
121}
122
123
124// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
125
126template<class ReactionThermo, class ThermoType>
129(
130 const word& modelType,
131 ReactionThermo& thermo,
133 const word& combustionProperties
134)
135:
136 ChemistryCombustion<ReactionThermo>
137 (
138 modelType,
139 thermo,
140 turb,
141 combustionProperties
142 ),
143 reactions_
144 (
145 dynamic_cast<const reactingMixture<ThermoType>&>(thermo)
146 ),
147 specieThermo_
148 (
149 dynamic_cast<const reactingMixture<ThermoType>&>(thermo).
150 speciesData()
151 ),
152 RijPtr_(reactions_.size()),
153 Ci_(reactions_.size(), 1.0),
154 fuelNames_(this->coeffs().lookup("fuels")),
155 oxidantNames_(this->coeffs().lookup("oxidants")),
156 qFuel_(reactions_.size()),
157 stoicRatio_(reactions_.size()),
158 s_(reactions_.size()),
159 YoxStream_(reactions_.size(), 0.23),
160 YfStream_(reactions_.size(), 1.0),
161 sigma_(reactions_.size(), 0.02),
162 oxidantRes_(this->coeffs().lookup("oxidantRes")),
163 ftCorr_(reactions_.size(), Zero),
164 alpha_(1),
165 laminarIgn_(false)
166{
167 init();
168}
169
170
171// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
172
173template<class ReactionThermo, class ThermoType>
174Foam::tmp<Foam::volScalarField> Foam::combustionModels::
175diffusionMulticomponent<ReactionThermo, ThermoType>::tc() const
176{
177 return this->chemistryPtr_->tc();
178}
179
180
181template<class ReactionThermo, class ThermoType>
184{
185 if (this->active())
186 {
187 typedef typename Reaction<ThermoType>::specieCoeffs specieCoeffs;
188 const speciesTable& species = this->thermo().composition().species();
189
190 const label nReactions = reactions_.size();
191
192 PtrList<volScalarField> RijlPtr(nReactions);
193
194 for (label k=0; k < nReactions; k++)
195 {
196 RijlPtr.set
197 (
198 k,
200 (
202 (
203 "Rijl" + Foam::name(k),
204 this->mesh_.time().timeName(),
205 this->mesh_,
208 false
209 ),
210 this->mesh_,
212 zeroGradientFvPatchScalarField::typeName
213 )
214 );
215
216 volScalarField& Rijl = RijlPtr[k];
217
218 // Obtain individual reaction rates for each reaction
219 const label fuelIndex = species.find(fuelNames_[k]);
220
221 if (laminarIgn_)
222 {
223 Rijl.ref() = -this->chemistryPtr_->calculateRR(k, fuelIndex);
224 }
225
226
227 // Look for the fuelStoic
228 const List<specieCoeffs>& rhs = reactions_[k].rhs();
229 const List<specieCoeffs>& lhs = reactions_[k].lhs();
230
231 // Set to zero RR's
232 forAll(lhs, l)
233 {
234 const label lIndex = lhs[l].index;
235 this->chemistryPtr_->RR(lIndex) =
237 }
238
239 forAll(rhs, l)
240 {
241 const label rIndex = rhs[l].index;
242 this->chemistryPtr_->RR(rIndex) =
244 }
245 }
246
247
248 for (label k=0; k < nReactions; k++)
249 {
250 const label fuelIndex = species.find(fuelNames_[k]);
251 const label oxidantIndex = species.find(oxidantNames_[k]);
252
253 const volScalarField& Yfuel =
254 this->thermo().composition().Y(fuelIndex);
255
256 const volScalarField& Yox =
257 this->thermo().composition().Y(oxidantIndex);
258
259 const volScalarField ft
260 (
261 "ft" + Foam::name(k),
262 (
263 s_[k]*Yfuel - (Yox - YoxStream_[k])
264 )
265 /(
266 s_[k]*YfStream_[k] + YoxStream_[k]
267 )
268 );
269
270 const scalar sigma = sigma_[k];
271
272 const scalar fStoich = 1.0/(1.0 + stoicRatio_[k]) + ftCorr_[k];
273
274 const volScalarField OAvailScaled
275 (
276 "OAvailScaled",
277 Yox/max(oxidantRes_[k], 1e-3)
278 );
279
280 const volScalarField preExp
281 (
282 "preExp" + Foam::name(k),
283 1.0 + sqr(OAvailScaled)
284 );
285
286 const volScalarField filter
287 (
288 (1.0/(sigma*sqrt(2.0*constant::mathematical::pi)))
289 * exp(-sqr(ft - fStoich)/(2*sqr(sigma)))
290 );
291
292 const volScalarField topHatFilter(pos(filter - 1e-3));
293
294 const volScalarField prob("prob" + Foam::name(k), preExp*filter);
295
296 const volScalarField RijkDiff
297 (
298 "RijkDiff",
299 Ci_[k]*this->turbulence().muEff()*prob*
300 (
301 mag(fvc::grad(Yfuel) & fvc::grad(Yox))
302 )
303 *pos(Yox)*pos(Yfuel)
304 );
305
306 volScalarField& Rijk = RijPtr_[k];
307
308 if (laminarIgn_)
309 {
310 Rijk =
311 min(RijkDiff, topHatFilter*RijlPtr[k]*pos(Yox)*pos(Yfuel));
312 }
313 else
314 {
315 Rijk = RijkDiff;
316 }
317
318 Rijk.relax(alpha_);
319
320 if (debug && this->mesh_.time().writeTime())
321 {
322 Rijk.write();
323 ft.write();
324 }
325
326 // Look for the fuelStoic
327 const List<specieCoeffs>& rhs = reactions_[k].rhs();
328 const List<specieCoeffs>& lhs = reactions_[k].lhs();
329
330 scalar fuelStoic = 1.0;
331 forAll(lhs, l)
332 {
333 const label lIndex = lhs[l].index;
334 if (lIndex == fuelIndex)
335 {
336 fuelStoic = lhs[l].stoichCoeff;
337 break;
338 }
339 }
340
341 const scalar MwFuel = specieThermo_[fuelIndex].W();
342
343 // Update left hand side species
344 forAll(lhs, l)
345 {
346 const label lIndex = lhs[l].index;
347
348 const scalar stoichCoeff = lhs[l].stoichCoeff;
349
350 this->chemistryPtr_->RR(lIndex) +=
351 -Rijk*stoichCoeff*specieThermo_[lIndex].W()/fuelStoic/MwFuel;
352
353 }
354
355 // Update right hand side species
356 forAll(rhs, r)
357 {
358 const label rIndex = rhs[r].index;
359
360 const scalar stoichCoeff = rhs[r].stoichCoeff;
361
362 this->chemistryPtr_->RR(rIndex) +=
363 Rijk*stoichCoeff*specieThermo_[rIndex].W()/fuelStoic/MwFuel;
364 }
365 }
366 }
367}
368
369
370template<class ReactionThermo, class ThermoType>
373(
375) const
376{
378
379 fvScalarMatrix& Su = tSu.ref();
380
381 if (this->active())
382 {
383 const label specieI =
384 this->thermo().composition().species().find(Y.member());
385
386 Su += this->chemistryPtr_->RR(specieI);
387 }
388
389 return tSu;
390}
391
392
393template<class ReactionThermo, class ThermoType>
397{
399 (
401 (
403 (
404 "Qdot",
405 this->mesh().time().timeName(),
406 this->mesh(),
409 false
410 ),
411 this->mesh(),
413 zeroGradientFvPatchScalarField::typeName
414 )
415 );
416
417 if (this->active())
418 {
419 volScalarField& Qdot = tQdot.ref();
420 Qdot = this->chemistryPtr_->Qdot();
421 }
422
423 return tQdot;
424}
425
426
427template<class ReactionThermo, class ThermoType>
430{
432 {
433 this->coeffs().readIfPresent("Ci", Ci_);
434 this->coeffs().readIfPresent("sigma", sigma_);
435 this->coeffs().readIfPresent("oxidantRes", oxidantRes_);
436 this->coeffs().readIfPresent("ftCorr", ftCorr_);
437 this->coeffs().readIfPresent("alpha", alpha_);
438 this->coeffs().readIfPresent("laminarIgn", laminarIgn_);
439 return true;
440 }
441
442 return false;
443}
444
445
446// ************************************************************************* //
label k
compressible::turbulenceModel & turb
Chemistry model wrapper for combustion models.
void relax(const scalar alpha)
Relax field (for steady-state solution).
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
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
const T * set(const label i) const
Definition: PtrList.H:138
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Diffusion based turbulent combustion model for multicomponent species.
virtual void correct()
Correct combustion rate.
virtual tmp< volScalarField > Qdot() const
Heat release rate calculated from fuel consumption rate matrix.
virtual bool read()
Update properties from given dictionary.
Abstract base class for turbulence models (RAS, LES and laminar).
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:121
A wordList with hashed named lookup, which can be faster in some situations than using the normal lis...
label find(const word &val) const
Find index of the value (searches the hash).
Lookup type of boundary radiation properties.
Definition: lookup.H:66
Foam::reactingMixture.
virtual bool write(const bool valid=true) const
Write using setting from DB.
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
T & ref() const
Definition: tmpI.H:227
A class for handling words, derived from Foam::string.
Definition: word.H:68
PtrList< volScalarField > & Y
dynamicFvMesh & mesh
compressible::turbulenceModel & turbulence
Calculate the gradient of the given field.
zeroField Su
Definition: alphaSuSp.H:1
word timeName
Definition: getTimeIndex.H:3
constexpr scalar pi(M_PI)
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensionedScalar pos(const dimensionedScalar &ds)
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)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimEnergy
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:45
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
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
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:60
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:51
scalar Qdot
Definition: solveChemistry.H:2
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
static const char *const typeName
The type name used in ensight case files.