EulerImplicit.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-2018 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 "EulerImplicit.H"
28 #include "simpleMatrix.H"
29 #include "Reaction.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 template<class ChemistryModel>
35 (
36  typename ChemistryModel::reactionThermo& thermo
37 )
38 :
40  coeffsDict_(this->subDict("EulerImplicitCoeffs")),
41  cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
42  eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter")),
43  cTp_(this->nEqns())
44 {}
45 
46 
47 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
48 
49 template<class ChemistryModel>
51 {}
52 
53 
54 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
55 
56 template<class ChemistryModel>
58 (
59  const label index,
60  const scalar pr,
61  const scalar pf,
62  const scalar corr,
63  const label lRef,
64  const label rRef,
65  const scalar p,
66  const scalar T,
68 ) const
69 {
71  this->reactions_[index];
72 
73  forAll(R.lhs(), s)
74  {
75  const label si = R.lhs()[s].index;
76  const scalar sl = R.lhs()[s].stoichCoeff;
77  RR[si][rRef] -= sl*pr*corr;
78  RR[si][lRef] += sl*pf*corr;
79  }
80 
81  forAll(R.rhs(), s)
82  {
83  const label si = R.rhs()[s].index;
84  const scalar sr = R.rhs()[s].stoichCoeff;
85  RR[si][lRef] -= sr*pf*corr;
86  RR[si][rRef] += sr*pr*corr;
87  }
88 }
89 
90 
91 template<class ChemistryModel>
93 (
94  scalarField& c,
95  scalar& T,
96  scalar& p,
97  scalar& deltaT,
98  scalar& subDeltaT
99 ) const
100 {
101  const label nSpecie = this->nSpecie();
102  simpleMatrix<scalar> RR(nSpecie, 0, 0);
103 
104  for (label i=0; i<nSpecie; i++)
105  {
106  c[i] = max(0, c[i]);
107  }
108 
109  // Calculate the absolute enthalpy
110  const scalar cTot = sum(c);
111  typename ChemistryModel::thermoType mixture
112  (
113  (this->specieThermo_[0].W()*c[0])*this->specieThermo_[0]
114  );
115  for (label i=1; i<nSpecie; i++)
116  {
117  mixture += (this->specieThermo_[i].W()*c[i])*this->specieThermo_[i];
118  }
119  const scalar ha = mixture.Ha(p, T);
120  const scalar deltaTEst = min(deltaT, subDeltaT);
121 
122  forAll(this->reactions(), i)
123  {
124  scalar pf, cf, pr, cr;
125  label lRef, rRef;
126 
127  const scalar omegai =
128  this->omegaI(i, c, T, p, pf, cf, lRef, pr, cr, rRef);
129 
130  scalar corr = 1;
131  if (eqRateLimiter_)
132  {
133  if (omegai < 0)
134  {
135  corr = 1/(1 + pr*deltaTEst);
136  }
137  else
138  {
139  corr = 1/(1 + pf*deltaTEst);
140  }
141  }
142 
143  updateRRInReactionI(i, pr, pf, corr, lRef, rRef, p, T, RR);
144  }
145 
146  // Calculate the stable/accurate time-step
147  scalar tMin = great;
148 
149  for (label i=0; i<nSpecie; i++)
150  {
151  scalar d = 0;
152  for (label j=0; j<nSpecie; j++)
153  {
154  d -= RR(i, j)*c[j];
155  }
156 
157  if (d < -small)
158  {
159  tMin = min(tMin, -(c[i] + small)/d);
160  }
161  else
162  {
163  d = max(d, small);
164  const scalar cm = max(cTot - c[i], 1e-5);
165  tMin = min(tMin, cm/d);
166  }
167  }
168 
169  subDeltaT = cTauChem_*tMin;
170  deltaT = min(deltaT, subDeltaT);
171 
172  // Add the diagonal and source contributions from the time-derivative
173  for (label i=0; i<nSpecie; i++)
174  {
175  RR(i, i) += 1/deltaT;
176  RR.source()[i] = c[i]/deltaT;
177  }
178 
179  // Solve for the new composition
180  c = RR.LUsolve();
181 
182  // Limit the composition
183  for (label i=0; i<nSpecie; i++)
184  {
185  c[i] = max(0, c[i]);
186  }
187 
188  // Update the temperature
189  mixture = (this->specieThermo_[0].W()*c[0])*this->specieThermo_[0];
190  for (label i=1; i<nSpecie; i++)
191  {
192  mixture += (this->specieThermo_[i].W()*c[i])*this->specieThermo_[i];
193  }
194  T = mixture.THa(ha, p, T);
195 }
196 
197 
198 // ************************************************************************* //
#define readScalar
Definition: doubleScalar.C:38
Field< Type > & source()
Return access to the source.
Definition: simpleMatrix.H:97
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
An abstract base class for solving chemistry.
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
const List< specieCoeffs > & lhs() const
Return the components of the left hand side.
Definition: ReactionI.H:58
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
rhoReactionThermo & thermo
Definition: createFields.H:28
const label nSpecie
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:55
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
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))
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:98
Field< Type > LUsolve() const
Solve the matrix using LU decomposition with pivoting.
Definition: simpleMatrix.C:86
A simple square matrix solver with scalar coefficients.
Definition: simpleMatrix.H:47
virtual void solve(scalarField &c, scalar &T, scalar &p, scalar &deltaT, scalar &subDeltaT) const
Update the concentrations and return the chemical time.
Definition: EulerImplicit.C:93
Info<< "Reading field p_rgh\"<< endl;volScalarField p_rgh(IOobject("p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Reading field U\"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Creating phaseChangeTwoPhaseMixture\"<< endl;autoPtr< phaseChangeTwoPhaseMixture > mixture
Definition: createFields.H:33
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
EulerImplicit(typename ChemistryModel::reactionThermo &thermo)
Construct from thermo.
Definition: EulerImplicit.C:35
#define R(A, B, C, D, E, F, K, M)
const scalarList W(::W(thermo))
const scalar RR
Universal gas constant (default in [J/(kmol K)])
void updateRRInReactionI(const label index, const scalar pr, const scalar pf, const scalar corr, const label lRef, const label rRef, const scalar p, const scalar T, simpleMatrix< scalar > &RR) const
Definition: EulerImplicit.C:58
virtual ~EulerImplicit()
Destructor.
Definition: EulerImplicit.C:50
const List< specieCoeffs > & rhs() const
Return the components of the right hand side.
Definition: ReactionI.H:66