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-2020 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"
27 #include "SubField.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class ChemistryModel>
34 (
35  const fluidReactionThermo& thermo
36 )
37 :
39  coeffsDict_(this->subDict("EulerImplicitCoeffs")),
40  cTauChem_(coeffsDict_.lookup<scalar>("cTauChem")),
41  cTp_(this->nEqns()),
42  R_(this->nEqns()),
43  J_(this->nEqns()),
44  E_(this->nEqns() - 2)
45 {}
46 
47 
48 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
49 
50 template<class ChemistryModel>
52 {}
53 
54 
55 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
56 
57 template<class ChemistryModel>
59 (
60  scalar& p,
61  scalar& T,
62  scalarField& c,
63  const label li,
64  scalar& deltaT,
65  scalar& subDeltaT
66 ) const
67 {
68  const label nSpecie = this->nSpecie();
69 
70  // Map the composition, temperature and pressure into cTp
71  for (int i=0; i<nSpecie; i++)
72  {
73  cTp_[i] = max(0, c[i]);
74  }
75  cTp_[nSpecie] = T;
76  cTp_[nSpecie + 1] = p;
77 
78  // Calculate the reaction rate and Jacobian
79  this->jacobian(0, cTp_, li, R_, J_);
80 
81  // Calculate the stable/accurate time-step
82  scalar tMin = great;
83  const scalar cTot = sum(c);
84 
85  for (label i=0; i<nSpecie; i++)
86  {
87  if (R_[i] < -small)
88  {
89  tMin = min(tMin, -(cTp_[i] + small)/R_[i]);
90  }
91  else
92  {
93  tMin = min
94  (
95  tMin,
96  max(cTot - cTp_[i], 1e-5)/max(R_[i], small)
97  );
98  }
99  }
100 
101  subDeltaT = cTauChem_*tMin;
102  deltaT = min(deltaT, subDeltaT);
103 
104  // Assemble the Euler implicit matrix for the composition
105  scalarField& source = E_.source();
106  for (label i=0; i<nSpecie; i++)
107  {
108  E_(i, i) = 1/deltaT - J_(i, i);
109  source[i] = R_[i] + E_(i, i)*cTp_[i];
110 
111  for (label j=0; j<nSpecie; j++)
112  {
113  if (i != j)
114  {
115  E_(i, j) = -J_(i, j);
116  source[i] += E_(i, j)*cTp_[j];
117  }
118  }
119  }
120 
121  // Solve for the new composition
122  scalarField::subField(cTp_, nSpecie) = E_.LUsolve();
123 
124  // Limit the composition and transfer back into c
125  for (label i=0; i<nSpecie; i++)
126  {
127  c[i] = max(0, cTp_[i]);
128  }
129 
130  // Euler explicit integrate the temperature.
131  // Separating the integration of temperature from composition
132  // is significantly more stable for exothermic systems
133  T += deltaT*R_[nSpecie];
134 }
135 
136 
137 // ************************************************************************* //
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
fluidReactionThermo & thermo
Definition: createFields.H:28
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Base-class for multi-component fluid thermodynamic properties.
const label nSpecie
Pre-declare related SubField type.
Definition: Field.H:60
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
virtual void solve(scalar &p, scalar &T, scalarField &c, const label li, scalar &deltaT, scalar &subDeltaT) const
Update the concentrations and return the chemical time.
Definition: EulerImplicit.C:59
EulerImplicit(const fluidReactionThermo &thermo)
Construct from thermo.
Definition: EulerImplicit.C:34
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const volScalarField & T
fvModels source(alpha1, mixture.thermo1().rho())
volScalarField & p
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105
virtual ~EulerImplicit()
Destructor.
Definition: EulerImplicit.C:51