diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C index 4619a4a76a..34e2bdc385 100644 --- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C +++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C @@ -74,6 +74,7 @@ Foam::surfaceFilmModels::standardPhaseChange::standardPhaseChange ) : phaseChangeModel(typeName, owner, dict), + Tb_(readScalar(coeffs_.lookup("Tb"))), deltaMin_(readScalar(coeffs_.lookup("deltaMin"))), L_(readScalar(coeffs_.lookup("L"))) {} @@ -116,48 +117,63 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct const scalarField hInf = film.htcs().h(); const scalarField hFilm = film.htcw().h(); const vectorField dU = film.UPrimary() - film.Us(); - const scalarField availableMass = delta*rho*magSf; + const scalarField availableMass = (delta - deltaMin_)*rho*magSf; - // Reynolds number - const scalarField Re = rho*mag(dU)*L_/mu; - - // molecular weight of vapour - const scalar Wvap = thermo.carrier().W(vapId); - - // molecular weight of liquid - const scalar Wliq = liq.W(); forAll(dMass, cellI) { if (delta[cellI] > deltaMin_) { - // cell pressure + // cell pressure [Pa] const scalar pc = pInf[cellI]; - // saturation pressure - const scalar pSat = liq.pv(pc, Ts[cellI]); + // local temperature - impose lower limit of 200 K for stability + const scalar Tloc = min(liq.Tc(), max(200.0, Ts[cellI])); - // latent heat - const scalar hVap = liq.hl(pc, Ts[cellI]); + // saturation pressure [Pa] + const scalar pSat = liq.pv(pc, Tloc); + + // latent heat [J/kg] + const scalar hVap = liq.hl(pc, Tloc); // calculate mass transfer - if (pSat > pc) + if (pSat >= pc) { // boiling const scalar qDotInf = hInf[cellI]*(TInf[cellI] - T[cellI]); const scalar qDotFilm = hFilm[cellI]*(T[cellI] - Tw[cellI]); - dMass[cellI] = dt*magSf[cellI]/hVap*(qDotInf - qDotFilm); + + const scalar cp = liq.cp(pc, Tloc); + const scalar qCorr = availableMass[cellI]*cp*(T[cellI] - Tb_); + dMass[cellI] = + dt*magSf[cellI]/hVap*(qDotInf + qDotFilm) + + qCorr/hVap; + } else { + // Reynolds number + const scalarField Re = rho*mag(dU)*L_/mu; + + // molecular weight of vapour [kg/kmol] + const scalar Wvap = thermo.carrier().W(vapId); + + // molecular weight of liquid [kg/kmol] + const scalar Wliq = liq.W(); + // vapour mass fraction at interface const scalar Ys = Wliq*pSat/(Wliq*pSat + Wvap*(pc - pSat)); - // bulk gas average density - const scalar rhoAve = pc/(specie::RR*Ts[cellI]); + // 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, Ts[cellI]); + const scalar Dab = liq.D(pc, Tloc); // Schmidt number const scalar Sc = mu[cellI]/(rho[cellI]*(Dab + ROOTVSMALL)); @@ -170,7 +186,7 @@ void Foam::surfaceFilmModels::standardPhaseChange::correct // add mass contribution to source dMass[cellI] = - dt*magSf[cellI]*rhoAve*hm*(Ys - YInf[cellI])/(1.0 - Ys); + dt*magSf[cellI]*rhoInf*hm*(Ys - YInf[cellI])/(1.0 - Ys); } dMass[cellI] = min(availableMass[cellI], max(0.0, dMass[cellI])); diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H index 726f36e2cc..d1633a6f93 100644 --- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H +++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H @@ -67,10 +67,13 @@ protected: // Protected data + //- Boiling temperature / [K] + const scalar Tb_; + //- Minimum film height for model to be active const scalar deltaMin_; - //- Length scalae / [m] + //- Length scale / [m] const scalar L_; diff --git a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.C b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.C index bbba0b5371..b7013da060 100644 --- a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.C +++ b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.C @@ -222,6 +222,10 @@ transferPrimaryRegionFields() rhoSp_.field() /= magSf_*deltaT; USp_.field() /= magSf_*deltaT; pSp_.field() /= magSf_*deltaT; + + // reset transfer to primary fields + massForPrimary_ == dimensionedScalar("zero", dimMass, 0.0); + diametersForPrimary_ == dimensionedScalar("zero", dimLength, -1.0); } @@ -271,9 +275,6 @@ Foam::surfaceFilmModels::kinematicSingleLayer::pp() void Foam::surfaceFilmModels::kinematicSingleLayer::correctDetachedFilm() { - massForPrimary_ == dimensionedScalar("zero", dimMass, 0.0); - diametersForPrimary_ == dimensionedScalar("zero", dimLength, -1.0); - const scalarField gNorm = this->gNorm(); forAll(gNorm, i) @@ -308,8 +309,7 @@ void Foam::surfaceFilmModels::kinematicSingleLayer::updateSubmodels() // Update source fields const dimensionedScalar deltaT = time_.deltaT(); - rhoSp_ -= massForPrimary_/magSf_/deltaT; - USp_ -= massForPrimary_*U_/magSf_/deltaT; + rhoSp_ -= (massForPrimary_ + massPhaseChangeForPrimary_)/magSf_/deltaT; } @@ -457,6 +457,12 @@ Foam::surfaceFilmModels::kinematicSingleLayer::solveMomentum USp_ + tau(U_) + fvc::grad(sigma_) + - fvm::Sp + ( + (massForPrimary_ + massPhaseChangeForPrimary_) + /magSf_/time_.deltaT(), + U_ + ) ); fvVectorMatrix& UEqn = tUEqn(); @@ -1239,7 +1245,7 @@ Foam::surfaceFilmModels::kinematicSingleLayer::Srho() const const label filmPatchI = filmBottomPatchIDs_[i]; scalarField patchMass = - massPhaseChangeForPrimary().boundaryField()[filmPatchI]; + massPhaseChangeForPrimary_.boundaryField()[filmPatchI]; distMap.distribute(patchMass); diff --git a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C index 3565cf7c12..103ab3bf6d 100644 --- a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C +++ b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C @@ -138,10 +138,8 @@ void Foam::surfaceFilmModels::noFilm::addSources const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::U() const { - FatalErrorIn - ( - "const volScalarField& noFilm::U() const" - ) << "U field not available for " << type() << abort(FatalError); + FatalErrorIn("const volScalarField& noFilm::U() const") + << "U field not available for " << type() << abort(FatalError); return volVectorField::null(); } @@ -149,10 +147,8 @@ const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::U() const const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::Us() const { - FatalErrorIn - ( - "const volScalarField& noFilm::Us() const" - ) << "Us field not available for " << type() << abort(FatalError); + FatalErrorIn("const volScalarField& noFilm::Us() const") + << "Us field not available for " << type() << abort(FatalError); return volVectorField::null(); } @@ -160,10 +156,8 @@ const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::Us() const const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::Uw() const { - FatalErrorIn - ( - "const volScalarField& noFilm::Uw() const" - ) << "Uw field not available for " << type() << abort(FatalError); + FatalErrorIn("const volScalarField& noFilm::Uw() const") + << "Uw field not available for " << type() << abort(FatalError); return volVectorField::null(); } @@ -171,10 +165,8 @@ const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::Uw() const const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::rho() const { - FatalErrorIn - ( - "const volScalarField& noFilm::rho() const" - ) << "rho field not available for " << type() << abort(FatalError); + FatalErrorIn("const volScalarField& noFilm::rho() const") + << "rho field not available for " << type() << abort(FatalError); return volScalarField::null(); } @@ -182,10 +174,8 @@ const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::rho() const const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::T() const { - FatalErrorIn - ( - "const Foam::volScalarField& Foam::noFilm::T() const" - ) << "T field not available for " << type() << abort(FatalError); + FatalErrorIn("const Foam::volScalarField& Foam::noFilm::T() const") + << "T field not available for " << type() << abort(FatalError); return volScalarField::null(); } @@ -193,10 +183,8 @@ const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::T() const const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::Ts() const { - FatalErrorIn - ( - "const Foam::volScalarField& Foam::noFilm::Ts() const" - ) << "Ts field not available for " << type() << abort(FatalError); + FatalErrorIn("const Foam::volScalarField& Foam::noFilm::Ts() const") + << "Ts field not available for " << type() << abort(FatalError); return volScalarField::null(); } @@ -204,10 +192,8 @@ const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::Ts() const const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::Tw() const { - FatalErrorIn - ( - "const Foam::volScalarField& Foam::noFilm::Tw() const" - ) << "Tw field not available for " << type() << abort(FatalError); + FatalErrorIn("const Foam::volScalarField& Foam::noFilm::Tw() const") + << "Tw field not available for " << type() << abort(FatalError); return volScalarField::null(); } @@ -215,10 +201,8 @@ const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::Tw() const const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::cp() const { - FatalErrorIn - ( - "const volScalarField& noFilm::cp() const" - ) << "cp field not available for " << type() << abort(FatalError); + FatalErrorIn("const volScalarField& noFilm::cp() const") + << "cp field not available for " << type() << abort(FatalError); return volScalarField::null(); } @@ -226,10 +210,8 @@ const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::cp() const const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::kappa() const { - FatalErrorIn - ( - "const volScalarField& noFilm::kappa() const" - ) << "kappa field not available for " << type() << abort(FatalError); + FatalErrorIn("const volScalarField& noFilm::kappa() const") + << "kappa field not available for " << type() << abort(FatalError); return volScalarField::null(); } @@ -238,10 +220,8 @@ const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::kappa() const const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::massForPrimary() const { - FatalErrorIn - ( - "const volScalarField& noFilm::massForPrimary() const" - ) << "massForPrimary field not available for " << type() + FatalErrorIn("const volScalarField& noFilm::massForPrimary() const") + << "massForPrimary field not available for " << type() << abort(FatalError); return volScalarField::null(); @@ -251,10 +231,8 @@ Foam::surfaceFilmModels::noFilm::massForPrimary() const const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::diametersForPrimary() const { - FatalErrorIn - ( - "const volScalarField& noFilm::diametersForPrimary() const" - ) << "diametersForPrimary field not available for " << type() + FatalErrorIn("const volScalarField& noFilm::diametersForPrimary() const") + << "diametersForPrimary field not available for " << type() << abort(FatalError); return volScalarField::null(); @@ -267,7 +245,7 @@ Foam::surfaceFilmModels::noFilm::massPhaseChangeForPrimary() const FatalErrorIn ( "const volScalarField& noFilm::massPhaseChangeForPrimary() const" - ) << "massPhaseChangeForPrimary field not available for " << type() + ) << "massPhaseChange field not available for " << type() << abort(FatalError); return volScalarField::null(); diff --git a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C index 480759783c..57f5a95ac5 100644 --- a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C +++ b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C @@ -175,8 +175,8 @@ void Foam::surfaceFilmModels::thermoSingleLayer::updateSubmodels() htcw_->correct(); // Update phase change - massPhaseChangeForPrimary_ == dimensionedScalar("zero", dimMass, 0.0); - energyPhaseChangeForPrimary_ == dimensionedScalar("zero", dimEnergy, 0.0); + massPhaseChangeForPrimary_.internalField() = 0.0; + energyPhaseChangeForPrimary_.internalField() = 0.0; phaseChange_->correct ( @@ -191,11 +191,7 @@ void Foam::surfaceFilmModels::thermoSingleLayer::updateSubmodels() kinematicSingleLayer::updateSubmodels(); // Update source fields - const dimensionedScalar deltaT = time_.deltaT(); - rhoSp_ -= massPhaseChangeForPrimary_/magSf_/deltaT; - USp_ -= massPhaseChangeForPrimary_*U_/magSf_/deltaT; - hsSp_ -= - (massForPrimary_*hs_ + energyPhaseChangeForPrimary_)/magSf_/deltaT; + hsSp_ -= energyPhaseChangeForPrimary_/magSf_/time_.deltaT(); } @@ -208,7 +204,7 @@ Foam::tmp Foam::surfaceFilmModels::thermoSingleLayer::q return ( - - fvm::Sp(htcs_->h()/cp_, hs) - htcs_->h()*(Tstd - Ts_) + - fvm::Sp(htcs_->h()/cp_, hs) - htcs_->h()*(Tstd - TPrimary_) - fvm::Sp(htcw_->h()/cp_, hs) - htcw_->h()*(Tstd - Tw_) ); } @@ -228,8 +224,9 @@ void Foam::surfaceFilmModels::thermoSingleLayer::solveEnergy() fvm::ddt(deltaRho_, hs_) + fvm::div(phi_, hs_) == - hsSp_ + fvm::Sp(hsSp_/hs_, hs_) + q(hs_) + - fvm::Sp(massForPrimary_/magSf_/time_.deltaT(), hs_) ); correctThermoFields(); @@ -627,7 +624,7 @@ Foam::surfaceFilmModels::thermoSingleLayer::Srho(const label i) const const directMappedWallPolyPatch& wpp = refCast ( - mesh_.boundaryMesh()[primaryPatchI] + mesh_.boundaryMesh()[primaryPatchI] ); const mapDistribute& distMap = wpp.map(); @@ -635,7 +632,7 @@ Foam::surfaceFilmModels::thermoSingleLayer::Srho(const label i) const const label filmPatchI = filmBottomPatchIDs_[i]; scalarField patchMass = - massPhaseChangeForPrimary().boundaryField()[filmPatchI]; + massPhaseChangeForPrimary_.boundaryField()[filmPatchI]; distMap.distribute(patchMass); @@ -672,7 +669,7 @@ Foam::surfaceFilmModels::thermoSingleLayer::Sh() const dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0) ) ); - +/* scalarField& Sh = tSh(); const scalarField& V = mesh_.V(); const scalar dt = time_.deltaTValue(); @@ -683,28 +680,25 @@ Foam::surfaceFilmModels::thermoSingleLayer::Sh() const const directMappedWallPolyPatch& wpp = refCast ( - mesh_.boundaryMesh()[primaryPatchI] + mesh_.boundaryMesh()[primaryPatchI] ); const mapDistribute& distMap = wpp.map(); const label filmPatchI = filmBottomPatchIDs_[i]; - scalarField patchMass = - massPhaseChangeForPrimary().boundaryField()[filmPatchI]; - distMap.distribute(patchMass); - - scalarField patchEnthalpy = hs_.boundaryField()[filmPatchI]; - distMap.distribute(patchEnthalpy); + scalarField patchEnergy = + energyPhaseChangeForPrimary_.boundaryField()[filmPatchI]; + distMap.distribute(patchEnergy); const unallocLabelList& cells = wpp.faceCells(); forAll(patchMass, j) { - Sh[cells[j]] += patchMass[j]*patchEnthalpy[j]/(V[cells[j]]*dt); + Sh[cells[j]] += patchEnergy[j]/(V[cells[j]]*dt); } } - +*/ return tSh; }