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