adjust update order for lagrangian parcels and add new evaporation model

src/lagrangian/Allwmake

@@ -7,5 +7,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 wmake $targetType intermediate
 wmake $targetType turbulence
 wmake $targetType spray
+wmake $targetType intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationSpalding
 
 #------------------------------------------------------------------------------

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

@@ -394,7 +394,7 @@ void Foam::ReactingParcel<ParcelType>::calc
     typedef typename TrackCloudType::reactingCloudType reactingCloudType;
     const CompositionModel<reactingCloudType>& composition =
         cloud.composition();
-
+    const SLGThermo& thermo = cloud.thermo();
 
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -404,13 +404,34 @@ void Foam::ReactingParcel<ParcelType>::calc
     const vector& U0 = this->U_;
     const scalar T0 = this->T_;
     const scalar mass0 = this->mass();
+    const scalar RR = 1000.0*constant::physicoChemical::R.value();
 
-
-    // Calc surface values
     scalar Ts, rhos, mus, Prs, kappas;
-    this->calcSurfaceValues(cloud, td, T0, Ts, rhos, mus, Prs, kappas);
-    scalar Res = this->Re(rhos, U0, td.Uc(), d0, mus);
 
+    // Calc carrier thermo properties
+    scalarField Yg(thermo.carrier().species().size());
+    scalar rhog = td.pc()*thermo.carrier().W()/(RR*td.Tc());
+    scalar mug = 0;
+    scalar kappag = 0;
+    scalar Cpg = 0;
+
+    thermo.carrier().calcMu(td.Tc(), td.pc());
+    thermo.carrier().calcCp(td.Tc(), td.pc());
+    forAll(Yg, i)
+    {
+        Yg[i] = thermo.carrier().Y(i)[this->cell()];
+        mug += Yg[i]*thermo.carrier().mu(i, td.pc(), td.Tc());
+        kappag += Yg[i]*thermo.carrier().kappa(i, td.pc(), td.Tc());
+        Cpg += Yg[i]*thermo.carrier().Cp(i, td.pc(), td.Tc());
+    }
+
+    rhog = max(rhog, rootVSmall);
+    mug = max(mug, rootVSmall);
+    kappag = max(kappag, rootVSmall);
+    Cpg = max(Cpg, rootVSmall);
+
+    scalar Prg = Cpg*mug/kappag;
+    Prg = max(Prg,rootVSmall);
 
     // Sources
     // ~~~~~~~
@@ -452,13 +473,21 @@ void Foam::ReactingParcel<ParcelType>::calc
     // Surface concentrations of emitted species
     scalarField Cs(composition.carrier().species().size(), 0.0);
 
+    this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
+
+    this->calcSurfaceValues(cloud, td, T0, Ts, rhos, mus, Prs, kappas);
+    // droplet Re calculated using carrier's rho and mu
+    scalar Red = this->Re(td.rhoc(), U0, td.Uc(), d0, mug);
+    // droplet Re calculated using surface value
+    scalar Res = this->Re(rhos, U0, td.Uc(), d0, mus);
+
     // Calc mass and enthalpy transfer due to phase change
     calcPhaseChange
     (
         cloud,
         td,
         dt,
-        Res,
+        Red,  //using new droplet Re
         Prs,
         Ts,
         mus/rhos,
@@ -524,16 +553,28 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     // Correct surface values due to emitted species
     correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
-    Res = this->Re(rhos, U0, td.Uc(), this->d(), mus);
+    Res = this->Re(rhos, U0, td.Uc(), this->d_, mus);
+    Red = this->Re(td.rhoc(), U0, td.Uc(), this->d_, mug);
+
+    // Motion
+    // ~~~~~~
+
+    // 3. Calculate new particle velocity
+    this->U_ =
+        this->calcVelocity(cloud, td, dt, Res, mus, mass1, Su, dUTrans, Spu);
+
+    this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
 
+    Red = this->Re(td.rhoc(), this->U_, td.Uc(), this->d_, mug);
+    Res = this->Re(rhos, this->U_, td.Uc(), this->d_, mus);
 
     // 3. Compute heat- and momentum transfers
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Heat transfer
     // ~~~~~~~~~~~~~
 
-    // Calculate new particle temperature
+    // 4. Calculate new particle temperature
     this->T_ =
         this->calcHeatTransfer
         (
@@ -549,18 +590,10 @@ void Foam::ReactingParcel<ParcelType>::calc
             Sph
         );
 
-    this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
-
-
-    // Motion
-    // ~~~~~~
-
-    // Calculate new particle velocity
-    this->U_ =
-        this->calcVelocity(cloud, td, dt, Res, mus, mass1, Su, dUTrans, Spu);
+    this->Cp_ = composition.Cp(0, Y_, td.pc(), this->T_);
 
 
-    // 4. Accumulate carrier phase source terms
+    // 5. Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     if (cloud.solution().coupled())
@@ -584,6 +617,8 @@ void Foam::ReactingParcel<ParcelType>::calc
 
         // Update sensible enthalpy transfer
         cloud.hsTrans()[this->cell()] += np0*dhsTrans;
+        // add additional term for work between droplets and carrier
+        cloud.hsTrans()[this->cell()] += U0&cloud.UTrans()[this->cell()]*dt;
         cloud.hsCoeff()[this->cell()] += np0*Sph;
 
         // Update radiation fields

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

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2019 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 "LiquidEvaporationSpalding.H"
+// #include "specie.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant::mathematical;
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::LiquidEvaporationSpalding<CloudType>::LiquidEvaporationSpalding
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    LiquidEvaporationBoil<CloudType>(dict, owner)
+{
+
+}
+
+
+template<class CloudType>
+Foam::LiquidEvaporationSpalding<CloudType>::LiquidEvaporationSpalding
+(
+    const LiquidEvaporationSpalding<CloudType>& pcm
+)
+:
+    LiquidEvaporationBoil<CloudType>(pcm)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::LiquidEvaporationSpalding<CloudType>::~LiquidEvaporationSpalding()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::LiquidEvaporationSpalding<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
+{
+    // immediately evaporate mass that has reached critical condition
+    // Info <<"Using new evaporation model"<< endl;
+    if ((this->liquids_.Tc(X) - T) < small)
+    {
+        if (debug)
+        {
+            WarningInFunction
+                << "Parcel reached critical conditions: "
+                << "evaporating all available mass" << endl;
+        }
+
+        forAll(this->activeLiquids_, i)
+        {
+            const label lid = this->liqToLiqMap_[i];
+            dMassPC[lid] = great;
+        }
+
+        return;
+    }
+
+    // droplet surface pressure assumed to surface vapour pressure
+    scalar ps = this->liquids_.pv(pc, Ts, X);
+
+    // vapour density at droplet surface [kg/m3]
+    const scalar RR = 1000.0*constant::physicoChemical::R.value(); // J/(kmol·k)
+    scalar rhos = ps*this->liquids_.W(X)/(RR*Ts);
+
+    // calculate mass transfer of each specie in liquid
+    forAll(this->activeLiquids_, i)
+    {
+        const label gid = this->liqToCarrierMap_[i];
+        const label lid = this->liqToLiqMap_[i];
+
+        // boiling temperature at cell pressure for liquid species lid [K]
+        const scalar TBoil = this->liquids_.properties()[lid].pvInvert(pc);
+        // limit droplet temperature to boiling/critical temperature
+        const scalar Td = min(T, 0.999*TBoil);
+        // saturation pressure for liquid species lid [Pa]
+        const scalar pSat = this->liquids_.properties()[lid].pv(pc, Td);
+        // surface molar fraction - Raoult's Law
+        const scalar Xs = X[lid]*pSat/pc;
+
+        // vapour diffusivity [m2/s]
+        const scalar Dab = this->liquids_.properties()[lid].D(ps, Ts);
+        // Schmidt number
+        const scalar Sc = nu/(Dab + rootVSmall);
+        // Sherwood number
+        const scalar Sh = this->Sh(Re, Sc);
+
+        // mixture molar weight
+        // scalar W_gas =this->owner().thermo().carrier().W();
+        scalar W_gas =this->owner().thermo().carrier().W(); // kg/kmol
+        // scalar W_gas =this->owner().thermo().carrier().cellMixture(celli).W(); // kg/kmol
+
+        const scalar Yc = this->owner().thermo().carrier().Y()[gid][celli];
+
+        // fuel molar weight [kg/kmol]
+        const scalar W_fuel = this->liquids_.properties()[lid].W();
+
+        const scalar Ys = Xs*W_fuel/W_gas/(1. + Xs*W_fuel/W_gas - Xs); // ref Hu.2017, Spalding.1953
+        const scalar Bm = (Ys - Yc)/max(small, 1. - Ys);
+
+        if (Bm > 0)
+        {
+            // mass transfer [kg]
+            dMassPC[lid] += pi*d*Sh*Dab*rhos*log(1.0 + Bm)*dt;
+        }
+    }
+}
+
+// ************************************************************************* //