Merge branch 'EvaporationModels' into 'trunk'

added convection-diffusion-controlled liquid evaporation model

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

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

@@ -96,7 +96,8 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         pc_,
         this->Tc_,
         X,
-        dMassPC
+        dMassPC,
+        this->rhoc_
     );
 
     // Limit phase change mass by availability of each specie

src/lagrangian/intermediate/parcels/include/makeReactingParcelPhaseChangeModels.H

@@ -30,6 +30,7 @@ License
 
 #include "NoPhaseChange.H"
 #include "LiquidEvaporation.H"
+#include "LiquidEvaporationConvDiff.H"
 #include "LiquidEvaporationBoil.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -39,6 +40,7 @@ License
     makePhaseChangeModel(CloudType);                                           \
     makePhaseChangeModelType(NoPhaseChange, CloudType);                        \
     makePhaseChangeModelType(LiquidEvaporation, CloudType);                    \
+    makePhaseChangeModelType(LiquidEvaporationConvDiff, CloudType);            \
     makePhaseChangeModelType(LiquidEvaporationBoil, CloudType);
 
 

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C

@@ -142,7 +142,8 @@ void Foam::LiquidEvaporation<CloudType>::calculate
     const scalar pc,
     const scalar Tc,
     const scalarField& X,
-    scalarField& dMassPC
+    scalarField& dMassPC,
+    const scalar rhoc
 ) const
 {
     // immediately evaporate mass that has reached critical condition
@@ -196,8 +197,8 @@ void Foam::LiquidEvaporation<CloudType>::calculate
         // vapour concentration at surface [kmol/m3] at film temperature
         const scalar Cs = pSat/(RR*Ts);
 
-        // vapour concentration in bulk gas [kmol/m3] at film temperature
+        // vapour concentration in bulk gas [kmol/m3] at bulk gas temperature
-        const scalar Cinf = Xc[gid]*pc/(RR*Ts);
+        const scalar Cinf = Xc[gid]*pc/(RR*Tc);
 
         // molar flux of vapour [kmol/m2/s]
         const scalar Ni = max(kc*(Cs - Cinf), 0.0);

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H

@@ -119,7 +119,8 @@ public:
             const scalar pc,
             const scalar Tc,
             const scalarField& X,
-            scalarField& dMassPC
+            scalarField& dMassPC,
+            const scalar rhoc
         ) const;
 
         //- Return the enthalpy per unit mass

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C

@@ -142,7 +142,8 @@ void Foam::LiquidEvaporationBoil<CloudType>::calculate
     const scalar pc,
     const scalar Tc,
     const scalarField& X,
-    scalarField& dMassPC
+    scalarField& dMassPC,
+    const scalar rhoc
 ) const
 {
     // immediately evaporate mass that has reached critical condition

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H

@@ -129,7 +129,8 @@ public:
             const scalar pc,
             const scalar Tc,
             const scalarField& X,
-            scalarField& dMassPC
+            scalarField& dMassPC,
+            const scalar rhoc
         ) const;
 
         //- Return the enthalpy per unit mass

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationConvDiff/LiquidEvaporationConvDiff.C

@@ -0,0 +1,287 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "LiquidEvaporationConvDiff.H"
+#include "specie.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant::mathematical;
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::tmp<Foam::scalarField> Foam::LiquidEvaporationConvDiff<CloudType>::calcXc
+(
+    const label celli
+) const
+{
+    scalarField Xc(this->owner().thermo().carrier().Y().size());
+
+    forAll(Xc, i)
+    {
+        Xc[i] =
+            this->owner().thermo().carrier().Y()[i][celli]
+           /this->owner().thermo().carrier().W(i);
+    }
+
+    return Xc/sum(Xc);
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporationConvDiff<CloudType>::Sh
+(
+    const scalar Re,
+    const scalar Sc
+) const
+{
+    return 2.0 + 0.6*Foam::sqrt(Re)*cbrt(Sc);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::LiquidEvaporationConvDiff<CloudType>::LiquidEvaporationConvDiff
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    PhaseChangeModel<CloudType>(dict, owner, typeName),
+    liquids_(owner.thermo().liquids()),
+    activeLiquids_(this->coeffDict().lookup("activeLiquids")),
+    liqToCarrierMap_(activeLiquids_.size(), -1),
+    liqToLiqMap_(activeLiquids_.size(), -1)
+{
+    if (activeLiquids_.size() == 0)
+    {
+        WarningInFunction
+            << "Evaporation model selected, but no active liquids defined"
+            << nl << endl;
+    }
+    else
+    {
+        Info<< "Participating liquid species:" << endl;
+
+        // Determine mapping between liquid and carrier phase species
+        forAll(activeLiquids_, i)
+        {
+            Info<< "    " << activeLiquids_[i] << endl;
+            liqToCarrierMap_[i] =
+                owner.composition().carrierId(activeLiquids_[i]);
+        }
+
+        // Determine mapping between model active liquids and global liquids
+        const label idLiquid = owner.composition().idLiquid();
+        forAll(activeLiquids_, i)
+        {
+            liqToLiqMap_[i] =
+                owner.composition().localId(idLiquid, activeLiquids_[i]);
+        }
+    }
+}
+
+
+template<class CloudType>
+Foam::LiquidEvaporationConvDiff<CloudType>::LiquidEvaporationConvDiff
+(
+    const LiquidEvaporationConvDiff<CloudType>& pcm
+)
+:
+    PhaseChangeModel<CloudType>(pcm),
+    liquids_(pcm.owner().thermo().liquids()),
+    activeLiquids_(pcm.activeLiquids_),
+    liqToCarrierMap_(pcm.liqToCarrierMap_),
+    liqToLiqMap_(pcm.liqToLiqMap_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::LiquidEvaporationConvDiff<CloudType>::~LiquidEvaporationConvDiff()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::LiquidEvaporationConvDiff<CloudType>::calculate
+(
+    const scalar dt,
+    const label celli,
+    const scalar Re,
+    const scalar Pr,
+    const scalar d,
+    const scalar nu,
+    const scalar T,
+    const scalar Ts,
+    const scalar pc,
+    const scalar Tc,
+    const scalarField& X,
+    scalarField& dMassPC, 
+    const scalar rhoc     //POSECH
+) const
+{
+    // immediately evaporate mass that has reached critical condition
+    if ((liquids_.Tc(X) - T) < SMALL)
+    {
+        if (debug)
+        {
+            WarningInFunction
+                << "Parcel reached critical conditions: "
+                << "evaporating all avaliable mass" << endl;
+        }
+
+        forAll(activeLiquids_, i)
+        {
+            const label lid = liqToLiqMap_[i];
+            dMassPC[lid] = GREAT;
+        }
+
+        return;
+    }
+
+    // construct carrier phase species volume fractions for cell, celli
+    const scalarField Xc(calcXc(celli));
+
+    // calculate mass transfer of each specie in liquid
+    forAll(activeLiquids_, i)
+    {
+        const label gid = liqToCarrierMap_[i];
+        const label lid = liqToLiqMap_[i];
+
+        // vapour diffusivity [m2/s]
+        const scalar Dab = liquids_.properties()[lid].D(pc, Ts);
+
+        // saturation pressure for species i [pa]
+        // - carrier phase pressure assumed equal to the liquid vapour pressure
+        //   close to the surface
+        // NOTE: if pSat > pc then particle is superheated
+        // calculated evaporation rate will be greater than that of a particle
+        // at boiling point, but this is not a boiling model
+        const scalar pSat = liquids_.properties()[lid].pv(pc, T);
+
+        // Schmidt number
+        const scalar Sc = nu/(Dab + ROOTVSMALL);
+
+        // Sherwood number
+        const scalar Sh = this->Sh(Re, Sc);
+
+        // mass transfer coefficient [m/s]
+        const scalar kc = Sh*Dab/(d + ROOTVSMALL);
+
+        // vapour concentration at surface [kmol/m3] at film temperature
+        const scalar Cs = pSat/(RR*Ts);
+
+        // vapour concentration in bulk gas [kmol/m3] at bulk gas temperature
+        const scalar Cinf = Xc[gid]*pc/(RR*Tc);
+
+        // molar flux of vapour [kmol/m2/s]
+        const scalar Ni = max(kc*(Cs - Cinf), 0.0);
+
+        // droplet surface pressure assumed to surface vapour pressure
+        const scalar ps = liquids_.pv(pc, Ts, X);
+
+        // vapour density at droplet surface [kg/m3]
+        const scalar rhos = ps*liquids_.W(X)/(RR*Ts);
+
+        // vapour mass fraction at surface
+        const scalar Ys = min((Cs*liquids_.W(X)/rhos),0.9);
+
+        // vapour mass fraction in bulk gas
+        const scalar Yinf = Cinf*liquids_.W(X)/rhos;
+
+        // Spalding mass number
+        const scalar Bm = max(((Ys - Yinf)/(1.0 - Ys)),0.0);
+
+        // mass transfer [kg]
+        dMassPC[lid] += kc*pi*sqr(d)*rhoc*log(1+Bm)*dt;
+
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporationConvDiff<CloudType>::dh
+(
+    const label idc,
+    const label idl,
+    const scalar p,
+    const scalar T
+) const
+{
+    scalar dh = 0;
+
+    typedef PhaseChangeModel<CloudType> parent;
+    switch (parent::enthalpyTransfer_)
+    {
+        case (parent::etLatentHeat):
+        {
+            dh = liquids_.properties()[idl].hl(p, T);
+            break;
+        }
+        case (parent::etEnthalpyDifference):
+        {
+            scalar hc = this->owner().composition().carrier().Ha(idc, p, T);
+            scalar hp = liquids_.properties()[idl].h(p, T);
+
+            dh = hc - hp;
+            break;
+        }
+        default:
+        {
+            FatalErrorInFunction
+                << "Unknown enthalpyTransfer type" << abort(FatalError);
+        }
+    }
+
+    return dh;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporationConvDiff<CloudType>::Tvap
+(
+    const scalarField& X
+) const
+{
+    return liquids_.Tpt(X);
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporationConvDiff<CloudType>::TMax
+(
+    const scalar p,
+    const scalarField& X
+) const
+{
+    return liquids_.pvInvert(p, X);
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationConvDiff/LiquidEvaporationConvDiff.H

@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::LiquidEvaporationConvDiff
+
+Description
+    Liquid evaporation model
+    - uses ideal gas assumption
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LiquidEvaporationConvDiff_H
+#define LiquidEvaporationConvDiff_H
+
+#include "PhaseChangeModel.H"
+#include "liquidMixtureProperties.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                     Class LiquidEvaporationConvDiff Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class LiquidEvaporationConvDiff
+:
+    public PhaseChangeModel<CloudType>
+{
+protected:
+
+    // Protected data
+
+        //- Global liquid properties data
+        const liquidMixtureProperties& liquids_;
+
+        //- List of active liquid names
+        List<word> activeLiquids_;
+
+        //- Mapping between liquid and carrier species
+        List<label> liqToCarrierMap_;
+
+        //- Mapping between local and global liquid species
+        List<label> liqToLiqMap_;
+
+
+    // Protected Member Functions
+
+        //- Sherwood number as a function of Reynolds and Schmidt numbers
+        scalar Sh(const scalar Re, const scalar Sc) const;
+
+        //- Calculate the carrier phase component volume fractions at celli
+        tmp<scalarField> calcXc(const label celli) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("liquidEvaporationConvDiff");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        LiquidEvaporationConvDiff(const dictionary& dict, CloudType& cloud);
+
+        //- Construct copy
+        LiquidEvaporationConvDiff(const LiquidEvaporationConvDiff<CloudType>& pcm);
+
+        //- Construct and return a clone
+        virtual autoPtr<PhaseChangeModel<CloudType>> clone() const
+        {
+            return autoPtr<PhaseChangeModel<CloudType>>
+            (
+                new LiquidEvaporationConvDiff<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~LiquidEvaporationConvDiff();
+
+
+    // Member Functions
+
+        //- Update model
+        virtual void calculate
+        (
+            const scalar dt,
+            const label celli,
+            const scalar Re,
+            const scalar Pr,
+            const scalar d,
+            const scalar nu,
+            const scalar T,
+            const scalar Ts,
+            const scalar pc,
+            const scalar Tc,
+            const scalarField& X,
+            scalarField& dMassPC,
+            const scalar rhoc
+        ) const;
+
+        //- Return the enthalpy per unit mass
+        virtual scalar dh
+        (
+            const label idc,
+            const label idl,
+            const scalar p,
+            const scalar T
+        ) const;
+
+        //- Return vapourisation temperature
+        virtual scalar Tvap(const scalarField& X) const;
+
+        //- Return maximum/limiting temperature
+        virtual scalar TMax(const scalar p, const scalarField& X) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "LiquidEvaporationConvDiff.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C

@@ -78,7 +78,8 @@ void Foam::NoPhaseChange<CloudType>::calculate
     const scalar pc,
     const scalar Tc,
     const scalarField& X,
-    scalarField& dMassPC
+    scalarField& dMassPC,
+    const scalar rhoc
 ) const
 {
     // Nothing to do...

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H

@@ -94,7 +94,8 @@ public:
             const scalar pc,
             const scalar Tc,
             const scalarField& X,
-            scalarField& dMassPC
+            scalarField& dMassPC,
+            const scalar rhoc
         ) const;
 };
 

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H

@@ -166,7 +166,8 @@ public:
             const scalar pc,
             const scalar Tc,
             const scalarField& X,
-            scalarField& dMassPC
+            scalarField& dMassPC,
+            const scalar rhoc
         ) const = 0;
 
         //- Return the enthalpy per unit mass