From 4ba032b2be067d1b73d7bee71b6d86baf873ad61 Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Mon, 7 Dec 2015 17:02:18 -0800
Subject: [PATCH] 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                |  11 +-
 .../relaxation/relaxation.C                   |  26 +-
 src/combustionModels/Make/files               |   5 +-
 .../combustionModel/combustionModel.H         |   4 +-
 .../combustionModel/combustionModelI.H        |   4 +-
 .../diffusionMulticomponent.C                 | 515 ++++++++++++++++++
 .../diffusionMulticomponent.H                 | 224 ++++++++
 .../diffusionMulticomponents.C                |  61 +++
 .../constant/chemistryProperties              |  29 +
 .../constant/combustionProperties             |  19 +
 .../constant/polyMesh/boundary                |   2 +-
 .../constant/radiationProperties              |   2 +-
 .../les/smallPoolFire2D/constant/reactions    |   2 +-
 .../constant/reactions.twoSteps               |  30 +
 .../constant/thermo.compressibleGas           |  23 +
 .../thermophysicalProperties.twoSteps         |  40 ++
 .../les/smallPoolFire2D/system/fvSchemes      |   3 +-
 17 files changed, 977 insertions(+), 23 deletions(-)
 create mode 100644 src/combustionModels/diffusionMulticomponent/diffusionMulticomponent.C
 create mode 100644 src/combustionModels/diffusionMulticomponent/diffusionMulticomponent.H
 create mode 100644 src/combustionModels/diffusionMulticomponent/diffusionMulticomponents.C
 create mode 100644 tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/chemistryProperties
 create mode 100644 tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactions.twoSteps
 create mode 100644 tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/thermophysicalProperties.twoSteps

diff --git a/src/combustionModels/FSD/FSD.C b/src/combustionModels/FSD/FSD.C
index 9d79904cd40..ed3cdb602c4 100644
--- a/src/combustionModels/FSD/FSD.C
+++ b/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 =
-        YO2.db().lookupObject<compressible::LESModel>("LESProperties");
+    // TODO: generalize delta for RAS and LES.
+    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
diff --git a/src/combustionModels/FSD/reactionRateFlameAreaModels/relaxation/relaxation.C b/src/combustionModels/FSD/reactionRateFlameAreaModels/relaxation/relaxation.C
index 9352dd83961..4294331f258 100644
--- a/src/combustionModels/FSD/reactionRateFlameAreaModels/relaxation/relaxation.C
+++ b/src/combustionModels/FSD/reactionRateFlameAreaModels/relaxation/relaxation.C
@@ -97,13 +97,23 @@ void Foam::reactionRateFlameAreaModels::relaxation::correct
         1e-4
     );
 
-    const compressible::LESModel& lesModel =
-        omega_.db().lookupObject<compressible::LESModel>("LESProperties");
+//     const compressible::LESModel& lesModel =
+//         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());
-    const surfaceScalarField phi(combModel_.phi());
+    tmp<surfaceScalarField> phi(combModel_.phi());
 
     solve
     (
-         fvm::ddt(rho, omega_)
+         fvm::ddt(omega_)
        + fvm::div(phi, omega_, "div(phi,omega)")
       ==
-         rho*Rc*omega0
-       - fvm::SuSp(rho*(tau + Rc), omega_)
+         Rc*omega0
+       - fvm::SuSp((tau + Rc), omega_)
     );
 
     omega_.min(omega0);
     omega_.max(0.0);
+
 }
 
 
diff --git a/src/combustionModels/Make/files b/src/combustionModels/Make/files
index b68fbb742f2..f1f87e9cfbb 100644
--- a/src/combustionModels/Make/files
+++ b/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
 
diff --git a/src/combustionModels/combustionModel/combustionModel.H b/src/combustionModels/combustionModel/combustionModel.H
index 5368d1af801..badd4848294 100644
--- a/src/combustionModels/combustionModel/combustionModel.H
+++ b/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
diff --git a/src/combustionModels/combustionModel/combustionModelI.H b/src/combustionModels/combustionModel/combustionModelI.H
index 65531556828..a69437b73a5 100644
--- a/src/combustionModels/combustionModel/combustionModelI.H
+++ b/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 "
diff --git a/src/combustionModels/diffusionMulticomponent/diffusionMulticomponent.C b/src/combustionModels/diffusionMulticomponent/diffusionMulticomponent.C
new file mode 100644
index 00000000000..bc60afe054a
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/combustionModels/diffusionMulticomponent/diffusionMulticomponent.H b/src/combustionModels/diffusionMulticomponent/diffusionMulticomponent.H
new file mode 100644
index 00000000000..2d6b9531f74
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/combustionModels/diffusionMulticomponent/diffusionMulticomponents.C b/src/combustionModels/diffusionMulticomponent/diffusionMulticomponents.C
new file mode 100644
index 00000000000..3aa6b794019
--- /dev/null
+++ b/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
+);
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/chemistryProperties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/chemistryProperties
new file mode 100644
index 00000000000..60eb1671344
--- /dev/null
+++ b/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;
+
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties
index 51e80f82d64..76a0324cf85 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties
+++ b/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;
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/polyMesh/boundary b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/polyMesh/boundary
index 6a17bae1f0d..0082e381a90 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/polyMesh/boundary
+++ b/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  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
index a25c81382a4..7220f7b9960 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
@@ -189,7 +189,7 @@ greyMeanAbsorptionEmissionCoeffs
 
 scatterModel    none;
 
-sootModel mixtureFractionSoot<gasHThermoPhysics>;
+sootModel none;//mixtureFractionSoot<gasHThermoPhysics>;
 
 mixtureFractionSootCoeffs
 {
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactions b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactions
index 67f2373c382..f8d7139b6b2 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactions
+++ b/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";
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactions.twoSteps b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactions.twoSteps
new file mode 100644
index 00000000000..360014f503f
--- /dev/null
+++ b/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;
+    }
+}
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/thermo.compressibleGas b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/thermo.compressibleGas
index db199425255..5e02e173a2c 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/thermo.compressibleGas
+++ b/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;
+    }
+}
+
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/thermophysicalProperties.twoSteps b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/thermophysicalProperties.twoSteps
new file mode 100644
index 00000000000..ae478175512
--- /dev/null
+++ b/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";
+
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSchemes b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSchemes
index 01608aaaa17..d8cde753f00 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSchemes
+++ b/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;
-- 
GitLab