ENH: Updated devolatilisation models to operate on a per-specie basis

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -287,9 +287,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         T0,
         mass0,
         this->mass0_,
-        idG,
+        YMix[GAS]*YGas_,
-        YMix[GAS],
-        YGas_,
         canCombust_,
         dMassDV,
         Sh,
@@ -494,9 +492,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     const scalar T,
     const scalar mass,
     const scalar mass0,
-    const label idVolatile,
+    const scalarField& YGasEff,
-    const scalar YVolatileTot,
-    const scalarField& YVolatile,
     bool& canCombust,
     scalarField& dMassDV,
     scalar& Sh,
@@ -510,30 +506,29 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     if
     (
         !td.cloud().devolatilisation().active()
-     || T < td.constProps().Tvap()
+     || T < td.cloud().constProps().Tvap()
     )
     {
         return;
     }
 
     // Total mass of volatiles evolved
-    const scalar dMassTot = td.cloud().devolatilisation().calculate
+    td.cloud().devolatilisation().calculate
     (
         dt,
         mass0,
         mass,
         T,
-        td.cloud().composition().YMixture0()[idVolatile],
+        YGasEff,
-        YVolatileTot,
+        canCombust,
-        canCombust
+        dMassDV
     );
 
-    // Volatile mass transfer - equal components of each volatile specie
+    scalar dMassTot = sum(dMassDV);
-    dMassDV = YVolatile*dMassTot;
 
     td.cloud().addToMassDevolatilisation(this->nParticle_*dMassTot);
 
-    Sh -= dMassTot*td.constProps().LDevol()/dt;
+    Sh -= dMassTot*td.cloud().constProps().LDevol()/dt;
 
     // Molar average molecular weight of carrier mix
     const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
@@ -546,8 +541,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
         const label id =
             td.cloud().composition().localToGlobalCarrierId(GAS, i);
-        const scalar Cp = td.cloud().thermo().carrier().Cp(id, Ts);
+        const scalar Cp = td.cloud().mcCarrierThermo().speciesData()[id].Cp(Ts);
-        const scalar W = td.cloud().thermo().carrier().W(id);
+        const scalar W = td.cloud().mcCarrierThermo().speciesData()[id].W();
         const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
         // Dab calc'd using API vapour mass diffusivity function

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H

@@ -225,9 +225,7 @@ protected:
             const scalar T,            // temperature
             const scalar mass,         // mass
             const scalar mass0,        // mass (initial on injection)
-            const label idVolatile,    // id of volatile phase
+            const scalarField& YGasEff,// Gas component mass fractions
-            const scalar YVolatileTot, // total volatile mass fraction
-            const scalarField& YVolatile, // volatile component mass fractions
             bool& canCombust,          // 'can combust' flag
             scalarField& dMassDV,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C

@@ -8,10 +8,10 @@
 License
     This file is part of OpenFOAM.
 
-    OpenFOAM is free software: you can redistribute it and/or modify it
+    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
+    under the terms of the GNU General Public License as published by the
-    the Free Software Foundation, either version 3 of the License, or
+    Free Software Foundation; either version 2 of the License, or (at your
-    (at your option) any later version.
+    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
@@ -19,7 +19,8 @@ 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/>.
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 \*---------------------------------------------------------------------------*/
 
@@ -35,12 +36,46 @@ Foam::ConstantRateDevolatilisation<CloudType>::ConstantRateDevolatilisation
 )
 :
     DevolatilisationModel<CloudType>(dict, owner, typeName),
-    A0_(dimensionedScalar(this->coeffDict().lookup("A0")).value()),
+    volatileData_(this->coeffDict().lookup("volatileData")),
-    volatileResidualCoeff_
+    YVolatile0_(volatileData_.size()),
+    volatileToGasMap_(volatileData_.size()),
+    residualCoeff_(readScalar(this->coeffDict().lookup("residualCoeff")))
+{
+    if (volatileData_.empty())
+    {
+        WarningIn
         (
-        readScalar(this->coeffDict().lookup("volatileResidualCoeff"))
+            "Foam::ConstantRateDevolatilisation<CloudType>::"
-    )
+            "ConstantRateDevolatilisation"
-{}
+            "("
+                "const dictionary& dict, "
+                "CloudType& owner"
+            ")"
+        )   << "Devolatilisation model selected, but no volatiles defined"
+            << nl << endl;
+    }
+    else
+    {
+        Info<< "Participating volatile species:" << endl;
+
+        // Determine mapping between active volatiles and cloud gas components
+        const label idGas = owner.composition().idGas();
+        const scalar YGasTot = owner.composition().YMixture0()[idGas];
+        const scalarField& YGas = owner.composition().Y0(idGas);
+        forAll(volatileData_, i)
+        {
+            const word& specieName = volatileData_[i].first();
+            Info<< "    " << specieName << endl;
+
+            const label id = owner.composition().localId(idGas, specieName);
+            volatileToGasMap_[i] = id;
+            YVolatile0_[i] = YGasTot*YGas[id];
+
+            Info<< "    " << specieName << ": particle mass fraction = "
+                << YVolatile0_[i] << endl;
+        }
+    }
+}
 
 
 template<class CloudType>
@@ -50,8 +85,10 @@ Foam::ConstantRateDevolatilisation<CloudType>::ConstantRateDevolatilisation
 )
 :
     DevolatilisationModel<CloudType>(dm),
-    A0_(dm.A0_),
+    volatileData_(dm.volatileData_),
-    volatileResidualCoeff_(dm.volatileResidualCoeff_)
+    YVolatile0_(dm.YVolatile0_),
+    volatileToGasMap_(dm.volatileToGasMap_),
+    residualCoeff_(dm.residualCoeff_)
 {}
 
 
@@ -65,29 +102,34 @@ Foam::ConstantRateDevolatilisation<CloudType>::~ConstantRateDevolatilisation()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
-Foam::scalar Foam::ConstantRateDevolatilisation<CloudType>::calculate
+void Foam::ConstantRateDevolatilisation<CloudType>::calculate
 (
     const scalar dt,
     const scalar mass0,
     const scalar mass,
     const scalar T,
-    const scalar YVolatile0,
+    const scalarField& YGasEff,
-    const scalar YVolatile,
+    bool& canCombust,
-    bool& canCombust
+    scalarField& dMassDV
 ) const
 {
-    const scalar massVolatile0 = YVolatile0*mass0;
+    forAll(volatileData_, i)
-    const scalar massVolatile  = YVolatile*mass;
-
-    if (massVolatile <= volatileResidualCoeff_*massVolatile0)
     {
-        canCombust = true;
+        const label id = volatileToGasMap_[i];
+        const scalar massVolatile0 = mass0*YVolatile0_[id];
+        const scalar massVolatile = mass*YGasEff[id];
+
+        // Combustion allowed once all volatile components evolved
+        canCombust =
+            canCombust
+         && (massVolatile <= residualCoeff_*massVolatile0);
+
+        // Model coefficients
+        const scalar A0 = volatileData_[i].second();
+
+        // Mass transferred from particle to carrier gas phase
+        dMassDV = min(dt*A0*massVolatile0, massVolatile);
     }
-
-    // Volatile devolatilisation from particle to carrier gas phase
-    const scalar dMass = min(dt*A0_*massVolatile0, massVolatile);
-
-    return dMass;
 }
 
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H

@@ -52,13 +52,19 @@ class ConstantRateDevolatilisation
 
         // Model constants
 
-            //- Rate constant (suggested default = 12) [1/s]
+            //- List of volatile data - (name A0)
-            const scalar A0_;
+            List<Tuple2<word, scalar> > volatileData_;
+
+            //- List of initial volatile mass fractions
+            List<scalar> YVolatile0_;
+
+            //- Mapping between local and cloud gaseous species
+            List<label> volatileToGasMap_;
 
             //- Volatile residual coefficient (0-1)
             //  When the fraction of volatiles are depleted below this
             //  threshold, combustion can occur
-            const scalar volatileResidualCoeff_;
+            const scalar residualCoeff_;
 
 
 public:
@@ -95,15 +101,15 @@ public:
     // Member Functions
 
         //- Update model
-        virtual scalar calculate
+        virtual void calculate
         (
             const scalar dt,
             const scalar mass0,
             const scalar mass,
             const scalar T,
-            const scalar YVolatile0,
+            const scalarField& YGasEff,
-            const scalar YVolatile,
+            bool& canCombust,
-            bool& canCombust
+            scalarField& dMassDV
         ) const;
 };
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.C

@@ -69,15 +69,15 @@ Foam::DevolatilisationModel<CloudType>::~DevolatilisationModel()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
-Foam::scalar Foam::DevolatilisationModel<CloudType>::calculate
+void Foam::DevolatilisationModel<CloudType>::calculate
 (
     const scalar,
     const scalar,
     const scalar,
     const scalar,
-    const scalar,
+    const scalarField&,
-    const scalar,
+    bool&,
-    bool&
+    scalarField&
 ) const
 {
     notImplemented
@@ -88,13 +88,11 @@ Foam::scalar Foam::DevolatilisationModel<CloudType>::calculate
             "const scalar, "
             "const scalar, "
             "const scalar, "
-            "const scalar, "
+            "const scalarField&, "
-            "const scalar, "
             "bool&, "
+            "scalarField&"
         ") const"
     );
-
-    return 0.0;
 }
 
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H

@@ -115,15 +115,15 @@ public:
     // Member Functions
 
         //- Update model
-        virtual scalar calculate
+        virtual void calculate
         (
             const scalar dt,
             const scalar mass0,
             const scalar mass,
             const scalar T,
-            const scalar YVolatile0,
+            const scalarField& YGasEff,
-            const scalar YVolatile,
+            bool& canCombust,
-            bool& canCombust
+            scalarField& dMassDV
         ) const;
 };
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C

@@ -65,21 +65,19 @@ bool Foam::NoDevolatilisation<CloudType>::active() const
 
 
 template<class CloudType>
-Foam::scalar Foam::NoDevolatilisation<CloudType>::calculate
+void Foam::NoDevolatilisation<CloudType>::calculate
 (
     const scalar,
     const scalar,
     const scalar,
     const scalar,
-    const scalar,
+    const scalarField&,
-    const scalar,
+    bool& canCombust,
-    bool& canCombust
+    scalarField&
 ) const
 {
     // Model does not stop combustion taking place
     canCombust = true;
-
-    return 0.0;
 }
 
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H

@@ -82,15 +82,15 @@ public:
         virtual bool active() const;
 
         //- Update model
-        virtual scalar calculate
+        virtual void calculate
         (
-            const scalar,
+            const scalar dt,
-            const scalar,
+            const scalar mass0,
-            const scalar,
+            const scalar mass,
-            const scalar,
+            const scalar T,
-            const scalar,
+            const scalarField& YGasEff,
-            const scalar,
+            bool& canCombust,
-            bool&
+            scalarField& dMassDV
         ) const;
 };
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C

@@ -8,10 +8,10 @@
 License
     This file is part of OpenFOAM.
 
-    OpenFOAM is free software: you can redistribute it and/or modify it
+    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
+    under the terms of the GNU General Public License as published by the
-    the Free Software Foundation, either version 3 of the License, or
+    Free Software Foundation; either version 2 of the License, or (at your
-    (at your option) any later version.
+    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
@@ -19,7 +19,8 @@ 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/>.
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 \*---------------------------------------------------------------------------*/
 
@@ -36,13 +37,44 @@ SingleKineticRateDevolatilisation
 )
 :
     DevolatilisationModel<CloudType>(dict, owner, typeName),
-    A1_(dimensionedScalar(this->coeffDict().lookup("A1")).value()),
+    volatileData_(this->coeffDict().lookup("volatileData")),
-    E_(dimensionedScalar(this->coeffDict().lookup("E")).value()),
+    YVolatile0_(volatileData_.size()),
-    volatileResidualCoeff_
+    volatileToGasMap_(volatileData_.size()),
+    residualCoeff_(readScalar(this->coeffDict().lookup("residualCoeff")))
+{
+    if (volatileData_.empty())
+    {
+        WarningIn
         (
-        readScalar(this->coeffDict().lookup("volatileResidualCoeff"))
+            "Foam::SingleKineticRateDevolatilisation<CloudType>::"
-    )
+            "SingleKineticRateDevolatilisation"
-{}
+            "("
+                "const dictionary& dict, "
+                "CloudType& owner"
+            ")"
+        )   << "Devolatilisation model selected, but no volatiles defined"
+            << nl << endl;
+    }
+    else
+    {
+        Info<< "Participating volatile species:" << endl;
+
+        // Determine mapping between active volatiles and cloud gas components
+        const label idGas = owner.composition().idGas();
+        const scalar YGasTot = owner.composition().YMixture0()[idGas];
+        const scalarField& YGas = owner.composition().Y0(idGas);
+        forAll(volatileData_, i)
+        {
+            const word& specieName = volatileData_[i].name();
+            const label id = owner.composition().localId(idGas, specieName);
+            volatileToGasMap_[i] = id;
+            YVolatile0_[i] = YGasTot*YGas[id];
+
+            Info<< "    " << specieName << ": particle mass fraction = "
+                << YVolatile0_[i] << endl;
+        }
+    }
+}
 
 
 template<class CloudType>
@@ -53,9 +85,10 @@ SingleKineticRateDevolatilisation
 )
 :
     DevolatilisationModel<CloudType>(dm),
-    A1_(dm.A1_),
+    volatileData_(dm.volatileData_),
-    E_(dm.E_),
+    YVolatile0_(dm.YVolatile0_),
-    volatileResidualCoeff_(dm.volatileResidualCoeff_)
+    volatileToGasMap_(dm.volatileToGasMap_),
+    residualCoeff_(dm.residualCoeff_)
 {}
 
 
@@ -70,32 +103,39 @@ Foam::SingleKineticRateDevolatilisation<CloudType>::
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
-Foam::scalar Foam::SingleKineticRateDevolatilisation<CloudType>::calculate
+void Foam::SingleKineticRateDevolatilisation<CloudType>::calculate
 (
     const scalar dt,
     const scalar mass0,
     const scalar mass,
     const scalar T,
-    const scalar YVolatile0,
+    const scalarField& YGasEff,
-    const scalar YVolatile,
+    bool& canCombust,
-    bool& canCombust
+    scalarField& dMassDV
 ) const
 {
-    const scalar massVolatile0 = YVolatile0*mass;
+    bool done = true;
-    const scalar massVolatile = YVolatile*mass;
+    forAll(volatileData_, i)
-
-    if (massVolatile <= volatileResidualCoeff_*massVolatile0)
     {
-        canCombust = true;
+        const label id = volatileToGasMap_[i];
-    }
+        const scalar massVolatile0 = mass0*YVolatile0_[id];
+        const scalar massVolatile = mass*YGasEff[id];
+
+        // Combustion allowed once all volatile components evolved
+        done = done && (massVolatile <= residualCoeff_*massVolatile0);
+
+        // Model coefficients
+        const scalar A1 = volatileData_[i].A1();
+        const scalar E = volatileData_[i].E();
 
         // Kinetic rate
-    const scalar kappa = A1_*exp(-E_/(specie::RR*T));
+        const scalar kappa = A1*exp(-E/(specie::RR*T));
 
-    // Volatile devolatilisation from particle to carrier gas phase
+        // Mass transferred from particle to carrier gas phase
-    const scalar dMass = min(dt*kappa*massVolatile, massVolatile);
+        dMassDV[id] = min(dt*kappa*massVolatile, massVolatile);
+    }
 
-    return dMass;
+    canCombust = done;
 }
 
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H

@@ -8,10 +8,10 @@
 License
     This file is part of OpenFOAM.
 
-    OpenFOAM is free software: you can redistribute it and/or modify it
+    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
+    under the terms of the GNU General Public License as published by the
-    the Free Software Foundation, either version 3 of the License, or
+    Free Software Foundation; either version 2 of the License, or (at your
-    (at your option) any later version.
+    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
@@ -19,13 +19,38 @@ 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/>.
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
     Foam::SingleKineticRateDevolatilisation
 
 Description
-    Single kinetic rate devolatisation model
+    Single kinetic rate devolatisation model.
+    - acts on a per-specie basis
+    - Rate given by Arrhenius eqn
+
+        kappa = A1.exp(- E/R.T)
+
+    Where:
+        kappa = rate constant
+        A1    = activation energy (user input)
+        E     = pre-exponential factor (user input)
+        R     = universal gas constant
+        T     = temperature
+
+    Usage:
+
+        SingleKineticRateDevolatilisationCoeffs
+        {
+            volatileData
+            (
+                (CH4     12     0.5)   // (name A1 E)
+                (CO2     12     0.5)   // (name A1 E)
+            );
+
+            volatileResidualCoeff 1e-6;
+        }
 
 \*---------------------------------------------------------------------------*/
 
@@ -47,20 +72,121 @@ class SingleKineticRateDevolatilisation
 :
     public DevolatilisationModel<CloudType>
 {
+    // Helper class to store specie-local volatile data
+    class volatileData
+    {
+        // Private data
+
+            //- Specie name
+            word name_;
+
+            //- Activation energy
+            scalar A1_;
+
+            //- Pre-exponential factor
+            scalar E_;
+
+
+    public:
+
+        // Constructors
+
+            //- Null constructor
+            volatileData()
+            :
+                name_(word::null),
+                A1_(0.0),
+                E_(0.0)
+            {}
+
+            //- Construct from Istream
+            volatileData(Istream& is)
+            :
+                name_(is),
+                A1_(readScalar(is)),
+                E_(readScalar(is))
+            {}
+
+            //- Construct as copy
+            volatileData(const volatileData& vd)
+            :
+                name_(vd.name_),
+                A1_(vd.A1_),
+                E_(vd.E_)
+            {}
+
+
+        //- Destructor
+        ~volatileData()
+        {}
+
+
+        // Public Member Functions
+
+            // Access
+
+                //- Return const access to the name
+                const word& name() const
+                {
+                    return name_;
+                }
+
+                //- Return const access to the activation energy
+                scalar A1() const
+                {
+                    return A1_;
+                }
+
+                //- Return const access to the pre-exponential factor
+                scalar E() const
+                {
+                    return E_;
+                }
+
+
+        // IOstream Operators
+
+            //- Read from Istream
+            friend Istream& operator>>(Istream& is, volatileData& vd)
+            {
+                is.readBeginList("volatileData");
+                is  >>  vd.name_  >> vd.A1_  >> vd.E_;
+                is.readEndList("volatileData");
+
+                return is;
+            }
+
+            //- Write to Ostream
+            friend Ostream& operator<<(Ostream& os, const volatileData& vd)
+            {
+                os  << token::BEGIN_LIST
+                    << vd.name_ << token::SPACE
+                    << vd.A1_ << token::SPACE
+                    << vd.E_
+                    << token::END_LIST;
+
+                return os;
+            }
+    };
+
+
     // Private data
 
         // Model constants
 
-            //- Activation energy
+            //- List of volatile data - (name A1 E)
-            const scalar A1_;
+            List<volatileData> volatileData_;
 
-            //- Pre-exoponential factor
+            //- List of initial volatile mass fractions
-            const scalar E_;
+            List<scalar> YVolatile0_;
+
+            //- Mapping between local and cloud gaseous species
+            List<label> volatileToGasMap_;
 
             //- Volatile residual coefficient (0-1)
             //  When the fraction of volatiles are depleted below this
             //  threshold, combustion can occur
-            const scalar volatileResidualCoeff_;
+            const scalar residualCoeff_;
 
 
 public:
@@ -101,15 +227,15 @@ public:
     // Member Functions
 
         //- Update model
-        virtual scalar calculate
+        virtual void calculate
         (
             const scalar dt,
             const scalar mass0,
             const scalar mass,
             const scalar T,
-            const scalar YVolatile0,
+            const scalarField& YGasEff,
-            const scalar YVolatile,
+            bool& canCombust,
-            bool& canCombust
+            scalarField& dMassDV
         ) const;
 };