ThermoCloud new constructor dynamic SLGThermo lookup

This reverts commit c0264dcdfa.

.gitignore

@@ -7,6 +7,8 @@
 *.bak
 *.bak[0-9][0-9]
 \#*\#
+.*.swp
+.*.swo
 
 # File-browser settings - anywhere
 .directory

.gitlab-ci.yml

@@ -1,23 +1,9 @@
+before_script:
+  - "[[ -d ../ThirdParty-4.x ]] || git clone https://github.com/OpenFOAM/ThirdParty-4.x.git ../ThirdParty-4.x"
+  - source etc/bashrc
+  - rm -rf platforms
+  - wclean all
+
 allmake:
-  image: park0d/of4builder
-  tags:
-    - openfoam4
-  before_script:
-    - "[[ -d ../ThirdParty-4.x ]] || git clone https://github.com/OpenFOAM/ThirdParty-4.x.git ../ThirdParty-4.x"
-    - source etc/bashrc
-    - rm -rf platforms
-    - wclean all
   script:
     - ./Allwmake -j
-  artifacts:
-    paths:
-      - ./
-
-cantera:
-    image: park0d/of4builder
-    tags:
-        - openfoam4
-    script:
-        - whereis libcantera
-    
-  

applications/solvers/lagrangian/ThermoMPPICFoam/EEqn.H

@@ -0,0 +1,40 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(alphac,rho, he) + mvConvection->fvmDiv(alphacPhi, he)
+      + fvc::ddt(alphac,rho, K) + fvc::div(alphacPhi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), U),
+                p,
+                "div(phiv,p)"
+            )
+          : -dpdt
+        )
+      - fvm::laplacian(alphac*turbulence->alphaEff(), he)
+     ==
+        alphac*rho*(U&g)
+      + parcels.Sh(he)
+      + radiation->Sh(thermo)
+      + combustion->Sh()
+      + fvOptions(alphac,rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/ThermoMPPICFoam/Make/files

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

applications/solvers/lagrangian/ThermoMPPICFoam/Make/options

@@ -0,0 +1,56 @@
+EXE_INC = \
+    -I. \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/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/chemistryModel/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)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lspecie \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lliquidProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lthermophysicalFunctions \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lradiationModels \
+    -lODE \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lcombustionModels \
+    -lfvOptions \
+    -lsampling

applications/solvers/lagrangian/ThermoMPPICFoam/ThermoMPPICFoam.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    ThermoMPPICFoam
+
+Description
+    Transient solver for compressible, turbulent flow with a reacting,
+    multiphase particle cloud, and optional sources/constraints.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulentFluidThermoModel.H"
+#include "basicThermoKinematicMPPICCloud.H"
+#include "rhoCombustionModel.H"
+#include "radiationModel.H"
+#include "fvOptions.H"
+#include "SLGThermo.H"
+#include "pimpleControl.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    turbulence->validate();
+
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        parcels.evolve();
+
+        // Update continuous phase volume fraction field
+        alphac = max(1.0 - parcels.theta(), alphacMin);
+        alphac.correctBoundaryConditions();
+        alphacf = fvc::interpolate(alphac);
+        alphacPhi = alphacf*phi;
+
+        fvVectorMatrix cloudSU(parcels.SU(U));
+        volVectorField cloudVolSUSu
+        (
+            IOobject
+            (
+                "cloudVolSUSu",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedVector
+            (
+                "0",
+                cloudSU.dimensions()/dimVolume,
+                vector::zero
+            ),
+            zeroGradientFvPatchVectorField::typeName
+        );
+
+        cloudVolSUSu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
+        cloudVolSUSu.correctBoundaryConditions();
+        cloudSU.source() = vector::zero;
+
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        rho = thermo.rho();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/ThermoMPPICFoam/UEqn.H

@@ -0,0 +1,43 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(alphac, rho, U)
+      + fvm::div(alphacPhi, U)
+      - fvm::laplacian(alphac*turbulence->muEff(),U)
+      - fvc::div(alphac*turbulence->muEff()*dev2(fvc::grad(U)().T()))
+      - fvm::Sp(fvc::ddt(alphac,rho) + fvc::div(alphacPhi), U)
+     ==
+      cloudSU
+      + fvOptions(alphac,rho, U)
+    );
+
+    UEqn.relax();
+
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+
+    surfaceScalarField phicForces
+    (
+        (fvc::interpolate(rAU*cloudVolSUSu) & mesh.Sf())
+        + alphacf*rhorAUf*(g & mesh.Sf())
+    );
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                phicForces/rAUf
+            )
+          - fvc::grad(p)
+        );
+
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
+    }

applications/solvers/lagrangian/ThermoMPPICFoam/YEqn.H

@@ -0,0 +1,53 @@
+tmp<fv::convectionScheme<scalar>> mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    combustion->correct();
+    dQ = combustion->dQ();
+    label inertIndex = -1;
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YEqn
+            (
+                fvm::ddt(alphac,rho, Yi)
+              + mvConvection->fvmDiv(alphacPhi, Yi)
+              - fvm::laplacian(alphac*turbulence->muEff(), Yi)
+             ==
+                combustion->R(Yi)
+              + fvOptions(alphac,rho, Yi)
+            );
+
+            YEqn.relax();
+
+            fvOptions.constrain(YEqn);
+
+            YEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}

applications/solvers/lagrangian/ThermoMPPICFoam/createClouds.H

@@ -0,0 +1,26 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicThermoKinematicMPPICCloud parcels
+(
+    "reactingCloud",
+    rho,
+    U,
+    slgThermo.thermo().mu(),
+    g
+);
+
+
+// Continuouse phase volume fraction lower limit
+scalar alphacMin
+(
+    1.0
+    - readScalar
+    (
+        parcels.particleProperties().subDict("constantProperties").lookup("alphaMax")
+    )
+);
+
+// Update alphac from the particle locations
+alphac = max(1.0 - parcels.theta(), alphacMin);
+alphac.correctBoundaryConditions();
+alphacf = fvc::interpolate(alphac);
+alphacPhi = alphacf*phi;

applications/solvers/lagrangian/ThermoMPPICFoam/createFieldRefs.H

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

applications/solvers/lagrangian/ThermoMPPICFoam/createFields.H

@@ -0,0 +1,158 @@
+#include "readGravitationalAcceleration.H"
+
+Info<< "Creating combustion model\n" << endl;
+
+autoPtr<combustionModels::rhoCombustionModel> combustion
+(
+    combustionModels::rhoCombustionModel::New(mesh)
+);
+
+rhoReactionThermo& thermo = combustion->thermo();
+thermo.validate(args.executable(), "h", "e");
+
+SLGThermo slgThermo(mesh, thermo);
+
+basicSpecieMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+const word inertSpecie(thermo.lookup("inertSpecie"));
+
+if (!composition.contains(inertSpecie))
+{
+    FatalErrorInFunction
+        << "Specified inert specie '" << inertSpecie << "' not found in "
+        << "species list. Available species:" << composition.species()
+        << exit(FatalError);
+}
+
+volScalarField& p = thermo.p();
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+mesh.setFluxRequired(p.name());
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        pimple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        pimple.dict(),
+        dimDensity,
+        0
+    )
+);
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+// Set the turbulence into the combustion 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));
+
+Info<< "Creating multi-variate interpolation scheme\n" << endl;
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+volScalarField dQ
+(
+    IOobject
+    (
+        "dQ",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+);
+
+volScalarField alphac
+(
+    IOobject
+    (
+        "alphac",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+        mesh,
+        dimensionedScalar("0", dimless, 1)
+);
+
+surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
+surfaceScalarField alphacPhi("alphacPhi", alphacf*phi);
+
+volScalarField alphacrho("alphacrho", alphac*rho);
+
+#include "createMRF.H"
+#include "createRadiationModel.H"
+#include "createClouds.H"

applications/solvers/lagrangian/ThermoMPPICFoam/pEqn.H

@@ -0,0 +1,81 @@
+{
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    rhorAUf= fvc::interpolate(rho*rAU);
+    surfaceScalarField rhof("rhof", fvc::interpolate(rho));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+         (fvc::interpolate(rho*HbyA) & mesh.Sf())
+         + rhorAUf*fvc::ddtCorr(rho, U, phi)
+         + phicForces*rhof
+
+        )
+    );
+
+       MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(alphac,rho) + alphac*psi*correction(fvm::ddt(p))
+      + fvc::div(alphacf*phiHbyA)
+     ==
+        fvOptions(alphac*psi, p, rho.name())
+    );
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            pDDtEqn
+          - fvm::laplacian(alphacf*rhorAUf, p)
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+
+           phi = phiHbyA + pEqn.flux()/alphacf;
+
+           alphacPhi = alphacf*phi;
+
+           p.relax();
+
+           thermo.rho() += psi*p;
+
+           #include "rhoEqn.H"
+           #include "compressibleContinuityErrs.H"
+
+           U = HbyA + rAU*fvc::reconstruct(phicForces/rAUf) - rAU*fvc::grad(p);
+
+         }
+    }
+
+    U.correctBoundaryConditions();
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+
+    if (thermo.dpdt())
+    {
+        dpdt = fvc::ddt(p);
+    }
+
+    rho = thermo.rho();
+    rho = max(rho, rhoMin);
+    rho = min(rho, rhoMax);
+
+    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+}

applications/solvers/lagrangian/ThermoMPPICFoam/rhoEqn.H

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(alphac,rho)
+      + fvc::div(alphacPhi)
+      ==
+       fvOptions(alphac,rho)
+    );
+
+    rhoEqn.solve();
+
+    fvOptions.correct(rho);
+
+    Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
+        << endl;
+
+    rho = max(rho, rhoMin);
+    rho = min(rho, rhoMax);
+}
+
+// ************************************************************************* //

applications/solvers/lagrangian/ThermoMPPICFoam/setRDeltaT.H

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    volScalarField& rDeltaT = trDeltaT.ref();
+
+    const dictionary& pimpleDict = pimple.dict();
+
+    // Maximum flow Courant number
+    scalar maxCo(readScalar(pimpleDict.lookup("maxCo")));
+
+    // Maximum time scale
+    scalar maxDeltaT(pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT));
+
+    // Smoothing parameter (0-1) when smoothing iterations > 0
+    scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    // Damping coefficient (1-0)
+    scalar rDeltaTDampingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 0.2)
+    );
+
+    // Maximum change in cell temperature per iteration
+    // (relative to previous value)
+    scalar alphaTemp(pimpleDict.lookupOrDefault("alphaTemp", 0.05));
+
+
+    Info<< "Time scales min/max:" << endl;
+
+    // Cache old reciprocal time scale field
+    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Flow time scale
+    {
+        rDeltaT.ref() =
+        (
+            fvc::surfaceSum(mag(phi))()()
+           /((2*maxCo)*mesh.V()*rho())
+        );
+
+        // Limit the largest time scale
+        rDeltaT.max(1/maxDeltaT);
+
+        Info<< "    Flow        = "
+            << gMin(1/rDeltaT.primitiveField()) << ", "
+            << gMax(1/rDeltaT.primitiveField()) << endl;
+    }
+
+    // Reaction source time scale
+    {
+        volScalarField::Internal rDeltaTT
+        (
+            mag
+            (
+               parcels.hsTrans()/(mesh.V()*runTime.deltaT())
+             + combustion->Sh()()
+            )
+           /(
+               alphaTemp
+              *rho()
+              *thermo.Cp()()()
+              *T()
+           )
+        );
+
+        Info<< "    Temperature = "
+            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
+            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
+
+        rDeltaT.ref() = max
+        (
+            rDeltaT(),
+            rDeltaTT
+        );
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    // Spatially smooth the time scale field
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+    }
+
+    Info<< "    Overall     = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+}
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/simpleCoalcombustionFoam/EEqn.H

@@ -0,0 +1,32 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        mvConvection->fvmDiv(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        parcels.Sh(he)
+      + radiation->Sh(thermo)
+      + combustion->Sh()
+      + fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/files

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/options

@@ -0,0 +1,57 @@
+EXE_INC = \
+    -I. \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/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/chemistryModel/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)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lcoalCombustion \
+    -llagrangianTurbulence \
+    -lspecie \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lliquidProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lthermophysicalFunctions \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lradiationModels \
+    -lODE \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lcombustionModels \
+    -lfvOptions \
+    -lsampling

applications/solvers/lagrangian/simpleCoalcombustionFoam/UEqn.H

@@ -0,0 +1,21 @@
+    MRF.correctBoundaryVelocity(U);
+
+    tmp<fvVectorMatrix> tUEqn
+    (
+        fvm::div(phi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho()*g
+      + parcels.SU(U)
+      + fvOptions(rho, U)
+    );
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    solve(UEqn == -fvc::grad(p));
+
+    fvOptions.correct(U);

applications/solvers/lagrangian/simpleCoalcombustionFoam/YEqn.H

@@ -0,0 +1,53 @@
+tmp<fv::convectionScheme<scalar>> mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    combustion->correct();
+    dQ = combustion->dQ();
+    label inertIndex = -1;
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YEqn
+            (
+                mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+             ==
+                parcels.Srho(i)
+              + combustion->R(Yi)
+              + fvOptions(rho, Yi)
+            );
+
+            YEqn.relax();
+
+            fvOptions.constrain(YEqn);
+
+            YEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}

applications/solvers/lagrangian/simpleCoalcombustionFoam/createClouds.H

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/createFieldRefs.H

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/createFields.H

@@ -0,0 +1,123 @@
+#include "readGravitationalAcceleration.H"
+
+Info<< "Creating combustion model\n" << endl;
+
+autoPtr<combustionModels::rhoCombustionModel> combustion
+(
+    combustionModels::rhoCombustionModel::New(mesh)
+);
+
+rhoReactionThermo& thermo = combustion->thermo();
+thermo.validate(args.executable(), "h", "e");
+
+SLGThermo slgThermo(mesh, thermo);
+
+basicSpecieMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+const word inertSpecie(thermo.lookup("inertSpecie"));
+
+if (!composition.contains(inertSpecie))
+{
+    FatalErrorInFunction
+        << "Specified inert specie '" << inertSpecie << "' not found in "
+        << "species list. Available species:" << composition.species()
+        << exit(FatalError);
+}
+
+volScalarField& p = thermo.p();
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+mesh.setFluxRequired(p.name());
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        simple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        simple.dict(),
+        dimDensity,
+        0
+    )
+);
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+// Set the turbulence into the combustion model
+combustion->setTurbulence(turbulence());
+
+Info<< "Creating multi-variate interpolation scheme\n" << endl;
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+volScalarField dQ
+(
+    IOobject
+    (
+        "dQ",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+);
+
+#include "createMRF.H"
+#include "createRadiationModel.H"
+#include "createClouds.H"

applications/solvers/lagrangian/simpleCoalcombustionFoam/pEqn.H

@@ -0,0 +1,57 @@
+{
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    tUEqn.clear();
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*fvc::flux(HbyA)
+    );
+
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUf, p)
+         ==
+            parcels.Srho()
+          + fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+
+    p.relax();
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    #include "compressibleContinuityErrs.H"
+
+    U = HbyA - rAU*fvc::grad(p);
+    U.correctBoundaryConditions();
+    fvOptions.correct(U);
+
+    rho = thermo.rho();
+    rho = max(rho, rhoMin);
+    rho = min(rho, rhoMax);
+    rho.relax();
+
+    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+}

applications/solvers/lagrangian/simpleCoalcombustionFoam/simpleCoalcombustionFoam.C

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    simpleCoalcombustionFoam
+
+Description
+    Steady state solver for compressible, turbulent flow with reacting,
+    multiphase particle clouds and optional sources/constraints.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulentFluidThermoModel.H"
+#include "coalCloudList.H"
+#include "rhoCombustionModel.H"
+#include "radiationModel.H"
+#include "IOporosityModelList.H"
+#include "fvOptions.H"
+#include "SLGThermo.H"
+#include "simpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    turbulence->validate();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        parcels.evolve();
+
+        // --- Pressure-velocity SIMPLE corrector loop
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "EEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H

@@ -34,6 +34,7 @@ License
 #include "COxidationKineticDiffusionLimitedRate.H"
 #include "COxidationHurtMitchell.H"
 #include "COxidationMurphyShaddix.H"
+#include "COxidationKD_CO.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,6 +48,7 @@ License
     );                                                                         \
     makeSurfaceReactionModelType(COxidationIntrinsicRate, CloudType);          \
     makeSurfaceReactionModelType(COxidationHurtMitchell, CloudType);           \
+    makeSurfaceReactionModelType(COxidationKD_CO, CloudType);                  \
     makeSurfaceReactionModelType(COxidationMurphyShaddix, CloudType);
 
 

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKD_CO/COxidationKD_CO.C

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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 "COxidationKD_CO.H"
+#include "mathematicalConstants.H"
+#include "thermodynamicConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationKD_CO<CloudType>::
+COxidationKD_CO
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    SurfaceReactionModel<CloudType>(dict, owner, typeName),
+    Sb_(readScalar(this->coeffDict().lookup("Sb"))),
+    C1_(readScalar(this->coeffDict().lookup("C1"))),
+    C2_(readScalar(this->coeffDict().lookup("C2"))),
+    E_(readScalar(this->coeffDict().lookup("E"))),
+    CsLocalId_(-1),
+    O2GlobalId_(owner.composition().carrierId("O2")),
+    COGlobalId_(owner.composition().carrierId("CO")),
+    WC_(0.0),
+    WO2_(0.0),
+    HcCO_(0.0)
+{
+    // Determine Cs ids
+    label idSolid = owner.composition().idSolid();
+    CsLocalId_ = owner.composition().localId(idSolid, "C");
+
+    // Set local copies of thermo properties
+    WO2_ = owner.thermo().carrier().W(O2GlobalId_);
+    const scalar WCO = owner.thermo().carrier().W(COGlobalId_);
+    WC_ = WCO - 0.5*WO2_;
+
+    HcCO_ = owner.thermo().carrier().Hc(COGlobalId_);
+
+    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
+    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
+    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
+}
+
+
+template<class CloudType>
+Foam::COxidationKD_CO<CloudType>::
+COxidationKD_CO
+(
+    const COxidationKD_CO<CloudType>& srm
+)
+:
+    SurfaceReactionModel<CloudType>(srm),
+    Sb_(srm.Sb_),
+    C1_(srm.C1_),
+    C2_(srm.C2_),
+    E_(srm.E_),
+    CsLocalId_(srm.CsLocalId_),
+    O2GlobalId_(srm.O2GlobalId_),
+    COGlobalId_(srm.COGlobalId_),
+    WC_(srm.WC_),
+    WO2_(srm.WO2_),
+    HcCO_(srm.HcCO_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationKD_CO<CloudType>::
+~COxidationKD_CO()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::COxidationKD_CO<CloudType>::calculate
+(
+    const scalar dt,
+    const label cellI,
+    const scalar d,
+    const scalar T,
+    const scalar Tc,
+    const scalar pc,
+    const scalar rhoc,
+    const scalar mass,
+    const scalarField& YGas,
+    const scalarField& YLiquid,
+    const scalarField& YSolid,
+    const scalarField& YMixture,
+    const scalar N,
+    scalarField& dMassGas,
+    scalarField& dMassLiquid,
+    scalarField& dMassSolid,
+    scalarField& dMassSRCarrier
+) const
+{
+    // Fraction of remaining combustible material
+    const label idSolid = CloudType::parcelType::SLD;
+    const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
+
+    // Surface combustion active combustible fraction is consumed
+    if (fComb < SMALL)
+    {
+        return 0.0;
+    }
+
+    const SLGThermo& thermo = this->owner().thermo();
+
+    // Local mass fraction of O2 in the carrier phase
+    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
+
+    // Diffusion rate coefficient
+    const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
+
+    // Kinetic rate
+    const scalar Rk = C2_*exp(-E_/(constant::thermodynamic::RR*T)); //
+
+    // Particle surface area
+    const scalar Ap = constant::mathematical::pi*sqr(d);
+
+    // Change in C mass [kg]
+    scalar dmC = Ap*rhoc*constant::thermodynamic::RR*Tc*YO2/WO2_*D0*Rk/(D0 + Rk)*dt;
+
+    // Limit mass transfer by availability of C
+    dmC = min(mass*fComb, dmC);
+
+    // Molar consumption
+    const scalar dOmega = dmC/WC_;
+
+    // Change in O2 mass [kg]
+    const scalar dmO2 = dOmega*Sb_*WO2_;
+
+    // Mass of newly created CO [kg]
+    const scalar dmCO = dOmega*(WC_ + Sb_*WO2_);
+
+    // Update local particle C mass
+    dMassSolid[CsLocalId_] += dOmega*WC_;
+
+    // Update carrier O2 and CO mass
+    dMassSRCarrier[O2GlobalId_] -= dmO2;
+    dMassSRCarrier[COGlobalId_] += dmCO;
+
+    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
+
+    // carrier sensible enthalpy exchange handled via change in mass
+
+    // Heat of reaction [J]
+//Info<<"Heat of reaction(Char)"<<tab<<(dmC*HsC - dmCO*HcCO_)/dmC<<endl;
+    return dmC*HsC - dmCO*HcCO_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKD_CO/COxidationKD_CO.H

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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
+    COxidationKD_CO
+
+Description
+    Kinetic/diffusion limited rate surface reaction model for coal parcels.
+    Limited to:
+
+        C(s) + Sb*O2 -> CO
+
+    where Sb is the stoichiometry of the reaction
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef COxidationKD_CO_H
+#define COxidationKD_CO_H
+
+#include "SurfaceReactionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+template<class CloudType>
+class COxidationKD_CO;
+
+/*---------------------------------------------------------------------------*\
+            Class COxidationKD_CO Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class COxidationKD_CO
+:
+    public SurfaceReactionModel<CloudType>
+{
+    // Private data
+
+        // Model constants
+
+            //- Stoichiometry of reaction
+            const scalar Sb_;
+
+            //- Mass diffusion limited rate constant, C1
+            const scalar C1_;
+
+            //- Kinetics limited rate pre-exponential constant, C2
+            const scalar C2_;
+
+            //- Kinetics limited rate activation energy
+            const scalar E_;
+
+
+        // Addressing
+
+            //- Cs positions in global/local lists
+            label CsLocalId_;
+
+            //- O2 position in global list
+            label O2GlobalId_;
+
+            //- CO2 positions in global list
+            label COGlobalId_;
+
+
+        // Local copies of thermo properties
+
+            //- Molecular weight of C [kg/kmol]
+            scalar WC_;
+
+            //- Molecular weight of O2 [kg/kmol]
+            scalar WO2_;
+
+            //- Formation enthalpy for CO2 [J/kg]
+            scalar HcCO_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("COxidationKD_CO");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        COxidationKD_CO
+        (
+            const dictionary& dict,
+            CloudType& owner
+        );
+
+        //- Construct copy
+        COxidationKD_CO
+        (
+            const COxidationKD_CO<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
+        {
+            return autoPtr<SurfaceReactionModel<CloudType> >
+            (
+                new COxidationKD_CO<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~COxidationKD_CO();
+
+
+    // Member Functions
+
+        //- Update surface reactions
+        virtual scalar calculate
+        (
+            const scalar dt,
+            const label cellI,
+            const scalar d,
+            const scalar T,
+            const scalar Tc,
+            const scalar pc,
+            const scalar rhoc,
+            const scalar mass,
+            const scalarField& YGas,
+            const scalarField& YLiquid,
+            const scalarField& YSolid,
+            const scalarField& YMixture,
+            const scalar N,
+            scalarField& dMassGas,
+            scalarField& dMassLiquid,
+            scalarField& dMassSolid,
+            scalarField& dMassSRCarrier
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "COxidationKD_CO.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/Make/files

@@ -50,6 +50,12 @@ $(KINEMATICMPPICPARCEL)/defineBasicKinematicMPPICParcel.C
 $(KINEMATICMPPICPARCEL)/makeBasicKinematicMPPICParcelSubmodels.C
 
 
+/* thermo kinematic MPPIC parcel sub-models */
+THERMOMPPICPARCEL=$(DERIVEDPARCELS)/basicThermoKinematicMPPICParcel
+$(THERMOMPPICPARCEL)/defineBasicThermoKinematicMPPICParcel.C
+$(THERMOMPPICPARCEL)/makeBasicThermoKinematicMPPICParcelSubmodels.C
+
+
 /* bolt-on models */
 RADIATION=submodels/addOns/radiation
 $(RADIATION)/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -213,6 +213,88 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
 }
 
 
+template<class CloudType>
+Foam::ThermoCloud<CloudType>::ThermoCloud
+(
+    const word& cloudName,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const volScalarField& mu,
+    const dimensionedVector& g,
+    bool readFields
+)
+:
+    CloudType
+    (
+        cloudName,
+        rho,
+        U,
+        mu,
+        g,
+        false
+    ),
+    thermoCloud(),
+    cloudCopyPtr_(NULL),
+    constProps_(this->particleProperties()),
+    thermo_(this->mesh().template lookupObject<SLGThermo>("SLGThermo")),
+    T_(thermo_.thermo().T()),
+    p_(thermo_.thermo().p()),
+    heatTransferModel_(NULL),
+    TIntegrator_(NULL),
+    radiation_(false),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
+    hsTrans_
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                this->name() + ":hsTrans",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimEnergy, 0.0)
+        )
+    ),
+    hsCoeff_
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                this->name() + ":hsCoeff",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimEnergy/dimTemperature, 0.0)
+        )
+    )
+{
+    if (this->solution().active())
+    {
+        setModels();
+
+        if (readFields)
+        {
+            parcelType::readFields(*this);
+        }
+    }
+
+    if (this->solution().resetSourcesOnStartup())
+    {
+        resetSourceTerms();
+    }
+}
+
+
 template<class CloudType>
 Foam::ThermoCloud<CloudType>::ThermoCloud
 (

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -180,6 +180,17 @@ public:
             bool readFields = true
         );
 
+        //- Construct given carrier gas fields
+        ThermoCloud
+        (
+            const word& cloudName,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const volScalarField& mu,
+            const dimensionedVector& g,
+            bool readFields = true
+        );
+
         //- Copy constructor with new name
         ThermoCloud(ThermoCloud<CloudType>& c, const word& name);
 

src/lagrangian/intermediate/clouds/derived/basicThermoKinematicMPPICCloud/basicThermoKinematicMPPICCloud.H

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

src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/basicThermoKinematicMPPICParcel.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::basicThermoKinematicMPPICParcel
+
+Description
+    Definition of basic Thermo kinematic MPPIC parcel
+
+SourceFiles
+    basicThermoKinematicMPPICParcel.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicThermoKinematicMPPICParcel_H
+#define basicThermoKinematicMPPICParcel_H
+
+#include "contiguous.H"
+#include "particle.H"
+#include "KinematicParcel.H"
+#include "MPPICParcel.H"
+
+#include "ThermoParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+    typedef MPPICParcel < ThermoParcel < KinematicParcel < particle > > > basicThermoKinematicMPPICParcel;
+
+    template<>
+    inline bool contiguous<basicThermoKinematicMPPICParcel>()
+    {
+        return true;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/defineBasicThermoKinematicMPPICParcel.C

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

src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/makeBasicThermoKinematicMPPICParcelSubmodels.C

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "basicThermoKinematicMPPICCloud.H"
+
+#include "makeParcelCloudFunctionObjects.H"
+
+// Kinematic sub-models
+#include "makeParcelForces.H"
+#include "makeParcelDispersionModels.H"
+#include "makeParcelInjectionModels.H"
+#include "makeParcelPatchInteractionModels.H"
+#include "makeParcelStochasticCollisionModels.H"
+#include "makeParcelSurfaceFilmModels.H"
+
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
+// Thermo sub-models
+#include "makeParcelHeatTransferModels.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeParcelCloudFunctionObjects(basicThermoKinematicMPPICCloud);
+
+// Kinematic sub-models
+makeParcelForces(basicThermoKinematicMPPICCloud);
+
+makeParcelDispersionModels(basicThermoKinematicMPPICCloud);
+makeParcelInjectionModels(basicThermoKinematicMPPICCloud);
+makeParcelPatchInteractionModels(basicThermoKinematicMPPICCloud);
+makeParcelStochasticCollisionModels(basicThermoKinematicMPPICCloud);
+makeParcelSurfaceFilmModels(basicThermoKinematicMPPICCloud);
+
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicThermoKinematicMPPICCloud);
+makeMPPICParcelIsotropyModels(basicThermoKinematicMPPICCloud);
+makeMPPICParcelPackingModels(basicThermoKinematicMPPICCloud);
+
+// Thermo sub-models
+makeParcelHeatTransferModels(basicThermoKinematicMPPICCloud);
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/Make/files

@@ -25,6 +25,7 @@ submodels/absorptionEmissionModel/binaryAbsorptionEmission/binaryAbsorptionEmiss
 submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
 submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
 submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
+submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.C
 
 /* Soot model */
 submodels/sootModel/sootModel/sootModel.C

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.C

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "wsggmAbsorptionEmission.H"
+#include "addToRunTimeSelectionTable.H"
+#include "unitConversion.H"
+#include "zeroGradientFvPatchFields.H"
+#include "basicSpecieMixture.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(wsggmAbsorptionEmission, 0);
+
+        addToRunTimeSelectionTable
+        (
+            absorptionEmissionModel,
+            wsggmAbsorptionEmission,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::wsggmAbsorptionEmission::wsggmAbsorptionEmission
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    absorptionEmissionModel(dict, mesh),
+    coeffsDict_((dict.subDict(typeName + "Coeffs"))),
+    thermo_(mesh.lookupObject<fluidThermo>("thermophysicalProperties")),
+    emissivityCoeffs_(coeffsDict_.lookup("emissivityCoeffs")),
+    fittingFactors_(coeffsDict_.lookup("fittingFactors")),
+    pathLength_(coeffsDict_.lookup("pathLength"))
+{
+    if (!isA<basicMultiComponentMixture>(thermo_))
+    {
+        FatalErrorIn
+        (
+            "radiation::wsggmAbsorptionEmission::wsggmAbsorptionEmission"
+            "("
+                "const dictionary&, "
+                "const fvMesh&"
+            ")"
+        )   << "Model requires a multi-component thermo package"
+            << abort(FatalError);
+    }
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::wsggmAbsorptionEmission::~wsggmAbsorptionEmission()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wsggmAbsorptionEmission::aCont(const label bandI) const
+{
+
+    const basicSpecieMixture& mixture =
+        dynamic_cast<const basicSpecieMixture&>(thermo_);
+
+    const volScalarField& T = thermo_.T();
+    const volScalarField& p = thermo_.p();
+
+    tmp<volScalarField> ta
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "aCont" + name(bandI),
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("a", dimless/dimLength, 0.0),
+            zeroGradientFvPatchVectorField::typeName
+        )
+    );
+
+    scalarField& a = ta.ref().primitiveFieldRef();
+
+    label indexCO2, indexH2O;
+    scalar wCO2, wH2O;
+    forAll(mixture.Y(), specieI)
+    {
+        if(mixture.Y(specieI).name() == "CO2")
+        {
+            indexCO2 = specieI;
+            wCO2 = mixture.W(indexCO2);
+        }
+        else if(mixture.Y(specieI).name() == "H2O")
+        {
+            indexH2O = specieI;
+            wH2O = mixture.W(indexH2O);
+        }
+    }
+
+    forAll(a,celli)
+    {
+        scalar invWt = 0.0;
+        forAll(mixture.Y(), specieI)
+        {
+            invWt += mixture.Y(specieI)[celli]/mixture.W(specieI);
+        }
+        scalar meanWt = 1.0/invWt; //Mean molecular weight [kg/kmole]
+
+        // Sum of partial pressures of all absorbing gases [atm]
+        scalar Pcell_atm =  paToAtm(meanWt*p[celli]*(mixture.Y(indexCO2)[celli]/wCO2 + mixture.Y(indexH2O)[celli]/wH2O));
+
+        // P*S at Eq.(1) of Ref. [atm*m]
+        scalar presPathLength = Pcell_atm*pathLength_.value();
+
+        // Limit cell temperature [K]
+        scalar Tcell = min(T[celli],2400.0);
+
+        scalar emissivity = 0.0;
+        forAll(emissivityCoeffs_, i) //k_i at Eq.(1) of Ref.
+        {
+            scalar weightingFactor = 0.0;
+            for(label j=0 ; j<4 ; j++)
+            {
+                //a_epsilon_i at Eq.(3) of Ref.
+                weightingFactor += fittingFactors_[i][j]*pow(Tcell, j);
+            }
+            //epsilon at Eq.(1) of Ref.
+            emissivity += weightingFactor * (1.0 - exp(-1.0*emissivityCoeffs_[i] * presPathLength));
+        }
+        emissivity = min(emissivity, 0.99999);
+
+        //Fluent theory guide(ver. 13) Eq.(5-106) [1/m]
+        a[celli] = max(1e-2, (-1.0)*log(1.0-emissivity)/pathLength_.value());
+    }
+
+    ta.ref().correctBoundaryConditions();
+
+    return ta;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wsggmAbsorptionEmission::eCont(const label bandI) const
+{
+   return aCont(bandI);
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::wsggmAbsorptionEmission::ECont(const label bandI) const
+{
+    tmp<volScalarField> E
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "ECont" + name(bandI),
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
+        )
+    );
+
+    return E;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.H

@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 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::wsggmAbsorptionEmission
+
+Description
+    wsggmAbsorptionEmission radiation absorption and emission coefficients
+    for continuous phase
+
+    Ref. T.F Smith, Z.F. Shen and J.N. Friedman
+         Evaluation of Coefficients for the Weighted Sum of Gray Gases Model,
+         Transactions of the ASME, Vol. 104, pp. 602-608
+
+    Valid when Ptotal = 1atm, 0.001 < P*S < 10.0atm.m, 600 < T < 2400K
+               (where P*S = sum of partial pressures of absorbing gases
+                            * mean beam length)
+
+    i.e. dictionary
+    \verbatim
+        wsggmAbsorptionEmissionCoeffs
+        {
+            pathLength pathLength [ 0 1  0  0  0  0  0] 0.251; //3.6V/S
+            // Pw/Pc = 2.0
+            emissivityCoeffs 3(0.4201 6.516 131.9);
+            fittingFactors
+            3
+            (
+                4(6.508e-1 -5.551e-4 3.029e-7  -5.353e-11)
+                4(-0.2504e-1  6.112e-4  -3.882e-7  6.528e-11)
+                4(2.718e-1 -3.118e-4  1.221e-7  -1.612e-11)
+            );
+        }
+    \endverbatim
+
+
+SourceFiles
+    wsggmAbsorptionEmission.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wsggmAbsorptionEmission_H
+#define wsggmAbsorptionEmission_H
+
+#include "interpolationLookUpTable.H"
+#include "absorptionEmissionModel.H"
+#include "HashTable.H"
+#include "absorptionCoeffs.H"
+#include "fluidThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class wsggmAbsorptionEmission Declaration
+\*---------------------------------------------------------------------------*/
+
+class wsggmAbsorptionEmission
+:
+    public absorptionEmissionModel
+{
+public:
+
+private:
+
+    // Private data
+
+    //- Absorption model dictionary
+    dictionary coeffsDict_;
+
+    //- SLG thermo package
+    const fluidThermo& thermo_;
+
+    //- absorption coefficitent for individual grey gases
+    scalarList emissivityCoeffs_;
+
+    //- fitting Factors
+    scalarListList fittingFactors_;
+
+    //characteristic path length [m] from dict
+    dimensionedScalar pathLength_;
+
+public:
+
+    //- Runtime type information
+    TypeName("wsggmAbsorptionEmission");
+
+
+    // Constructors
+
+        //- Construct from components
+        wsggmAbsorptionEmission(const dictionary& dict, const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~wsggmAbsorptionEmission();
+
+
+    // Member Functions
+
+        // Access
+
+            // Absorption coefficient
+
+                //- Absorption coefficient for continuous phase
+                tmp<volScalarField> aCont(const label bandI = 0) const;
+
+
+            // Emission coefficient
+
+                //- Emission coefficient for continuous phase
+                tmp<volScalarField> eCont(const label bandI = 0) const;
+
+
+            // Emission contribution
+
+                //- Emission contribution for continuous phase
+                tmp<volScalarField> ECont(const label bandI = 0) const;
+
+
+    // Member Functions
+
+        inline bool isGrey() const
+        {
+            return true;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L

@@ -148,7 +148,7 @@ irreversibleReactionDelimiter {space}"=>"{space}
 startPDependentSpecie {space}"("{space}"+"{space}
 pDependentSpecie      {specieName}")"{space}
 reactionCoeffs        {space}{floatNum}{some_space}{floatNum}{some_space}{floatNum}{space}
-reactionKeyword       {space}[A-Za-z](([A-Za-z0-9)*-])|("("[^+]))*((\+({some_space}))?){space}
+reactionKeyword       {space}[A-Za-z](([A-Za-z0-9)*-])|("("[^+]))*{space}
 reactionKeywordSlash  {reactionKeyword}"/"{space}
 thirdBodyEfficiency   {space}{floatNum}{space}"/"{space}
 startReactionCoeffs   {space}"/"{space}
@@ -637,8 +637,7 @@ bool finishReaction = false;
                     (
                         currentSpecieName,
                         1.0,
-                        molecularWeight(currentSpecieComposition),
-                        chargeNumber(currentSpecieComposition)
+                        molecularWeight(currentSpecieComposition)
                     ),
                     currentLowT,
                     currentHighT,

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C

@@ -144,28 +144,6 @@ Foam::scalar Foam::chemkinReader::molecularWeight
 }
 
 
-Foam::scalar Foam::chemkinReader::chargeNumber
-(
-    const List<specieElement>& specieComposition
-) const
-{
-    scalar nElemCharges = 0.0;
-
-    forAll(specieComposition, i)
-    {
-        label nAtoms = specieComposition[i].nAtoms;
-        const word& elementName = specieComposition[i].elementName;
-
-        if (elementName[0] == 'e')
-        {
-            nElemCharges -= nAtoms;
-        }
-    }
-
-    return nElemCharges;
-}
-
-
 void Foam::chemkinReader::checkCoeffs
 (
     const scalarList& reactionCoeffs,

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.H

@@ -254,11 +254,6 @@ private:
             const List<specieElement>& specieComposition
         ) const;
 
-        scalar chargeNumber
-        (
-            const List<specieElement>& specieComposition
-        ) const;
-
         void finishElements(labelList& currentAtoms);
 
         void checkCoeffs

src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C

@@ -68,70 +68,6 @@ Foam::scalar Foam::SpecieMixture<MixtureType>::W
 }
 
 
-template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::z
-(
-    const label speciei
-) const
-{
-    return this->getLocalThermo(speciei).z();
-}
-
-
-template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::Qc
-(
-    const label speciei
-) const
-{
-    return this->getLocalThermo(speciei).Qc();
-}
-
-
-template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::Qc2
-(
-    const label speciei
-) const
-{
-    return this->getLocalThermo(speciei).Qc2();
-}
-
-
-template<class MixtureType>
-void Foam::SpecieMixture<MixtureType>::Cp
-(
-    scalarField& Cps,
-    const scalar p,
-    const scalar T
-) const
-{
-    forAll(this->species(), speciei)
-    {
-        Cps[speciei] = this->getLocalThermo(speciei).Cp(p, T);
-    }
-
-    return;
-}
-
-
-template<class MixtureType>
-void Foam::SpecieMixture<MixtureType>::Cv
-(
-    scalarField& Cvs,
-    const scalar p,
-    const scalar T
-) const
-{
-    forAll(this->species(), speciei)
-    {
-        Cvs[speciei] = this->getLocalThermo(speciei).Cv(p, T);
-    }
-
-    return;
-}
-
-
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::Cp
 (

src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H

@@ -81,34 +81,6 @@ public:
             //- Molecular weight of the given specie [kg/kmol]
             virtual scalar W(const label speciei) const;
 
-            //- Number of charges of the given specie []
-            virtual scalar z(const label speciei) const;
-
-            //- Specific charge of the given specie []
-            virtual scalar Qc(const label specieI) const;
-
-            //- Absolute specific charge of the given specie []
-            virtual scalar Qc2(const label specieI) const;
-
-
-        // All species thermo properties
-
-            //- Heat capacity at constant pressure [J/(kg K)]
-            virtual void Cp
-            (
-                scalarField& Cps,
-                const scalar p,
-                const scalar T
-            ) const;
-
-            //- Heat capacity at constant volume [J/(kg K)]
-            virtual void Cv
-            (
-                scalarField& Cvs,
-                const scalar p,
-                const scalar T
-            ) const;
-
 
         // Per specie thermo properties
 

src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.C

@@ -81,72 +81,4 @@ Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::W() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::Qc() const
-{
-    const PtrList<volScalarField>& Y(basicMultiComponentMixture::Y());
-
-    tmp<volScalarField> tQc
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName("Qc", Y[0].group()),
-                Y[0].time().timeName(),
-                Y[0].mesh()
-            ),
-            Y[0].mesh(),
-            dimensionedScalar("zero", dimCurrent*dimTime/dimMass, 0)
-        )
-    );
-
-    volScalarField& Qc = tQc.ref();
-
-    dimensionedScalar Qci("Qci", dimCurrent*dimTime/dimMass, 0);
-
-    forAll(Y, i)
-    {
-        Qci.value() = this->Qc(i);
-        Qc += Y[i]*Qci;
-        // Qc += Y[i]*this->Qc(i);
-    }
-
-    return tQc;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::Qc2() const
-{
-    const PtrList<volScalarField>& Y(basicMultiComponentMixture::Y());
-
-    tmp<volScalarField> tQc
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName("Qc2", Y[0].group()),
-                Y[0].time().timeName(),
-                Y[0].mesh()
-            ),
-            Y[0].mesh(),
-            dimensionedScalar("zero", sqr(dimCurrent*dimTime)/dimMass, 0)
-        )
-    );
-
-    volScalarField& Qc = tQc.ref();
-
-    dimensionedScalar Qc2i("Qc2i", sqr(dimCurrent*dimTime)/dimMass, 0);
-
-    forAll(Y, i)
-    {
-        Qc2i.value() = this->Qc2(i);
-        Qc += Y[i]*Qc2i;
-        //Qc += Y[i]*this->Qc2(i);
-    }
-
-    return tQc;
-}
-
-
 // ************************************************************************* //

src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H

@@ -88,43 +88,9 @@ public:
             //- Molecular weight of the given specie [kg/kmol]
             virtual scalar W(const label speciei) const = 0;
 
-            //- Number of charges of the given specie []
-            virtual scalar z(const label speciei) const = 0;
-
-            //- Specific charge of the given specie [C/kg]
-            virtual scalar Qc(const label speciei) const = 0;
-
-            //- Absolute specific charge of the given specie [C/kg]
-            virtual scalar Qc2(const label speciei) const = 0;
-
             //- Molecular weight of the mixture [kg/kmol]
             tmp<volScalarField> W() const;
 
-            //- Specific charge of the mixture [C/kg]
-            tmp<volScalarField> Qc() const;
-
-            //- Absolute specific charge of the mixture [C^2/kg]
-            tmp<volScalarField> Qc2() const;
-
-
-        // All species thermo properties
-
-            //- Heat capacity at constant pressure [J/(kg K)]
-            virtual void Cp
-            (
-                scalarField& Cps,
-                const scalar p,
-                const scalar T
-            ) const = 0;
-
-            //- Heat capacity at constant volume [J/(kg K)]
-            virtual void Cv
-            (
-                scalarField& Cvs,
-                const scalar p,
-                const scalar T
-            ) const = 0;
-
 
         // Per specie thermo properties
 

src/thermophysicalModels/specie/specie/specie.C

@@ -29,24 +29,6 @@ License
 
 /* * * * * * * * * * * * * * * public constants  * * * * * * * * * * * * * * */
 
-//- Universal gas constant (default in [J/(kmol K)])
-const Foam::scalar Foam::specie::RR = constant::physicoChemical::R.value()*1000;
-
-//- Standard pressure (default in [Pa])
-const Foam::scalar Foam::specie::Pstd = constant::standard::Pstd.value();
-
-//- Standard temperature (default in [K])
-const Foam::scalar Foam::specie::Tstd = constant::standard::Tstd.value();
-
-//- Elementary charge (default in [C])
-const Foam::scalar Foam::specie::e = constant::electromagnetic::e.value();
-
-//- Avogadro number (default in [1/mol])
-const Foam::scalar Foam::specie::NA = constant::physicoChemical::NA.value()*1000;
-
-//- Boltzmann constant (default in [J/K])
-const Foam::scalar Foam::specie::k = constant::physicoChemical::k.value();
-
 namespace Foam
 {
     defineTypeNameAndDebug(specie, 0);
@@ -59,8 +41,7 @@ Foam::specie::specie(Istream& is)
 :
     name_(is),
     nMoles_(readScalar(is)),
-    molWeight_(readScalar(is)),
-    nCharges_(readScalar(is))
+    molWeight_(readScalar(is))
 {
     is.check("specie::specie(Istream& is)");
 }
@@ -70,8 +51,7 @@ Foam::specie::specie(const dictionary& dict)
 :
     name_(dict.dictName()),
     nMoles_(readScalar(dict.subDict("specie").lookup("nMoles"))),
-    molWeight_(readScalar(dict.subDict("specie").lookup("molWeight"))),
-    nCharges_(dict.subDict("specie").lookupOrDefault("nCharges", 0.0))
+    molWeight_(readScalar(dict.subDict("specie").lookup("molWeight")))
 {}
 
 
@@ -82,7 +62,6 @@ void Foam::specie::write(Ostream& os) const
     dictionary dict("specie");
     dict.add("nMoles", nMoles_);
     dict.add("molWeight", molWeight_);
-    dict.add("nCharges", nCharges_);
     os  << indent << dict.dictName() << dict;
 }
 
@@ -93,8 +72,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const specie& st)
 {
     os  << st.name_ << tab
         << st.nMoles_ << tab
-        << st.molWeight_ << tab
-        << st.nCharges_;
+        << st.molWeight_;
 
     os.check("Ostream& operator<<(Ostream& os, const specie& st)");
     return os;

src/thermophysicalModels/specie/specie/specie.H

@@ -77,9 +77,6 @@ class specie
         //- Molecular weight of specie [kg/kmol]
         scalar molWeight_;
 
-        //- Number of elementary charges of specie
-        scalar nCharges_;
-
 
 public:
 
@@ -87,29 +84,6 @@ public:
     ClassName("specie");
 
 
-    // Public constants
-
-        // Thermodynamic constants
-
-            //- Universal gas constant [J/(kmol K)]
-            static const scalar RR;
-
-            //- Standard pressure [Pa]
-            static const scalar Pstd;
-
-            //- Standard temperature [K]
-            static const scalar Tstd;
-
-            //- Elementary charge [C]
-            static const scalar e;
-
-            //- Avogadro number [1/kmol]
-            static const scalar NA;
-
-            //- Boltzmann constant [J/K]
-            static const scalar k;
-
-
     // Constructors
 
 
@@ -124,23 +98,6 @@ public:
             const scalar molWeight
         );
 
-        //- Construct from components without name
-        inline specie
-        (
-            const scalar nMoles,
-            const scalar molWeight,
-            const scalar nCharges
-        );
-
-        //- Construct from components with name
-        inline specie
-        (
-            const word& name,
-            const scalar nMoles,
-            const scalar molWeight,
-            const scalar nCharges
-        );
-
         //- Construct as copy
         inline specie(const specie&);
 
@@ -170,27 +127,6 @@ public:
             //- Gas constant [J/(kg K)]
             inline scalar R() const;
 
-            //- Charge number
-            inline scalar z() const;
-
-
-        // Electric charge function
-
-
-            // Mole specific properties
-
-                //- Electric charge [C/kmol]
-                inline scalar qc() const;
-
-
-            // Mass specific properties
-
-                //- Electric charge [C/kg]
-                inline scalar Qc() const;
-
-                //- Electric charge squared [C^2/kg]
-                inline scalar Qc2() const;
-
 
         // I-O
 

src/thermophysicalModels/specie/specie/specieI.H

@@ -41,8 +41,7 @@ inline specie::specie
 :
     name_(name),
     nMoles_(nMoles),
-    molWeight_(molWeight),
-    nCharges_(0.0)
+    molWeight_(molWeight)
 {}
 
 
@@ -53,36 +52,7 @@ inline specie::specie
 )
 :
     nMoles_(nMoles),
-    molWeight_(molWeight),
-    nCharges_(0.0)
-{}
-
-
-inline specie::specie
-(
-    const word& name,
-    const scalar nMoles,
-    const scalar molWeight,
-    const scalar nCharges
-)
-:
-    name_(name),
-    nMoles_(nMoles),
-    molWeight_(molWeight),
-    nCharges_(nCharges)
-{}
-
-
-inline specie::specie
-(
-    const scalar nMoles,
-    const scalar molWeight,
-    const scalar nCharges
-)
-:
-    nMoles_(nMoles),
-    molWeight_(molWeight),
-    nCharges_(nCharges)
+    molWeight_(molWeight)
 {}
 
 
@@ -92,8 +62,7 @@ inline specie::specie(const specie& st)
 :
     name_(st.name_),
     nMoles_(st.nMoles_),
-    molWeight_(st.molWeight_),
-    nCharges_(st.nCharges_)
+    molWeight_(st.molWeight_)
 {}
 
 
@@ -101,8 +70,7 @@ inline specie::specie(const word& name, const specie& st)
 :
     name_(name),
     nMoles_(st.nMoles_),
-    molWeight_(st.molWeight_),
-    nCharges_(st.nCharges_)
+    molWeight_(st.molWeight_)
 {}
 
 
@@ -132,30 +100,6 @@ inline scalar specie::R() const
 }
 
 
-inline scalar specie::z() const
-{
-    return nCharges_;
-}
-
-
-inline scalar specie::qc() const
-{
-    return nCharges_ * NA * e;
-}
-
-
-inline scalar specie::Qc() const
-{
-    return nCharges_ * NA * e / molWeight_;
-}
-
-
-inline scalar specie::Qc2() const
-{
-    return nCharges_ * nCharges_ * NA * e * e / molWeight_;
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline void specie::operator=(const specie& st)
@@ -163,7 +107,6 @@ inline void specie::operator=(const specie& st)
     //name_ = st.name_;
     nMoles_ = st.nMoles_;
     molWeight_ = st.molWeight_;
-    nCharges_ = st.nCharges_;
 }
 
 
@@ -171,10 +114,6 @@ inline void specie::operator+=(const specie& st)
 {
     scalar sumNmoles = max(nMoles_ + st.nMoles_, SMALL);
 
-    nCharges_ =
-        nMoles_/sumNmoles*nCharges_
-      + st.nMoles_/sumNmoles*st.nCharges_;
-
     molWeight_ =
         nMoles_/sumNmoles*molWeight_
       + st.nMoles_/sumNmoles*st.molWeight_;
@@ -191,10 +130,6 @@ inline void specie::operator-=(const specie& st)
         diffnMoles = SMALL;
     }
 
-    nCharges_ =
-        nMoles_/diffnMoles*nCharges_
-      - st.nMoles_/diffnMoles*st.nCharges_;
-
     molWeight_ =
         nMoles_/diffnMoles*molWeight_
       - st.nMoles_/diffnMoles*st.molWeight_;
@@ -219,9 +154,7 @@ inline specie operator+(const specie& st1, const specie& st2)
     (
         sumNmoles,
         st1.nMoles_/sumNmoles*st1.molWeight_
-      + st2.nMoles_/sumNmoles*st2.molWeight_,
-        st1.nMoles_/sumNmoles*st1.nCharges_
-      + st2.nMoles_/sumNmoles*st2.nCharges_
+      + st2.nMoles_/sumNmoles*st2.molWeight_
     );
 }
 
@@ -238,9 +171,7 @@ inline specie operator-(const specie& st1, const specie& st2)
     (
         diffNmoles,
         st1.nMoles_/diffNmoles*st1.molWeight_
-      - st2.nMoles_/diffNmoles*st2.molWeight_,
-        st1.nMoles_/diffNmoles*st1.nCharges_
-      - st2.nMoles_/diffNmoles*st2.nCharges_
+      - st2.nMoles_/diffNmoles*st2.molWeight_
     );
 }
 
@@ -250,8 +181,7 @@ inline specie operator*(const scalar s, const specie& st)
     return specie
     (
         s*st.nMoles_,
-        st.molWeight_,
-        st.nCharges_
+        st.molWeight_
     );
 }
 

src/thermophysicalModels/specie/thermo/thermo/thermo.H

@@ -336,7 +336,6 @@ public:
             ) const;
 
 
-
         // I-O
 
             //- Write to Ostream