Commit 68af83db authored by Andrew Heather's avatar Andrew Heather
Browse files

ENH: Updates to surface film phase change modelling

parent ecbdccf2
......@@ -108,11 +108,11 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
const scalarField& YInf = film.YPrimary()[vapId];
const scalarField& pInf = film.pPrimary();
const scalarField& T = film.T();
const scalarField& Ts = film.Ts();
const scalarField& Tw = film.Tw();
const scalarField& TInf = film.TPrimary();
const scalarField& rho = film.rho();
const scalarField& mu = film.mu();
const scalarField& TInf = film.TPrimary();
const scalarField& rhoInf = film.rhoPrimary();
const scalarField& muInf = film.muPrimary();
const scalarField& magSf = film.magSf();
const scalarField hInf = film.htcs().h();
const scalarField hFilm = film.htcw().h();
......@@ -128,7 +128,7 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
const scalar pc = pInf[cellI];
// local temperature - impose lower limit of 200 K for stability
const scalar Tloc = min(liq.Tc(), max(200.0, Ts[cellI]));
const scalar Tloc = min(liq.Tc(), max(200.0, T[cellI]));
// saturation pressure [Pa]
const scalar pSat = liq.pv(pc, Tloc);
......@@ -137,7 +137,7 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
const scalar hVap = liq.hl(pc, Tloc);
// calculate mass transfer
if (pSat >= pc)
if (pSat >= 0.95*pc)
{
// boiling
const scalar qDotInf = hInf[cellI]*(TInf[cellI] - T[cellI]);
......@@ -152,8 +152,14 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
}
else
{
// Primary region density [kg/m3]
const scalar rhoInfc = rhoInf[cellI];
// Primary region viscosity [Pa.s]
const scalar muInfc = muInf[cellI];
// Reynolds number
const scalarField Re = rho*mag(dU)*L_/mu;
const scalar Re = rhoInfc*mag(dU[cellI])*L_/muInfc;
// molecular weight of vapour [kg/kmol]
const scalar Wvap = thermo.carrier().W(vapId);
......@@ -164,29 +170,21 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct
// vapour mass fraction at interface
const scalar Ys = Wliq*pSat/(Wliq*pSat + Wvap*(pc - pSat));
// bulk gas average density [kg/m3]
scalar Winf = 0;
forAll(thermo.carrier().Y(), i)
{
Winf += film.YPrimary()[i][cellI]*thermo.carrier().W(i);
}
const scalar rhoInf = Winf*pc/(specie::RR*Tloc);
// vapour diffusivity [m2/s]
const scalar Dab = liq.D(pc, Tloc);
// Schmidt number
const scalar Sc = mu[cellI]/(rho[cellI]*(Dab + ROOTVSMALL));
const scalar Sc = muInfc/(rhoInfc*(Dab + ROOTVSMALL));
// Sherwood number
const scalar Sh = this->Sh(Re[cellI], Sc);
const scalar Sh = this->Sh(Re, Sc);
// mass transfer coefficient [m/s]
const scalar hm = Sh*Dab/L_;
// add mass contribution to source
dMass[cellI] =
dt*magSf[cellI]*rhoInf*hm*(Ys - YInf[cellI])/(1.0 - Ys);
dt*magSf[cellI]*rhoInfc*hm*(Ys - YInf[cellI])/(1.0 - Ys);
}
dMass[cellI] = min(availableMass[cellI], max(0.0, dMass[cellI]));
......
......@@ -203,10 +203,12 @@ resetPrimaryRegionSourceTerms()
void Foam::surfaceFilmModels::kinematicSingleLayer::
transferPrimaryRegionFields()
{
// Update pressure and velocity from primary region via direct mapped
// Update fields from primary region via direct mapped
// (coupled) boundary conditions
UPrimary_.correctBoundaryConditions();
pPrimary_.correctBoundaryConditions();
rhoPrimary_.correctBoundaryConditions();
muPrimary_.correctBoundaryConditions();
// Retrieve the source fields from the primary region via direct mapped
// (coupled) boundary conditions
......@@ -930,6 +932,34 @@ Foam::surfaceFilmModels::kinematicSingleLayer::kinematicSingleLayer
),
filmRegion_
),
rhoPrimary_
(
IOobject
(
"rho", // must have same name as rho to enable mapping
time_.timeName(),
filmRegion_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
filmRegion_,
dimensionedScalar("zero", dimDensity, 0.0),
pPrimary_.boundaryField().types()
),
muPrimary_
(
IOobject
(
"mu", // must have same name as mu to enable mapping
time_.timeName(),
filmRegion_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
filmRegion_,
dimensionedScalar("zero", dimPressure*dimTime, 0.0),
pPrimary_.boundaryField().types()
),
injection_(injectionModel::New(*this, coeffs_)),
......@@ -1193,15 +1223,15 @@ void Foam::surfaceFilmModels::kinematicSingleLayer::info() const
// Output Courant number for info only - does not change time step
CourantNumber();
Info<< indent << "added mass = "
Info<< indent << "added mass = "
<< returnReduce<scalar>(addedMass_, sumOp<scalar>()) << nl
<< indent << "current mass = "
<< indent << "current mass = "
<< gSum((deltaRho_*magSf_)()) << nl
<< indent << "detached mass = "
<< indent << "detached mass = "
<< returnReduce<scalar>(detachedMass_, sumOp<scalar>()) << nl
<< indent << "min/max(mag(U)) = " << min(mag(U_)).value() << ", "
<< indent << "min/max(mag(U)) = " << min(mag(U_)).value() << ", "
<< max(mag(U_)).value() << nl
<< indent << "min/max(delta) = " << min(delta_).value() << ", "
<< indent << "min/max(delta) = " << min(delta_).value() << ", "
<< max(delta_).value() << nl;
}
......
......@@ -212,6 +212,12 @@ protected:
//- Pressure / [Pa]
volScalarField pPrimary_;
//- Density / [kg/m3]
volScalarField rhoPrimary_;
//- Viscosity / [Pa.s]
volScalarField muPrimary_;
// Sub-models
......@@ -483,6 +489,12 @@ public:
//- Pressure / [Pa]
inline const volScalarField& pPrimary() const;
//- Density / [kg/m3]
inline const volScalarField& rhoPrimary() const;
//- Viscosity / [Pa.s]
inline const volScalarField& muPrimary() const;
// Sub-models
......
......@@ -171,6 +171,20 @@ Foam::surfaceFilmModels::kinematicSingleLayer::pPrimary() const
}
inline const Foam::volScalarField&
Foam::surfaceFilmModels::kinematicSingleLayer::rhoPrimary() const
{
return rhoPrimary_;
}
inline const Foam::volScalarField&
Foam::surfaceFilmModels::kinematicSingleLayer::muPrimary() const
{
return muPrimary_;
}
inline Foam::surfaceFilmModels::injectionModel&
Foam::surfaceFilmModels::kinematicSingleLayer::injection()
{
......
......@@ -126,19 +126,7 @@ void Foam::surfaceFilmModels::thermoSingleLayer::updateSurfaceTemperatures()
Tw_.correctBoundaryConditions();
// Update film surface temperature
dimensionedScalar deltaSmall("SMALL", dimLength, SMALL);
volScalarField kappaDeltaBy2 = kappa_/(0.5*delta_ + deltaSmall);
Ts_ =
(
// qRad
- energyPhaseChangeForPrimary_/(time_.deltaT()*magSf_)
+ TPrimary_*htcs_->h()
+ kappaDeltaBy2*T_
)
/(
htcs_->h()
+ kappaDeltaBy2
);
Ts_ = T_;
Ts_.correctBoundaryConditions();
}
......@@ -581,10 +569,12 @@ void Foam::surfaceFilmModels::thermoSingleLayer::info() const
{
kinematicSingleLayer::info();
Info<< indent << "min/max(T) = " << min(T_).value() << ", "
Info<< indent << "min/max(T) = " << min(T_).value() << ", "
<< max(T_).value() << nl
<< indent << "mass phase change = "
<< returnReduce(totalMassPhaseChange_, sumOp<scalar>()) << nl;
<< indent << "mass phase change = "
<< returnReduce(totalMassPhaseChange_, sumOp<scalar>()) << nl
<< indent << "vapourisation rate = "
<< sum(massPhaseChangeForPrimary_).value()/time_.deltaTValue() << nl;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment