diff --git a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C index 698a76b81f..02ba442f2e 100644 --- a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C +++ b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C @@ -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(), diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H index ccc13b8035..9efad2770e 100644 --- a/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H +++ b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H @@ -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; } diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C index aab07b4d7d..9b4beff0ee 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C @@ -511,7 +511,7 @@ void Foam::ReactingMultiphaseParcel::calcDevolatilisation if ( !td.cloud().devolatilisation().active() - || T < td.cloud().constProps().Tvap() + || T < td.cloud().constProps().TDevol() || canCombust == -1 ) { diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H index ff13d240b9..5f1093b36d 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H @@ -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; diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H index 4a4a97eddf..0b062a6426 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H @@ -30,6 +30,7 @@ inline Foam::ReactingMultiphaseParcel::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 +inline Foam::ReactingMultiphaseParcel::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 inline Foam::ReactingMultiphaseParcel::ReactingMultiphaseParcel ( @@ -148,6 +196,14 @@ inline Foam::ReactingMultiphaseParcel::ReactingMultiphaseParcel // * * * * * * * * * constantProperties Member Functions * * * * * * * * * * // +template +inline Foam::scalar +Foam::ReactingMultiphaseParcel::constantProperties::TDevol() const +{ + return TDevol_; +} + + template inline Foam::scalar Foam::ReactingMultiphaseParcel::constantProperties::LDevol() const diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C index 9c424896a0..9d9a4cdbe8 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C @@ -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::calcPhaseChange scalarField& Cs ) { - if - ( - !td.cloud().phaseChange().active() - || T < td.cloud().constProps().Tvap() - || YPhase < SMALL - ) + typedef typename TrackData::cloudType::reactingCloudType reactingCloudType; + PhaseChangeModel& 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& 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::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& 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; } diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H index 854e2f5730..2edffd92be 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H @@ -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; }; diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H index 1329741c45..fd5700f05b 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H @@ -31,8 +31,7 @@ Foam::ReactingParcel::constantProperties::constantProperties() : ParcelType::constantProperties(), pMin_(0.0), - constantVolume_(false), - Tvap_(0.0) + constantVolume_(false) {} @@ -44,8 +43,7 @@ inline Foam::ReactingParcel::constantProperties::constantProperties : ParcelType::constantProperties(cp), pMin_(cp.pMin_), - constantVolume_(cp.constantVolume_), - Tvap_(cp.Tvap_) + constantVolume_(cp.constantVolume_) {} @@ -58,8 +56,7 @@ inline Foam::ReactingParcel::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::constantProperties::constantProperties } this->dict().lookup("constantVolume") >> constantVolume_; - this->dict().lookup("Tvap") >> Tvap_; } } @@ -91,8 +87,7 @@ inline Foam::ReactingParcel::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::constantProperties::constantProperties Pr ), pMin_(pMin), - constantVolume_(constantVolume), - Tvap_(Tvap) + constantVolume_(constantVolume) {} @@ -198,14 +192,6 @@ Foam::ReactingParcel::constantProperties::constantVolume() const } -template -inline Foam::scalar -Foam::ReactingParcel::constantProperties::Tvap() const -{ - return Tvap_; -} - - // * * * * * * * * * * ThermoParcel Member Functions * * * * * * * * * * * * // template diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C index 8bd8a4d140..61dfb6c6ee 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C +++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C @@ -151,6 +151,44 @@ void Foam::LiquidEvaporation::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::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::dh template -Foam::scalar Foam::LiquidEvaporation::TMax(const scalar pIn) const +Foam::scalar Foam::LiquidEvaporation::Tvap +( + const scalarField& Y +) const { - return liquids_.TMax(pIn); + const scalarField X(liquids_.X(Y)); + + return liquids_.Tpt(X); +} + + +template +Foam::scalar Foam::LiquidEvaporation::TMax +( + const scalar p, + const scalarField& Y +) const +{ + const scalarField X(liquids_.X(Y)); + + return liquids_.pvInvert(p, X); } diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H index fb2e2cdd9c..3d5334cd88 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H +++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H @@ -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; }; diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C index 6436d39605..d78d6b110b 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C +++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C @@ -151,18 +151,53 @@ void Foam::LiquidEvaporationBoil::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::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::dh template -Foam::scalar Foam::LiquidEvaporationBoil::TMax +Foam::scalar Foam::LiquidEvaporationBoil::Tvap ( - const scalar pIn + const scalarField& Y ) const { - return liquids_.TMax(pIn); + const scalarField X(liquids_.X(Y)); + + return liquids_.Tpt(X); +} + + +template +Foam::scalar Foam::LiquidEvaporationBoil::TMax +( + const scalar p, + const scalarField& Y +) const +{ + const scalarField X(liquids_.X(Y)); + + return liquids_.pvInvert(p, X); } diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H index 2d11f55f32..f3b532dff5 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H +++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H @@ -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; }; diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C index 3ba13f27e1..9ec9307278 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C +++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C @@ -179,12 +179,23 @@ Foam::scalar Foam::PhaseChangeModel::dh template -Foam::scalar Foam::PhaseChangeModel::TMax(const scalar) const +Foam::scalar Foam::PhaseChangeModel::TMax +( + const scalar, + const scalarField& +) const { return GREAT; } +template +Foam::scalar Foam::PhaseChangeModel::Tvap(const scalarField& Y) const +{ + return -GREAT; +} + + template void Foam::PhaseChangeModel::addToPhaseChangeMass(const scalar dMass) { diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H index f61031c758..6f3a171708 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H +++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H @@ -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); diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C index 876fa675a1..198dc717f8 100644 --- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C +++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C @@ -64,47 +64,62 @@ void Foam::SprayParcel::calc const label cellI ) { + typedef typename TrackData::cloudType::reactingCloudType reactingCloudType; + const CompositionModel& 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; diff --git a/src/thermophysicalModels/properties/liquidMixtureProperties/liquidMixtureProperties/liquidMixtureProperties.C b/src/thermophysicalModels/properties/liquidMixtureProperties/liquidMixtureProperties/liquidMixtureProperties.C index 0360731e21..ac26ec1390 100644 --- a/src/thermophysicalModels/properties/liquidMixtureProperties/liquidMixtureProperties/liquidMixtureProperties.C +++ b/src/thermophysicalModels/properties/liquidMixtureProperties/liquidMixtureProperties/liquidMixtureProperties.C @@ -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; diff --git a/src/thermophysicalModels/properties/liquidMixtureProperties/liquidMixtureProperties/liquidMixtureProperties.H b/src/thermophysicalModels/properties/liquidMixtureProperties/liquidMixtureProperties/liquidMixtureProperties.H index 0ca2cd0c3e..2e56b02a2f 100644 --- a/src/thermophysicalModels/properties/liquidMixtureProperties/liquidMixtureProperties/liquidMixtureProperties.H +++ b/src/thermophysicalModels/properties/liquidMixtureProperties/liquidMixtureProperties/liquidMixtureProperties.H @@ -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; diff --git a/src/thermophysicalModels/properties/liquidProperties/liquidProperties/liquidProperties.H b/src/thermophysicalModels/properties/liquidProperties/liquidProperties/liquidProperties.H index 7cb4aa0710..2210d567a3 100644 --- a/src/thermophysicalModels/properties/liquidProperties/liquidProperties/liquidProperties.H +++ b/src/thermophysicalModels/properties/liquidProperties/liquidProperties/liquidProperties.H @@ -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;