src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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 "addToRunTimeSelectionTable.H"
  28 #include "simpleMatrix.H"
  29
  30 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  31
  32 template<class ODEChemistryType>
  33 Foam::EulerImplicit<ODEChemistryType>::EulerImplicit
  34 (
  35     const fvMesh& mesh,
  36     const word& ODEModelName,
  37     const word& thermoType
  38 )
  39 :
  40     chemistrySolver<ODEChemistryType>(mesh, ODEModelName, thermoType),
  41     coeffsDict_(this->subDict("EulerImplicitCoeffs")),
  42     cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
  43     eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter"))
  44 {}
  45
  46
  47 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  48
  49 template<class ODEChemistryType>
  50 Foam::EulerImplicit<ODEChemistryType>::~EulerImplicit()
  51 {}
  52
  53
  54 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  55
  56 template<class ODEChemistryType>
  57 Foam::scalar Foam::EulerImplicit<ODEChemistryType>::solve
  58 (
  59     scalarField &c,
  60     const scalar T,
  61     const scalar p,
  62     const scalar t0,
  63     const scalar dt
  64 ) const
  65 {
  66     scalar pf, cf, pr, cr;
  67     label lRef, rRef;
  68
  69     const label nSpecie = this->nSpecie();
  70     simpleMatrix<scalar> RR(nSpecie, 0, 0);
  71
  72     for (label i = 0; i < nSpecie; i++)
  73     {
  74         c[i] = max(0.0, c[i]);
  75     }
  76
  77     for (label i = 0; i < nSpecie; i++)
  78     {
  79         RR.source()[i] = c[i]/dt;
  80     }
  81
  82     forAll(this->reactions(), i)
  83     {
  84         scalar omegai = this->omegaI(i, c, T, p, pf, cf, lRef, pr, cr, rRef);
  85
  86         scalar corr = 1.0;
  87         if (eqRateLimiter_)
  88         {
  89             if (omegai < 0.0)
  90             {
  91                 corr = 1.0/(1.0 + pr*dt);
  92             }
  93             else
  94             {
  95                 corr = 1.0/(1.0 + pf*dt);
  96             }
  97         }
  98
  99         this->updateRRInReactionI(i, pr, pf, corr, lRef, rRef, RR);
 100     }
 101
 102
 103     for (label i = 0; i < nSpecie; i++)
 104     {
 105         RR[i][i] += 1.0/dt;
 106     }
 107
 108     c = RR.LUsolve();
 109     for (label i = 0; i < nSpecie; i++)
 110     {
 111         c[i] = max(0.0, c[i]);
 112     }
 113
 114     // estimate the next time step
 115     scalar tMin = GREAT;
 116     const label nEqns = this->nEqns();
 117     scalarField c1(nEqns, 0.0);
 118
 119     for (label i = 0; i < nSpecie; i++)
 120     {
 121         c1[i] = c[i];
 122     }
 123     c1[nSpecie] = T;
 124     c1[nSpecie+1] = p;
 125
 126     scalarField dcdt(nEqns, 0.0);
 127     this->derivatives(0.0, c1, dcdt);
 128
 129     const scalar sumC = sum(c);
 130
 131     for (label i = 0; i < nSpecie; i++)
 132     {
 133         scalar d = dcdt[i];
 134         if (d < -SMALL)
 135         {
 136             tMin = min(tMin, -(c[i] + SMALL)/d);
 137         }
 138         else
 139         {
 140             d = max(d, SMALL);
 141             scalar cm = max(sumC - c[i], 1.0e-5);
 142             tMin = min(tMin, cm/d);
 143         }
 144     }
 145
 146     return cTauChem_*tMin;
 147 }
 148
 149
 150 // ************************************************************************* //