FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/simpleMatrix.H>
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class CompType, class ThermoType>
34 (
36  const word& modelName
37 )
38 :
39  chemistrySolver<CompType, ThermoType>(model, modelName),
40  coeffsDict_(model.subDict(modelName + "Coeffs")),
41  cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
42  equil_(coeffsDict_.lookup("equilibriumRateLimiter"))
43 {}
44 
45 
46 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
47 
48 template<class CompType, class ThermoType>
50 {}
51 
52 
53 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
54 
55 template<class CompType, class ThermoType>
57 (
58  scalarField &c,
59  const scalar T,
60  const scalar p,
61  const scalar t0,
62  const scalar dt
63 ) const
64 {
65  scalar pf, cf, pr, cr;
66  label lRef, rRef;
67 
68  label nSpecie = this->model_.nSpecie();
69  simpleMatrix<scalar> RR(nSpecie, 0.0, 0.0);
70 
71  for (label i=0; i<nSpecie; i++)
72  {
73  c[i] = max(0.0, c[i]);
74  }
75 
76  for (label i=0; i<nSpecie; i++)
77  {
78  RR.source()[i] = c[i]/dt;
79  }
80 
81  for (label i=0; i<this->model_.reactions().size(); i++)
82  {
83  const Reaction<ThermoType>& R = this->model_.reactions()[i];
84 
85  scalar omegai = this->model_.omega
86  (
87  R, c, T, p, pf, cf, lRef, pr, cr, rRef
88  );
89 
90  scalar corr = 1.0;
91  if (equil_)
92  {
93  if (omegai<0.0)
94  {
95  corr = 1.0/(1.0 + pr*dt);
96  }
97  else
98  {
99  corr = 1.0/(1.0 + pf*dt);
100  }
101  }
102 
103  for (label s=0; s<R.lhs().size(); s++)
104  {
105  label si = R.lhs()[s].index;
106  scalar sl = R.lhs()[s].stoichCoeff;
107  RR[si][rRef] -= sl*pr*corr;
108  RR[si][lRef] += sl*pf*corr;
109  }
110 
111  for (label s=0; s<R.rhs().size(); s++)
112  {
113  label si = R.rhs()[s].index;
114  scalar sr = R.rhs()[s].stoichCoeff;
115  RR[si][lRef] -= sr*pf*corr;
116  RR[si][rRef] += sr*pr*corr;
117  }
118 
119  } // end for(label i...
120 
121 
122  for (label i=0; i<nSpecie; i++)
123  {
124  RR[i][i] += 1.0/dt;
125  }
126 
127  c = RR.LUsolve();
128  for (label i=0; i<nSpecie; i++)
129  {
130  c[i] = max(0.0, c[i]);
131  }
132 
133  // estimate the next time step
134  scalar tMin = GREAT;
135  label nEqns = this->model_.nEqns();
136  scalarField c1(nEqns, 0.0);
137 
138  for (label i=0; i<nSpecie; i++)
139  {
140  c1[i] = c[i];
141  }
142  c1[nSpecie] = T;
143  c1[nSpecie+1] = p;
144 
145  scalarField dcdt(nEqns, 0.0);
146  this->model_.derivatives(0.0, c1, dcdt);
147 
148  scalar sumC = sum(c);
149 
150  for (label i=0; i<nSpecie; i++)
151  {
152  scalar d = dcdt[i];
153  if (d < -SMALL)
154  {
155  tMin = min(tMin, -(c[i] + SMALL)/d);
156  }
157  else
158  {
159  d = max(d, SMALL);
160  scalar cm = max(sumC - c[i], 1.0e-5);
161  tMin = min(tMin, cm/d);
162  }
163  }
164 
165  return cTauChem_*tMin;
166 }
167 
168 
169 // ************************ vim: set sw=4 sts=4 et: ************************ //