ENH: Adding diffusionMulticomponent combustion model.

Adding optional files to smallPoolFire2D to run using this model.
Taking out of the compilation of FSD combustion. It needs futher work to run using the new turbulent framework

src/combustionModels/FSD/FSD.C

@@ -56,7 +56,7 @@ FSD<CombThermoType, ThermoType>::FSD
         (
             this->coeffs(),
             this->mesh(),
-            *this
+           *this
         )
     ),
     ft_
@@ -188,18 +188,17 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
     volScalarField& omegaFuelBar = tomegaFuel();
 
     // Calculation of the mixture fraction variance (ftVar)
-    const compressible::LESModel& lesModel =
+    // TODO: generalize delta for RAS and LES.
-        YO2.db().lookupObject<compressible::LESModel>("LESProperties");
+    const volScalarField& delta =
+        refCast<const compressible::LESModel>(this->turbulence()).delta();
 
-    const volScalarField& delta = lesModel.delta();
     const volScalarField ftVar(Cv_*sqr(delta)*sqr(mgft));
 
     // Thickened flame (average flame thickness for counterflow configuration
     // is 1.5 mm)
-
     volScalarField  deltaF
     (
-        lesModel.delta()/dimensionedScalar("flame", dimLength, 1.5e-3)
+        delta/dimensionedScalar("flame", dimLength, 1.5e-3)
     );
 
     // Linear correlation between delta and flame thickness

src/combustionModels/FSD/reactionRateFlameAreaModels/relaxation/relaxation.C

@@ -97,13 +97,23 @@ void Foam::reactionRateFlameAreaModels::relaxation::correct
         1e-4
     );
 
-    const compressible::LESModel& lesModel =
+//     const compressible::LESModel& lesModel =
-        omega_.db().lookupObject<compressible::LESModel>("LESProperties");
+//         omega_.db().lookupObject<compressible::LESModel>("LESProperties");
 
     // Total strain : resolved and sub-grid (just LES for now)
+//     const volScalarField sigmaTotal
+//     (
+//         sigma + alpha_*lesModel.epsilon()/(lesModel.k() + lesModel.kMin())
+//     );
+
     const volScalarField sigmaTotal
     (
-        sigma + alpha_*lesModel.epsilon()/(lesModel.k() + lesModel.kMin())
+        sigma
+      + alpha_*combModel_.turbulence().epsilon()
+      / (
+            combModel_.turbulence().k()
+          + dimensionedScalar("kMin", sqr(dimVelocity), SMALL)
+        )
     );
 
     const volScalarField omegaInf(correlation_.omega0Sigma(sigmaTotal));
@@ -118,20 +128,20 @@ void Foam::reactionRateFlameAreaModels::relaxation::correct
        /(sqr(omega0 - omegaInf) + sqr(omegaMin))
     );
 
-    const volScalarField rho(combModel_.rho());
+    tmp<surfaceScalarField> phi(combModel_.phi());
-    const surfaceScalarField phi(combModel_.phi());
 
     solve
     (
-         fvm::ddt(rho, omega_)
+         fvm::ddt(omega_)
        + fvm::div(phi, omega_, "div(phi,omega)")
       ==
-         rho*Rc*omega0
+         Rc*omega0
-       - fvm::SuSp(rho*(tau + Rc), omega_)
+       - fvm::SuSp((tau + Rc), omega_)
     );
 
     omega_.min(omega0);
     omega_.max(0.0);
+
 }
 
 

src/combustionModels/Make/files

@@ -17,12 +17,15 @@ PaSR/PaSRs.C
 
 laminar/laminars.C
 
+/*
 FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
 FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameArea.C
 FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameAreaNew.C
 FSD/reactionRateFlameAreaModels/relaxation/relaxation.C
-
 FSD/FSDs.C
+*/
+
+diffusionMulticomponent/diffusionMulticomponents.C
 
 noCombustion/noCombustions.C
 

src/combustionModels/combustionModel/combustionModel.H

@@ -112,9 +112,9 @@ public:
             inline const fvMesh& mesh() const;
 
             //- Return const access to phi
-            inline const surfaceScalarField& phi() const;
+            inline tmp<surfaceScalarField> phi() const;
 
-            //- Return const access to rho
+            //- Returns rho
             virtual tmp<volScalarField> rho() const = 0;
 
             //- Return access to turbulence

src/combustionModels/combustionModel/combustionModelI.H

@@ -31,7 +31,7 @@ inline const Foam::fvMesh& Foam::combustionModel::mesh() const
 }
 
 
-inline const Foam::surfaceScalarField& Foam::combustionModel::phi() const
+inline Foam::tmp<Foam::surfaceScalarField> Foam::combustionModel::phi() const
 {
     if (turbulencePtr_)
     {
@@ -41,7 +41,7 @@ inline const Foam::surfaceScalarField& Foam::combustionModel::phi() const
     {
         FatalErrorIn
         (
-            "const Foam::compressibleTurbulenceModel& "
+            "const Foam::tmp<Foam::surfaceScalarField> "
             "Foam::combustionModel::turbulence() const "
         )   << "turbulencePtr_ is empty. Please use "
             << "combustionModel::setTurbulence "

src/combustionModels/diffusionMulticomponent/diffusionMulticomponent.C

@@ -0,0 +1,515 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd
+     \\/     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 "diffusionMulticomponent.H"
+#include "fvcGrad.H"
+#include "reactingMixture.H"
+#include "fvCFD.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+void Foam::combustionModels::
+diffusionMulticomponent<CombThermoType, ThermoType>::init()
+{
+    // Load default values
+    if (this->coeffs().found("Ci"))
+    {
+        this->coeffs().lookup("Ci") >> Ci_;
+    }
+
+    if (this->coeffs().found("YoxStream"))
+    {
+        this->coeffs().lookup("YoxStream") >> YoxStream_;
+    }
+
+    if (this->coeffs().found("YfStream"))
+    {
+        this->coeffs().lookup("YfStream") >> YfStream_;
+    }
+
+    if (this->coeffs().found("sigma"))
+    {
+        this->coeffs().lookup("sigma") >> sigma_;
+    }
+
+    if (this->coeffs().found("ftCorr"))
+    {
+        this->coeffs().lookup("ftCorr") >> ftCorr_;
+    }
+
+    if (this->coeffs().found("alpha"))
+    {
+        alpha_ = readScalar(this->coeffs().lookup("alpha"));
+    }
+
+    if (this->coeffs().found("laminarIgn"))
+    {
+        this->coeffs().lookup("laminarIgn") >> laminarIgn_;
+    }
+
+    typedef typename Reaction<ThermoType>::specieCoeffs specieCoeffs;
+
+    const speciesTable& species = this->thermo().composition().species();
+
+    scalarList specieStoichCoeffs(species.size());
+    const label nReactions = reactions_.size();
+
+    for (label k=0; k < nReactions; k++)
+    {
+        RijPtr_.set
+        (
+            k,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "Rijk" + name(k),
+                    this->mesh_.time().timeName(),
+                    this->mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                this->mesh_,
+                dimensionedScalar("Rij", dimMass/dimTime/dimVolume, 0.0),
+                zeroGradientFvPatchScalarField::typeName
+            )
+        );
+
+        RijPtr_[k].storePrevIter();
+
+        const List<specieCoeffs>& lhs = reactions_[k].lhs();
+        const List<specieCoeffs>& rhs = reactions_[k].rhs();
+
+        const label fuelIndex = species[fuelNames_[k]];
+        const label oxydantIndex = species[oxydantNames_[k]];
+
+        const scalar Wu = specieThermo_[fuelIndex].W();
+        const scalar Wox = specieThermo_[oxydantIndex].W();
+
+        forAll(lhs, i)
+        {
+            const label specieI = lhs[i].index;
+            specieStoichCoeffs[specieI] = -lhs[i].stoichCoeff;
+            qFuel_[k] +=
+                specieThermo_[specieI].hc()*lhs[i].stoichCoeff/Wu;
+        }
+
+        forAll(rhs, i)
+        {
+            const label specieI = rhs[i].index;
+            specieStoichCoeffs[specieI] = rhs[i].stoichCoeff;
+            qFuel_[k] -=
+                specieThermo_[specieI].hc()*rhs[i].stoichCoeff/Wu;
+        }
+
+        Info << "Fuel heat of combustion : " << qFuel_[k] << endl;
+
+        s_[k] =
+            (Wox*mag(specieStoichCoeffs[oxydantIndex]))
+          / (Wu*mag(specieStoichCoeffs[fuelIndex]));
+
+        Info << "stoichiometric oxygen-fuel ratio : " << s_[k] << endl;
+
+        stoicRatio_[k] = s_[k]*YfStream_[k]/YoxStream_[k];
+
+        Info << "stoichiometric air-fuel ratio : " << stoicRatio_[k] << endl;
+
+        const scalar fStoich = 1.0/(1.0 + stoicRatio_[k]);
+
+        Info << "stoichiometric mixture fraction : " << fStoich << endl;
+
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+Foam::combustionModels::diffusionMulticomponent<CombThermoType, ThermoType>::
+diffusionMulticomponent
+(
+    const word& modelType, const fvMesh& mesh, const word& phaseName
+)
+:
+    CombThermoType(modelType, mesh, phaseName),
+    reactions_
+    (
+        dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo())
+    ),
+    specieThermo_
+    (
+        dynamic_cast<const reactingMixture<ThermoType>&>
+            (this->thermo()).speciesData()
+    ),
+    RijPtr_(reactions_.size()),
+    Ci_(reactions_.size(), 1.0),
+    fuelNames_(this->coeffs().lookup("fuels")),
+    oxydantNames_(this->coeffs().lookup("oxydants")),
+    qFuel_(reactions_.size()),
+    stoicRatio_(reactions_.size()),
+    s_(reactions_.size()),
+    YoxStream_(reactions_.size(), 0.23),
+    YfStream_(reactions_.size(), 1.0),
+    sigma_(reactions_.size(), 0.02),
+    oxydantRes_(this->coeffs().lookup("oxydantRes")),
+    ftCorr_(reactions_.size(), 0.0),
+    alpha_(1),
+    laminarIgn_(false)
+{
+    init();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+Foam::combustionModels::diffusionMulticomponent<CombThermoType, ThermoType>::
+~diffusionMulticomponent()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField> Foam::combustionModels::
+diffusionMulticomponent<CombThermoType, ThermoType>::tc() const
+{
+    return this->chemistryPtr_->tc();
+}
+
+
+template<class CombThermoType, class ThermoType>
+void Foam::combustionModels::
+diffusionMulticomponent<CombThermoType, ThermoType>::correct()
+{
+    if (this->active())
+    {
+        typedef typename Reaction<ThermoType>::specieCoeffs specieCoeffs;
+        const speciesTable& species = this->thermo().composition().species();
+
+        const label nReactions = reactions_.size();
+
+        PtrList<volScalarField> RijlPtr(nReactions);
+
+        for (label k=0; k < nReactions; k++)
+        {
+
+            RijlPtr.set
+            (
+                k,
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "Rijl" + word(k),
+                        this->mesh_.time().timeName(),
+                        this->mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE,
+                        false
+                    ),
+                    this->mesh_,
+                    dimensionedScalar("Rijl", dimMass/dimTime/dimVolume, 0.0),
+                    zeroGradientFvPatchScalarField::typeName
+                )
+            );
+
+            volScalarField& Rijl = RijlPtr[k];
+
+            // Obtain individual reaction rates for each reaction
+            const label fuelIndex = species[fuelNames_[k]];
+
+            if (laminarIgn_)
+            {
+                Rijl.dimensionedInternalField() =
+                    -this->chemistryPtr_->calculateRR(k, fuelIndex);
+            }
+
+
+             // Look for the fuelStoic
+            const List<specieCoeffs>& rhs = reactions_[k].rhs();
+            const List<specieCoeffs>& lhs = reactions_[k].lhs();
+
+            // Set to zero RR's
+            forAll (lhs, l)
+            {
+                const label lIndex = lhs[l].index;
+                this->chemistryPtr_->RR(lIndex) =
+                    dimensionedScalar("zero", dimMass/dimTime/dimVolume, 0.0);
+            }
+
+            forAll (rhs, l)
+            {
+                const label rIndex = rhs[l].index;
+                this->chemistryPtr_->RR(rIndex) =
+                    dimensionedScalar("zero", dimMass/dimTime/dimVolume, 0.0);
+            }
+        }
+
+
+        for (label k=0; k < nReactions; k++)
+        {
+            const label fuelIndex = species[fuelNames_[k]];
+            const label oxydantIndex = species[oxydantNames_[k]];
+
+            const volScalarField& Yfuel =
+                this->thermo().composition().Y(fuelIndex);
+
+            const volScalarField& Yox =
+                this->thermo().composition().Y(oxydantIndex);
+
+            const volScalarField ft
+            (
+                "ft" + name(k),
+                (
+                    s_[k]*Yfuel
+                  - (Yox - YoxStream_[k])
+                )
+                /
+                (
+                    s_[k]*YfStream_[k] + YoxStream_[k]
+                )
+            );
+
+            const scalar sigma = sigma_[k];
+
+            const scalar fStoich = 1.0/(1.0 + stoicRatio_[k]) + ftCorr_[k];
+
+            const volScalarField OAvailScaled
+            (
+                "OAvailScaled",
+                Yox/max(oxydantRes_[k], 1e-3)
+            );
+
+            const volScalarField preExp
+            (
+                "preExp" + name(k),
+                 1.0  + sqr(OAvailScaled)
+            );
+
+            const volScalarField filter
+            (
+                (1.0/(sigma*sqrt(2.0*constant::mathematical::pi)))
+              * exp(-sqr(ft - fStoich)/(2*sqr(sigma)))
+            );
+
+            const volScalarField topHatFilter
+            (
+                pos(filter - 1e-3)
+            );
+
+            const volScalarField prob
+            (
+                "prob" + name(k), preExp*filter
+            );
+
+            const volScalarField RijkDiff
+            (
+               "RijkDiff",
+                Ci_[k]*this->turbulence().muEff()*prob*
+                (
+                    mag(fvc::grad(Yfuel) & fvc::grad(Yox))
+                )
+               *pos(Yox)*pos(Yfuel)
+            );
+
+            volScalarField& Rijk = RijPtr_[k];
+
+            if (laminarIgn_)
+            {
+                Rijk =
+                    min(RijkDiff, topHatFilter*RijlPtr[k]*pos(Yox)*pos(Yfuel));
+            }
+            else
+            {
+                Rijk = RijkDiff;
+            }
+
+            Rijk.relax(alpha_);
+
+            if (this->mesh_.time().outputTime() && debug)
+            {
+                Rijk.write();
+                ft.write();
+            }
+
+            // Look for the fuelStoic
+            const List<specieCoeffs>& rhs = reactions_[k].rhs();
+            const List<specieCoeffs>& lhs = reactions_[k].lhs();
+
+            scalar fuelStoic = 1.0;
+            forAll (lhs, l)
+            {
+                const label lIndex = lhs[l].index;
+                if (lIndex == fuelIndex)
+                {
+                    fuelStoic = lhs[l].stoichCoeff;
+                    break;
+                }
+            }
+
+            const scalar MwFuel = specieThermo_[fuelIndex].W();
+
+            // Update left hand side species
+            forAll (lhs, l)
+            {
+                const label lIndex = lhs[l].index;
+
+                const scalar stoichCoeff = lhs[l].stoichCoeff;
+
+                this->chemistryPtr_->RR(lIndex) +=
+                    -Rijk*stoichCoeff*specieThermo_[lIndex].W()
+                   /fuelStoic/MwFuel;
+
+            }
+
+            // Update right hand side species
+            forAll (rhs, r)
+            {
+                const label rIndex = rhs[r].index;
+
+                const scalar stoichCoeff = rhs[r].stoichCoeff;
+
+                this->chemistryPtr_->RR(rIndex) +=
+                    Rijk*stoichCoeff*specieThermo_[rIndex].W()
+                   /fuelStoic/MwFuel;
+            }
+        }
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::combustionModels::diffusionMulticomponent<CombThermoType, ThermoType>::
+R(volScalarField& Y) const
+{
+    tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
+
+    fvScalarMatrix& Su = tSu();
+
+    if (this->active())
+    {
+        const label specieI = this->thermo().composition().species()[Y.name()];
+        Su += this->chemistryPtr_->RR(specieI);
+    }
+
+    return tSu;
+}
+
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::diffusionMulticomponent<CombThermoType, ThermoType>::
+dQ() const
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "dQ",
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    if (this->active())
+    {
+        volScalarField& dQ = tdQ();
+        dQ = this->chemistryPtr_->dQ();
+    }
+
+    return tdQ;
+}
+
+
+template<class CombThermoType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::diffusionMulticomponent<CombThermoType, ThermoType>::
+Sh() const
+{
+    tmp<volScalarField> tSh
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Sh",
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    if (this->active())
+    {
+        scalarField& Sh = tSh();
+        Sh = this->chemistryPtr_->Sh();
+    }
+
+    return tSh;
+}
+
+
+template<class CombThermoType, class ThermoType>
+bool Foam::combustionModels::
+diffusionMulticomponent<CombThermoType, ThermoType>::read()
+{
+    if (CombThermoType::read())
+    {
+        this->coeffs().readIfPresent("Ci", Ci_);
+        this->coeffs().readIfPresent("sigma", sigma_);
+        this->coeffs().readIfPresent("oxydantRes", oxydantRes_);
+        this->coeffs().readIfPresent("ftCorr", ftCorr_);
+        this->coeffs().readIfPresent("alpha", alpha_);
+        this->coeffs().readIfPresent("laminarIgn", laminarIgn_);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/combustionModels/diffusionMulticomponent/diffusionMulticomponent.H

@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::diffusionMulticomponent
+
+Description
+
+    Diffusion based turbulent combustion model for multicomponent species.
+
+    The model calculates the laminar finite rate source terms based on
+    the kinetic for each reaction in order to begin the combustion and
+    evaluates the minimum between this and the cross diffusion rate term
+    defined as D*prob*muEff*mag(grad(Yi)*grad(Yj)) if laminarIgn is true.
+
+    where:
+
+    D :     is a model constant.
+    muEff : is the effective turbulent diffusivity
+    prob :  is a Gaussian shaped distribution around the stoichiometric value
+            of each reaction. The distribtion has the input parameter 'sigma'
+            for standard deviation.
+
+    The variable prob is multiplied by the factor: 1  + pow2(O2/oxydantRes),
+    where oxydantRes is the residual oxydant specified for each reaction.
+
+    In the combustion properties dictionary:
+
+    diffusionMulticomponentCoeffs
+    {
+        Ci          (1.0 1.0);      // Default to 1
+        fuels       (CH4 CO);
+        oxydants    (O2 O2);
+        YoxStream   (0.23 0.23);    // Default to 0.23
+        YfStream    (1.0 1.0);      // Default to 1.0
+        sigma       (0.02 0.02);    // Default to 0.02
+        oxydantRes  (0.025 0.005);
+        ftCorr      (0.0 0.0);      // Default to 0.0
+        laminarIgn  false;          // Default false
+    }
+
+    ftCorr is used to shift the location of the flame and corrects the
+    stochimetric mixture value when the mesh resolution is not enough
+    to resolve the combustion zone.
+
+    NOTE: Optionally the switch 'laminarIgn' can be used to ignite the mixture
+    using laminar combustion.
+
+
+SourceFiles
+    diffusionMulticomponent.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef diffusionMulticomponent_H
+#define diffusionMulticomponent_H
+
+#include "scalarList.H"
+#include "tmp.H"
+#include "Reaction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class diffusionMulticomponent Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType, class ThermoType>
+class diffusionMulticomponent
+:
+    public CombThermoType
+{
+    // Private data
+
+        //- Reactions
+        const PtrList<Reaction<ThermoType> >& reactions_;
+
+        //- Thermodynamic data of the species
+        const PtrList<ThermoType>& specieThermo_;
+
+        //- Pointer list of source terms
+        PtrList<volScalarField> RijPtr_;
+
+        //- Model constants
+        scalarList Ci_;
+
+        //- List of fuels for each reaction
+        wordList fuelNames_;
+
+        //- List of oxydants for each reaction
+        wordList oxydantNames_;
+
+        //- Heat of combustion [J/Kg]
+        scalarList qFuel_;
+
+        //- Stoichiometric air-fuel mass ratio
+        scalarList stoicRatio_;
+
+        //- Stoichiometric oxygen-fuel mass ratio
+        scalarList s_;
+
+        //- Oxydaser sream mass concentrations
+        scalarList YoxStream_;
+
+        //- Fuel stream mass concentrations
+        scalarList YfStream_;
+
+        //- Mean distribution for gaussian probabililty
+        scalarList sigma_;
+
+        //- Residual oxydaser
+        scalarList oxydantRes_;
+
+        //- ft stochiometric correction
+        scalarList ftCorr_;
+
+        //- Relaxatnio factor on total source
+        scalar alpha_;
+
+        //- Switch on to laminar combustion for ignition
+        bool laminarIgn_;
+
+
+    // Private Member Functions
+
+        //- Return the chemical time scale
+        tmp<volScalarField> tc() const;
+
+        //- Initialize
+        void init();
+
+        //- Disallow copy construct
+        diffusionMulticomponent(const diffusionMulticomponent&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const diffusionMulticomponent&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("diffusionMulticomponent");
+
+
+    // Constructors
+
+        //- Construct from components
+        diffusionMulticomponent
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~diffusionMulticomponent();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+            //- Fuel consumption rate matrix.
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
+
+            //- Heat release rate calculated from fuel consumption rate matrix
+            virtual tmp<volScalarField> dQ() const;
+
+            //-  Return source for enthalpy equation [kg/m/s3]
+            virtual tmp<volScalarField> Sh() const;
+
+
+    // IO
+
+            //- Update properties from given dictionary
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "diffusionMulticomponent.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/diffusionMulticomponent/diffusionMulticomponents.C

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

tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/chemistryProperties

@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      chemistryProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+chemistryType
+{
+    chemistrySolver   noChemistrySolver;
+    chemistryThermo   psi;
+}
+
+chemistry           on;
+
+initialChemicalTimeStep 1e-07;
+
+
+// ************************************************************************* //

tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties

@@ -16,7 +16,12 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 combustionModel  infinitelyFastChemistry<psiThermoCombustion,gasHThermoPhysics>;
+
 //combustionModel  FSD<psiThermoCombustion,gasHThermoPhysics>;
+//combustionModel  diffusionMulticomponent<psiChemistryCombustion,gasHThermoPhysics>;
+//NOTE: To use diffusionMulticomponent combustion model you need to rename files:
+// reactions.twoSteps -> reactions
+// thermophysicalProperties.twoSteps -> thermophysicalProperties
 
 active  true;
 
@@ -26,8 +31,22 @@ infinitelyFastChemistryCoeffs
     C           5.0;
 }
 
+diffusionMulticomponentCoeffs
+{
+    Ci          (1.0 1.5);
+    fuels       (CH4 CO);
+    oxydants    (O2 O2);
+    YoxStream   (0.23 0.23);
+    YfStream    (1.0 1.0);
+    sigma       (0.02 0.02);
+    oxydantRes  (0.015 0.005);
+    ftCorr      (0  0);
+    laminarIgn  false;
+}
+
 FSDCoeffs
 {
+    semiImplicit no;
     Cv          0.1;
     ftVarMin    1e-2;
 

tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/polyMesh/boundary

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
+|  \\    /   O peration     | Version:  dev-OpenCFD                           |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/

tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties

@@ -189,7 +189,7 @@ greyMeanAbsorptionEmissionCoeffs
 
 scatterModel    none;
 
-sootModel mixtureFractionSoot<gasHThermoPhysics>;
+sootModel none;//mixtureFractionSoot<gasHThermoPhysics>;
 
 mixtureFractionSootCoeffs
 {

tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactions

@@ -9,7 +9,7 @@ species
 
 reactions
 {
-    propaneReaction
+    methaneReaction
     {
         type         irreversibleinfiniteReaction;
         reaction     "CH4 + 2O2 + 7.5N2 = CO2 + 2H2O + 7.5N2";

tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactions.twoSteps

@@ -0,0 +1,30 @@
+species
+(
+    O2
+    H2O
+    CH4
+    CO2
+    CO
+    N2
+);
+
+reactions
+{
+    methaneReaction
+    {
+        type        irreversibleArrheniusReaction;
+        reaction    "CH4^0.9 + 2O2^1.1 = CO + 2H2O";
+        A           2.119e12;
+        beta        0;
+        Ta          17676;
+    }
+
+    COReaction
+    {
+        type         irreversibleArrheniusReaction;
+        reaction     "CO + 0.5O2^0.25 = CO2";
+        A            1e6;
+        beta         0;
+        Ta           6060;
+    }
+}

tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/thermo.compressibleGas

@@ -124,3 +124,26 @@ N2
         Ts              170.672;
     }
 }
+
+CO
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       28.0106;
+    }
+    thermodynamics
+    {
+        Tlow            200;
+        Thigh           6000;
+        Tcommon         1000;
+        highCpCoeffs    ( 3.04849 0.00135173 -4.85794e-07 7.88536e-11 -4.69807e-15 -14266.1 6.0171 );
+        lowCpCoeffs     ( 3.57953 -0.000610354 1.01681e-06 9.07006e-10 -9.04424e-13 -14344.1 3.50841 );
+    }
+    transport
+    {
+        As              1.67212e-06;
+        Ts              170.672;
+    }
+}
+

tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/thermophysicalProperties.twoSteps

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.2.2                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            hePsiThermo;
+    mixture         reactingMixture;
+    transport       sutherland;
+    thermo          janaf;
+    energy          sensibleEnthalpy;
+    equationOfState perfectGas;
+    specie          specie;
+}
+
+inertSpecie N2;
+
+fuel        CH4;
+
+chemistryReader foamChemistryReader;
+
+foamChemistryFile "$FOAM_CASE/constant/reactions";
+
+foamChemistryThermoFile "$FOAM_CASE/constant/thermo.compressibleGas";
+
+
+// ************************************************************************* //

tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSchemes

@@ -37,10 +37,11 @@ divSchemes
         N2              limitedLinear01 1;
         H2O             limitedLinear01 1;
         CO2             limitedLinear01 1;
+        CO              limitedLinear01 1;
         h               limitedLinear 1;
     };
     div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
-    div(phi,omega)  Gauss limitedLinear 1;
+    div(phi,omega)  Gauss upwind;
     div(phi,K)      Gauss limitedLinear 1;
     div(U)          Gauss linear;
     div(Ji,Ii_h)    Gauss upwind;