Commit 2d2787bd authored by Henry's avatar Henry
Browse files

thermodynamics: Added pressure as an addition argument to all primitive thermodynamic functions

Added additional layer of templating to reactingMixture to support specie functions in a generic manner.
parent 0cdb34db
......@@ -87,7 +87,7 @@
forAll(Y, i)
{
Y[i] = Y0[i];
hs0 += Y0[i]*specieData[i].Hs(T0);
hs0 += Y0[i]*specieData[i].Hs(p[i], T0);
}
hs = dimensionedScalar("h", dimEnergy/dimMass, hs0);
......
......@@ -211,7 +211,7 @@ int main(int argc, char *argv[])
co = co2*
min
(
CO2BreakUp.Kn(equilibriumFlameTemperature, P, N)
CO2BreakUp.Kn(P, equilibriumFlameTemperature, N)
/::sqrt(max(ores, 0.001)),
1.0
);
......@@ -219,7 +219,7 @@ int main(int argc, char *argv[])
h2 = h2o*
min
(
H2OBreakUp.Kn(equilibriumFlameTemperature, P, N)
H2OBreakUp.Kn(P, equilibriumFlameTemperature, N)
/::sqrt(max(ores, 0.001)),
1.0
);
......
......@@ -424,7 +424,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
{
scalar dm = np0*dMassGas[i];
label gid = composition.localToGlobalCarrierId(GAS, i);
scalar hs = composition.carrier().Hs(gid, T0);
scalar hs = composition.carrier().Hs(gid, pc, T0);
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().hsTrans()[cellI] += dm*hs;
......@@ -433,7 +433,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
{
scalar dm = np0*dMassLiquid[i];
label gid = composition.localToGlobalCarrierId(LIQ, i);
scalar hs = composition.carrier().Hs(gid, T0);
scalar hs = composition.carrier().Hs(gid, pc, T0);
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().hsTrans()[cellI] += dm*hs;
......@@ -444,7 +444,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
{
scalar dm = np0*dMassSolid[i];
label gid = composition.localToGlobalCarrierId(SLD, i);
scalar hs = composition.carrier().Hs(gid, T0);
scalar hs = composition.carrier().Hs(gid, pc, T0);
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().hsTrans()[cellI] += dm*hs;
......@@ -453,7 +453,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
forAll(dMassSRCarrier, i)
{
scalar dm = np0*dMassSRCarrier[i];
scalar hs = composition.carrier().Hs(i, T0);
scalar hs = composition.carrier().Hs(i, pc, T0);
td.cloud().rhoTrans(i)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
td.cloud().hsTrans()[cellI] += dm*hs;
......@@ -541,7 +541,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
forAll(dMassDV, i)
{
const label id = composition.localToGlobalCarrierId(GAS, i);
const scalar Cp = composition.carrier().Cp(id, Ts);
const scalar Cp = composition.carrier().Cp(id, this->pc_, Ts);
const scalar W = composition.carrier().W(id);
const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
......
......@@ -104,7 +104,12 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
forAll(td.cloud().rhoTrans(), i)
{
scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
CpEff += Y*td.cloud().composition().carrier().Cp(i, this->Tc_);
CpEff += Y*td.cloud().composition().carrier().Cp
(
i,
this->pc_,
this->Tc_
);
}
const scalar Cpc = td.CpInterp().psi()[cellI];
......@@ -206,9 +211,9 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
const scalar cbrtW = cbrt(W);
rhos += Xs[i]*W;
mus += Ys[i]*sqrtW*thermo.carrier().mu(i, T);
kappas += Ys[i]*cbrtW*thermo.carrier().kappa(i, T);
Cps += Xs[i]*thermo.carrier().Cp(i, T);
mus += Ys[i]*sqrtW*thermo.carrier().mu(i, pc_, T);
kappas += Ys[i]*cbrtW*thermo.carrier().kappa(i, pc_, T);
Cps += Xs[i]*thermo.carrier().Cp(i, pc_, T);
sumYiSqrtW += Ys[i]*sqrtW;
sumYiCbrtW += Ys[i]*cbrtW;
......@@ -378,7 +383,7 @@ void Foam::ReactingParcel<ParcelType>::calc
{
scalar dmi = dm*Y_[i];
label gid = composition.localToGlobalCarrierId(0, i);
scalar hs = composition.carrier().Hs(gid, T0);
scalar hs = composition.carrier().Hs(gid, pc_, T0);
td.cloud().rhoTrans(gid)[cellI] += dmi;
td.cloud().hsTrans()[cellI] += dmi*hs;
......@@ -439,7 +444,7 @@ void Foam::ReactingParcel<ParcelType>::calc
{
scalar dm = np0*dMass[i];
label gid = composition.localToGlobalCarrierId(0, i);
scalar hs = composition.carrier().Hs(gid, T0);
scalar hs = composition.carrier().Hs(gid, pc_, T0);
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().UTrans()[cellI] += dm*U0;
......@@ -543,7 +548,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
const label idc = composition.localToGlobalCarrierId(idPhase, i);
const label idl = composition.globalIds(idPhase)[i];
const scalar Cp = composition.carrier().Cp(idc, Ts);
const scalar Cp = composition.carrier().Cp(idc, pc_, Ts);
const scalar W = composition.carrier().W(idc);
const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
......
......@@ -396,7 +396,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::H
forAll(Y, i)
{
label gid = props.globalIds()[i];
HMixture += Y[i]*thermo_.carrier().Hs(gid, T);
HMixture += Y[i]*thermo_.carrier().Hs(gid, p, T);
}
break;
}
......@@ -460,7 +460,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::Hs
forAll(Y, i)
{
label gid = props.globalIds()[i];
HsMixture += Y[i]*thermo_.carrier().Hs(gid, T);
HsMixture += Y[i]*thermo_.carrier().Hs(gid, p, T);
}
break;
}
......@@ -584,7 +584,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::Cp
forAll(Y, i)
{
label gid = props.globalIds()[i];
CpMixture += Y[i]*thermo_.carrier().Cp(gid, T);
CpMixture += Y[i]*thermo_.carrier().Cp(gid, p, T);
}
break;
}
......
......@@ -216,7 +216,7 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
}
case (parent::etEnthalpyDifference):
{
scalar hc = this->owner().composition().carrier().Hs(idc, T);
scalar hc = this->owner().composition().carrier().Hs(idc, p, T);
scalar hp = liquids_.properties()[idl].h(p, T);
dh = hc - hp;
......
......@@ -171,10 +171,10 @@ void Foam::LiquidEvaporationBoil<CloudType>::calculate
forAll(this->owner().thermo().carrier().Y(), i)
{
scalar Yc = this->owner().thermo().carrier().Y()[i][cellI];
Hc += Yc*this->owner().thermo().carrier().Hs(i, Tc);
Hsc += Yc*this->owner().thermo().carrier().Hs(i, Ts);
Cpc += Yc*this->owner().thermo().carrier().Cp(i, Ts);
kappac += Yc*this->owner().thermo().carrier().kappa(i, Ts);
Hc += Yc*this->owner().thermo().carrier().Hs(i, pc, Tc);
Hsc += Yc*this->owner().thermo().carrier().Hs(i, ps, Ts);
Cpc += Yc*this->owner().thermo().carrier().Cp(i, ps, Ts);
kappac += Yc*this->owner().thermo().carrier().kappa(i, ps, Ts);
}
// calculate mass transfer of each specie in liquid
......@@ -315,7 +315,7 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
}
case (parent::etEnthalpyDifference):
{
scalar hc = this->owner().composition().carrier().Hs(idc, TDash);
scalar hc = this->owner().composition().carrier().Hs(idc, p, TDash);
scalar hp = liquids_.properties()[idl].h(p, TDash);
dh = hc - hp;
......
......@@ -146,7 +146,8 @@ void reactingOneDim::updatePhiGas()
forAll(gasTable, gasI)
{
tmp<volScalarField> tHsiGas = solidChemistry_->gasHs(T_, gasI);
tmp<volScalarField> tHsiGas =
solidChemistry_->gasHs(p, T_, gasI);
tmp<volScalarField> tRRiGas = solidChemistry_->RRg(gasI);
const volScalarField& HsiGas = tHsiGas();
......
......@@ -153,6 +153,7 @@ public:
//- Enthalpy/Internal energy for cell-set [J/kg]
virtual tmp<scalarField> he
(
const scalarField& p,
const scalarField& T,
const labelList& cells
) const = 0;
......@@ -160,6 +161,7 @@ public:
//- Enthalpy/Internal energy for patch [J/kg]
virtual tmp<scalarField> he
(
const scalarField& p,
const scalarField& T,
const label patchi
) const = 0;
......@@ -171,6 +173,7 @@ public:
virtual tmp<scalarField> THE
(
const scalarField& h,
const scalarField& p,
const scalarField& T0, // starting temperature
const labelList& cells
) const = 0;
......@@ -179,6 +182,7 @@ public:
virtual tmp<scalarField> THE
(
const scalarField& h,
const scalarField& p,
const scalarField& T0, // starting temperature
const label patchi
) const = 0;
......@@ -195,6 +199,7 @@ public:
//- Heat capacity at constant pressure for patch [J/kg/K]
virtual tmp<scalarField> Cp
(
const scalarField& p,
const scalarField& T,
const label patchi
) const = 0;
......@@ -205,6 +210,7 @@ public:
//- Heat capacity at constant volume for patch [J/kg/K]
virtual tmp<scalarField> Cv
(
const scalarField& p,
const scalarField& T,
const label patchi
) const = 0;
......@@ -215,6 +221,7 @@ public:
//- gamma = Cp/Cv for patch []
virtual tmp<scalarField> gamma
(
const scalarField& p,
const scalarField& T,
const label patchi
) const = 0;
......@@ -225,6 +232,7 @@ public:
//- Heat capacity at constant pressure/volume for patch [J/kg/K]
virtual tmp<scalarField> Cpv
(
const scalarField& p,
const scalarField& T,
const label patchi
) const = 0;
......@@ -235,6 +243,7 @@ public:
//- Heat capacity ratio for patch []
virtual tmp<scalarField> CpByCpv
(
const scalarField& p,
const scalarField& T,
const label patchi
) const = 0;
......
......@@ -110,6 +110,7 @@ void Foam::energyJumpFvPatchScalarField::updateCoeffs()
label patchID = patch().index();
const scalarField& pp = thermo.p().boundaryField()[patchID];
const temperatureJumpFvPatchScalarField& TbPatch =
refCast<const temperatureJumpFvPatchScalarField>
(
......@@ -124,7 +125,7 @@ void Foam::energyJumpFvPatchScalarField::updateCoeffs()
const labelUList& faceCells = this->patch().faceCells();
jump_ = thermo.he(jumpTb, faceCells);
jump_ = thermo.he(pp, jumpTb, faceCells);
}
fixedJumpFvPatchField<scalar>::updateCoeffs();
......
......@@ -104,10 +104,11 @@ void Foam::fixedEnergyFvPatchScalarField::updateCoeffs()
const label patchi = patch().index();
const scalarField& pw = thermo.p().boundaryField()[patchi];
fvPatchScalarField& Tw =
const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
Tw.evaluate();
operator==(thermo.he(Tw, patchi));
operator==(thermo.he(pw, Tw, patchi));
fixedValueFvPatchScalarField::updateCoeffs();
}
......
......@@ -104,16 +104,17 @@ void Foam::gradientEnergyFvPatchScalarField::updateCoeffs()
const label patchi = patch().index();
const scalarField& pw = thermo.p().boundaryField()[patchi];
fvPatchScalarField& Tw =
const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
Tw.evaluate();
gradient() = thermo.Cpv(Tw, patchi)*Tw.snGrad()
gradient() = thermo.Cpv(pw, Tw, patchi)*Tw.snGrad()
+ patch().deltaCoeffs()*
(
thermo.he(Tw, patchi)
- thermo.he(Tw, patch().faceCells())
thermo.he(pw, Tw, patchi)
- thermo.he(pw, Tw, patch().faceCells())
);
fixedGradientFvPatchScalarField::updateCoeffs();
......
......@@ -109,6 +109,7 @@ void Foam::mixedEnergyFvPatchScalarField::updateCoeffs()
const label patchi = patch().index();
const scalarField& pw = thermo.p().boundaryField()[patchi];
mixedFvPatchScalarField& Tw = refCast<mixedFvPatchScalarField>
(
const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi])
......@@ -117,13 +118,13 @@ void Foam::mixedEnergyFvPatchScalarField::updateCoeffs()
Tw.evaluate();
valueFraction() = Tw.valueFraction();
refValue() = thermo.he(Tw.refValue(), patchi);
refValue() = thermo.he(pw, Tw.refValue(), patchi);
refGrad() =
thermo.Cpv(Tw, patchi)*Tw.refGrad()
thermo.Cpv(pw, Tw, patchi)*Tw.refGrad()
+ patch().deltaCoeffs()*
(
thermo.he(Tw, patchi)
- thermo.he(Tw, patch().faceCells())
thermo.he(pw, Tw, patchi)
- thermo.he(pw, Tw, patch().faceCells())
);
mixedFvPatchScalarField::updateCoeffs();
......
......@@ -156,8 +156,9 @@ void Foam::wallHeatTransferFvPatchScalarField::updateCoeffs()
const label patchi = patch().index();
const scalarField& pw = thermo.p().boundaryField()[patchi];
const scalarField& Tw = thermo.T().boundaryField()[patchi];
const scalarField Cpw(thermo.Cp(Tw, patchi));
const scalarField Cpw(thermo.Cp(pw, Tw, patchi));
valueFraction() =
1.0/
......
......@@ -49,17 +49,23 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo(const fvMesh& mesh)
)
{
scalarField& heCells = he_.internalField();
const scalarField& pCells = this->p_.internalField();
const scalarField& TCells = this->T_.internalField();
forAll(heCells, celli)
{
heCells[celli] = this->cellMixture(celli).HE(TCells[celli]);
heCells[celli] =
this->cellMixture(celli).HE(pCells[celli], TCells[celli]);
}
forAll(he_.boundaryField(), patchi)
{
he_.boundaryField()[patchi] ==
he(this->T_.boundaryField()[patchi], patchi);
he_.boundaryField()[patchi] == he
(
this->p_.boundaryField()[patchi],
this->T_.boundaryField()[patchi],
patchi
);
}
this->heBoundaryCorrection(he_);
......@@ -78,6 +84,7 @@ Foam::heThermo<BasicThermo, MixtureType>::~heThermo()
template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
(
const scalarField& p,
const scalarField& T,
const labelList& cells
) const
......@@ -87,7 +94,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
forAll(T, celli)
{
he[celli] = this->cellMixture(cells[celli]).HE(T[celli]);
he[celli] = this->cellMixture(cells[celli]).HE(p[celli], T[celli]);
}
return the;
......@@ -97,6 +104,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
(
const scalarField& p,
const scalarField& T,
const label patchi
) const
......@@ -106,7 +114,8 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
forAll(T, facei)
{
he[facei] = this->patchFaceMixture(patchi, facei).HE(T[facei]);
he[facei] =
this->patchFaceMixture(patchi, facei).HE(p[facei], T[facei]);
}
return the;
......@@ -161,6 +170,7 @@ Foam::heThermo<BasicThermo, MixtureType>::hc() const
template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cp
(
const scalarField& p,
const scalarField& T,
const label patchi
) const
......@@ -170,7 +180,8 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cp
forAll(T, facei)
{
cp[facei] = this->patchFaceMixture(patchi, facei).Cp(T[facei]);
cp[facei] =
this->patchFaceMixture(patchi, facei).Cp(p[facei], T[facei]);
}
return tCp;
......@@ -204,17 +215,20 @@ Foam::heThermo<BasicThermo, MixtureType>::Cp() const
forAll(this->T_, celli)
{
cp[celli] = this->cellMixture(celli).Cp(this->T_[celli]);
cp[celli] =
this->cellMixture(celli).Cp(this->p_[celli], this->T_[celli]);
}
forAll(this->T_.boundaryField(), patchi)
{
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
const fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
fvPatchScalarField& pCp = cp.boundaryField()[patchi];
forAll(pT, facei)
{
pCp[facei] = this->patchFaceMixture(patchi, facei).Cp(pT[facei]);
pCp[facei] =
this->patchFaceMixture(patchi, facei).Cp(pp[facei], pT[facei]);
}
}
......@@ -226,6 +240,7 @@ template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::scalarField>
Foam::heThermo<BasicThermo, MixtureType>::Cv
(
const scalarField& p,
const scalarField& T,
const label patchi
) const
......@@ -235,7 +250,8 @@ Foam::heThermo<BasicThermo, MixtureType>::Cv
forAll(T, facei)
{
cv[facei] = this->patchFaceMixture(patchi, facei).Cv(T[facei]);
cv[facei] =
this->patchFaceMixture(patchi, facei).Cv(p[facei], T[facei]);
}
return tCv;
......@@ -269,13 +285,18 @@ Foam::heThermo<BasicThermo, MixtureType>::Cv() const
forAll(this->T_, celli)
{
cv[celli] = this->cellMixture(celli).Cv(this->T_[celli]);
cv[celli] =
this->cellMixture(celli).Cv(this->p_[celli], this->T_[celli]);
}
forAll(this->T_.boundaryField(), patchi)
{
cv.boundaryField()[patchi] =
Cv(this->T_.boundaryField()[patchi], patchi);
cv.boundaryField()[patchi] = Cv
(
this->p_.boundaryField()[patchi],
this->T_.boundaryField()[patchi],
patchi
);
}
return tCv;
......@@ -285,6 +306,7 @@ Foam::heThermo<BasicThermo, MixtureType>::Cv() const
template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::gamma
(
const scalarField& p,
const scalarField& T,
const label patchi
) const
......@@ -294,7 +316,8 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::gamma
forAll(T, facei)
{
cpv[facei] = this->patchFaceMixture(patchi, facei).gamma(T[facei]);
cpv[facei] =
this->patchFaceMixture(patchi, facei).gamma(p[facei], T[facei]);
}
return tgamma;
......@@ -328,18 +351,23 @@ Foam::heThermo<BasicThermo, MixtureType>::gamma() const
forAll(this->T_, celli)
{
cpv[celli] = this->cellMixture(celli).gamma(this->T_[celli]);
cpv[celli] =
this->cellMixture(celli).gamma(this->p_[celli], this->T_[celli]);
}
forAll(this->T_.boundaryField(), patchi)
{
const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
const fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
fvPatchScalarField& pgamma = cpv.boundaryField()[patchi];
forAll(pT, facei)
{
pgamma[facei] =
this->patchFaceMixture(patchi, facei).gamma(pT[facei]);
pgamma[facei] = this->patchFaceMixture(patchi, facei).gamma
(
pp[facei],
pT[facei]
);
}
}
......@@ -350,6 +378,7 @@ Foam::heThermo<BasicThermo, MixtureType>::gamma() const
template<class BasicThermo, class MixtureType>
Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cpv
(
const scalarField& p,
const scalarField& T,
const label patchi
) const
......@@ -359,7 +388,8 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cpv