Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "EulerImplicit.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <OpenFOAM/simpleMatrix.H>
00029
00030
00031
00032 template<class CompType, class ThermoType>
00033 Foam::EulerImplicit<CompType, ThermoType>::EulerImplicit
00034 (
00035 ODEChemistryModel<CompType, ThermoType>& model,
00036 const word& modelName
00037 )
00038 :
00039 chemistrySolver<CompType, ThermoType>(model, modelName),
00040 coeffsDict_(model.subDict(modelName + "Coeffs")),
00041 cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
00042 equil_(coeffsDict_.lookup("equilibriumRateLimiter"))
00043 {}
00044
00045
00046
00047
00048 template<class CompType, class ThermoType>
00049 Foam::EulerImplicit<CompType, ThermoType>::~EulerImplicit()
00050 {}
00051
00052
00053
00054
00055 template<class CompType, class ThermoType>
00056 Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
00057 (
00058 scalarField &c,
00059 const scalar T,
00060 const scalar p,
00061 const scalar t0,
00062 const scalar dt
00063 ) const
00064 {
00065 scalar pf, cf, pr, cr;
00066 label lRef, rRef;
00067
00068 label nSpecie = this->model_.nSpecie();
00069 simpleMatrix<scalar> RR(nSpecie, 0.0, 0.0);
00070
00071 for (label i=0; i<nSpecie; i++)
00072 {
00073 c[i] = max(0.0, c[i]);
00074 }
00075
00076 for (label i=0; i<nSpecie; i++)
00077 {
00078 RR.source()[i] = c[i]/dt;
00079 }
00080
00081 for (label i=0; i<this->model_.reactions().size(); i++)
00082 {
00083 const Reaction<ThermoType>& R = this->model_.reactions()[i];
00084
00085 scalar omegai = this->model_.omega
00086 (
00087 R, c, T, p, pf, cf, lRef, pr, cr, rRef
00088 );
00089
00090 scalar corr = 1.0;
00091 if (equil_)
00092 {
00093 if (omegai<0.0)
00094 {
00095 corr = 1.0/(1.0 + pr*dt);
00096 }
00097 else
00098 {
00099 corr = 1.0/(1.0 + pf*dt);
00100 }
00101 }
00102
00103 for (label s=0; s<R.lhs().size(); s++)
00104 {
00105 label si = R.lhs()[s].index;
00106 scalar sl = R.lhs()[s].stoichCoeff;
00107 RR[si][rRef] -= sl*pr*corr;
00108 RR[si][lRef] += sl*pf*corr;
00109 }
00110
00111 for (label s=0; s<R.rhs().size(); s++)
00112 {
00113 label si = R.rhs()[s].index;
00114 scalar sr = R.rhs()[s].stoichCoeff;
00115 RR[si][lRef] -= sr*pf*corr;
00116 RR[si][rRef] += sr*pr*corr;
00117 }
00118
00119 }
00120
00121
00122 for (label i=0; i<nSpecie; i++)
00123 {
00124 RR[i][i] += 1.0/dt;
00125 }
00126
00127 c = RR.LUsolve();
00128 for (label i=0; i<nSpecie; i++)
00129 {
00130 c[i] = max(0.0, c[i]);
00131 }
00132
00133
00134 scalar tMin = GREAT;
00135 label nEqns = this->model_.nEqns();
00136 scalarField c1(nEqns, 0.0);
00137
00138 for (label i=0; i<nSpecie; i++)
00139 {
00140 c1[i] = c[i];
00141 }
00142 c1[nSpecie] = T;
00143 c1[nSpecie+1] = p;
00144
00145 scalarField dcdt(nEqns, 0.0);
00146 this->model_.derivatives(0.0, c1, dcdt);
00147
00148 scalar sumC = sum(c);
00149
00150 for (label i=0; i<nSpecie; i++)
00151 {
00152 scalar d = dcdt[i];
00153 if (d < -SMALL)
00154 {
00155 tMin = min(tMin, -(c[i] + SMALL)/d);
00156 }
00157 else
00158 {
00159 d = max(d, SMALL);
00160 scalar cm = max(sumC - c[i], 1.0e-5);
00161 tMin = min(tMin, cm/d);
00162 }
00163 }
00164
00165 return cTauChem_*tMin;
00166 }
00167
00168
00169