diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C index 5f4e839508..a7f92224bd 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C @@ -378,8 +378,11 @@ updateCoeffs() // TvNbr: vapour Tp scalarField c(TcNbr*KDeltaLiqNbr + TvNbr*KDeltaVapNbr); - valueFraction() = KDeltaNbr/(KDeltaNbr + KDelta); - refValue() = c/KDeltaNbr; + //valueFraction() = KDeltaNbr/(KDeltaNbr + KDelta); + //refValue() = c/KDeltaNbr; + scalarField KDeltaLiqVapNbr(KDeltaLiqNbr + KDeltaVapNbr); + valueFraction() = KDeltaLiqVapNbr/(KDeltaLiqVapNbr + KDelta); + refValue() = c/KDeltaLiqVapNbr; refGrad() = (qr + qrNbr)/kappa(Tp); if (debug) @@ -393,7 +396,7 @@ updateCoeffs() << " walltemperature " << " min:" << gMin(Tp) << " max:" << gMax(Tp) - << " avg:" << gAverage(Tp) + << " avg:" << gAverage(Tp) << nl << endl; } } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/createFluidFields.H index 6d1e3b9bcc..00e0f30b23 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/createFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/createFluidFields.H @@ -21,7 +21,178 @@ const uniformDimensionedVectorField& g = meshObjects::gravity::New(runTime); PtrList pimpleFluid(fluidRegions.size()); -// Populate fluid field pointer lists + +//Debug Fields +/* +PtrList faceRegimesFluid(fluidRegions.size()); +PtrList qcFluid(fluidRegions.size()); +PtrList qFilmFluid(fluidRegions.size()); +PtrList htcFilmBoilingFluid(fluidRegions.size()); +PtrList qtbFluid(fluidRegions.size()); +PtrList qSubCoolFluid(fluidRegions.size()); +PtrList CHFtotalFluid(fluidRegions.size()); +PtrList TdnbFluid(fluidRegions.size()); +PtrList phiFluid(fluidRegions.size()); + +forAll(fluidRegions, i) +{ + faceRegimesFluid.set + ( + i, + new volScalarField + ( + IOobject + ( + "faceRegimes", + runTime.timeName(), + fluidRegions[i], + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluidRegions[i], + dimensionedScalar(dimless, Zero) + + ) + ); + qcFluid.set + ( + i, + new volScalarField + ( + IOobject + ( + "qc", + runTime.timeName(), + fluidRegions[i], + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluidRegions[i], + dimensionedScalar(dimless, Zero) + ) + ); + qFilmFluid.set + ( + i, + new volScalarField + ( + IOobject + ( + "qFilm", + runTime.timeName(), + fluidRegions[i], + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluidRegions[i], + dimensionedScalar(dimless, Zero) + ) + ); + htcFilmBoilingFluid.set + ( + i, + new volScalarField + ( + IOobject + ( + "htcFilmBoiling", + runTime.timeName(), + fluidRegions[i], + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluidRegions[i], + dimensionedScalar(dimless, Zero) + ) + ); + qtbFluid.set + ( + i, + new volScalarField + ( + IOobject + ( + "qtb", + runTime.timeName(), + fluidRegions[i], + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluidRegions[i], + dimensionedScalar(dimless, Zero) + ) + ); + qSubCoolFluid.set + ( + i, + new volScalarField + ( + IOobject + ( + "qSubCool", + runTime.timeName(), + fluidRegions[i], + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluidRegions[i], + dimensionedScalar(dimless, Zero) + ) + ); + CHFtotalFluid.set + ( + i, + new volScalarField + ( + IOobject + ( + "CHFtotal", + runTime.timeName(), + fluidRegions[i], + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluidRegions[i], + dimensionedScalar(dimless, Zero) + ) + ); + TdnbFluid.set + ( + i, + new volScalarField + ( + IOobject + ( + "Tdnb", + runTime.timeName(), + fluidRegions[i], + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluidRegions[i], + dimensionedScalar(dimless, Zero) + ) + ); + phiFluid.set + ( + i, + new volScalarField + ( + IOobject + ( + "phiTb", + runTime.timeName(), + fluidRegions[i], + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluidRegions[i], + dimensionedScalar(dimless, Zero) + ) + ); +} +*/ + forAll(fluidRegions, i) { Info<< "*** Reading fluid mesh thermophysical properties for region " diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H index 7b55ecf529..7543293492 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H @@ -56,3 +56,16 @@ pimpleControl& pimple = pimpleFluid[i]; const dimensionedScalar& pMin = pMinFluid[i]; + + // Debug fields +/* + volScalarField& faceRegimes = faceRegimesFluid[i]; + volScalarField& qc = qcFluid[i]; + volScalarField& qFilm = qFilmFluid[i]; + volScalarField& htcFilmBoiling = htcFilmBoilingFluid[i]; + volScalarField& qtb = qtbFluid[i]; + volScalarField& qSubCool = qSubCoolFluid[i]; + volScalarField& CHFtotal = CHFtotalFluid[i]; + volScalarField& Tdnb = TdnbFluid[i]; + volScalarField& phiTb = phiFluid[i]; +*/ diff --git a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C index bb8858fb01..be52a40152 100644 --- a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C +++ b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C @@ -74,7 +74,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF), otherPhaseName_("vapor"), phaseType_(liquidPhase), - relax_(0.5), + relax_(), AbyV_(p.size(), 0), alphatConv_(p.size(), 0), dDep_(p.size(), 1e-5), @@ -112,7 +112,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF, dict), otherPhaseName_(dict.lookup("otherPhase")), phaseType_(phaseTypeNames_.read(dict.lookup("phaseType"))), - relax_(dict.lookupOrDefault("relax", 0.5)), + relax_(Function1::New("relax", dict)), AbyV_(p.size(), 0), alphatConv_(p.size(), 0), dDep_(p.size(), 1e-5), @@ -155,20 +155,6 @@ alphatWallBoilingWallFunctionFvPatchScalarField const dictionary* LeidenfrostDict = dict.findDict("LeidenfrostModel"); - if (LeidenfrostDict) - { - LeidenfrostModel_ = - wallBoilingModels::LeidenfrostModel::New(*LeidenfrostDict); - } - - const dictionary* filmDict = dict.findDict("filmBoilingModel"); - - if (filmDict) - { - filmBoilingModel_ = - wallBoilingModels::filmBoilingModel::New(*filmDict); - } - dmdt_ = 0; break; @@ -304,7 +290,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField ), otherPhaseName_(psf.otherPhaseName_), phaseType_(psf.phaseType_), - relax_(psf.relax_), + relax_(psf.relax_.clone()), AbyV_(psf.AbyV_), alphatConv_(psf.alphatConv_, mapper), dDep_(psf.dDep_, mapper), @@ -332,7 +318,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf), otherPhaseName_(psf.otherPhaseName_), phaseType_(psf.phaseType_), - relax_(psf.relax_), + relax_(psf.relax_.clone()), AbyV_(psf.AbyV_), alphatConv_(psf.alphatConv_), dDep_(psf.dDep_), @@ -361,7 +347,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf, iF), otherPhaseName_(psf.otherPhaseName_), phaseType_(psf.phaseType_), - relax_(psf.relax_), + relax_(psf.relax_.clone()), AbyV_(psf.AbyV_), alphatConv_(psf.alphatConv_), dDep_(psf.dDep_), @@ -457,6 +443,9 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs() const label patchi = patch().index(); + const scalar t = this->db().time().timeOutputValue(); + scalar relax = relax_->value(t); + switch (phaseType_) { case vaporPhase: @@ -469,102 +458,9 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs() const fvPatchScalarField& hewv = vapor.thermo().he().boundaryField()[patchi]; - const phaseModel& liquid(fluid.phases()[otherPhaseName_]); - - const phaseCompressibleTurbulenceModel& turbModel = - db().lookupObject - ( - IOobject::groupName - ( - turbulenceModel::propertiesName, - liquid.name() - ) - ); - - const phaseCompressibleTurbulenceModel& vaporTurbModel = - db().lookupObject - ( - IOobject::groupName - ( - turbulenceModel::propertiesName, - vapor.name() - ) - ); - - const fvPatchScalarField& rhoVaporw = - vaporTurbModel.rho().boundaryField()[patchi]; - // Vapor Liquid phase fraction at the wall const scalarField vaporw(vapor.boundaryField()[patchi]); - const fvPatchScalarField& Tw = - liquid.thermo().T().boundaryField()[patchi]; - const scalarField Tc(Tw.patchInternalField()); - - // Saturation temperature - const tmp tTsat = - satModel.Tsat(liquid.thermo().p()); - const volScalarField& Tsat = tTsat(); - const fvPatchScalarField& Tsatw(Tsat.boundaryField()[patchi]); - const scalarField Tsatc(Tsatw.patchInternalField()); - - const fvPatchScalarField& hewl = - liquid.thermo().he().boundaryField()[patchi]; - - const fvPatchScalarField& pw = - liquid.thermo().p().boundaryField()[patchi]; - - const fvPatchScalarField& rhow = - turbModel.rho().boundaryField()[patchi]; - - const scalarField hw - ( - liquid.thermo().he().member() == "e" - ? hewl.patchInternalField() + pw/rhow.patchInternalField() - : hewl.patchInternalField() - ); - - const scalarField L - ( - vapor.thermo().he().member() == "e" - ? vapor.thermo().he(pw, Tsatc, patchi) + pw/rhoVaporw - hw - : vapor.thermo().he(pw, Tsatc, patchi) - hw - ); - - // Film boiling models - - scalarField htcFilmBoiling(this->size(), 0); - scalarField TLeiden(this->size(), GREAT); - - if (filmBoilingModel_.valid() && LeidenfrostModel_.valid()) - { - // htc for film boiling - htcFilmBoiling = - filmBoilingModel_->htcFilmBoil - ( - liquid, - vapor, - patchi, - Tc, - Tsatw, - L - ); - - // Leidenfrost Temperature - TLeiden = - LeidenfrostModel_->TLeid - ( - liquid, - vapor, - patchi, - Tc, - Tsatw, - L - ); - } - - const scalarField qFilm(htcFilmBoiling*max(Tw - Tsatw, scalar(0))); - // NOTE! Assumes 1-thisPhase for liquid fraction in // multiphase simulations const scalarField fLiquid @@ -582,30 +478,11 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs() forAll (*this, i) { - if (Tw[i] > TLeiden[i]) - { - this->operator[](i) = - ( - max - ( - (1 - fLiquid[i]) - *( - (qFilm[i]/heSnGrad[i]) - /max(vaporw[i], scalar(1e-8)) - - alphaw[i] - ), - -alphaw[i] - ) - ); - } - else - { - this->operator[](i) = - ( - (1 - fLiquid[i])*(alphatv[i]) - /max(vaporw[i], scalar(1e-8)) - ); - } + this->operator[](i) = + ( + (1 - fLiquid[i])*(alphatv[i]) + /max(vaporw[i], scalar(1e-8)) + ); } if (debug) @@ -759,132 +636,36 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs() const scalarField fLiquid(partitioningModel_->fLiquid(liquidw)); - for (label i=0; i<10; i++) + // Liquid temperature at y+=250 is estimated from logarithmic + // thermal wall function (Koncar, Krepper & Egorov, 2005) + const scalarField Tplus_y250 + ( + Prt_*(log(E_*250)/kappa_ + P) + ); + const scalarField Tplus(Prt_*(log(E_*yPlus)/kappa_ + P)); + scalarField Tl(Tw - (Tplus_y250/Tplus)*(Tw - Tc)); + Tl = max(Tc - 40, Tl); + + // Film, transient boiling regimes + scalarField Qtb(this->size(), 0); + scalarField tDNB(this->size(), GREAT); + scalarField TLeiden(this->size(), GREAT); + scalarField htcFilmBoiling(this->size(), 0); + + if + ( + CHFModel_.valid() + && CHFSoobModel_.valid() + && TDNBModel_.valid() + && MHFModel_.valid() + && LeidenfrostModel_.valid() + && filmBoilingModel_.valid() + ) { - // Liquid temperature at y+=250 is estimated from logarithmic - // thermal wall function (Koncar, Krepper & Egorov, 2005) - const scalarField Tplus_y250(Prt_*(log(E_*250)/kappa_ + P)); - const scalarField Tplus(Prt_*(log(E_*yPlus)/kappa_ + P)); - scalarField Tl(Tw - (Tplus_y250/Tplus)*(Tw - Tc)); - Tl = max(Tc - 40, Tl); - // Film, transient boiling regimes - scalarField Qtb(this->size(), 0); - scalarField tDNB(this->size(), GREAT); - scalarField TLeiden(this->size(), GREAT); - scalarField htcFilmBoiling(this->size(), 0); - - if + const scalarField CHF ( - CHFModel_.valid() - && CHFSoobModel_.valid() - && TDNBModel_.valid() - && MHFModel_.valid() - && LeidenfrostModel_.valid() - && filmBoilingModel_.valid() - ) - { - - const scalarField CHF - ( - CHFModel_->CHF - ( - liquid, - vapor, - patchi, - Tl, - Tsatw, - L - ) - ); - - // Effect of sub-cooling to the CHF in saturated conditions - const scalarField CHFSubCool - ( - CHFSoobModel_->CHFSubCool - ( - liquid, - vapor, - patchi, - Tl, - Tsatw, - L - ) - ); - - const scalarField CHFtotal(CHF*CHFSubCool); - - tDNB = - TDNBModel_->TDNB - ( - liquid, - vapor, - patchi, - Tl, - Tsatw, - L - ); - - const scalarField MHF - ( - MHFModel_->MHF - ( - liquid, - vapor, - patchi, - Tl, - Tsatw, - L - ) - ); - - TLeiden = - LeidenfrostModel_->TLeid - ( - liquid, - vapor, - patchi, - Tl, - Tsatw, - L - ); - - // htc for film boiling - htcFilmBoiling = - filmBoilingModel_->htcFilmBoil - ( - liquid, - vapor, - patchi, - Tl, - Tsatw, - L - ); - - // htc for film transition boiling - - // Indicator between CHF (phi = 0) and MHF (phi = 1) - const scalarField phi - ( - min - ( - max - ( - wp_*(Tw - tDNB)/(TLeiden - tDNB), - scalar(0) - ) - , scalar(1) - ) - ); - - Qtb = CHFtotal*(1 - phi) + phi*MHF; - } - - - // Sub-cool boiling Nucleation - const scalarField N - ( - nucleationSiteModel_->N + CHFModel_->CHF ( liquid, vapor, @@ -895,8 +676,123 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs() ) ); - // Bubble departure diameter: - dDep_ = departureDiamModel_->dDeparture + // Effect of sub-cooling to the CHF in saturated conditions + const scalarField CHFSubCool + ( + CHFSoobModel_->CHFSubCool + ( + liquid, + vapor, + patchi, + Tl, + Tsatw, + L + ) + ); + + const scalarField CHFtotal(CHF*CHFSubCool); + + tDNB = + TDNBModel_->TDNB + ( + liquid, + vapor, + patchi, + Tl, + Tsatw, + L + ); + + const scalarField MHF + ( + MHFModel_->MHF + ( + liquid, + vapor, + patchi, + Tl, + Tsatw, + L + ) + ); + + TLeiden = + LeidenfrostModel_->TLeid + ( + liquid, + vapor, + patchi, + Tl, + Tsatw, + L + ); + + // htc for film boiling + htcFilmBoiling = + filmBoilingModel_->htcFilmBoil + ( + liquid, + vapor, + patchi, + Tl, + Tsatw, + L + ); + + // htc for film transition boiling + + // Indicator between CHF (phi = 0) and MHF (phi = 1) + const scalarField phi + ( + min + ( + max + ( + wp_*(Tw - tDNB)/(TLeiden - tDNB), + scalar(0) + ), + scalar(1) + ) + ); + + Qtb = CHFtotal*(1 - phi) + phi*MHF; + + + if (debug & 2) + { + + scalarField& phiFluid = + const_cast + ( + patch().lookupPatchField + < + volScalarField, + scalar + >("phiTb") + ); + + phiFluid = phi; + + scalarField& CHFFluid = + const_cast + ( + patch().lookupPatchField + < + volScalarField, + scalar + >("CHFtotal") + ); + + CHFFluid = CHFtotal; + } + + } + + + // Sub-cool boiling Nucleation + const scalarField N + ( + nucleationSiteModel_->N ( liquid, vapor, @@ -904,369 +800,435 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs() Tl, Tsatw, L - ); + ) + ); - // Bubble departure frequency: - const scalarField fDep + // Bubble departure diameter: + dDep_ = departureDiamModel_->dDeparture + ( + liquid, + vapor, + patchi, + Tl, + Tsatw, + L + ); + + // Bubble departure frequency: + const scalarField fDep + ( + departureFreqModel_->fDeparture ( - departureFreqModel_->fDeparture - ( - liquid, - vapor, - patchi, - dDep_ - ) - ); + liquid, + vapor, + patchi, + dDep_ + ) + ); - // Convective thermal diffusivity for single phase - alphatConv_ = calcAlphat(alphatConv_); + // Convective thermal diffusivity for single phase + alphatConv_ = calcAlphat(alphatConv_); - // Convective heat transfer area for Sub-cool boiling - scalarField A1(this->size(), 0); + // Convective heat transfer area for Sub-cool boiling + scalarField A1(this->size(), 0); - const scalarField hewSn(hew.snGrad()); + const scalarField hewSn(hew.snGrad()); - scalarField alphaFilm(this->size(), 0); + scalarField alphaFilm(this->size(), 0); - // Use to identify regimes per face - labelField regimeTypes(A1.size(), -1); + // Use to identify regimes per face + labelField regimeTypes(A1.size(), -1); - forAll (*this, i) + forAll (*this, i) + { + if (Tw[i] > Tsatw[i]) { - if (Tw[i] > Tsatw[i]) + // Sub-cool boiling + if (Tw[i] < tDNB[i]) { // Sub-cool boiling - if (Tw[i] < tDNB[i]) - { - // Sub-cool boiling - regimeTypes[i] = regimeType::subcool; + regimeTypes[i] = regimeType::subcool; - Tl = (Tw - (Tplus_y250/Tplus)*(Tw - Tc)); - Tl = max(Tc - 40, Tl); + // Del Valle & Kenning (1985) + const scalar Ja = + rhow[i]*Cpw[i]*(Tsatw[i] - Tl[i]) + /(rhoVaporw[i]*L[i]); - // Area fractions: + const scalar Al = + fLiquid[i]*4.8*exp(min(-Ja/80, log(VGREAT))); - /* - // Del Valle & Kenning (1985) - const scalarField Ja + const scalar A2 + ( + min(pi*sqr(dDep_[i])*N[i]*Al/4, scalar(1)) + ); + + A1[i] = max(1 - A2, scalar(1e-4)); + + // Following Bowring(1962) + const scalar A2E + ( + min ( - rhoLiquidw*Cpw*(Tsatw - Tl)/(rhoVaporw*L) + pi*sqr(dDep_[i])*N[i]*Al/4, + scalar(5) + ) + ); + + // Volumetric mass source in the near wall cell due + // to the wall boiling + dmdt_[i] = + ( + (1 - relax)*dmdt_[i] + + relax*(1.0/6.0)*A2E*dDep_[i]*rhoVaporw[i] + * fDep[i]*AbyV_[i] ); - const scalarField Al - ( - fLiquid*4.8*exp(min(-Ja/80, log(VGREAT))) - ); - */ + // Volumetric source in the near wall cell due to + // the wall boiling + mDotL_[i] = dmdt_[i]*L[i]; - // More simple method to calculate area affected by - // bubbles - const scalar A2 + // Quenching heat transfer coefficient + const scalar hQ + ( + 2*(alphaw[i]*Cpw[i])*fDep[i] + *sqrt ( - min - ( - fLiquid[i]*pi*sqr(dDep_[i])*N[i]*K_/4, - scalar(1) - ) + (0.8/max(fDep[i], SMALL))/(pi*alphaw[i]/rhow[i]) + ) + ); + + // Quenching heat flux in Sub-cool boiling + qq_[i] = + ( + (1 - relax)*qq_[i] + + relax*A2*hQ*max(Tw[i] - Tl[i], scalar(0)) ); - A1[i] = max(1 - A2, 0.0); - - // Following Bowring(1962) - const scalar A2E - ( - min - ( - fLiquid[i]*pi*sqr(dDep_[i])*N[i], - scalar(5) - ) - ); - - // Volumetric mass source in the near wall cell due - // to the wall boiling - dmdt_[i] = - ( - (1 - relax_)*dmdt_[i] - + relax_*(1.0/6.0)*A2E*dDep_[i]*rhoVaporw[i] - * fDep[i]*AbyV_[i] - ); - - // Volumetric source in the near wall cell due to - // the wall boiling - mDotL_[i] = dmdt_[i]*L[i]; - - // Quenching heat transfer coefficient - const scalar hQ - ( - 2*(alphaw[i]*Cpw[i])*fDep[i] - *sqrt - ( - (0.8/fDep[i])/(pi*alphaw[i]/rhow[i]) - ) - ); - - // Quenching heat flux in Sub-cool boiling - qq_[i] = - ( - (1 - relax_)*qq_[i] - + relax_*A2*hQ*max(Tw[i] - Tl[i], scalar(0)) - ); - - this->operator[](i) = - ( - max - ( - A1[i]*alphatConv_[i] - + ( - (qq_[i] + mDotL_[i]/AbyV_[i]) - / max(hewSn[i], scalar(1e-16)) - ) - /max(liquidw[i], scalar(1e-8)), - 1e-8 - ) - ); - } - else if (Tw[i] > tDNB[i] && Tw[i] < TLeiden[i]) - { - // transient boiling - regimeTypes[i] = regimeType::transient; - - // No convective heat tranfer - alphatConv_[i] = 0.0; - - // transient boiling - dmdt_[i] = - fLiquid[i] - *( - relax_*Qtb[i]*AbyV_[i]/L[i] - + (1 - relax_)*dmdt_[i] - ); - - mDotL_[i] = dmdt_[i]*L[i]; - - // No quenching flux - qq_[i] = 0.0; - - this->operator[](i) = - ( - max - ( - ( - mDotL_[i]/AbyV_[i] - /max(hewSn[i], scalar(1e-16)) - )/max(liquidw[i], scalar(1e-8)), - 1e-8 - ) - ); - } - else if (Tw[i] > TLeiden[i]) - { - regimeTypes[i] = regimeType::film; // film boiling - - // No convective heat tranfer - alphatConv_[i] = 0.0; - - // Film boiling - dmdt_[i] = - fLiquid[i] - *( - relax_*htcFilmBoiling[i] - *max(Tw[i] - Tsatw[i], 0)*AbyV_[i]/L[i] - + (1 - relax_)*dmdt_[i] - ); - - - mDotL_[i] = dmdt_[i]*L[i]; - - // No quenching flux - qq_[i] = 0.0; - - alphaFilm[i] = - ( - mDotL_[i]/AbyV_[i]/max(hewSn[i], scalar(1e-16)) - ); - - // alphat is added alphal and multiplied by phase - // alphaFilm is lower than alphal. Then we substract - // alpha and divide by phase to get alphaFilm - this->operator[](i) = - ( - (alphaFilm[i]/max(liquidw[i], scalar(1e-8)) - - alphaw[i]) - ); - } - } - else - { - // Tw below Tsat. No boiling single phase convection - // Single phase - regimeTypes[i] = regimeType::nonBoiling; - A1[i] = 1.0; - qq_[i] = 0.0; - mDotL_[i] = 0.0; - dmdt_[i] = 0.0; - - // Turbulente thermal diffusivity for single phase. this->operator[](i) = ( max ( - fLiquid[i]*A1[i]*(alphatConv_[i]) - /max(liquidw[i], scalar(1e-8)), - 1e-16 + A1[i]*alphatConv_[i] + + ( + (qq_[i] + mDotL_[i]/AbyV_[i]) + / max(hewSn[i], scalar(1e-16)) + ) + /max(liquidw[i], scalar(1e-8)), + 1e-8 ) ); } + else if (Tw[i] > tDNB[i] && Tw[i] < TLeiden[i]) + { + // transient boiling + regimeTypes[i] = regimeType::transient; + + // No convective heat tranfer + alphatConv_[i] = 0.0; + + // transient boiling + dmdt_[i] = + fLiquid[i] + *( + relax*Qtb[i]*AbyV_[i]/L[i] + + (1 - relax)*dmdt_[i] + ); + + mDotL_[i] = dmdt_[i]*L[i]; + + + // No quenching flux + qq_[i] = 0.0; + + this->operator[](i) = + max + ( + ( + mDotL_[i]/AbyV_[i] + /max(hewSn[i], scalar(1e-16)) + )/max(liquidw[i], scalar(1e-8)), + 1e-8 + ); + + } + else if (Tw[i] > TLeiden[i]) + { + regimeTypes[i] = regimeType::film; // film boiling + + // No convective heat tranfer + alphatConv_[i] = 0.0; + + // Film boiling + dmdt_[i] = + fLiquid[i] + *( + relax*htcFilmBoiling[i] + *max(Tw[i] - Tsatw[i], 0) + *AbyV_[i]/L[i] + + (1 - relax)*dmdt_[i] + ); + + + mDotL_[i] = dmdt_[i]*L[i]; + + // No quenching flux + qq_[i] = 0.0; + + alphaFilm[i] = + ( + mDotL_[i]/AbyV_[i]/max(hewSn[i], scalar(1e-16)) + ); + + // alphat is added alphal and multiplied by phase + // alphaFilm in the coupled BC. We substract + // alpha and divide by phase to get a net alphaFilm + this->operator[](i) = + ( + alphaFilm[i]/max(liquidw[i], scalar(1e-8)) + - alphaw[i] + ); + } } - - scalarField TsupPrev(max((Tw - Tsatw), scalar(0))); - - // NOTE: lagging Tw update. - //const_cast(Tw).evaluate(); - scalarField TsupNew(max((Tw - Tsatw), scalar(0))); - - scalar maxErr(max(mag(TsupPrev - TsupNew))); - - if (debug) + else { - const scalarField qEff + // Tw below Tsat. No boiling single phase convection + // Single phase + regimeTypes[i] = regimeType::nonBoiling; + + qq_[i] = 0.0; + mDotL_[i] = 0.0; + dmdt_[i] = 0.0; + + // Turbulente thermal diffusivity for single phase. + this->operator[](i) = ( - liquidw*(*this + alphaw)*hew.snGrad() + max + ( + fLiquid[i]*(alphatConv_[i]) + /max(liquidw[i], scalar(1e-8)), + 1e-8 + ) + ); + } + } + + if (debug) + { + const scalarField qEff + ( + liquidw*(*this + alphaw)*hew.snGrad() + ); + + Info<< "alphat for liquid: " << nl << endl; + + Info<< " alphatl: " << gMin((*this)) << " - " + << gMax((*this)) << endl; + + Info<< " dmdt: " << gMin((dmdt_)) << " - " + << gMax((dmdt_)) << endl; + + Info<< " alphatlEff: " << gMin(liquidw*(*this + alphaw)) + << " - " << gMax(liquidw*(*this + alphaw)) << endl; + + scalar Qeff = gSum(qEff*patch().magSf()); + Info<< " Effective heat transfer rate to liquid: " << Qeff + << endl << nl; + + if (debug & 2) + { + /* + scalarField& faceRegimes = + const_cast + ( + patch().lookupPatchField + < + volScalarField, + scalar + >("faceRegimes") + ); + */ + + scalar nSubCool(0); + scalar nTransient(0); + scalar nFilm(0); + scalar nNonBoiling(0); + + scalarField nSubCools(this->size(), 0); + scalarField nTransients(this->size(), 0); + scalarField nFilms(this->size(), 0); + scalarField nNonBoilings(this->size(), 0); + + forAll (*this, i) + { + //faceRegimes[i] = regimeTypes[i]; + switch (regimeTypes[i]) + { + case regimeType::subcool: + nSubCool++; + nSubCools[i] = 1; + break; + + case regimeType::transient: + nTransient++; + nTransients[i] = 1; + break; + + case regimeType::film: + nFilm++; + nFilms[i] = 1; + break; + + case regimeType::nonBoiling: + nNonBoiling++; + nNonBoilings[i] = 1; + break; + } + } + + Info<< "Faces regime : " << nl << endl; + + Info<< " sub Cool faces : " << nSubCool << endl; + Info<< " transient faces : " << nTransient << endl; + Info<< " film faces : " << nFilm << endl; + Info<< " non-Boiling faces : " << nNonBoiling << endl; + Info<< " total faces : " << this->size() << endl << nl; + + const scalarField qc + ( + nNonBoilings*fLiquid*A1*(alphatConv_ + alphaw) + *hew.snGrad() + ); + /* + scalarField& qcFluid = + const_cast + ( + patch().lookupPatchField + < + volScalarField, + scalar + >("qc") + ); + + qcFluid = qc; + */ + + scalar Qc = gSum(qc*patch().magSf()); + Info<< " Convective heat transfer: " << Qc << endl; + + const scalarField qFilm + ( + relax*fLiquid*nFilms*htcFilmBoiling*(Tw - Tsatw) + ); + /* + scalarField& qFilmFluid = + const_cast + ( + patch().lookupPatchField + < + volScalarField, + scalar + >("qFilm") + ); + + qFilmFluid = qFilm; + */ + + /* + scalarField& htcFilm = + const_cast + ( + patch().lookupPatchField + < + volScalarField, + scalar + >("htcFilmBoiling") ); - Info<< "alphat for liquid: " << nl << endl; + htcFilm = htcFilmBoiling; + */ + scalar QFilm = gSum(qFilm*patch().magSf()); + Info<< " Film boiling heat transfer: " << QFilm << endl; - Info<< " alphatl: " << gMin((*this)) << " - " - << gMax((*this)) << endl; - - Info<< " dmdt: " << gMin((dmdt_)) << " - " - << gMax((dmdt_)) << endl; - - Info<< " alphatlEff: " << gMin(liquidw*(*this + alphaw)) - << " - " << gMax(liquidw*(*this + alphaw)) << endl; - - scalar Qeff = gSum(qEff*patch().magSf()); - Info<< " Effective heat transfer rate to liquid: " << Qeff - << endl << nl; - - if (debug & 2) - { - scalar nSubCool(0); - scalar nTransient(0); - scalar nFilm(0); - scalar nNonBoiling(0); - - scalarField nSubCools(this->size(), 0); - scalarField nTransients(this->size(), 0); - scalarField nFilms(this->size(), 0); - scalarField nNonBoilings(this->size(), 0); - - forAll (*this, i) - { - switch (regimeTypes[i]) - { - case regimeType::subcool: - nSubCool++; - nSubCools[i] = 1; - break; - - case regimeType::transient: - nTransient++; - nTransients[i] = 1; - break; - - case regimeType::film: - nFilm++; - nFilms[i] = 1; - break; - - case regimeType::nonBoiling: - nNonBoiling++; - nNonBoilings[i] = 1; - break; - } - } - - Info<< "Faces regime : " << nl << endl; - - Info<< " sub Cool faces : " << nSubCool << endl; - Info<< " transient faces : " << nTransient << endl; - Info<< " film faces : " << nFilm << endl; - Info<< " non-Boiling faces : " << nNonBoiling << endl; - Info<< " total faces : " << this->size() << endl << nl; - - const scalarField qc + Info<< " Htc Film Boiling coeff: " + << gMin(nFilms*htcFilmBoiling) + << " - " + << gMax(nFilms*htcFilmBoiling) << endl; + /* + scalarField& qtbFluid = + const_cast ( - nNonBoilings*fLiquid*A1*(alphatConv_ + alphaw) - *hew.snGrad() + patch().lookupPatchField + < + volScalarField, + scalar + >("qtb") + ); + qtbFluid = fLiquid*Qtb; + */ + + scalar Qtbtot = + gSum(fLiquid*nTransients*Qtb*patch().magSf()); + Info<< " Transient boiling heat transfer:" << Qtbtot + << endl; + + /* + scalarField& TdnbFluid = + const_cast + ( + patch().lookupPatchField + < + volScalarField, + scalar + >("Tdnb") + ); + TdnbFluid = tDNB; + */ + + Info<< " TDNB: " << gMin(tDNB) << " - " << gMax(tDNB) + << endl; + + const scalarField qSubCool + ( + fLiquid*nSubCools* + ( + A1*alphatConv_*hew.snGrad() + + qe() + qq() + ) + ); + + /* + scalarField& qSubCoolFluid = + const_cast + ( + patch().lookupPatchField + < + volScalarField, + scalar + >("qSubCool") ); - scalar Qc = gSum(qc*patch().magSf()); - Info<< " Convective heat transfer: " << Qc << endl; + qSubCoolFluid = qSubCool; + */ + + scalar QsubCool = gSum(qSubCool*patch().magSf()); - const scalarField qFilm - ( - relax_*fLiquid*nFilms*htcFilmBoiling*(Tw - Tsatw) - ); + Info<< " Sub Cool boiling heat transfer: " << QsubCool + << endl; - scalar QFilm = gSum(qFilm*patch().magSf()); - Info<< " Film boiling heat transfer: " << QFilm << endl; - - Info<< " Htc Film Boiling coeff: " - << gMin(nFilms*htcFilmBoiling) - << " - " - << gMax(nFilms*htcFilmBoiling) << endl; - - scalar Qtbtot = - gSum(fLiquid*nTransients*Qtb*patch().magSf()); - Info<< " Transient boiling heat transfer:" << Qtbtot - << endl; - - Info<< " TDNB: " << gMin(tDNB) << " - " << gMax(tDNB) - << endl; - - const scalarField qSubCool - ( - fLiquid*nSubCools* - ( - A1*alphatConv_*hew.snGrad() - + qe() + qq() - ) - ); - - scalar QsubCool = gSum(qSubCool*patch().magSf()); - - Info<< " Sub Cool boiling heat transfer: " << QsubCool - << endl; - - Info<< " N: " << gMin(nSubCools*N) << " - " - << gMax(nSubCools*N) << endl; - Info<< " dDep: " << gMin(nSubCools*dDep_) << " - " - << gMax(nSubCools*dDep_) << endl; - Info<< " fDep: " << gMin(nSubCools*fDep) << " - " - << gMax(nSubCools*fDep) << endl; - Info<< " A1: " << gMin(nSubCools*A1) << " - " - << gMax(nSubCools*A1) << endl; - - Info<< nl; - - } + Info<< " N: " << gMin(nSubCools*N) << " - " + << gMax(nSubCools*N) << endl; + Info<< " dDep: " << gMin(nSubCools*dDep_) << " - " + << gMax(nSubCools*dDep_) << endl; + Info<< " fDep: " << gMin(nSubCools*fDep) << " - " + << gMax(nSubCools*fDep) << endl; + Info<< " A1: " << gMin(nSubCools*A1) << " - " + << gMax(nSubCools*A1) << endl; + Info<< nl; } - - if (maxErr < 1e-1) - { - if (i > 0) - { - Info<< "Wall boiling wall function iterations: " - << i + 1 << endl; - } - break; - } - } - break; } + break; default: { FatalErrorInFunction @@ -1286,7 +1248,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const os.writeKeyword("phaseType") << phaseTypeNames_[phaseType_] << token::END_STATEMENT << nl; - os.writeKeyword("relax") << relax_ << token::END_STATEMENT << nl; + relax_->writeData(os); switch (phaseType_) { diff --git a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H index aab97c57d6..4a39dc5ce0 100644 --- a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H +++ b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H @@ -224,6 +224,8 @@ SourceFiles #ifndef compressible_alphatWallBoilingWallFunctionFvPatchScalarField_H #define compressible_alphatWallBoilingWallFunctionFvPatchScalarField_H +#include "Function1.H" + #include "alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H" #include "partitioningModel.H" #include "nucleationSiteModel.H" @@ -287,7 +289,7 @@ private: phaseType phaseType_; //- dmdt relaxationFactor - scalar relax_; + autoPtr> relax_; //- Patch face area by cell volume scalarField AbyV_; diff --git a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/Schroeder/Schroeder.C b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/Schroeder/Schroeder.C index 057cc7f908..128a7cdfd8 100644 --- a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/Schroeder/Schroeder.C +++ b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/Schroeder/Schroeder.C @@ -58,7 +58,7 @@ Foam::wallBoilingModels::TDNBModels::Schroeder::Schroeder ) : TDNBModel(), - kg_(dict.lookupOrDefault("kg", 5/3)) + kg_(dict.lookupOrDefault("kg", 1.666)) {} diff --git a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.C b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.C index 86d253b16e..8439d3f274 100644 --- a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.C +++ b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.C @@ -26,6 +26,9 @@ License #include "Bromley.H" #include "addToRunTimeSelectionTable.H" #include "uniformDimensionedFields.H" +#include "constants.H" + +using namespace Foam::constant; // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -55,7 +58,9 @@ Foam::wallBoilingModels::filmBoilingModels::Bromley::Bromley ) : filmBoilingModel(), - Cn_(dict.lookupOrDefault("Cn", 0.62)) + Cn_(dict.lookupOrDefault("Cn", 0.62)), + emissivity_(dict.lookupOrDefault("emissivity", 1)), + L_(dict.get("L")) {} @@ -84,7 +89,11 @@ Foam::wallBoilingModels::filmBoilingModels::Bromley::htcFilmBoil const uniformDimensionedVectorField& g = liquid.mesh().time().lookupObject("g"); - const scalarField rhoVapor(vapor.thermo().rho(patchi)); + const fvPatchScalarField& rhoVaporw + ( + vapor.thermo().rho()().boundaryField()[patchi] + ); + const scalarField rhoLiq(liquid.thermo().rho(patchi)); const scalarField kappaVapor(vapor.kappa(patchi)); @@ -93,17 +102,23 @@ Foam::wallBoilingModels::filmBoilingModels::Bromley::htcFilmBoil const scalarField& CpVapor = Cp.boundaryField()[patchi]; const scalarField muVapor(vapor.mu(patchi)); - const scalarField dbVapor(vapor.d()().boundaryField()[patchi]); + //const scalarField dbVapor(vapor.d()().boundaryField()[patchi]); + + const scalarField htcRad + ( + emissivity_*physicoChemical::sigma.value()*(pow4(Tw) - pow4(Tsatw)) + / max((Tw - Tsatw), scalar(1e-4)) + ); return Cn_*pow ( pow3(kappaVapor) - *rhoVapor*(rhoLiq - rhoVapor)*mag(g.value()) + *rhoVaporw*(rhoLiq - rhoVaporw)*mag(g.value()) *(L + 0.4*CpVapor*max((Tw-Tsatw), scalar(0))) - /(dbVapor*muVapor*max((Tw-Tsatw), scalar(1e-4))), + /(L_*muVapor*max((Tw-Tsatw), scalar(1e-4))), 0.25 - ); + ) + 0.75*htcRad; } @@ -114,6 +129,7 @@ void Foam::wallBoilingModels::filmBoilingModels::Bromley::write { filmBoilingModel::write(os); os.writeKeyword("Cn") << Cn_ << token::END_STATEMENT << nl; + os.writeKeyword("emissivity") << emissivity_ << token::END_STATEMENT << nl; } diff --git a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.H b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.H index 782c7e4aab..8427e752aa 100644 --- a/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.H +++ b/src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.H @@ -66,6 +66,12 @@ class Bromley //- Coefficient for nucleation site density scalar Cn_; + //- Wall emissivity + scalar emissivity_; + + //- Characteristic lenght-scale + scalar L_; + public: //- Runtime type information diff --git a/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.gas b/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.gas index e502ff6f6e..b9829fa82b 100644 --- a/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.gas +++ b/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.gas @@ -43,23 +43,13 @@ boundaryField Cmu 0.09; kappa 0.41; E 9.8; + relax constant 1; partitioningModel { type Lavieville; alphaCrit 0.2; } - - LeidenfrostModel - { - type Spiegler; - Tcrit 647; - } - - filmBoilingModel - { - type Bromley; - } - value uniform 0; + value uniform 1e-8; } maxY { diff --git a/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.liquid b/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.liquid index a2a87eeec3..633f99ec62 100644 --- a/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.liquid +++ b/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.liquid @@ -43,6 +43,7 @@ boundaryField Cmu 0.09; kappa 0.41; E 9.8; + relax constant 1; dmdt uniform 0; partitioningModel { @@ -88,6 +89,7 @@ boundaryField filmBoilingModel { type Bromley; + L 1e-3; } value uniform 1e-8; }