separate soot models require turbulence model

src/Allwmake

@@ -61,11 +61,10 @@ parallel/Allwmake $targetType $*
 wmake $targetType ODE
 wmake $targetType randomProcesses
 
-wmakeLnIncludeAll TurbulenceModels
-
 transportModels/Allwmake $targetType $*
 thermophysicalModels/Allwmake $targetType $*
 TurbulenceModels/Allwmake $targetType $*
+wmake $targetType thermophysicalModels/radiationSootModel
 wmake $targetType combustionModels
 regionModels/Allwmake $targetType $*
 lagrangian/Allwmake $targetType $*

src/thermophysicalModels/radiation/Make/files

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

src/thermophysicalModels/radiationSootModel/Make/files

@@ -0,0 +1,5 @@
+/* Soot model */
+khanGreeveSoot/khanGreeveSoot.C
+MossBrookesSoot/MossBrookesSoot.C
+
+LIB = $(FOAM_LIBBIN)/libradiationSootModels

src/thermophysicalModels/radiationSootModel/Make/options

@@ -0,0 +1,30 @@
+EXE_INC = \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/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/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+LIB_LIBS = \
+    -lradiationModels \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lsolidThermo \
+    -lSLGThermo \
+    -lsolidMixtureProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lliquidProperties \
+    -lfiniteVolume \
+    -lmeshTools

src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.C

@@ -0,0 +1,487 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "MossBrookesSoot.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvm.H"
+#include "basicSpecieMixture.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(MossBrookesSoot, 0);
+
+        addToRunTimeSelectionTable
+        (
+            sootModel,
+            MossBrookesSoot,
+            dictionary
+        );
+    }
+}
+
+
+const Foam::basicSpecieMixture&
+Foam::radiation::MossBrookesSoot::checkThermo
+(
+    const fluidThermo& thermo
+)
+{
+    if (isA<basicSpecieMixture>(thermo))
+    {
+        return dynamic_cast<const basicSpecieMixture& >
+        (
+            thermo
+        );
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Inconsistent thermo package for " << thermo.type()
+            << "Please select a thermo package based on "
+            << "basicSpecieMixture" << exit(FatalError);
+
+        return dynamic_cast<const basicSpecieMixture& >
+        (
+            thermo
+        );
+    }
+
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::MossBrookesSoot::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)),
+    mixture_(checkThermo(thermo_)),
+    turbulence_(mesh.lookupObject<compressibleTurbulenceModel>(turbulenceModel::propertiesName)),
+    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)
+    ),
+    aggParam_
+    (
+        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"))),
+    radiationActive_(readBool(coeffsDict_.lookup("Radiationeffect"))),
+    cliqueNuc(readScalar(coeffsDict_.lookup("NucExpCorrectConstant"))),
+    cliqueMassGrow(readScalar(coeffsDict_.lookup("MassGrowCorrectConstant")))
+{
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::MossBrookesSoot::~MossBrookesSoot()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::MossBrookesSoot::a(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 (radiationActive_ == 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;
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::MossBrookesSoot::e(const label bandI) const
+{
+
+    return a(bandI);
+}
+
+void Foam::radiation::MossBrookesSoot::correct()
+{
+
+    calcSource();
+
+    const volScalarField rho_ = thermo_.rho();
+
+    const volScalarField &nut_ = turbulence().nut();
+
+    const surfaceScalarField phi_ (turbulence().phi());
+
+    fvScalarMatrix SootEqn
+    (
+        fvm::div(phi_, soot_)
+      - fvm::laplacian((rho_*nut_/PrtSoot), soot_)
+     ==
+        Snet_
+    );
+
+    SootEqn.relax();
+
+    SootEqn.solve();
+
+    fvScalarMatrix NucConcEqn
+    (
+        fvm::div(phi_, NucConc_)
+      - fvm::laplacian((rho_*nut_/PrtSoot), NucConc_)
+     ==
+        NSnet_
+    );
+
+    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;
+    }
+}
+
+
+void Foam::radiation::MossBrookesSoot::calcSource()
+{
+    const volScalarField rho_ = thermo_.rho();
+    const volScalarField T_ = thermo_.T();
+    const volScalarField p_ = thermo_.p();
+
+    const scalar precI = mixture_.species()[coeffsDict_.lookup("Precursor")];
+    const scalar growI = mixture_.species()[coeffsDict_.lookup("SurfaceGrowSpecie")];
+    const scalar oxidI = mixture_.species()[coeffsDict_.lookup("SootOxidizerSpecie")];
+
+    const volScalarField &Yprec_ = mixture_.Y(precI);
+    const volScalarField &Ygrow_ = mixture_.Y(growI);
+    const volScalarField &Yoxid_ = mixture_.Y(oxidI);
+
+    const volScalarField XPrec_(mixture_.W()*Yprec_/mixture_.W(precI));
+    const volScalarField XSgs_ (mixture_.W()*Ygrow_/mixture_.W(growI));
+    const volScalarField XOxid_(mixture_.W()*Yoxid_/mixture_.W(oxidI));
+
+    forAll(rho_, celli)
+    {
+        const scalar XPreci = XPrec_[celli];
+        const scalar XSgsi = XSgs_[celli];
+        const scalar XOxidi = XOxid_[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 = (XPreci*pi)/(8314.5*Ti);
+        const scalar MolarConcSgsi = (XSgsi*pi)/(8314.5*Ti);
+        const scalar MolarConcOxidi = (XOxidi*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);
+
+        aggParam_[celli] = 24.0*Boltzmann*Ti/rhoSoot;
+
+        NSagg_[celli] = (Cbeta*pow((aggParam_[celli]),0.5)*pow(SootMeanDiameteri,0.5)*pow(NumDeni,2.0))*(1e-15);
+
+        NSnet_[celli] = NSnuc_[celli] - NSagg_[celli];
+    }
+}
+
+
+const Foam::compressibleTurbulenceModel&
+Foam::radiation::MossBrookesSoot::turbulence() const
+{
+    return turbulence_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.H

@@ -0,0 +1,283 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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;
+    }
+
+SourceFiles
+    MossBrookesSoot.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MossBrookesSoot_H
+#define MossBrookesSoot_H
+
+#include "sootModel.H"
+#include "HashTable.H"
+
+#include "fluidThermo.H"
+#include "basicSpecieMixture.H"
+#include "turbulentFluidThermoModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class MossBrookesSoot Declaration
+\*---------------------------------------------------------------------------*/
+class MossBrookesSoot
+:
+    public sootModel
+
+{
+
+    // Static functions
+
+        //- Check mixture in thermo
+        static const basicSpecieMixture& checkThermo
+        (
+            const fluidThermo&
+        );
+
+
+    // 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_;
+
+        //- Composition
+        const basicSpecieMixture& mixture_;
+
+        //- Reference to the turbulence model
+        const compressibleTurbulenceModel &turbulence_;
+
+        //- 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 aggParam_;
+
+        //- 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 radiationActive_;
+
+        //- 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
+
+    // Protected Member Functions
+
+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_;
+            }
+
+            //- Return access to turbulence
+            const compressibleTurbulenceModel& turbulence() const;
+
+            // Soot absorption coefficient
+
+                //- Soot absorption coefficient (net)
+                virtual tmp<volScalarField> a(const label bandI = 0) const;
+
+            // Soot emission coefficient
+
+                //- Soot emission coefficient (net)
+                virtual tmp<volScalarField> e(const label bandI = 0) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiationSootModel/khanGreeveSoot/khanGreeveSoot.C

@@ -0,0 +1,374 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "khanGreeveSoot.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvm.H"
+#include "basicSpecieMixture.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(khanGreeveSoot, 0);
+
+        addToRunTimeSelectionTable
+        (
+            sootModel,
+            khanGreeveSoot,
+            dictionary
+        );
+    }
+}
+
+
+const Foam::basicSpecieMixture&
+Foam::radiation::khanGreeveSoot::checkThermo
+(
+    const fluidThermo& thermo
+)
+{
+    if (isA<basicSpecieMixture>(thermo))
+    {
+        return dynamic_cast<const basicSpecieMixture& >
+        (
+            thermo
+        );
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Inconsistent thermo package for " << thermo.type()
+            << "Please select a thermo package based on "
+            << "basicSpecieMixture" << exit(FatalError);
+
+        return dynamic_cast<const basicSpecieMixture& >
+        (
+            thermo
+        );
+    }
+
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::khanGreeveSoot::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)),
+    mixture_(checkThermo(thermo_)),
+    turbulence_(mesh.lookupObject<compressibleTurbulenceModel>(turbulenceModel::propertiesName)),
+    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"))),
+    radiationActive_(readBool(coeffsDict_.lookup("Radiationeffect")))
+{
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::khanGreeveSoot::~khanGreeveSoot()
+
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot::a(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("aSoot", dimless/dimLength, 0.0),
+            extrapolatedCalculatedFvPatchVectorField::typeName
+        )
+    );
+
+    scalarField& as = tas.ref().primitiveFieldRef();
+    if (radiationActive_)
+    {
+        forAll(as, celli)
+        {
+            as[celli] = (1232.4)*(1800.0)*soot_[celli]*(1.0+(4.8e-4)*(T_[celli]-2000.0));
+        }
+    }
+    tas.ref().correctBoundaryConditions();
+
+    return tas;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::khanGreeveSoot::e(const label bandI) const
+{
+
+    return a(bandI);
+}
+
+
+void Foam::radiation::khanGreeveSoot::correct()
+{
+
+    calcSource();
+
+    const volScalarField &rho_ = thermo_.rho();
+
+    const volScalarField &nut_ = turbulence().nut();
+
+    const surfaceScalarField phi_ (turbulence().phi());
+
+    fvScalarMatrix SootEqn
+    (
+        fvm::ddt(rho_, soot_)
+      + 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;
+
+    }
+
+}
+
+
+void Foam::radiation::khanGreeveSoot::calcSource()
+{
+    const volScalarField rho_ = thermo_.rho();
+    const volScalarField T_ = thermo_.T();
+    const volScalarField p_ = thermo_.p();
+
+    const scalar oxI = mixture_.species()[coeffsDict_.lookup("Oxidizer")];
+
+    const scalar fuelI = mixture_.species()[coeffsDict_.lookup("Fuel")];
+
+    const volScalarField &Yox_ = mixture_.Y(oxI);
+
+    const volScalarField &Yfuel_ = mixture_.Y(fuelI);
+
+    const volScalarField Xfuel_(mixture_.W()*Yfuel_/mixture_.W(fuelI));
+
+    const volScalarField &epsilon_ = turbulence().epsilon();
+
+    const volScalarField &k_ = turbulence().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 Xfueli = Xfuel_[celli];
+        const scalar Ti = T_[celli];
+
+        Sform_[celli] = Cs_*pi*Xfueli*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];
+    }
+
+}
+
+
+const Foam::compressibleTurbulenceModel&
+Foam::radiation::khanGreeveSoot::turbulence() const
+{
+    return turbulence_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiationSootModel/khanGreeveSoot/khanGreeveSoot.H

@@ -0,0 +1,230 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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;
+    }
+
+SourceFiles
+    khanGreeveSoot.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef khanGreeveSoot_H
+#define khanGreeveSoot_H
+
+#include "sootModel.H"
+#include "HashTable.H"
+
+#include "fluidThermo.H"
+#include "basicSpecieMixture.H"
+#include "turbulentFluidThermoModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class khanGreeveSoot Declaration
+\*---------------------------------------------------------------------------*/
+class khanGreeveSoot
+:
+    public sootModel
+
+{
+
+    // Static functions
+
+        //- Check mixture in thermo
+        static const basicSpecieMixture& checkThermo
+        (
+            const fluidThermo&
+        );
+
+
+    // Private data
+
+         //- Soot mass fraction
+        volScalarField soot_;
+
+        //- Soot volume fraction
+        volScalarField sootVF_;
+
+        //- Soot model dictionary
+        dictionary coeffsDict_;
+
+        //- Thermo package
+        const fluidThermo& thermo_;
+
+        //- Composition
+        const basicSpecieMixture& mixture_;
+
+        //- Reference to the turbulence model
+        const compressibleTurbulenceModel &turbulence_;
+
+        //- 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 radiationActive_;
+
+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_;
+            }
+
+            //- Return access to turbulence
+            const compressibleTurbulenceModel& turbulence() const;
+
+            // Soot absorption coefficient
+
+                //- Soot absorption coefficient (net)
+                virtual tmp<volScalarField> a(const label bandI = 0) const;
+
+            // Soot emission coefficient
+
+                //- Soot emission coefficient (net)
+                virtual tmp<volScalarField> e(const label bandI = 0) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace radiation
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //