EulerImplicit.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 | Copyright (C) 2011-2015 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 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class ChemistryModel>
34 (
35  const fvMesh& mesh,
36  const word& phaseName
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  label si = R.lhs()[s].index;
76  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  label si = R.rhs()[s].index;
84  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.0, c[i]);
107  }
108 
109  // Calculate the absolute enthalpy
110  scalar cTot = sum(c);
111  typename ChemistryModel::thermoType mixture
112  (
113  (c[0]/cTot)*this->specieThermo_[0]
114  );
115  for (label i=1; i<nSpecie; i++)
116  {
117  mixture += (c[i]/cTot)*this->specieThermo_[i];
118  }
119  scalar ha = mixture.Ha(p, T);
120 
121  scalar deltaTEst = min(deltaT, subDeltaT);
122 
123  forAll(this->reactions(), i)
124  {
125  scalar pf, cf, pr, cr;
126  label lRef, rRef;
127 
128  scalar omegai = this->omegaI(i, c, T, p, pf, cf, lRef, pr, cr, rRef);
129 
130  scalar corr = 1.0;
131  if (eqRateLimiter_)
132  {
133  if (omegai < 0.0)
134  {
135  corr = 1.0/(1.0 + pr*deltaTEst);
136  }
137  else
138  {
139  corr = 1.0/(1.0 + 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  scalar cm = max(cTot - c[i], 1.0e-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.0/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.0, c[i]);
186  }
187 
188  // Update the temperature
189  cTot = sum(c);
190  mixture = (c[0]/cTot)*this->specieThermo_[0];
191  for (label i=1; i<nSpecie; i++)
192  {
193  mixture += (c[i]/cTot)*this->specieThermo_[i];
194  }
195  T = mixture.THa(ha, p, T);
196 
197  /*
198  for (label i=0; i<nSpecie; i++)
199  {
200  cTp_[i] = c[i];
201  }
202  cTp_[nSpecie] = T;
203  cTp_[nSpecie+1] = p;
204  */
205 }
206 
207 
208 // ************************************************************************* //
Field< Type > & source()
Return access to the source.
Definition: simpleMatrix.H:97
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
A simple square matrix solver with scalar coefficients.
Definition: simpleMatrix.H:47
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 ))
label nSpecie
#define readScalar
Definition: doubleScalar.C:38
virtual ~EulerImplicit()
Destructor.
Definition: EulerImplicit.C:50
#define R(A, B, C, D, E, F, K, M)
EulerImplicit(const fvMesh &mesh, const word &phaseName)
Construct from mesh and phase name.
Definition: EulerImplicit.C:34
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
A class for handling words, derived from string.
Definition: word.H:59
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
Info<< "Reading field p_rgh\n"<< endl;volScalarField p_rgh(IOobject("p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Reading field U\n"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Creating phaseChangeTwoPhaseMixture\n"<< endl;autoPtr< phaseChangeTwoPhaseMixture > mixture
Definition: createFields.H:33
An abstract base class for solving chemistry.
dynamicFvMesh & mesh
const List< specieCoeffs > & rhs() const
Definition: ReactionI.H:59
Field< Type > LUsolve() const
Solve the matrix using LU decomposition with pivoting.
Definition: simpleMatrix.C:86
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
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Simple extension of ReactionThermo to handle reaction kinetics in addition to the equilibrium thermod...
Definition: Reaction.H:53
const List< specieCoeffs > & lhs() const
Definition: ReactionI.H:51
const scalar RR
Universal gas constant (default in [J/(kmol K)])
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)