Merge branch 'EDM' into 'trunk'

EDM,SLFM,ThermoMPPICFoam

See merge request combustion-postech/OpenFOAM-4.x!17

.gitignore

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

applications/solvers/combustion/SLFMFoam/AMC/AMC.C

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "AMC.H"
+
+#include "error.H"
+
+// * * * * * * * * * * * * * Functions * * * * * * * * * * * * //
+
+//Amplitude Mapping Closure (from KIVA)
+Foam::tmp<Foam::scalarField> Foam::CMC::AMC(const scalarField& eta)
+{
+    tmp<scalarField> tRes(new scalarField(eta.size(), 0.0));
+
+    scalarField& Res = tRes.ref();
+
+    for(label i=0 ; i<eta.size() ; i++)
+    {
+        Res[i] = AMC(eta[i]);
+    }
+
+    return tRes;
+}
+
+
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/AMC/AMC.H

@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Function
+    Foam::CMC::AMC
+
+Description
+
+SourceFiles
+    AMC.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef AMC_H
+#define AMC_H
+
+#include "tmp.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace CMC
+{
+
+/*---------------------------------------------------------------------------*\
+                         Function AMC Declaration
+\*---------------------------------------------------------------------------*/
+
+    const scalar pi = 3.141592;
+
+    const scalar spi = Foam::sqrt(pi);
+
+    //Amplitude Mapping Closure (from KIVA)
+    tmp<scalarField> AMC(const scalarField& eta);
+
+
+    //Amplitude Mapping Closure (from KIVA)
+    //Define exp(-2*(erf^-1(2*eta - 1))^2)
+    inline scalar AMC(const scalar eta);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CMC
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "AMCI.H"
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/AMC/AMCI.H

@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+//Amplitude Mapping Closure (from KIVA)
+//Define exp(-2*(erf^-1(2*eta - 1))^2)
+inline Foam::scalar Foam::CMC::AMC(const scalar eta)
+{
+    const scalar a0 = (2.0*eta - 1.0);
+    scalar a1 = 0.5;
+    scalar slope;
+    scalar da = GREAT;
+
+    while(mag(da) > SMALL)
+    {
+        slope = (2.0/spi)*Foam::exp(-1.0*Foam::pow(a1,2.0));
+        da = (a0 - Foam::erf(a1))/slope;
+        a1 = a1 + da;
+    }
+
+    return Foam::exp(-2.0*Foam::pow(a1,2.0));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunction.C

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BetaFunction.H"
+
+#include "error.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// const scalar Foam::CMC::BetaFunction::staticData();
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::CMC::BetaFunction::_setParameters()
+{
+        scalar gamma = mf_*(1.0 - mf_)/(mfVar_ + SMALL) - 1.0;
+
+        if (gamma < SMALL || (mfVar_)/(mf_*(1 - mf_) + SMALL) < 0.001)
+        {
+            fdelta_ = true;
+        }
+        else
+        {
+            a_ = max(mf_*gamma, 0.0);
+            b_ = max((1.0 - mf_)*gamma, 0.0);
+
+            if(a_ < 1.0)
+            {
+                delta0_ = true;
+            }
+            if(b_ < 1.0)
+            {
+                delta1_ = true;
+            }
+
+            _limitAB();
+        }
+}
+
+
+void Foam::CMC::BetaFunction::_limitAB()
+{
+    scalar fmax = 1.0/(1.0 + (b_ - 1.0)/(a_ - 1.0));
+
+    if(a_ > 500.0)
+    {
+        a_ = 500.0;
+        b_ = (a_ - 1.0 - fmax*(a_ - 2.0))/fmax;
+    }
+    else if(b_ > 500.0)
+    {
+        b_ = 500.0;
+        a_ = (1.0 + fmax*(b_ - 2.0))/(1.0 - fmax);
+    }
+}
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CMC::BetaFunction::BetaFunction(const scalar mf, const scalar mfVar)
+:
+    mf_(mf),
+    mfVar_(mfVar),
+    a_(0.0),
+    b_(0.0),
+    fdelta_(false),
+    delta0_(false),
+    delta1_(false)
+{
+    _setParameters();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::CMC::BetaFunction::~BetaFunction()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunction.H

@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::BetaFunction
+
+Description
+
+SourceFiles
+    BetaFunctionI.H
+    BetaFunction.C
+    BetaFunctionIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BetaFunction_H
+#define BetaFunction_H
+
+#include "scalar.H"
+#include "Switch.H"
+#include "scalarField.H"
+#include "tmp.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+namespace CMC
+{
+    class BetaFunction;
+}
+Istream& operator>>(Istream&, CMC::BetaFunction&);
+Ostream& operator<<(Ostream&, const CMC::BetaFunction&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class BetaFunction Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace CMC
+{
+
+class BetaFunction
+{
+    // Private data
+
+        //- Mean mixture fraction
+        scalar mf_;
+
+        //- Mixture fraction variance
+        scalar mfVar_;
+
+        //- Beta pdf parameter alpha
+        scalar a_;
+
+        //- Beta pdf parameter beta
+        scalar b_;
+
+        //- Flag beta pdf special case: forced-delta function
+        Switch fdelta_;
+
+        //- delta function at pure oxidizer stream
+        Switch delta0_;
+
+        //- delta function at pure fuel stream
+        Switch delta1_;
+
+
+    // Private Member Functions
+
+        //- Calculate alpha and beta and set special case flags
+        void _setParameters();
+
+        //- Limit too large value of beta pdf parameters
+        void _limitAB();
+
+        //- Disallow default bitwise copy construct
+        // BetaFunction(const BetaFunction&);
+
+        //- Disallow default bitwise assignment
+        // void operator=(const BetaFunction&);
+
+
+public:
+
+    // Static data members
+
+        //- Static data staticData
+        // static const scalar staticData;
+
+
+    // Constructors
+
+        //- Construct from components
+        BetaFunction(const scalar mf, const scalar mfVar);
+
+        //- Construct from Istream
+        BetaFunction(Istream&);
+
+        //- Construct as copy
+        // BetaFunction(const BetaFunction&);
+
+
+    //- Destructor
+    ~BetaFunction();
+
+
+    // Member Functions
+
+        tmp<scalarField> etaFunc(const UList<scalar>& eta) const
+        {
+            return pow(eta, a_-1.0)*pow(1.0-eta, b_-1.0);
+        }
+
+        scalar etaFunc(const scalar eta) const
+        {
+            return Foam::pow(eta, a_-1.0)*Foam::pow(1.0-eta, b_-1.0);
+        }
+
+
+        // Access
+
+            Switch fdelta() const {return fdelta_;}
+
+            Switch delta0() const {return delta0_;}
+
+            Switch delta1() const {return delta1_;}
+
+            scalar mf() const {return mf_;}
+
+            scalar mfVar() const {return mfVar_;}
+
+            scalar alpha() const {return a_;}
+
+            scalar beta() const {return b_;}
+
+        // Check
+
+        // Edit
+
+        // Write
+
+
+    // Member Operators
+
+        // void operator=(const BetaFunction&);
+
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+        friend Istream& ::Foam::operator>>(Istream&, BetaFunction&);
+        friend Ostream& ::Foam::operator<<(Ostream&, const BetaFunction&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CMC
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "BetaFunctionI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunctionI.H

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

applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunctionIO.C

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

applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGrid.C

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BetaGrid.H"
+
+#include "error.H"
+#include "dictionary.H"
+#include "SubList.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// const scalar Foam::CMC::BetaGrid::staticData();
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CMC::BetaGrid::BetaGrid(const dictionary& dict)
+:
+    etaCut_(dict.lookup("detailedEta")),
+    nSpacings_(dict.lookup("detailedN")),
+    etaSpace_(sum(nSpacings_)+3)
+{
+    // Check for input list size
+    if (etaCut_.size() != nSpacings_.size() + 1)
+    {
+        FatalErrorInFunction
+            << "Number of grid interval points does not match number of grid sizes"
+            << abort(FatalError);
+    }
+
+    // Check for evenness of number of spacings in each interval
+    forAll (nSpacings_, i)
+    {
+        if ((nSpacings_[i] % 2) == 1)
+        {
+            FatalErrorInFunction
+                << "Number of spacings in intervals should be even"
+                << abort(FatalError);
+        }
+    }
+
+    forAll(nSpacings_, i)
+    {
+        const label ns = nSpacings_[i];
+        const labelList::subList prev(nSpacings_, i);
+        const label baseI = sum(prev) + 1;
+
+        const scalar low = etaCut_[i];
+        const scalar upp = etaCut_[i+1];
+        const scalar delta = (upp - low) / scalar(ns);
+
+        for (label j = 0; j < ns; j++)
+        {
+            etaSpace_[baseI+j] = low + delta * j;
+        }
+    }
+    etaSpace_[0] = 0.0;
+    etaSpace_[etaSpace_.size()-2] = etaCut_.last();
+    etaSpace_[etaSpace_.size()-1] = 1.0;
+}
+
+
+Foam::CMC::BetaGrid::BetaGrid(const BetaGrid&)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::CMC::BetaGrid::~BetaGrid()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+void Foam::CMC::BetaGrid::operator=(const BetaGrid& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+}
+
+// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGrid.H

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::BetaGrid
+
+Description
+
+SourceFiles
+    BetaGridI.H
+    BetaGrid.C
+    BetaGridIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BetaGrid_H
+#define BetaGrid_H
+
+#include "scalarField.H"
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+namespace CMC
+{
+    class BetaGrid;
+}
+class dictionary;
+
+Istream& operator>>(Istream&, CMC::BetaGrid&);
+Ostream& operator<<(Ostream&, const CMC::BetaGrid&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class BetaGrid Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace CMC
+{
+
+class BetaGrid
+{
+    // Private data
+
+        //- Piecewise uniform interval points
+        scalarField etaCut_;
+
+        //- Sizes of piecewise uniform intervals
+        labelList nSpacings_;
+
+        //- Mixture fraction eta grid points
+        scalarField etaSpace_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        BetaGrid(const BetaGrid&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const BetaGrid&);
+
+
+public:
+
+    // Static data members
+
+        //- Static data staticData
+        // static const scalar staticData;
+
+
+    // Constructors
+
+        //- Construct from components
+        BetaGrid(const dictionary& dict);
+
+        //- Construct from Istream
+        BetaGrid(Istream&);
+
+        //- Construct as copy
+        // BetaGrid(const BetaGrid&);
+
+
+    //- Destructor
+    ~BetaGrid();
+
+
+    // Member Functions
+
+        // Access
+        const labelList& N() const {return nSpacings_;}
+
+        const scalarField& etaCut() const {return etaCut_;}
+
+        const scalarField& etaSpace() const {return etaSpace_;}
+
+        // Check
+
+        // Edit
+
+        // Write
+
+
+    // Member Operators
+
+        // void operator=(const BetaGrid&);
+
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+        friend Istream& ::Foam::operator>>(Istream&, BetaGrid&);
+        friend Ostream& ::Foam::operator<<(Ostream&, const BetaGrid&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CMC
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "BetaGridI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGridI.H

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

applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGridIO.C

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

applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegrator.C

@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BetaIntegrator.H"
+
+
+#include "error.H"
+#include "SubList.H"
+#include "SubField.H"
+#include "interpolateXY.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// const scalar Foam::CMC::BetaIntegrator::staticData();
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+//extended Simpson's rule (Numerical recipes, 2nd Ed. p.128)
+//for equally spaced and even N intervals (or odd N+1 points)
+Foam::scalar Foam::CMC::BetaIntegrator::simps
+(const scalar xl, const scalar xh, const label n, const UList<scalar>& fx)
+    const
+{
+    scalar evensum(0.0), oddsum(0.0), sum(0.0);
+
+    scalar h = (xh - xl)/scalar(n);
+
+    for(label i=0 ; i<fx.size() ; i++)
+    {
+        if(i%2 == 0)
+        {
+            evensum += fx[i];
+        }
+        else
+        {
+            oddsum += fx[i];
+        }
+    }
+
+    sum = -1.0*fx.first() + 2.0*evensum + 4.0*oddsum - 1.0*fx.last();
+
+    return sum*(h/3.0);
+}
+
+
+Foam::scalar Foam::CMC::BetaIntegrator::sumSimps(const UList<scalar>& f) const
+{
+    scalar total = 0.0;
+
+    const labelList &N(bg_.N());
+    const scalarField &etaCut(bg_.etaCut());
+
+    forAll(N, i)
+    {
+        const labelList::subList prev(N, i);
+        const label baseI = sum(prev) + 1;
+
+        const scalarField::subField subInterval(f, N[i]+1, baseI);
+
+        total += simps(etaCut[i], etaCut[i+1], N[i], subInterval);
+    }
+
+    return total;
+}
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CMC::BetaIntegrator::BetaIntegrator(const BetaFunction& bf, const BetaGrid &bg)
+:
+    etaSpace_(bg.etaSpace()),
+    pdfNum_(bg.etaSpace().size(), 0.0),
+    pdfDen_(0.0),
+    bf_(bf),
+    bg_(bg)
+{
+    const scalarField::subField domain(etaSpace_, etaSpace_.size()-2, 1);
+    scalarField::subField range(pdfNum_, etaSpace_.size()-2, 1);
+
+    if (bf_.fdelta())
+    {
+        pdfDen_ = 1.0;
+    }
+    else
+    {
+        range = bf_.etaFunc(domain);
+
+        if (!bf_.delta0())
+        {
+            pdfNum_.first() = bf_.etaFunc(etaSpace_.first());
+        }
+
+        if (!bf_.delta1())
+        {
+            pdfNum_.last() = bf_.etaFunc(etaSpace_.last());
+        }
+
+        pdfDen_ = sumSimps(pdfNum_)
+            + Foam::pow(etaSpace_[1], bf_.alpha())/(bf_.alpha() + SMALL)
+            + Foam::pow(etaSpace_[1], bf_.beta())/(bf_.beta() + SMALL);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::CMC::BetaIntegrator::~BetaIntegrator()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+// beta-pdf weighted integration for given mf, mfVar and f
+Foam::scalar Foam::CMC::BetaIntegrator::betaIntegrate(const scalarField& f) const
+{
+    scalar result = 0.0;
+
+    if(bf_.fdelta())
+    {
+        result = interpolateXY(bf_.mf(), etaSpace_, f);
+    }
+    else
+    {
+        scalar num = sumSimps (f * pdfNum_)
+            + f.first()*Foam::pow(etaSpace_[1], bf_.alpha())/(bf_.alpha() + SMALL)
+            + f.last()*Foam::pow(etaSpace_[1], bf_.beta())/(bf_.beta() + SMALL);
+        result = num/pdfDen_;
+    }
+
+    return result;
+}
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+void Foam::CMC::BetaIntegrator::operator=(const BetaIntegrator& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+}
+
+// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegrator.H

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::CMC::BetaIntegrator
+
+Description
+
+SourceFiles
+    BetaIntegratorI.H
+    BetaIntegrator.C
+    BetaIntegratorIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BetaIntegrator_H
+#define BetaIntegrator_H
+
+#include "scalarField.H"
+#include "BetaFunction.H"
+#include "BetaGrid.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+namespace CMC
+{
+    class BetaIntegrator;
+}
+Istream& operator>>(Istream&, CMC::BetaIntegrator&);
+Ostream& operator<<(Ostream&, const CMC::BetaIntegrator&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class BetaIntegrator Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace CMC
+{
+
+class BetaIntegrator
+{
+    // Private data
+
+        //- Quadrature points
+        const scalarField &etaSpace_;
+
+        //- Beta pdf numerators
+        scalarField pdfNum_;
+
+        //- Beta pdf denominator
+        scalar pdfDen_;
+
+        //- Description of data_
+        BetaFunction bf_;
+
+        //- Description of data_
+        const BetaGrid &bg_;
+
+
+    // Private Member Functions
+
+        //- Calculate integral using Simpson's rule
+        scalar simps(const scalar xl, const scalar xh, const label n, const UList<scalar>& fx) const;
+
+        //- Calculate integral using Simpson's rule
+        scalar sumSimps(const UList<scalar>& f) const;
+
+        //- Disallow default bitwise copy construct
+        BetaIntegrator(const BetaIntegrator&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const BetaIntegrator&);
+
+
+public:
+
+    // Static data members
+
+        //- Static data staticData
+        // static const scalar staticData;
+
+
+    // Constructors
+
+        //- Construct from components
+        BetaIntegrator(const BetaFunction& bf, const BetaGrid &bg);
+
+        //- Construct as copy
+        // BetaIntegrator(const BetaIntegrator&);
+
+
+    //- Destructor
+    ~BetaIntegrator();
+
+
+    // Member Functions
+
+        scalar betaIntegrate(const scalarField& f) const;
+
+        // Access
+
+        scalar pdfDen () const {return pdfDen_;}
+
+        // Check
+
+        // Edit
+
+        // Write
+
+
+    // Member Operators
+
+        // void operator=(const BetaIntegrator&);
+
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+        friend Istream& ::Foam::operator>>(Istream&, BetaIntegrator&);
+        friend Ostream& ::Foam::operator<<(Ostream&, const BetaIntegrator&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CMC
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "BetaIntegratorI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegratorI.H

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

applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegratorIO.C

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

applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructure.C

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "FlameStructure.H"
+
+#include "fileName.H"
+#include "IFstream.H"
+#include "hashedWordList.H"
+#include "interpolateXY.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::SLFM::FlameStructure::FlameStructure(const string& NstDir, scalarList NstList, const scalarField &etaGrid, const hashedWordList &spTable)
+:
+    NstDir_(NstDir),
+    NstList_(NstList),
+    eta_(etaGrid),
+    species_(spTable),
+    Y_(spTable.size(), scalarFieldArray1d(NstList_.size(), scalarField(eta_.size(), 0.0))),
+    W_(spTable.size(), scalarFieldArray1d(NstList_.size(), scalarField(eta_.size(), 0.0))),
+    T_(NstList_.size(), scalarField(eta_.size(), 0.0)),
+    Q_(NstList_.size(), scalarField(eta_.size(), 0.0)),
+    h_(NstList_.size(), scalarField(eta_.size(), 0.0)),
+    rho_(NstList_.size(), scalarField(eta_.size(), 0.0)),
+    Rgas_(NstList_.size(), scalarField(eta_.size(), 0.0))
+{
+    for (label n = 0; n < NstList_.size(); n++)
+    {
+        scalar N = NstList_[n];
+
+        fileName fname = "sdr"+Foam::name(N)+".inp";
+        IFstream fin(NstDir_/fname);
+
+        string gbg;
+        label etamax_SLFM, NoSpecies;
+
+        //INPUT FILE...
+        fin.getLine(gbg);
+
+        //No.GRID...
+        fin.getLine(gbg);
+
+        //Number of eta-point and species
+        fin >> etamax_SLFM >> NoSpecies;
+
+        //blank line
+        fin.getLine(gbg);
+
+        //PRESSURE...
+        fin.getLine(gbg);
+
+        //1.00...
+        fin.getLine(gbg);
+
+        //MIXTURE...
+        fin.getLine(gbg);
+
+        //read eta space (from SLFM library)
+        scalarField etaValue_SLFM(etamax_SLFM, 0.0);
+        for(label j = 0 ; j < etamax_SLFM ; j++)
+        {
+            fin>>etaValue_SLFM[j];
+        }
+
+        //blank line
+        fin.getLine(gbg);
+
+        //INITIAL...
+        fin.getLine(gbg);
+
+        //species loop
+        scalarField temp(etamax_SLFM, 0.0);
+        for(label i = 0 ; i < NoSpecies ; i++)
+        {
+            word spName;
+            //H2... (species name)
+            fin >> spName;
+            const label spI (species_[spName]);
+
+            // Info << spI << spName << endl;
+
+            //read species mass fraction (from SLFM library)
+            for(label j=0 ; j<etamax_SLFM ; j++)
+            {
+                scalar Yj = 0.0;
+                fin >> Yj;
+                temp[j] = max(0.0, Yj);
+            }
+
+            Y_[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
+
+            //read reaction rate (from SLFM library)
+            for(label j = 0 ; j < etamax_SLFM ; j++)
+            {
+                fin>>temp[j];
+            }
+
+            W_[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
+
+            //blank line
+            fin.getLine(gbg);
+        }
+
+        //read temperature
+        {
+            //H2... (species name)
+            fin.getLine(gbg);
+
+            //read temperature
+            for(label j = 0 ; j < etamax_SLFM ; j++)
+            {
+                scalar Yj = 0.0;
+                fin >> Yj;
+                temp[j] = max(0.0, Yj);
+            }
+
+            T_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
+
+            //read heat source
+            for(label j = 0 ; j < etamax_SLFM ; j++)
+            {
+                fin>>temp[j];
+            }
+
+            Q_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
+
+            //blank line
+            fin.getLine(gbg);
+        }
+
+        //DENSITY...
+        fin.getLine(gbg);
+
+        //read density (from SLFM library)
+        for(label j = 0 ; j < etamax_SLFM ; j++)
+        {
+            fin>>temp[j];
+        }
+        //Info<<rho_temp<<endl;
+
+        rho_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
+
+        //blank line
+        fin.getLine(gbg);
+
+        //ENTHALPY...
+        fin.getLine(gbg);
+
+        //read enthalpy (from SLFM library)
+        for(label j = 0 ; j < etamax_SLFM ; j++)
+        {
+            fin>>temp[j];
+        }
+
+        h_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::SLFM::FlameStructure::~FlameStructure()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructure.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::SLFM::FlameStructure
+
+Description
+
+SourceFiles
+    FlameStructureI.H
+    FlameStructure.C
+    FlameStructureIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef FlameStructure_H
+#define FlameStructure_H
+
+#include "error.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+class hashedWordList;
+
+namespace SLFM
+{
+    class FlameStructure;
+}
+
+// Forward declaration of friend functions and operators
+Istream& operator>>(Istream&, SLFM::FlameStructure&);
+Ostream& operator<<(Ostream&, const SLFM::FlameStructure&);
+
+namespace SLFM
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class FlameStructure Declaration
+\*---------------------------------------------------------------------------*/
+
+class FlameStructure
+{
+    // Private data
+        typedef List<scalarField> scalarFieldArray1d;
+
+        typedef List<scalarFieldArray1d> scalarFieldArray2d;
+
+        typedef List<scalarFieldArray2d> scalarFieldArray3d;
+
+        //- Species mass fraction
+        word NstDir_;
+
+        //- Species mass fraction
+        const scalarList NstList_;
+
+        //- Species mass fraction
+        const scalarField eta_;
+
+        //- Species mass fraction
+        const hashedWordList &species_;
+
+        //- Species mass fraction
+        scalarFieldArray2d Y_;
+
+        //- Species production rate
+        scalarFieldArray2d W_;
+
+        //- Temperature
+        scalarFieldArray1d T_;
+
+        //- Heat source
+        scalarFieldArray1d Q_;
+
+        //- Enthalpy
+        scalarFieldArray1d h_;
+
+        //- Density
+        scalarFieldArray1d rho_;
+
+        //- Specific gas constant for mean W
+        scalarFieldArray1d Rgas_;
+
+
+    // Private Member Functions
+
+        //- Disallow Construct null
+        FlameStructure();
+
+        //- Disallow default bitwise copy construct
+        // FlameStructure(const FlameStructure&);
+
+        //- Disallow default bitwise assignment
+        // void operator=(const FlameStructure&);
+
+
+public:
+
+    // Static data members
+
+        //- Static data staticData
+        // static const dataType staticData;
+
+
+    // Constructors
+
+        //- Construct from components
+        FlameStructure(const string& NstDir, scalarList NstList, const scalarField &etaGrid, const hashedWordList& spTable);
+
+        //- Construct from Istream
+        // FlameStructure(Istream&);
+
+        //- Construct as copy
+        // FlameStructure(const FlameStructure&);
+
+
+    // Selectors
+
+        //- Select null constructed
+        // static autoPtr<FlameStructure> New();
+
+
+    //- Destructor
+    ~FlameStructure();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Select null constructed
+            scalarFieldArray2d& Y() {return Y_;}
+
+            //- Select null constructed
+            scalarFieldArray2d& W() {return W_;}
+
+            //- Select null constructed
+            scalarFieldArray1d& T() {return T_;}
+
+            //- Select null constructed
+            scalarFieldArray1d& Q() {return Q_;}
+
+            //- Select null constructed
+            scalarFieldArray1d& rho() {return rho_;}
+
+            //- Select null constructed
+            scalarFieldArray1d& h() {return h_;}
+
+            //- Select null constructed
+            scalarFieldArray1d& Rgas() {return Rgas_;}
+
+        // Check
+
+        // Edit
+
+        // Write
+
+
+    // Member Operators
+
+        // void operator=(const FlameStructure&);
+
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, FlameStructure&);
+        friend Ostream& operator<<(Ostream&, const FlameStructure&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace SLFM
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "FlameStructureI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructureI.H

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

applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructureIO.C

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

applications/solvers/combustion/SLFMFoam/Make/files

@@ -0,0 +1,17 @@
+FlameStructure/FlameStructure.C
+FlameStructure/FlameStructureIO.C
+
+./AMC/AMC.C
+
+./BetaFunction/BetaFunction.C
+./BetaFunction/BetaFunctionIO.C
+
+./BetaGrid/BetaGrid.C
+./BetaGrid/BetaGridIO.C
+
+./BetaIntegrator/BetaIntegrator.C
+./BetaIntegrator/BetaIntegratorIO.C
+
+SLFMFoam.C
+
+EXE = $(FOAM_APPBIN)/SLFMFoam

applications/solvers/combustion/SLFMFoam/Make/options

@@ -0,0 +1,34 @@
+EXE_INC = \
+    -I. \
+    -I./AMC \
+    -I./BetaFunction \
+    -I./BetaGrid \
+    -I./BetaIntegrator \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lspecie \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lreactionThermophysicalModels \
+    -lchemistryModel \
+    -lODE \
+    -lcombustionModels \
+    -lfvOptions \
+    -lsampling

applications/solvers/combustion/SLFMFoam/Mixturefraction.H

@@ -0,0 +1,10 @@
+// mixture fraction equation
+fvScalarMatrix mfEqn
+(
+    fvm::div(phi, mf)
+  - fvm::laplacian((1/Sc)*turbulence->mut(), mf) //let 1/Sc = 1.47
+);
+
+mfEqn.relax();
+
+mfEqn.solve();

applications/solvers/combustion/SLFMFoam/MixturefractionVar.H

@@ -0,0 +1,17 @@
+//mixture fraction variance equation
+//cf) mfVar should not be ZERO
+SDR = turbulence->epsilon() * mfVar / turbulence->k();
+volVectorField Gradmf = fvc::grad(mf);
+
+fvScalarMatrix mfVarEqn
+(
+    fvm::div(phi, mfVar)
+  - fvm::laplacian((1/Sc) * turbulence->mut(), mfVar)
+ ==
+    2 * ((1/Sc) * turbulence->mut()) * (Gradmf & Gradmf)
+  - 2 * rho * SDR
+);
+
+mfVarEqn.relax();
+
+mfVarEqn.solve();

applications/solvers/combustion/SLFMFoam/SLFMFoam.C

@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    SLFMFoam
+
+Description
+    Steady Laminar Flamelet Model(SLFM) solver for turbulent combustion. SLFM
+    assumes a turbulent flame composed of thin stretched laminar flamelets.
+    Each flamelet library is generated by external program and imported by
+    the solver.
+
+References
+    A.Y. Klimenko, R.W. Bilger, Progress in Energy and Combustion Science 25 (1999) 595-687
+    N. Peters, Turbulent Combustion, Cambridge University Press (2000)
+
+Contact
+    POSTECH combustion lab.
+    huh@postech.ac.kr
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulentFluidThermoModel.H"
+#include "rhoCombustionModel.H"
+#include "fvOptions.H"
+#include "simpleControl.H"
+#include "FlameStructure/FlameStructure.H"
+#include "interpolateXY.H"
+
+#include "AMC.H"
+#include "BetaFunction.H"
+#include "BetaGrid.H"
+#include "BetaIntegrator.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    turbulence->validate();
+
+    //SLFM
+    Info<<"Read SLFM library"<<endl;
+    #include "readSLFMlib.H" //read SLFM library
+
+    Info<<"Construct Beta-PDF"<<endl;
+    #include "preIntegration.H" //AMC C1 coefficient pdf integration
+
+    Info<< "\nStarting time loop\n" << endl;
+
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity SIMPLE corrector loop
+        {
+            #include "UEqn.H"
+            #include "Mixturefraction.H" //mean mixture fraction
+            #include "MixturefractionVar.H" //mean mixture fraction variance
+            #include "setSDR.H" //calculate scalar dissipation rate
+            #include "updateT_RHO.H" //update temperature and density
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        if (runTime.write())
+        {
+            #include "updateYi.H" //update species
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/SLFMFoam/UEqn.H

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

applications/solvers/combustion/SLFMFoam/createFieldRefs.H

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

applications/solvers/combustion/SLFMFoam/createFields.H

@@ -0,0 +1,374 @@
+#include "readGravitationalAcceleration.H"
+
+Info<< "Creating combustion model\n" << endl;
+
+autoPtr<combustionModels::rhoCombustionModel> combustion
+(
+    combustionModels::rhoCombustionModel::New(mesh)
+);
+
+rhoReactionThermo& thermo = combustion->thermo();
+thermo.validate(args.executable(), "h", "e");
+
+basicSpecieMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+const word inertSpecie(thermo.lookup("inertSpecie"));
+
+if (!composition.contains(inertSpecie))
+{
+    FatalErrorInFunction
+        << "Specified inert specie '" << inertSpecie << "' not found in "
+        << "species list. Available species:" << composition.species()
+        << exit(FatalError);
+}
+
+volScalarField& p = thermo.p();
+volScalarField& T = thermo.T();
+
+wordList rhoBoundaryTypes
+(
+	T.boundaryField().size(),
+	zeroGradientFvPatchScalarField::typeName
+);
+
+forAll(composition.Y(inertSpecie).boundaryField(), patchi)
+{
+	rhoBoundaryTypes[patchi]
+        = composition.Y(inertSpecie).boundaryField()[patchi].type();
+}
+
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho(),
+    rhoBoundaryTypes
+);
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+mesh.setFluxRequired(p.name());
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        simple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        simple.dict(),
+        dimDensity,
+        0
+    )
+);
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+// Set the turbulence into the combustion model
+combustion->setTurbulence(turbulence());
+
+Info<< "Creating multi-variate interpolation scheme\n" << endl;
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+
+// SLFM-related fields
+
+Info<<"Creating field mf, mfVar, SDR\n"<<endl;
+
+volScalarField mf //mixture fraction
+(
+    IOobject
+    (
+        "mf",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+volScalarField mfVar //mixture fraction variance
+(
+    IOobject
+    (
+        "mfVar",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+volScalarField SDR //scalar dissipation rate
+(
+    IOobject
+    (
+        "SDR",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("SDR", dimEnergy/dimTime/dimEnergy, 0.0)
+);
+
+volScalarField jlc //eta-index, lower
+(
+    IOobject
+    (
+        "jlc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("jlc", dimless, 0.0)
+);
+
+volScalarField jhc //eta-index, higher
+(
+    IOobject
+    (
+        "jhc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("jhc", dimless, 0.0)
+);
+
+volScalarField jfc //eta-factor
+(
+    IOobject
+    (
+        "jfc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("jfc", dimless, 0.0)
+);
+
+volScalarField vlc //var-index, lower
+(
+    IOobject
+    (
+        "vlc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("vlc", dimless, 0.0)
+);
+
+volScalarField vhc //var-index, higher
+(
+    IOobject
+    (
+        "vhc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("vhc", dimless, 0.0)
+);
+
+volScalarField vfc //var-factor
+(
+    IOobject
+    (
+        "vfc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("vfc", dimless, 0.0)
+);
+
+volScalarField nlc //sdr-index, lower
+(
+    IOobject
+    (
+        "nlc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("nlc", dimless, 0.0)
+);
+
+volScalarField nhc //sdr-index, higher
+(
+    IOobject
+    (
+        "nhc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("nhc", dimless, 0.0)
+);
+
+volScalarField nfc //sdr-factor
+(
+    IOobject
+    (
+        "nfc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("nfc", dimless, 0.0)
+);
+
+//SLFM property read
+
+IOdictionary SLFMdict //SLFM properties dictionary
+(
+    IOobject
+    (
+        "SLFMdict",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+const wordList postSpecieNames(SLFMdict.lookup("postSpecieNames"));
+labelList postSpecieIndices(postSpecieNames.size(),-1);
+forAll(postSpecieNames, i)
+{
+    const word &spName (postSpecieNames[i]);
+    const label spI (composition.species()[spName]);
+    postSpecieIndices[i] = spI;
+}
+scalar stoiMF(readScalar(SLFMdict.lookup("stoiMF")));
+scalar stretchFac(readScalar(SLFMdict.lookup("stretchFactor")));
+scalar lowerN(readScalar(SLFMdict.lookup("lowerN")));
+scalar upperN(readScalar(SLFMdict.lookup("upperN")));
+scalar Sc(readScalar(SLFMdict.lookup("Sc")));
+const label etamax(lowerN+upperN+1);
+word outletName(SLFMdict.lookup("outletName"));
+
+labelList nIntervals (SLFMdict.lookup("detailedN"));
+
+CMC::BetaGrid bg(SLFMdict);
+const scalarField &etaGrid(bg.etaSpace());
+const scalarField amcGridValue(CMC::AMC(etaGrid));
+
+Info << etaGrid << endl;
+
+scalarField etaValue(etamax,0);
+scalarField etaIndex(etamax,0);
+scalarField deta_h(etamax,0);
+scalarField deta_l(etamax,0);
+scalarField deta_m(etamax,0);
+
+#include "setEtaSpace.H"
+
+Info<<"\nNumber of eta-point = "<<etamax<<endl;
+
+
+scalar NVar(readScalar(SLFMdict.lookup("NVar")));
+scalarField varValue(NVar+1,0);
+scalarField varIndex(NVar+1,0);
+
+#include "setVarSpace.H"
+
+
+
+PtrList<volScalarField> Neta(etamax);
+
+for(label j=0 ; j<etamax ; j++)
+{
+
+
+    word Netai = "Neta_"+Foam::name(j);
+    Neta.set
+    (
+        j,
+        new volScalarField
+        (
+            IOobject
+            (
+                Netai,
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(Netai, dimless, 0.0)
+        )
+    );
+}
+
+#include "createMRF.H"

applications/solvers/combustion/SLFMFoam/pEqn.H

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

applications/solvers/combustion/SLFMFoam/preIntegration.H

@@ -0,0 +1,78 @@
+scalarFieldArray1d C1table //AMC C1 coefficient table
+(
+    etamax,
+    scalarField(NVar+1, 0)
+);
+
+scalarFieldArray2d Ttable
+(
+    etamax,
+    scalarFieldArray1d
+    (
+        NVar+1,
+        scalarField(NstList.size(), 0)
+    )
+);
+
+scalarFieldArray2d Rtable
+(
+    etamax,
+    scalarFieldArray1d
+    (
+        NVar+1,
+        scalarField(NstList.size(), 0)
+    )
+);
+
+scalarFieldArray3d Ytable
+(
+    Y.size(),
+    scalarFieldArray2d
+    (
+        etamax,
+        scalarFieldArray1d
+        (
+            NVar+1,
+            scalarField(NstList.size(), 0)
+        )
+    )
+);
+
+Info<<"Construct C1 coefficient table"<<endl;
+for(label j = 1; j < etamax - 1; j++)
+{
+    for(label v = 1; v < varValue.size() - 1; v++)
+    {
+        const scalar mf = etaValue[j];
+        const scalar mfVar = varValue[v]*mf*(1.0-mf);
+
+        CMC::BetaFunction bf(mf, mfVar);
+        CMC::BetaIntegrator bi(bf, bg);
+
+        C1table[j][v] = 1.0/bi.betaIntegrate(amcGridValue);
+    }
+}
+
+Info<<"Construct Temp and Rgas table"<<endl;
+for(label j=0 ; j<etamax ; j++)
+{
+    for(label v=0 ; v<=NVar ; v++)
+    {
+        const scalar mf = etaValue[j];
+        const scalar mfVar = varValue[v]*mf*(1.0-mf);
+
+        CMC::BetaFunction bf(mf, mfVar);
+        CMC::BetaIntegrator bi(bf, bg);
+
+        for(label n=0 ; n<NstList.size() ; n++)
+        {
+            Ttable[j][v][n] = bi.betaIntegrate(T_SLFM[n]);
+            Rtable[j][v][n] = bi.betaIntegrate(Rgas_SLFM[n]);
+            forAll(Y, yi)
+            {
+                Ytable[yi][j][v][n] = bi.betaIntegrate(Y_SLFM[yi][n]);
+            }
+        }
+    }
+    Info<<scalar(j+1)/etamax*100<<"%"<<endl;
+}

applications/solvers/combustion/SLFMFoam/readSLFMlib.H

@@ -0,0 +1,33 @@
+Switch Equilibrium = SLFMdict.lookup("Equilibrium");
+word NstFolder = SLFMdict.lookup("NstFolder");
+scalarField NstList = SLFMdict.lookup("NstList");
+scalarField NstIndex(NstList.size(),0);
+
+forAll(NstList, n)
+{
+    NstIndex[n] = n;
+}
+
+typedef List<scalarField> scalarFieldArray1d;
+typedef List<scalarFieldArray1d> scalarFieldArray2d;
+typedef List<scalarFieldArray2d> scalarFieldArray3d;
+
+SLFM::FlameStructure slfmLibrary(runTime.constant()/NstFolder, NstList, etaGrid, composition.species());
+
+scalarFieldArray2d& Y_SLFM (slfmLibrary.Y());
+scalarFieldArray2d& W_SLFM (slfmLibrary.W());
+scalarFieldArray1d& T_SLFM (slfmLibrary.T());
+scalarFieldArray1d& Q_SLFM (slfmLibrary.Q());
+scalarFieldArray1d& rho_SLFM (slfmLibrary.rho());
+scalarFieldArray1d& h_SLFM (slfmLibrary.h());
+scalarFieldArray1d& Rgas_SLFM (slfmLibrary.Rgas());
+
+for(label n=0 ; n<NstList.size() ; n++)
+{
+    for(label i=0 ; i<Y.size() ; i++)
+    {
+        Rgas_SLFM[n] += Y_SLFM[i][n]/composition.W(i); // = 1/meanMW
+    }
+    Rgas_SLFM[n] = 8314.46*Rgas_SLFM[n]; //universal gas constant = 8314.46 [J/K.Kmole]
+
+}

applications/solvers/combustion/SLFMFoam/setEtaSpace.H

@@ -0,0 +1,51 @@
+//set eta-space
+scalar upperSum(0), upperDelta(0);
+scalar lowerSum(0), lowerDelta(0);
+
+//summation of geometric progression (eta > stoi)
+upperSum = Foam::pow(stretchFac,upperN)-1.0;
+upperSum = upperSum/(stretchFac-1.0);
+
+//summation of geometric progression (eta < stoi)
+lowerSum = Foam::pow(stretchFac,lowerN)-1.0;
+lowerSum = lowerSum/(stretchFac-1.0);
+
+//set reference delta eta
+upperDelta = (1.0-stoiMF)/upperSum;
+lowerDelta = stoiMF/lowerSum;
+
+//set boundary and stoichiometric value of eta
+etaValue[0] = 0.0;
+etaValue[lowerN] = stoiMF;
+etaValue.last() = 1.0;
+
+//set internal value of eta (eta > stoi)
+for(label j=0 ; j<upperN-1 ; j++)
+{
+    etaValue[lowerN+j+1] = etaValue[lowerN+j] + upperDelta*pow(stretchFac,j);
+}
+
+//set internal value of eta (eta < stoi)
+for(label j=0 ; j<lowerN-1 ; j++)
+{
+    etaValue[lowerN-j-1] = etaValue[lowerN-j] - lowerDelta*pow(stretchFac,j);
+}
+
+Info<<etaValue<<endl;
+
+//calculate delta eta
+for(label j=1 ; j<etamax-1 ; j++)
+{
+    deta_h[j] = etaValue[j+1] - etaValue[j];
+    deta_l[j] = etaValue[j] - etaValue[j-1];
+    deta_m[j] = 0.5*(etaValue[j+1] - etaValue[j-1]);
+}
+
+//Info<<deta_h<<endl;
+//Info<<deta_l<<endl;
+//Info<<deta_m<<endl;
+
+for(label j=0 ; j<etamax ; j++)
+{
+    etaIndex[j] = j;
+}

applications/solvers/combustion/SLFMFoam/setSDR.H

@@ -0,0 +1,79 @@
+scalarField AMCcoeff(CMC::AMC(etaValue));
+
+//cell loop
+
+forAll(mf, cellI)
+{
+
+
+    scalar jl(0), jh(0), vl(0), vh(0);
+    scalar jfac(0), vfac(0);
+    scalar C1coeff(0);
+
+    //find eta-index, factor
+    if(mf[cellI] < etaValue[1])
+    {
+        jl = 0;
+        jh = 1;
+    }
+    else if(mf[cellI] > etaValue[etamax-2])
+    {
+        jl = etamax-2;
+        jh = etamax-1;
+    }
+    else
+    {
+        jl = label( interpolateXY(mf[cellI], etaValue, etaIndex) );
+        jh = jl+1;
+    }
+    jfac = max(0.0, (mf[cellI]-etaValue[jl])/(etaValue[jh]-etaValue[jl]));
+
+    //find var-index, factor
+    scalar scaledVar = min(0.99999, mfVar[cellI]/(mf[cellI]*(1.0-mf[cellI])+SMALL));
+    if(scaledVar < varValue[1])
+    {
+        vl = 0;
+        vh = 1;
+    }
+    else if(scaledVar > varValue[NVar-1])
+    {
+        vl = NVar-1;
+        vh = NVar;
+    }
+    else
+    {
+        vl = label( interpolateXY(scaledVar, varValue, varIndex) );
+        vh = vl+1;
+    }
+    vfac = max(0.0, (scaledVar-varValue[vl])/(varValue[vh]-varValue[vl]));
+
+
+    //Bi-linear interpolation on j and v
+    //Numerical recipes, 2nd Ed. p.117
+    C1coeff
+    = C1table[jl][vl]*(1-jfac)*(1-vfac)
+    + C1table[jh][vl]*(jfac)*(1-vfac)
+    + C1table[jh][vh]*(jfac)*(vfac)
+    + C1table[jl][vh]*(1-jfac)*(vfac);
+
+    jlc[cellI] = jl;
+    jhc[cellI] = jh;
+    jfc[cellI] = jfac;
+    vlc[cellI] = vl;
+    vhc[cellI] = vh;
+    vfc[cellI] = vfac;
+
+
+    scalar Coeff = AMCcoeff[lowerN];
+
+    if (Equilibrium == true)
+    {
+        Neta[lowerN][cellI] = 0; //Equilibrium
+    }
+    else
+    {
+        Neta[lowerN][cellI] = SDR[cellI]*Coeff*C1coeff; //CSDR at stoichiometic m.f.
+    }
+}
+
+

applications/solvers/combustion/SLFMFoam/setVarSpace.H

@@ -0,0 +1,23 @@
+//set normalized variance space
+{
+    scalar sum(0), del(0), fac(1.075);
+
+    sum = Foam::pow(fac,NVar) - 1.0;
+    sum = sum/(fac-1.0);
+
+    del = 1.0/sum;
+
+    varValue[0] = 0.0;
+    varValue[NVar] = 1.0;
+
+    for(label v=0 ; v<NVar-1 ; v++)
+    {
+        varValue[v+1] = varValue[v] + del*Foam::pow(fac,v);
+    }
+
+    for(label v=0 ; v<=NVar ; v++)
+    {
+        varIndex[v] = v;
+    }
+    Info<<varValue<<endl;
+}

applications/solvers/combustion/SLFMFoam/updateT_RHO.H

@@ -0,0 +1,60 @@
+//force the saving of the old-time values
+//important!!
+
+rho.oldTime();
+
+forAll(rho, cellI)
+{
+    scalar& Nst = Neta[lowerN][cellI];
+
+    //find Nst index
+    if(Nst < NstList.first())
+    {
+        nlc[cellI] = 0;
+        nfc[cellI] = 0;
+    }
+    else if(Nst > NstList[NstList.size()-2])
+    {
+        nhc[cellI] = NstList.size()-2; //just before extinction
+        nfc[cellI] = 1;
+    }
+    else
+    {
+        nlc[cellI] = label( interpolateXY(Nst, NstList, NstIndex) );
+        nhc[cellI] = nlc[cellI]+1;
+        nfc[cellI] = (Nst-NstList[nlc[cellI]])
+                     /(NstList[nhc[cellI]]-NstList[nlc[cellI]]);
+    }
+
+    //update Favre mean temp. and rho(=1/RT)
+    scalar jl = jlc[cellI], jh = jhc[cellI];
+    scalar vl = vlc[cellI], vh = vhc[cellI];
+    scalar nl = nlc[cellI], nh = nhc[cellI];
+    scalar jfac = jfc[cellI], vfac = vfc[cellI], nfac = nfc[cellI];
+
+    scalar T00 = Ttable[jl][vl][nl]*(1-jfac)+Ttable[jh][vl][nl]*jfac;
+    scalar T01 = Ttable[jl][vl][nh]*(1-jfac)+Ttable[jh][vl][nh]*jfac;
+    scalar T10 = Ttable[jl][vh][nl]*(1-jfac)+Ttable[jh][vh][nl]*jfac;
+    scalar T11 = Ttable[jl][vh][nh]*(1-jfac)+Ttable[jh][vh][nh]*jfac;
+
+    scalar T0 = T00*(1-vfac)+T10*vfac;
+    scalar T1 = T01*(1-vfac)+T11*vfac;
+
+    T[cellI] = T0*(1-nfac)+T1*nfac;
+
+    scalar R00 = Rtable[jl][vl][nl]*(1-jfac)+Rtable[jh][vl][nl]*jfac;
+    scalar R01 = Rtable[jl][vl][nh]*(1-jfac)+Rtable[jh][vl][nh]*jfac;
+    scalar R10 = Rtable[jl][vh][nl]*(1-jfac)+Rtable[jh][vh][nl]*jfac;
+    scalar R11 = Rtable[jl][vh][nh]*(1-jfac)+Rtable[jh][vh][nh]*jfac;
+
+    scalar R0 = R00*(1-vfac)+R10*vfac;
+    scalar R1 = R01*(1-vfac)+R11*vfac;
+
+    scalar Rgas = R0*(1-nfac)+R1*nfac;
+    rho[cellI] = p[cellI]/(Rgas*T[cellI]);
+}
+
+//correct processor boundary value of T and rho
+T.correctBoundaryConditions();
+rho.correctBoundaryConditions();
+

applications/solvers/combustion/SLFMFoam/updateYi.H

@@ -0,0 +1,56 @@
+//force the saving of the old-time values
+//important!!
+
+forAll(rho, cellI)
+{
+    scalar& Nst = Neta[lowerN][cellI];
+
+    //find Nst index
+    if(Nst < NstList.first())
+    {
+        nlc[cellI] = 0;
+        nfc[cellI] = 0;
+    }
+    else if(Nst > NstList[NstList.size()-2])
+    {
+        nhc[cellI] = NstList.size()-2; //just before extinction
+        nfc[cellI] = 1;
+    }
+    else
+    {
+        nlc[cellI] = label( interpolateXY(Nst, NstList, NstIndex) );
+        nhc[cellI] = nlc[cellI]+1;
+        nfc[cellI] = (Nst-NstList[nlc[cellI]])
+                     /(NstList[nhc[cellI]]-NstList[nlc[cellI]]);
+    }
+
+    //update Favre mean temp. and rho(=1/RT)
+    scalar jl = jlc[cellI], jh = jhc[cellI];
+    scalar vl = vlc[cellI], vh = vhc[cellI];
+    scalar nl = nlc[cellI], nh = nhc[cellI];
+    scalar jfac = jfc[cellI], vfac = vfc[cellI], nfac = nfc[cellI];
+
+
+    forAll(Y, yi)
+    {
+        scalar y00 = Ytable[yi][jl][vl][nl]*(1-jfac)+Ytable[yi][jh][vl][nl]*jfac;
+        scalar y01 = Ytable[yi][jl][vl][nh]*(1-jfac)+Ytable[yi][jh][vl][nh]*jfac;
+        scalar y10 = Ytable[yi][jl][vh][nl]*(1-jfac)+Ytable[yi][jh][vh][nl]*jfac;
+        scalar y11 = Ytable[yi][jl][vh][nh]*(1-jfac)+Ytable[yi][jh][vh][nh]*jfac;
+
+        scalar y0 = y00*(1-vfac)+y10*vfac;
+        scalar y1 = y01*(1-vfac)+y11*vfac;
+
+        scalar ygas = y0*(1-nfac)+y1*nfac;
+
+        Y[yi][cellI] = ygas;
+    }
+}
+
+//correct processor boundary value of T and rho
+forAll(Y, yi)
+{
+    Y[yi].correctBoundaryConditions();
+}
+
+

applications/solvers/combustion/postSootFoam/Make/files

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

applications/solvers/lagrangian/ThermoMPPICFoam/EEqn.H

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

applications/solvers/lagrangian/ThermoMPPICFoam/Make/files

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

applications/solvers/lagrangian/ThermoMPPICFoam/Make/options

@@ -1,7 +1,7 @@
 EXE_INC = \
+    -I. \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
@@ -18,38 +18,39 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/properties/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/combustionModels/lnInclude
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
 
 EXE_LIBS = \
-    -lchemistryModel \
-    -lcombustionModels \
-    -lradiationModels \
-    -lradiationSootModels \
-    -llagrangian \
-    -llagrangianIntermediate \
-    -llagrangianTurbulence \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -lreactionThermophysicalModels \
     -lfiniteVolume \
     -lmeshTools \
     -lturbulenceModels \
     -lcompressibleTurbulenceModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lspecie \
     -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
     -lliquidProperties \
     -lliquidMixtureProperties \
     -lsolidProperties \
     -lsolidMixtureProperties \
     -lthermophysicalFunctions \
+    -lreactionThermophysicalModels \
     -lSLGThermo \
+    -lchemistryModel \
+    -lradiationModels \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
-    -lfvOptions
+    -lcombustionModels \
-
+    -lfvOptions \
+    -lsampling

applications/solvers/lagrangian/ThermoMPPICFoam/ThermoMPPICFoam.C

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

applications/solvers/lagrangian/ThermoMPPICFoam/UEqn.H

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

applications/solvers/lagrangian/ThermoMPPICFoam/YEqn.H

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

applications/solvers/lagrangian/ThermoMPPICFoam/createClouds.H

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

applications/solvers/lagrangian/ThermoMPPICFoam/createFieldRefs.H

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

applications/solvers/lagrangian/ThermoMPPICFoam/createFields.H

@@ -1,18 +1,31 @@
-Info<< "\nCreating thermophysical model for Fluid mesh\n" << endl;
 #include "readGravitationalAcceleration.H"
 
+Info<< "Creating combustion model\n" << endl;
+
 autoPtr<combustionModels::rhoCombustionModel> combustion
 (
     combustionModels::rhoCombustionModel::New(mesh)
 );
 
 rhoReactionThermo& thermo = combustion->thermo();
-thermo.validate(args.executable(), "ha");
+thermo.validate(args.executable(), "h", "e");
+
 SLGThermo slgThermo(mesh, thermo);
+
 basicSpecieMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 
-word inertSpecie(thermo.lookup("inertSpecie"));
+const word inertSpecie(thermo.lookup("inertSpecie"));
+
+if (!composition.contains(inertSpecie))
+{
+    FatalErrorInFunction
+        << "Specified inert specie '" << inertSpecie << "' not found in "
+        << "species list. Available species:" << composition.species()
+        << exit(FatalError);
+}
+
+volScalarField& p = thermo.p();
 
 volScalarField rho
 (
@@ -27,7 +40,7 @@ volScalarField rho
     thermo.rho()
 );
 
-Info<< "Reading field U\n" << endl;
+Info<< "\nReading field U\n" << endl;
 volVectorField U
 (
     IOobject
@@ -41,16 +54,16 @@ volVectorField U
     mesh
 );
 
-
-volScalarField& p = thermo.p();
-
 #include "compressibleCreatePhi.H"
+
+mesh.setFluxRequired(p.name());
+
 dimensionedScalar rhoMax
 (
     dimensionedScalar::lookupOrDefault
     (
         "rhoMax",
-        simple.dict(),
+        pimple.dict(),
         dimDensity,
         GREAT
     )
@@ -61,13 +74,13 @@ dimensionedScalar rhoMin
     dimensionedScalar::lookupOrDefault
     (
         "rhoMin",
-        simple.dict(),
+        pimple.dict(),
         dimDensity,
         0
     )
 );
 
-Info << "Creating turbulence model.\n" << nl;
+Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
     compressible::turbulenceModel::New
@@ -79,7 +92,7 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-// Set the turbulence into the reaction model
+// Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
 Info<< "Creating field dpdt\n" << endl;
@@ -98,6 +111,7 @@ volScalarField dpdt
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
+Info<< "Creating multi-variate interpolation scheme\n" << endl;
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
 forAll(Y, i)
@@ -106,37 +120,39 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-/*
+volScalarField dQ
-volScalarField Hha //absolute Enthalpy
 (
     IOobject
     (
-        "Hha",
+        "dQ",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("Hha", dimEnergy, 0.0)
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
-volScalarField Hla //absolute Enthalpy
-(
-    IOobject
-    (
-        "Hla",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("Hla", dimEnergy, 0.0)
-);
-*/
 
-#include "createRadiationModel.H"
+volScalarField alphac
+(
+    IOobject
+    (
+        "alphac",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+        mesh,
+        dimensionedScalar("0", dimless, 1)
+);
+
+surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
+surfaceScalarField alphacPhi("alphacPhi", alphacf*phi);
+
+volScalarField alphacrho("alphacrho", alphac*rho);
 
 #include "createMRF.H"
+#include "createRadiationModel.H"
 #include "createClouds.H"
-

applications/solvers/lagrangian/ThermoMPPICFoam/pEqn.H

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

applications/solvers/lagrangian/ThermoMPPICFoam/rhoEqn.H

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

applications/solvers/lagrangian/ThermoMPPICFoam/setRDeltaT.H

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/EEqn.H

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/files

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/options

@@ -1,30 +1,57 @@
 EXE_INC = \
+    -I. \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
-LIB_LIBS = \
+
-    -lradiationModels \
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lcoalCombustion \
+    -llagrangianTurbulence \
+    -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lspecie \
+    -lliquidProperties \
-    -lsolidThermo \
-    -lSLGThermo \
-    -lsolidMixtureProperties \
     -lliquidMixtureProperties \
     -lsolidProperties \
-    -lliquidProperties \
+    -lsolidMixtureProperties \
-    -lfiniteVolume \
+    -lthermophysicalFunctions \
-    -lmeshTools
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lradiationModels \
+    -lODE \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lcombustionModels \
+    -lfvOptions \
+    -lsampling

applications/solvers/lagrangian/simpleCoalcombustionFoam/UEqn.H

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/YEqn.H

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/createClouds.H

@@ -1,7 +1,6 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingMultiphaseCloud parcels
+coalCloudList parcels
 (
-    "reactingCloud1",
     rho,
     U,
     g,

applications/solvers/lagrangian/simpleCoalcombustionFoam/createFieldRefs.H

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/createFields.H

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/pEqn.H

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

applications/solvers/lagrangian/simpleCoalcombustionFoam/simpleCoalcombustionFoam.C

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    simpleReactingParcelFoam
+    simpleCoalcombustionFoam
 
 Description
     Steady state solver for compressible, turbulent flow with reacting,
@@ -32,9 +32,10 @@ Description
 
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicReactingMultiphaseCloud.H"
+#include "coalCloudList.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
+#include "IOporosityModelList.H"
 #include "fvOptions.H"
 #include "SLGThermo.H"
 #include "simpleControl.H"
@@ -43,24 +44,38 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "postProcess.H"
 
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-
     #include "createControl.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
 
     turbulence->validate();
 
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
     Info<< "\nStarting time loop\n" << endl;
 
     while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        radiation->correct();
+        parcels.evolve();
+
+        // --- Pressure-velocity SIMPLE corrector loop
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "EEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
 
         runTime.write();
 

src/Allwmake

@@ -64,7 +64,6 @@ wmake $targetType randomProcesses
 transportModels/Allwmake $targetType $*
 thermophysicalModels/Allwmake $targetType $*
 TurbulenceModels/Allwmake $targetType $*
-wmake $targetType thermophysicalModels/radiationSootModel
 wmake $targetType combustionModels
 regionModels/Allwmake $targetType $*
 lagrangian/Allwmake $targetType $*

src/combustionModels/EDM/EDM.C

@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "EDM.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+template<class Type>
+Foam::combustionModels::EDM<Type>::EDM
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const word& phaseName
+)
+:
+    Type(modelType, mesh, phaseName),
+    finiteRate_ (this->coeffs().lookupOrDefault("finiteRate", false)),
+    A_ (this->coeffs().lookupOrDefault("A", 4.0)),
+    B_ (this->coeffs().lookupOrDefault("B", 0.5))
+{
+    if (finiteRate_)
+    {
+        Info<< "    using Finite-rate/Eddy Dissipation Model" << endl;
+    }
+    else
+    {
+        Info<< "    using Eddy Dissipation Model" << endl;
+    }
+    Info<< "    A = " << A_ << endl;
+    Info<< "    B = " << B_ << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::combustionModels::EDM<Type>::~EDM()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::EDM<Type>::tc() const
+{
+    return this->chemistryPtr_->tc();
+}
+
+
+template<class Type>
+void Foam::combustionModels::EDM<Type>::correct()
+{
+    if (this->active())
+    {
+        tmp<volScalarField> tk(this->turbulence().k());
+        const volScalarField& k = tk();
+
+        tmp<volScalarField> tepsilon(this->turbulence().epsilon());
+        const volScalarField& epsilon = tepsilon();
+
+        this->chemistryPtr_->edm(epsilon/k, A_, B_);
+
+        if (finiteRate_)
+        {
+            // Copy edm rate stored in chemistryPtr_->RR(i)
+            PtrList<DimensionedField<scalar, volMesh>> edmRates
+            (
+                this->thermo().composition().Y().size()
+            );
+            forAll(edmRates, i)
+            {
+                edmRates.set
+                (
+                    i,
+                    new DimensionedField<scalar, volMesh>
+                    (
+                        this->chemistryPtr_->RR(i)
+                    )
+                );
+            }
+
+            this->chemistryPtr_->calculate();
+
+            forAll(edmRates, i)
+            {
+                DimensionedField<scalar, volMesh> &finiteRatei
+                (
+                    this->chemistryPtr_->RR(i)
+                );
+
+                forAll(finiteRatei, celli)
+                {
+                    finiteRatei[celli] = min(finiteRatei[celli], edmRates[i][celli]);
+                    finiteRatei[celli] =
+                        finiteRatei[celli]*edmRates[i][celli]
+                        / (finiteRatei[celli] + edmRates[i][celli] + SMALL);
+                }
+
+            }
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::combustionModels::EDM<Type>::R(volScalarField& Y) const
+{
+    tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
+
+    fvScalarMatrix& Su = tSu.ref();
+
+    if (this->active())
+    {
+        const label specieI =
+            this->thermo().composition().species()[Y.member()];
+
+        Su += this->chemistryPtr_->RR(specieI);
+    }
+
+    return tSu;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::EDM<Type>::dQ() const
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName(typeName + ":dQ", this->phaseName_),
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+        )
+    );
+
+    if (this->active())
+    {
+        tdQ.ref() = this->chemistryPtr_->dQ();
+    }
+
+    return tdQ;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::EDM<Type>::Sh() const
+{
+    tmp<volScalarField> tSh
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName(typeName + ":Sh", this->phaseName_),
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
+        )
+    );
+
+    if (this->active())
+    {
+        tSh.ref() = this->chemistryPtr_->Sh();
+    }
+
+    return tSh;
+}
+
+
+template<class Type>
+bool Foam::combustionModels::EDM<Type>::read()
+{
+    if (Type::read())
+    {
+        this->coeffs().lookup("finiteRate") >> finiteRate_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/combustionModels/EDM/EDM.H

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::EDM
+
+Description
+    Eddy dissipation model with finite-rate chemistry.
+                                                                    28.Feb.2017
+                                                        Combustion Lab. POSTECH
+                                                                      Karam Han
+                                              (Version upgraded by Jinwoo Park)
+
+SourceFiles
+    EDM.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef EDM_H
+#define EDM_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class EDM Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class EDM
+:
+    public Type
+{
+    // Private data
+
+        //- finite-rate/EDM using Arrhenius formula
+        bool finiteRate_;
+
+        //- EDM Parameter A
+        scalar A_;
+
+        //- EDM Parameter B
+        scalar B_;
+
+protected:
+
+    // Protected Member Functions
+
+        //- Return the chemical time scale
+        tmp<volScalarField> tc() const;
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        EDM(const EDM&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const EDM&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("EDM");
+
+
+    // Constructors
+
+        //- Construct from components
+        EDM
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const word& phaseName
+        );
+
+
+    //- Destructor
+    virtual ~EDM();
+
+
+    // 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 "EDM.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/EDM/EDMs.C

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

src/combustionModels/Make/files

@@ -13,6 +13,8 @@ rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
 diffusion/diffusions.C
 infinitelyFastChemistry/infinitelyFastChemistrys.C
 
+EDM/EDMs.C
+
 PaSR/PaSRs.C
 
 laminar/laminars.C

src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H

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

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

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "COxidationKD_CO.H"
+#include "mathematicalConstants.H"
+#include "thermodynamicConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationKD_CO<CloudType>::
+COxidationKD_CO
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    SurfaceReactionModel<CloudType>(dict, owner, typeName),
+    Sb_(readScalar(this->coeffDict().lookup("Sb"))),
+    C1_(readScalar(this->coeffDict().lookup("C1"))),
+    C2_(readScalar(this->coeffDict().lookup("C2"))),
+    E_(readScalar(this->coeffDict().lookup("E"))),
+    CsLocalId_(-1),
+    O2GlobalId_(owner.composition().carrierId("O2")),
+    COGlobalId_(owner.composition().carrierId("CO")),
+    WC_(0.0),
+    WO2_(0.0),
+    HcCO_(0.0)
+{
+    // Determine Cs ids
+    label idSolid = owner.composition().idSolid();
+    CsLocalId_ = owner.composition().localId(idSolid, "C");
+
+    // Set local copies of thermo properties
+    WO2_ = owner.thermo().carrier().W(O2GlobalId_);
+    const scalar WCO = owner.thermo().carrier().W(COGlobalId_);
+    WC_ = WCO - 0.5*WO2_;
+
+    HcCO_ = owner.thermo().carrier().Hc(COGlobalId_);
+
+    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
+    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
+    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
+}
+
+
+template<class CloudType>
+Foam::COxidationKD_CO<CloudType>::
+COxidationKD_CO
+(
+    const COxidationKD_CO<CloudType>& srm
+)
+:
+    SurfaceReactionModel<CloudType>(srm),
+    Sb_(srm.Sb_),
+    C1_(srm.C1_),
+    C2_(srm.C2_),
+    E_(srm.E_),
+    CsLocalId_(srm.CsLocalId_),
+    O2GlobalId_(srm.O2GlobalId_),
+    COGlobalId_(srm.COGlobalId_),
+    WC_(srm.WC_),
+    WO2_(srm.WO2_),
+    HcCO_(srm.HcCO_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationKD_CO<CloudType>::
+~COxidationKD_CO()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::COxidationKD_CO<CloudType>::calculate
+(
+    const scalar dt,
+    const label cellI,
+    const scalar d,
+    const scalar T,
+    const scalar Tc,
+    const scalar pc,
+    const scalar rhoc,
+    const scalar mass,
+    const scalarField& YGas,
+    const scalarField& YLiquid,
+    const scalarField& YSolid,
+    const scalarField& YMixture,
+    const scalar N,
+    scalarField& dMassGas,
+    scalarField& dMassLiquid,
+    scalarField& dMassSolid,
+    scalarField& dMassSRCarrier
+) const
+{
+    // Fraction of remaining combustible material
+    const label idSolid = CloudType::parcelType::SLD;
+    const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
+
+    // Surface combustion active combustible fraction is consumed
+    if (fComb < SMALL)
+    {
+        return 0.0;
+    }
+
+    const SLGThermo& thermo = this->owner().thermo();
+
+    // Local mass fraction of O2 in the carrier phase
+    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
+
+    // Diffusion rate coefficient
+    const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
+
+    // Kinetic rate
+    const scalar Rk = C2_*exp(-E_/(constant::thermodynamic::RR*T)); //
+
+    // Particle surface area
+    const scalar Ap = constant::mathematical::pi*sqr(d);
+
+    // Change in C mass [kg]
+    scalar dmC = Ap*rhoc*constant::thermodynamic::RR*Tc*YO2/WO2_*D0*Rk/(D0 + Rk)*dt;
+
+    // Limit mass transfer by availability of C
+    dmC = min(mass*fComb, dmC);
+
+    // Molar consumption
+    const scalar dOmega = dmC/WC_;
+
+    // Change in O2 mass [kg]
+    const scalar dmO2 = dOmega*Sb_*WO2_;
+
+    // Mass of newly created CO [kg]
+    const scalar dmCO = dOmega*(WC_ + Sb_*WO2_);
+
+    // Update local particle C mass
+    dMassSolid[CsLocalId_] += dOmega*WC_;
+
+    // Update carrier O2 and CO mass
+    dMassSRCarrier[O2GlobalId_] -= dmO2;
+    dMassSRCarrier[COGlobalId_] += dmCO;
+
+    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
+
+    // carrier sensible enthalpy exchange handled via change in mass
+
+    // Heat of reaction [J]
+//Info<<"Heat of reaction(Char)"<<tab<<(dmC*HsC - dmCO*HcCO_)/dmC<<endl;
+    return dmC*HsC - dmCO*HcCO_;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    COxidationKD_CO
+
+Description
+    Kinetic/diffusion limited rate surface reaction model for coal parcels.
+    Limited to:
+
+        C(s) + Sb*O2 -> CO
+
+    where Sb is the stoichiometry of the reaction
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef COxidationKD_CO_H
+#define COxidationKD_CO_H
+
+#include "SurfaceReactionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+template<class CloudType>
+class COxidationKD_CO;
+
+/*---------------------------------------------------------------------------*\
+            Class COxidationKD_CO Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class COxidationKD_CO
+:
+    public SurfaceReactionModel<CloudType>
+{
+    // Private data
+
+        // Model constants
+
+            //- Stoichiometry of reaction
+            const scalar Sb_;
+
+            //- Mass diffusion limited rate constant, C1
+            const scalar C1_;
+
+            //- Kinetics limited rate pre-exponential constant, C2
+            const scalar C2_;
+
+            //- Kinetics limited rate activation energy
+            const scalar E_;
+
+
+        // Addressing
+
+            //- Cs positions in global/local lists
+            label CsLocalId_;
+
+            //- O2 position in global list
+            label O2GlobalId_;
+
+            //- CO2 positions in global list
+            label COGlobalId_;
+
+
+        // Local copies of thermo properties
+
+            //- Molecular weight of C [kg/kmol]
+            scalar WC_;
+
+            //- Molecular weight of O2 [kg/kmol]
+            scalar WO2_;
+
+            //- Formation enthalpy for CO2 [J/kg]
+            scalar HcCO_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("COxidationKD_CO");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        COxidationKD_CO
+        (
+            const dictionary& dict,
+            CloudType& owner
+        );
+
+        //- Construct copy
+        COxidationKD_CO
+        (
+            const COxidationKD_CO<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
+        {
+            return autoPtr<SurfaceReactionModel<CloudType> >
+            (
+                new COxidationKD_CO<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~COxidationKD_CO();
+
+
+    // Member Functions
+
+        //- Update surface reactions
+        virtual scalar calculate
+        (
+            const scalar dt,
+            const label cellI,
+            const scalar d,
+            const scalar T,
+            const scalar Tc,
+            const scalar pc,
+            const scalar rhoc,
+            const scalar mass,
+            const scalarField& YGas,
+            const scalarField& YLiquid,
+            const scalarField& YSolid,
+            const scalarField& YMixture,
+            const scalar N,
+            scalarField& dMassGas,
+            scalarField& dMassLiquid,
+            scalarField& dMassSolid,
+            scalarField& dMassSRCarrier
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "COxidationKD_CO.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/Make/files

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

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

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

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

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

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

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

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,17 +23,16 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "basicThermoKinematicMPPICParcel.H"
+#include "Cloud.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+namespace Foam
-
-inline const Foam::compressibleTurbulenceModel&
-Foam::radiation::sootModel::turbulence() const
 {
-    return turbulence_;
+    defineTemplateTypeNameAndDebug(basicThermoKinematicMPPICParcel, 0);
+    defineTemplateTypeNameAndDebug(Cloud<basicThermoKinematicMPPICParcel>, 0);
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 
 // ************************************************************************* //

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

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

src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C

@@ -86,4 +86,17 @@ Foam::basicChemistryModel::~basicChemistryModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::basicChemistryModel::edm
+(
+    const scalarField& mixing,
+    const scalar A,
+    const scalar B
+)
+{
+    NotImplemented;
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H

@@ -159,6 +159,14 @@ public:
                 //- Calculates the reaction rates
                 virtual void calculate() = 0;
 
+                //- Calculates the eddy dissipation reaction rates
+                virtual void edm
+                (
+                    const scalarField& mixing,
+                    const scalar A,
+                    const scalar B
+                );
+
                 //- Solve the reaction system for the given time step
                 //  and return the characteristic time
                 virtual scalar solve(const scalar deltaT) = 0;

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C

@@ -756,6 +756,104 @@ void Foam::chemistryModel<CompType, ThermoType>::calculate()
 }
 
 
+template<class CompType, class ThermoType>
+void Foam::chemistryModel<CompType, ThermoType>::edm
+(
+    const scalarField& mixing,
+    const scalar A,
+    const scalar B
+)
+{
+    if (!this->chemistry_)
+    {
+        return;
+    }
+
+    const volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            this->time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->thermo().rho()
+    );
+
+    forAll(rho, celli)
+    {
+        const scalar rhoi = rho[celli];
+        const scalar mixingi = mixing[celli];
+
+        scalarField Yi(nSpecie_, 0.0);
+        for (label k=0; k<nSpecie_; k++)
+        {
+            Yi[k] = Y_[k][celli];
+        }
+
+        scalarField om(nSpecie_, 0.0);
+
+        scalar minYR = 1.0;
+        scalar YP = 0.0;
+        forAll(reactions_, m)
+        {
+            const Reaction<ThermoType>& R = reactions_[m];
+
+            scalar num = 0.0;
+            scalar den = 0.0;
+
+            forAll(R.lhs(), s)
+            {
+                const label k = R.lhs()[s].index;
+                const scalar vk = R.lhs()[s].stoichCoeff + SMALL;
+                const scalar Wk = specieThermo_[k].W();
+
+                minYR = min(minYR, Yi[k]/Wk/vk);
+            }
+
+            forAll(R.rhs(), s)
+            {
+                const label k = R.rhs()[s].index;
+                const scalar vk = R.rhs()[s].stoichCoeff;
+                const scalar Wk = specieThermo_[k].W();
+
+                num += Yi[k];
+                den += vk*Wk;
+            }
+            YP = num/den;
+
+            scalar omegam = min (A*rhoi*mixingi*minYR, A*B*rhoi*mixingi*YP);
+
+            forAll(R.lhs(), s)
+            {
+                const label k = R.lhs()[s].index;
+                const scalar vk = R.lhs()[s].stoichCoeff;
+                const scalar Wk = specieThermo_[k].W();
+
+                om[k] -= vk*Wk*omegam;
+            }
+
+            forAll(R.rhs(), s)
+            {
+                const label k = R.rhs()[s].index;
+                const scalar vk = R.rhs()[s].stoichCoeff;
+                const scalar Wk = specieThermo_[k].W();
+
+                om[k] += vk*Wk*omegam;
+            }
+        }
+
+        for (label k=0; k<nSpecie_; k++)
+        {
+            RR_[k][celli] = om[k];
+        }
+    }
+}
+
+
 template<class CompType, class ThermoType>
 template<class DeltaTType>
 Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H

@@ -191,6 +191,14 @@ public:
         //- Calculates the reaction rates
         virtual void calculate();
 
+        //- Calculates the eddy dissipation reaction rates
+        virtual void edm
+        (
+            const scalarField& mixing,
+            const scalar A,
+            const scalar B
+        );
+
 
         // Chemistry model functions (overriding abstract functions in
         // basicChemistryModel.H)

src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModels.C

@@ -111,23 +111,6 @@ namespace Foam
         rhoChemistryModel,
         icoPoly8EThermoPhysics
     );
-
-    // Chemistry moldels based on absoluteEnthalpy
-    makeChemistryModel
-    (
-        chemistryModel,
-        rhoChemistryModel,
-        gasHaThermoPhysics
-    );
-
-    // Chemistry moldels based on absoluteInternalEnergy
-    makeChemistryModel
-    (
-        chemistryModel,
-        rhoChemistryModel,
-        gasEaThermoPhysics
-    );
-
 }
 
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolvers.C

@@ -88,12 +88,6 @@ namespace Foam
         incompressibleGasEThermoPhysics
     );
     makeChemistrySolverTypes(rhoChemistryModel, icoPoly8EThermoPhysics);
-
-    // Chemistry solvers based on absoluteEnthalpy
-    makeChemistrySolverTypes(rhoChemistryModel, gasHaThermoPhysics);
-
-    // Chemistry solvers based on absoluteInternalEnergy
-    makeChemistrySolverTypes(rhoChemistryModel, gasEaThermoPhysics);
 }
 
 

src/thermophysicalModels/radiation/Make/files

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

src/thermophysicalModels/radiation/Make/options

@@ -9,8 +9,6 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 

src/thermophysicalModels/radiation/radiationModels/P1/P1.C

@@ -28,7 +28,6 @@ License
 #include "fvmSup.H"
 #include "absorptionEmissionModel.H"
 #include "scatterModel.H"
-#include "sootModel.H"
 #include "constants.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -89,19 +88,6 @@ Foam::radiation::P1::P1(const volScalarField& T)
         mesh_,
         dimensionedScalar("a", dimless/dimLength, 0.0)
     ),
-    aSoot_
-    (
-        IOobject
-        (
-            "aSoot",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("aSoot", dimless/dimLength, 0.0)
-    ),
     e_
     (
         IOobject
@@ -172,19 +158,6 @@ Foam::radiation::P1::P1(const dictionary& dict, const volScalarField& T)
         mesh_,
         dimensionedScalar("a", dimless/dimLength, 0.0)
     ),
-    aSoot_
-    (
-        IOobject
-        (
-            "aSoot",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("aSoot", dimless/dimLength, 0.0)
-    ),
     e_
     (
         IOobject
@@ -239,26 +212,9 @@ bool Foam::radiation::P1::read()
 
 void Foam::radiation::P1::calculate()
 {
-    a_ = absorptionEmission_->a() + soot_->a();
+    a_ = absorptionEmission_->a();
-    aSoot_ = soot_->a();
+    e_ = absorptionEmission_->e();
-    e_ = absorptionEmission_->e() + soot_->e();
     E_ = absorptionEmission_->E();
-
-    scalar sootradiationChecker = 0.0;
-
-    forAll(aSoot_, celli)
-    {
-        if (aSoot_[celli] != 0.0)
-        {
-            sootradiationChecker = 1.0;
-        }
-    }
-
-    if (sootradiationChecker == 0.0)
-    {
-        Info<<"soot radiation effects are unsuccessfully considered!"<<endl;
-    }
-
     const volScalarField sigmaEff(scatter_->sigmaEff());
 
     const dimensionedScalar a0 ("a0", a_.dimensions(), ROOTVSMALL);
@@ -316,7 +272,7 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::P1::Rp() const
                 IOobject::NO_WRITE,
                 false
             ),
-            4.0*(absorptionEmission_->eCont()+soot_->e())*physicoChemical::sigma
+            4.0*absorptionEmission_->eCont()*physicoChemical::sigma
         )
     );
 }
@@ -330,7 +286,7 @@ Foam::radiation::P1::Ru() const
     const DimensionedField<scalar, volMesh> E =
         absorptionEmission_->ECont()()();
     const DimensionedField<scalar, volMesh> a =
-        absorptionEmission_->aCont()()()+soot_->a()()();
+        absorptionEmission_->aCont()()();
 
     return a*G - E;
 }

src/thermophysicalModels/radiation/radiationModels/P1/P1.H

@@ -70,8 +70,6 @@ class P1
         //- Absorption coefficient
         volScalarField a_;
 
-        volScalarField aSoot_;
-
         //- Emission coefficient
         volScalarField e_;
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C

@@ -167,46 +167,6 @@ void Foam::radiation::fvDOM::initialise()
 
 
     // Construct absorption field for each wavelength
-    forAll(agasLambda_, lambdaI)
-    {
-        agasLambda_.set
-        (
-            lambdaI,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "agasLambda_" + Foam::name(lambdaI) ,
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                agas_
-            )
-        );
-    }
-
-    forAll(aSootLambda_, lambdaI)
-    {
-        aSootLambda_.set
-        (
-            lambdaI,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "aSootLambda_" + Foam::name(lambdaI) ,
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                aSoot_
-            )
-        );
-    }
-
     forAll(aLambda_, lambdaI)
     {
         aLambda_.set
@@ -227,7 +187,6 @@ void Foam::radiation::fvDOM::initialise()
         );
     }
 
-
     Info<< "fvDOM : Allocated " << IRay_.size()
         << " rays with average orientation:" << nl;
 
@@ -327,32 +286,6 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         mesh_,
         dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
     ),
-    agas_
-    (
-        IOobject
-        (
-            "agas",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("agas", dimless/dimLength, 0.0)
-    ),
-    aSoot_
-    (
-        IOobject
-        (
-            "aSoot",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("aSoot", dimless/dimLength, 0.0)
-    ),
     a_
     (
         IOobject
@@ -370,8 +303,6 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     nPhi_(readLabel(coeffs_.lookup("nPhi"))),
     nRay_(0),
     nLambda_(absorptionEmission_->nBands()),
-    agasLambda_(nLambda_),
-    aSootLambda_(nLambda_),
     aLambda_(nLambda_),
     blackBody_(nLambda_, T),
     IRay_(0),
@@ -444,32 +375,6 @@ Foam::radiation::fvDOM::fvDOM
         mesh_,
         dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
     ),
-    agas_
-    (
-        IOobject
-        (
-            "agas",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("agas", dimless/dimLength, 0.0)
-    ),
-    aSoot_
-    (
-        IOobject
-        (
-            "aSoot",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("aSoot", dimless/dimLength, 0.0)
-    ),
     a_
     (
         IOobject
@@ -487,8 +392,6 @@ Foam::radiation::fvDOM::fvDOM
     nPhi_(readLabel(coeffs_.lookup("nPhi"))),
     nRay_(0),
     nLambda_(absorptionEmission_->nBands()),
-    agasLambda_(nLambda_),
-    aSootLambda_(nLambda_),
     aLambda_(nLambda_),
     blackBody_(nLambda_, T),
     IRay_(0),
@@ -529,17 +432,7 @@ bool Foam::radiation::fvDOM::read()
 
 void Foam::radiation::fvDOM::calculate()
 {
-    absorptionEmission_->correct(agas_, agasLambda_);
+    absorptionEmission_->correct(a_, aLambda_);
-
-    soot_->correct(aSoot_, aSootLambda_);
-
-//    aLambda_ = agasLambda_ + aSootLambda_;
-
-    forAll(aLambda_, celli)
-    {
-        a_[celli] = agas_[celli] + aSoot_[celli];
-        aLambda_[celli] = agasLambda_[celli] + aSootLambda_[celli];
-    }
 
     updateBlackBodyEmission();
 
@@ -590,7 +483,7 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
                 false
             ),
             // Only include continuous phase emission
-            4*(absorptionEmission_->aCont() + soot_->a())*physicoChemical::sigma
+            4*absorptionEmission_->aCont()*physicoChemical::sigma
         )
     );
 }
@@ -608,7 +501,7 @@ Foam::radiation::fvDOM::Ru() const
 
     // Only include continuous phase absorption
     const DimensionedField<scalar, volMesh> a =
-        absorptionEmission_->aCont()()() + soot_->a()()();
+        absorptionEmission_->aCont()()();
 
     return a*G - E;
 }

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H

@@ -101,10 +101,6 @@ class fvDOM
         volScalarField Qin_;
 
         //- Total absorption coefficient [1/m]
-        volScalarField agas_;
-
-        volScalarField aSoot_;
-
         volScalarField a_;
 
         //- Number of solid angles in theta
@@ -120,10 +116,6 @@ class fvDOM
         label nLambda_;
 
         //- Wavelength total absorption coefficient [1/m]
-        PtrList<volScalarField> agasLambda_;
-
-        PtrList<volScalarField> aSootLambda_;
-
         PtrList<volScalarField> aLambda_;
 
         //- Black body

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H

@@ -37,7 +37,6 @@ SourceFiles
 
 #include "absorptionEmissionModel.H"
 #include "blackBodyEmission.H"
-#include "sootModel.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

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

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

@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::radiation::wsggmAbsorptionEmission
+
+Description
+    wsggmAbsorptionEmission radiation absorption and emission coefficients
+    for continuous phase
+
+    Ref. T.F Smith, Z.F. Shen and J.N. Friedman
+         Evaluation of Coefficients for the Weighted Sum of Gray Gases Model,
+         Transactions of the ASME, Vol. 104, pp. 602-608
+
+    Valid when Ptotal = 1atm, 0.001 < P*S < 10.0atm.m, 600 < T < 2400K
+               (where P*S = sum of partial pressures of absorbing gases
+                            * mean beam length)
+
+    i.e. dictionary
+    \verbatim
+        wsggmAbsorptionEmissionCoeffs
+        {
+            pathLength pathLength [ 0 1  0  0  0  0  0] 0.251; //3.6V/S
+            // Pw/Pc = 2.0
+            emissivityCoeffs 3(0.4201 6.516 131.9);
+            fittingFactors
+            3
+            (
+                4(6.508e-1 -5.551e-4 3.029e-7  -5.353e-11)
+                4(-0.2504e-1  6.112e-4  -3.882e-7  6.528e-11)
+                4(2.718e-1 -3.118e-4  1.221e-7  -1.612e-11)
+            );
+        }
+    \endverbatim
+
+
+SourceFiles
+    wsggmAbsorptionEmission.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wsggmAbsorptionEmission_H
+#define wsggmAbsorptionEmission_H
+
+#include "interpolationLookUpTable.H"
+#include "absorptionEmissionModel.H"
+#include "HashTable.H"
+#include "absorptionCoeffs.H"
+#include "fluidThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class wsggmAbsorptionEmission Declaration
+\*---------------------------------------------------------------------------*/
+
+class wsggmAbsorptionEmission
+:
+    public absorptionEmissionModel
+{
+public:
+
+private:
+
+    // Private data
+
+    //- Absorption model dictionary
+    dictionary coeffsDict_;
+
+    //- SLG thermo package
+    const fluidThermo& thermo_;
+
+    //- absorption coefficitent for individual grey gases
+    scalarList emissivityCoeffs_;
+
+    //- fitting Factors
+    scalarListList fittingFactors_;
+
+    //characteristic path length [m] from dict
+    dimensionedScalar pathLength_;
+
+public:
+
+    //- Runtime type information
+    TypeName("wsggmAbsorptionEmission");
+
+
+    // Constructors
+
+        //- Construct from components
+        wsggmAbsorptionEmission(const dictionary& dict, const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~wsggmAbsorptionEmission();
+
+
+    // Member Functions
+
+        // Access
+
+            // Absorption coefficient
+
+                //- Absorption coefficient for continuous phase
+                tmp<volScalarField> aCont(const label bandI = 0) const;
+
+
+            // Emission coefficient
+
+                //- Emission coefficient for continuous phase
+                tmp<volScalarField> eCont(const label bandI = 0) const;
+
+
+            // Emission contribution
+
+                //- Emission contribution for continuous phase
+                tmp<volScalarField> ECont(const label bandI = 0) const;
+
+
+    // Member Functions
+
+        inline bool isGrey() const
+        {
+            return true;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/sootModel/MossBrookesSoot/MossBrookesSoot.C

@@ -1,487 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "MossBrookesSoot.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvm.H"
-#include "basicSpecieMixture.H"
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-namespace Foam
-{
-    namespace radiation
-    {
-        defineTypeNameAndDebug(MossBrookesSoot, 0);
-
-        addToRunTimeSelectionTable
-        (
-            sootModel,
-            MossBrookesSoot,
-            dictionary
-        );
-    }
-}
-
-
-const Foam::basicSpecieMixture&
-Foam::radiation::MossBrookesSoot::checkThermo
-(
-    const fluidThermo& thermo
-)
-{
-    if (isA<basicSpecieMixture>(thermo))
-    {
-        return dynamic_cast<const basicSpecieMixture& >
-        (
-            thermo
-        );
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Inconsistent thermo package for " << thermo.type()
-            << "Please select a thermo package based on "
-            << "basicSpecieMixture" << exit(FatalError);
-
-        return dynamic_cast<const basicSpecieMixture& >
-        (
-            thermo
-        );
-    }
-
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::radiation::MossBrookesSoot::MossBrookesSoot
-(
-    const dictionary& dict,
-    const fvMesh& mesh,
-    const word& modelType
-)
-:
-    sootModel(dict, mesh, modelType),
-    soot_
-    (
-        IOobject
-        (
-            "soot",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    sootVF_
-    (
-        IOobject
-        (
-            "sootVF",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    NucConc_
-    (
-        IOobject
-        (
-            "NucConc",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    NumDen_
-    (
-        IOobject
-        (
-            "NumDen",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    coeffsDict_(dict.subOrEmptyDict(modelType + "Coeffs")),
-    thermo_(mesh.lookupObject<fluidThermo>(basicThermo::dictName)),
-    mixture_(checkThermo(thermo_)),
-    turbulence_(mesh.lookupObject<compressibleTurbulenceModel>(turbulenceModel::propertiesName)),
-    Snet_
-    (
-        IOobject
-        (
-            "Snet",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Snet", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Snuc_
-    (
-        IOobject
-        (
-            "Snuc",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Snuc", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Smassgrow_
-    (
-        IOobject
-        (
-            "Smassgrow",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Smassgrow", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Soxid_
-    (
-        IOobject
-        (
-            "Soxid",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Soxid", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    NSnuc_
-    (
-        IOobject
-        (
-            "NSnuc",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("NSnuc", dimless/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    NSagg_
-    (
-        IOobject
-        (
-            "NSagg",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("NSagg", dimless/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    NSnet_
-    (
-        IOobject
-        (
-            "NSnet",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("NSnet", dimless/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    aggParam_
-    (
-        IOobject
-        (
-            "ParameterAgg",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("ParameterAgg", dimless, 0.0)
-    ),
-    Sf_
-    (
-        IOobject
-        (
-            "Sf",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Sf", dimless, 0.0)
-    ),
-    sootMeanDiameter_
-    (
-        IOobject
-        (
-            "SootMeanDiameter",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("SootMeanDiameter", dimLength, 0.0)
-    ),
-    PrtSoot(readScalar(coeffsDict_.lookup("TurbulentPrandtlNumberSoot"))),
-    rhoSoot(readScalar(coeffsDict_.lookup("SootDensity"))),
-    Calpha(readScalar(coeffsDict_.lookup("Calpha"))),
-    Mp(readScalar(coeffsDict_.lookup("Mp"))),
-    Cgamma(readScalar(coeffsDict_.lookup("Cgamma"))),
-    Coxid(readScalar(coeffsDict_.lookup("Coxid"))),
-    Cw(readScalar(coeffsDict_.lookup("Cw"))),
-    etaColl(readScalar(coeffsDict_.lookup("etaColl"))),
-    Talpha(readScalar(coeffsDict_.lookup("Talpha"))),
-    Tgamma(readScalar(coeffsDict_.lookup("Tgamma"))),
-    Navog(readScalar(coeffsDict_.lookup("AvogadroNumber"))),
-    Cbeta(readScalar(coeffsDict_.lookup("Cbeta"))),
-    Boltzmann(readScalar(coeffsDict_.lookup("BoltzmannConst"))),
-    radiationActive_(readBool(coeffsDict_.lookup("Radiationeffect"))),
-    cliqueNuc(readScalar(coeffsDict_.lookup("NucExpCorrectConstant"))),
-    cliqueMassGrow(readScalar(coeffsDict_.lookup("MassGrowCorrectConstant")))
-{
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::radiation::MossBrookesSoot::~MossBrookesSoot()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::MossBrookesSoot::a(const label bandI) const
-{
-
-    const volScalarField T_ = thermo_.T();
-
-    tmp<volScalarField> tas
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "aSoot" + name(bandI),
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("as", dimless/dimLength, 0.0),
-            extrapolatedCalculatedFvPatchVectorField::typeName
-        )
-    );
-
-    scalarField& as = tas.ref().primitiveFieldRef();
-
-    forAll(as, celli)
-    {
-        if (radiationActive_ == true)
-        {
-            as[celli] = (1232.4)*(1800.0)*soot_[celli]*(1.0+(4.8e-4)*(T_[celli]-2000.0));
-        }
-        else
-        {
-            as[celli] = 0.0;
-        }
-    }
-    tas.ref().correctBoundaryConditions();
-    return tas;
-}
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::MossBrookesSoot::e(const label bandI) const
-{
-
-    return a(bandI);
-}
-
-void Foam::radiation::MossBrookesSoot::correct()
-{
-
-    calcSource();
-
-    const volScalarField rho_ = thermo_.rho();
-
-    const volScalarField &nut_ = turbulence().nut();
-
-    const surfaceScalarField phi_ (turbulence().phi());
-
-    fvScalarMatrix SootEqn
-    (
-        fvm::div(phi_, soot_)
-      - fvm::laplacian((rho_*nut_/PrtSoot), soot_)
-     ==
-        Snet_
-    );
-
-    SootEqn.relax();
-
-    SootEqn.solve();
-
-    fvScalarMatrix NucConcEqn
-    (
-        fvm::div(phi_, NucConc_)
-      - fvm::laplacian((rho_*nut_/PrtSoot), NucConc_)
-     ==
-        NSnet_
-    );
-
-    NucConcEqn.relax();
-
-    NucConcEqn.solve();
-
-    scalar minusindicator1 = 0.0;
-    scalar minusindicator2 = 0.0;
-
-    forAll(rho_, celli)
-    {
-        if (soot_[celli] < 0.0) {
-            soot_[celli] = 0.0;
-            minusindicator1 = 1.0;
-        }
-        if (NucConc_[celli] < 0.0) {
-            NucConc_[celli] =0.0;
-            minusindicator2 = 1.0;
-        }
-
-        NumDen_[celli] = rho_[celli]*(NucConc_[celli]*1e+15);
-
-        sootVF_[celli] = soot_[celli]*rho_[celli]/rhoSoot;
-
-    }
-    if (minusindicator1 == 1.0 and minusindicator2 == 0.0) {
-        Info<< "Minus soot value! corrected to zero" <<endl;
-    }
-    else if (minusindicator1 == 0.0 and minusindicator2 ==1.0) {
-        Info<< "Minus NucConc value! corrected to zero" <<endl;
-    }
-    else if (minusindicator1 == 1.0 and minusindicator2 == 1.0) {
-        Info<< "Minus soot and NucConc value! corrected to zero" << endl;
-    }
-}
-
-
-void Foam::radiation::MossBrookesSoot::calcSource()
-{
-    const volScalarField rho_ = thermo_.rho();
-    const volScalarField T_ = thermo_.T();
-    const volScalarField p_ = thermo_.p();
-
-    const scalar precI = mixture_.species()[coeffsDict_.lookup("Precursor")];
-    const scalar growI = mixture_.species()[coeffsDict_.lookup("SurfaceGrowSpecie")];
-    const scalar oxidI = mixture_.species()[coeffsDict_.lookup("SootOxidizerSpecie")];
-
-    const volScalarField &Yprec_ = mixture_.Y(precI);
-    const volScalarField &Ygrow_ = mixture_.Y(growI);
-    const volScalarField &Yoxid_ = mixture_.Y(oxidI);
-
-    const volScalarField XPrec_(mixture_.W()*Yprec_/mixture_.W(precI));
-    const volScalarField XSgs_ (mixture_.W()*Ygrow_/mixture_.W(growI));
-    const volScalarField XOxid_(mixture_.W()*Yoxid_/mixture_.W(oxidI));
-
-    forAll(rho_, celli)
-    {
-        const scalar XPreci = XPrec_[celli];
-        const scalar XSgsi = XSgs_[celli];
-        const scalar XOxidi = XOxid_[celli];
-        const scalar Ti = T_[celli];
-        const scalar pi = p_[celli];
-        const scalar rhoi = rho_[celli];
-        const scalar sooti = soot_[celli];
-        const scalar NumDeni = NumDen_[celli];
-
-        const scalar MolarConcPreci = (XPreci*pi)/(8314.5*Ti);
-        const scalar MolarConcSgsi = (XSgsi*pi)/(8314.5*Ti);
-        const scalar MolarConcOxidi = (XOxidi*pi)/(8314.5*Ti);
-
-        scalar SootMeanDiameteri = 0.0;
-
-        if (NumDeni == 0.0) {
-            sootMeanDiameter_[celli] = 0.0;
-            Sf_[celli] = 0.0;
-        }
-        else {
-            const scalar parameter = (6.0/3.1416)*(rhoi/rhoSoot)*(sooti/NumDeni);
-            sootMeanDiameter_[celli] = pow(parameter,0.33333);
-            SootMeanDiameteri = sootMeanDiameter_[celli];
-            Sf_[celli] = (3.1416)*pow(SootMeanDiameteri,2.0)*NumDeni;
-        }
-
-        Snuc_[celli] = Mp*Calpha*exp(-cliqueNuc*Talpha/Ti)*MolarConcPreci;
-
-        Smassgrow_[celli] = Cgamma*MolarConcSgsi*exp(-cliqueMassGrow*Tgamma/Ti)*Sf_[celli];
-
-        Soxid_[celli] = Coxid*Cw*etaColl*MolarConcOxidi*pow(Ti,0.5)*Sf_[celli];
-
-        Snet_[celli] = Snuc_[celli] + Smassgrow_[celli] - Soxid_[celli];
-
-        NSnuc_[celli] = (Calpha*Navog*MolarConcPreci*exp(-cliqueNuc*Talpha/Ti))*(1e-15);
-
-        aggParam_[celli] = 24.0*Boltzmann*Ti/rhoSoot;
-
-        NSagg_[celli] = (Cbeta*pow((aggParam_[celli]),0.5)*pow(SootMeanDiameteri,0.5)*pow(NumDeni,2.0))*(1e-15);
-
-        NSnet_[celli] = NSnuc_[celli] - NSagg_[celli];
-    }
-}
-
-
-const Foam::compressibleTurbulenceModel&
-Foam::radiation::MossBrookesSoot::turbulence() const
-{
-    return turbulence_;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiation/submodels/sootModel/MossBrookesSoot/MossBrookesSoot.H

@@ -1,283 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::radiation::MossBrookesSoot
-
-Description
-    This soot model solves two transport equation for soot mass fraction
-    and soot normalized nuclei concentration.
-    The detail of the model strictly follows Moss & Brookes model in fluent.
-
-    Parameter settings are available in constant/radiationProperties.
-
-    Setting example:
-
-    sootModel MossBrookesSoot<rhoChemistryModel, gasHaThermoPhysics>;
-
-    MossBrookesSootCoeffs
-    {
-        SootDensity            1800.0;
-        Calpha                 54.0;
-        Mp                     144;
-        Cgamma                 11700;
-        Coxid                  1.0;
-        Cw                     105.8125;
-        etaColl                0.04;
-        Talpha                 21100.0;
-        Tgamma                 12100.0;
-        AvogadroNumber         6.022e+26;
-        Cbeta                  1.0;
-        BoltzmannConst         1.38e-23;
-        Radiationeffect        true;
-
-        Precursor              C2H2;
-        SurfaceGrowSpecie      C2H2;
-        SootOxidizerSpecie     OH;
-    }
-
-SourceFiles
-    MossBrookesSoot.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef MossBrookesSoot_H
-#define MossBrookesSoot_H
-
-#include "sootModel.H"
-#include "HashTable.H"
-
-#include "fluidThermo.H"
-#include "basicSpecieMixture.H"
-#include "turbulentFluidThermoModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-namespace radiation
-{
-
-/*---------------------------------------------------------------------------*\
-                    Class MossBrookesSoot Declaration
-\*---------------------------------------------------------------------------*/
-class MossBrookesSoot
-:
-    public sootModel
-
-{
-
-    // Static functions
-
-        //- Check mixture in thermo
-        static const basicSpecieMixture& checkThermo
-        (
-            const fluidThermo&
-        );
-
-
-    // Private data
-
-        //- Soot mass fraction
-        volScalarField soot_;
-
-        //- Soot volume fraction
-        volScalarField sootVF_;
-
-        //- Normalized nuclei concentration of soot
-        volScalarField NucConc_;
-
-        //- Number density of soot
-        volScalarField NumDen_;
-
-        //- Soot model dictionary
-        dictionary coeffsDict_;
-
-        //- Thermo package
-        const fluidThermo& thermo_;
-
-        //- Composition
-        const basicSpecieMixture& mixture_;
-
-        //- Reference to the turbulence model
-        const compressibleTurbulenceModel &turbulence_;
-
-        //- Net rate of soot generation
-        volScalarField Snet_;
-
-        //- The rate of soot formation by nucleation
-        volScalarField Snuc_;
-
-        //- The rate of soot formation by mass growth
-        volScalarField Smassgrow_;
-
-        //- The rate of soot formation by soot oxidation
-        volScalarField Soxid_;
-
-        //- The rate of soot number density formation by nucleation
-        volScalarField NSnuc_;
-
-        //- The rate of soot number density formation by agglomeration
-        volScalarField NSagg_;
-
-        //- Net rate of soot numbe density formation
-        volScalarField NSnet_;
-
-        //- Parameter for Agglomeration source term
-        volScalarField aggParam_;
-
-        //- Available surface area of soot oxidation per volume
-        volScalarField Sf_;
-
-        //- Soot particle mean diameter
-        volScalarField sootMeanDiameter_;
-
-        //- Soot turbulent Prandtl number
-        scalar PrtSoot;
-
-        //- Density of soot particle
-        scalar rhoSoot;
-
-        //- Model constant for soot inception rate
-        scalar Calpha;
-
-        //- Mass of an incipient soot particle
-        scalar Mp;
-
-        //- Surface growth rate scaling factor
-        scalar Cgamma;
-
-        //- Oxidation rate scaling parameter
-        scalar Coxid;
-
-        //- Oxidation model constant
-        scalar Cw;
-
-        //- Collisional efficiency parameter
-        scalar etaColl;
-
-        //- Activation temperature for soot inception
-        scalar Talpha;
-
-        //- Activation temperature for soot growth rate
-        scalar Tgamma;
-
-        //- Avogadro's number
-        scalar Navog;
-
-        //- Model constant for coagulation rate
-        scalar Cbeta;
-
-        //- Boltzmann constant
-        scalar Boltzmann;
-
-        //- Soot radiative effect on/off flag
-        bool radiationActive_;
-
-        //- Correction constant for Exponent in Nucleation
-        scalar cliqueNuc;
-
-        //- Correction constant for Exponent in Soot Mass Grow
-        scalar cliqueMassGrow;
-
-        //- Soot precursor specie name
-        const word precSpecie;
-
-        //- Soot surface growth specie name
-        const word sgsSpecie;
-
-        //- Soot oxidation specie name
-        const word oxSpecie;
-
-protected:
-
-    // Protected data
-
-    // Protected Member Functions
-
-public:
-
-    //- Runtime type information
-    TypeName("MossBrookesSoot");
-
-
-    // Constructors
-
-        //- Construct from components
-        MossBrookesSoot
-        (
-            const dictionary& dict,
-            const fvMesh& mesh,
-            const word& modelType
-        );
-
-
-    //- Destructor
-    virtual ~MossBrookesSoot();
-
-
-        // Edit
-
-            //- Main update/correction routine
-            virtual void correct();
-
-            //- Calculation of the source term routine
-            virtual void calcSource();
-
-
-        // Access
-
-            //- Return Ysoot
-            const volScalarField& soot() const
-            {
-                return soot_;
-            }
-
-            //- Return access to turbulence
-            const compressibleTurbulenceModel& turbulence() const;
-
-            // Soot absorption coefficient
-
-                //- Soot absorption coefficient (net)
-                virtual tmp<volScalarField> a(const label bandI = 0) const;
-
-            // Soot emission coefficient
-
-                //- Soot emission coefficient (net)
-                virtual tmp<volScalarField> e(const label bandI = 0) const;
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-} // End namespace radiation
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/sootModel/khanGreeveSoot/khanGreeveSoot.C

@@ -1,374 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "khanGreeveSoot.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvm.H"
-#include "basicSpecieMixture.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-namespace Foam
-{
-    namespace radiation
-    {
-        defineTypeNameAndDebug(khanGreeveSoot, 0);
-
-        addToRunTimeSelectionTable
-        (
-            sootModel,
-            khanGreeveSoot,
-            dictionary
-        );
-    }
-}
-
-
-const Foam::basicSpecieMixture&
-Foam::radiation::khanGreeveSoot::checkThermo
-(
-    const fluidThermo& thermo
-)
-{
-    if (isA<basicSpecieMixture>(thermo))
-    {
-        return dynamic_cast<const basicSpecieMixture& >
-        (
-            thermo
-        );
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Inconsistent thermo package for " << thermo.type()
-            << "Please select a thermo package based on "
-            << "basicSpecieMixture" << exit(FatalError);
-
-        return dynamic_cast<const basicSpecieMixture& >
-        (
-            thermo
-        );
-    }
-
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::radiation::khanGreeveSoot::khanGreeveSoot
-
-(
-    const dictionary& dict,
-    const fvMesh& mesh,
-    const word& modelType
-)
-:
-    sootModel(dict, mesh, modelType),
-    soot_
-    (
-        IOobject
-        (
-            "soot",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    sootVF_
-    (
-        IOobject
-        (
-            "sootVF",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    coeffsDict_(dict.subOrEmptyDict(modelType + "Coeffs")),
-    thermo_(mesh.lookupObject<fluidThermo>(basicThermo::dictName)),
-    mixture_(checkThermo(thermo_)),
-    turbulence_(mesh.lookupObject<compressibleTurbulenceModel>(turbulenceModel::propertiesName)),
-    Snet_
-    (
-        IOobject
-        (
-            "Snet",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Snet", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Sform_
-    (
-        IOobject
-        (
-            "Sform",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Sform", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Scomb_
-    (
-        IOobject
-        (
-            "Scomb",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Scomb", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    S1_
-    (
-        IOobject
-        (
-            "S1",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("S1", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    S2_
-    (
-        IOobject
-        (
-            "S2",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("S2", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    PrtSoot(readScalar(coeffsDict_.lookup("TurbulentPrandtlNumberSoot"))),
-    rhoSoot(readScalar(coeffsDict_.lookup("SootDensity"))),
-    Cs_(readScalar(coeffsDict_.lookup("Cs"))),
-    r_(readScalar(coeffsDict_.lookup("r"))),
-    Ta_(readScalar(coeffsDict_.lookup("ActivationTemperature"))),
-    A_(readScalar(coeffsDict_.lookup("A"))),
-    nuFuel_(readScalar(coeffsDict_.lookup("nufuel"))),
-    nuSoot_(readScalar(coeffsDict_.lookup("nusoot"))),
-    eqv_
-    (
-        IOobject
-        (
-            "eqv",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("eqv", dimless, 0.0)
-    ),
-    eqvMax(readScalar(coeffsDict_.lookup("EquivalenceRatioMax"))),
-    eqvMin(readScalar(coeffsDict_.lookup("EquivalenceRatioMin"))),
-    radiationActive_(readBool(coeffsDict_.lookup("Radiationeffect")))
-{
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::radiation::khanGreeveSoot::~khanGreeveSoot()
-
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::khanGreeveSoot::a(const label bandI) const
-{
-
-    const volScalarField T_ = thermo_.T();
-
-    tmp<volScalarField> tas
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "aSoot" + name(bandI),
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("aSoot", dimless/dimLength, 0.0),
-            extrapolatedCalculatedFvPatchVectorField::typeName
-        )
-    );
-
-    scalarField& as = tas.ref().primitiveFieldRef();
-    if (radiationActive_)
-    {
-        forAll(as, celli)
-        {
-            as[celli] = (1232.4)*(1800.0)*soot_[celli]*(1.0+(4.8e-4)*(T_[celli]-2000.0));
-        }
-    }
-    tas.ref().correctBoundaryConditions();
-
-    return tas;
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::khanGreeveSoot::e(const label bandI) const
-{
-
-    return a(bandI);
-}
-
-
-void Foam::radiation::khanGreeveSoot::correct()
-{
-
-    calcSource();
-
-    const volScalarField &rho_ = thermo_.rho();
-
-    const volScalarField &nut_ = turbulence().nut();
-
-    const surfaceScalarField phi_ (turbulence().phi());
-
-    fvScalarMatrix SootEqn
-    (
-        fvm::ddt(rho_, soot_)
-      + fvm::div(phi_, soot_)
-      - fvm::laplacian((rho_*nut_/PrtSoot), soot_)
-     ==
-        Snet_
-    );
-
-    SootEqn.relax();
-
-    SootEqn.solve();
-
-    forAll(rho_, celli)
-    {
-        if (soot_[celli] < 1e-30){
-            soot_[celli] = 0.0;
-        }
-
-        sootVF_[celli] = soot_[celli]*rho_[celli]/rhoSoot;
-
-    }
-
-}
-
-
-void Foam::radiation::khanGreeveSoot::calcSource()
-{
-    const volScalarField rho_ = thermo_.rho();
-    const volScalarField T_ = thermo_.T();
-    const volScalarField p_ = thermo_.p();
-
-    const scalar oxI = mixture_.species()[coeffsDict_.lookup("Oxidizer")];
-
-    const scalar fuelI = mixture_.species()[coeffsDict_.lookup("Fuel")];
-
-    const volScalarField &Yox_ = mixture_.Y(oxI);
-
-    const volScalarField &Yfuel_ = mixture_.Y(fuelI);
-
-    const volScalarField Xfuel_(mixture_.W()*Yfuel_/mixture_.W(fuelI));
-
-    const volScalarField &epsilon_ = turbulence().epsilon();
-
-    const volScalarField &k_ = turbulence().k();
-
-    forAll(rho_, celli)
-    {
-        const scalar Yoxi = Yox_[celli];
-        const scalar Yfueli = Yfuel_[celli];
-
-        eqv_[celli] = (Yfueli/Yoxi)*nuFuel_;
-
-        if (eqv_[celli] > eqvMax){
-            eqv_[celli] = 0.0;
-        }
-        else if (eqv_[celli] < eqvMin){
-            eqv_[celli] = 0.0;
-        }
-
-        const scalar pi = p_[celli];
-        const scalar Xfueli = Xfuel_[celli];
-        const scalar Ti = T_[celli];
-
-        Sform_[celli] = Cs_*pi*Xfueli*pow(eqv_[celli],r_)*exp(-Ta_/Ti);
-
-        const scalar rhoi = rho_[celli];
-        const scalar sooti = soot_[celli];
-        const scalar epsiloni = epsilon_[celli];
-        const scalar ki = k_[celli];
-
-        S1_[celli] = A_*rhoi*sooti*epsiloni/ki;
-
-        if (sooti == 0.0 and Yfueli == 0.0 ){
-            S2_[celli] = 0.0;
-        }
-        else {
-            S2_[celli] = A_*rhoi*(Yoxi/nuSoot_)*(sooti*nuSoot_/(sooti*nuSoot_+(Yfueli)*nuFuel_))*epsiloni/ki;
-
-        }
-
-        Scomb_[celli] = min(S1_[celli], S2_[celli]);
-
-        Snet_[celli] = Sform_[celli] - Scomb_[celli];
-    }
-
-}
-
-
-const Foam::compressibleTurbulenceModel&
-Foam::radiation::khanGreeveSoot::turbulence() const
-{
-    return turbulence_;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiation/submodels/sootModel/khanGreeveSoot/khanGreeveSoot.H

@@ -1,230 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::radiation::khanGreeveSoot
-
-Description
-    This soot model solves transport equation for soot mass fraction.
-    The detail of the model can be describe in the fluent theory guide.
-
-    Parameter setup is available in constant/radiationProperties file.
-
-    Setup example:
-
-    sootModel khanGreeveSoot<rhoChemistryModel,gasHaThermoPhysics>;
-
-    khanGreeveSootCoeffs
-    {
-        SootDensity                1800.0;
-        Cs                         1.5;
-        r                          3;
-        ActivationTemperature      20000;
-        A                          4;
-        nufuel                     3.2;
-        nusoot                     2.6667;
-        EquivalenceRatioMax        3.0;
-        EquivalenceRatioMin        1.67;
-        Radiationeffect            true;
-
-        Fuel                       PHC3H7;
-        Oxidizer                   O2;
-    }
-
-SourceFiles
-    khanGreeveSoot.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef khanGreeveSoot_H
-#define khanGreeveSoot_H
-
-#include "sootModel.H"
-#include "HashTable.H"
-
-#include "fluidThermo.H"
-#include "basicSpecieMixture.H"
-#include "turbulentFluidThermoModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-namespace radiation
-{
-
-/*---------------------------------------------------------------------------*\
-                    Class khanGreeveSoot Declaration
-\*---------------------------------------------------------------------------*/
-class khanGreeveSoot
-:
-    public sootModel
-
-{
-
-    // Static functions
-
-        //- Check mixture in thermo
-        static const basicSpecieMixture& checkThermo
-        (
-            const fluidThermo&
-        );
-
-
-    // Private data
-
-         //- Soot mass fraction
-        volScalarField soot_;
-
-        //- Soot volume fraction
-        volScalarField sootVF_;
-
-        //- Soot model dictionary
-        dictionary coeffsDict_;
-
-        //- Thermo package
-        const fluidThermo& thermo_;
-
-        //- Composition
-        const basicSpecieMixture& mixture_;
-
-        //- Reference to the turbulence model
-        const compressibleTurbulenceModel &turbulence_;
-
-        //- Net rate of soot generation
-        volScalarField Snet_;
-
-        //- The rate of soot formation
-        volScalarField Sform_;
-
-        //- The rate of soot combustion
-        volScalarField Scomb_;
-
-        //- The 1st rate of soot combustion
-        volScalarField S1_;
-
-        //- The 2nd rate of soot combustion
-        volScalarField S2_;
-
-        //- Soot turbulent Prandtl number
-        scalar PrtSoot;
-
-        //- Soot density
-        scalar rhoSoot;
-
-        //- Soot formation constant
-        scalar Cs_;
-
-        //- Equivalence ratio exponent
-        scalar r_;
-
-        //- Activation temperature
-        scalar Ta_;
-
-        //- Magnussen model constant
-        scalar A_;
-
-        //- Mass stoichiometry for the fuel oxidation
-        scalar nuFuel_;
-
-        //- Mass stoichiometry for the soot oxidation
-        scalar nuSoot_;
-
-        //- Equivalence ratio for the fuel
-        volScalarField eqv_;
-
-        //- Maximum equivalence ratio
-        scalar eqvMax;
-
-        //- Minimum equivalence ratio
-        scalar eqvMin;
-
-        //- Soot radiative effect on/off flag
-        bool radiationActive_;
-
-public:
-
-    //- Runtime type information
-    TypeName("khanGreeveSoot");
-
-
-    // Constructors
-
-        //- Construct from components
-        khanGreeveSoot
-        (
-            const dictionary& dict,
-            const fvMesh& mesh,
-            const word& modelType
-        );
-
-
-    //- Destructor
-    virtual ~khanGreeveSoot();
-
-
-        // Edit
-
-            //- Main update/correction routine
-            virtual void correct();
-
-            //- Calculation of the source term routine
-            virtual void calcSource();
-
-
-        // Access
-
-            //- Return Ysoot
-            const volScalarField& soot() const
-            {
-                return soot_;
-            }
-
-            //- Return access to turbulence
-            const compressibleTurbulenceModel& turbulence() const;
-
-            // Soot absorption coefficient
-
-                //- Soot absorption coefficient (net)
-                virtual tmp<volScalarField> a(const label bandI = 0) const;
-
-            // Soot emission coefficient
-
-                //- Soot emission coefficient (net)
-                virtual tmp<volScalarField> e(const label bandI = 0) const;
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-} // End namespace radiation
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/sootModel/sootModel/makeSootTypes.H

@@ -31,33 +31,6 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#define makeSootTypesCombThermo(sootModelType, Comp, Thermo)                   \
-    typedef Foam::radiation::sootModelType                                     \
-        <Foam::Comp, Foam::Thermo>                                             \
-        sootModelType##Comp##Thermo;                                           \
-                                                                               \
-    defineTemplateTypeNameAndDebugWithName                                     \
-    (                                                                          \
-        sootModelType##Comp##Thermo,                                           \
-        #sootModelType"<"#Comp","#Thermo">",                                   \
-        0                                                                      \
-    );                                                                         \
-                                                                               \
-    namespace Foam                                                             \
-    {                                                                          \
-        namespace radiation                                                    \
-        {                                                                      \
-            typedef sootModelType<Comp, Thermo> sootModelType##Comp##Thermo;   \
-                                                                               \
-            addToRunTimeSelectionTable                                         \
-            (                                                                  \
-                sootModel,                                                     \
-                sootModelType##Comp##Thermo,                                   \
-                dictionary                                                     \
-            );                                                                 \
-        }                                                                      \
-    }
-
 #define makeSootTypesThermo(sootModelType, Thermo)                             \
     typedef Foam::radiation::sootModelType<Foam::Thermo>                       \
         sootModelType##Thermo;                                                 \

src/thermophysicalModels/radiation/submodels/sootModel/sootModel/sootModel.C

@@ -58,58 +58,5 @@ Foam::radiation::sootModel::~sootModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::sootModel::a(const label bandI) const
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "aSoot",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimless/dimLength, 0.0)
-        )
-    );
-}
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::sootModel::e(const label bandI) const
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "eSoot",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimless/dimLength, 0.0)
-        )
-    );
-}
-
-void Foam::radiation::sootModel::correct
-(
-    volScalarField& a,
-    PtrList<volScalarField>& aj
-) const
-{
-    a = this->a();
-    aj[0] = a;
-}
 
 // ************************************************************************* //

src/thermophysicalModels/radiation/submodels/sootModel/sootModel/sootModel.H

@@ -138,26 +138,6 @@ public:
 
                     //- Return const reference to soot
                     virtual const volScalarField& soot() const = 0;
-
-                // Soot absorption coefficient
-
-                    //- Soot absorption coefficient (net)
-                    virtual tmp<volScalarField> a(const label bandI = 0) const;
-
-                // Soot emission coefficient
-
-                    //- Soot emission coefficient (net)
-                    virtual tmp<volScalarField> e(const label bandI = 0) const;
-
-                // Correct absorption coefficients
-
-                    virtual void correct
-                    (
-                        volScalarField& a,
-                        PtrList<volScalarField>& aj
-                    ) const;
-
-
 };
 
 

src/thermophysicalModels/radiationSootModel/Make/files

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

src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.C

@@ -1,485 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "MossBrookesSoot.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvm.H"
-#include "constants.H"
-#include "basicSpecieMixture.H"
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-namespace Foam
-{
-    namespace radiation
-    {
-        defineTypeNameAndDebug(MossBrookesSoot, 0);
-
-        addToRunTimeSelectionTable
-        (
-            sootModel,
-            MossBrookesSoot,
-            dictionary
-        );
-    }
-}
-
-
-const Foam::basicSpecieMixture&
-Foam::radiation::MossBrookesSoot::checkThermo
-(
-    const fluidThermo& thermo
-)
-{
-    if (isA<basicSpecieMixture>(thermo))
-    {
-        return dynamic_cast<const basicSpecieMixture& >
-        (
-            thermo
-        );
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Inconsistent thermo package for " << thermo.type()
-            << "Please select a thermo package based on "
-            << "basicSpecieMixture" << exit(FatalError);
-
-        return dynamic_cast<const basicSpecieMixture& >
-        (
-            thermo
-        );
-    }
-
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::radiation::MossBrookesSoot::MossBrookesSoot
-(
-    const dictionary& dict,
-    const fvMesh& mesh,
-    const word& modelType
-)
-:
-    sootModel(dict, mesh, modelType),
-    soot_
-    (
-        IOobject
-        (
-            "soot",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    sootVF_
-    (
-        IOobject
-        (
-            "sootVF",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    NucConc_
-    (
-        IOobject
-        (
-            "NucConc",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    NumDen_
-    (
-        IOobject
-        (
-            "NumDen",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    coeffsDict_(dict.subOrEmptyDict(modelType + "Coeffs")),
-    thermo_(mesh.lookupObject<fluidThermo>(basicThermo::dictName)),
-    mixture_(checkThermo(thermo_)),
-    turbulence_(mesh.lookupObject<compressibleTurbulenceModel>(turbulenceModel::propertiesName)),
-    Snet_
-    (
-        IOobject
-        (
-            "Snet",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Snet", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Snuc_
-    (
-        IOobject
-        (
-            "Snuc",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Snuc", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Smassgrow_
-    (
-        IOobject
-        (
-            "Smassgrow",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Smassgrow", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Soxid_
-    (
-        IOobject
-        (
-            "Soxid",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Soxid", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    NSnuc_
-    (
-        IOobject
-        (
-            "NSnuc",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("NSnuc", dimless/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    NSagg_
-    (
-        IOobject
-        (
-            "NSagg",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("NSagg", dimless/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    NSnet_
-    (
-        IOobject
-        (
-            "NSnet",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("NSnet", dimless/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    aggParam_
-    (
-        IOobject
-        (
-            "ParameterAgg",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("ParameterAgg", dimless, 0.0)
-    ),
-    Sf_
-    (
-        IOobject
-        (
-            "Sf",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Sf", dimless, 0.0)
-    ),
-    sootMeanDiameter_
-    (
-        IOobject
-        (
-            "SootMeanDiameter",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("SootMeanDiameter", dimLength, 0.0)
-    ),
-    Nnorm("Nnorm", dimless, 1e15),
-    PrtSoot(readScalar(coeffsDict_.lookup("TurbulentPrandtlNumberSoot"))),
-    rhoSoot(readScalar(coeffsDict_.lookup("SootDensity"))),
-    Calpha(readScalar(coeffsDict_.lookup("Calpha"))),
-    Mp(readScalar(coeffsDict_.lookup("Mp"))),
-    Cgamma(readScalar(coeffsDict_.lookup("Cgamma"))),
-    Coxid(readScalar(coeffsDict_.lookup("Coxid"))),
-    Cw(readScalar(coeffsDict_.lookup("Cw"))),
-    etaColl(readScalar(coeffsDict_.lookup("etaColl"))),
-    Talpha(readScalar(coeffsDict_.lookup("Talpha"))),
-    Tgamma(readScalar(coeffsDict_.lookup("Tgamma"))),
-    Cbeta(readScalar(coeffsDict_.lookup("Cbeta"))),
-    radiationActive_(readBool(coeffsDict_.lookup("Radiationeffect"))),
-    cliqueNuc(readScalar(coeffsDict_.lookup("NucExpCorrectConstant"))),
-    cliqueMassGrow(readScalar(coeffsDict_.lookup("MassGrowCorrectConstant")))
-{
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::radiation::MossBrookesSoot::~MossBrookesSoot()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::MossBrookesSoot::a(const label bandI) const
-{
-
-    const volScalarField T_ = thermo_.T();
-
-    tmp<volScalarField> tas
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "aSoot" + name(bandI),
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("as", dimless/dimLength, 0.0),
-            extrapolatedCalculatedFvPatchVectorField::typeName
-        )
-    );
-
-    scalarField& as = tas.ref().primitiveFieldRef();
-
-    forAll(as, celli)
-    {
-        if (radiationActive_ == true)
-        {
-            as[celli] = (1232.4)*(1800.0)*soot_[celli]*(1.0+(4.8e-4)*(T_[celli]-2000.0));
-        }
-        else
-        {
-            as[celli] = 0.0;
-        }
-    }
-    tas.ref().correctBoundaryConditions();
-    return tas;
-}
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::MossBrookesSoot::e(const label bandI) const
-{
-
-    return a(bandI);
-}
-
-void Foam::radiation::MossBrookesSoot::correct()
-{
-
-    calcSource();
-
-    const volScalarField &rho_ = thermo_.rho();
-
-    const volVectorField &U_ = turbulence().U();
-
-    volScalarField Dsoot (turbulence().mut()/PrtSoot);
-
-    const surfaceScalarField phi_
-    (
-        IOobject
-        (
-            "phiSoot",
-            mesh().time().timeName(),
-            mesh()
-        ),
-        linearInterpolate(rho_*U_) & mesh().Sf()
-    );
-
-    fvScalarMatrix SootEqn
-    (
-        fvm::div(phi_, soot_)
-      - fvm::laplacian(Dsoot, soot_)
-     ==
-        Snet_
-    );
-
-    SootEqn.relax();
-
-    SootEqn.solve();
-
-    fvScalarMatrix NucConcEqn
-    (
-        fvm::div(phi_, NucConc_)
-      - fvm::laplacian(Dsoot, NucConc_)
-     ==
-        NSnet_
-    );
-
-    NucConcEqn.relax();
-
-    NucConcEqn.solve();
-
-
-    soot_.max(0.0);
-    NucConc_.max(0.0);
-
-    NumDen_ = rho_*NucConc_*Nnorm;
-
-    sootVF_ = soot_*rho_/dimensionedScalar("rhoSoot", dimDensity, rhoSoot);
-}
-
-
-void Foam::radiation::MossBrookesSoot::calcSource()
-{
-    const volScalarField &rho_ = thermo_.rho();
-    const volScalarField &T_ = thermo_.T();
-    const volScalarField &p_ = thermo_.p();
-
-    const scalar precI = mixture_.species()[coeffsDict_.lookup("Precursor")];
-    const scalar growI = mixture_.species()[coeffsDict_.lookup("SurfaceGrowSpecie")];
-    const scalar oxidI = mixture_.species()[coeffsDict_.lookup("SootOxidizerSpecie")];
-
-    const volScalarField &Yprec_ = mixture_.Y(precI);
-    const volScalarField &Ygrow_ = mixture_.Y(growI);
-    const volScalarField &Yoxid_ = mixture_.Y(oxidI);
-
-    const volScalarField XPrec_(mixture_.W()*Yprec_/mixture_.W(precI));
-    const volScalarField XSgs_ (mixture_.W()*Ygrow_/mixture_.W(growI));
-    const volScalarField XOxid_(mixture_.W()*Yoxid_/mixture_.W(oxidI));
-
-    const dimensionedScalar NA (Foam::constant::physicoChemical::NA);
-    const dimensionedScalar kB (Foam::constant::physicoChemical::k);
-    const dimensionedScalar R (Foam::constant::physicoChemical::R);
-
-    const scalar pi314 = Foam::constant::mathematical::pi;
-
-    forAll(rho_, celli)
-    {
-        const scalar XPreci = XPrec_[celli];
-        const scalar XSgsi = XSgs_[celli];
-        const scalar XOxidi = XOxid_[celli];
-        const scalar Ti = T_[celli];
-        const scalar pi = p_[celli];
-        const scalar rhoi = rho_[celli];
-        const scalar sooti = soot_[celli];
-        const scalar NumDeni = NumDen_[celli];
-
-        // concentrations in [kmol/m3]
-        const scalar MolarConcPreci = (XPreci*pi)/(1e3*R.value()*Ti);
-        const scalar MolarConcSgsi  = (XSgsi *pi)/(1e3*R.value()*Ti);
-        const scalar MolarConcOxidi = (XOxidi*pi)/(1e3*R.value()*Ti);
-
-        if (NumDeni == 0.0)
-        {
-            sootMeanDiameter_[celli] = 0.0;
-            Sf_[celli] = 0.0;
-        }
-        else
-        {
-            const scalar parameter = (6.0/pi314)*(rhoi/rhoSoot)*(sooti/NumDeni);
-            sootMeanDiameter_[celli] = pow(parameter, 1./3.);
-            Sf_[celli] = pi314*pow(sootMeanDiameter_[celli],2.0)*NumDeni;
-        }
-
-        const scalar sootMeanDiameteri = sootMeanDiameter_[celli];
-
-        Snuc_[celli] = Mp*Calpha*exp(-cliqueNuc*Talpha/Ti)*MolarConcPreci;
-
-        Smassgrow_[celli] = Cgamma*MolarConcSgsi*exp(-cliqueMassGrow*Tgamma/Ti)*Sf_[celli];
-
-        Soxid_[celli] = Coxid*Cw*etaColl*MolarConcOxidi*sqrt(Ti)*Sf_[celli];
-
-        Snet_[celli] = Snuc_[celli] + Smassgrow_[celli] - Soxid_[celli];
-
-        NSnuc_[celli] = (Calpha*NA.value()*MolarConcPreci*exp(-cliqueNuc*Talpha/Ti))/Nnorm.value();
-
-        aggParam_[celli] = 24.0*kB.value()*Ti/rhoSoot;
-
-        NSagg_[celli] = (Cbeta*sqrt(aggParam_[celli])*sqrt(sootMeanDiameteri)*pow(NumDeni,2.0))/Nnorm.value();
-
-        NSnet_[celli] = NSnuc_[celli] - NSagg_[celli];
-    }
-}
-
-
-const Foam::compressibleTurbulenceModel&
-Foam::radiation::MossBrookesSoot::turbulence() const
-{
-    return turbulence_;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/radiationSootModel/MossBrookesSoot/MossBrookesSoot.H

@@ -1,280 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::radiation::MossBrookesSoot
-
-Description
-    This soot model solves two transport equation for soot mass fraction
-    and soot normalized nuclei concentration.
-    The detail of the model strictly follows Moss & Brookes model in fluent.
-
-    Parameter settings are available in constant/radiationProperties.
-
-    Setting example:
-
-    sootModel MossBrookesSoot<rhoChemistryModel, gasHaThermoPhysics>;
-
-    MossBrookesSootCoeffs
-    {
-        SootDensity            1800.0;
-        Calpha                 54.0;
-        Mp                     144;
-        Cgamma                 11700;
-        Coxid                  1.0;
-        Cw                     105.8125;
-        etaColl                0.04;
-        Talpha                 21100.0;
-        Tgamma                 12100.0;
-        AvogadroNumber         6.022e+26;
-        Cbeta                  1.0;
-        BoltzmannConst         1.38e-23;
-        Radiationeffect        true;
-
-        Precursor              C2H2;
-        SurfaceGrowSpecie      C2H2;
-        SootOxidizerSpecie     OH;
-    }
-
-SourceFiles
-    MossBrookesSoot.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef MossBrookesSoot_H
-#define MossBrookesSoot_H
-
-#include "sootModel.H"
-#include "HashTable.H"
-
-#include "fluidThermo.H"
-#include "basicSpecieMixture.H"
-#include "turbulentFluidThermoModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-namespace radiation
-{
-
-/*---------------------------------------------------------------------------*\
-                    Class MossBrookesSoot Declaration
-\*---------------------------------------------------------------------------*/
-class MossBrookesSoot
-:
-    public sootModel
-
-{
-
-    // Static functions
-
-        //- Check mixture in thermo
-        static const basicSpecieMixture& checkThermo
-        (
-            const fluidThermo&
-        );
-
-
-    // Private data
-
-        //- Soot mass fraction
-        volScalarField soot_;
-
-        //- Soot volume fraction
-        volScalarField sootVF_;
-
-        //- Normalized nuclei concentration of soot
-        volScalarField NucConc_;
-
-        //- Number density of soot
-        volScalarField NumDen_;
-
-        //- Soot model dictionary
-        dictionary coeffsDict_;
-
-        //- Thermo package
-        const fluidThermo& thermo_;
-
-        //- Composition
-        const basicSpecieMixture& mixture_;
-
-        //- Reference to the turbulence model
-        const compressibleTurbulenceModel &turbulence_;
-
-        //- Net rate of soot generation
-        volScalarField Snet_;
-
-        //- The rate of soot formation by nucleation
-        volScalarField Snuc_;
-
-        //- The rate of soot formation by mass growth
-        volScalarField Smassgrow_;
-
-        //- The rate of soot formation by soot oxidation
-        volScalarField Soxid_;
-
-        //- The rate of soot number density formation by nucleation
-        volScalarField NSnuc_;
-
-        //- The rate of soot number density formation by agglomeration
-        volScalarField NSagg_;
-
-        //- Net rate of soot numbe density formation
-        volScalarField NSnet_;
-
-        //- Parameter for Agglomeration source term
-        volScalarField aggParam_;
-
-        //- Available surface area of soot oxidation per volume
-        volScalarField Sf_;
-
-        //- Soot particle mean diameter
-        volScalarField sootMeanDiameter_;
-
-        //- soot number density normalization factor
-        const dimensionedScalar Nnorm;
-
-        //- Soot turbulent Prandtl number
-        scalar PrtSoot;
-
-        //- Density of soot particle
-        scalar rhoSoot;
-
-        //- Model constant for soot inception rate
-        scalar Calpha;
-
-        //- Mass of an incipient soot particle
-        scalar Mp;
-
-        //- Surface growth rate scaling factor
-        scalar Cgamma;
-
-        //- Oxidation rate scaling parameter
-        scalar Coxid;
-
-        //- Oxidation model constant
-        scalar Cw;
-
-        //- Collisional efficiency parameter
-        scalar etaColl;
-
-        //- Activation temperature for soot inception
-        scalar Talpha;
-
-        //- Activation temperature for soot growth rate
-        scalar Tgamma;
-
-        //- Model constant for coagulation rate
-        scalar Cbeta;
-
-        //- Soot radiative effect on/off flag
-        bool radiationActive_;
-
-        //- Correction constant for Exponent in Nucleation
-        scalar cliqueNuc;
-
-        //- Correction constant for Exponent in Soot Mass Grow
-        scalar cliqueMassGrow;
-
-        //- Soot precursor specie name
-        const word precSpecie;
-
-        //- Soot surface growth specie name
-        const word sgsSpecie;
-
-        //- Soot oxidation specie name
-        const word oxSpecie;
-
-protected:
-
-    // Protected data
-
-    // Protected Member Functions
-
-public:
-
-    //- Runtime type information
-    TypeName("MossBrookesSoot");
-
-
-    // Constructors
-
-        //- Construct from components
-        MossBrookesSoot
-        (
-            const dictionary& dict,
-            const fvMesh& mesh,
-            const word& modelType
-        );
-
-
-    //- Destructor
-    virtual ~MossBrookesSoot();
-
-
-        // Edit
-
-            //- Main update/correction routine
-            virtual void correct();
-
-            //- Calculation of the source term routine
-            virtual void calcSource();
-
-
-        // Access
-
-            //- Return Ysoot
-            const volScalarField& soot() const
-            {
-                return soot_;
-            }
-
-            //- Return access to turbulence
-            const compressibleTurbulenceModel& turbulence() const;
-
-            // Soot absorption coefficient
-
-                //- Soot absorption coefficient (net)
-                virtual tmp<volScalarField> a(const label bandI = 0) const;
-
-            // Soot emission coefficient
-
-                //- Soot emission coefficient (net)
-                virtual tmp<volScalarField> e(const label bandI = 0) const;
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-} // End namespace radiation
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/radiationSootModel/khanGreeveSoot/khanGreeveSoot.C

@@ -1,384 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "khanGreeveSoot.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvm.H"
-#include "basicSpecieMixture.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-namespace Foam
-{
-    namespace radiation
-    {
-        defineTypeNameAndDebug(khanGreeveSoot, 0);
-
-        addToRunTimeSelectionTable
-        (
-            sootModel,
-            khanGreeveSoot,
-            dictionary
-        );
-    }
-}
-
-
-const Foam::basicSpecieMixture&
-Foam::radiation::khanGreeveSoot::checkThermo
-(
-    const fluidThermo& thermo
-)
-{
-    if (isA<basicSpecieMixture>(thermo))
-    {
-        return dynamic_cast<const basicSpecieMixture& >
-        (
-            thermo
-        );
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Inconsistent thermo package for " << thermo.type()
-            << "Please select a thermo package based on "
-            << "basicSpecieMixture" << exit(FatalError);
-
-        return dynamic_cast<const basicSpecieMixture& >
-        (
-            thermo
-        );
-    }
-
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::radiation::khanGreeveSoot::khanGreeveSoot
-
-(
-    const dictionary& dict,
-    const fvMesh& mesh,
-    const word& modelType
-)
-:
-    sootModel(dict, mesh, modelType),
-    soot_
-    (
-        IOobject
-        (
-            "soot",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    sootVF_
-    (
-        IOobject
-        (
-            "sootVF",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    coeffsDict_(dict.subOrEmptyDict(modelType + "Coeffs")),
-    thermo_(mesh.lookupObject<fluidThermo>(basicThermo::dictName)),
-    mixture_(checkThermo(thermo_)),
-    turbulence_(mesh.lookupObject<compressibleTurbulenceModel>(turbulenceModel::propertiesName)),
-    Snet_
-    (
-        IOobject
-        (
-            "Snet",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Snet", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Sform_
-    (
-        IOobject
-        (
-            "Sform",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Sform", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    Scomb_
-    (
-        IOobject
-        (
-            "Scomb",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Scomb", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    S1_
-    (
-        IOobject
-        (
-            "S1",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("S1", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    S2_
-    (
-        IOobject
-        (
-            "S2",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("S2", dimMass/(pow3(dimLength)*dimTime), 0.0)
-    ),
-    PrtSoot(readScalar(coeffsDict_.lookup("TurbulentPrandtlNumberSoot"))),
-    rhoSoot(readScalar(coeffsDict_.lookup("SootDensity"))),
-    Cs_(readScalar(coeffsDict_.lookup("Cs"))),
-    r_(readScalar(coeffsDict_.lookup("r"))),
-    Ta_(readScalar(coeffsDict_.lookup("ActivationTemperature"))),
-    A_(readScalar(coeffsDict_.lookup("A"))),
-    nuFuel_(readScalar(coeffsDict_.lookup("nufuel"))),
-    nuSoot_(readScalar(coeffsDict_.lookup("nusoot"))),
-    eqv_
-    (
-        IOobject
-        (
-            "eqv",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("eqv", dimless, 0.0)
-    ),
-    eqvMax(readScalar(coeffsDict_.lookup("EquivalenceRatioMax"))),
-    eqvMin(readScalar(coeffsDict_.lookup("EquivalenceRatioMin"))),
-    radiationActive_(readBool(coeffsDict_.lookup("Radiationeffect")))
-{
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::radiation::khanGreeveSoot::~khanGreeveSoot()
-
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::khanGreeveSoot::a(const label bandI) const
-{
-
-    const volScalarField T_ = thermo_.T();
-
-    tmp<volScalarField> tas
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "aSoot" + name(bandI),
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("aSoot", dimless/dimLength, 0.0),
-            extrapolatedCalculatedFvPatchVectorField::typeName
-        )
-    );
-
-    scalarField& as = tas.ref().primitiveFieldRef();
-    if (radiationActive_)
-    {
-        forAll(as, celli)
-        {
-            as[celli] = (1232.4)*(1800.0)*soot_[celli]*(1.0+(4.8e-4)*(T_[celli]-2000.0));
-        }
-    }
-    tas.ref().correctBoundaryConditions();
-
-    return tas;
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::khanGreeveSoot::e(const label bandI) const
-{
-
-    return a(bandI);
-}
-
-
-void Foam::radiation::khanGreeveSoot::correct()
-{
-
-    calcSource();
-
-    const volScalarField &rho_ = thermo_.rho();
-
-    const volVectorField &U_ = turbulence().U();
-
-    const surfaceScalarField phi_
-    (
-        IOobject
-        (
-            "phiSoot",
-            mesh().time().timeName(),
-            mesh()
-        ),
-        linearInterpolate(rho_*U_) & mesh().Sf()
-    );
-
-    fvScalarMatrix SootEqn
-    (
-        fvm::ddt(rho_, soot_)
-      + fvm::div(phi_, soot_)
-      - fvm::laplacian((turbulence().mut()/PrtSoot), soot_)
-     ==
-        Snet_
-    );
-
-    SootEqn.relax();
-
-    SootEqn.solve();
-
-    forAll(rho_, celli)
-    {
-        if (soot_[celli] < 1e-30){
-            soot_[celli] = 0.0;
-        }
-
-        sootVF_[celli] = soot_[celli]*rho_[celli]/rhoSoot;
-
-    }
-
-}
-
-
-void Foam::radiation::khanGreeveSoot::calcSource()
-{
-    const volScalarField rho_ = thermo_.rho();
-    const volScalarField T_ = thermo_.T();
-    const volScalarField p_ = thermo_.p();
-
-    const scalar oxI = mixture_.species()[coeffsDict_.lookup("Oxidizer")];
-
-    const scalar fuelI = mixture_.species()[coeffsDict_.lookup("Fuel")];
-
-    const volScalarField &Yox_ = mixture_.Y(oxI);
-
-    const volScalarField &Yfuel_ = mixture_.Y(fuelI);
-
-    const volScalarField Xfuel_(mixture_.W()*Yfuel_/mixture_.W(fuelI));
-
-    const volScalarField &epsilon_ = turbulence().epsilon();
-
-    const volScalarField &k_ = turbulence().k();
-
-    forAll(rho_, celli)
-    {
-        const scalar Yoxi = Yox_[celli];
-        const scalar Yfueli = Yfuel_[celli];
-
-        eqv_[celli] = (Yfueli/Yoxi)*nuFuel_;
-
-        if (eqv_[celli] > eqvMax){
-            eqv_[celli] = 0.0;
-        }
-        else if (eqv_[celli] < eqvMin){
-            eqv_[celli] = 0.0;
-        }
-
-        const scalar pi = p_[celli];
-        const scalar Xfueli = Xfuel_[celli];
-        const scalar Ti = T_[celli];
-
-        Sform_[celli] = Cs_*pi*Xfueli*pow(eqv_[celli],r_)*exp(-Ta_/Ti);
-
-        const scalar rhoi = rho_[celli];
-        const scalar sooti = soot_[celli];
-        const scalar epsiloni = epsilon_[celli];
-        const scalar ki = k_[celli];
-
-        S1_[celli] = A_*rhoi*sooti*epsiloni/ki;
-
-        if (sooti == 0.0 and Yfueli == 0.0 )
-        {
-            S2_[celli] = 0.0;
-        }
-        else 
-        {
-            S2_[celli] = A_*rhoi*(Yoxi/nuSoot_)*(sooti*nuSoot_/(sooti*nuSoot_+(Yfueli)*nuFuel_))*epsiloni/ki;
-        }
-
-        Scomb_[celli] = min(S1_[celli], S2_[celli]);
-
-        Snet_[celli] = Sform_[celli] - Scomb_[celli];
-    }
-
-}
-
-
-const Foam::compressibleTurbulenceModel&
-Foam::radiation::khanGreeveSoot::turbulence() const
-{
-    return turbulence_;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //