Commit 11d17fec authored by sergio's avatar sergio Committed by Andrew Heather
Browse files

ENH: Adding evaporation-condensation lagragian model for solution

1) Adding LiquidEvapFuchsKnudsen model for lagrangian evaporation.
   This models is based on a diffusion type of evaporation/
   condensation on particles composed of solution (liquid + solid).

2) Adding modes of calculating the particle rho and volume change.
   The new keyword in constantProperties is 'volumeUpdateMethod'
   which three options:
        a) constantRho
        b) constantVolume
        c) updateRhoAndVol

   The old keyword 'constantVolume' true/face is still valid

3) The entry rho0 is now optional for multicomponent parcels.
   If defined , it is used, but if it is not the actual mixture
   provided is used to calculate rho0 of the particle.
   T0 is still used as initial T and Cp0 is over-written in the
   multicomponent cloud but still required.

4) Adding tutorial for evaporation/condensation model
parent 973e2d4e
......@@ -562,7 +562,11 @@ void Foam::KinematicCloud<CloudType>::setParcelThermoProperties
const scalar lagrangianDt
)
{
parcel.rho() = constProps_.rho0();
// If rho0 is given in the const properties
if (constProps_.rho0() != -1)
{
parcel.rho() = constProps_.rho0();
}
}
......@@ -581,6 +585,14 @@ void Foam::KinematicCloud<CloudType>::checkParcelProperties
{
parcel.typeId() = constProps_.parcelTypeId();
}
if (parcel.rho() == -1)
{
FatalErrorInFunction
<< "The kinematic cloud needs rho0 in the constantProperties "
<< " dictionary. " << nl
<< abort(FatalError);
}
}
......
......@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2018-2019 OpenCFD Ltd.
Copyright (C) 2018-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -136,6 +136,25 @@ void Foam::ReactingHeterogeneousCloud<CloudType>::setParcelThermoProperties
// Set the parcel to combust
parcel.canCombust() = 1;
// If rho0 was given in constProp use it. If not use the composition
// to set tho
if (this->constProps_.rho0() == -1)
{
const label idGas = this->composition().idGas();
const label idLiquid = this->composition().idLiquid();
const label idSolid = this->composition().idSolid();
const scalarField& Ygas = this->composition().Y0(idGas);
const scalarField& Yliq = this->composition().Y0(idLiquid);
const scalarField& Ysol = this->composition().Y0(idSolid);
const scalar p0 =
this->composition().thermo().thermo().p()[parcel.cell()];
const scalar T0 = this->constProps_.T0();
parcel.rho() = this->composition().rho(Ygas, Yliq, Ysol, T0, p0);
}
}
......
......@@ -174,6 +174,21 @@ void Foam::ReactingMultiphaseCloud<CloudType>::setParcelThermoProperties
parcel.YGas() = this->composition().Y0(idGas);
parcel.YLiquid() = this->composition().Y0(idLiquid);
parcel.YSolid() = this->composition().Y0(idSolid);
// If rho0 was given in constProp use it. If not use the composition
// to set tho
if (constProps_.rho0() == -1)
{
const scalarField& Ygas = this->composition().Y0(idGas);
const scalarField& Yliq = this->composition().Y0(idLiquid);
const scalarField& Ysol = this->composition().Y0(idSolid);
const scalar p0 =
this->composition().thermo().thermo().p()[parcel.cell()];
const scalar T0 = constProps_.T0();
parcel.rho() = this->composition().rho(Ygas, Yliq, Ysol, T0, p0);
}
}
......
......@@ -67,7 +67,7 @@ inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
dict_(parentDict.subOrEmptyDict("constantProperties")),
parcelTypeId_(dict_, "parcelTypeId", -1),
rhoMin_(dict_, "rhoMin", 1e-15),
rho0_(dict_, "rho0"),
rho0_(dict_, "rho0", -1),
minParcelMass_(dict_, "minParcelMass", 1e-15)
{}
......
......@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2018-2019 OpenCFD Ltd.
Copyright (C) 2018-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -80,6 +80,30 @@ Foam::scalar Foam::ReactingHeterogeneousParcel<ParcelType>::LEff
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
template<class ParcelType>
template<class TrackCloudType>
Foam::scalar Foam::ReactingHeterogeneousParcel<ParcelType>::updatedDeltaVolume
(
TrackCloudType& cloud,
const scalarField& dMass,
const scalar p,
const scalar T
)
{
const auto& composition = cloud.composition();
scalarField dVolSolid(dMass.size(), Zero);
forAll(dVolSolid, i)
{
dVolSolid[i] =
dMass[i]/composition.solids().properties()[i].rho();
}
return sum(dVolSolid);
}
template<class ParcelType>
template<class TrackCloudType>
void Foam::ReactingHeterogeneousParcel<ParcelType>::calc
......@@ -191,16 +215,52 @@ void Foam::ReactingHeterogeneousParcel<ParcelType>::calc
(void)this->updateMassFraction(mass0, dMassSolid, this->Y_);
// Update particle density or diameter
if (cloud.constProps().constantVolume())
if
(
cloud.constProps().volUpdateType()
== constantProperties::volumeUpadteType::mUndefined
)
{
this->rho_ = mass1/this->volume();
if (cloud.constProps().constantVolume())
{
this->rho_ = mass1/this->volume();
}
else
{
this->d_ = cbrt(mass1/this->rho_*6/pi);
}
}
else
{
this->d_ = cbrt(mass1/this->rho_*6.0/pi);
switch (cloud.constProps().volUpdateType())
{
case constantProperties::volumeUpadteType::mConstRho :
{
this->d_ = cbrt(mass1/this->rho_*6/pi);
break;
}
case constantProperties::volumeUpadteType::mConstVol :
{
this->rho_ = mass1/this->volume();
break;
}
case constantProperties::volumeUpadteType::mUpdateRhoAndVol :
{
scalar deltaVol =
updatedDeltaVolume
(
cloud,
dMassSolid,
pc,
T0
);
this->rho_ = mass1/(this->volume() + deltaVol);
this->d_ = cbrt(mass1/this->rho_*6.0/pi);
break;
}
}
}
// Correct surface values due to emitted species
this->correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
Res = this->Re(rhos, U0, td.Uc(), this->d_, mus);
......
......@@ -168,6 +168,16 @@ protected:
// Protected Member Functions
//- Return change of volume due to mass exchange
template<class TrackCloudType>
scalar updatedDeltaVolume
(
TrackCloudType& cloud,
const scalarField& dMass,
const scalar p,
const scalar T
);
//- Calculate surface reactions
template<class TrackCloudType>
......
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -139,6 +139,48 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
}
template<class ParcelType>
template<class TrackCloudType>
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updatedDeltaVolume
(
TrackCloudType& cloud,
const scalarField& dMassGas,
const scalarField& dMassLiquid,
const scalarField& dMassSolid,
const label idG,
const label idL,
const label idS,
const scalar p,
const scalar T
)
{
const auto& props = cloud.composition().phaseProps()[idG];
const auto& thermo = cloud.composition().thermo();
scalarField dVolGas(dMassGas.size(), Zero);
forAll(dMassGas, i)
{
label cid = props.carrierIds()[i];
dVolGas[i] = -dMassGas[i]/thermo.carrier().rho(cid, p, T);
}
scalarField dVolLiquid(dMassLiquid.size(), Zero);
forAll(dMassLiquid, i)
{
dVolLiquid[i] =
-dMassLiquid[i]/thermo.liquids().properties()[i].rho(p, T);
}
scalarField dVolSolid(dMassSolid.size(), Zero);
forAll(dMassSolid, i)
{
dVolSolid[i] = -dMassSolid[i]/thermo.solids().properties()[i].rho();
}
return (sum(dVolGas) + sum(dVolLiquid) + sum(dMassSolid));
}
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
template<class ParcelType>
......@@ -188,6 +230,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const vector& U0 = this->U_;
const scalar T0 = this->T_;
const scalar mass0 = this->mass();
const scalar rho0 = this->rho_;
const scalar pc = td.pc();
......@@ -256,9 +300,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
d0,
T0,
mass0,
rho0,
idL,
YMix[LIQ],
YLiquid_,
YMix[SLD]*YSolid_,
dMassPC,
Sh,
Ne,
......@@ -382,16 +428,57 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
(void)updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
// Update particle density or diameter
if (cloud.constProps().constantVolume())
if
(
cloud.constProps().volUpdateType()
== constantProperties::volumeUpadteType::mUndefined
)
{
this->rho_ = mass1/this->volume();
if (cloud.constProps().constantVolume())
{
this->rho_ = mass1/this->volume();
}
else
{
this->d_ = cbrt(mass1/this->rho_*6/pi);
}
}
else
{
this->d_ = cbrt(mass1/this->rho_*6.0/pi);
switch (cloud.constProps().volUpdateType())
{
case constantProperties::volumeUpadteType::mConstRho :
{
this->d_ = cbrt(mass1/this->rho_*6/pi);
break;
}
case constantProperties::volumeUpadteType::mConstVol :
{
this->rho_ = mass1/this->volume();
break;
}
case constantProperties::volumeUpadteType::mUpdateRhoAndVol :
{
scalar deltaVol =
updatedDeltaVolume
(
cloud,
dMassGas,
dMassLiquid,
dMassSolid,
idG,
idL,
idS,
pc,
T0
);
this->rho_ = mass1/(this->volume() + deltaVol);
this->d_ = cbrt(mass1/this->rho_*6/pi);
break;
}
}
}
// Correct surface values due to emitted species
this->correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
Res = this->Re(rhos, U0, td.Uc(), this->d_, mus);
......
......@@ -209,6 +209,22 @@ protected:
// Protected Member Functions
//- Return change of volume due to mass exchange
template<class TrackCloudType>
scalar updatedDeltaVolume
(
TrackCloudType& cloud,
const scalarField& dMassGas,
const scalarField& dMassLiquid,
const scalarField& dMassSolid,
const label idG,
const label idL,
const label idS,
const scalar p,
const scalar T
);
//- Calculate Devolatilisation
template<class TrackCloudType>
void calcDevolatilisation
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -35,6 +36,29 @@ using namespace Foam::constant::mathematical;
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
template<class ParcelType>
template<class TrackCloudType>
Foam::scalar Foam::ReactingParcel<ParcelType>::updatedDeltaVolume
(
TrackCloudType& cloud,
const scalarField& dMass,
const scalar p,
const scalar T
)
{
const auto& composition = cloud.composition();
scalarField dVolLiquid(dMass.size(), Zero);
forAll(dVolLiquid, i)
{
dVolLiquid[i] =
dMass[i]/composition.liquids().properties()[i].rho(p, T);
}
return sum(dVolLiquid);
}
template<class ParcelType>
template<class TrackCloudType>
void Foam::ReactingParcel<ParcelType>::calcPhaseChange
......@@ -49,9 +73,11 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
const scalar d,
const scalar T,
const scalar mass,
const scalar rho,
const label idPhase,
const scalar YPhase,
const scalarField& Y,
const scalarField& Ysol,
scalarField& dMassPC,
scalar& Sh,
scalar& N,
......@@ -62,7 +88,8 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
const auto& composition = cloud.composition();
auto& phaseChange = cloud.phaseChange();
if (!phaseChange.active() || (YPhase < SMALL))
// Some models allow evaporation and condensation
if (!phaseChange.active())
{
return;
}
......@@ -70,7 +97,6 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
scalarField X(composition.liquids().X(Y));
scalar Tvap = phaseChange.Tvap(X);
if (T < Tvap)
{
return;
......@@ -89,16 +115,27 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
Pr,
d,
nus,
rho,
Tdash,
Tsdash,
td.pc(),
td.Tc(),
X,
X, // components molar fractions purely in the liquid
Ysol*mass, // total solid mass
YPhase*Y*mass, // total liquid mass
dMassPC
);
// Limit phase change mass by availability of each specie
dMassPC = min(mass*YPhase*Y, dMassPC);
forAll (Y, i)
{
// evaporation
if (dMassPC[i] > 0)
{
dMassPC[i] = min(mass*YPhase*Y[i], dMassPC[i]);
}
// condensation Do something?
}
const scalar dMassTot = sum(dMassPC);
......@@ -394,6 +431,7 @@ void Foam::ReactingParcel<ParcelType>::calc
const vector& U0 = this->U_;
const scalar T0 = this->T_;
const scalar mass0 = this->mass();
const scalar rho0 = this->rho_;
// Calc surface values
......@@ -455,9 +493,11 @@ void Foam::ReactingParcel<ParcelType>::calc
d0,
T0,
mass0,
rho0,
0,
1.0,
Y_,
scalarField(0),
dMassPC,
Sh,
Ne,
......@@ -474,14 +514,52 @@ void Foam::ReactingParcel<ParcelType>::calc
this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
// Update particle density or diameter
if (cloud.constProps().constantVolume())
if
(
cloud.constProps().volUpdateType()
== constantProperties::volumeUpadteType::mUndefined
)
{
this->rho_ = mass1/this->volume();
// Update particle density or diameter
if (cloud.constProps().constantVolume())
{
this->rho_ = mass1/this->volume();
}
else
{
this->d_ = cbrt(mass1/this->rho_*6/pi);
}
}
else
{
this->d_ = cbrt(mass1/this->rho_*6.0/pi);
switch (cloud.constProps().volUpdateType())
{
case constantProperties::volumeUpadteType::mConstRho :
{
this->d_ = cbrt(mass1/this->rho_*6/pi);
break;
}
case constantProperties::volumeUpadteType::mConstVol :
{
this->rho_ = mass1/this->volume();
break;
}
case constantProperties::volumeUpadteType::mUpdateRhoAndVol :
{
scalar deltaVol =
updatedDeltaVolume
(
cloud,
dMass,
td.pc(),
T0
);
this->rho_ = mass1/(this->volume() - deltaVol);
this->d_ = cbrt(mass1/this->rho_*6/pi);
break;
}
}
}
// Remove the particle when mass falls below minimum threshold
......
......@@ -79,19 +79,39 @@ public:
static const std::size_t sizeofFields;
//- Class to hold reacting parcel constant properties
class constantProperties
:
public ParcelType::constantProperties
{
// Private data
public:
//- Type of activity coefficient models
enum volumeUpadteType
{
mConstRho,
mConstVol,
mUpdateRhoAndVol,
mUndefined
};
private:
// Private data
//- Minimum pressure [Pa]