Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
mattijs
2013-08-13 12:09:54 +01:00
18 changed files with 355 additions and 78 deletions

View File

@ -79,8 +79,13 @@ void writeWeights(const polyMesh& mesh)
if (cpp.owner())
{
Info<< "Calculating AMI weights between owner patch: "
<< cpp.name() << " and neighbour patch: "
<< cpp.neighbPatch().name() << endl;
const AMIPatchToPatchInterpolation& ami =
cpp.AMI();
writeWeights
(
ami.tgtWeightsSum(),

View File

@ -28,7 +28,23 @@ Group
grpRASTurbulence
Description
Standard k-epsilon turbulence model
Mixture k-epsilon turbulence model for two-phase gas-liquid systems
The basic structure of the model is based on:
\verbatim
"Modelling of dispersed bubble and droplet ow at high phase fractions"
A. Behzadi, R.I. Issa , H. Rusche,
Chemical Engineering Science (59) 2004 pp.759-770.
\endverbatim
but an effective density for the gas is used in the averaging and an
alternative model for bubble-generated turbulence from:
\verbatim
"The simulation of multidimensional multiphase flows",
Lahey R.T.,
Nucl. Eng. & Design
(235) 2005 pp.1043-1060.
\endverbatim
The default model coefficients correspond to the following:
\verbatim
@ -37,7 +53,7 @@ Description
Cmu 0.09;
C1 1.44;
C2 1.92;
C3 -0.33;
C3 C2;
sigmak 1.0;
sigmaEps 1.3;
}

View File

@ -511,7 +511,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
if
(
!td.cloud().devolatilisation().active()
|| T < td.cloud().constProps().Tvap()
|| T < td.cloud().constProps().TDevol()
|| canCombust == -1
)
{

View File

@ -81,6 +81,9 @@ public:
{
// Private data
//- Devolatilisation activation temperature [K]
scalar TDevol_;
//- Latent heat of devolatilisation [J/kg]
scalar LDevol_;
@ -106,9 +109,32 @@ public:
const bool readFields = true
);
//- Construct from components
constantProperties
(
const label parcelTypeId,
const scalar rhoMin,
const scalar rho0,
const scalar minParticleMass,
const scalar youngsModulus,
const scalar poissonsRatio,
const scalar T0,
const scalar TMin,
const scalar TMax,
const scalar Cp0,
const scalar epsilon0,
const scalar f0,
const scalar Pr,
const scalar pMin,
const Switch& constantVolume,
const scalar TDevol
);
// Access
//- Return const access to the devolatilisation temperature
inline scalar TDevol() const;
//- Return const access to the latent heat of devolatilisation
inline scalar LDevol() const;

View File

@ -30,6 +30,7 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
constantProperties()
:
ParcelType::constantProperties(),
TDevol_(0.0),
LDevol_(0.0),
hRetentionCoeff_(0.0)
{}
@ -43,6 +44,7 @@ constantProperties
)
:
ParcelType::constantProperties(cp),
TDevol_(cp.TDevol_),
LDevol_(cp.LDevol_),
hRetentionCoeff_(cp.hRetentionCoeff_)
{}
@ -57,11 +59,13 @@ constantProperties
)
:
ParcelType::constantProperties(parentDict, readFields),
TDevol_(0.0),
LDevol_(0.0),
hRetentionCoeff_(0.0)
{
if (readFields)
{
this->dict().lookup("TDevol") >> TDevol_;
this->dict().lookup("LDevol") >> LDevol_;
this->dict().lookup("hRetentionCoeff") >> hRetentionCoeff_;
@ -80,6 +84,50 @@ constantProperties
}
template<class ParcelType>
inline Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
constantProperties
(
const label parcelTypeId,
const scalar rhoMin,
const scalar rho0,
const scalar minParticleMass,
const scalar youngsModulus,
const scalar poissonsRatio,
const scalar T0,
const scalar TMin,
const scalar TMax,
const scalar Cp0,
const scalar epsilon0,
const scalar f0,
const scalar Pr,
const scalar pMin,
const Switch& constantVolume,
const scalar TDevol
)
:
ParcelType::constantProperties
(
parcelTypeId,
rhoMin,
rho0,
minParticleMass,
youngsModulus,
poissonsRatio,
T0,
TMin,
TMax,
Cp0,
epsilon0,
f0,
Pr,
pMin,
constantVolume
),
TDevol_(TDevol)
{}
template<class ParcelType>
inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
(
@ -148,6 +196,14 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
// * * * * * * * * * constantProperties Member Functions * * * * * * * * * * //
template<class ParcelType>
inline Foam::scalar
Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::TDevol() const
{
return TDevol_;
}
template<class ParcelType>
inline Foam::scalar
Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::LDevol() const

View File

@ -26,6 +26,7 @@ License
#include "ReactingParcel.H"
#include "specie.H"
#include "CompositionModel.H"
#include "PhaseChangeModel.H"
#include "mathematicalConstants.H"
using namespace Foam::constant::mathematical;
@ -56,26 +57,22 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
scalarField& Cs
)
{
if
(
!td.cloud().phaseChange().active()
|| T < td.cloud().constProps().Tvap()
|| YPhase < SMALL
)
typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
PhaseChangeModel<reactingCloudType>& phaseChange = td.cloud().phaseChange();
scalar Tvap = phaseChange.Tvap(YComponents);
if (!phaseChange.active() || T < Tvap || YPhase < SMALL)
{
return;
}
typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
const CompositionModel<reactingCloudType>& composition =
td.cloud().composition();
const scalar TMax = td.cloud().phaseChange().TMax(pc_);
const scalar TMax = phaseChange.TMax(pc_, YComponents);
const scalar Tdash = min(T, TMax);
const scalar Tsdash = min(Ts, TMax);
// Calculate mass transfer due to phase change
td.cloud().phaseChange().calculate
phaseChange.calculate
(
dt,
cellI,
@ -97,14 +94,17 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
const scalar dMassTot = sum(dMassPC);
// Add to cumulative phase change mass
td.cloud().phaseChange().addToPhaseChangeMass(this->nParticle_*dMassTot);
phaseChange.addToPhaseChangeMass(this->nParticle_*dMassTot);
const CompositionModel<reactingCloudType>& composition =
td.cloud().composition();
forAll(dMassPC, i)
{
const label idc = composition.localToGlobalCarrierId(idPhase, i);
const label idl = composition.globalIds(idPhase)[i];
const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, Tdash);
const scalar dh = phaseChange.dh(idc, idl, pc_, Tdash);
Sh -= dMassPC[i]*dh/dt;
}

View File

@ -80,9 +80,6 @@ public:
//- Constant volume flag - e.g. during mass transfer
Switch constantVolume_;
//- Vaporisation temperature [K]
scalar Tvap_;
public:
@ -118,8 +115,7 @@ public:
const scalar f0,
const scalar Pr,
const scalar pMin,
const Switch& constantVolume,
const scalar Tvap
const Switch& constantVolume
);
@ -130,9 +126,6 @@ public:
//- Return const access to the constant volume flag
inline Switch constantVolume() const;
//- Return const access to the vaporisation temperature
inline scalar Tvap() const;
};

View File

@ -31,8 +31,7 @@ Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties()
:
ParcelType::constantProperties(),
pMin_(0.0),
constantVolume_(false),
Tvap_(0.0)
constantVolume_(false)
{}
@ -44,8 +43,7 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
:
ParcelType::constantProperties(cp),
pMin_(cp.pMin_),
constantVolume_(cp.constantVolume_),
Tvap_(cp.Tvap_)
constantVolume_(cp.constantVolume_)
{}
@ -58,8 +56,7 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
:
ParcelType::constantProperties(parentDict, readFields),
pMin_(1000.0),
constantVolume_(false),
Tvap_(0.0)
constantVolume_(false)
{
if (readFields)
{
@ -69,7 +66,6 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
}
this->dict().lookup("constantVolume") >> constantVolume_;
this->dict().lookup("Tvap") >> Tvap_;
}
}
@ -91,8 +87,7 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
const scalar f0,
const scalar Pr,
const scalar pMin,
const Switch& constantVolume,
const scalar Tvap
const Switch& constantVolume
)
:
ParcelType::constantProperties
@ -112,8 +107,7 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
Pr
),
pMin_(pMin),
constantVolume_(constantVolume),
Tvap_(Tvap)
constantVolume_(constantVolume)
{}
@ -198,14 +192,6 @@ Foam::ReactingParcel<ParcelType>::constantProperties::constantVolume() const
}
template<class ParcelType>
inline Foam::scalar
Foam::ReactingParcel<ParcelType>::constantProperties::Tvap() const
{
return Tvap_;
}
// * * * * * * * * * * ThermoParcel Member Functions * * * * * * * * * * * * //
template<class ParcelType>

View File

@ -151,6 +151,44 @@ void Foam::LiquidEvaporation<CloudType>::calculate
scalarField& dMassPC
) const
{
// liquid volume fraction
const scalarField X(liquids_.X(Yl));
// immediately evaporate mass that has reached critical condition
if (mag(T - liquids_.Tc(X)) < SMALL)
{
if (debug)
{
WarningIn
(
"void Foam::LiquidEvaporation<CloudType>::calculate"
"("
"const scalar, "
"const label, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalarField&, "
"scalarField&"
") const"
) << "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));
@ -242,9 +280,27 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
template<class CloudType>
Foam::scalar Foam::LiquidEvaporation<CloudType>::TMax(const scalar pIn) const
Foam::scalar Foam::LiquidEvaporation<CloudType>::Tvap
(
const scalarField& Y
) const
{
return liquids_.TMax(pIn);
const scalarField X(liquids_.X(Y));
return liquids_.Tpt(X);
}
template<class CloudType>
Foam::scalar Foam::LiquidEvaporation<CloudType>::TMax
(
const scalar p,
const scalarField& Y
) const
{
const scalarField X(liquids_.X(Y));
return liquids_.pvInvert(p, X);
}

View File

@ -131,8 +131,11 @@ public:
const scalar T
) const;
//- Return vapourisation temperature
virtual scalar Tvap(const scalarField& Y) const;
//- Return maximum/limiting temperature
virtual scalar TMax(const scalar pIn) const;
virtual scalar TMax(const scalar p, const scalarField& Y) const;
};

View File

@ -151,18 +151,53 @@ void Foam::LiquidEvaporationBoil<CloudType>::calculate
scalarField& dMassPC
) const
{
// construct carrier phase species volume fractions for cell, cellI
const scalarField XcMix(calcXc(cellI));
// liquid volume fraction
const scalarField X(liquids_.X(Yl));
// immediately evaporate mass that has reached critical condition
if (mag(T - liquids_.Tc(X)) < SMALL)
{
if (debug)
{
WarningIn
(
"void Foam::LiquidEvaporationBoil<CloudType>::calculate"
"("
"const scalar, "
"const label, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalarField&, "
"scalarField&"
") const"
) << "Parcel reached critical conditions: "
<< "evaporating all avaliable mass" << endl;
}
forAll(activeLiquids_, i)
{
const label lid = liqToLiqMap_[i];
dMassPC[lid] = GREAT;
}
return;
}
// droplet surface pressure assumed to surface vapour pressure
scalar ps = liquids_.pv(pc, Ts, X);
// vapour density at droplet surface [kg/m3]
scalar rhos = ps*liquids_.W(X)/(specie::RR*Ts);
// construct carrier phase species volume fractions for cell, cellI
const scalarField XcMix(calcXc(cellI));
// carrier thermo properties
scalar Hsc = 0.0;
scalar Hc = 0.0;
@ -341,12 +376,27 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
template<class CloudType>
Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::TMax
Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::Tvap
(
const scalar pIn
const scalarField& Y
) const
{
return liquids_.TMax(pIn);
const scalarField X(liquids_.X(Y));
return liquids_.Tpt(X);
}
template<class CloudType>
Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::TMax
(
const scalar p,
const scalarField& Y
) const
{
const scalarField X(liquids_.X(Y));
return liquids_.pvInvert(p, X);
}

View File

@ -141,8 +141,11 @@ public:
const scalar T
) const;
//- Return vapourisation temperature
virtual scalar Tvap(const scalarField& Y) const;
//- Return maximum/limiting temperature
virtual scalar TMax(const scalar pIn) const;
virtual scalar TMax(const scalar p, const scalarField& Y) const;
};

View File

@ -179,12 +179,23 @@ Foam::scalar Foam::PhaseChangeModel<CloudType>::dh
template<class CloudType>
Foam::scalar Foam::PhaseChangeModel<CloudType>::TMax(const scalar) const
Foam::scalar Foam::PhaseChangeModel<CloudType>::TMax
(
const scalar,
const scalarField&
) const
{
return GREAT;
}
template<class CloudType>
Foam::scalar Foam::PhaseChangeModel<CloudType>::Tvap(const scalarField& Y) const
{
return -GREAT;
}
template<class CloudType>
void Foam::PhaseChangeModel<CloudType>::addToPhaseChangeMass(const scalar dMass)
{

View File

@ -184,8 +184,11 @@ public:
const scalar T
) const;
//- Return vapourisation temperature
virtual scalar Tvap(const scalarField& Y) const;
//- Return maximum/limiting temperature
virtual scalar TMax(const scalar pIn) const;
virtual scalar TMax(const scalar p, const scalarField& Y) const;
//- Add to phase change mass
void addToPhaseChangeMass(const scalar dMass);

View File

@ -64,47 +64,62 @@ void Foam::SprayParcel<ParcelType>::calc
const label cellI
)
{
typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
const CompositionModel<reactingCloudType>& composition =
td.cloud().composition();
bool coupled = td.cloud().solution().coupled();
// check if parcel belongs to liquid core
if (liquidCore() > 0.5)
{
// liquid core parcels should not interact with the gas
if (td.cloud().solution().coupled())
{
td.cloud().solution().coupled() = false;
}
td.cloud().solution().coupled() = false;
}
// get old mixture composition
const scalarField& Y0(this->Y());
scalarField X0(composition.liquids().X(Y0));
// check if we have critical or boiling conditions
scalar TMax = composition.liquids().Tc(X0);
const scalar T0 = this->T();
const scalar pc0 = this->pc_;
if (composition.liquids().pv(pc0, T0, X0) >= pc0*0.999)
{
// set TMax to boiling temperature
TMax = composition.liquids().pvInvert(pc0, X0);
}
// set the maximum temperature limit
const scalar TMax = td.cloud().composition().liquids().TMax(this->pc_);
td.cloud().constProps().TMax() = TMax;
// store the parcel properties
const scalarField& Y(this->Y());
scalarField X(td.cloud().composition().liquids().X(Y));
scalarField X(composition.liquids().X(Y));
scalar T0 = this->T();
this->Cp() = td.cloud().composition().liquids().Cp(this->pc_, T0, X);
sigma_ = td.cloud().composition().liquids().sigma(this->pc_, T0, X);
scalar rho0 = td.cloud().composition().liquids().rho(this->pc_, T0, X);
this->Cp() = composition.liquids().Cp(this->pc_, T0, X);
sigma_ = composition.liquids().sigma(this->pc_, T0, X);
scalar rho0 = composition.liquids().rho(this->pc_, T0, X);
this->rho() = rho0;
ParcelType::calc(td, dt, cellI);
if (td.keepParticle)
{
// update Cp, sigma, diameter and density due to change in temperature
// update Cp, sigma, density and diameter due to change in temperature
// and/or composition
scalar T1 = this->T();
const scalarField& Y1(this->Y());
scalarField X1(td.cloud().composition().liquids().X(Y1));
scalarField X1(composition.liquids().X(Y1));
this->Cp() = td.cloud().composition().liquids().Cp(this->pc_, T1, X1);
sigma_ = td.cloud().composition().liquids().sigma(this->pc_, T1, X);
this->Cp() = composition.liquids().Cp(this->pc_, T1, X1);
scalar rho1 = td.cloud().composition().liquids().rho(this->pc_, T1, X1);
sigma_ = composition.liquids().sigma(this->pc_, T1, X);
scalar rho1 = composition.liquids().rho(this->pc_, T1, X1);
this->rho() = rho1;
scalar d1 = this->d()*cbrt(rho0/rho1);
this->d() = d1;

View File

@ -103,12 +103,62 @@ Foam::scalar Foam::liquidMixtureProperties::Tc(const scalarField& x) const
}
Foam::scalar Foam::liquidMixtureProperties::TMax(const scalar p) const
Foam::scalar Foam::liquidMixtureProperties::Tpt(const scalarField& x) const
{
scalar T = -GREAT;
scalar Tpt = 0.0;
forAll(properties_, i)
{
T = max(T, properties_[i].pvInvert(p));
Tpt += x[i]*properties_[i].Tt();
}
return Tpt;
}
Foam::scalar Foam::liquidMixtureProperties::pvInvert
(
const scalar p,
const scalarField& x
) const
{
// Set upper and lower bounds
scalar Thi = Tc(x);
scalar Tlo = Tpt(x);
// Check for critical and solid phase conditions
if (p >= pv(p, Thi, x))
{
return Thi;
}
else if (p < pv(p, Tlo, x))
{
WarningIn
(
"Foam::scalar Foam::liquidMixtureProperties::pvInvert"
"("
" const scalar,"
" const scalarField&"
") const"
) << "Pressure below triple point pressure: "
<< "p = " << p << " < Pt = " << pv(p, Tlo, x) << nl << endl;
return -1;
}
// Set initial guess
scalar T = (Thi + Tlo)*0.5;
while ((Thi - Tlo) > 1.0e-4)
{
if ((pv(p, T, x) - p) <= 0.0)
{
Tlo = T;
}
else
{
Thi = T;
}
T = (Thi + Tlo)*0.5;
}
return T;

View File

@ -141,8 +141,9 @@ public:
//- Calculate the critical temperature of mixture
scalar Tc(const scalarField& x) const;
//- Calculate the maximum temperature of mixture at given pressure
scalar TMax(const scalar p) const;
//- Invert the vapour pressure relationship to retrieve the boiling
// temperature of the mixture as a function of pressure
scalar pvInvert(const scalar p, const scalarField& x) const;
//- Return pseudocritical temperature according to Kay's rule
scalar Tpc(const scalarField& x) const;
@ -150,6 +151,9 @@ public:
//- Return pseudocritical pressure (modified Prausnitz and Gunn)
scalar Ppc(const scalarField& x) const;
//- Return pseudo triple point temperature (mole averaged formulation)
scalar Tpt(const scalarField& x) const;
//- Return mixture accentric factor
scalar omega(const scalarField& x) const;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -252,7 +252,7 @@ public:
virtual scalar D(scalar p, scalar T, scalar Wb) const;
//- Invert the vapour pressure relationship to retrieve the
// boiling temperuture as a function of temperature
// boiling temperuture as a function of pressure
virtual scalar pvInvert(scalar p) const;