FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

EulerImplicit.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "EulerImplicit.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <OpenFOAM/simpleMatrix.H>
00029 
00030 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
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 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00047 
00048 template<class CompType, class ThermoType>
00049 Foam::EulerImplicit<CompType, ThermoType>::~EulerImplicit()
00050 {}
00051 
00052 
00053 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
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     } // end for(label i...
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     // estimate the next time step
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 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines