Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-dev-OpenCFD into develop

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;