postSootFoam, 2 soot models, absolute H/E gas thermos, old EDM

applications/solvers/combustion/postSootFoam/Make/files

@@ -0,0 +1,3 @@
+postSootFoam.C
+
+EXE = $(FOAM_APPBIN)/postSootFoam

applications/solvers/combustion/postSootFoam/Make/options

@@ -0,0 +1,56 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIb_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/properties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/combustionModels/lnInclude
+
+EXE_LIBS = \
+    -lchemistryModel \
+    -lcombustionModels \
+    -lradiationModels \
+    -lthermoAdd \
+    -lliquidProperties_kerosene \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lreactionThermophysicalModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lcompressibleTransportModels \
+    -lliquidProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lthermophysicalFunctions \
+    -lSLGThermo \
+    -lODE \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lfvOptions
+

applications/solvers/combustion/postSootFoam/createClouds.H

@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicReactingMultiphaseCloud parcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    slgThermo
+);

applications/solvers/combustion/postSootFoam/createFieldRefs.H

@@ -0,0 +1 @@
+const volScalarField& psi = thermo.psi();

applications/solvers/combustion/postSootFoam/createFields.H

@@ -0,0 +1,141 @@
+Info<< "\nCreating thermophysical model for Fluid mesh\n" << endl;
+#include "readGravitationalAcceleration.H"
+
+combustionModels::rhoCombustionModel* pRxn = (combustionModels::rhoCombustionModel::New(mesh)).ptr();
+autoPtr<combustionModels::rhoChemistryCombustion> combustion
+(
+    dynamic_cast <combustionModels::rhoChemistryCombustion*> (pRxn)
+);
+autoPtr<rhoChemistryModel> chemistry = combustion->chem();
+rhoReactionThermo& thermo = chemistry->thermo();
+thermo.validate(args.executable(), "ha");
+SLGThermo slgThermo(mesh, thermo);
+basicSpecieMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+word inertSpecie(thermo.lookup("inertSpecie"));
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+
+volScalarField& p = thermo.p();
+
+#include "compressibleCreatePhi.H"
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        simple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        simple.dict(),
+        dimDensity,
+        0
+    )
+);
+
+Info << "Creating turbulence model.\n" << nl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+// Set the turbulence into the reaction model
+combustion->setTurbulence(turbulence());
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+/*
+volScalarField Hha //absolute Enthalpy
+(
+    IOobject
+    (
+        "Hha",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Hha", dimEnergy, 0.0)
+);
+volScalarField Hla //absolute Enthalpy
+(
+    IOobject
+    (
+        "Hla",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Hla", dimEnergy, 0.0)
+);
+*/
+
+#include "createRadiationModel.H"
+
+#include "createMRF.H"
+#include "createClouds.H"
+

applications/solvers/combustion/postSootFoam/postSootFoam.C

@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Application
+    simpleReactingParcelFoam
+
+Description
+    Steady state solver for compressible, turbulent flow with reacting,
+    multiphase particle clouds and optional sources/constraints.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulentFluidThermoModel.H"
+#include "basicReactingMultiphaseCloud.H"
+#include "rhoChemistryCombustion.H"
+#include "radiationModel.H"
+#include "fvOptions.H"
+#include "SLGThermo.H"
+#include "simpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+
+    turbulence->validate();
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        radiation->correct();
+
+        if(runTime.write() == true)
+        {
+            rho.write();
+            psi.write();
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/combustionModels/EDM/EDM.C

@@ -0,0 +1,195 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "EDM.H"
+#include "fvmSup.H"
+#include "localEulerDdtScheme.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+template<class Type>
+Foam::combustionModels::EDM<Type>::EDM
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    Type(modelType, mesh, phaseName),
+    finiteRate_
+    (
+        this->coeffs().lookupOrDefault("finiteRate", false)
+    ),
+    A_
+    (
+        this->coeffs().lookupOrDefault("A", 4.0)
+    ),
+    B_
+    (
+        this->coeffs().lookupOrDefault("B",0.5)
+    )
+{
+    if (finiteRate_)
+    {
+        Info<< "    using Finite-rate/Eddy Dissipation Model" << endl;
+        Info<< "    A = "<<A_<< endl;
+        Info<< "    B = "<<B_<< endl;
+    }
+    else
+    {
+        Info<< "    using Eddy Dissipation Model" << endl;
+        Info<< "    A = "<<A_<< endl;
+        Info<< "    B = "<<B_<< endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::combustionModels::EDM<Type>::~EDM()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::EDM<Type>::tc() const
+{
+    return this->chemistryPtr_->tc();
+}
+
+
+template<class Type>
+void Foam::combustionModels::EDM<Type>::correct()
+{
+    if (this->active())
+    {
+        this->chemistryPtr_->calculateEDM(finiteRate_, A_, B_);
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::combustionModels::EDM<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::EDM<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::EDM<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::EDM<Type>::read()
+{
+    if (Type::read())
+    {
+        this->coeffs().lookup("finiteRate")
+            >> finiteRate_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/combustionModels/EDM/EDM.H

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::EDM
+
+Description
+    Eddy dissipation model with finite-rate chemistry.
+                                                                    28.Feb.2017
+                                                        Combustion Lab. POSTECH
+                                                                      Karam Han
+                                              (Version upgraded by Jinwoo Park)
+
+SourceFiles
+    EDM.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef EDM_H
+#define EDM_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class EDM Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class EDM
+:
+    public Type
+{
+    // Private data
+
+        //- Select EDM or finite-rate/EDM
+        bool finiteRate_;
+        scalar A_;
+        scalar B_;
+
+protected:
+
+    // Protected Member Functions
+
+        //- Return the chemical time scale
+        tmp<volScalarField> tc() const;
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        EDM(const EDM&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const EDM&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("EDM");
+
+
+    // Constructors
+
+        //- Construct from components
+        EDM
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~EDM();
+
+
+    // 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 "EDM.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/EDM/EDMs.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 "EDM.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeCombustionTypes(EDM, psiChemistryCombustion, psiCombustionModel);
+makeCombustionTypes(EDM, rhoChemistryCombustion, rhoCombustionModel);
+
+
+// ************************************************************************* //

src/combustionModels/Make/files

@@ -13,6 +13,8 @@ rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
 diffusion/diffusions.C
 infinitelyFastChemistry/infinitelyFastChemistrys.C
 
+EDM/EDMs.C
+
 PaSR/PaSRs.C
 
 laminar/laminars.C

src/combustionModels/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.C

@@ -47,6 +47,12 @@ Foam::combustionModels::psiChemistryCombustion::~psiChemistryCombustion()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+Foam::autoPtr<Foam::psiChemistryModel>
+Foam::combustionModels::psiChemistryCombustion::chem()
+{
+    return chemistryPtr_;
+}
+
 Foam::psiReactionThermo&
 Foam::combustionModels::psiChemistryCombustion::thermo()
 {

src/combustionModels/psiCombustionModel/psiChemistryCombustion/psiChemistryCombustion.H

@@ -90,6 +90,8 @@ public:
 
     // Member Functions
 
+        autoPtr<psiChemistryModel> chem();
+
         //- Return access to the thermo package
         virtual psiReactionThermo& thermo();
 

src/combustionModels/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C

@@ -46,6 +46,11 @@ Foam::combustionModels::rhoChemistryCombustion::~rhoChemistryCombustion()
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+Foam::autoPtr<Foam::rhoChemistryModel>
+Foam::combustionModels::rhoChemistryCombustion::chem()
+{
+    return chemistryPtr_;
+}
 
 Foam::rhoReactionThermo&
 Foam::combustionModels::rhoChemistryCombustion::thermo()

src/combustionModels/rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.H

@@ -90,6 +90,8 @@ public:
 
     // Member Functions
 
+        autoPtr<rhoChemistryModel> chem();
+
         //- Return access to the thermo package
         virtual rhoReactionThermo& thermo();
 

src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H

@@ -175,6 +175,49 @@ public:
 
                 //- Return the heat release, i.e. enthalpy/sec [m2/s3]
                 virtual tmp<volScalarField> dQ() const = 0;
+
+                //- Return the mass fraction of the *specieName
+                //- 28.Nov.2017 Jinwoo Park
+                virtual tmp<volScalarField> Yspecie(const word *specieName) const = 0;
+
+                //- Return the mass fraction of the *specieName
+                //- 28.Nov.2017 Jinwoo Park
+                virtual tmp<volScalarField> MoleFracSpecie(const word *specieName) const = 0;
+
+                //- Return the epsilon
+                virtual tmp<volScalarField> epsilon() const = 0;
+
+                //- Return the k
+                virtual tmp<volScalarField> k() const = 0;
+
+                //- Return the U field
+                virtual tmp<volVectorField> U() const = 0;
+
+                //- Return the turbulent kinematic diffusivity of chemical species (Le=1 is assumed)
+                virtual tmp<volScalarField> nut() const = 0;
+
+                virtual scalar solveCMCchem
+                (
+                    scalar deltaT, scalar& rho, scalar& T,
+                    scalar& p, scalarField& Qi, scalar& Qh,
+                    scalarField& RRCMC, scalar& ChemDeltaT
+                ) = 0; //CMC
+
+                virtual scalar calculateTCMC(scalar& Qh, scalarField& Qi, scalar& Told, scalar& rho, scalar& p) = 0; //CMC
+                virtual scalar calculateRHOCMC(scalarField& Qi, scalar& T, scalar& p) = 0; //CMC
+                virtual scalar calculateHCMC(scalarField& Qi, scalar& T, scalar& rho, scalar& p) = 0; //CMC
+                virtual scalar calculateShCMC(scalarField& RRCMC, label i) = 0; //CMC
+            virtual void correction
+            (
+                scalarField& b, scalarField& Ta, scalarField& Wa, scalarField& Wb,
+                scalarField& Gab, scalarField& Gat, scalarField& Gbt,
+                scalarField& Gaa, scalarField& Gbb, scalar& Gtt, scalar& rhorho
+            ) = 0;
+
+                //- Calculate the reaction rate (EDM or finite-rate/EDM)
+                //- 11.Apr.2017 Karam Han (updated by Jinwoo Park)
+                virtual void calculateEDM(bool finiteRate, const scalar A, const scalar B) = 0;
+
 };
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C

@@ -27,6 +27,7 @@ License
 #include "reactingMixture.H"
 #include "UniformField.H"
 #include "extrapolatedCalculatedFvPatchFields.H"
+#include "wallFvPatch.H"                                             //jinwoo, implemented according to the of231 version of simpleCoalCombustionfoam from karam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -49,11 +50,12 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
         dynamic_cast<const reactingMixture<ThermoType>&>
             (this->thermo()).speciesData()
     ),
-
     nSpecie_(Y_.size()),
     nReaction_(reactions_.size()),
     Treact_(CompType::template lookupOrDefault<scalar>("Treact", 0.0)),
-    RR_(nSpecie_)
+    RR_(nSpecie_),
+    edmRR_(reactions_.size()),
+    kineticRR_(reactions_.size())
 {
     // create the fields for the chemistry sources
     forAll(RR_, fieldi)
@@ -77,6 +79,45 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
         );
     }
 
+    forAll(edmRR_, reactioni)
+    {
+        edmRR_.set
+        (
+            reactioni,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "edmRR." + reactions_[reactioni].name(),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            )
+        );
+
+        kineticRR_.set
+        (
+            reactioni,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "kineticRR." + reactions_[reactioni].name(),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            )
+        );
+    }
+
     Info<< "chemistryModel: Number of species = " << nSpecie_
         << " and reactions = " << nReaction_ << endl;
 }
@@ -585,6 +626,249 @@ Foam::chemistryModel<CompType, ThermoType>::Sh() const
     return tSh;
 }
 
+// * * * * * * * quoted by JW 20170320   (13) * * * * * * * * * * * * //
+
+
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::chemistryModel<CompType, ThermoType>::Yspecie(const word *specieName) const
+{
+    tmp<volScalarField> tYspecie
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Yspecie",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimless, 0.0)
+        )
+    );
+
+    scalarField& Yspecie = tYspecie.ref();
+
+    scalar initializationchecker=0.0;
+
+    forAll(Yspecie, celli)
+    {
+        if (Yspecie[celli] !=0.0){
+            initializationchecker = 1.0;
+        }
+    }
+
+    if (initializationchecker != 0.0){
+        Info<<"mass fraction field failed to be initialized!"<<endl;
+    }
+
+    char SpecieIndex = dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo()).species()[*specieName];
+
+
+    forAll(Yspecie, celli)
+    {
+        Yspecie[celli] = Y_[SpecieIndex][celli];
+    }
+
+    return tYspecie;
+}
+
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::chemistryModel<CompType, ThermoType>::MoleFracSpecie(const word *specieName) const
+{
+    tmp<volScalarField> tMoleFracSpecie
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "MoleFracSpecie",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimless, 0.0)
+        )
+    );
+
+    scalarField& MoleFracSpecie = tMoleFracSpecie.ref();
+
+    scalar initializationchecker=0.0;
+
+    forAll(MoleFracSpecie, celli)
+    {
+        if (MoleFracSpecie[celli] !=0.0){
+            initializationchecker = 1.0;
+        }
+    }
+
+    if (initializationchecker != 0.0){
+        Info<<"mole fraction field failed to be initialized!"<<endl;
+    }
+
+    char SpecieIndex = dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo()).species()[*specieName];
+
+    forAll(MoleFracSpecie, celli)
+    {
+
+        scalar totMol(0.0);
+
+        for (label i=0; i<nSpecie_; i++)
+        {
+            totMol += Y_[i][celli]/specieThermo_[i].W();
+        }
+
+        MoleFracSpecie[celli] = (Y_[SpecieIndex][celli]/specieThermo_[SpecieIndex].W())/totMol;
+    }
+
+    return tMoleFracSpecie;
+}
+
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::chemistryModel<CompType, ThermoType>::epsilon() const
+{
+    tmp<volScalarField> tepsilon
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "epsilon_temp",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimless, 0.0)
+        )
+    );
+
+    scalarField& epsilon = tepsilon.ref();
+
+    const volScalarField& epsilon_ref = this->db().objectRegistry::lookupObject<volScalarField>("epsilon");
+
+    forAll(epsilon, celli)
+    {
+        epsilon[celli] = epsilon_ref[celli];
+    }
+
+    return tepsilon;
+}
+
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::chemistryModel<CompType, ThermoType>::k() const
+{
+    tmp<volScalarField> tk
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "k_temp",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimless, 0.0)
+        )
+    );
+
+    scalarField& k = tk.ref();
+
+    const volScalarField& k_ref = this->db().objectRegistry::lookupObject<volScalarField>("k");
+
+    forAll(k, celli)
+    {
+        k[celli] = k_ref[celli];
+    }
+
+    return tk;
+}
+
+// * * * * * quoted by JW 20170331(15) * * * * * * //
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::volVectorField>
+Foam::chemistryModel<CompType, ThermoType>::U() const
+{
+    tmp<volVectorField> tU
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                "U_temp",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedVector("U_temp", dimLength/dimTime, vector(0.0, 0.0, 0.0))
+        )
+    );
+
+    vectorField& U = tU.ref();
+
+    const volVectorField& U_ref = this->db().objectRegistry::lookupObject<volVectorField>("U");
+
+    forAll(U, celli)
+    {
+        U[celli] = U_ref[celli];
+    }
+
+    return tU;
+}
+
+
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::chemistryModel<CompType, ThermoType>::nut() const
+{
+    tmp<volScalarField> tnut
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "nut_temp",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", pow(dimLength,2)/dimTime, 0.0)
+        )
+    );
+
+    scalarField& nut = tnut.ref();
+
+    const volScalarField& nut_ref = this->db().objectRegistry::lookupObject<volScalarField>("nut");
+
+    forAll(nut, celli)
+    {
+        nut[celli] = nut_ref[celli];
+    }
+
+    return tnut;
+}
 
 template<class CompType, class ThermoType>
 Foam::tmp<Foam::volScalarField>
@@ -877,5 +1161,313 @@ void Foam::chemistryModel<CompType, ThermoType>::solve
     NotImplemented;
 }
 
+template<class CompType, class ThermoType>
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solveCMCchem
+(
+    scalar deltaT, scalar& rho, scalar& T, scalar& p, scalarField& Qi, scalar& Qh, scalarField& RRCMC, scalar& ChemDeltaT
+)
+{
+    scalar deltaTMin = GREAT;
+    //if(!this->chemistry_)
+    //{
+    //    return deltaTMin;
+    //}
+
+    scalar rhoi = rho;
+    scalar Ti = T;
+    scalar pi = p;
+
+    scalarField c(nSpecie_);
+    scalarField c0(nSpecie_);
+
+    for(label i=0 ; i<nSpecie_ ; i++)
+    {
+        c[i] = rhoi*Qi[i]/specieThermo_[i].W();
+        c0[i] = c[i];
+    }
+
+    scalar timeLeft = deltaT;
+
+    while(timeLeft > SMALL)
+    {
+        scalar dt = timeLeft;
+        this->solve(c, Ti, pi, dt, ChemDeltaT);
+        timeLeft -= dt;
+    }
+    deltaTMin = min(ChemDeltaT, deltaTMin);
+
+    for(label i=0 ; i<nSpecie_; i++)
+    {
+        RRCMC[i] = (c[i]-c0[i])*specieThermo_[i].W()/deltaT;
+    }
+    return deltaTMin;
+}
+
+template<class CompType, class ThermoType>
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateTCMC
+(
+    scalar& Qh, scalarField& Qi, scalar& Told, scalar& rho, scalar& p
+)
+{
+    scalar rhoi = rho;
+    scalar Ti = Told;
+    scalar pi = p;
+
+    scalarField c(nSpecie_);
+
+    for (label i=0; i<nSpecie_; i++)
+    {
+        c[i] = rhoi*Qi[i]/specieThermo_[i].W();
+    }
+
+    scalar cTot = 0.0;
+
+    // update the temperature
+    cTot = sum(c);
+    ThermoType mixture(0.0*specieThermo_[0]);
+    for (label i=0; i<nSpecie_; i++)
+    {
+        mixture += (c[i]/cTot)*specieThermo_[i];
+    }
+    scalar TCMC = mixture.THa(Qh, pi,  Ti);
+    //scalar TCMCnew = mixture.THs(Qh, Ti);
+
+    return TCMC;//
+}
+
+template<class CompType, class ThermoType>
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateRHOCMC
+(
+    scalarField& Qi, scalar& T, scalar& p
+)
+{
+    scalar Ti = T; //conditional temperature [K]
+    scalar pi = p; //pressure [Pa]
+
+    scalar MeanMW(0); //mean molecular weight [kg/kmol]
+    scalar InvMeanMW(0); //inverse of MeanMW [kmol/kg]
+    scalar RHOCMC(0); //conditional density [kg/kmol]
+
+    for (label i=0; i<nSpecie_; i++) //get mean molecular weight respect to Qi
+    {
+        InvMeanMW += (Qi[i]/specieThermo_[i].W());
+    }
+    MeanMW = 1.0/InvMeanMW ;
+
+    RHOCMC = pi/(8314.4621/MeanMW*Ti); //universial gas constant 8314.51[J/kmol K] = [kg.m2/kmol.K.sec2]
+
+    return RHOCMC;//
+}
+
+template<class CompType, class ThermoType>
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateHCMC
+(
+    scalarField& Qi, scalar& T, scalar& rho, scalar& p
+)
+{
+    scalar rhoi = rho; //conditional density [kg/kmol]
+    scalar pi = p;
+    scalar Ti = T; //conditional temperature [K]
+    scalar HCMC(0); //conditional total enthalpy [J/kg]
+    scalarField c(nSpecie_);
+
+    for (label i=0; i<nSpecie_; i++)
+    {
+        c[i] = rhoi*Qi[i]/specieThermo_[i].W();
+    }
+
+    scalar cTot = 0.0;
+
+    // update the temperature
+    cTot = sum(c);
+    ThermoType mixture(0.0*specieThermo_[0]);
+    for (label i=0; i<nSpecie_; i++)
+    {
+        mixture += (c[i]/cTot)*specieThermo_[i];
+    }
+    HCMC = mixture.Ha(pi, Ti);
+
+    return HCMC;
+}
+
+template<class CompType, class ThermoType>
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateShCMC
+(
+    scalarField& RRCMC, label i
+)
+{
+    scalar tSh(0);
+    if (this->chemistry_)
+    {
+        scalar hi = specieThermo_[i].Hc();
+        //Info<<"hi = "<<hi<<endl;
+        //Info<<"RRCMC[i] = "<<RRCMC[i]<<endl;
+        tSh = hi*RRCMC[i];
+        //Info<<"tSh = "<<tSh<<endl;
+    }
+    return tSh;
+}
+
+template<class CompType, class ThermoType>
+void Foam::chemistryModel<CompType, ThermoType>::correction
+(
+        scalarField& b, scalarField& Ta, scalarField& Wa, scalarField& Wb,
+        scalarField& Gab, scalarField& Gat, scalarField& Gbt,
+        scalarField& Gaa, scalarField& Gbb, scalar& Gtt, scalar& rhorho
+)
+{
+        for(label i=0 ; i<3 ; i++)
+        {
+                b_[i] = b[i];
+                Ta_[i] = Ta[i];
+                Wa_[i] = Wa[i];
+                Wb_[i] = Wb[i];
+
+                Gab_[i] = Gab[i];
+                Gat_[i] = Gat[i];
+                Gbt_[i] = Gbt[i];
+                Gaa_[i] = Gaa[i];
+                Gbb_[i] = Gbb[i];
+        }
+        Gtt_ = Gtt;
+        rhorho_ = rhorho;
+}
+
+/*
+template<class CompType, class ThermoType>
+void Foam::chemistryModel<CompType, ThermoType>::specieMoleFrac
+(
+    char specieName
+)
+{
+}
+*/
+
+template<class CompType, class ThermoType>
+void Foam::chemistryModel<CompType, ThermoType>::calculateEDM
+(
+    bool finiteRate, const scalar A, const scalar B
+)
+{
+    if (!this->chemistry_)
+    {
+        Info<<"Reaction rates are not calculated!!"<<endl;
+        Info<<"Please check constant/chemistryProperties"<<endl;
+        return;
+    }
+
+    const volScalarField& rho = this->thermo().rho();
+    const volScalarField& k = this->db().objectRegistry::lookupObject<volScalarField>("k");
+    const volScalarField& epsilon = this->db().objectRegistry::lookupObject<volScalarField>("epsilon");
+
+    const scalarField& T = this->thermo().T();
+    const scalarField& p = this->thermo().p();
+
+    scalar EDMrate(0.0);
+
+    forAll(rho, celli)
+    {
+        scalarField RRedm(nSpecie_, 0.0);
+        scalarField RRkinetic(nSpecie_, 0.0);
+
+        forAll(reactions_, reactioni)
+        {
+            const Reaction<ThermoType>& R = reactions_[reactioni];
+
+            scalar Ymin(1.0);
+            forAll(R.lhs(), lhs_speciei)
+            {
+                const label speciei = R.lhs()[lhs_speciei].index;
+                const scalar stoiCoeff = R.lhs()[lhs_speciei].stoichCoeff;
+                if(stoiCoeff != 0)
+                {
+                    Ymin = min(Ymin, Y_[speciei][celli]/stoiCoeff/specieThermo_[speciei].W());
+                }
+            }
+
+            scalar YPmin(0.0);
+            scalar Num(0.0);
+            scalar Den(0.0);
+            forAll(R.rhs(), rhs_speciei)
+            {
+                const label speciei = R.rhs()[rhs_speciei].index;
+                const scalar stoiCoeff = R.rhs()[rhs_speciei].stoichCoeff;
+                if(stoiCoeff != 0)
+                {
+                    Num += Y_[speciei][celli];
+                    Den += stoiCoeff*specieThermo_[speciei].W();
+                }
+            }
+            YPmin = Num / Den;
+
+            Ymin = min(Ymin, B*YPmin);
+
+            EDMrate = A*rho[celli]*epsilon[celli]/k[celli]*Ymin;
+
+            scalar omegai = 0.0;
+            if(finiteRate == true)
+            {
+                scalar pf, cf, pr, cr;
+                label lRef, rRef;
+
+                scalarField c(nSpecie_, 0.0);
+                for(label i=0 ; i<nSpecie_ ; i++)
+                {
+                    const scalar Yi = Y_[i][celli];
+                    c[i] = rho[celli]*Yi/specieThermo_[i].W();
+                }
+                omegai = omega
+                (
+                    R, c, T[celli], p[celli], pf, cf, lRef, pr, cr, rRef
+                );
+            }
+
+            forAll(R.lhs(), lhs_speciei)
+            {
+                const label speciei = R.lhs()[lhs_speciei].index;
+                const scalar stoiCoeff = R.lhs()[lhs_speciei].stoichCoeff;
+
+                RRedm[speciei] -= stoiCoeff*specieThermo_[speciei].W()*EDMrate;
+                RRkinetic[speciei] -= stoiCoeff*specieThermo_[speciei].W()*omegai;
+
+                if(lhs_speciei == 0)
+                {
+                    //for post-processing
+                    edmRR_[reactioni][celli] = stoiCoeff*specieThermo_[speciei].W()*EDMrate;
+                    kineticRR_[reactioni][celli] = stoiCoeff*specieThermo_[speciei].W()*omegai;
+                }
+            }
+            forAll(R.rhs(), rhs_speciei)
+            {
+                const label speciei = R.rhs()[rhs_speciei].index;
+                const scalar stoiCoeff = R.rhs()[rhs_speciei].stoichCoeff;
+
+                RRedm[speciei] += stoiCoeff*specieThermo_[speciei].W()*EDMrate;
+                RRkinetic[speciei] += stoiCoeff*specieThermo_[speciei].W()*omegai;
+
+            }
+        }
+
+        for (label i=0; i<nSpecie_; i++)
+        {
+            if(finiteRate == true)
+            {
+                if(RRedm[i] == 0.0 || RRkinetic[i] == 0.0)
+                {
+                    RR_[i][celli] = 0.0;
+                }
+                else
+                {
+                    RR_[i][celli] = (RRedm[i]*RRkinetic[i])/(RRedm[i]+RRkinetic[i]);
+                }
+            }
+            else
+            {
+                RR_[i][celli] = RRedm[i];
+            }
+        }
+    }
+}
+
 
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H

@@ -103,6 +103,13 @@ protected:
         //- List of reaction rate per specie [kg/m3/s]
         PtrList<DimensionedField<scalar, volMesh>> RR_;
 
+        scalarField b_, Ta_, Wa_, Wb_, Gab_, Gat_, Gbt_, Gaa_, Gbb_;
+
+        scalar Gtt_, rhorho_;
+
+        PtrList<DimensionedField<scalar, volMesh>> edmRR_;
+
+        PtrList<DimensionedField<scalar, volMesh>> kineticRR_;
 
     // Protected Member Functions
 
@@ -222,6 +229,24 @@ public:
             //  and return the characteristic time
             virtual scalar solve(const scalarField& deltaT);
 
+            virtual scalar solveCMCchem
+            (
+                scalar deltaT, scalar& rho, scalar& T,
+                scalar& p, scalarField& Qi, scalar& Qh,
+                scalarField& RRCMC, scalar& ChemDeltaT
+            );
+
+            virtual scalar calculateTCMC(scalar& Qh, scalarField& Qi, scalar& Told, scalar& rho, scalar& p);
+            virtual scalar calculateRHOCMC(scalarField& Qi, scalar& T, scalar& p);
+            virtual scalar calculateHCMC(scalarField& Qi, scalar& T, scalar& rho, scalar& p);
+            virtual scalar calculateShCMC(scalarField& RRCMC, label i);
+            virtual void correction
+            (
+                scalarField& b, scalarField& Ta, scalarField& Wa, scalarField& Wb,
+                scalarField& Gab, scalarField& Gat, scalarField& Gbt,
+                scalarField& Gaa, scalarField& Gbb, scalar& Gtt, scalar& rhorho
+            );
+
             //- Return the chemical time scale
             virtual tmp<volScalarField> tc() const;
 
@@ -231,6 +256,33 @@ public:
             //- Return the heat release, i.e. enthalpy/sec [kg/m2/s3]
             virtual tmp<volScalarField> dQ() const;
 
+            //- Return mass fraction of *specieName
+            //- 28.Nov.2017 Jinwoo Park
+            virtual tmp<volScalarField> Yspecie(const word *specieName) const;
+
+            //- Return mole fraction of *specieName
+            //- 28.Nov.2017 Jinwoo Park
+            virtual tmp<volScalarField> MoleFracSpecie(const word *specieName) const;
+
+            //- Return the epsilon
+            //- 28.Nov.2017 Jinwoo Park
+            virtual tmp<volScalarField> epsilon() const;
+
+            //- Return the turbulent kinetic energy
+            //- 28.Nov.2017 Jinwoo Park
+            virtual tmp<volScalarField> k() const;
+
+            //- Return the velocity field
+            //- 28.Nov.2017 Jinwoo Park
+            virtual tmp<volVectorField> U() const;
+
+            //- Return the turbulent kinematic diffusivity of chemical species (Le=1)
+            //- 28.Nov.2017 Jinwoo Park
+            virtual tmp<volScalarField> nut() const;
+
+            //- Calculate the reaction rates (EDM or finite-rate EDM)
+            //- 11.Apr.2017 Karam Han (updated by Jinwoo Park)
+            virtual void calculateEDM(bool finiteRate, const scalar A, const scalar B);
 
         // ODE functions (overriding abstract functions in ODE.H)
 

src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModels.C

@@ -111,6 +111,23 @@ namespace Foam
         rhoChemistryModel,
         icoPoly8EThermoPhysics
     );
+
+    // Chemistry moldels based on absoluteEnthalpy
+    makeChemistryModel
+    (
+        chemistryModel,
+        rhoChemistryModel,
+        gasHaThermoPhysics
+    );
+
+    // Chemistry moldels based on absoluteInternalEnergy
+    makeChemistryModel
+    (
+        chemistryModel,
+        rhoChemistryModel,
+        gasEaThermoPhysics
+    );
+
 }
 
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolvers.C

@@ -88,6 +88,14 @@ namespace Foam
         incompressibleGasEThermoPhysics
     );
     makeChemistrySolverTypes(rhoChemistryModel, icoPoly8EThermoPhysics);
+
+    // Chemistry solvers based on absoluteEnthalpy
+    makeChemistrySolverTypes(rhoChemistryModel, gasHaThermoPhysics);
+    // makeChemistrySolverTypes(rhoChemistryModel, gasHaThermoPhysics);
+
+    // Chemistry solvers based on absoluteInternalEnergy
+    makeChemistrySolverTypes(rhoChemistryModel, gasEaThermoPhysics);
+    // makeChemistrySolverTypes(rhoChemistryModel, gasEaThermoPhysics);
 }
 
 

src/thermophysicalModels/radiation/Make/files

@@ -31,6 +31,8 @@ submodels/sootModel/sootModel/sootModel.C
 submodels/sootModel/sootModel/sootModelNew.C
 submodels/sootModel/mixtureFractionSoot/mixtureFractionSoots.C
 submodels/sootModel/noSoot/noSoot.C
+submodels/sootModel/khanGreeveSoot/khanGreeveSoots.C
+submodels/sootModel/MossBrookesSoot/MossBrookesSoots.C
 
 /* Boundary conditions */
 derivedFvPatchFields/MarshakRadiation/MarshakRadiationFvPatchScalarField.C

src/thermophysicalModels/radiation/Make/options

@@ -9,6 +9,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
@@ -22,5 +23,6 @@ LIB_LIBS = \
     -lliquidMixtureProperties \
     -lsolidProperties \
     -lliquidProperties \
+    -lchemistryModel \
     -lfiniteVolume \
     -lmeshTools

src/thermophysicalModels/radiation/radiationModels/P1/P1.C

@@ -28,6 +28,7 @@ License
 #include "fvmSup.H"
 #include "absorptionEmissionModel.H"
 #include "scatterModel.H"
+#include "sootModel.H"
 #include "constants.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -88,6 +89,19 @@ Foam::radiation::P1::P1(const volScalarField& T)
         mesh_,
         dimensionedScalar("a", dimless/dimLength, 0.0)
     ),
+    aSoot_
+    (
+        IOobject
+        (
+            "aSoot",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("aSoot", dimless/dimLength, 0.0)
+    ),
     e_
     (
         IOobject
@@ -158,6 +172,19 @@ Foam::radiation::P1::P1(const dictionary& dict, const volScalarField& T)
         mesh_,
         dimensionedScalar("a", dimless/dimLength, 0.0)
     ),
+    aSoot_
+    (
+        IOobject
+        (
+            "aSoot",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("aSoot", dimless/dimLength, 0.0)
+    ),
     e_
     (
         IOobject
@@ -212,9 +239,28 @@ bool Foam::radiation::P1::read()
 
 void Foam::radiation::P1::calculate()
 {
-    a_ = absorptionEmission_->a();
+    a_ = absorptionEmission_->a()+ soot_->aSoot();
-    e_ = absorptionEmission_->e();
+//    a_ = absorptionEmission_->a();
+    aSoot_ = soot_->aSoot();
+    e_ = absorptionEmission_->e()+ soot_->eSoot();
+//    e_ = absorptionEmission_->e();
     E_ = absorptionEmission_->E();
+
+    scalar sootradiationChecker = 0.0;
+
+    forAll(aSoot_, celli)
+    {
+        if (aSoot_[celli] != 0.0)
+        {
+            sootradiationChecker = 1.0;
+        }
+    }
+
+    if (sootradiationChecker == 0.0)
+    {
+        Info<<"soot radiation effects are unsuccessfully considered!"<<endl;
+    }
+
     const volScalarField sigmaEff(scatter_->sigmaEff());
 
     const dimensionedScalar a0 ("a0", a_.dimensions(), ROOTVSMALL);
@@ -272,7 +318,7 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::P1::Rp() const
                 IOobject::NO_WRITE,
                 false
             ),
-            4.0*absorptionEmission_->eCont()*physicoChemical::sigma
+            4.0*(absorptionEmission_->eCont()+soot_->eSoot())*physicoChemical::sigma
         )
     );
 }
@@ -286,7 +332,7 @@ Foam::radiation::P1::Ru() const
     const DimensionedField<scalar, volMesh> E =
         absorptionEmission_->ECont()()();
     const DimensionedField<scalar, volMesh> a =
-        absorptionEmission_->aCont()()();
+        absorptionEmission_->aCont()()()+soot_->aSoot()()();
 
     return a*G - E;
 }

src/thermophysicalModels/radiation/radiationModels/P1/P1.H

@@ -70,6 +70,8 @@ class P1
         //- Absorption coefficient
         volScalarField a_;
 
+        volScalarField aSoot_;
+
         //- Emission coefficient
         volScalarField e_;
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C

@@ -167,6 +167,46 @@ void Foam::radiation::fvDOM::initialise()
 
 
     // Construct absorption field for each wavelength
+    forAll(agasLambda_, lambdaI)
+    {
+        agasLambda_.set
+        (
+            lambdaI,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "agasLambda_" + Foam::name(lambdaI) ,
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                agas_
+            )
+        );
+    }
+
+    forAll(aSootLambda_, lambdaI)
+    {
+        aSootLambda_.set
+        (
+            lambdaI,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "aSootLambda_" + Foam::name(lambdaI) ,
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                aSoot_
+            )
+        );
+    }
+
     forAll(aLambda_, lambdaI)
     {
         aLambda_.set
@@ -187,6 +227,7 @@ void Foam::radiation::fvDOM::initialise()
         );
     }
 
+
     Info<< "fvDOM : Allocated " << IRay_.size()
         << " rays with average orientation:" << nl;
 
@@ -286,6 +327,32 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         mesh_,
         dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
     ),
+    agas_
+    (
+        IOobject
+        (
+            "agas",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("agas", dimless/dimLength, 0.0)
+    ),
+    aSoot_
+    (
+        IOobject
+        (
+            "aSoot",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("aSoot", dimless/dimLength, 0.0)
+    ),
     a_
     (
         IOobject
@@ -303,6 +370,8 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     nPhi_(readLabel(coeffs_.lookup("nPhi"))),
     nRay_(0),
     nLambda_(absorptionEmission_->nBands()),
+    agasLambda_(nLambda_),
+    aSootLambda_(nLambda_),
     aLambda_(nLambda_),
     blackBody_(nLambda_, T),
     IRay_(0),
@@ -375,6 +444,32 @@ Foam::radiation::fvDOM::fvDOM
         mesh_,
         dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
     ),
+    agas_
+    (
+        IOobject
+        (
+            "agas",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("agas", dimless/dimLength, 0.0)
+    ),
+    aSoot_
+    (
+        IOobject
+        (
+            "aSoot",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("aSoot", dimless/dimLength, 0.0)
+    ),
     a_
     (
         IOobject
@@ -392,6 +487,8 @@ Foam::radiation::fvDOM::fvDOM
     nPhi_(readLabel(coeffs_.lookup("nPhi"))),
     nRay_(0),
     nLambda_(absorptionEmission_->nBands()),
+    agasLambda_(nLambda_),
+    aSootLambda_(nLambda_),
     aLambda_(nLambda_),
     blackBody_(nLambda_, T),
     IRay_(0),
@@ -432,7 +529,17 @@ bool Foam::radiation::fvDOM::read()
 
 void Foam::radiation::fvDOM::calculate()
 {
-    absorptionEmission_->correct(a_, aLambda_);
+    absorptionEmission_->correct(agas_, agasLambda_);
+
+    soot_->correct(aSoot_, aSootLambda_);
+
+//    aLambda_ = agasLambda_ + aSootLambda_;
+
+    forAll(aLambda_, celli)
+    {
+        a_[celli] = agas_[celli] + aSoot_[celli];
+        aLambda_[celli] = agasLambda_[celli] + aSootLambda_[celli];
+    }
 
     updateBlackBodyEmission();
 
@@ -483,7 +590,7 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
                 false
             ),
             // Only include continuous phase emission
-            4*absorptionEmission_->aCont()*physicoChemical::sigma
+            4*(absorptionEmission_->aCont() + soot_->aSoot())*physicoChemical::sigma
         )
     );
 }
@@ -501,7 +608,7 @@ Foam::radiation::fvDOM::Ru() const
 
     // Only include continuous phase absorption
     const DimensionedField<scalar, volMesh> a =
-        absorptionEmission_->aCont()()();
+        absorptionEmission_->aCont()()() + soot_->aSoot()()();
 
     return a*G - E;
 }

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H

@@ -101,6 +101,10 @@ class fvDOM
         volScalarField Qin_;
 
         //- Total absorption coefficient [1/m]
+        volScalarField agas_;
+
+        volScalarField aSoot_;
+
         volScalarField a_;
 
         //- Number of solid angles in theta
@@ -116,6 +120,10 @@ class fvDOM
         label nLambda_;
 
         //- Wavelength total absorption coefficient [1/m]
+        PtrList<volScalarField> agasLambda_;
+
+        PtrList<volScalarField> aSootLambda_;
+
         PtrList<volScalarField> aLambda_;
 
         //- Black body

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H

@@ -37,6 +37,7 @@ SourceFiles
 
 #include "absorptionEmissionModel.H"
 #include "blackBodyEmission.H"
+#include "sootModel.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiation/submodels/sootModel/MossBrookesSoot/MossBrookesSoot.C

@@ -0,0 +1,530 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "fvmLaplacian.H"
+#include "fvmDiv.H"
+#include "linear.H"
+#include "MossBrookesSoot.H"
+#include "reactingMixture.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::MossBrookesSoot
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const word& modelType
+)
+:
+    sootModel(dict, mesh, modelType),
+    soot_
+    (
+        IOobject
+        (
+            "soot",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    sootVF_
+    (
+        IOobject
+        (
+            "sootVF",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    NucConc_
+    (
+        IOobject
+        (
+            "NucConc",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    NumDen_
+    (
+        IOobject
+        (
+            "NumDen",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    coeffsDict_(dict.subOrEmptyDict(modelType + "Coeffs")),
+    thermo_(mesh.lookupObject<fluidThermo>(basicThermo::dictName)),
+    Snet_
+    (
+        IOobject
+        (
+            "Snet",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Snet", dimMass/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    Snuc_
+    (
+        IOobject
+        (
+            "Snuc",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Snuc", dimMass/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    Smassgrow_
+    (
+        IOobject
+        (
+            "Smassgrow",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Smassgrow", dimMass/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    Soxid_
+    (
+        IOobject
+        (
+            "Soxid",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Soxid", dimMass/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    NSnuc_
+    (
+        IOobject
+        (
+            "NSnuc",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("NSnuc", dimless/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    NSagg_
+    (
+        IOobject
+        (
+            "NSagg",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("NSagg", dimless/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    NSnet_
+    (
+        IOobject
+        (
+            "NSnet",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("NSnet", dimless/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    ParameterAgg_
+    (
+        IOobject
+        (
+            "ParameterAgg",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("ParameterAgg", dimless, 0.0)
+    ),
+    Sf_
+    (
+        IOobject
+        (
+            "Sf",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Sf", dimless, 0.0)
+    ),
+    SootMeanDiameter_
+    (
+        IOobject
+        (
+            "SootMeanDiameter",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("SootMeanDiameter", dimLength, 0.0)
+    ),
+    Prtsoot(readScalar(coeffsDict_.lookup("TurbulentPrandtlNumberSoot"))),
+    rhosoot(readScalar(coeffsDict_.lookup("SootDensity"))),
+    Calpha(readScalar(coeffsDict_.lookup("Calpha"))),
+    Mp(readScalar(coeffsDict_.lookup("Mp"))),
+    Cgamma(readScalar(coeffsDict_.lookup("Cgamma"))),
+    Coxid(readScalar(coeffsDict_.lookup("Coxid"))),
+    Cw(readScalar(coeffsDict_.lookup("Cw"))),
+    etaColl(readScalar(coeffsDict_.lookup("etaColl"))),
+    Talpha(readScalar(coeffsDict_.lookup("Talpha"))),
+    Tgamma(readScalar(coeffsDict_.lookup("Tgamma"))),
+    Navog(readScalar(coeffsDict_.lookup("AvogadroNumber"))),
+    Cbeta(readScalar(coeffsDict_.lookup("Cbeta"))),
+    Boltzmann(readScalar(coeffsDict_.lookup("BoltzmannConst"))),
+    EffectonRadiation_(readBool(coeffsDict_.lookup("Radiationeffect"))),
+    cliqueNuc(readScalar(coeffsDict_.lookup("NucExpCorrectConstant"))),
+    cliqueMassGrow(readScalar(coeffsDict_.lookup("MassGrowCorrectConstant"))),
+    chemistryPtr_(CombThermoType::New(mesh))
+{
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::~MossBrookesSoot()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volVectorField>
+Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::U() const
+{
+    return this->chemistryPtr_->U();
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::nut() const
+{
+    return this->chemistryPtr_->nut();
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::MoleFracPrec() const
+{
+    word precSpecie(coeffsDict_.lookup("Precursor"));
+    const word *precspecie;
+    precspecie = &precSpecie;
+    return this->chemistryPtr_->MoleFracSpecie(precspecie);
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::MoleFracSgs() const
+{
+    word sgsSpecie(coeffsDict_.lookup("SurfaceGrowSpecie"));
+    const word *sgsspecie;
+    sgsspecie = &sgsSpecie;
+    return this->chemistryPtr_->MoleFracSpecie(sgsspecie);
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::MoleFracOxid() const
+{
+    word oxSpecie(coeffsDict_.lookup("SootOxidizerSpecie"));
+    const word *oxspecie;
+    oxspecie = &oxSpecie;
+    return this->chemistryPtr_->MoleFracSpecie(oxspecie);
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::aSoot(const label bandI) const
+{
+
+    const volScalarField T_ = thermo_.T();
+
+    tmp<volScalarField> tas
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "aSoot" + name(bandI),
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("as", dimless/dimLength, 0.0),
+            extrapolatedCalculatedFvPatchVectorField::typeName
+        )
+    );
+
+    scalarField& as = tas.ref().primitiveFieldRef();
+
+    forAll(as, celli)
+    {
+        if (EffectonRadiation_ == true)
+        {
+            as[celli] = (1232.4)*(1800.0)*soot_[celli]*(1.0+(4.8e-4)*(T_[celli]-2000.0));
+        }
+        else
+        {
+            as[celli] = 0.0;
+        }
+    }
+    tas.ref().correctBoundaryConditions();
+    return tas;
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::eSoot(const label bandI) const
+{
+
+    return aSoot(bandI);
+}
+
+template<class CombThermoType, class ThermoType>
+void Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::correct()
+{
+
+    calcSource();
+
+    const volScalarField rho_ = thermo_.rho();
+
+    const volVectorField U_ = U();
+
+    const volScalarField nut_ = nut();
+
+    const surfaceScalarField phi_
+    (
+        IOobject
+        (
+            "phi_for_soot",
+            mesh().time().timeName(),
+            mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        linearInterpolate(rho_*U_) & mesh().Sf()
+    );
+
+    fvScalarMatrix SootEqn
+    (
+        fvm::div(phi_, soot_)
+        - fvm::laplacian((rho_*nut_/Prtsoot), soot_)
+        ==
+        Snet_
+    );
+
+    fvScalarMatrix NucConcEqn
+    (
+        fvm::div(phi_, NucConc_)
+        - fvm::laplacian((rho_*nut_/Prtsoot), NucConc_)
+        ==
+        NSnet_
+    );
+
+    SootEqn.relax();
+
+    SootEqn.solve();
+
+    NucConcEqn.relax();
+
+    NucConcEqn.solve();
+
+    scalar minusindicator1 = 0.0;
+    scalar minusindicator2 = 0.0;
+
+    forAll(rho_, celli)
+    {
+        if (soot_[celli] < 0.0) {
+            soot_[celli] = 0.0;
+            minusindicator1 = 1.0;
+        }
+        if (NucConc_[celli] < 0.0) {
+            NucConc_[celli] =0.0;
+            minusindicator2 = 1.0;
+        }
+
+        NumDen_[celli] = rho_[celli]*(NucConc_[celli]*1e+15);
+
+        sootVF_[celli] = soot_[celli]*rho_[celli]/rhosoot;
+
+    }
+    if (minusindicator1 == 1.0 and minusindicator2 == 0.0) {
+        Info<< "Minus soot value! corrected to zero" <<endl;
+    }
+    else if (minusindicator1 == 0.0 and minusindicator2 ==1.0) {
+        Info<< "Minus NucConc value! corrected to zero" <<endl;
+    }
+    else if (minusindicator1 == 1.0 and minusindicator2 == 1.0) {
+        Info<< "Minus soot and NucConc value! corrected to zero" << endl;
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+void Foam::radiation::MossBrookesSoot<CombThermoType, ThermoType>::calcSource()
+{
+    const volScalarField rho_ = thermo_.rho();
+    const volScalarField T_ = thermo_.T();
+    const volScalarField p_ = thermo_.p();
+
+    const volScalarField MoleFracPrec_ = MoleFracPrec();
+
+    const volScalarField MoleFracSgs_ = MoleFracSgs();
+
+    const volScalarField MoleFracOxid_ = MoleFracOxid();
+
+    scalar ErrorIndex_p = 0;
+    scalar ErrorIndex_rho = 0;
+    scalar ErrorIndex_T = 0;
+    scalar ErrorIndex_precursor = 0;
+    scalar ErrorIndex_oxidizer = 0;
+
+    forAll(rho_, celli)
+    {
+        if (p_[celli] < 0) {
+            ErrorIndex_p = 1;
+        }
+        if (rho_[celli] < 0) {
+            ErrorIndex_rho = 1;
+        }
+        if (T_[celli] < 0) {
+            ErrorIndex_T = 1;
+        }
+        if (MoleFracPrec_[celli] < 0) {
+            ErrorIndex_precursor = 1;
+        }
+        if (MoleFracOxid_[celli] < 0) {
+            ErrorIndex_oxidizer = 1;
+        }
+    }
+
+    if (ErrorIndex_p == 1){
+        Info<<"Warning! Improper pressure field value has been fed into the soot calculation!"<<endl;
+    }
+    if (ErrorIndex_rho == 1){
+        Info<<"Warning! Improper rho field value has been fed into the soot calculation!"<<endl;
+    }
+    if (ErrorIndex_T == 1){
+        Info<<"Warning! Improper T field value has been fed into the soot calculation!"<<endl;
+    }
+    if (ErrorIndex_precursor == 1){
+        Info<<"Warning! Improper soot precursor mole fraction field value has been fed into the soot calculation!"<<endl;
+    }
+    if (ErrorIndex_oxidizer == 1){
+        Info<<"Warning! Improper soot oxidizer mole fraction field value has been fed into the soot calculation!"<<endl;
+    }
+
+    forAll(rho_, celli)
+    {
+        const scalar MoleFracPreci = MoleFracPrec_[celli];
+        const scalar MoleFracSgsi = MoleFracSgs_[celli];
+        const scalar MoleFracOxidi = MoleFracOxid_[celli];
+        const scalar Ti = T_[celli];
+        const scalar pi = p_[celli];
+        const scalar rhoi = rho_[celli];
+        const scalar sooti = soot_[celli];
+        const scalar NumDeni = NumDen_[celli];
+
+        const scalar MolarConcPreci = (MoleFracPreci*pi)/(8314.5*Ti);
+        const scalar MolarConcSgsi = (MoleFracSgsi*pi)/(8314.5*Ti);
+        const scalar MolarConcOxidi = (MoleFracOxidi*pi)/(8314.5*Ti);
+
+        scalar SootMeanDiameteri = 0.0;
+
+        if (NumDeni == 0.0) {
+            SootMeanDiameter_[celli] = 0.0;
+            Sf_[celli] = 0.0;
+        }
+        else {
+            const scalar parameter = (6.0/3.1416)*(rhoi/rhosoot)*(sooti/NumDeni);
+            SootMeanDiameter_[celli] = pow(parameter,0.33333);
+            SootMeanDiameteri = SootMeanDiameter_[celli];
+            Sf_[celli] = (3.1416)*pow(SootMeanDiameteri,2.0)*NumDeni;
+        }
+
+        Snuc_[celli] = Mp*Calpha*exp(-cliqueNuc*Talpha/Ti)*MolarConcPreci;
+
+        Smassgrow_[celli] = Cgamma*MolarConcSgsi*exp(-cliqueMassGrow*Tgamma/Ti)*Sf_[celli];
+
+        Soxid_[celli] = Coxid*Cw*etaColl*MolarConcOxidi*pow(Ti,0.5)*Sf_[celli];
+
+        Snet_[celli] = Snuc_[celli] + Smassgrow_[celli] - Soxid_[celli];
+
+        NSnuc_[celli] = (Calpha*Navog*MolarConcPreci*exp(-cliqueNuc*Talpha/Ti))*(1e-15);
+
+        ParameterAgg_[celli] = 24.0*Boltzmann*Ti/rhosoot;
+
+        NSagg_[celli] = (Cbeta*pow((ParameterAgg_[celli]),0.5)*pow(SootMeanDiameteri,0.5)*pow(NumDeni,2.0))*(1e-15);
+
+        NSnet_[celli] = NSnuc_[celli] - NSagg_[celli];
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiation/submodels/sootModel/MossBrookesSoot/MossBrookesSoot.H

@@ -0,0 +1,294 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::radiation::MossBrookesSoot
+
+Description
+    This soot model solves two transport equation for soot mass fraction and soot normalized nuclei concentration.
+    The detail of the model strictly follows Moss & Brookes model in fluent.
+
+    Parameter settings are available in constant/radiationProperties.
+
+    Setting example:
+
+    sootModel MossBrookesSoot<rhoChemistryModel, gasHaThermoPhysics>;
+
+    MossBrookesSootCoeffs
+    {
+        SootDensity            1800.0;
+        Calpha                 54.0;
+        Mp                     144;
+        Cgamma                 11700;
+        Coxid                  1.0;
+        Cw                     105.8125;
+        etaColl                0.04;
+        Talpha                 21100.0;
+        Tgamma                 12100.0;
+        AvogadroNumber         6.022e+26;
+        Cbeta                  1.0;
+        BoltzmannConst         1.38e-23;
+        Radiationeffect        true;
+
+        Precursor              C2H2;
+        SurfaceGrowSpecie      C2H2;
+        SootOxidizerSpecie     OH;
+    }
+
+                                                                                                2017.11.28
+                                                                                    POSTECH Combustion Lab
+                                                                                               Jinwoo Park
+
+SourceFiles
+    MossBrookesSoot.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MossBrookesSoot_H
+#define MossBrookesSoot_H
+
+#include "interpolationLookUpTable.H"
+#include "sootModel.H"
+#include "HashTable.H"
+
+#include "fluidThermo.H"
+#include "Reaction.H"
+#include "autoPtr.H"
+#include "psiChemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class MossBrookesSoot Declaration
+\*---------------------------------------------------------------------------*/
+template<class CombThermoType, class ThermoType>
+class MossBrookesSoot
+:
+    public sootModel
+
+{
+
+    // Private data
+
+        //- Soot mass fraction
+        volScalarField soot_;
+
+        //- Soot volume fraction
+        volScalarField sootVF_;
+
+        //- Normalized nuclei concentration of soot
+        volScalarField NucConc_;
+
+        //- Number density of soot
+        volScalarField NumDen_;
+
+        //- Soot model dictionary
+        dictionary coeffsDict_;
+
+        //- Thermo package
+        const fluidThermo& thermo_;
+
+        //- Net rate of soot generation
+        volScalarField Snet_;
+
+        //- The rate of soot formation by nucleation
+        volScalarField Snuc_;
+
+        //- The rate of soot formation by mass growth
+        volScalarField Smassgrow_;
+
+        //- The rate of soot formation by soot oxidation
+        volScalarField Soxid_;
+
+        //- The rate of soot number density formation by nucleation
+        volScalarField NSnuc_;
+
+        //- The rate of soot number density formation by agglomeration
+        volScalarField NSagg_;
+
+        //- Net rate of soot numbe density formation
+        volScalarField NSnet_;
+
+        //- Parameter for Agglomeration source term
+        volScalarField ParameterAgg_;
+
+        //- Available surface area of soot oxidation per volume
+        volScalarField Sf_;
+
+        //- Soot particle mean diameter
+        volScalarField SootMeanDiameter_;
+
+        //- Soot turbulent Prandtl number
+        scalar Prtsoot;
+
+        //- Density of soot particle
+        scalar rhosoot;
+
+        //- Model constant for soot inception rate
+        scalar Calpha;
+
+        //- Mass of an incipient soot particle
+        scalar Mp;
+
+        //- Surface growth rate scaling factor
+        scalar Cgamma;
+
+        //- Oxidation rate scaling parameter
+        scalar Coxid;
+
+        //- Oxidation model constant
+        scalar Cw;
+
+        //- Collisional efficiency parameter
+        scalar etaColl;
+
+        //- Activation temperature for soot inception
+        scalar Talpha;
+
+        //- Activation temperature for soot growth rate
+        scalar Tgamma;
+
+        //- Avogadro's number
+        scalar Navog;
+
+        //- Model constant for coagulation rate
+        scalar Cbeta;
+
+        //- Boltzmann constant
+        scalar Boltzmann;
+
+        //- Soot radiative effect on/off flag
+        bool EffectonRadiation_;
+
+        //- Correction constant for Exponent in Nucleation
+        scalar cliqueNuc;
+
+        //- Correction constant for Exponent in Soot Mass Grow
+        scalar cliqueMassGrow;
+
+        //- Soot precursor specie name
+        const word precSpecie;
+
+        //- Soot surface growth specie name
+        const word sgsSpecie;
+
+        //- Soot oxidation specie name
+        const word oxSpecie;
+
+protected:
+
+    // Protected data
+
+        //- Pointer to chemistry model
+        autoPtr<CombThermoType> chemistryPtr_;
+
+    // Protected Member Functions
+
+        //- Return velocity field U
+        tmp<volVectorField> U() const;
+
+        //- Return nut
+        tmp<volScalarField> nut() const;
+
+        //- Return mole fraction of soot precursor
+        tmp<volScalarField> MoleFracPrec() const;
+
+        //- Return mole fraction of soot surface growth specie
+        tmp<volScalarField> MoleFracSgs() const;
+
+        //- Return mole fraction of soot oxidation specie
+        tmp<volScalarField> MoleFracOxid() const;
+
+public:
+
+    //- Runtime type information
+    TypeName("MossBrookesSoot");
+
+
+    // Constructors
+
+        //- Construct from components
+        MossBrookesSoot
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const word& modelType
+        );
+
+
+    //- Destructor
+    virtual ~MossBrookesSoot();
+
+
+        // Edit
+
+            //- Main update/correction routine
+            virtual void correct();
+
+            //- Calculation of the source term routine
+            virtual void calcSource();
+
+
+        // Access
+
+            //- Return Ysoot
+            const volScalarField& soot() const
+            {
+                return soot_;
+            }
+
+            // Soot absorption coefficient
+
+                //- Soot absorption coefficient (net)
+                tmp<volScalarField> aSoot(const label bandI = 0) const;
+
+            // Soot emission coefficient
+
+                //- Soot emission coefficient (net)
+                tmp<volScalarField> eSoot(const label bandI = 0) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "MossBrookesSoot.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/sootModel/MossBrookesSoot/MossBrookesSoots.C

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "MossBrookesSoot.H"
+#include "makeSootTypes.H"
+#include "thermoPhysicsTypes.H"
+#include "psiChemistryModel.H"
+#include "rhoChemistryModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makeSootTypesCombThermo(MossBrookesSoot, psiChemistryModel, gasHThermoPhysics);
+makeSootTypesCombThermo(MossBrookesSoot, psiChemistryModel, gasEThermoPhysics);
+makeSootTypesCombThermo(MossBrookesSoot, psiChemistryModel, gasHaThermoPhysics);
+makeSootTypesCombThermo(MossBrookesSoot, psiChemistryModel, gasEaThermoPhysics);
+
+makeSootTypesCombThermo(MossBrookesSoot, rhoChemistryModel, gasHThermoPhysics);
+makeSootTypesCombThermo(MossBrookesSoot, rhoChemistryModel, gasEThermoPhysics);
+makeSootTypesCombThermo(MossBrookesSoot, rhoChemistryModel, gasHaThermoPhysics);
+makeSootTypesCombThermo(MossBrookesSoot, rhoChemistryModel, gasEaThermoPhysics);
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/sootModel/khanGreeveSoot/khanGreeveSoot.C

@@ -0,0 +1,394 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "fvmLaplacian.H"
+#include "fvmDiv.H"
+#include "linear.H"
+#include "khanGreeveSoot.H"
+#include "reactingMixture.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::khanGreeveSoot
+
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const word& modelType
+)
+:
+    sootModel(dict, mesh, modelType),
+    soot_
+    (
+        IOobject
+        (
+            "soot",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    sootVF_
+    (
+        IOobject
+        (
+            "sootVF",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    coeffsDict_(dict.subOrEmptyDict(modelType + "Coeffs")),
+    thermo_(mesh.lookupObject<fluidThermo>(basicThermo::dictName)),
+    Snet_
+    (
+        IOobject
+        (
+            "Snet",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Snet", dimMass/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    Sform_
+    (
+        IOobject
+        (
+            "Sform",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Sform", dimMass/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    Scomb_
+    (
+        IOobject
+        (
+            "Scomb",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Scomb", dimMass/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    S1_
+    (
+        IOobject
+        (
+            "S1",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("S1", dimMass/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    S2_
+    (
+        IOobject
+        (
+            "S2",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("S2", dimMass/(pow3(dimLength)*dimTime), 0.0)
+    ),
+    Prtsoot(readScalar(coeffsDict_.lookup("TurbulentPrandtlNumberSoot"))),
+    rhosoot(readScalar(coeffsDict_.lookup("SootDensity"))),
+    Cs_(readScalar(coeffsDict_.lookup("Cs"))),
+    r_(readScalar(coeffsDict_.lookup("r"))),
+    Ta_(readScalar(coeffsDict_.lookup("ActivationTemperature"))),
+    A_(readScalar(coeffsDict_.lookup("A"))),
+    nufuel_(readScalar(coeffsDict_.lookup("nufuel"))),
+    nusoot_(readScalar(coeffsDict_.lookup("nusoot"))),
+    Eqv_
+    (
+        IOobject
+        (
+            "Eqv",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Eqv", dimless, 0.0)
+    ),
+    EqvMax(readScalar(coeffsDict_.lookup("EquivalenceRatioMax"))),
+    EqvMin(readScalar(coeffsDict_.lookup("EquivalenceRatioMin"))),
+    EffectonRadiation_(readBool(coeffsDict_.lookup("Radiationeffect"))),
+    chemistryPtr_(CombThermoType::New(mesh))
+{
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::~khanGreeveSoot()
+
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::Yox() const
+{
+    word oxidSpecie(coeffsDict_.lookup("Oxidizer"));
+    const word *oxidspecie;
+    oxidspecie = &oxidSpecie;
+    return this->chemistryPtr_->Yspecie(oxidspecie);
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::Yfuel() const
+{
+    word fuelSpecie(coeffsDict_.lookup("Fuel"));
+    const word *fuelspecie;
+    fuelspecie = &fuelSpecie;
+    return this->chemistryPtr_->Yspecie(fuelspecie);
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::MoleFracfuel() const
+{
+    word fuelSpecie(coeffsDict_.lookup("Fuel"));
+    const word *fuelspecie;
+    fuelspecie = &fuelSpecie;
+    return this->chemistryPtr_->MoleFracSpecie(fuelspecie);
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::epsilon() const
+{
+    return this->chemistryPtr_->epsilon();
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::k() const
+{
+    return this->chemistryPtr_->k();
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volVectorField>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::U() const
+{
+    return this->chemistryPtr_->U();
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::nut() const
+{
+    return this->chemistryPtr_->nut();
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::aSoot(const label bandI) const
+{
+
+    const volScalarField T_ = thermo_.T();
+
+    tmp<volScalarField> tas
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "aSoot" + name(bandI),
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("as", dimless/dimLength, 0.0),
+            extrapolatedCalculatedFvPatchVectorField::typeName
+        )
+    );
+
+    scalarField& as = tas.ref().primitiveFieldRef();
+
+    forAll(as, celli)
+    {
+        if (EffectonRadiation_ == true)
+        {
+            as[celli] = (1232.4)*(1800.0)*soot_[celli]*(1.0+(4.8e-4)*(T_[celli]-2000.0));
+        }
+        else
+        {
+            as[celli] = 0.0;
+        }
+    }
+    tas.ref().correctBoundaryConditions();
+    return tas;
+}
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::eSoot(const label bandI) const
+{
+
+    return aSoot(bandI);
+}
+
+template<class CombThermoType, class ThermoType>
+void Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::correct()
+{
+
+    calcSource();
+
+    const volScalarField rho_ = thermo_.rho();
+
+    const volVectorField U_ = U();
+
+    const volScalarField nut_ = nut();
+
+    const surfaceScalarField phi_
+    (
+        IOobject
+        (
+            "phi_for_soot",
+            mesh().time().timeName(),
+            mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        linearInterpolate(rho_*U_) & mesh().Sf()
+    );
+
+    fvScalarMatrix SootEqn
+    (
+        fvm::div(phi_, soot_)
+        - fvm::laplacian((rho_*nut_/Prtsoot), soot_)
+        ==
+        Snet_
+    );
+
+    SootEqn.relax();
+
+    SootEqn.solve();
+
+    forAll(rho_, celli)
+    {
+        if (soot_[celli] < 1e-30){
+            soot_[celli] = 0.0;
+        }
+
+        sootVF_[celli] = soot_[celli]*rho_[celli]/rhosoot;
+
+    }
+
+}
+
+
+template<class CombThermoType, class ThermoType>
+void Foam::radiation::khanGreeveSoot<CombThermoType, ThermoType>::calcSource()
+{
+    const volScalarField rho_ = thermo_.rho();
+    const volScalarField T_ = thermo_.T();
+    const volScalarField p_ = thermo_.p();
+
+    const volScalarField Yox_ = Yox();
+
+    const volScalarField Yfuel_ = Yfuel();
+
+    const volScalarField MoleFracfuel_ = MoleFracfuel();
+
+    const volScalarField epsilon_ = epsilon();
+
+    const volScalarField k_ = k();
+
+    forAll(rho_, celli)
+    {
+        const scalar Yoxi = Yox_[celli];
+        const scalar Yfueli = Yfuel_[celli];
+
+        Eqv_[celli] = (Yfueli/Yoxi)*nufuel_;
+
+        if (Eqv_[celli] > EqvMax){
+            Eqv_[celli] = 0.0;
+        }
+        else if (Eqv_[celli] < EqvMin){
+            Eqv_[celli] = 0.0;
+        }
+
+        const scalar pi = p_[celli];
+        const scalar MoleFracfueli = MoleFracfuel_[celli];
+        const scalar Ti = T_[celli];
+
+        Sform_[celli] = Cs_*pi*MoleFracfueli*pow(Eqv_[celli],r_)*exp(-Ta_/Ti);
+
+        const scalar rhoi = rho_[celli];
+        const scalar sooti = soot_[celli];
+        const scalar epsiloni = epsilon_[celli];
+        const scalar ki = k_[celli];
+
+        S1_[celli] = A_*rhoi*sooti*epsiloni/ki;
+
+        if (sooti == 0.0 and Yfueli == 0.0 ){
+            S2_[celli] = 0.0;
+        }
+        else {
+            S2_[celli] = A_*rhoi*(Yoxi/nusoot_)*(sooti*nusoot_/(sooti*nusoot_+(Yfueli)*nufuel_))*epsiloni/ki;
+
+        }
+
+        Scomb_[celli] = min(S1_[celli], S2_[celli]);
+
+        Snet_[celli] = Sform_[celli] - Scomb_[celli];
+    }
+
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiation/submodels/sootModel/khanGreeveSoot/khanGreeveSoot.H

@@ -0,0 +1,254 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::radiation::khanGreeveSoot
+
+Description
+    This soot model solves transport equation for soot mass fraction.
+    The detail of the model can be describe in the fluent theory guide.
+
+    Parameter setup is available in constant/radiationProperties file.
+
+    Setup example:
+
+    sootModel khanGreeveSoot<rhoChemistryModel,gasHaThermoPhysics>;
+
+    khanGreeveSootCoeffs
+    {
+        SootDensity                1800.0;
+        Cs                         1.5;
+        r                          3;
+        ActivationTemperature      20000;
+        A                          4;
+        nufuel                     3.2;
+        nusoot                     2.6667;
+        EquivalenceRatioMax        3.0;
+        EquivalenceRatioMin        1.67;
+        Radiationeffect            true;
+
+        Fuel                       PHC3H7;
+        Oxidizer                   O2;
+    }
+
+                                                                                                2017.11.28
+                                                                                    POSTECH Combustion Lab
+                                                                                               Jinwoo Park
+
+SourceFiles
+    khanGreeveSoot.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef khanGreeveSoot_H
+#define khanGreeveSoot_H
+
+#include "interpolationLookUpTable.H"
+#include "sootModel.H"
+#include "HashTable.H"
+
+#include "fluidThermo.H"
+#include "Reaction.H"
+#include "autoPtr.H"
+#include "psiChemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class khanGreeveSoot Declaration
+\*---------------------------------------------------------------------------*/
+template<class CombThermoType, class ThermoType>
+class khanGreeveSoot
+:
+    public sootModel
+
+{
+
+    // Private data
+
+         //- Soot mass fraction
+        volScalarField soot_;
+
+        //- Soot volume fraction
+        volScalarField sootVF_;
+
+        //- Soot model dictionary
+        dictionary coeffsDict_;
+
+        //- Thermo package
+        const fluidThermo& thermo_;
+
+        //- Net rate of soot generation
+        volScalarField Snet_;
+
+        //- The rate of soot formation
+        volScalarField Sform_;
+
+        //- The rate of soot combustion
+        volScalarField Scomb_;
+
+        //- The 1st rate of soot combustion
+        volScalarField S1_;
+
+        //- The 2nd rate of soot combustion
+        volScalarField S2_;
+
+        //- Soot turbulent Prandtl number
+        scalar Prtsoot;
+
+        //- Soot density
+        scalar rhosoot;
+
+        //- Soot formation constant
+        scalar Cs_;
+
+        //- Equivalence ratio exponent
+        scalar r_;
+
+        //- Activation temperature
+        scalar Ta_;
+
+        //- Magnussen model constant
+        scalar A_;
+
+        //- Mass stoichiometry for the fuel oxidation
+        scalar nufuel_;
+
+        //- Mass stoichiometry for the soot oxidation
+        scalar nusoot_;
+
+        //- Equivalence ratio for the fuel
+        volScalarField Eqv_;
+
+        //- Maximum equivalence ratio
+        scalar EqvMax;
+
+        //- Minimum equivalence ratio
+        scalar EqvMin;
+
+        //- Soot radiative effect on/off flag
+        bool EffectonRadiation_;
+
+protected:
+
+    // Protected data
+
+        //- Pointer to chemistry model
+        autoPtr<CombThermoType> chemistryPtr_;
+
+    // Protected Member Functions
+
+        //- Return the oxidizer mass fraction
+        tmp<volScalarField> Yox() const;
+
+        //- Return the fuel mass fraction
+        tmp<volScalarField> Yfuel() const;
+
+        //- Return the mole fraction of the fuel
+        tmp<volScalarField> MoleFracfuel() const;
+
+        //- Return the epsilon
+        tmp<volScalarField> epsilon() const;
+
+        //- Return the k
+        tmp<volScalarField> k() const;
+
+        //- Return the velocity field U
+        tmp<volVectorField> U() const;
+
+        //- Return the nut
+        tmp<volScalarField> nut() const;
+
+public:
+
+    //- Runtime type information
+    TypeName("khanGreeveSoot");
+
+
+    // Constructors
+
+        //- Construct from components
+        khanGreeveSoot
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const word& modelType
+        );
+
+
+    //- Destructor
+    virtual ~khanGreeveSoot();
+
+
+        // Edit
+
+            //- Main update/correction routine
+            virtual void correct();
+
+            //- Calculation of the source term routine
+            virtual void calcSource();
+
+
+        // Access
+
+            //- Return Ysoot
+            const volScalarField& soot() const
+            {
+                return soot_;
+            }
+
+            // Soot absorption coefficient
+
+                //- Soot absorption coefficient (net)
+                tmp<volScalarField> aSoot(const label bandI = 0) const;
+
+            // Soot emission coefficient
+
+                //- Soot emission coefficient (net)
+                tmp<volScalarField> eSoot(const label bandI = 0) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "khanGreeveSoot.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/sootModel/khanGreeveSoot/khanGreeveSoots.C

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "khanGreeveSoot.H"
+#include "makeSootTypes.H"
+#include "thermoPhysicsTypes.H"
+#include "psiChemistryModel.H"                                                           //JW 20170317   (12)
+#include "rhoChemistryModel.H"                                                           //JW 20170317   (12)
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makeSootTypesCombThermo(khanGreeveSoot, psiChemistryModel, gasHThermoPhysics);                    //JW 20170317   (12)
+makeSootTypesCombThermo(khanGreeveSoot, psiChemistryModel, gasEThermoPhysics);                    //JW 20170317   (12)
+makeSootTypesCombThermo(khanGreeveSoot, psiChemistryModel, gasHaThermoPhysics);
+makeSootTypesCombThermo(khanGreeveSoot, psiChemistryModel, gasEaThermoPhysics);
+
+makeSootTypesCombThermo(khanGreeveSoot, rhoChemistryModel, gasHThermoPhysics);                    //JW 20170317   (12)
+makeSootTypesCombThermo(khanGreeveSoot, rhoChemistryModel, gasEThermoPhysics);                    //JW 20170317   (12)
+makeSootTypesCombThermo(khanGreeveSoot, rhoChemistryModel, gasHaThermoPhysics);
+makeSootTypesCombThermo(khanGreeveSoot, rhoChemistryModel, gasEaThermoPhysics);
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/sootModel/sootModel/makeSootTypes.H

@@ -31,6 +31,33 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#define makeSootTypesCombThermo(sootModelType, Comp, Thermo)                   \
+    typedef Foam::radiation::sootModelType                                     \
+        <Foam::Comp, Foam::Thermo>                                             \
+        sootModelType##Comp##Thermo;                                           \
+                                                                               \
+    defineTemplateTypeNameAndDebugWithName                                     \
+    (                                                                          \
+        sootModelType##Comp##Thermo,                                           \
+        #sootModelType"<"#Comp","#Thermo">",                                   \
+        0                                                                      \
+    );                                                                         \
+                                                                               \
+    namespace Foam                                                             \
+    {                                                                          \
+        namespace radiation                                                    \
+        {                                                                      \
+            typedef sootModelType<Comp, Thermo> sootModelType##Comp##Thermo;   \
+                                                                               \
+            addToRunTimeSelectionTable                                         \
+            (                                                                  \
+                sootModel,                                                     \
+                sootModelType##Comp##Thermo,                                   \
+                dictionary                                                     \
+            );                                                                 \
+        }                                                                      \
+    }
+
 #define makeSootTypesThermo(sootModelType, Thermo)                             \
     typedef Foam::radiation::sootModelType<Foam::Thermo>                       \
         sootModelType##Thermo;                                                 \

src/thermophysicalModels/radiation/submodels/sootModel/sootModel/sootModel.C

@@ -58,5 +58,58 @@ Foam::radiation::sootModel::~sootModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::sootModel::aSoot(const label bandI) const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "aSoot",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
+        )
+    );
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::sootModel::eSoot(const label bandI) const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "eSoot",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
+        )
+    );
+}
+
+void Foam::radiation::sootModel::correct
+(
+    volScalarField& a,
+    PtrList<volScalarField>& aj
+) const
+{
+    a = this->aSoot();
+    aj[0] = a;
+}
 
 // ************************************************************************* //

src/thermophysicalModels/radiation/submodels/sootModel/sootModel/sootModel.H

@@ -138,6 +138,25 @@ public:
 
                     //- Return const reference to soot
                     virtual const volScalarField& soot() const = 0;
+
+                // Soot absorption coefficient
+
+                    //- Soot absorption coefficient (net)
+                    virtual tmp<volScalarField> aSoot(const label bandI = 0) const;
+
+                // Soot emission coefficient
+
+                    //- Soot emission coefficient (net)
+                    virtual tmp<volScalarField> eSoot(const label bandI = 0) const;
+
+                // Correct absorption coefficients
+
+                    virtual void correct
+                    (
+                        volScalarField& a,
+                        PtrList<volScalarField>& aj
+                    ) const;
+
 };
 
 

src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/makeChemistryReaders.C

@@ -86,6 +86,19 @@ makeChemistryReaderType(foamChemistryReader, hConstSolidThermoPhysics);
 makeChemistryReaderType(foamChemistryReader, hPowerSolidThermoPhysics);
 makeChemistryReaderType(foamChemistryReader, hExpKappaConstSolidThermoPhysics);
 
+
+// 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

src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermos.C

@@ -212,6 +212,47 @@ makeReactionThermo
     specie
 );
 
+//reactionthermo using absoluteEnthalpy
+
+makeReactionThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    homogeneousMixture,
+    sutherlandTransport,
+    absoluteEnthalpy,
+    janafThermo,
+    perfectGas,
+    specie
+);
+
+makeReactionThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    inhomogeneousMixture,
+    sutherlandTransport,
+    absoluteEnthalpy,
+    janafThermo,
+    perfectGas,
+    specie
+);
+
+makeReactionThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    veryInhomogeneousMixture,
+    sutherlandTransport,
+    absoluteEnthalpy,
+    janafThermo,
+    perfectGas,
+    specie
+);
+
 
 // Multi-component thermo for internal energy
 
@@ -260,6 +301,29 @@ makeReactionMixtureThermo
     icoPoly8EThermoPhysics
 );
 
+// Multi-component thermo for absolute enthalpy
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    gasHaThermoPhysics
+);
+
+
+// Multi-component thermo for internal energy
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    gasEaThermoPhysics
+);
+
 
     // Multi-component reaction thermo
 
@@ -424,6 +488,47 @@ makeReactionMixtureThermo
     gasHThermoPhysics
 );
 
+// Multi-component reaction thermo for absolute enthalpy
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    reactingMixture,
+    gasHaThermoPhysics
+);
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    singleStepReactingMixture,
+    gasHaThermoPhysics
+);
+
+
+// Multi-component reaction thermo for internal energy
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    reactingMixture,
+    gasEaThermoPhysics
+);
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    singleStepReactingMixture,
+    gasEaThermoPhysics
+);
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/thermophysicalModels/specie/include/reactionTypes.H

@@ -68,6 +68,12 @@ namespace Foam
         incompressibleGasEReaction;
 
     typedef Reaction<icoPoly8EThermoPhysics> icoPoly8EReaction;
+
+    // absolute enthalpy based reactions
+    typedef Reaction<gasHaThermoPhysics> gasHaReaction;
+
+    // absolute internal energy based reactions
+    typedef Reaction<gasEaThermoPhysics> gasEaReaction;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/specie/include/thermoPhysicsTypes.H

@@ -40,6 +40,8 @@ Description
 
 #include "sensibleEnthalpy.H"
 #include "sensibleInternalEnergy.H"
+#include "absoluteEnthalpy.H"
+#include "absoluteInternalEnergy.H"
 #include "thermo.H"
 #include "sutherlandTransport.H"
 #include "constTransport.H"
@@ -188,6 +190,40 @@ namespace Foam
         >,
         8
     > icoPoly8EThermoPhysics;
+
+
+    // 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;
+
+
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/specie/reaction/reactions/makeReactions.C

@@ -103,6 +103,10 @@ namespace Foam
     )
     makeReactions(incompressibleGasEThermoPhysics, incompressibleGasEReaction)
     makeReactions(icoPoly8EThermoPhysics, icoPoly8EReaction)
+
+    // absolute enthalpy based reactions
+    makeReactions(gasHaThermoPhysics, gasHaReaction)
+    makeReactions(gasEaThermoPhysics, gasEaReaction)
 }
 
 // ************************************************************************* //