ENH: autoLayerDriver: better layering information message
[OpenFOAM-2.0.x.git] / src / thermophysicalModels / reactionThermo / mixtures / singleStepReactingMixture / singleStepReactingMixture.C
blob10a317f1804ea8da5a99726904d7e0979a4b5e7e
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
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.
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.
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/>.
24 \*---------------------------------------------------------------------------*/
26 #include "singleStepReactingMixture.H"
27 #include "fvMesh.H"
29 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
31 template<class ThermoType>
32 void Foam::singleStepReactingMixture<ThermoType>::calculateqFuel()
34     const Reaction<ThermoType>& reaction = this->operator[](0);
35     const  scalar Wu = this->speciesData()[fuelIndex_].W();
37     forAll(reaction.lhs(), i)
38     {
39         const label specieI = reaction.lhs()[i].index;
40         const scalar stoichCoeff = reaction.lhs()[i].stoichCoeff;
41         specieStoichCoeffs_[specieI] = -stoichCoeff;
42         qFuel_.value() += this->speciesData()[specieI].hc()*stoichCoeff/Wu;
43     }
45     forAll(reaction.rhs(), i)
46     {
47         const label specieI = reaction.rhs()[i].index;
48         const scalar stoichCoeff = reaction.rhs()[i].stoichCoeff;
49         specieStoichCoeffs_[specieI] = stoichCoeff;
50         qFuel_.value() -= this->speciesData()[specieI].hc()*stoichCoeff/Wu;
51         specieProd_[specieI] = -1;
52     }
54     Info << "Fuel heat of combustion :" << qFuel_.value() << endl;
58 template<class ThermoType>
59 void Foam::singleStepReactingMixture<ThermoType>::massAndAirStoichRatios()
61     const label O2Index = this->species()["O2"];
62     const scalar Wu = this->speciesData()[fuelIndex_].W();
64     stoicRatio_ =
65        (this->speciesData()[inertIndex_].W()
66       * specieStoichCoeffs_[inertIndex_]
67       + this->speciesData()[O2Index].W()
68       * mag(specieStoichCoeffs_[O2Index]))
69       / (Wu*mag(specieStoichCoeffs_[fuelIndex_]));
71     s_ =
72         (this->speciesData()[O2Index].W()
73       * mag(specieStoichCoeffs_[O2Index]))
74       / (Wu*mag(specieStoichCoeffs_[fuelIndex_]));
76     Info << "stoichiometric air-fuel ratio :" << stoicRatio_.value() << endl;
78     Info << "stoichiometric oxygen-fuel ratio :" << s_.value() << endl;
82 template<class ThermoType>
83 void Foam::singleStepReactingMixture<ThermoType>::calculateMaxProducts()
85     const Reaction<ThermoType>& reaction = this->operator[](0);
87     scalar Wm = 0.0;
88     scalar totalMol = 0.0;
89     forAll(reaction.rhs(), i)
90     {
91         label specieI = reaction.rhs()[i].index;
92         totalMol += mag(specieStoichCoeffs_[specieI]);
93     }
95     scalarList Xi(reaction.rhs().size());
97     forAll(reaction.rhs(), i)
98     {
99         const label specieI = reaction.rhs()[i].index;
100         Xi[i] = mag(specieStoichCoeffs_[specieI])/totalMol;
102         Wm += Xi[i]*this->speciesData()[specieI].W();
103     }
105     forAll(reaction.rhs(), i)
106     {
107         const label specieI = reaction.rhs()[i].index;
108         Yprod0_[specieI] =  this->speciesData()[specieI].W()/Wm*Xi[i];
109     }
111     // Normalize the stoichiometric coeff to mass
112     forAll(specieStoichCoeffs_, i)
113     {
114         specieStoichCoeffs_[i] =
115             specieStoichCoeffs_[i]
116           * this->speciesData()[i].W()
117           / (this->speciesData()[fuelIndex_].W()
118           * mag(specieStoichCoeffs_[fuelIndex_]));
119     }
123 template<class ThermoType>
124 void Foam::singleStepReactingMixture<ThermoType>::fresCorrect()
126     const Reaction<ThermoType>& reaction = this->operator[](0);
128     label O2Index = this->species()["O2"];
129     const volScalarField& YFuel = this->Y()[fuelIndex_];
130     const volScalarField& YO2 = this->Y()[O2Index];
132     // reactants
133     forAll(reaction.lhs(), i)
134     {
135         const label specieI = reaction.lhs()[i].index;
136         if (specieI == fuelIndex_)
137         {
138             fres_[specieI] =  max(YFuel - YO2/s_, scalar(0.0));
139         }
140         else if (specieI == O2Index)
141         {
142             fres_[specieI] =  max(YO2 - YFuel*s_, scalar(0.0));
143         }
144     }
147     // products
148     forAll(reaction.rhs(), i)
149     {
150         const label specieI = reaction.rhs()[i].index;
151         if (specieI != inertIndex_)
152         {
153             forAll(fres_[specieI], cellI)
154             {
155                 if (fres_[fuelIndex_][cellI] > 0.0)
156                 {
157                     // rich mixture
158                     fres_[specieI][cellI] =
159                         Yprod0_[specieI]
160                       * (1.0 + YO2[cellI]/s_.value() - YFuel[cellI]);
161                 }
162                 else
163                 {
164                     // lean mixture
165                     fres_[specieI][cellI] =
166                         Yprod0_[specieI]
167                       * (
168                             1.0
169                           - YO2[cellI]/s_.value()*stoicRatio_.value()
170                           + YFuel[cellI]*stoicRatio_.value()
171                         );
172                 }
173             }
174         }
175     }
179 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
181 template<class ThermoType>
182 Foam::singleStepReactingMixture<ThermoType>::singleStepReactingMixture
184     const dictionary& thermoDict,
185     const fvMesh& mesh
188     reactingMixture<ThermoType>(thermoDict, mesh),
189     stoicRatio_(dimensionedScalar("stoicRatio", dimless, 0.0)),
190     s_(dimensionedScalar("s", dimless, 0.0)),
191     qFuel_(dimensionedScalar("qFuel", sqr(dimVelocity), 0.0)),
192     specieStoichCoeffs_(this->species_.size(), 0.0),
193     Yprod0_(this->species_.size(), 0.0),
194     fres_(Yprod0_.size()),
195     inertIndex_(this->species()[thermoDict.lookup("inertSpecie")]),
196     fuelIndex_(this->species()[thermoDict.lookup("fuel")]),
197     specieProd_(Yprod0_.size(), 1)
199     if (this->size() == 1)
200     {
201         forAll(fres_, fresI)
202         {
203             IOobject header
204             (
205                 "fres_" + this->species()[fresI],
206                 mesh.time().timeName(),
207                 mesh,
208                 IOobject::NO_READ,
209                 IOobject::NO_WRITE
210             );
212             fres_.set
213             (
214                 fresI,
215                 new volScalarField
216                 (
217                     header,
218                     mesh,
219                     dimensionedScalar("fres" + name(fresI), dimless, 0.0)
220                 )
221             );
222         }
224         calculateqFuel();
226         massAndAirStoichRatios();
228         calculateMaxProducts();
230         autoPtr<chemistryReader<ThermoType> >::clear();
231     }
232     else
233     {
234         FatalErrorIn
235         (
236             "singleStepReactingMixture::<ThermoType>::singleStepReactingMixture"
237             "("
238                 "const dictionary&, "
239                 "const fvMesh&"
240             ")"
241         )   << "Only one reaction required for single step reaction"
242             << exit(FatalError);
243     }
247 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
249 template<class ThermoType>
250 void Foam::singleStepReactingMixture<ThermoType>::read
252     const dictionary& thermoDict
257 // ************************************************************************* //