.gitignore

@@ -82,20 +82,3 @@ doc/Doxygen/DTAGS
 
 # Ignore the test directory
 /tutorialsTest
-
-.*.swp
-*.gz
-log.*
-postProcessing
-*~
-processor*
-[0-9]*
-[0-9]*.[0-9]*
-!0.org
-constant/extendedFeatureEdgeMesh
-*.eMesh
-constant/polyMesh/*
-!constant/polyMesh/blockMeshDict
-constant/cellToRegion
-*.OpenFOAM
-*.foam

examples/FlameD_LagrangianCMC/constant/thermophysicalProperties

@@ -21,7 +21,7 @@ thermoType
     mixture         reactingMixture;
     transport       sutherland;
     thermo          janaf;
-    energy          absoluteEnthalpy;
+    energy          sensibleEnthalpy;
     equationOfState perfectGas;
     specie          specie;
 }

examples/FlameD_LagrangianCMC/system/fvSchemes

@@ -33,7 +33,7 @@ divSchemes
     div(phi,U)      Gauss linearUpwindV cellLimited Gauss linear 1;
     div(phi,Yi_h)   Gauss multivariateSelection
     {
-	ha linearUpwind cellLimited Gauss linear 1;
+	h linearUpwind cellLimited Gauss linear 1;
 	H2 linearUpwind cellLimited Gauss linear 1;
 	H linearUpwind cellLimited Gauss linear 1;
 	O linearUpwind cellLimited Gauss linear 1;

examples/FlameD_LagrangianCMC/system/fvSolution

@@ -46,14 +46,14 @@ solvers
         relTol           0.0;
     }
 
-    "(U|Yi|ha|k|epsilon|mf|mfVar|Fk)"
+    "(U|Yi|h|k|epsilon|mf|mfVar|Fk)"
     {
         solver          PBiCG;
         preconditioner  DILU;
         tolerance       1e-08;
         relTol          0;
     }
-    "(U|Yi|ha|k|epsilon|mf|mfVar|Fk)Final"
+    "(U|Yi|h|k|epsilon|mf|mfVar|Fk)Final"
     {
         solver          PBiCG;
         preconditioner  DILU;

libs/Allwmake

@@ -5,6 +5,6 @@ targetType=libso
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
 
 wmake $targetType chemistryModel_POSTECH
-wmake $targetType thermos
+wmake $targetType combustionModels_POSTECH
 
 #------------------------------------------------------------------------------

libs/combustionModels_POSTECH/FSD/FSD.C

@@ -0,0 +1,370 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "FSD.H"
+#include "addToRunTimeSelectionTable.H"
+#include "LESModel.H"
+#include "fvcGrad.H"
+#include "fvcDiv.H"
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+FSD<CombThermoType, ThermoType>::FSD
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    singleStepCombustion<CombThermoType, ThermoType>
+    (
+        modelType,
+        mesh,
+        phaseName
+    ),
+    reactionRateFlameArea_
+    (
+        reactionRateFlameArea::New
+        (
+            this->coeffs(),
+            this->mesh(),
+            *this
+        )
+    ),
+    ft_
+    (
+        IOobject
+        (
+            IOobject::groupName("ft", phaseName),
+            this->mesh().time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh(),
+        dimensionedScalar("zero", dimless, 0.0)
+    ),
+    YFuelFuelStream_(dimensionedScalar("YFuelStream", dimless, 1.0)),
+    YO2OxiStream_(dimensionedScalar("YOxiStream", dimless, 0.23)),
+    Cv_(readScalar(this->coeffs().lookup("Cv"))),
+    C_(5.0),
+    ftMin_(0.0),
+    ftMax_(1.0),
+    ftDim_(300),
+    ftVarMin_(readScalar(this->coeffs().lookup("ftVarMin")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+FSD<CombThermoType, ThermoType>::~FSD()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
+{
+    this->singleMixturePtr_->fresCorrect();
+
+    const label fuelI = this->singleMixturePtr_->fuelIndex();
+
+    const volScalarField& YFuel = this->thermoPtr_->composition().Y()[fuelI];
+
+    const volScalarField& YO2 = this->thermoPtr_->composition().Y("O2");
+
+    const dimensionedScalar s = this->singleMixturePtr_->s();
+
+    ft_ =
+        (s*YFuel - (YO2 - YO2OxiStream_))/(s*YFuelFuelStream_ + YO2OxiStream_);
+
+
+    volVectorField nft(fvc::grad(ft_));
+
+    volScalarField mgft(mag(nft));
+
+    surfaceVectorField SfHat(this->mesh().Sf()/this->mesh().magSf());
+
+    volScalarField cAux(scalar(1) - ft_);
+
+    dimensionedScalar dMgft = 1.0e-3*
+        (ft_*cAux*mgft)().weightedAverage(this->mesh().V())
+       /((ft_*cAux)().weightedAverage(this->mesh().V()) + SMALL)
+      + dimensionedScalar("ddMgft", mgft.dimensions(), SMALL);
+
+    mgft += dMgft;
+
+    nft /= mgft;
+
+    const volVectorField& U = YO2.db().lookupObject<volVectorField>("U");
+
+    const volScalarField sigma
+    (
+        (nft & nft)*fvc::div(U) - (nft & fvc::grad(U) & nft)
+    );
+
+    reactionRateFlameArea_->correct(sigma);
+
+    const volScalarField& omegaFuel = reactionRateFlameArea_->omega();
+
+
+    const scalar ftStoich =
+        YO2OxiStream_.value()
+       /(
+            s.value()*YFuelFuelStream_.value() + YO2OxiStream_.value()
+        );
+
+    tmp<volScalarField> tPc
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("Pc", this->phaseName_),
+                U.time().timeName(),
+                U.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U.mesh(),
+            dimensionedScalar("Pc", dimless, 0)
+        )
+    );
+
+    volScalarField& pc = tPc.ref();
+
+    tmp<volScalarField> tomegaFuel
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("omegaFuelBar", this->phaseName_),
+                U.time().timeName(),
+                U.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U.mesh(),
+            dimensionedScalar
+            (
+                "omegaFuelBar",
+                omegaFuel.dimensions(),
+                0
+            )
+        )
+    );
+
+    volScalarField& omegaFuelBar = tomegaFuel.ref();
+
+    // Calculation of the mixture fraction variance (ftVar)
+    const compressible::LESModel& lesModel =
+        YO2.db().lookupObject<compressible::LESModel>
+        (
+            turbulenceModel::propertiesName
+        );
+
+    const volScalarField& delta = lesModel.delta();
+    const volScalarField ftVar(Cv_*sqr(delta)*sqr(mgft));
+
+    // Thickened flame (average flame thickness for counterflow configuration
+    // is 1.5 mm)
+
+    volScalarField  deltaF
+    (
+        lesModel.delta()/dimensionedScalar("flame", dimLength, 1.5e-3)
+    );
+
+    // Linear correlation between delta and flame thickness
+    volScalarField omegaF(max(deltaF*(4.0/3.0) + (2.0/3.0), scalar(1)));
+
+    scalar deltaFt = 1.0/ftDim_;
+
+    forAll(ft_, celli)
+    {
+        if (ft_[celli] > ftMin_ && ft_[celli] < ftMax_)
+        {
+            scalar ftCell = ft_[celli];
+
+            if (ftVar[celli] > ftVarMin_) //sub-grid beta pdf of ft_
+            {
+                scalar ftVarc = ftVar[celli];
+                scalar a =
+                    max(ftCell*(ftCell*(1.0 - ftCell)/ftVarc - 1.0), 0.0);
+                scalar b = max(a/ftCell - a, 0.0);
+
+                for (int i=1; i<ftDim_; i++)
+                {
+                    scalar ft = i*deltaFt;
+                    pc[celli] += pow(ft, a-1.0)*pow(1.0 - ft, b - 1.0)*deltaFt;
+                }
+
+                for (int i=1; i<ftDim_; i++)
+                {
+                    scalar ft = i*deltaFt;
+                    omegaFuelBar[celli] +=
+                        omegaFuel[celli]/omegaF[celli]
+                       *exp
+                        (
+                           -sqr(ft - ftStoich)
+                           /(2.0*sqr(0.01*omegaF[celli]))
+                        )
+                       *pow(ft, a - 1.0)
+                       *pow(1.0 - ft, b - 1.0)
+                       *deltaFt;
+                }
+                omegaFuelBar[celli] /= max(pc[celli], 1e-4);
+            }
+            else
+            {
+                omegaFuelBar[celli] =
+                   omegaFuel[celli]/omegaF[celli]
+                  *exp(-sqr(ftCell - ftStoich)/(2.0*sqr(0.01*omegaF[celli])));
+            }
+        }
+        else
+        {
+            omegaFuelBar[celli] = 0.0;
+        }
+    }
+
+
+    // Combustion progress variable, c
+
+    List<label> productsIndex(2, label(-1));
+    {
+        label i = 0;
+        forAll(this->singleMixturePtr_->specieProd(), specieI)
+        {
+            if (this->singleMixturePtr_->specieProd()[specieI] < 0)
+            {
+                productsIndex[i] = specieI;
+                i++;
+            }
+        }
+    }
+
+
+    // Flamelet probability of the progress c based on IFC (reuse pc)
+    scalar YprodTotal = 0;
+    forAll(productsIndex, j)
+    {
+        YprodTotal += this->singleMixturePtr_->Yprod0()[productsIndex[j]];
+    }
+
+    forAll(ft_, celli)
+    {
+        if (ft_[celli] < ftStoich)
+        {
+            pc[celli] = ft_[celli]*(YprodTotal/ftStoich);
+        }
+        else
+        {
+            pc[celli] = (1.0 - ft_[celli])*(YprodTotal/(1.0 - ftStoich));
+        }
+    }
+
+    tmp<volScalarField> tproducts
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("products", this->phaseName_),
+                U.time().timeName(),
+                U.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U.mesh(),
+            dimensionedScalar("products", dimless, 0)
+        )
+    );
+
+    volScalarField& products = tproducts.ref();
+
+    forAll(productsIndex, j)
+    {
+        label specieI = productsIndex[j];
+        const volScalarField& Yp = this->thermoPtr_->composition().Y()[specieI];
+        products += Yp;
+    }
+
+    volScalarField c
+    (
+        max(scalar(1) - products/max(pc, scalar(1e-5)), scalar(0))
+    );
+
+    pc = min(C_*c, scalar(1));
+
+    const volScalarField fres(this->singleMixturePtr_->fres(fuelI));
+
+    this->wFuel_ == mgft*pc*omegaFuelBar;
+}
+
+
+template<class CombThermoType, class ThermoType>
+void FSD<CombThermoType, ThermoType>::correct()
+{
+    this->wFuel_ ==
+        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);
+
+    if (this->active())
+    {
+        calculateSourceNorm();
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+bool FSD<CombThermoType, ThermoType>::read()
+{
+    if (singleStepCombustion<CombThermoType, ThermoType>::read())
+    {
+        this->coeffs().lookup("Cv") >> Cv_ ;
+        this->coeffs().lookup("ftVarMin") >> ftVarMin_;
+        reactionRateFlameArea_->read(this->coeffs());
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs/combustionModels_POSTECH/FSD/FSD.H

@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::FSD
+
+Description
+
+    Flame Surface Dennsity (FDS) combustion model.
+
+    The fuel source term is given by mgft*pc*omegaFuelBar.
+
+    where:
+          mgft: filtered flame area.
+          pc:   probability of the combustion progress.
+          omegaFuelBar: filtered consumption speed per unit of flame area.
+
+    pc is considered from the IFC solution.
+    omegaFuelBar is calculated solving a relaxation equation which tends to
+    omegaEq. This omegaEq is obtained from the flamelet solution for
+    different strain rates and fit using a expential distribution.
+
+    The spacial distribution of the consumption speed (omega) is obtained also
+    from a strained flamelet solution and it is assumed to have a guassian
+    distribution.
+
+    If the grid resolution is not enough to resolve the flame, the consumption
+    speed distribution is linearly thickened conserving the overall heat
+    release.
+
+    If the turbulent fluctuation of the mixture fraction at the sub-grid level
+    is large (>1e-04) then a beta pdf is used for filtering.
+
+    At the moment the flame area combustion model is only fit to work in a LES
+    frame work. In RAS the subgrid fluctuation has to be solved by an extra
+    transport equation.
+
+SourceFiles
+    FSD.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef FSD_H
+#define FSD_H
+
+#include "singleStepCombustion.H"
+#include "reactionRateFlameArea.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class FSD Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType, class ThermoType>
+class FSD
+:
+    public singleStepCombustion <CombThermoType, ThermoType>
+{
+    // Private data
+
+        //- Auto pointer to consumption speed per unit of flame area model
+        autoPtr<reactionRateFlameArea> reactionRateFlameArea_;
+
+        //- Mixture fraction
+        volScalarField ft_;
+
+        //- Fuel mass concentration on the fuel stream
+        dimensionedScalar YFuelFuelStream_;
+
+        //- Oxygen mass concentration on the oxydizer stream
+        dimensionedScalar YO2OxiStream_;
+
+        //- Similarity constant for the sub-grid ft fluctuations
+        scalar Cv_;
+
+        //- Model constant
+        scalar C_;
+
+        //- Lower flammability limit
+        scalar ftMin_;
+
+        //- Upper flammability limit
+        scalar ftMax_;
+
+        //- Dimension of the ft space. Used to integrate the beta-pdf
+        scalar ftDim_;
+
+        //- Minimum mixture freaction variance to calculate pdf
+        scalar ftVarMin_;
+
+
+    // Private Member Functions
+
+        //- Calculate the normalised fuel source term
+        void calculateSourceNorm();
+
+        //- Disallow copy construct
+        FSD(const FSD&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const FSD&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("FSD");
+
+
+    // Constructors
+
+        //- Construct from components
+        FSD(const word& modelType, const fvMesh& mesh, const word& phaseName);
+
+
+    // Destructor
+        virtual ~FSD();
+
+
+    // Evolution
+
+        //- Correct combustion rate
+        virtual void correct();
+
+
+    // IO
+
+        //- Update properties
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "FSD.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/FSD/FSDs.C

@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "thermoPhysicsTypes.H"
+#include "psiCombustionModel.H"
+#include "psiThermoCombustion.H"
+#include "rhoCombustionModel.H"
+#include "rhoThermoCombustion.H"
+#include "FSD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Combustion models based on sensibleEnthalpy
+makeCombustionTypesThermo
+(
+    FSD,
+    psiThermoCombustion,
+    gasHThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    FSD,
+    psiThermoCombustion,
+    constGasHThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    FSD,
+    rhoThermoCombustion,
+    gasHThermoPhysics,
+    rhoCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    FSD,
+    rhoThermoCombustion,
+    constGasHThermoPhysics,
+    rhoCombustionModel
+);
+
+// Combustion models based on sensibleInternalEnergy
+makeCombustionTypesThermo
+(
+    FSD,
+    psiThermoCombustion,
+    gasEThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    FSD,
+    psiThermoCombustion,
+    constGasEThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    FSD,
+    rhoThermoCombustion,
+    gasEThermoPhysics,
+    rhoCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    FSD,
+    rhoThermoCombustion,
+    constGasEThermoPhysics,
+    rhoCombustionModel
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs/combustionModels_POSTECH/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "consumptionSpeed.H"
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+namespace Foam
+{
+    defineTypeNameAndDebug(consumptionSpeed, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::consumptionSpeed::consumptionSpeed
+(
+    const dictionary& dict
+)
+:   omega0_(readScalar(dict.lookup("omega0"))),
+    eta_(readScalar(dict.lookup("eta"))),
+    sigmaExt_(readScalar(dict.lookup("sigmaExt"))),
+    omegaMin_(readScalar(dict.lookup("omegaMin")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::consumptionSpeed::~consumptionSpeed()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::consumptionSpeed::omega0Sigma
+(
+    scalar sigma,
+    scalar a
+) const
+{
+    scalar omega0 = 0.0;
+
+    if (sigma < sigmaExt_)
+    {
+        omega0 = max
+        (
+            a*omega0_*(1.0 - exp(eta_*max(sigma, 0.0))),
+            omegaMin_
+        ) ;
+    }
+
+    return omega0;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::consumptionSpeed::omega0Sigma
+(
+    const volScalarField& sigma
+)
+{
+    tmp<volScalarField> tomega0
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "omega0",
+                sigma.time().timeName(),
+                sigma.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            sigma.mesh(),
+            dimensionedScalar
+            (
+                "omega0",
+                dimensionSet(1, -2, -1, 0, 0, 0, 0),
+                0
+            )
+        )
+    );
+
+    volScalarField& omega0 = tomega0.ref();
+
+    volScalarField::Internal& iomega0 = omega0;
+
+    forAll(iomega0, celli)
+    {
+        iomega0[celli] = omega0Sigma(sigma[celli], 1.0);
+    }
+
+    volScalarField::Boundary& bomega0 = omega0.boundaryFieldRef();
+
+    forAll(bomega0, patchi)
+    {
+        forAll(bomega0[patchi], facei)
+        {
+            bomega0[patchi][facei] =
+                omega0Sigma
+                (
+                    sigma.boundaryField()[patchi][facei],
+                    1.0
+                );
+        }
+    }
+
+    return tomega0;
+}
+
+
+void  Foam::consumptionSpeed::read(const dictionary& dict)
+{
+    dict.lookup("omega0") >> omega0_ ;
+    dict.lookup("eta") >> eta_ ;
+    dict.lookup("sigmaExt") >> sigmaExt_;
+    dict.lookup("omegaMin") >> omegaMin_;
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.H

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::reactionRateFlameAreaModels::consumptionSpeed
+
+Description
+    Correlation function for laminar consumption speed obtained from flamelet
+    solution at increasing strain rates.
+
+SourceFiles
+    consumptionSpeed.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef consumptionSpeed_H
+#define consumptionSpeed_H
+
+#include "IOdictionary.H"
+#include "volFields.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class consumptionSpeed Declaration
+\*---------------------------------------------------------------------------*/
+
+class consumptionSpeed
+{
+    // Private Data
+
+        //- Maximum consumption speed
+        scalar omega0_;
+
+        //- Exponential factor
+        scalar eta_;
+
+        //- Extinction strain
+        scalar sigmaExt_;
+
+        //- Minimum consumption speed
+        scalar omegaMin_;
+
+
+    // Private member functions
+
+        //- Return consumption rate
+        scalar omega0Sigma(scalar sigma, scalar a) const;
+
+        //- Disallow copy construct
+        consumptionSpeed(const consumptionSpeed&);
+
+         //- Disallow default bitwise assignment
+        void operator=(const consumptionSpeed&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("consumptionSpeed");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        consumptionSpeed(const dictionary& dict);
+
+
+    //- Destructor
+        virtual ~consumptionSpeed();
+
+
+    // Member functions
+
+        //- Return speed consumption rate temp
+        tmp<volScalarField> omega0Sigma(const volScalarField& sigma);
+
+
+        // Access functions
+
+            scalar omega0() const
+            {
+                return omega0_;
+            }
+
+            scalar eta() const
+            {
+                return eta_;
+            }
+
+            scalar sigmaExt() const
+            {
+                return sigmaExt_;
+            }
+
+            scalar omegaMin() const
+            {
+                return omegaMin_;
+            }
+
+
+    // IO
+
+        //- Update properties
+        void read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameArea.C

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "reactionRateFlameArea.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(reactionRateFlameArea, 0);
+    defineRunTimeSelectionTable(reactionRateFlameArea, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::reactionRateFlameArea::reactionRateFlameArea
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const combustionModel& combModel
+)
+:
+    coeffDict_(dictionary::null),
+    mesh_(mesh),
+    combModel_(combModel),
+    fuel_(dict.lookup("fuel")),
+    omega_
+    (
+        IOobject
+        (
+            "omega",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    )
+{}
+
+
+Foam::reactionRateFlameArea::reactionRateFlameArea
+(
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const combustionModel& combModel
+)
+:
+    coeffDict_(dict.subDict(modelType + "Coeffs")),
+    mesh_(mesh),
+    combModel_(combModel),
+    fuel_(dict.lookup("fuel")),
+    omega_
+    (
+        IOobject
+        (
+            "omega",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::reactionRateFlameArea::~reactionRateFlameArea()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::reactionRateFlameArea::read(const dictionary& dict)
+{
+    dict.lookup("fuel") >> fuel_;
+
+    return true;
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameArea.H

@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::reactionRateFlameArea
+
+Description
+    Abstract class for reaction rate per flame area unit
+
+SourceFiles
+    reactionRateFlameArea.C
+    reactionRateFlameAreaNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef reactionRateFlameArea_H
+#define reactionRateFlameArea_H
+
+#include "runTimeSelectionTables.H"
+#include "dictionary.H"
+#include "autoPtr.H"
+#include "volFields.H"
+#include "combustionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class reactionRateFlameArea Declaration
+\*---------------------------------------------------------------------------*/
+
+class reactionRateFlameArea
+{
+
+protected:
+
+    // Protected data
+
+        //- Dictionary
+        dictionary coeffDict_;
+
+        //- Mesh reference
+        const fvMesh& mesh_;
+
+        //- Combstion model owner
+        const combustionModel& combModel_;
+
+        //- Fuel name
+        word fuel_;
+
+        //- Fuel consumption rate per unit of flame area
+        volScalarField omega_;
+
+
+private:
+
+    // Private member functions
+
+         //- Disallow copy construct
+        reactionRateFlameArea(const reactionRateFlameArea&);
+
+         //- Disallow default bitwise assignment
+        void operator=(const reactionRateFlameArea&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("reactionRateFlameArea");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            reactionRateFlameArea,
+            dictionary,
+            (
+                const word modelType,
+                const dictionary& dict,
+                const fvMesh& mesh,
+                const combustionModel& combModel
+            ),
+            (modelType, dict, mesh, combModel)
+        );
+
+
+    // Constructors
+
+        //- Construct from dictionary and psiReactionThermo
+        reactionRateFlameArea
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const combustionModel& combModel
+        );
+
+        //- Construct from components
+        reactionRateFlameArea
+        (
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const combustionModel& combModel
+        );
+
+
+    // Selector
+
+        static autoPtr<reactionRateFlameArea> New
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const combustionModel& combModel
+        );
+
+
+    // Destructor
+
+        virtual ~reactionRateFlameArea();
+
+
+    // Member functions
+
+        //- Access functions
+
+            //- Return omega
+            const volScalarField& omega() const
+            {
+                return omega_;
+            }
+
+
+        //- Correct omega
+        virtual void correct(const volScalarField& sigma) = 0;
+
+        //- Update from dictionary
+        virtual bool read(const dictionary& dictProperties);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameAreaNew.C

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "reactionRateFlameArea.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::reactionRateFlameArea> Foam::reactionRateFlameArea::New
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const combustionModel& combModel
+)
+{
+    word reactionRateFlameAreaType
+    (
+        dict.lookup("reactionRateFlameArea")
+    );
+
+    Info<< "Selecting reaction rate flame area correlation "
+        << reactionRateFlameAreaType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(reactionRateFlameAreaType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalIOErrorInFunction
+        (
+            dict
+        )   << "Unknown reactionRateFlameArea type "
+            << reactionRateFlameAreaType << endl << endl
+            << "Valid reaction rate flame area types are :" << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalIOError);
+    }
+
+    const label tempOpen = reactionRateFlameAreaType.find('<');
+
+    const word className = reactionRateFlameAreaType(0, tempOpen);
+
+    return autoPtr<reactionRateFlameArea>
+        (cstrIter()(className, dict, mesh, combModel));
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/FSD/reactionRateFlameAreaModels/relaxation/relaxation.C

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "relaxation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvm.H"
+#include "LESModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace reactionRateFlameAreaModels
+{
+    defineTypeNameAndDebug(relaxation, 0);
+    addToRunTimeSelectionTable
+    (
+        reactionRateFlameArea,
+        relaxation,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::reactionRateFlameAreaModels::relaxation::relaxation
+(
+    const word modelType,
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const combustionModel& combModel
+)
+:
+    reactionRateFlameArea(modelType, dict, mesh, combModel),
+    correlation_(dict.subDict(typeName + "Coeffs").subDict(fuel_)),
+    C_(readScalar(dict.subDict(typeName + "Coeffs").lookup("C"))),
+    alpha_(readScalar(dict.subDict(typeName + "Coeffs").lookup("alpha")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::reactionRateFlameAreaModels::relaxation::~relaxation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::reactionRateFlameAreaModels::relaxation::correct
+(
+    const volScalarField& sigma
+)
+{
+    dimensionedScalar omega0
+    (
+        "omega0",
+        dimensionSet(1, -2, -1, 0, 0, 0, 0),
+        correlation_.omega0()
+    );
+
+    dimensionedScalar sigmaExt
+    (
+        "sigmaExt",
+        dimensionSet(0, 0, -1, 0, 0, 0, 0),
+        correlation_.sigmaExt()
+    );
+
+    dimensionedScalar omegaMin
+    (
+        "omegaMin",
+        omega0.dimensions(),
+        1e-4
+    );
+
+    dimensionedScalar kMin
+    (
+        "kMin",
+        sqr(dimVelocity),
+        SMALL
+    );
+
+    const compressibleTurbulenceModel& turbulence = combModel_.turbulence();
+
+    // Total strain
+    const volScalarField sigmaTotal
+    (
+        sigma + alpha_*turbulence.epsilon()/(turbulence.k() + kMin)
+    );
+
+    const volScalarField omegaInf(correlation_.omega0Sigma(sigmaTotal));
+
+    dimensionedScalar sigma0("sigma0", sigma.dimensions(), 0.0);
+
+    const volScalarField tau(C_*mag(sigmaTotal));
+
+    volScalarField Rc
+    (
+        (tau*omegaInf*(omega0 - omegaInf) + sqr(omegaMin)*sigmaExt)
+       /(sqr(omega0 - omegaInf) + sqr(omegaMin))
+    );
+
+    const volScalarField rho(combModel_.rho());
+    const surfaceScalarField phi(combModel_.phi());
+
+    solve
+    (
+         fvm::ddt(rho, omega_)
+       + fvm::div(phi, omega_, "div(phi,omega)")
+      ==
+         rho*Rc*omega0
+       - fvm::SuSp(rho*(tau + Rc), omega_)
+    );
+
+    omega_.min(omega0);
+    omega_.max(0.0);
+}
+
+
+bool  Foam::reactionRateFlameAreaModels::relaxation::read
+(
+    const dictionary& dict
+)
+{
+    if (reactionRateFlameArea::read(dict))
+    {
+        coeffDict_ = dict.subDict(typeName + "Coeffs");
+        coeffDict_.lookup("C") >> C_;
+        coeffDict_.lookup("alpha") >> alpha_;
+        correlation_.read
+        (
+            coeffDict_.subDict(fuel_)
+        );
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/FSD/reactionRateFlameAreaModels/relaxation/relaxation.H

@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::reactionRateFlameAreaModels::relaxation
+
+Description
+    Consumption rate per unit of flame area obtained from a relaxation equation
+
+SourceFiles
+    relaxation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef relaxation_H
+#define relaxation_H
+
+#include "reactionRateFlameArea.H"
+#include "consumptionSpeed.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace reactionRateFlameAreaModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class relaxation Declaration
+\*---------------------------------------------------------------------------*/
+
+class relaxation
+:
+    public reactionRateFlameArea
+{
+    // Private Data
+
+        //- Correlation
+        consumptionSpeed correlation_;
+
+        //- Proportionality constant for time scale in the relaxation Eq.
+        scalar C_;
+
+        //- Proportionality constant for sub-grid strain
+        scalar alpha_;
+
+
+     // Private Member Functions
+
+        //- Disallow copy construct
+        relaxation(const relaxation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const relaxation&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("relaxation");
+
+
+    // Constructors
+
+        //- Construct from dictionary and psiReactionThermo
+        relaxation
+        (
+            const word modelType,
+            const dictionary& dictCoeffs,
+            const fvMesh& mesh,
+            const combustionModel& combModel
+        );
+
+
+    // Destructor
+
+        virtual ~relaxation();
+
+
+    // Member functions
+
+        //- Correct omega
+        virtual void correct(const volScalarField& sigma);
+
+
+    // IO
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& dictProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End reactionRateFlameAreaModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/Make/files

@@ -0,0 +1,29 @@
+combustionModel/combustionModel.C
+
+psiCombustionModel/psiCombustionModel/psiCombustionModel.C
+psiCombustionModel/psiCombustionModel/psiCombustionModelNew.C
+psiCombustionModel/psiThermoCombustion/psiThermoCombustion.C
+psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.C
+
+rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.C
+rhoCombustionModel/rhoCombustionModel/rhoCombustionModelNew.C
+rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.C
+rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
+
+diffusion/diffusions.C
+infinitelyFastChemistry/infinitelyFastChemistrys.C
+
+PaSR/PaSRs.C
+
+laminar/laminars.C
+
+FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
+FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameArea.C
+FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameAreaNew.C
+FSD/reactionRateFlameAreaModels/relaxation/relaxation.C
+
+FSD/FSDs.C
+
+noCombustion/noCombustions.C
+
+LIB = $(FOAM_USER_LIBBIN)/libcombustionModels_POSTECH

libs/combustionModels_POSTECH/Make/options

@@ -0,0 +1,18 @@
+EXE_INC = \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+LIB_LIBS = \
+    -lcompressibleTransportModels \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lchemistryModel \
+    -lfiniteVolume \
+    -lmeshTools

libs/combustionModels_POSTECH/PaSR/PaSR.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PaSR.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::combustionModels::PaSR<Type>::PaSR
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    laminar<Type>(modelType, mesh, phaseName),
+    Cmix_(readScalar(this->coeffs().lookup("Cmix"))),
+    turbulentReaction_(this->coeffs().lookup("turbulentReaction")),
+    kappa_
+    (
+        IOobject
+        (
+            IOobject::groupName("PaSR:kappa", phaseName),
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("kappa", dimless, 0.0)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::combustionModels::PaSR<Type>::~PaSR()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::combustionModels::PaSR<Type>::correct()
+{
+    if (this->active())
+    {
+        laminar<Type>::correct();
+
+        if (turbulentReaction_)
+        {
+            tmp<volScalarField> tepsilon(this->turbulence().epsilon());
+            const volScalarField& epsilon = tepsilon();
+            tmp<volScalarField> tmuEff(this->turbulence().muEff());
+            const volScalarField& muEff = tmuEff();
+            tmp<volScalarField> ttc(this->tc());
+            const volScalarField& tc = ttc();
+            tmp<volScalarField> trho(this->rho());
+            const volScalarField& rho = trho();
+
+            forAll(epsilon, i)
+            {
+                scalar tk =
+                    Cmix_*sqrt(max(muEff[i]/rho[i]/(epsilon[i] + SMALL), 0));
+
+                if (tk > SMALL)
+                {
+                    kappa_[i] = tc[i]/(tc[i] + tk);
+                }
+                else
+                {
+                    kappa_[i] = 1.0;
+                }
+            }
+        }
+        else
+        {
+            kappa_ = 1.0;
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::combustionModels::PaSR<Type>::R(volScalarField& Y) const
+{
+    return kappa_*laminar<Type>::R(Y);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::PaSR<Type>::dQ() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject::groupName("PaSR:dQ", this->phaseName_),
+            kappa_*laminar<Type>::dQ()
+        )
+    );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::PaSR<Type>::Sh() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject::groupName("PaSR:Sh", this->phaseName_),
+            kappa_*laminar<Type>::Sh()
+        )
+    );
+}
+
+
+template<class Type>
+bool Foam::combustionModels::PaSR<Type>::read()
+{
+    if (laminar<Type>::read())
+    {
+        this->coeffs().lookup("Cmix") >> Cmix_;
+        this->coeffs().lookup("turbulentReaction") >> turbulentReaction_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/PaSR/PaSR.H

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::PaSR
+
+Description
+    Partially stirred reactor combustion model.  The model calculates a finite
+    rate, based on both turbulence and chemistry time scales.  Depending on
+    mesh resolution, the Cmix parameter can be used to scale the turbulence
+    mixing time scale.
+
+SourceFiles
+    PaSR.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PaSR_H
+#define PaSR_H
+
+#include "../laminar/laminar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class PaSR Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class PaSR
+:
+    public laminar<Type>
+{
+    // Private data
+
+        //- Mixing constant
+        scalar Cmix_;
+
+        //- Turbulent reaction switch
+        Switch turbulentReaction_;
+
+        //- Mixing parameter
+        volScalarField kappa_;
+
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        PaSR(const PaSR&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const PaSR&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("PaSR");
+
+
+    // Constructors
+
+        //- Construct from components
+        PaSR(const word& modelType, const fvMesh& mesh, const word& phaseName);
+
+
+    //- Destructor
+    virtual ~PaSR();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+            //- Fuel consumption rate matrix.
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
+
+            //- Heat release rate calculated from fuel consumption rate matrix
+            virtual tmp<volScalarField> dQ() const;
+
+            //-  Return source for enthalpy equation [kg/m/s3]
+            virtual tmp<volScalarField> Sh() const;
+
+
+    // IO
+
+            //- Update properties from given dictionary
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "PaSR.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/PaSR/PaSRs.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "psiChemistryCombustion.H"
+#include "rhoChemistryCombustion.H"
+#include "PaSR.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeCombustionTypes(PaSR, psiChemistryCombustion, psiCombustionModel);
+makeCombustionTypes(PaSR, rhoChemistryCombustion, rhoCombustionModel);
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/combustionModel/combustionModel.C

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "combustionModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(combustionModel, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::combustionModel::combustionModel
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            IOobject::groupName("combustionProperties", phaseName),
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+    turbulencePtr_(),
+    mesh_(mesh),
+    active_(lookupOrDefault<Switch>("active", true)),
+    coeffs_(subDict(modelType + "Coeffs")),
+    modelType_(modelType),
+    phaseName_(phaseName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::combustionModel::~combustionModel()
+{
+    if (turbulencePtr_)
+    {
+        turbulencePtr_ = 0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+bool Foam::combustionModel::read()
+{
+    if (regIOobject::read())
+    {
+        this->lookup("active") >> active_;
+        coeffs_ = subDict(modelType_ + "Coeffs");
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::combustionModel::Sh() const
+{
+    return tmp<Foam::volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("Sh", phaseName_),
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
+        )
+    );
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/combustionModel/combustionModel.H

@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModel
+
+Description
+    Base class for combustion models
+
+SourceFiles
+    combustionModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef combustionModel_H
+#define combustionModel_H
+
+#include "IOdictionary.H"
+#include "turbulentFluidThermoModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class combustionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class combustionModel
+:
+    public IOdictionary
+{
+    // Private Member Functions
+
+        //- Disallow copy construct
+        combustionModel(const combustionModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const combustionModel&);
+
+
+protected:
+
+    // Protected data
+
+        //- Reference to the turbulence model
+        compressibleTurbulenceModel* turbulencePtr_;
+
+        //- Reference to the mesh database
+        const fvMesh& mesh_;
+
+        //- Active
+        Switch active_;
+
+        //- Dictionary of the model
+        dictionary coeffs_;
+
+        //- Model type
+        const word modelType_;
+
+        //- Phase name
+        const word phaseName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("combustionModel");
+
+
+    // Constructors
+
+        //- Construct from components
+        combustionModel
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName=word::null
+        );
+
+
+    //- Destructor
+    virtual ~combustionModel();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to the mesh database
+            inline const fvMesh& mesh() const;
+
+            //- Return const access to phi
+            inline const surfaceScalarField& phi() const;
+
+            //- Return const access to rho
+            virtual tmp<volScalarField> rho() const = 0;
+
+            //- Return access to turbulence
+            inline const compressibleTurbulenceModel& turbulence() const;
+
+            //- Set turbulence
+            inline void setTurbulence
+            (
+                compressibleTurbulenceModel& turbModel
+            );
+
+            //- Is combustion active?
+            inline const Switch& active() const;
+
+            //- Return const dictionary of the model
+            inline const dictionary& coeffs() const;
+
+
+    // Evolution
+
+        //- Correct combustion rate
+        virtual void correct() = 0;
+
+        //- Fuel consumption rate matrix, i.e. source term for fuel equation
+        virtual tmp<fvScalarMatrix> R(volScalarField& Y) const = 0;
+
+        //- Heat release rate calculated from fuel consumption rate matrix
+        virtual tmp<volScalarField> dQ() const = 0;
+
+        //-  Return source for enthalpy equation [kg/m/s3]
+        virtual tmp<volScalarField> Sh() const;
+
+
+    // IO
+
+        //- Update properties from given dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "combustionModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/combustionModel/combustionModelI.H

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::fvMesh& Foam::combustionModel::mesh() const
+{
+    return mesh_;
+}
+
+
+inline const Foam::surfaceScalarField& Foam::combustionModel::phi() const
+{
+    if (turbulencePtr_)
+    {
+        return turbulencePtr_->phi();
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "turbulencePtr_ is empty. Please use "
+            << "combustionModel::setTurbulence "
+            << "(compressibleTurbulenceModel& )"
+            << abort(FatalError);
+
+        return turbulencePtr_->phi();
+    }
+}
+
+
+inline const Foam::compressibleTurbulenceModel&
+Foam::combustionModel::turbulence() const
+{
+    if (turbulencePtr_)
+    {
+        return *turbulencePtr_;
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "turbulencePtr_ is empty. Please use "
+            << "combustionModel::setTurbulence "
+            << "(compressibleTurbulenceModel& )"
+            << abort(FatalError);
+
+        return *turbulencePtr_;
+    }
+}
+
+
+inline const Foam::Switch& Foam::combustionModel::active() const
+{
+    return active_;
+}
+
+
+inline void Foam::combustionModel::setTurbulence
+(
+    compressibleTurbulenceModel& turbModel
+)
+{
+    turbulencePtr_ = &turbModel;
+}
+
+
+inline const Foam::dictionary& Foam::combustionModel::coeffs() const
+{
+    return coeffs_;
+}
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/combustionModel/makeCombustionTypes.H

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeCombustionTypes_H
+#define makeCombustionTypes_H
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeCombustionTypesThermo(CombModel, CombType, Thermo, Table)          \
+                                                                               \
+    typedef Foam::combustionModels::CombModel                                  \
+        <Foam::combustionModels::CombType, Foam::Thermo>                       \
+        CombModel##CombType##Thermo;                                           \
+                                                                               \
+    defineTemplateTypeNameAndDebugWithName                                     \
+    (                                                                          \
+        CombModel##CombType##Thermo,                                           \
+        #CombModel"<"#CombType","#Thermo">",                                   \
+        0                                                                      \
+    );                                                                         \
+                                                                               \
+    namespace Foam                                                             \
+    {                                                                          \
+        namespace combustionModels                                             \
+        {                                                                      \
+            typedef CombModel<CombType, Thermo> CombModel##CombType##Thermo;   \
+            addToRunTimeSelectionTable                                         \
+            (                                                                  \
+                Table,                                                         \
+                CombModel##CombType##Thermo,                                   \
+                dictionary                                                     \
+            );                                                                 \
+        }                                                                      \
+    }
+
+
+#define makeCombustionTypes(CombModel, CombType, Table)                        \
+                                                                               \
+    typedef Foam::combustionModels::CombModel                                  \
+        <Foam::combustionModels::CombType>                                     \
+        CombModel##CombType;                                                   \
+                                                                               \
+    defineTemplateTypeNameAndDebugWithName                                     \
+    (                                                                          \
+        CombModel##CombType,                                                   \
+        #CombModel"<"#CombType">",                                             \
+        0                                                                      \
+    );                                                                         \
+                                                                               \
+    namespace Foam                                                             \
+    {                                                                          \
+        namespace combustionModels                                             \
+        {                                                                      \
+            typedef CombModel<CombType> CombModel##CombType;                   \
+                                                                               \
+            addToRunTimeSelectionTable                                         \
+            (                                                                  \
+                Table,                                                         \
+                CombModel##CombType,                                           \
+                dictionary                                                     \
+            );                                                                 \
+        }                                                                      \
+    }
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/diffusion/diffusion.C

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "diffusion.H"
+#include "fvcGrad.H"
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+diffusion<CombThermoType, ThermoType>::diffusion
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    singleStepCombustion<CombThermoType, ThermoType>
+    (
+        modelType,
+        mesh,
+        phaseName
+    ),
+    C_(readScalar(this->coeffs().lookup("C"))),
+    oxidantName_(this->coeffs().template lookupOrDefault<word>("oxidant", "O2"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+diffusion<CombThermoType, ThermoType>::~diffusion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+void diffusion<CombThermoType, ThermoType>::correct()
+{
+    this->wFuel_ ==
+        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);
+
+    if (this->active())
+    {
+        this->singleMixturePtr_->fresCorrect();
+
+        const label fuelI = this->singleMixturePtr_->fuelIndex();
+
+        const volScalarField& YFuel =
+            this->thermoPtr_->composition().Y()[fuelI];
+
+        if (this->thermoPtr_->composition().contains(oxidantName_))
+        {
+            const volScalarField& YO2 =
+                this->thermoPtr_->composition().Y(oxidantName_);
+
+            this->wFuel_ ==
+                C_*this->turbulence().muEff()
+               *mag(fvc::grad(YFuel) & fvc::grad(YO2))
+               *pos(YFuel)*pos(YO2);
+        }
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+bool diffusion<CombThermoType, ThermoType>::read()
+{
+    if (singleStepCombustion<CombThermoType, ThermoType>::read())
+    {
+        this->coeffs().lookup("C") >> C_ ;
+        this->coeffs().readIfPresent("oxidant", oxidantName_);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs/combustionModels_POSTECH/diffusion/diffusion.H

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::diffusion
+
+Description
+    Simple diffusion-based combustion model based on the principle mixed is
+    burnt. Additional parameter C is used to distribute the heat release rate
+    in time.
+
+SourceFiles
+    diffusion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef diffusion_H
+#define diffusion_H
+
+#include "singleStepCombustion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class diffusion Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType, class ThermoType>
+class diffusion
+:
+    public singleStepCombustion<CombThermoType, ThermoType>
+{
+    // Private data
+
+        //- Model constant
+        scalar C_;
+
+        //- Name of oxidant - default is "O2"
+        word oxidantName_;
+
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        diffusion(const diffusion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const diffusion&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("diffusion");
+
+
+    // Constructors
+
+        //- Construct from components
+        diffusion
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~diffusion();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+
+        // IO
+
+            //- Update properties
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "diffusion.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/diffusion/diffusions.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "thermoPhysicsTypes.H"
+#include "psiThermoCombustion.H"
+#include "rhoThermoCombustion.H"
+#include "diffusion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Combustion models based on sensibleEnthalpy
+makeCombustionTypesThermo
+(
+    diffusion,
+    psiThermoCombustion,
+    gasHThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    diffusion,
+    psiThermoCombustion,
+    constGasHThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    diffusion,
+    rhoThermoCombustion,
+    gasHThermoPhysics,
+    rhoCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    diffusion,
+    rhoThermoCombustion,
+    constGasHThermoPhysics,
+    rhoCombustionModel
+);
+
+ // Combustion models based on sensibleInternalEnergy
+
+makeCombustionTypesThermo
+(
+    diffusion,
+    psiThermoCombustion,
+    gasEThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    diffusion,
+    psiThermoCombustion,
+    constGasEThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    diffusion,
+    rhoThermoCombustion,
+    gasEThermoPhysics,
+    rhoCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    diffusion,
+    rhoThermoCombustion,
+    constGasEThermoPhysics,
+    rhoCombustionModel
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs/combustionModels_POSTECH/infinitelyFastChemistry/infinitelyFastChemistry.C

@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "infinitelyFastChemistry.H"
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+infinitelyFastChemistry<CombThermoType, ThermoType>::infinitelyFastChemistry
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    singleStepCombustion<CombThermoType, ThermoType>
+    (
+        modelType,
+        mesh,
+        phaseName
+    ),
+    C_(readScalar(this->coeffs().lookup("C")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+infinitelyFastChemistry<CombThermoType, ThermoType>::~infinitelyFastChemistry()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+void infinitelyFastChemistry<CombThermoType, ThermoType>::correct()
+{
+    this->wFuel_ ==
+        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);
+
+    if (this->active())
+    {
+        this->singleMixturePtr_->fresCorrect();
+
+        const label fuelI = this->singleMixturePtr_->fuelIndex();
+
+        const volScalarField& YFuel =
+            this->thermoPtr_->composition().Y()[fuelI];
+
+        const dimensionedScalar s = this->singleMixturePtr_->s();
+
+        if (this->thermoPtr_->composition().contains("O2"))
+        {
+            const volScalarField& YO2 = this->thermoPtr_->composition().Y("O2");
+
+            this->wFuel_ ==
+                this->rho()/(this->mesh().time().deltaT()*C_)
+               *min(YFuel, YO2/s.value());
+        }
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+bool infinitelyFastChemistry<CombThermoType, ThermoType>::read()
+{
+    if (singleStepCombustion<CombThermoType, ThermoType>::read())
+    {
+        this->coeffs().lookup("C") >> C_ ;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs/combustionModels_POSTECH/infinitelyFastChemistry/infinitelyFastChemistry.H

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::infinitelyFastChemistry
+
+Description
+    Simple infinitely fast chemistry combustion model based on the principle
+    mixed is burnt. Additional parameter C is used to distribute the heat
+    release rate.in time
+
+SourceFiles
+    infinitelyFastChemistry.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef infinitelyFastChemistry_H
+#define infinitelyFastChemistry_H
+
+#include "singleStepCombustion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class infinitelyFastChemistry Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType, class ThermoType>
+class infinitelyFastChemistry
+:
+    public singleStepCombustion<CombThermoType, ThermoType>
+{
+    // Private data
+
+        //- Model constant
+        scalar C_;
+
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        infinitelyFastChemistry(const infinitelyFastChemistry&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const infinitelyFastChemistry&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("infinitelyFastChemistry");
+
+
+    // Constructors
+
+        //- Construct from components
+        infinitelyFastChemistry
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~infinitelyFastChemistry();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+
+        // IO
+
+            //- Update properties
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "infinitelyFastChemistry.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/infinitelyFastChemistry/infinitelyFastChemistrys.C

@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "thermoPhysicsTypes.H"
+#include "psiThermoCombustion.H"
+#include "rhoThermoCombustion.H"
+#include "infinitelyFastChemistry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Combustion models based on sensibleEnthalpy
+
+makeCombustionTypesThermo
+(
+    infinitelyFastChemistry,
+    psiThermoCombustion,
+    gasHThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    infinitelyFastChemistry,
+    psiThermoCombustion,
+    constGasHThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    infinitelyFastChemistry,
+    rhoThermoCombustion,
+    gasHThermoPhysics,
+    rhoCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    infinitelyFastChemistry,
+    rhoThermoCombustion,
+    constGasHThermoPhysics,
+    rhoCombustionModel
+);
+
+// Combustion models based on sensibleInternalEnergy
+
+makeCombustionTypesThermo
+(
+    infinitelyFastChemistry,
+    psiThermoCombustion,
+    gasEThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    infinitelyFastChemistry,
+    psiThermoCombustion,
+    constGasEThermoPhysics,
+    psiCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    infinitelyFastChemistry,
+    rhoThermoCombustion,
+    gasEThermoPhysics,
+    rhoCombustionModel
+);
+
+makeCombustionTypesThermo
+(
+    infinitelyFastChemistry,
+    rhoThermoCombustion,
+    constGasEThermoPhysics,
+    rhoCombustionModel
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs/combustionModels_POSTECH/laminar/laminar.C

@@ -0,0 +1,214 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "laminar.H"
+#include "fvmSup.H"
+#include "localEulerDdtScheme.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::combustionModels::laminar<Type>::laminar
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    Type(modelType, mesh, phaseName),
+    integrateReactionRate_
+    (
+        this->coeffs().lookupOrDefault("integrateReactionRate", true)
+    )
+{
+    if (integrateReactionRate_)
+    {
+        Info<< "    using integrated reaction rate" << endl;
+    }
+    else
+    {
+        Info<< "    using instantaneous reaction rate" << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::combustionModels::laminar<Type>::~laminar()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::laminar<Type>::tc() const
+{
+    return this->chemistryPtr_->tc();
+}
+
+
+template<class Type>
+void Foam::combustionModels::laminar<Type>::correct()
+{
+    if (this->active())
+    {
+        if (integrateReactionRate_)
+        {
+            if (fv::localEulerDdt::enabled(this->mesh()))
+            {
+                const scalarField& rDeltaT =
+                    fv::localEulerDdt::localRDeltaT(this->mesh());
+
+                if (this->coeffs().found("maxIntegrationTime"))
+                {
+                    scalar maxIntegrationTime
+                    (
+                        readScalar(this->coeffs().lookup("maxIntegrationTime"))
+                    );
+
+                    this->chemistryPtr_->solve
+                    (
+                        min(1.0/rDeltaT, maxIntegrationTime)()
+                    );
+                }
+                else
+                {
+                    this->chemistryPtr_->solve((1.0/rDeltaT)());
+                }
+            }
+            else
+            {
+                this->chemistryPtr_->solve(this->mesh().time().deltaTValue());
+            }
+        }
+        else
+        {
+            this->chemistryPtr_->calculate();
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::combustionModels::laminar<Type>::R(volScalarField& Y) const
+{
+    tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
+
+    fvScalarMatrix& Su = tSu.ref();
+
+    if (this->active())
+    {
+        const label specieI =
+            this->thermo().composition().species()[Y.member()];
+
+        Su += this->chemistryPtr_->RR(specieI);
+    }
+
+    return tSu;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::laminar<Type>::dQ() const
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName(typeName + ":dQ", this->phaseName_),
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+        )
+    );
+
+    if (this->active())
+    {
+        tdQ.ref() = this->chemistryPtr_->dQ();
+    }
+
+    return tdQ;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::laminar<Type>::Sh() const
+{
+    tmp<volScalarField> tSh
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName(typeName + ":Sh", this->phaseName_),
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
+        )
+    );
+
+    if (this->active())
+    {
+        tSh.ref() = this->chemistryPtr_->Sh();
+    }
+
+    return tSh;
+}
+
+
+template<class Type>
+bool Foam::combustionModels::laminar<Type>::read()
+{
+    if (Type::read())
+    {
+        this->coeffs().lookup("integrateReactionRate")
+            >> integrateReactionRate_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/laminar/laminar.H

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::laminar
+
+Description
+    Laminar combustion model.
+
+SourceFiles
+    laminar.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef laminar_H
+#define laminar_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class laminar Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class laminar
+:
+    public Type
+{
+    // Private data
+
+        //- Integrate reaction rate over the time-step
+        //  using the selected ODE solver
+        bool integrateReactionRate_;
+
+protected:
+
+    // Protected Member Functions
+
+        //- Return the chemical time scale
+        tmp<volScalarField> tc() const;
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        laminar(const laminar&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const laminar&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("laminar");
+
+
+    // Constructors
+
+        //- Construct from components
+        laminar
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~laminar();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+            //- Fuel consumption rate matrix.
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
+
+            //- Heat release rate calculated from fuel consumption rate matrix
+            virtual tmp<volScalarField> dQ() const;
+
+            //-  Return source for enthalpy equation [kg/m/s3]
+            virtual tmp<volScalarField> Sh() const;
+
+
+    // IO
+
+            //- Update properties from given dictionary
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "laminar.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/laminar/laminars.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "psiChemistryCombustion.H"
+#include "rhoChemistryCombustion.H"
+#include "laminar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeCombustionTypes(laminar, psiChemistryCombustion, psiCombustionModel);
+makeCombustionTypes(laminar, rhoChemistryCombustion, rhoCombustionModel);
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/noCombustion/noCombustion.C

@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "noCombustion.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType>
+Foam::combustionModels::noCombustion<CombThermoType>::noCombustion
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    CombThermoType(modelType, mesh, phaseName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CombThermoType>
+Foam::combustionModels::noCombustion<CombThermoType>::~noCombustion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType>
+void Foam::combustionModels::noCombustion<CombThermoType>::correct()
+{
+//  Do Nothing
+}
+
+
+template<class CombThermoType>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::combustionModels::noCombustion<CombThermoType>::R
+(
+    volScalarField& Y
+) const
+{
+    tmp<fvScalarMatrix> tSu
+    (
+        new fvScalarMatrix(Y, dimMass/dimTime)
+    );
+
+    return tSu;
+}
+
+
+template<class CombThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::noCombustion<CombThermoType>::dQ() const
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("dQ", this->phaseName_),
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+        )
+    );
+
+    return tdQ;
+}
+
+
+template<class CombThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::noCombustion<CombThermoType>::Sh() const
+{
+    tmp<volScalarField> tSh
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("Sh", this->phaseName_),
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
+        )
+    );
+
+    return tSh;
+}
+
+
+template<class CombThermoType>
+bool Foam::combustionModels::noCombustion<CombThermoType>::read()
+{
+    if (CombThermoType::read())
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/noCombustion/noCombustion.H

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::noCombustion
+
+Description
+    Dummy combustion model for 'no combustion'
+
+SourceFiles
+    noCombustion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noCombustion_H
+#define noCombustion_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class noCombustion Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType>
+class noCombustion
+:
+    public CombThermoType
+{
+
+    //- Disallow copy construct
+    noCombustion(const noCombustion&);
+
+    //- Disallow default bitwise assignment
+    void operator=(const noCombustion&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("noCombustion");
+
+
+    // Constructors
+
+        //- Construct from components
+        noCombustion
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~noCombustion();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+            //- Fuel consumption rate matrix
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
+
+            //- Heat release rate calculated from fuel consumption rate matrix
+            virtual tmp<volScalarField> dQ() const;
+
+            //-  Return source for enthalpy equation [kg/m/s3]
+            virtual tmp<volScalarField> Sh() const;
+
+    // IO
+
+            //- Update properties from given dictionary
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "noCombustion.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/noCombustion/noCombustions.C

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "psiCombustionModel.H"
+#include "rhoCombustionModel.H"
+#include "psiThermoCombustion.H"
+#include "rhoThermoCombustion.H"
+
+#include "noCombustion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeCombustionTypes
+(
+    noCombustion,
+    psiThermoCombustion,
+    psiCombustionModel
+);
+
+makeCombustionTypes
+(
+    noCombustion,
+    rhoThermoCombustion,
+    rhoCombustionModel
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs/combustionModels_POSTECH/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.C

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "psiChemistryCombustion.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::combustionModels::psiChemistryCombustion::psiChemistryCombustion
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    psiCombustionModel(modelType, mesh, phaseName),
+    chemistryPtr_(psiChemistryModel::New(mesh, phaseName))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::psiChemistryCombustion::~psiChemistryCombustion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::psiChemistryModel>
+Foam::combustionModels::psiChemistryCombustion::chem()
+{
+    return chemistryPtr_;
+}
+
+Foam::psiReactionThermo&
+Foam::combustionModels::psiChemistryCombustion::thermo()
+{
+    return chemistryPtr_->thermo();
+}
+
+
+const Foam::psiReactionThermo&
+Foam::combustionModels::psiChemistryCombustion::thermo() const
+{
+    return chemistryPtr_->thermo();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::psiChemistryCombustion::rho() const
+{
+    return chemistryPtr_->thermo().rho();
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.H

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::psiChemistryCombustion
+
+Description
+    Compressibility-based chemistry model wrapper for combustion models
+
+SourceFiles
+    psiChemistryCombustion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef psiChemistryCombustion_H
+#define psiChemistryCombustion_H
+
+#include "autoPtr.H"
+#include "psiCombustionModel.H"
+#include "psiChemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    class psiChemistryCombustion Declaration
+\*---------------------------------------------------------------------------*/
+
+class psiChemistryCombustion
+:
+    public psiCombustionModel
+{
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        psiChemistryCombustion(const psiChemistryCombustion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const psiChemistryCombustion&);
+
+
+protected:
+
+    // Protected data
+
+        //- Pointer to chemistry model
+        autoPtr<psiChemistryModel> chemistryPtr_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components and thermo
+        psiChemistryCombustion
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~psiChemistryCombustion();
+
+
+    // Member Functions
+        autoPtr<psiChemistryModel> chem(); 
+        //- Return access to the thermo package
+        virtual psiReactionThermo& thermo();
+
+        //- Return const access to the thermo package
+        virtual const psiReactionThermo& thermo() const;
+
+        //- Return const access to the density field
+        virtual tmp<volScalarField> rho() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/psiCombustionModel/psiCombustionModel/psiCombustionModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,49 +23,50 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "FlameStructure.H"
+#include "psiCombustionModel.H"
-#include "IOstreams.H"
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+namespace Foam
+{
+namespace combustionModels
+{
+    defineTypeNameAndDebug(psiCombustionModel, 0);
+    defineRunTimeSelectionTable(psiCombustionModel, dictionary);
+}
+}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::CMC::FlameStructure::FlameStructure(Istream& is)
+Foam::combustionModels::psiCombustionModel::psiCombustionModel
-    /*
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
 :
-    base1(is),
+    combustionModel(modelType, mesh, phaseName)
-    base2(is),
+{}
-    member1(is),
+
-    member2(is)
+
-    */
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::psiCombustionModel::~psiCombustionModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::combustionModels::psiCombustionModel::read()
 {
-    // Check state of Istream
+    if (combustionModel::read())
-    is.check("Foam::CMC::FlameStructure::FlameStructure(Foam::Istream&)");
+    {
-}
+        return true;
-
+    }
-
+    else
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+    {
-
+        return false;
-Foam::Istream& Foam::operator>>(Istream& is, CMC::FlameStructure&)
+    }
-{
-    // Check state of Istream
-    is.check
-    (
-        "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::CMC::FlameStructure&)"
-    );
-
-    return is;
-}
-
-
-Foam::Ostream& Foam::operator<<(Ostream& os, const CMC::FlameStructure&)
-{
-    // Check state of Ostream
-    os.check
-    (
-        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
-        "const Foam::CMC::FlameStructure&)"
-    );
-
-    return os;
 }
 
 

libs/combustionModels_POSTECH/psiCombustionModel/psiCombustionModel/psiCombustionModel.H

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::psiCombustionModel
+
+Description
+    Combustion models for compressibility-based thermodynamics
+
+SourceFiles
+    psiCombustionModelI.H
+    psiCombustionModel.C
+    psiCombustionModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef psiCombustionModel_H
+#define psiCombustionModel_H
+
+#include "combustionModel.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "psiReactionThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     class psiCombustionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class psiCombustionModel
+:
+    public combustionModel
+{
+
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        psiCombustionModel(const psiCombustionModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const psiCombustionModel&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("psiCombustionModel");
+
+
+    //- Declare run-time constructor selection tables
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        psiCombustionModel,
+        dictionary,
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        ),
+        (modelType, mesh, phaseName)
+    );
+
+
+    // Constructors
+
+        //- Construct from components
+        psiCombustionModel
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+
+    //- Selector
+    static autoPtr<psiCombustionModel> New
+    (
+        const fvMesh& mesh,
+        const word& phaseName=word::null
+    );
+
+
+    //- Destructor
+    virtual ~psiCombustionModel();
+
+
+    // Member Functions
+
+        //- Return access to the thermo package
+        virtual psiReactionThermo& thermo() = 0;
+
+        //- Return const access to the thermo package
+        virtual const psiReactionThermo& thermo() const = 0;
+
+        //- Return reaction rate as volScalarField ,karam
+        //- defined at laminar.H and .C
+        virtual tmp<volScalarField> RR(label i) const
+        {
+            return tmp<volScalarField> (NULL);
+        }
+ 
+        //- Return tmp of rho
+        virtual tmp<volScalarField> rho() const = 0;
+
+
+    // IO
+
+        //- Update properties from given dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/psiCombustionModel/psiCombustionModel/psiCombustionModelNew.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "psiCombustionModel.H"
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::combustionModels::psiCombustionModel>
+Foam::combustionModels::psiCombustionModel::New
+(
+    const fvMesh& mesh,
+    const word& phaseName
+)
+{
+    const word combModelName
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                IOobject::groupName("combustionProperties", phaseName),
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).lookup("combustionModel")
+    );
+
+    Info<< "Selecting combustion model " << combModelName << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(combModelName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown psiCombustionModel type "
+            << combModelName << endl << endl
+            << "Valid  combustionModels are : " << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    const label tempOpen = combModelName.find('<');
+
+    const word className = combModelName(0, tempOpen);
+
+    return autoPtr<psiCombustionModel>(cstrIter()(className, mesh, phaseName));
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/psiCombustionModel/psiThermoCombustion/psiThermoCombustion.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,36 +23,49 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include "psiThermoCombustion.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::combustionModels::psiThermoCombustion::psiThermoCombustion
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    psiCombustionModel(modelType, mesh, phaseName),
+    thermoPtr_(psiReactionThermo::New(mesh, phaseName))
+{}
+
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-
+Foam::combustionModels::psiThermoCombustion::~psiThermoCombustion()
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+{}
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::psiReactionThermo&
+Foam::combustionModels::psiThermoCombustion::thermo()
+{
+    return thermoPtr_();
+}
 
 
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+const Foam::psiReactionThermo&
+Foam::combustionModels::psiThermoCombustion::thermo() const
+{
+    return thermoPtr_();
+}
 
 
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+Foam::tmp<Foam::volScalarField>
-
+Foam::combustionModels::psiThermoCombustion::rho() const
-
+{
-// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+    return thermoPtr_->rho();
-
+}
-
-// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
 // ************************************************************************* //

libs/combustionModels_POSTECH/psiCombustionModel/psiThermoCombustion/psiThermoCombustion.H

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::psiThermoCombustion
+
+Description
+    Compressibility-based thermo model wrapper for combustion models
+
+SourceFiles
+    psiThermoCombustion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef psiThermoCombustion_H
+#define psiThermoCombustion_H
+
+#include "autoPtr.H"
+#include "psiCombustionModel.H"
+#include "psiReactionThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    class psiThermoCombustion Declaration
+\*---------------------------------------------------------------------------*/
+
+class psiThermoCombustion
+:
+    public psiCombustionModel
+{
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        psiThermoCombustion(const psiThermoCombustion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const psiThermoCombustion&);
+
+
+protected:
+
+    // Protected data
+
+        //- Pointer to chemistry model
+        autoPtr<psiReactionThermo> thermoPtr_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        psiThermoCombustion
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~psiThermoCombustion();
+
+
+    // Member Functions
+
+        //- Return access to the thermo package
+        virtual psiReactionThermo& thermo();
+
+        //- Return const access to the thermo package
+        virtual const psiReactionThermo& thermo() const;
+
+        //- Return const access to the density field
+        virtual tmp<volScalarField> rho() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rhoChemistryCombustion.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoChemistryCombustion::rhoChemistryCombustion
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    rhoCombustionModel(modelType, mesh, phaseName),
+    chemistryPtr_(rhoChemistryModel::New(mesh, phaseName))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoChemistryCombustion::~rhoChemistryCombustion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+Foam::autoPtr<Foam::rhoChemistryModel>
+Foam::combustionModels::rhoChemistryCombustion::chem()
+{
+    return chemistryPtr_;
+}
+
+Foam::rhoReactionThermo&
+Foam::combustionModels::rhoChemistryCombustion::thermo()
+{
+    return chemistryPtr_->thermo();
+}
+
+
+const Foam::rhoReactionThermo&
+Foam::combustionModels::rhoChemistryCombustion::thermo() const
+{
+    return chemistryPtr_->thermo();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::rhoChemistryCombustion::rho() const
+{
+    return chemistryPtr_->thermo().rho();
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.H

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::rhoChemistryCombustion
+
+Description
+    Density-based chemistry model wrapper for combustion models
+
+SourceFiles
+    rhoChemistryCombustion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rhoChemistryCombustion_H
+#define rhoChemistryCombustion_H
+
+#include "autoPtr.H"
+#include "rhoCombustionModel.H"
+#include "rhoChemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    class rhoChemistryCombustion Declaration
+\*---------------------------------------------------------------------------*/
+
+class rhoChemistryCombustion
+:
+    public rhoCombustionModel
+{
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        rhoChemistryCombustion(const rhoChemistryCombustion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const rhoChemistryCombustion&);
+
+
+protected:
+
+    // Protected data
+
+        //- Pointer to chemistry model
+        autoPtr<rhoChemistryModel> chemistryPtr_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components and thermo
+        rhoChemistryCombustion
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~rhoChemistryCombustion();
+
+
+    // Member Functions
+        autoPtr<rhoChemistryModel> chem(); 
+        //- Return access to the thermo package
+        virtual rhoReactionThermo& thermo();
+
+        //- Return const access to the thermo package
+        virtual const rhoReactionThermo& thermo() const;
+
+        //- Return const access to the density field
+        virtual tmp<volScalarField> rho() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.C

@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rhoCombustionModel.H"
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+namespace Foam
+{
+namespace combustionModels
+{
+    defineTypeNameAndDebug(rhoCombustionModel, 0);
+    defineRunTimeSelectionTable(rhoCombustionModel, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoCombustionModel::rhoCombustionModel
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    combustionModel(modelType, mesh, phaseName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoCombustionModel::~rhoCombustionModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::combustionModels::rhoCombustionModel::read()
+{
+    if (combustionModel::read())
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/rhoCombustionModel/rhoCombustionModel/rhoCombustionModel.H

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::rhoCombustionModel
+
+Description
+    Combustion models for rho-based thermodynamics
+
+SourceFiles
+    rhoCombustionModelI.H
+    rhoCombustionModel.C
+    rhoCombustionModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rhoCombustionModel_H
+#define rhoCombustionModel_H
+
+#include "combustionModel.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "rhoReactionThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     class rhoCombustionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class rhoCombustionModel
+:
+    public combustionModel
+{
+
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        rhoCombustionModel(const rhoCombustionModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const rhoCombustionModel&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("rhoCombustionModel");
+
+
+    //- Declare run-time constructor selection tables
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        rhoCombustionModel,
+        dictionary,
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        ),
+        (modelType, mesh, phaseName)
+    );
+
+
+    // Constructors
+
+
+        //- Construct from components
+        rhoCombustionModel
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+
+    //- Selector
+    static autoPtr<rhoCombustionModel> New
+    (
+        const fvMesh& mesh,
+        const word& phaseName=word::null
+    );
+
+
+    //- Destructor
+    virtual ~rhoCombustionModel();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Access combustion dict
+             inline const dictionary& coeff() const;
+
+
+        //- Return access to the thermo package
+        virtual rhoReactionThermo& thermo() = 0;
+
+        //- Return const access to the thermo package
+        virtual const rhoReactionThermo& thermo() const = 0;
+
+        //- Return tmp of rho
+        virtual tmp<volScalarField> rho() const = 0;
+
+
+     // IO
+
+        //- Update properties from given dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/rhoCombustionModel/rhoCombustionModel/rhoCombustionModelNew.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rhoCombustionModel.H"
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::combustionModels::rhoCombustionModel>
+Foam::combustionModels::rhoCombustionModel::New
+(
+    const fvMesh& mesh,
+    const word& phaseName
+)
+{
+    const word combTypeName
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                IOobject::groupName("combustionProperties", phaseName),
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).lookup("combustionModel")
+    );
+
+    Info<< "Selecting combustion model " << combTypeName << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(combTypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown rhoCombustionModel type "
+            << combTypeName << endl << endl
+            << "Valid  combustionModels are : " << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    const label tempOpen = combTypeName.find('<');
+
+    const word className = combTypeName(0, tempOpen);
+
+    return autoPtr<rhoCombustionModel> (cstrIter()(className, mesh, phaseName));
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.C

@@ -0,0 +1,71 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rhoThermoCombustion.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoThermoCombustion::rhoThermoCombustion
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    rhoCombustionModel(modelType, mesh, phaseName),
+    thermoPtr_(rhoReactionThermo::New(mesh, phaseName))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::rhoThermoCombustion::~rhoThermoCombustion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::rhoReactionThermo&
+Foam::combustionModels::rhoThermoCombustion::thermo()
+{
+    return thermoPtr_();
+}
+
+
+const Foam::rhoReactionThermo&
+Foam::combustionModels::rhoThermoCombustion::thermo() const
+{
+    return thermoPtr_();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::rhoThermoCombustion::rho() const
+{
+    return thermoPtr_().rho();
+}
+
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.H

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::rhoThermoCombustion
+
+Description
+    Density-based thermo model wrapper for combustion models
+
+SourceFiles
+    rhoThermoCombustion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rhoThermoCombustion_H
+#define rhoThermoCombustion_H
+
+#include "autoPtr.H"
+#include "rhoCombustionModel.H"
+#include "rhoChemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    class rhoThermoCombustion Declaration
+\*---------------------------------------------------------------------------*/
+
+class rhoThermoCombustion
+:
+    public rhoCombustionModel
+{
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        rhoThermoCombustion(const rhoThermoCombustion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const rhoThermoCombustion&);
+
+
+protected:
+
+    // Protected data
+
+        //- Pointer to thermo model
+        autoPtr<rhoReactionThermo> thermoPtr_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        rhoThermoCombustion
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~rhoThermoCombustion();
+
+
+    // Member Functions
+
+        //- Return access to the thermo package
+        virtual rhoReactionThermo& thermo();
+
+        //- Return const access to the thermo package
+        virtual const rhoReactionThermo& thermo() const;
+
+        //- Return const access to the density field
+        virtual tmp<volScalarField> rho() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/combustionModels_POSTECH/singleStepCombustion/singleStepCombustion.C

@@ -0,0 +1,189 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "singleStepCombustion.H"
+#include "fvmSup.H"
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    CombThermoType(modelType, mesh, phaseName),
+    singleMixturePtr_(NULL),
+    wFuel_
+    (
+        IOobject
+        (
+            IOobject::groupName("wFuel", phaseName),
+            this->mesh().time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh(),
+        dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+    ),
+    semiImplicit_(readBool(this->coeffs_.lookup("semiImplicit")))
+{
+    if (isA<singleStepReactingMixture<ThermoType>>(this->thermo()))
+    {
+        singleMixturePtr_ =
+            &dynamic_cast<singleStepReactingMixture<ThermoType>&>
+            (
+                this->thermo()
+            );
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Inconsistent thermo package for " << this->type() << " model:\n"
+            << "    " << this->thermo().type() << nl << nl
+            << "Please select a thermo package based on "
+            << "singleStepReactingMixture" << exit(FatalError);
+    }
+
+    if (semiImplicit_)
+    {
+        Info<< "Combustion mode: semi-implicit" << endl;
+    }
+    else
+    {
+        Info<< "Combustion mode: explicit" << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+singleStepCombustion<CombThermoType, ThermoType>::~singleStepCombustion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+tmp<fvScalarMatrix> singleStepCombustion<CombThermoType, ThermoType>::R
+(
+    volScalarField& Y
+) const
+{
+    const label specieI =
+        this->thermoPtr_->composition().species()[Y.member()];
+
+    volScalarField wSpecie
+    (
+        wFuel_*singleMixturePtr_->specieStoichCoeffs()[specieI]
+    );
+
+    if (semiImplicit_)
+    {
+        const label fNorm = singleMixturePtr_->specieProd()[specieI];
+        const volScalarField fres(singleMixturePtr_->fres(specieI));
+        wSpecie /= max(fNorm*(Y - fres), scalar(1e-2));
+
+        return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y);
+    }
+    else
+    {
+        return wSpecie + fvm::Sp(0.0*wSpecie, Y);
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+tmp<volScalarField>
+singleStepCombustion<CombThermoType, ThermoType>::Sh() const
+{
+    const label fuelI = singleMixturePtr_->fuelIndex();
+    volScalarField& YFuel =
+        const_cast<volScalarField&>(this->thermoPtr_->composition().Y(fuelI));
+
+    return -singleMixturePtr_->qFuel()*(R(YFuel) & YFuel);
+}
+
+
+template<class CombThermoType, class ThermoType>
+tmp<volScalarField>
+singleStepCombustion<CombThermoType, ThermoType>::dQ() const
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("dQ", this->phaseName_),
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+        )
+    );
+
+    if (this->active())
+    {
+        volScalarField& dQ = tdQ.ref();
+        dQ.ref() = this->mesh().V()*Sh()();
+    }
+    return tdQ;
+}
+
+
+template<class CombThermoType, class ThermoType>
+bool singleStepCombustion<CombThermoType, ThermoType>::read()
+{
+    if (CombThermoType::read())
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs/combustionModels_POSTECH/singleStepCombustion/singleStepCombustion.H

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::singleStepCombustion
+
+Description
+    Base class for combustion models using singleStepReactingMixture.
+
+SourceFiles
+    singleStepCombustion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef singleStepCombustion_H
+#define singleStepCombustion_H
+
+#include "singleStepReactingMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class singleStepCombustion Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType, class ThermoType>
+class singleStepCombustion
+:
+    public CombThermoType
+{
+    // Private Member Functions
+
+        //- Disallow copy construct
+        singleStepCombustion(const singleStepCombustion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const singleStepCombustion&);
+
+
+protected:
+
+    // Protected data
+
+        //- Pointer to singleStepReactingMixture mixture
+        singleStepReactingMixture<ThermoType>* singleMixturePtr_;
+
+        //- Fuel consumption rate
+        volScalarField wFuel_;
+
+        //- Semi-implicit (true) or explicit (false) treatment
+        bool semiImplicit_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        singleStepCombustion
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~singleStepCombustion();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Fuel consumption rate matrix
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
+
+            //- Heat release rate calculated from fuel consumption rate matrix
+            virtual tmp<volScalarField> dQ() const;
+
+            //- Sensible enthalpy source term
+            virtual tmp<volScalarField> Sh() const;
+
+
+        // IO
+
+            //- Update properties from given dictionary
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "singleStepCombustion.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

libs/thermos/Make/files

@@ -1,3 +0,0 @@
-thermoAdd.C
-
-LIB = $(FOAM_USER_LIBBIN)/libthermoAdd

libs/thermos/Make/options

@@ -1,29 +0,0 @@
-OFP_LIB_SRC = ../
-
-EXE_INC = \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(OFP_LIB_SRC)//chemistryModel_POSTECH/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/functions/Polynomial \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-LIB_LIBS = \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lreactionThermophysicalModels \
-    -lspecie \
-    -lthermophysicalFunctions \
-    -lODE \
-    -lfiniteVolume \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -lsolidSpecie \
-    -lfiniteVolume \
-    -lOpenFOAM

libs/thermos/thermoAdd.C

@@ -1,384 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Typedefs
-    Foam::thermoPhysicsTypes
-
-Description
-    Type definitions for thermo-physics models
-
-\*---------------------------------------------------------------------------*/
-
-#include "specie.H"
-#include "perfectGas.H"
-#include "janafThermo.H"
-
-#include "absoluteEnthalpy.H"
-#include "absoluteInternalEnergy.H"
-#include "thermo.H"
-#include "sutherlandTransport.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    // thermo physics types based on absoluteEnthalpy
-
-    typedef
-    sutherlandTransport
-    <
-        species::thermo
-        <
-            janafThermo
-            <
-                perfectGas<specie>
-            >,
-            absoluteEnthalpy
-        >
-    > gasHaThermoPhysics;
-
-
-    // thermo physics types based on absoluteInternalEnergy
-
-    typedef
-    sutherlandTransport
-    <
-        species::thermo
-        <
-            janafThermo
-            <
-                perfectGas<specie>
-            >,
-            absoluteInternalEnergy
-        >
-    > gasEaThermoPhysics;
-
-} // End namespace Foam
-
-/*---------------------------------------------------------------------------*\
-Typedefs
-    Foam::reactionTypes
-
-Description
-    Type definitions for reactions
-\*---------------------------------------------------------------------------*/
-
-#include "Reaction.H"
-
-namespace Foam
-{
-    // absolute enthalpy based reactions
-    typedef Reaction<gasHaThermoPhysics> gasHaReaction;
-
-    // absolute internal ennergy based reactions
-    typedef Reaction<gasEaThermoPhysics> gasEaReaction;
-
-}
-
-// ************************************************************************* //
-
-#include "makeReactionThermo.H"
-
-#include "chemistryReader.H"
-#include "foamChemistryReader.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Gas chemistry readers based on absoluteEnthalpy
-
-makeChemistryReader(gasHaThermoPhysics);
-
-makeChemistryReaderType(foamChemistryReader, gasHaThermoPhysics);
-
-
-
-// Gas chemistry readers based on absoluteInternalEnergy
-
-makeChemistryReader(gasEaThermoPhysics);
-
-makeChemistryReaderType(foamChemistryReader, gasEaThermoPhysics);
-
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "makeReaction.H"
-
-#include "ArrheniusReactionRate.H"
-#include "infiniteReactionRate.H"
-#include "LandauTellerReactionRate.H"
-#include "thirdBodyArrheniusReactionRate.H"
-
-#include "ChemicallyActivatedReactionRate.H"
-#include "JanevReactionRate.H"
-#include "powerSeriesReactionRate.H"
-
-#include "FallOffReactionRate.H"
-#include "LindemannFallOffFunction.H"
-#include "SRIFallOffFunction.H"
-#include "TroeFallOffFunction.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#define makeReactions(Thermo, Reaction)                                        \
-                                                                               \
-    defineTemplateTypeNameAndDebug(Reaction, 0);                               \
-    defineTemplateRunTimeSelectionTable(Reaction, Istream);                    \
-    defineTemplateRunTimeSelectionTable(Reaction, dictionary);                 \
-                                                                               \
-    makeIRNReactions(Thermo, ArrheniusReactionRate)                            \
-    makeIRNReactions(Thermo, infiniteReactionRate)                             \
-    makeIRNReactions(Thermo, LandauTellerReactionRate)                         \
-    makeIRNReactions(Thermo, thirdBodyArrheniusReactionRate)                   \
-                                                                               \
-    makeIRReactions(Thermo, JanevReactionRate)                                 \
-    makeIRReactions(Thermo, powerSeriesReactionRate)                           \
-                                                                               \
-    makePressureDependentReactions                                             \
-    (                                                                          \
-       Thermo,                                                                 \
-       ArrheniusReactionRate,                                                  \
-       LindemannFallOffFunction                                                \
-    )                                                                          \
-                                                                               \
-    makePressureDependentReactions                                             \
-    (                                                                          \
-       Thermo,                                                                 \
-       ArrheniusReactionRate,                                                  \
-       TroeFallOffFunction                                                     \
-    )                                                                          \
-                                                                               \
-    makePressureDependentReactions                                             \
-    (                                                                          \
-       Thermo,                                                                 \
-       ArrheniusReactionRate,                                                  \
-       SRIFallOffFunction                                                      \
-    )
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    // absolute enthalpy based reactions
-    makeReactions(gasHaThermoPhysics, gasHaReaction)
-    makeReactions(gasEaThermoPhysics, gasEaReaction)
-}
-
-// ************************************************************************* //
-
-#include "makeReactionThermo.H"
-
-#include "psiReactionThermo.H"
-#include "hePsiThermo.H"
-
-#include "specie.H"
-#include "perfectGas.H"
-#include "janafThermo.H"
-#include "absoluteEnthalpy.H"
-#include "thermo.H"
-#include "sutherlandTransport.H"
-
-#include "homogeneousMixture.H"
-#include "inhomogeneousMixture.H"
-#include "veryInhomogeneousMixture.H"
-#include "multiComponentMixture.H"
-#include "reactingMixture.H"
-#include "singleStepReactingMixture.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// sutherlandTransport, janafThermo
-
-makeReactionThermo
-(
-    psiThermo,
-    psiReactionThermo,
-    hePsiThermo,
-    homogeneousMixture,
-    sutherlandTransport,
-    absoluteEnthalpy,
-    janafThermo,
-    perfectGas,
-    specie
-);
-
-makeReactionThermo
-(
-    psiThermo,
-    psiReactionThermo,
-    hePsiThermo,
-    inhomogeneousMixture,
-    sutherlandTransport,
-    absoluteEnthalpy,
-    janafThermo,
-    perfectGas,
-    specie
-);
-
-makeReactionThermo
-(
-    psiThermo,
-    psiReactionThermo,
-    hePsiThermo,
-    veryInhomogeneousMixture,
-    sutherlandTransport,
-    absoluteEnthalpy,
-    janafThermo,
-    perfectGas,
-    specie
-);
-
-
-// Multi-component thermo for absolute enthalpy
-
-makeReactionMixtureThermo
-(
-    psiThermo,
-    psiReactionThermo,
-    hePsiThermo,
-    multiComponentMixture,
-    gasHaThermoPhysics
-);
-
-
-// Multi-component thermo for internal energy
-
-makeReactionMixtureThermo
-(
-    psiThermo,
-    psiReactionThermo,
-    hePsiThermo,
-    multiComponentMixture,
-    gasEaThermoPhysics
-);
-
-
-// Multi-component reaction thermo for absolute enthalpy
-
-makeReactionMixtureThermo
-(
-    psiThermo,
-    psiReactionThermo,
-    hePsiThermo,
-    reactingMixture,
-    gasHaThermoPhysics
-);
-
-makeReactionMixtureThermo
-(
-    psiThermo,
-    psiReactionThermo,
-    hePsiThermo,
-    singleStepReactingMixture,
-    gasHaThermoPhysics
-);
-
-
-// Multi-component reaction thermo for internal energy
-
-makeReactionMixtureThermo
-(
-    psiThermo,
-    psiReactionThermo,
-    hePsiThermo,
-    reactingMixture,
-    gasEaThermoPhysics
-);
-
-makeReactionMixtureThermo
-(
-    psiThermo,
-    psiReactionThermo,
-    hePsiThermo,
-    singleStepReactingMixture,
-    gasEaThermoPhysics
-);
-
-} // End namespace Foam
-
-/*---------------------------------------------------------------------------*\
-InClass
-    Foam::psiChemistryModel
-
-Description
-    Creates chemistry model instances templated on the type of thermodynamics
-\*---------------------------------------------------------------------------*/
-
-#include "makeChemistryModel.H"
-
-#include "psiChemistryModel.H"
-#include "chemistryModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    // Chemistry moldels based on absoluteEnthalpy
-    makeChemistryModel
-    (
-        chemistryModel,
-        psiChemistryModel,
-        gasHaThermoPhysics
-    );
-
-    // Chemistry moldels based on absoluteInternalEnergy
-    makeChemistryModel
-    (
-        chemistryModel,
-        psiChemistryModel,
-        gasEaThermoPhysics
-    );
-}
-
-// ************************************************************************* //
-
-#include "makeChemistrySolverTypes.H"
-
-#include "psiChemistryModel.H"
-#include "rhoChemistryModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    // Chemistry solvers based on absoluteEnthalpy
-    makeChemistrySolverTypes(psiChemistryModel, gasHaThermoPhysics);
-    // makeChemistrySolverTypes(rhoChemistryModel, gasHaThermoPhysics);
-
-    // Chemistry solvers based on absoluteInternalEnergy
-    makeChemistrySolverTypes(psiChemistryModel, gasEaThermoPhysics);
-    // makeChemistrySolverTypes(rhoChemistryModel, gasEaThermoPhysics);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-

solvers_post/LagrangianCMCFoam/CMCintegration.H

@@ -1,37 +1,50 @@
-
+                
         forAll(rho, celli)
         {
-            if(F_total[celli] > 0)
+                if(F_total[celli] > 0)
-            {
-                scalarField Ytemp(Ysize, 0.0);
-
-                for(label i = 0 ; i < Ysize ; i++) //mass fraction integration over an eta-space for species i
                 {
-                    for(label k=0 ; k<=group ; k++)
+			scalarField Ytemp(Ysize, 0.0);
-                    {
+			
-                        for(label j=0 ; j<=etamax ; j++)
+                        for(label i=0 ; i < Ysize ; i++) //reaction rate integration over an eta-space for species i
                         {
-                            pdf[j] = Peta[j][celli];
+        			for(label k=0 ; k<=group ; k++)
-                            f[j] = QiCMC[(k*(etamax+1)+j)*Ysize + i];
+        			{
-                        }
+					for(label j=0 ; j<= etamax ; j++)
-                        Ytemp[i] += Ffrac[k][celli]*integration(deltaftn[celli], MFcut, Neta, pdf, f);
+					{
-                    }
+						pdf[j] = Peta[j][celli];
-                }
+						f[j] = QiCMC[(k*(etamax+1)+j)*Ysize + i];
+					}
+					Ytemp[i] += Ffrac[k][celli]*integration(deltaftn[celli], MFcut, Neta, pdf, f);
+                		}
+        		}
 
-                Switch ignoreOFIncrease = true;
+                        if( Ytemp[o2index] > Y[o2index][celli] || Ytemp[fuelindex] > Y[fuelindex][celli])
-                Switch oxygenIncrease =  Ytemp[o2index] > Y[o2index][celli];
+                        {
-                Switch fuelIncrease = Ytemp[fuelindex] > Y[fuelindex][celli];
+                                SUMSh[celli] = 0;
-                if( ignoreOFIncrease && (oxygenIncrease || fuelIncrease))
+                                //keep Yi value
-                {
+                        }
-                    // Do not update Yi
+                        else
+                        {
+                                scalarField ReactionRate(Ysize, 0.0);
+                                scalar tShtemp(0);
+
+                                for(label i = 0 ; i<Ysize ; i++)
+                                {
+                                        ReactionRate[i] = rho[celli] * (Ytemp[i] - Y[i][celli]) / runTime.deltaT().value();
+                                        tShtemp -= chemistry->calculateShCMC( ReactionRate , i);
+                                }
+                                SUMSh[celli] = tShtemp;
+
+                                for(label i = 0 ; i <Ysize ; i++)
+                                {
+                                        Y[i][celli] = Ytemp[i];
+                                }
+
+                        }
                 }
                 else
                 {
-                    for(label i = 0 ; i < Ysize ; i++)
+                        SUMSh[celli] = 0;
-                    {
-                        Y[i][celli] = Ytemp[i];
-                    }
                 }
-            }
         }
-
+ 

solvers_post/LagrangianCMCFoam/EEqn.H

@@ -6,7 +6,7 @@
         fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
       + fvc::ddt(rho, K) + fvc::div(phi, K)
       + (
-            he.name() == "ea"
+            he.name() == "e"
           ? fvc::div
             (
                 fvc::absolute(phi/fvc::interpolate(rho), U),
@@ -17,7 +17,8 @@
         )
       - fvm::laplacian(1.47*turbulence->mut(), he)
      ==
-        fvOptions(rho, he)
+        SUMSh
+       +fvOptions(rho, he)
     );
 
     EEqn.relax();

solvers_post/LagrangianCMCFoam/Flamegroup.H

@@ -19,20 +19,3 @@
         F_total += F[k];
     }
     
-    forAll(rho , celli)
-    {
-        if(F_total[celli] > 0)
-        {
-            for(label k=0 ; k<=group ; k++)
-            {
-                Ffrac[k][celli] = F[k][celli]/F_total[celli];
-            }
-        }
-        else if(F_total[celli] <= 0)
-        {
-            for(label k=0 ; k<=group ; k++)
-            {
-                Ffrac[k][celli] = 0.0;
-            }
-        }
-    }

solvers_post/LagrangianCMCFoam/InjNewFG.H

@@ -1,7 +1,7 @@
 if(fvc::domainIntegrate(F[group]*rho).value() >= fuelperFG) 
 {
 	NewFGindex = NewFGindex + 1;
-	if(group < (nf-1))
+	if(group < (nf-1.5))
 	{
 		group++;
         	forAll(mf.boundaryFieldRef(), patchi)
@@ -17,7 +17,7 @@ if(fvc::domainIntegrate(F[group]*rho).value() >= fuelperFG)
         	}
 		eventRecordFile<<runTime.value()<<tab<<"New flame group is injected"<<tab<<group<<endl;		
 	}
-	else if(group >= (nf-1))
+	else if(group > (nf-1.5))
 	{
 		//retain group, BC
 		forAll(rho, celli)

solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C

@@ -39,7 +39,7 @@ Contact
 
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "psiChemistryModel.H"
+#include "psiChemistryCombustion.H"
 #include "multivariateScheme.H"
 #include "pisoControl.H"
 #include "fvOptions.H"
@@ -48,9 +48,6 @@ Contact
 #include "IFstream.H"
 #include "OFstream.H"
 #include "Math.H" //Mathematical functions for CMC (AMC, gammaln, TDMA)
-
-#include "LinearMesh.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -76,8 +73,6 @@ int main(int argc, char *argv[])
     #include "startSummary.H" //Make logSummary file for CMC
     #include "readCMCProperties.H" //Read and set fields for CMC calculation
 
-    #include "createLinearMesh.H"
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -92,8 +87,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "QiEqn.H"
-
         #include "rhoEqn.H"
         #include "UEqn.H"
         if(init_start_CMC == true)
@@ -103,6 +96,7 @@ int main(int argc, char *argv[])
                 init_start_CMC = false;
         }
         #include "CMCequation.H" //Solve CMC equations
+        #include "YEqn.H"
         #include "CMCintegration.H" //Conditional field integration and get reaction rate
         #include "EEqn.H"
 

solvers_post/LagrangianCMCFoam/LinearMesh.H

@@ -1,337 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::LinearMesh
-
-Description
-
-SourceFiles
-    LinearMeshI.H
-    LinearMesh.C
-    LinearMeshIO.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef LinearMesh_H
-#define LinearMesh_H
-
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class Istream;
-class Ostream;
-
-// Forward declaration of friend functions and operators
-class LinearMesh;
-Istream& operator>>(Istream&, LinearMesh&);
-Ostream& operator<<(Ostream&, const LinearMesh&);
-
-
-/*---------------------------------------------------------------------------*\
-                         Class LinearMesh Declaration
-\*---------------------------------------------------------------------------*/
-
-class LinearMesh
-{
-    // Private data
-
-        //- Lower end value of 1-D domain
-        scalar x0_;
-
-        //- Upper end value of 1-D domain
-        scalar x1_;
-
-        //- Total number of cells
-        label nCells_;
-
-        //- Total number of cells
-        pointField points_;
-
-        //- Total number of cells
-        faceList faces_;
-
-        //- Total number of cells
-        labelList owner_;
-
-        //- Total number of cells
-        labelList neighbour_;
-
-        //- Total number of cells
-        label nPatches_;
-
-        //- Total number of cells
-        wordList patchNames_;
-
-        //- Total number of cells
-        labelList patchStartIndices_;
-
-        //- Total number of cells
-        labelList patchSizes_;
-
-
-    // Private Member Functions
-
-        //- Disallow Construct null
-        LinearMesh();
-
-        //- Disallow default bitwise copy construct
-        LinearMesh(const LinearMesh&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const LinearMesh&);
-
-        //- Generate points and faces
-        void init(const scalarField &etaValue)
-        {
-            for(label i = 0; i < (nCells_+1); i++)
-            {
-                points_[4*i+0] = vector(etaValue[i], 0.0, 0.0);
-                points_[4*i+1] = vector(etaValue[i], 0.1, 0.0);
-                points_[4*i+2] = vector(etaValue[i], 0.1, 0.1);
-                points_[4*i+3] = vector(etaValue[i], 0.0, 0.1);
-            }
-
-            const label nBoundaryFaces = 2;
-            const label nInternalFaces = nCells_-1;
-            const label nEmptyFaces = 4*nCells_;
-
-            for(label i = 0; i < nInternalFaces; i++)
-            {
-                IStringStream faceFormatStream ("4(4 5 6 7)");
-                labelField faceFormat (faceFormatStream);
-                faces_[i] = face(faceFormat + 4*i);
-            }
-
-            {
-                IStringStream faceFormatStream ("4(0 3 2 1)");
-                faces_[nInternalFaces] = face(faceFormatStream);
-            }
-
-            {
-                IStringStream faceFormatStream ("4(4 5 6 7)");
-                labelField faceFormat (faceFormatStream);
-                faces_[nInternalFaces+1] = face(faceFormat + 4*nInternalFaces);
-            }
-
-            for(label i = 0; i < nCells_; i++)
-            {
-                IStringStream faceFormatStream ("4(4(0 1 5 4)4(1 2 6 5)4(2 3 7 6)4(3 0 4 7))");
-                List<labelField> faceFormats (faceFormatStream);
-
-                label fi = nInternalFaces + nBoundaryFaces + 4*i;
-                faces_[fi+0] = face(faceFormats[0] + 4*i);
-                faces_[fi+1] = face(faceFormats[1] + 4*i);
-                faces_[fi+2] = face(faceFormats[2] + 4*i);
-                faces_[fi+3] = face(faceFormats[3] + 4*i);
-            }
-
-
-            for(label i = 0; i < nInternalFaces; i++)
-            {
-                owner_[i] = i;
-                neighbour_[i] = i+1;
-            }
-
-            {
-                owner_[nInternalFaces] = 0;
-            }
-
-            {
-                owner_[nInternalFaces+1] = nCells_ - 1;
-            }
-
-            for(label i = 0; i < nCells_; i++)
-            {
-                label fi = nInternalFaces + nBoundaryFaces + 4*i;
-                owner_[fi+0] = i;
-                owner_[fi+1] = i;
-                owner_[fi+2] = i;
-                owner_[fi+3] = i;
-            }
-
-
-            nPatches_ = 3;
-
-            IStringStream sInd
-                ( 
-                 word("3(")
-                 + name(nInternalFaces) + " "
-                 + name(nInternalFaces+1) + " "
-                 + name(nInternalFaces+2) + " "
-                 + word(")")
-                );
-            patchStartIndices_ = labelList(sInd);
-
-            IStringStream sSize(word("3(1 1 ")+name(4*nCells_)+word(")"));
-            patchSizes_ = labelList(sSize);
-
-            IStringStream sName("3(lowerEnd upperEnd Sides)");
-            patchNames_ = wordList(sName);
-        }
-
-public:
-
-    // Static data members
-
-        //- Static data staticData
-        // static const dataType staticData;
-
-
-    // Constructors
-
-        //- Construct from components
-        LinearMesh(const scalar x0, const scalar x1, const label nCells)
-            :
-                x0_(x0),
-                x1_(x1),
-                nCells_(nCells),
-                points_(4*(nCells+1)),
-                faces_(5*nCells+1),
-                owner_(faces_.size()),
-                neighbour_(nCells-1)
-        {
-
-            scalarField etaValue(nCells+1, 0.0);
-            scalar delta = (x1 - x0) / scalar(nCells);
-            forAll(etaValue, i)
-            {
-                etaValue[i] = x0 + delta*i;
-            }
-
-
-            init(etaValue);
-
-        }
-
-        //- Construct from components
-        LinearMesh(const scalarField &x)
-            :
-                x0_(x[0]),
-                x1_(x.last()),
-                nCells_(x.size()-1),
-                points_(4*(nCells_+1)),
-                faces_(5*nCells_+1),
-                owner_(faces_.size()),
-                neighbour_(nCells_-1)
-        {
-
-            init(x);
-
-        }
-
-        //- Construct from Istream
-        LinearMesh(Istream&);
-
-        //- Construct as copy
-        // LinearMesh(const LinearMesh&);
-
-
-    // Selectors
-
-        //- Select null constructed
-        // static autoPtr<LinearMesh> New();
-
-
-    //- Destructor
-    ~LinearMesh() {}
-
-
-    // Member Functions
-
-        // Access
-            //- Lower end value of 1-D domain
-            scalar x0() const {return x0_;};
-
-            //- Upper end value of 1-D domain
-            scalar x1() const {return x1_;};
-
-            //- Total number of cells
-            label nCells() const {return nCells_;};
-
-            //- Total number of cells
-            const pointField &points() const {return points_;};
-
-            //- Total number of cells
-            const faceList &faces() const {return faces_;};
-
-            //- Total number of cells
-            const labelList &owner() const {return owner_;};
-
-            //- Total number of cells
-            const labelList &neighbour() const {return neighbour_;};
-
-            //- Total number of cells
-            label nPatches () const {return nPatches_;};
-
-            //- Total number of cells
-            word patchName (const label i) const {return patchNames_[i];};
-
-            //- Total number of cells
-            label patchStartIndex (const label i) const {return patchStartIndices_[i];};
-
-            //- Total number of cells
-            label patchSize (const label i) const {return patchSizes_[i];};
-
-
-        // Check
-
-        // Edit
-
-        // Write
-
-
-    // Member Operators
-
-        // void operator=(const LinearMesh&);
-
-
-    // Friend Functions
-
-    // Friend Operators
-
-    // IOstream Operators
-
-        friend Istream& operator>>(Istream&, LinearMesh&);
-        friend Ostream& operator<<(Ostream&, const LinearMesh&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// #include "LinearMeshI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

solvers_post/LagrangianCMCFoam/Make/options

@@ -11,12 +11,13 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(DEV_PATH)/chemistryModel_POSTECH/lnInclude
+    -I$(DEV_PATH)/chemistryModel_POSTECH/lnInclude \
+    -I$(DEV_PATH)/combustionModels_POSTECH/lnInclude
 
 EXE_LIBS = \
     -L$(FOAM_USER_LIBBIN) \
         -lchemistryModel_POSTECH \
-        -lthermoAdd \
+        -lcombustionModels_POSTECH \
     -lfiniteVolume \
     -lfvOptions \
     -lmeshTools \

solvers_post/LagrangianCMCFoam/QCMC_init.H

@@ -1,94 +1,97 @@
-IFstream slfmFile (SLFMinit);
+		IFstream slfmFile
+		(
+		        SLFMinit
+		);
 
-word garbage, yname;
+		word garbage, yname;
-scalar neta(0), ny(0);
+		scalar neta(0), ny(0);
 
-slfmFile.getLine(garbage);//INPUT
+	        slfmFile.getLine(garbage);//INPUT
-slfmFile.getLine(garbage);//NO.
+	        slfmFile.getLine(garbage);//NO.
-slfmFile>>neta>>ny;//number of eta and species
+	        slfmFile>>neta>>ny;//number of eta and species
 
-slfmFile.getLine(garbage);//blank line
+	        slfmFile.getLine(garbage);//blank line
-slfmFile.getLine(garbage);//PRESSURE
+	        slfmFile.getLine(garbage);//PRESSURE
-slfmFile.getLine(garbage);//1.0
+	        slfmFile.getLine(garbage);//1.0
-slfmFile.getLine(garbage);//MIXTURE
+	        slfmFile.getLine(garbage);//MIXTURE
 
-scalarField eta(neta, 0);
+	        scalarField eta(neta, 0);
-scalarField yi_temp(neta*(ny+1), 0), yi((etamax+1)*(ny+1),0), wi(neta*(ny+1), 0);
+	        scalarField yi_temp(neta*(ny+1), 0), yi((etamax+1)*(ny+1),0), wi(neta*(ny+1), 0);
 
-for(label j=0 ; j<neta ; j++)
+	        for(label j=0 ; j<neta ; j++)
-{
+	        {
-    slfmFile>>eta[j];
+	                slfmFile>>eta[j];
-}
+	        }
 
-slfmFile.getLine(garbage); //blank line
+	        slfmFile.getLine(garbage); //blank line
-slfmFile.getLine(garbage); //INITIAL
+	        slfmFile.getLine(garbage); //INITIAL
 
-for(label i=0; i<ny ; i++) //species
+	        for(label i=0; i<ny ; i++) //species
-{
+	        {
-    slfmFile>>yname;
+	                slfmFile>>yname;
-    Info<<yname<<endl;
+			Info<<yname<<endl;
-    for(label j=0 ; j<neta ; j++)
+			for(label j=0 ; j<neta ; j++)
-    {
+			{
-        slfmFile>>yi_temp[j*(ny+1)+i];
+				slfmFile>>yi_temp[j*(ny+1)+i];
-    }
+			}
 
-    for(label j=0 ; j<neta ; j++)
+			for(label j=0 ; j<neta ; j++)
-    {
+			{
-        slfmFile>>wi[j*(ny+1)+i]; //reaction rates
+				slfmFile>>wi[j*(ny+1)+i]; //reaction rates
-    }
+			}
-    slfmFile.getLine(garbage);//blank line
+			slfmFile.getLine(garbage);//blank line
-}
+		}
-//Temperature
+		//Temperature
 
-slfmFile>>yname;
+		slfmFile>>yname;
-Info<<yname<<endl;
+		Info<<yname<<endl;
-for(label j=0 ; j<neta ; j++)
+		for(label j=0 ; j<etamax+1 ; j++)
-{
+		{
-    slfmFile>>yi_temp[j*(ny+1)+ny];
+			slfmFile>>yi_temp[j*(ny+1)+ny];
 
-}
+		}
 
-Interpolation(etaValue, eta, 1, ny+1, yi_temp, yi);
+		Interpolation(etaValue, eta, 1, ny+1, yi_temp, yi);
 
-for(label j = 0 ; j<=etamax ; j++)
+		for(label j = 0 ; j<=etamax ; j++)
-{
+		{
-    for(label i = 0 ; i<Ysize ; i++)
+			for(label i = 0 ; i<Ysize ; i++)
-    {
+			{
-        QiCMC[j*Ysize+i] = yi[j*(ny+1)+i];
+				QiCMC[j*Ysize+i] = yi[j*(ny+1)+i];
-    }
+			}
-    TCMC[j] = yi[j*(ny+1)+ny];
+			TCMC[j] = yi[j*(ny+1)+ny];
-}
+		}
 
-for(label k=1 ; k<nf ; k++)
+		for(label k=1 ; k<nf ; k++)
-{
+		{
-    for(label j = 0 ; j<=etamax ; j++)
+			for(label j = 0 ; j<=etamax ; j++)
-    {
+			{
-        for(label i = 0 ; i<Ysize ; i++)
+				for(label i = 0 ; i<Ysize ; i++)
-        {
+				{
-            QiCMC[(k*(etamax+1)+j)*Ysize+i] = QiCMC[j*Ysize+i];
+					QiCMC[(k*(etamax+1)+j)*Ysize+i] = QiCMC[j*Ysize+i];
-        }
+				}
-        TCMC[j+k*(etamax+1)] = TCMC[j];
+				TCMC[j+k*(etamax+1)] = TCMC[j];
-    }
+			}
-}
+		}
 
 
 
 for(label j=0 ; j <= etamax ; j++) //get conditional enthalpy 
 {
-    for(label k = 0 ; k<nf ; k++)
+        for(label k = 0 ; k<nf ; k++)
-    {
-        scalar Tin = TCMC[ j+k*(etamax+1) ];
-        scalarField QiCMCin(Ysize, 0);
-
-        for(label i=0 ; i<Ysize ; i++)
         {
-            QiCMCin[i] = QiCMC[(k*(etamax+1)+j)*Ysize + i];
+                scalar Tin = TCMC[ j+k*(etamax+1) ];
+                scalarField QiCMCin(Ysize, 0);
+
+                for(label i=0 ; i<Ysize ; i++)
+                {
+                        QiCMCin[i] = QiCMC[(k*(etamax+1)+j)*Ysize + i];
+                }
+
+                rhoCMC[j+k*(etamax+1)] = chemistry->calculateRHOCMC(QiCMCin , Tin, Pin);
+
+                QhCMC[j+k*(etamax+1)] = chemistry->calculateHCMC(QiCMCin , Tin, rhoCMC[ j+k*(etamax+1)], Pin);
+
         }
-
-        rhoCMC[j+k*(etamax+1)] = chemistry->calculateRHOCMC(QiCMCin , Tin, Pin);
-
-        QhCMC[j+k*(etamax+1)] = chemistry->calculateHCMC(QiCMCin , Tin, rhoCMC[ j+k*(etamax+1)], Pin);
-
-    }
 }
 

solvers_post/LagrangianCMCFoam/QiEqn.H

@@ -1,29 +0,0 @@
-{
-    fvScalarMatrix testSource(testEtaMesh, testEtaMesh.dimensions()*dimVol/dimTime);
-
-    testSource.source() = 1.0;
-    
-    fvScalarMatrix QiEqn
-    (
-        fvm::ddt(testEtaMesh)
-//      + mvConvection->fvmDiv(phi, he)
-//      + fvc::ddt(rho, K) + fvc::div(phi, K)
-//      + (
-//            he.name() == "ea"
-//          ? fvc::div
-//            (
-//                fvc::absolute(phi/fvc::interpolate(rho), U),
-//                p,
-//                "div(phiv,p)"
-//            )
-//          : -dpdt
-//        )
-      - fvm::laplacian(DQi, testEtaMesh)
-      ==
-      testSource
-    );
-
-    QiEqn.relax();
-
-    QiEqn.solve();
-}

solvers_post/LagrangianCMCFoam/createFields.H

@@ -1,10 +1,16 @@
 Info<< "Creating reaction model\n" << endl;
 
-autoPtr<psiChemistryModel> chemistry(psiChemistryModel::New(mesh));
+combustionModels::psiCombustionModel* pRxnCMC = (combustionModels::psiCombustionModel::New(mesh)).ptr();
 
+autoPtr<combustionModels::psiChemistryCombustion> reaction
+(
+    dynamic_cast <combustionModels::psiChemistryCombustion*> (pRxnCMC)
+);
+
+autoPtr<psiChemistryModel> chemistry = reaction->chem();
 psiReactionThermo& thermo = chemistry->thermo();
 
-thermo.validate(args.executable(), "ha");
+thermo.validate(args.executable(), "h", "e");
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -77,6 +83,9 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
+// Set the turbulence into the reaction model
+reaction->setTurbulence(turbulence());
+
 
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt

solvers_post/LagrangianCMCFoam/createLinearMesh.H

@@ -1,145 +0,0 @@
-Info<< "Constructing single cell mesh" << nl << endl;
-
-LinearMesh lMesh(0, 1, 100);
-
-pointField CMCPoints(lMesh.points());
-faceList CMCFaces(lMesh.faces());
-labelList CMCOwners(lMesh.owner());
-labelList CMCNeighbors(lMesh.neighbour());
-
-fvMesh mfMesh
-(
-    IOobject
-    (
-        "mfSpace",
-        runTime.timeName(),
-        runTime,
-        IOobject::READ_IF_PRESENT
-    ),
-    xferMove<Field<vector>>(CMCPoints),
-    CMCFaces.xfer(),
-    CMCOwners.xfer(),
-    CMCNeighbors.xfer()
-);
-
-if (mfMesh.boundaryMesh().empty())
-{
-    List<polyPatch*> mfPatches(3);
-
-    mfPatches[0] = new polyPatch
-    (
-        lMesh.patchName(0),
-        lMesh.patchSize(0),
-        lMesh.patchStartIndex(0),
-        0,
-        mfMesh.boundaryMesh(),
-        polyPatch::typeName
-    );
-
-    mfPatches[1] = new polyPatch
-    (
-        lMesh.patchName(1),
-        lMesh.patchSize(1),
-        lMesh.patchStartIndex(1),
-        1,
-        mfMesh.boundaryMesh(),
-        polyPatch::typeName
-    );
-
-    mfPatches[2] = new emptyPolyPatch
-    (
-        lMesh.patchName(2),
-        lMesh.patchSize(2),
-        lMesh.patchStartIndex(2),
-        2,
-        mfMesh.boundaryMesh(),
-        emptyPolyPatch::typeName
-    );
-
-    mfMesh.addFvPatches(mfPatches);
-}
-
-IStringStream tempStream("3(fixedValue fixedValue empty)");
-wordList testEtaPatchTypeNames(tempStream);
-
-volScalarField testEtaMesh
-(
-    IOobject
-    (
-        "testEtaMesh",
-        runTime.timeName(),
-        mfMesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mfMesh,
-    dimensionedScalar("testEtaMesh", dimless, 1.0),
-    testEtaPatchTypeNames
-);
-
-const volVectorField &cEtaMesh(mfMesh.C());
-forAll(cEtaMesh, celli)
-{
-    const scalar x = cEtaMesh[celli].x();
-    if (x > 0.5)
-    {
-        testEtaMesh[celli] = 1.0;
-    }
-    else
-    {
-        testEtaMesh[celli] = 0.0;
-    }
-    
-}
-
-volScalarField DQi
-(
-    IOobject
-    (
-        "DQi",
-        runTime.timeName(),
-        mfMesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    mfMesh,
-    dimensionedScalar("DQi", dimArea/dimTime, 0.01),
-    testEtaPatchTypeNames
-);
-
-PtrList<PtrList<volScalarField> > FStruct (nf);
-for(label k = 0 ; k < nf ; k++)
-{
-    FStruct.set (k, new PtrList<volScalarField> (Ysize+3));
-
-    PtrList<volScalarField> &Fk(FStruct[k]);
-    Info << Fk.size() << endl;
-    for (label i = 0; i < Ysize + 3; i++)
-    {
-        Fk.set (i,
-            new volScalarField 
-            (
-                IOobject
-                (
-                    "test",
-                    runTime.timeName(),
-                    mfMesh,
-                    IOobject::NO_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                mfMesh,
-                dimensionedScalar("test", dimless, 0.0),
-                testEtaPatchTypeNames
-            )
-        );
-
-    }
-}
-
-for(label k = 0 ; k < nf ; k++)
-{
-    for (label i = 0; i < Ysize + 3; i++)
-    {
-        FStruct[k][i].rename("test."+name(k)+"."+name(i));
-    }
-}

solvers_post/LagrangianCMCFoam/readCMCProperties.H

@@ -81,7 +81,7 @@ word fuelName(CMCdict.lookup("fuel_name"));
 word CMCresults_Time(CMCdict.lookup("CMCresults_Time"));
 scalar previousnf(readScalar(CMCdict.lookup("previousnf")));
 scalar fuelperFG(readScalar(CMCdict.lookup("fuelperFG")));
-label nf(readLabel(CMCdict.lookup("nf")));
+scalar nf(readScalar(CMCdict.lookup("nf")));
 
 scalar etamax( readScalar(CMCdict.lookup("etamax")) );
 scalarField etaValue(etamax+1, 0.0);
@@ -149,7 +149,7 @@ else if(uni_eta == false)
 }
 
 scalarField ChemDeltaT( (etamax+1)*nf, 0.0000001); //initial chemical time step;
-label group(0); //flame group number to add evaporation source (default = 0), Flamegroup.H line# 29
+scalar group(0); //flame group number to add evaporation source (default = 0), Flamegroup.H line# 29
 
 
 //Peta field

solvers_post/LagrangianCMCFoam/setSDRCMC.H

@@ -31,6 +31,20 @@
                                 CSDRCELL[j][celli] = C2 * AMC(eta);
                         }
                 }
+		if(F_total[celli] > 0)
+               	{
+			for(label k=0 ; k<=group ; k++)
+			{
+                		Ffrac[k][celli] = F[k][celli]/F_total[celli];
+			}
+              	}
+		else if(F_total[celli] <= 0)
+		{
+                        for(label k=0 ; k<=group ; k++)
+                        {
+                                Ffrac[k][celli] = 0.0;
+                        }			
+		}
         }
 
 //get density weighted averaged value of conditional SDR for each eta and flame-group 

solvers_post/LagrangianCMCSprayFoam/LagrangianCMCSprayFoam.C

@@ -34,7 +34,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "basicSprayCloud.H"
 //#include "psiCombustionModel.H"
-#include "psiChemistryModel.H"
+#include "psiChemistryCombustion.H" //karam
 #include "radiationModel.H"
 #include "SLGThermo.H"
 //#include "pimpleControl.H"

solvers_post/LagrangianCMCSprayFoam/createFields.H

@@ -2,8 +2,14 @@
 
 Info<< "Creating combustion model\n" << endl;
 
-autoPtr<psiChemistryModel> chemistry(psiChemistryModel::New(mesh));
+combustionModels::psiCombustionModel* pRxnCMC = (combustionModels::psiCombustionModel::New(mesh)).ptr();
 
+autoPtr<combustionModels::psiChemistryCombustion> reaction
+(
+    dynamic_cast <combustionModels::psiChemistryCombustion*> (pRxnCMC)
+);
+
+autoPtr<psiChemistryModel> chemistry = reaction->chem();
 psiReactionThermo& thermo = chemistry->thermo();
 
 thermo.validate(args.executable(), "h", "e");
@@ -88,6 +94,9 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
+// Set the turbulence into the combustion model
+reaction->setTurbulence(turbulence());
+
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (

solvers_post/LagrangianCMCdieselEngineFoam4x/LagrangianCMCdieselEngineFoam.C

@@ -35,7 +35,7 @@ Description
 #include "engineMesh.H"
 #include "turbulentFluidThermoModel.H"
 #include "basicSprayCloud.H"
-#include "psiChemistryModel.H"
+#include "psiChemistryCombustion.H" //karam
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pisoControl.H"

solvers_post/LagrangianCMCdieselEngineFoam4x/Make/options

@@ -25,11 +25,13 @@ EXE_INC = \
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(DEV_PATH)/chemistryModel_POSTECH/lnInclude
+    -I$(DEV_PATH)/chemistryModel_POSTECH/lnInclude \
+    -I$(DEV_PATH)/combustionModels_POSTECH/lnInclude
 
 EXE_LIBS = \
     -L$(FOAM_USER_LIBBIN) \
         -lchemistryModel_POSTECH \
+        -lcombustionModels_POSTECH \
     -lfiniteVolume \
     -lmeshTools \
     -lsampling \

solvers_post/LagrangianCMCdieselEngineFoam4x/createFields.H

@@ -1,8 +1,14 @@
 #include "readGravitationalAcceleration.H"
 Info<< "Creating reaction model\n" << endl;
 
-autoPtr<psiChemistryModel> chemistry(psiChemistryModel::New(mesh));
+combustionModels::psiCombustionModel* pRxnCMC = (combustionModels::psiCombustionModel::New(mesh)).ptr();
 
+autoPtr<combustionModels::psiChemistryCombustion> reaction
+(
+    dynamic_cast <combustionModels::psiChemistryCombustion*> (pRxnCMC)
+);
+
+autoPtr<psiChemistryModel> chemistry = reaction->chem();
 psiReactionThermo& thermo = chemistry->thermo();
 
 thermo.validate(args.executable(), "h", "e");
@@ -87,6 +93,9 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
+// Set the turbulence into the reaction model
+reaction->setTurbulence(turbulence());
+
 
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt

solvers_post/SLFMFoam/BetaPDF.H

@@ -3,27 +3,15 @@ class BetaPDF
 
     //private variable 
 
-    //raw parameter mixture fraction and mf variance
-    scalar mf_, mfVar_;
-
     //beta-pdf parameter alpha and beta
     scalar alpha_, beta_;
 
-    //cutting points
+    //cutting point and number of eta-space
-    scalarField etaCut_;
+    scalarField etaCut_, N_;
-
-    //number of eta-space
-    labelList N_;
 
     //detailed integration space
     scalarField etaSpace_;
 
-    //pdf numerators
-    scalarField pdfNum_;
-
-    //pdf denominator
-    scalar pdfDen_;
-
     //detailed integration space, part
     typedef List<scalarField> scalarFieldArray1d;
     scalarFieldArray1d etaPart_;
@@ -48,6 +36,8 @@ public:
     etaPart_(N_.size())
     {
         //Ref. F.Liu et al., INT. J. THERM. SCI. 41 (2002) 763-772.
+        Info<<"Construct Beta-PDF"<<endl;
+
         scalar del(0);
         label cnt(0);
 
@@ -69,9 +59,6 @@ public:
         }
 
         AMCfine_ = AMC(etaSpace_);
-
-        pdfNum_ = scalarField(etaSpace_.size(), 0.0);
-        pdfDen_ = 1.0;
     }
 
     // Destructor
@@ -80,11 +67,8 @@ public:
     // Member functions
 
     // set beta-pdf parameter alpha_ and beta_
-    void setParameter(const scalar mf, const scalar mfVar)
+    void setParameter(scalar& mf, scalar& mfVar)
     {
-        mf_ = mf;
-        mfVar_ = mfVar;
-
         fdelta_ = false;
         delta_ox = false;
         delta_fu = false;
@@ -108,17 +92,10 @@ public:
             {
                 delta_fu = true;
             }
-
-            limit_ab();
         }
-
-        pdfNum_ = etaFunc(alpha_, beta_, etaSpace_);
-        pdfDen_ = integrate(pdfNum_)
-            + Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL)
-            + Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
     }
     // beta-pdf weighted integration for given mf, mfVar and f
-    scalar evaluate(const scalar mf, const scalar mfVar, const scalarField& etaValue, const scalarField& f)
+    scalar evaluate(scalar& mf, scalar& mfVar, scalarField& etaValue, scalarField& f)
     {
         scalar result(0);
 
@@ -126,6 +103,8 @@ public:
 
         if(fdelta_ == false)
         {
+            limit_ab();
+
             scalar num(0), den(0);
             for(label i=0 ; i<etaCut_.size()-1 ; i++)
             {
@@ -151,45 +130,7 @@ public:
 
         return result;
     }
-
+    scalar value(scalar& mf, scalar& mfVar, scalar& etaValue)
-    // beta-pdf weighted integration for given mf, mfVar and f
-    scalar betaIntegrate(const scalarField& f) const
-    {
-        scalar result = 0.0;
-
-        if(fdelta_)
-        {
-            result = interpolateXY(mf_, etaSpace_, f);
-        }
-        else
-        {
-            scalar num = integrate (f * pdfNum_)
-                + f.first()*Foam::pow(etaSpace_[0], alpha_)/(alpha_ + SMALL)
-                + f.last()*Foam::pow(etaSpace_[0], beta_)/(beta_ + SMALL);
-            result = num/pdfDen_;
-        }
-
-        return result;
-    }
-
-    scalar integrate(const scalarField& f) const
-    {
-        scalar total = 0.0;
-
-        forAll(N_, i)
-        {
-            const labelList::subList prev(N_, i);
-            const label baseI = sum(prev);
-
-            const scalarField::subField subInterval(f, N_[i]+1, baseI);
-
-            total += simps(etaCut_[i], etaCut_[i+1], N_[i], subInterval);
-        }
-
-        return total;
-    }
-
-    scalar value(const scalar mf, const scalar mfVar, const scalar etaValue)
     {
         scalar result(0);
 
@@ -197,6 +138,8 @@ public:
         
         if(fdelta_ == false)
         {
+            limit_ab();
+
             scalar den(0);
             for(label i=0 ; i<etaCut_.size()-1 ; i++)
             {
@@ -229,21 +172,21 @@ public:
             alpha_ = (1.0+fmax*(beta_-2.0))/(1.0-fmax);
         }
     }
-    scalarField etaFunc(const scalar a, const scalar b, const scalarField& eta) const
+    scalarField etaFunc(scalar& a, scalar& b, scalarField& eta)
     {
         return pow(eta, a-1.0)*pow(1.0-eta, b-1.0);
     }
-    scalar etaFunc(const scalar a, const scalar b, const scalar eta) const
+    scalar etaFunc(scalar& a, scalar& b, scalar& eta)
     {
         return Foam::pow(eta, a-1.0)*Foam::pow(1.0-eta, b-1.0);
     }
     //extended Simpson's rule (Numerical recipes, 2nd Ed. p.128)
     //for equally spaced and even N intervals (or odd N+1 points) 
-    scalar simps(const scalar xl, const scalar xh, const label N, const UList<scalar>& fx) const
+    scalar simps(scalar& xl, scalar& xh, scalar& N, scalarField& fx)
     {
         scalar evensum(0.0), oddsum(0.0), sum(0.0);
 
-        scalar h = (xh - xl)/scalar(N);
+        scalar h = (xh - xl)/N;
 
         for(label i=0 ; i<fx.size() ; i++)
         {
@@ -263,7 +206,7 @@ public:
     }
     //Amplitude Mapping Closure (from KIVA)
     //Define exp(-2*(erf^-1(2*eta - 1))^2)
-    scalarField AMC(const scalarField& eta) const
+    scalarField AMC(scalarField& eta)
     {
         const scalar pi = 3.141592;
         const scalar spi = Foam::sqrt(pi);
@@ -285,11 +228,11 @@ public:
 
         return result;
     }
-    scalar AMC(const scalar eta) const
+    scalar AMC(scalar& eta)
     {
         return interpolateXY(eta, etaSpace_, AMCfine_);
     }
-    void C1coeff(const scalar mf, const scalarField& varValue, scalarField& C1table)
+    void C1coeff(scalar& mf, scalarField& varValue, scalarField& C1table)
     {
         scalar maxVar = mf*(1.0-mf);
         scalarField x, fx;

solvers_post/SLFMFoam/FlameStructure/FlameStructure.C

@@ -1,105 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "FlameStructure.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-// const Foam::CMC::dataType Foam::CMC::FlameStructure::staticData();
-
-
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::CMC::FlameStructure::FlameStructure(const dictionary &dict, const label nEta, const label nY)
-:
-    // baseClassName(),
-    NstDir_(dict.lookup("NstFolder")),
-    NstList_(dict.lookup("NstList")),
-    Y_(nY, scalarFieldArray1d(NstList_.size(), scalarField(nEta, 0.0))),
-    W_(nY, scalarFieldArray1d(NstList_.size(), scalarField(nEta, 0.0))),
-    T_(NstList_.size(), scalarField(nEta, 0.0)),
-    Q_(NstList_.size(), scalarField(nEta, 0.0)),
-    h_(NstList_.size(), scalarField(nEta, 0.0)),
-    rho_(NstList_.size(), scalarField(nEta, 0.0)),
-    Rgas_(NstList_.size(), scalarField(nEta, 0.0))
-{}
-
-
-Foam::CMC::FlameStructure::FlameStructure(const FlameStructure&)
-// :
-    // baseClassName(),
-    // data_()
-{}
-
-
-// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
-
-/*
-Foam::autoPtr<Foam::CMC::FlameStructure>
-Foam::CMC::FlameStructure::New()
-{
-    return autoPtr<FlameStructure>(new FlameStructure);
-}
-*/
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::CMC::FlameStructure::~FlameStructure()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
-
-void Foam::CMC::FlameStructure::operator=(const FlameStructure& rhs)
-{
-    // Check for assignment to self
-    if (this == &rhs)
-    {
-        FatalErrorInFunction
-            << "Attempted assignment to self"
-            << abort(FatalError);
-    }
-}
-
-// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
-
-
-// ************************************************************************* //

solvers_post/SLFMFoam/FlameStructure/FlameStructure.H

@@ -1,210 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::CMC::FlameStructure
-
-Description
-
-SourceFiles
-    FlameStructureI.H
-    FlameStructure.C
-    FlameStructureIO.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef FlameStructure_H
-#define FlameStructure_H
-
-#include "error.H"
-#include "scalarField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class Istream;
-class Ostream;
-
-namespace CMC
-{
-class FlameStructure;
-}
-
-// Forward declaration of friend functions and operators
-Istream& operator>>(Istream&, CMC::FlameStructure&);
-Ostream& operator<<(Ostream&, const CMC::FlameStructure&);
-
-namespace CMC
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class FlameStructure Declaration
-\*---------------------------------------------------------------------------*/
-
-class FlameStructure
-{
-    // Private data
-        typedef List<scalarField> scalarFieldArray1d;
-
-        typedef List<scalarFieldArray1d> scalarFieldArray2d;
-
-        typedef List<scalarFieldArray2d> scalarFieldArray3d;
-
-        //- Species mass fraction
-        word NstDir_;
-
-        //- Species mass fraction
-        scalarList NstList_;
-
-        //- Species mass fraction
-        scalarFieldArray2d Y_;
-
-        //- Species production rate
-        scalarFieldArray2d W_;
-
-        //- Temperature
-        scalarFieldArray1d T_;
-
-        //- Heat source
-        scalarFieldArray1d Q_;
-
-        //- Enthalpy
-        scalarFieldArray1d h_;
-
-        //- Density
-        scalarFieldArray1d rho_;
-
-        //- Specific gas constant for mean W
-        scalarFieldArray1d Rgas_;
-
-
-    // Private Member Functions
-
-        //- Disallow Construct null
-        FlameStructure();
-
-        //- Disallow default bitwise copy construct
-        FlameStructure(const FlameStructure&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const FlameStructure&);
-
-        //- limit value of alpha and beta
-        void limitAB();
-
-
-public:
-
-    // Static data members
-
-        //- Static data staticData
-        // static const dataType staticData;
-
-
-    // Constructors
-
-        //- Construct from components
-        FlameStructure(const dictionary &dict, const label nEta, const label nY);
-
-        //- Construct from Istream
-        FlameStructure(Istream&);
-
-        //- Construct as copy
-        // FlameStructure(const FlameStructure&);
-
-
-    // Selectors
-
-        //- Select null constructed
-        // static autoPtr<FlameStructure> New();
-
-
-    //- Destructor
-    ~FlameStructure();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Select null constructed
-            scalarFieldArray2d& Y() {return Y_;}
-
-            //- Select null constructed
-            scalarFieldArray2d& W() {return W_;}
-
-            //- Select null constructed
-            scalarFieldArray1d& T() {return T_;}
-
-            //- Select null constructed
-            scalarFieldArray1d& Q() {return Q_;}
-
-            //- Select null constructed
-            scalarFieldArray1d& rho() {return rho_;}
-
-            //- Select null constructed
-            scalarFieldArray1d& h() {return h_;}
-
-            //- Select null constructed
-            scalarFieldArray1d& Rgas() {return Rgas_;}
-
-        // Check
-
-        // Edit
-
-        // Write
-
-
-    // Member Operators
-
-        // void operator=(const FlameStructure&);
-
-
-    // Friend Functions
-
-    // Friend Operators
-
-    // IOstream Operators
-
-        friend Istream& operator>>(Istream&, FlameStructure&);
-        friend Ostream& operator<<(Ostream&, const FlameStructure&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace CMC
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "FlameStructureI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

solvers_post/SLFMFoam/Make/files

@@ -1,5 +1,3 @@
-FlameStructure/FlameStructure.C
-FlameStructure/FlameStructureIO.C
 SLFMFoam.C
 
 EXE = $(FOAM_USER_APPBIN)/SLFMFoam

solvers_post/SLFMFoam/Make/options

@@ -11,12 +11,13 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(DEV_PATH)/chemistryModel_POSTECH/lnInclude \
-    -I$(LIB_SRC)/combustionModels/lnInclude
+    -I$(DEV_PATH)/combustionModels_POSTECH/lnInclude
 
 EXE_LIBS = \
-    -lchemistryModel \
+    -L$(FOAM_USER_LIBBIN) \
-    -lcombustionModels \
+        -lchemistryModel_POSTECH \
+        -lcombustionModels_POSTECH \
     -lfiniteVolume \
     -lfvOptions \
     -lmeshTools \

solvers_post/SLFMFoam/SLFMFoam.C

@@ -45,8 +45,8 @@ Contact
 #include "rhoChemistryCombustion.H" 
 #include "fvOptions.H"
 #include "simpleControl.H"
-#include "FlameStructure/FlameStructure.H"
 #include "interpolateXY.H" 
+#include "Math.H" //Mathematical functions for CMC (AMC, gammaln, TDMA)
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -68,7 +68,6 @@ int main(int argc, char *argv[])
 
     //SLFM
     #include "readSLFMlib.H" //read SLFM library 
-    Info<<"Construct Beta-PDF"<<endl;
     BetaPDF bpdf(SLFMdict);
     #include "preIntegration.H" //AMC C1 coefficient pdf integration
 
@@ -95,10 +94,16 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
+            #include "updateYi.H" //update species
         if(runTime.write() == true)
         {
-            #include "updateYi.H" //update species
             rho.write();
+            #include "updateYi.H" //update species
+            forAll(postSpecieIndices, yi)
+            {
+                const label y = postSpecieIndices[yi];
+                Y[y].write();
+            }
         }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

solvers_post/SLFMFoam/createFields.H

@@ -1,16 +1,21 @@
 Info<< "\nCreating thermophysical model for Fluid mesh\n" << endl;
 #include "readGravitationalAcceleration.H"
 
-autoPtr<combustionModels::rhoCombustionModel> reaction
+combustionModels::rhoCombustionModel* pRxn = (combustionModels::rhoCombustionModel::New(mesh)).ptr();
+autoPtr<combustionModels::rhoChemistryCombustion> combustion
 (
-    combustionModels::rhoCombustionModel::New(mesh)
+    dynamic_cast <combustionModels::rhoChemistryCombustion*> (pRxn)
 );
-
+autoPtr<rhoChemistryModel> chemistry = combustion->chem();
-rhoReactionThermo& thermo = reaction->thermo();
+rhoReactionThermo& thermo = chemistry->thermo();
 thermo.validate(args.executable(), "h", "e");
 
 basicSpecieMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
+forAll (Y, i)
+{
+    Y[i].writeOpt() = IOobject::NO_WRITE;
+}
 
 word inertSpecie(thermo.lookup("inertSpecie"));
 
@@ -79,8 +84,7 @@ autoPtr<compressible::turbulenceModel> turbulence
 );
 
 // Set the turbulence into the reaction model
-reaction->setTurbulence(turbulence());
+combustion->setTurbulence(turbulence());
-
 
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
@@ -98,7 +102,6 @@ volScalarField dpdt
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
-
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
 forAll(Y, i)
@@ -291,58 +294,33 @@ IOdictionary SLFMdict //SLFM properties dictionary
     )
 );
 
-const label Ysize = Y.size();
+label Ysize = Y.size();
-const wordList postSpecieNames(SLFMdict.lookup("postSpecieNames"));
+wordList postSpecieNames(SLFMdict.lookup("postSpecieNames"));
 labelList postSpecieIndices(postSpecieNames.size(),-1);
 forAll(postSpecieNames, i)
 {
     const word &spName (postSpecieNames[i]);
     const label spI (composition.species()[spName]);
     postSpecieIndices[i] = spI;
+    Y[spI].writeOpt() = IOobject::AUTO_WRITE;
 }
 scalar stoiMF(readScalar(SLFMdict.lookup("stoiMF")));
 scalar stretchFac(readScalar(SLFMdict.lookup("stretchFactor")));
 scalar lowerN(readScalar(SLFMdict.lookup("lowerN")));
 scalar upperN(readScalar(SLFMdict.lookup("upperN")));
 scalar Sc(readScalar(SLFMdict.lookup("Sc")));
-const label etamax(lowerN+upperN+1);
+scalar etamax(lowerN+upperN);
 word outletName(SLFMdict.lookup("outletName"));
 
-labelList nIntervals (SLFMdict.lookup("detailedN"));
+scalarField etaValue(etamax+1,0);
-scalarField etaGrid(sum(nIntervals)+1, 0.0);
+scalarField etaIndex(etamax+1,0);
-
+scalarField deta_h(etamax+1,0);
-{
+scalarField deta_l(etamax+1,0);
-    scalarList etas (SLFMdict.lookup("detailedEta"));
+scalarField deta_m(etamax+1,0);
-
-    forAll(nIntervals, i)
-    {
-        const label itv = nIntervals[i];
-        const labelList::subList prev(nIntervals, i);
-        const label baseI = sum(prev);
-
-        const scalar low = etas[i];
-        const scalar upp = etas[i+1];
-        const scalar delta = (upp - low) / scalar(itv);
-
-        for (label j = 0; j < itv; j++)
-        {
-            etaGrid[baseI+j] = low + delta * j;
-        }
-    }
-    etaGrid.last() = etas.last();
-}
-
-Info << etaGrid << endl;
-
-scalarField etaValue(etamax,0);
-scalarField etaIndex(etamax,0);
-scalarField deta_h(etamax,0);
-scalarField deta_l(etamax,0);
-scalarField deta_m(etamax,0);
 
 #include "setEtaSpace.H"
 
-Info<<"\nNumber of eta-point = "<<etamax<<endl;
+Info<<"\nNumber of eta-point = "<<etamax+1<<endl;
 
 
 scalar NVar(readScalar(SLFMdict.lookup("NVar")));
@@ -353,9 +331,9 @@ scalarField varIndex(NVar+1,0);
 
 
 
-PtrList<volScalarField> Neta(etamax);
+PtrList<volScalarField> Neta(etamax+1);
 
-for(label j=0 ; j<etamax ; j++)
+for(label j=0 ; j<=etamax ; j++)
 {
 
 

solvers_post/SLFMFoam/preIntegration.H

@@ -1,12 +1,12 @@
 scalarFieldArray1d C1table //AMC C1 coefficient table
 (
-    etamax,
+    etamax+1,
     scalarField(NVar+1, 0)
 );
 
 scalarFieldArray2d Ttable
 (
-    etamax,
+    etamax+1,
     scalarFieldArray1d 
     (
         NVar+1,
@@ -16,7 +16,7 @@ scalarFieldArray2d Ttable
 
 scalarFieldArray2d Rtable
 (
-    etamax,
+    etamax+1,
     scalarFieldArray1d
     (
         NVar+1,
@@ -26,10 +26,10 @@ scalarFieldArray2d Rtable
 
 scalarFieldArray3d Ytable
 (
-    Y.size(),
+    postSpecieNames.size(),
     scalarFieldArray2d
-    (
+        (
-        etamax,
+        etamax+1,
         scalarFieldArray1d
         (
             NVar+1,
@@ -39,30 +39,28 @@ scalarFieldArray3d Ytable
 );
 
 Info<<"Construct C1 coefficient table"<<endl;
-for(label j=1 ; j<etamax-1; j++)
+for(label j=1 ; j<etamax; j++)
 {
     bpdf.C1coeff(etaValue[j], varValue, C1table[j]);
 }
 
 Info<<"Construct Temp and Rgas table"<<endl;
-for(label j=0 ; j<etamax ; j++)
+for(label j=0 ; j<=etamax ; j++)
 {
-    const scalar mf = etaValue[j];
+    scalar mf = etaValue[j];
     for(label v=0 ; v<=NVar ; v++)
     {
-        const scalar mfVar = varValue[v]*mf*(1.0-mf);
+        scalar mfVar = varValue[v]*mf*(1.0-mf);
-
-        bpdf.setParameter(mf, mfVar);
-
         for(label n=0 ; n<NstList.size() ; n++)
         {
-            Ttable[j][v][n] = bpdf.betaIntegrate(T_SLFM[n]);
+            Ttable[j][v][n] = bpdf.evaluate(mf, mfVar, etaValue, Y_SLFM[n][Ysize]);
-            Rtable[j][v][n] = bpdf.betaIntegrate(Rgas_SLFM[n]);
+            Rtable[j][v][n] = bpdf.evaluate(mf, mfVar, etaValue, Rgas_SLFM[n]);
-            forAll(Y, yi)
+            forAll(postSpecieIndices, yi)
             {
-                Ytable[yi][j][v][n] = bpdf.betaIntegrate(Y_SLFM[yi][n]);
+                const label y = postSpecieIndices[yi];
+                Ytable[yi][j][v][n] = bpdf.evaluate(mf, mfVar, etaValue, Y_SLFM[n][y]);
             }
         }
     }
-    Info<<scalar(j+1)/etamax*100<<"%"<<endl;
+    Info<<j/etamax*100<<"%"<<endl;
 }

solvers_post/SLFMFoam/readSLFMlib.H

@@ -12,15 +12,39 @@ typedef List<scalarField> scalarFieldArray1d;
 typedef List<scalarFieldArray1d> scalarFieldArray2d;
 typedef List<scalarFieldArray2d> scalarFieldArray3d;
 
-CMC::FlameStructure slfmLibrary(SLFMdict, etaGrid.size(), Y.size());
+scalarFieldArray2d Y_SLFM //species mass frac. [-]
-
+(
-scalarFieldArray2d& Y_SLFM (slfmLibrary.Y());
+    NstList.size(),
-scalarFieldArray2d& W_SLFM (slfmLibrary.W());
+    scalarFieldArray1d
-scalarFieldArray1d& T_SLFM (slfmLibrary.T());
+    (
-scalarFieldArray1d& Q_SLFM (slfmLibrary.Q());
+        Ysize+1,
-scalarFieldArray1d& rho_SLFM (slfmLibrary.rho());
+        scalarField(etamax+1,0)
-scalarFieldArray1d& h_SLFM (slfmLibrary.h());
+    )
-scalarFieldArray1d& Rgas_SLFM (slfmLibrary.Rgas());
+);
+scalarFieldArray2d W_SLFM //reaction rate [1/sec]
+(
+    NstList.size(),
+    scalarFieldArray1d
+    (
+        Ysize+1,
+        scalarField(etamax+1,0)
+    )
+);
+scalarFieldArray1d rho_SLFM //density [g/cm3]
+(
+    NstList.size(),
+    scalarField(etamax+1, 0)
+);
+scalarFieldArray1d h_SLFM //enthalpy [erg/gm]
+(
+    NstList.size(),
+    scalarField(etamax+1, 0)
+);
+scalarFieldArray1d Rgas_SLFM //specific gas constant for meanMW
+(
+    NstList.size(),
+    scalarField(etamax+1, 0)
+);
 
 for(label n=0 ; n<NstList.size() ; n++)
 {
@@ -30,7 +54,7 @@ for(label n=0 ; n<NstList.size() ; n++)
     IFstream fin(runTime.constant()/NstFolder/fname);
 
     string gbg;
-    label etamax_SLFM, NoSpecies;
+    scalar etamax_SLFM, NoSpecies;
 
     //INPUT FILE...
     fin.getLine(gbg);
@@ -40,6 +64,7 @@ for(label n=0 ; n<NstList.size() ; n++)
 
     //Number of eta-point and species
     fin>>etamax_SLFM>>NoSpecies;
+    etamax_SLFM = etamax_SLFM-1;
 
     //blank line
     fin.getLine(gbg);
@@ -54,8 +79,8 @@ for(label n=0 ; n<NstList.size() ; n++)
     fin.getLine(gbg);
 
     //read eta space (from SLFM library)
-    scalarField etaValue_SLFM(etamax_SLFM, 0.0);
+    scalarField etaValue_SLFM(etamax_SLFM+1, 0.0);
-    for(label j=0 ; j<etamax_SLFM ; j++)
+    for(label j=0 ; j<=etamax_SLFM ; j++)
     {
         fin>>etaValue_SLFM[j];
     }
@@ -66,64 +91,30 @@ for(label n=0 ; n<NstList.size() ; n++)
     //INITIAL...
     fin.getLine(gbg);
 
-    //species loop
+    //species and temperature loop
-    scalarField Y_temp(etamax_SLFM, 0.0);
+    scalarField Y_temp(etamax_SLFM+1, 0.0);
-    scalarField W_temp(etamax_SLFM, 0.0);
+    scalarField W_temp(etamax_SLFM+1, 0.0);
-    for(label i=0 ; i<NoSpecies ; i++)
+    for(label i=0 ; i<=NoSpecies ; i++)
     {
         //H2... (species name)
-        string spLine;
+        fin.getLine(gbg);
-        fin.getLine(spLine);
-        IStringStream spStream(spLine);
-        word spName;
-        spStream >> spName;
-        const label spI (composition.species()[spName]);
-        // Info<<spI << spName << endl;
 
         //read species mass fraction (from SLFM library)
-        for(label j=0 ; j<etamax_SLFM ; j++)
+        for(label j=0 ; j<=etamax_SLFM ; j++)
         {
-            scalar Yj = 0.0;
+            fin>>Y_temp[j];
-            fin >> Yj;
-            Y_temp[j] = max(0.0, Yj);
         }
+        Y_temp = max(0.0, Y_temp);
 
-        Y_SLFM[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, Y_temp);
+        Y_SLFM[n][i] = interpolateXY(etaValue, etaValue_SLFM, Y_temp);
 
         //read reaction rate (from SLFM library)
-        for(label j=0 ; j<etamax_SLFM ; j++)
+        for(label j=0 ; j<=etamax_SLFM ; j++)
         {
             fin>>W_temp[j];
         }
 
-        W_SLFM[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, W_temp);
+        W_SLFM[n][i] = interpolateXY(etaValue, etaValue_SLFM, W_temp);
-
-        //blank line
-        fin.getLine(gbg);
-    }
-
-    //read temperature
-    {
-        //H2... (species name)
-        fin.getLine(gbg);
-
-        //read temperature
-        for(label j=0 ; j<etamax_SLFM ; j++)
-        {
-            scalar Yj = 0.0;
-            fin >> Yj;
-            Y_temp[j] = max(0.0, Yj);
-        }
-
-        T_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, Y_temp);
-
-        //read heat source
-        for(label j=0 ; j<etamax_SLFM ; j++)
-        {
-            fin>>W_temp[j];
-        }
-
-        Q_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, W_temp);
 
         //blank line
         fin.getLine(gbg);
@@ -133,14 +124,14 @@ for(label n=0 ; n<NstList.size() ; n++)
     fin.getLine(gbg);
 
     //read density (from SLFM library)
-    scalarField rho_temp(etamax_SLFM, 0.0);
+    scalarField rho_temp(etamax_SLFM+1, 0.0);
-    for(label j=0 ; j<etamax_SLFM ; j++)
+    for(label j=0 ; j<=etamax_SLFM ; j++)
     {
         fin>>rho_temp[j];
     }
     //Info<<rho_temp<<endl;
 
-    rho_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, rho_temp);
+    rho_SLFM[n] = interpolateXY(etaValue, etaValue_SLFM, rho_temp);
 
     //blank line
     fin.getLine(gbg);
@@ -149,17 +140,17 @@ for(label n=0 ; n<NstList.size() ; n++)
     fin.getLine(gbg);
 
     //read enthalpy (from SLFM library)
-    scalarField h_temp(etamax_SLFM, 0.0);
+    scalarField h_temp(etamax_SLFM+1, 0.0);
-    for(label j=0 ; j<etamax_SLFM ; j++)
+    for(label j=0 ; j<=etamax_SLFM ; j++)
     {
         fin>>h_temp[j];
     }
 
-    h_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, h_temp);
+    h_SLFM[n] = interpolateXY(etaValue, etaValue_SLFM, h_temp);
 
-    for(label i=0 ; i<Y.size() ; i++)
+    for(label i=0 ; i<NoSpecies ; i++)
     {
-        Rgas_SLFM[n] += Y_SLFM[i][n]/composition.W(i); // = 1/meanMW
+        Rgas_SLFM[n] += Y_SLFM[n][i]/composition.W(i); // = 1/meanMW
     }
     Rgas_SLFM[n] = 8314.46*Rgas_SLFM[n]; //universal gas constant = 8314.46 [J/K.Kmole]
 

solvers_post/SLFMFoam/setEtaSpace.H

@@ -17,7 +17,7 @@ lowerDelta = stoiMF/lowerSum;
 //set boundary and stoichiometric value of eta
 etaValue[0] = 0.0;
 etaValue[lowerN] = stoiMF;
-etaValue.last() = 1.0;
+etaValue[etamax] = 1.0;
 
 //set internal value of eta (eta > stoi)
 for(label j=0 ; j<upperN-1 ; j++)
@@ -34,7 +34,7 @@ for(label j=0 ; j<lowerN-1 ; j++)
 Info<<etaValue<<endl;
 
 //calculate delta eta
-for(label j=1 ; j<etamax-1 ; j++)
+for(label j=1 ; j<etamax ; j++)
 {
     deta_h[j] = etaValue[j+1] - etaValue[j];
     deta_l[j] = etaValue[j] - etaValue[j-1];
@@ -45,7 +45,7 @@ for(label j=1 ; j<etamax-1 ; j++)
 //Info<<deta_l<<endl;
 //Info<<deta_m<<endl;
 
-for(label j=0 ; j<etamax ; j++)
+for(label j=0 ; j<=etamax ; j++)
 {
     etaIndex[j] = j;
 }

solvers_post/SLFMFoam/setSDR.H

@@ -16,10 +16,10 @@ forAll(mf, cellI)
         jl = 0;
         jh = 1;
     }
-    else if(mf[cellI] > etaValue[etamax-2])
+    else if(mf[cellI] > etaValue[etamax-1])
     {
-        jl = etamax-2;
+        jl = etamax-1;
-        jh = etamax-1;
+        jh = etamax;
     }
     else
     {

solvers_post/SLFMFoam/updateYi.H

@@ -31,8 +31,9 @@ forAll(rho, cellI)
     scalar jfac = jfc[cellI], vfac = vfc[cellI], nfac = nfc[cellI];
 
 
-    forAll(Y, yi)
+    forAll(postSpecieIndices, yi)
     {
+        const label y = postSpecieIndices[yi];
         scalar y00 = Ytable[yi][jl][vl][nl]*(1-jfac)+Ytable[yi][jh][vl][nl]*jfac;
         scalar y01 = Ytable[yi][jl][vl][nh]*(1-jfac)+Ytable[yi][jh][vl][nh]*jfac;
         scalar y10 = Ytable[yi][jl][vh][nl]*(1-jfac)+Ytable[yi][jh][vh][nl]*jfac;
@@ -43,14 +44,15 @@ forAll(rho, cellI)
 
         scalar ygas = y0*(1-nfac)+y1*nfac;
 
-        Y[yi][cellI] = ygas;
+        Y[y][cellI] = ygas;
     }
 }
 
 //correct processor boundary value of T and rho
-forAll(Y, yi)
+forAll(postSpecieIndices, yi)
 {
-    Y[yi].correctBoundaryConditions();
+    const label y = postSpecieIndices[yi];
+    Y[y].correctBoundaryConditions();
 }