From a2ad7167615523538f13ac99f613b33c007407c9 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Sun, 30 Apr 2023 10:19:25 +0100 Subject: [PATCH] isothermalFilm/fvModels/filmCloudTransfer: New models to transfer Lagrangian parcels to film The parcel transfer occurs from the cloudFilmTransfer surfaceFilmModel specified in the region constant//cloudProperties dictionary: . . . libs ("libfilmCloudTransfer.so"); . . . surfaceFilmModel cloudFilmTransfer; and the film filmCloudTransfer specified in the region constant//fvModels dictionary: . . . filmCloudTransfer { type filmCloudTransfer; libs ("libfilmCloudTransfer.so"); } For an example of cloud->film->VoF transfer see the tutorials/modules/multiRegion/film/cylinder tutorial case. Note that parcel transfer from film to Lagrangian cloud is not yet supported, this will be added soon. --- .../modules/film/thermophysicalPredictor.C | 4 + .../solvers/modules/isothermalFilm/Allwmake | 1 + .../modules/isothermalFilm/correctAlpha.C | 4 +- .../fvModels/filmCloudTransfer/Make/files | 9 + .../fvModels/filmCloudTransfer/Make/options | 22 + .../CloudFilmTransfer/CloudFilmTransfer.C | 743 ++++++++++++++++++ .../CloudFilmTransfer/CloudFilmTransfer.H | 319 ++++++++ .../CloudFilmTransfer/CloudFilmTransferBase.C | 48 ++ ...eactingMultiphaseParcelCloudFilmTransfer.C | 33 + .../makeReactingParcelCloudFilmTransfer.C | 33 + .../makeSprayParcelCloudFilmTransfer.C | 33 + .../makeThermoParcelCloudFilmTransfer.C | 33 + .../filmCloudTransfer/filmCloudTransfer.C | 284 +++++++ .../filmCloudTransfer/filmCloudTransfer.H | 212 +++++ .../fvModels/filmVoFTransfer/Make/files | 4 +- .../fvModels/filmVoFTransfer/Make/options | 2 +- ...{VoFtoFilmTransfer.C => VoFFilmTransfer.C} | 82 +- ...{VoFtoFilmTransfer.H => VoFFilmTransfer.H} | 39 +- ...{filmToVoFTransfer.C => filmVoFTransfer.C} | 86 +- ...{filmToVoFTransfer.H => filmVoFTransfer.H} | 37 +- .../modules/isothermalFilm/isothermalFilm.C | 19 +- .../modules/isothermalFilm/isothermalFilm.H | 13 +- .../isothermalFilm/momentumPredictor.C | 6 +- .../modules/isothermalFilm/predictAlpha.C | 20 +- src/OpenFOAM/dimensionSet/dimensionSets.C | 1 + src/OpenFOAM/dimensionSet/dimensionSets.H | 3 +- .../momentumSurfaceFilm/momentumSurfaceFilm.C | 20 +- .../film/VoFToFilm/constant/VoF/fvModels | 4 +- .../film/VoFToFilm/constant/film/fvModels | 4 +- .../multiRegion/film/cylinderVoF/0/VoF/T | 42 + .../multiRegion/film/cylinderVoF/0/VoF/T.air | 46 ++ .../film/cylinderVoF/0/VoF/T.liquid | 46 ++ .../multiRegion/film/cylinderVoF/0/VoF/U | 42 + .../film/cylinderVoF/0/VoF/alpha.liquid | 39 + .../multiRegion/film/cylinderVoF/0/VoF/p | 42 + .../multiRegion/film/cylinderVoF/0/VoF/p_rgh | 40 + .../multiRegion/film/cylinderVoF/0/film/T | 46 ++ .../multiRegion/film/cylinderVoF/0/film/U | 47 ++ .../multiRegion/film/cylinderVoF/0/film/delta | 53 ++ .../multiRegion/film/cylinderVoF/0/film/p | 45 ++ .../multiRegion/film/cylinderVoF/Allrun | 32 + .../cylinderVoF/constant/VoF/cloudProperties | 133 ++++ .../film/cylinderVoF/constant/VoF/fvModels | 42 + .../film/cylinderVoF/constant/VoF/g | 21 + .../constant/VoF/momentumTransport | 20 + .../constant/VoF/parcelInjectionProperties | 24 + .../cylinderVoF/constant/VoF/phaseProperties | 21 + .../constant/VoF/physicalProperties.air | 59 ++ .../constant/VoF/physicalProperties.liquid | 31 + .../cylinderVoF/constant/VoF/speciesThermo | 72 ++ .../constant/VoF/surfaceFilmProperties | 76 ++ .../film/cylinderVoF/constant/film/fvModels | 36 + .../film/cylinderVoF/constant/film/g | 21 + .../constant/film/momentumTransport | 20 + .../constant/film/physicalProperties | 39 + .../film/cylinderVoF/system/VoF/blockMeshDict | 143 ++++ .../cylinderVoF/system/VoF/decomposeParDict | 47 ++ .../film/cylinderVoF/system/VoF/fvSchemes | 58 ++ .../film/cylinderVoF/system/VoF/fvSolution | 89 +++ .../film/cylinderVoF/system/controlDict | 60 ++ .../film/cylinderVoF/system/extrudeMeshDict | 36 + .../cylinderVoF/system/film/blockMeshDict | 143 ++++ .../film/cylinderVoF/system/film/fvSchemes | 47 ++ .../film/cylinderVoF/system/film/fvSolution | 55 ++ .../film/cylinderVoF/system/fvSchemes | 30 + .../film/cylinderVoF/system/fvSolution | 21 + .../film/rivuletBox/system/film/fvSchemes | 4 +- 67 files changed, 3834 insertions(+), 152 deletions(-) create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/Make/files create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/Make/options create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.C create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.H create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransferBase.C create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/makeReactingMultiphaseParcelCloudFilmTransfer.C create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/makeReactingParcelCloudFilmTransfer.C create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/makeSprayParcelCloudFilmTransfer.C create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/makeThermoParcelCloudFilmTransfer.C create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.C create mode 100644 applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/filmCloudTransfer.H rename applications/solvers/modules/isothermalFilm/fvModels/filmVoFTransfer/{VoFtoFilmTransfer.C => VoFFilmTransfer.C} (81%) rename applications/solvers/modules/isothermalFilm/fvModels/filmVoFTransfer/{VoFtoFilmTransfer.H => VoFFilmTransfer.H} (90%) rename applications/solvers/modules/isothermalFilm/fvModels/filmVoFTransfer/{filmToVoFTransfer.C => filmVoFTransfer.C} (78%) rename applications/solvers/modules/isothermalFilm/fvModels/filmVoFTransfer/{filmToVoFTransfer.H => filmVoFTransfer.H} (89%) create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/VoF/T create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/VoF/T.air create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/VoF/T.liquid create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/VoF/U create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/VoF/alpha.liquid create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/VoF/p create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/VoF/p_rgh create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/film/T create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/film/U create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/film/delta create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/0/film/p create mode 100755 tutorials/modules/multiRegion/film/cylinderVoF/Allrun create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/cloudProperties create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/fvModels create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/g create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/momentumTransport create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/parcelInjectionProperties create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/phaseProperties create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/physicalProperties.air create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/physicalProperties.liquid create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/speciesThermo create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/VoF/surfaceFilmProperties create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/film/fvModels create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/film/g create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/film/momentumTransport create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/constant/film/physicalProperties create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/VoF/blockMeshDict create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/VoF/decomposeParDict create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/VoF/fvSchemes create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/VoF/fvSolution create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/controlDict create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/extrudeMeshDict create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/film/blockMeshDict create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/film/fvSchemes create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/film/fvSolution create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/fvSchemes create mode 100644 tutorials/modules/multiRegion/film/cylinderVoF/system/fvSolution diff --git a/applications/solvers/modules/film/thermophysicalPredictor.C b/applications/solvers/modules/film/thermophysicalPredictor.C index 03a1bd601a..2942d01e80 100644 --- a/applications/solvers/modules/film/thermophysicalPredictor.C +++ b/applications/solvers/modules/film/thermophysicalPredictor.C @@ -36,6 +36,7 @@ void Foam::solvers::film::thermophysicalPredictor() fvScalarMatrix heEqn ( fvm::ddt(alpha, rho, he) + fvm::div(alphaRhoPhi, he) + - fvm::Sp(contErr(), he) + thermophysicalTransport->divq(he) == fvModels().source(alpha, rho, he) @@ -50,6 +51,9 @@ void Foam::solvers::film::thermophysicalPredictor() fvConstraints().constrain(he); thermo_.correct(); + + Info<< max(alpha) << " " << min(alpha) << endl; + Info<< max(thermo.T()) << " " << min(thermo.T()) << endl; } diff --git a/applications/solvers/modules/isothermalFilm/Allwmake b/applications/solvers/modules/isothermalFilm/Allwmake index ebeddb9bb4..a22fc08c01 100755 --- a/applications/solvers/modules/isothermalFilm/Allwmake +++ b/applications/solvers/modules/isothermalFilm/Allwmake @@ -7,5 +7,6 @@ cd ${0%/*} || exit 1 # Run from this directory wmake $targetType filmCompressibleMomentumTransportModels wmake $targetType wmake $targetType fvModels/filmVoFTransfer +wmake $targetType fvModels/filmCloudTransfer #------------------------------------------------------------------------------ diff --git a/applications/solvers/modules/isothermalFilm/correctAlpha.C b/applications/solvers/modules/isothermalFilm/correctAlpha.C index ec6af6b9b4..0b7121e9e5 100644 --- a/applications/solvers/modules/isothermalFilm/correctAlpha.C +++ b/applications/solvers/modules/isothermalFilm/correctAlpha.C @@ -66,7 +66,7 @@ void Foam::solvers::isothermalFilm::correctAlpha() const surfaceScalarField alphaf ( - constrainedField(max(fvc::interpolate(alpha), scalar(0))) + constrainedField(fvc::interpolate(alpha)) ); const surfaceScalarField alpharAUf @@ -76,7 +76,7 @@ void Foam::solvers::isothermalFilm::correctAlpha() const surfaceScalarField phig("phig", phip + pbByAlphaGradRhof*alphaf); - phi = constrainPhiHbyA(fvc::flux(HbyA) - alpharAUf*phig, U, alpha); + phi = constrainedField(fvc::flux(HbyA) - alpharAUf*phig); const surfaceScalarField phid("phid", rhof*phi); diff --git a/applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/Make/files b/applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/Make/files new file mode 100644 index 0000000000..5ce92efe65 --- /dev/null +++ b/applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/Make/files @@ -0,0 +1,9 @@ +filmCloudTransfer.C + +cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransferBase.C +cloudFilmTransfer/makeThermoParcelCloudFilmTransfer.C +cloudFilmTransfer/makeReactingParcelCloudFilmTransfer.C +cloudFilmTransfer/makeReactingMultiphaseParcelCloudFilmTransfer.C +cloudFilmTransfer/makeSprayParcelCloudFilmTransfer.C + +LIB = $(FOAM_LIBBIN)/libfilmCloudTransfer diff --git a/applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/Make/options b/applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/Make/options new file mode 100644 index 0000000000..c955a4c391 --- /dev/null +++ b/applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/Make/options @@ -0,0 +1,22 @@ +EXE_INC = \ + -I$(FOAM_SOLVERS)/modules/isothermalFilm/lnInclude \ + -I$(FOAM_SOLVERS)/modules/isothermalFilm/filmCompressibleMomentumTransportModels/lnInclude \ + -I$(LIB_SRC)/physicalProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \ + -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \ + \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/lagrangian/parcel/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \ + \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + +LIB_LIBS = \ + -lisothermalFilm \ + -llagrangian \ + -llagrangianParcel diff --git a/applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.C b/applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.C new file mode 100644 index 0000000000..d386b5d5d3 --- /dev/null +++ b/applications/solvers/modules/isothermalFilm/fvModels/filmCloudTransfer/cloudFilmTransfer/CloudFilmTransfer/CloudFilmTransfer.C @@ -0,0 +1,743 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "CloudFilmTransfer.H" +#include "filmCloudTransfer.H" +#include "mappedPatchBase.H" +#include "ThermoCloud.H" +#include "meshTools.H" +#include "mathematicalConstants.H" +#include "addToRunTimeSelectionTable.H" + +using namespace Foam::constant::mathematical; + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +template +Foam::vector Foam::CloudFilmTransfer::splashDirection +( + const vector& tanVec1, + const vector& tanVec2, + const vector& nf +) const +{ + // Azimuthal angle [rad] + const scalar phiSi = twoPi*rndGen_.sample01(); + + // Ejection angle [rad] + const scalar thetaSi = pi/180.0*(rndGen_.sample01()*(50 - 5) + 5); + + // Direction vector of new parcel + const scalar alpha = sin(thetaSi); + const scalar dcorr = cos(thetaSi); + const vector normal = alpha*(tanVec1*cos(phiSi) + tanVec2*sin(phiSi)); + vector dirVec = dcorr*nf; + dirVec += normal; + + return dirVec/mag(dirVec); +} + + +template +void Foam::CloudFilmTransfer::absorbInteraction +( + fv::filmCloudTransfer& filmCloudTransfer, + const parcelType& p, + const polyPatch& pp, + const label facei, + const scalar mass, + bool& keepParticle +) +{ + if (debug) + { + Info<< "Parcel " << p.origId() << " absorbInteraction" << endl; + } + + const fluidThermo& carrierThermo = + static_cast&>(this->owner()) + .carrierThermo(); + + const parcelThermo& thermo = + static_cast&>(this->owner()).thermo(); + + // Patch face normal + const vector& nf = pp.faceNormals()[facei]; + + // Patch velocity + const vector& Up = this->owner().U().boundaryField()[pp.index()][facei]; + + // Relative parcel velocity + const vector Urel = p.U() - Up; + + // Parcel normal velocity + const vector Un = nf*(Urel & nf); + + // Parcel tangential velocity + const vector Ut = Urel - Un; + + const liquidProperties& liq = thermo.liquids().properties()[0]; + + // Local pressure + const scalar pc = carrierThermo.p()[p.cell()]; + + filmCloudTransfer.parcelFromCloud + ( + facei, + mass, // mass + mass*Ut, // tangential momentum + mass*mag(Un), // impingement pressure + mass*liq.Hs(pc, p.T()) // energy + ); + + this->nParcelsTransferred()++; + + keepParticle = false; +} + + +template +void Foam::CloudFilmTransfer::bounceInteraction +( + parcelType& p, + const polyPatch& pp, + const label facei, + bool& keepParticle +) const +{ + if (debug) + { + Info<< "Parcel " << p.origId() << " bounceInteraction" << endl; + } + + // Patch face normal + const vector& nf = pp.faceNormals()[facei]; + + // Patch velocity + const vector& Up = this->owner().U().boundaryField()[pp.index()][facei]; + + // Relative parcel velocity + const vector Urel = p.U() - Up; + + // Flip parcel normal velocity component + p.U() -= 2.0*nf*(Urel & nf); + + keepParticle = true; +} + + +template +void Foam::CloudFilmTransfer::drySplashInteraction +( + fv::filmCloudTransfer& filmCloudTransfer, + const parcelType& p, + const polyPatch& pp, + const label facei, + bool& keepParticle +) +{ + if (debug) + { + Info<< "Parcel " << p.origId() << " drySplashInteraction" << endl; + } + + const fluidThermo& carrierThermo = + static_cast&>(this->owner()) + .carrierThermo(); + + const parcelThermo& thermo = + static_cast&>(this->owner()).thermo(); + + const liquidProperties& liq = thermo.liquids().properties()[0]; + + // Patch face velocity and normal + const vector& Up = this->owner().U().boundaryField()[pp.index()][facei]; + const vector& nf = pp.faceNormals()[facei]; + + // Local pressure + const scalar pc = carrierThermo.p()[p.cell()]; + + // Retrieve parcel properties + const scalar m = p.mass()*p.nParticle(); + const scalar rho = p.rho(); + const scalar d = p.d(); + const scalar sigma = liq.sigma(pc, p.T()); + const scalar mu = liq.mu(pc, p.T()); + const vector Urel = p.U() - Up; + const vector Un = nf*(Urel & nf); + + // Laplace number + const scalar La = rho*sigma*d/sqr(mu); + + // Weber number + const scalar We = rho*magSqr(Un)*d/sigma; + + // Critical Weber number + const scalar Wec = Adry_*pow(La, -0.183); + + if (We < Wec) // Adhesion - assume absorb + { + absorbInteraction(filmCloudTransfer, p, pp, facei, m, keepParticle); + } + else // Splash + { + // Ratio of incident mass to splashing mass + const scalar mRatio = 0.2 + 0.6*rndGen_.sample01(); + splashInteraction + ( + filmCloudTransfer, + p, + pp, + facei, + mRatio, + We, + Wec, + sigma, + keepParticle + ); + } +} + + +template +void Foam::CloudFilmTransfer::wetSplashInteraction +( + fv::filmCloudTransfer& filmCloudTransfer, + parcelType& p, + const polyPatch& pp, + const label facei, + bool& keepParticle +) +{ + if (debug) + { + Info<< "Parcel " << p.origId() << " wetSplashInteraction" << endl; + } + + const fluidThermo& carrierThermo = + static_cast&>(this->owner()) + .carrierThermo(); + + const parcelThermo& thermo = + static_cast&>(this->owner()).thermo(); + + const liquidProperties& liq = thermo.liquids().properties()[0]; + + // Patch face velocity and normal + const vector& Up = this->owner().U().boundaryField()[pp.index()][facei]; + const vector& nf = pp.faceNormals()[facei]; + + // Local pressure + const scalar pc = carrierThermo.p()[p.cell()]; + + // Retrieve parcel properties + const scalar m = p.mass()*p.nParticle(); + const scalar rho = p.rho(); + const scalar d = p.d(); + vector& U = p.U(); + const scalar sigma = liq.sigma(pc, p.T()); + const scalar mu = liq.mu(pc, p.T()); + const vector Urel = p.U() - Up; + const vector Un = nf*(Urel & nf); + const vector Ut = Urel - Un; + + // Laplace number + const scalar La = rho*sigma*d/sqr(mu); + + // Weber number + const scalar We = rho*magSqr(Un)*d/sigma; + + // Critical Weber number + const scalar Wec = Awet_*pow(La, -0.183); + + if (We < 2) // Adhesion - assume absorb + { + absorbInteraction(filmCloudTransfer, p, pp, facei, m, keepParticle); + } + else if ((We >= 2) && (We < 20)) // Bounce + { + // Incident angle of impingement + const scalar theta = pi/2 - acos(U/mag(U) & nf); + + // Restitution coefficient + const scalar epsilon = 0.993 - theta*(1.76 - theta*(1.56 - theta*0.49)); + + // Update parcel velocity + U = -epsilon*(Un) + 5.0/7.0*(Ut); + + keepParticle = true; + return; + } + else if ((We >= 20) && (We < Wec)) // Spread - assume absorb + { + absorbInteraction(filmCloudTransfer, p, pp, facei, m, keepParticle); + } + else // Splash + { + // Ratio of incident mass to splashing mass + // splash mass can be > incident mass due to film entrainment + const scalar mRatio = 0.2 + 0.9*rndGen_.sample01(); + splashInteraction + ( + filmCloudTransfer, + p, + pp, + facei, + mRatio, + We, + Wec, + sigma, + keepParticle + ); + } +} + + +template +void Foam::CloudFilmTransfer::splashInteraction +( + fv::filmCloudTransfer& filmCloudTransfer, + const parcelType& p, + const polyPatch& pp, + const label facei, + const scalar mRatio, + const scalar We, + const scalar Wec, + const scalar sigma, + bool& keepParticle +) +{ + // Patch face velocity and normal + const fvMesh& mesh = this->owner().mesh(); + const vector& Up = this->owner().U().boundaryField()[pp.index()][facei]; + const vector& nf = pp.faceNormals()[facei]; + + // Determine direction vectors tangential to patch normal + const vector tanVec1 = normalised(perpendicular(nf)); + const vector tanVec2 = nf^tanVec1; + + // Retrieve parcel properties + const scalar np = p.nParticle(); + const scalar m = p.mass()*np; + const scalar d = p.d(); + const vector Urel = p.U() - Up; + const vector Un = nf*(Urel & nf); + const vector Ut = Urel - Un; + const vector& posC = mesh.C()[p.cell()]; + const vector& posCf = mesh.Cf().boundaryField()[pp.index()][facei]; + + // Total mass of (all) splashed parcels + const scalar mSplash = m*mRatio; + + // Number of splashed particles per incoming particle + const scalar Ns = 5.0*(We/Wec - 1.0); + + // Average diameter of splashed particles + const scalar dBarSplash = 1/cbrt(6.0)*cbrt(mRatio/Ns)*d + rootVSmall; + + // Cumulative diameter splash distribution + const scalar dMax = 0.9*cbrt(mRatio)*d; + const scalar dMin = 0.1*dMax; + const scalar K = exp(-dMin/dBarSplash) - exp(-dMax/dBarSplash); + + // Surface energy of secondary parcels [J] + scalar ESigmaSec = 0; + + // Sample splash distribution to determine secondary parcel diameters + scalarList dNew(parcelsPerSplash_); + scalarList npNew(parcelsPerSplash_); + forAll(dNew, i) + { + const scalar y = rndGen_.sample01(); + dNew[i] = -dBarSplash*log(exp(-dMin/dBarSplash) - y*K); + npNew[i] = mRatio*np*pow3(d)/pow3(dNew[i])/parcelsPerSplash_; + ESigmaSec += npNew[i]*sigma*p.areaS(dNew[i]); + } + + // Incident kinetic energy [J] + const scalar EKIn = 0.5*m*magSqr(Un); + + // Incident surface energy [J] + const scalar ESigmaIn = np*sigma*p.areaS(d); + + // Dissipative energy + const scalar Ed = max(0.8*EKIn, np*Wec/12*pi*sigma*sqr(d)); + + // Total energy [J] + const scalar EKs = EKIn + ESigmaIn - ESigmaSec - Ed; + + // Switch to absorb if insufficient energy for splash + if (EKs <= 0) + { + absorbInteraction(filmCloudTransfer, p, pp, facei, m, keepParticle); + return; + } + + // Helper variables to calculate magUns0 + const scalar logD = log(d); + const scalar coeff2 = log(dNew[0]) - logD + rootVSmall; + scalar coeff1 = 0.0; + forAll(dNew, i) + { + coeff1 += sqr(log(dNew[i]) - logD); + } + + // Magnitude of the normal velocity of the first splashed parcel + const scalar magUns0 = + sqrt(2.0*parcelsPerSplash_*EKs/mSplash/(1.0 + coeff1/sqr(coeff2))); + + // Set splashed parcel properties + forAll(dNew, i) + { + const vector dirVec = splashDirection(tanVec1, tanVec2, -nf); + + // Create a new parcel by copying source parcel + parcelType* pPtr = new parcelType(p); + + pPtr->origId() = pPtr->getNewParticleID(); + + pPtr->origProc() = Pstream::myProcNo(); + + if (splashParcelType_ >= 0) + { + pPtr->typeId() = splashParcelType_; + } + + // Perturb new parcels towards the owner cell centre + pPtr->track(mesh, 0.5*rndGen_.sample01()*(posC - posCf), 0); + + pPtr->nParticle() = npNew[i]; + + pPtr->d() = dNew[i]; + + pPtr->U() = dirVec*(mag(Cf_*Ut) + magUns0*(log(dNew[i]) - logD)/coeff2); + + // Apply correction to velocity for 2-D cases + meshTools::constrainDirection(mesh, mesh.solutionD(), pPtr->U()); + + // Add the new parcel + this->owner().addParticle(pPtr); + + nParcelsSplashed_++; + } + + // Transfer remaining part of parcel to film 0 - splashMass can be -ve + // if entraining from the film + const scalar mDash = m - mSplash; + absorbInteraction(filmCloudTransfer, p, pp, facei, mDash, keepParticle); +} + + +template +Foam::UPtrList& +Foam::CloudFilmTransfer::filmTransferPtrs() const +{ + if (filmTransfers_.empty()) + { + const polyBoundaryMesh& patches = this->owner().mesh().boundaryMesh(); + + label nFilms = 0; + + forAll(patches, patchi) + { + if (isA(patches[patchi])) + { + const mappedPatchBase& mpb = + refCast(patches[patchi]); + + Foam::fvModels& fvModels + ( + fvModels::New + ( + refCast(mpb.nbrMesh()) + ) + ); + + fv::filmCloudTransfer* filmCloudPtr = nullptr; + + forAll(fvModels, i) + { + if (isType(fvModels[i])) + { + filmCloudPtr = + &refCast(fvModels[i]); + } + } + + if (filmCloudPtr) + { + Info<< "Found filmCloudTransfer fvModel " + "for the film region " << mpb.nbrMesh().name() + << endl; + + filmTransfers_.resize(nFilms + 1); + filmPatches_.resize(nFilms + 1); + + filmTransfers_.set(nFilms, filmCloudPtr); + filmPatches_[nFilms] = patchi; + + nFilms++; + } + } + } + } + + return filmTransfers_; +} + + +template +const Foam::labelList& Foam::CloudFilmTransfer::filmPatches() const +{ + // Ensure filmPatches_ has been initialise + filmTransferPtrs(); + + return filmPatches_; +} + + +template +void Foam::CloudFilmTransfer::cacheFilmFields(const label filmi) +{ + // Film->Cloud transfer Not implemented yet + + fv::filmCloudTransfer& filmCloudTransfer = filmTransferPtrs()[filmi]; + + filmCloudTransfer.resetFromCloudFields(); + + this->massParcelPatch_.setSize + ( + this->owner().mesh().boundary()[filmPatches_[filmi]].size(), + 0.0 + ); +} + + +template +void Foam::CloudFilmTransfer::setParcelProperties +( + parcelType& p, + const label filmFacei +) const +{ + // Set parcel properties + const scalar vol = + mathematical::pi/6.0*pow3(this->diameterParcelPatch_[filmFacei]); + p.d() = this->diameterParcelPatch_[filmFacei]; + p.U() = UFilmPatch_[filmFacei]; + p.rho() = rhoFilmPatch_[filmFacei]; + + p.nParticle() = this->massParcelPatch_[filmFacei]/p.rho()/vol; + + if (this->ejectedParcelType_ >= 0) + { + p.typeId() = this->ejectedParcelType_; + } + + // Set parcel properties + p.T() = TFilmPatch_[filmFacei]; + p.Cp() = CpFilmPatch_[filmFacei]; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::CloudFilmTransfer::CloudFilmTransfer +( + const dictionary& dict, + CloudType& owner +) +: + SurfaceFilmModel(dict, owner, typeName), + rndGen_(owner.rndGen()), + UFilmPatch_(0), + rhoFilmPatch_(0), + TFilmPatch_(0), + CpFilmPatch_(0), + interactionType_ + ( + interactionTypeNames_.read(this->coeffDict().lookup("interactionType")) + ), + deltaWet_(0.0), + splashParcelType_(0), + parcelsPerSplash_(0), + Adry_(0.0), + Awet_(0.0), + Cf_(0.0), + nParcelsSplashed_(0) +{ + Info<< " Applying " << interactionTypeNames_[interactionType_] + << " interaction model" << endl; + + if (interactionType_ == interactionType::splashBai) + { + this->coeffDict().lookup("deltaWet") >> deltaWet_; + splashParcelType_ = + this->coeffDict().lookupOrDefault("splashParcelType", -1); + parcelsPerSplash_ = + this->coeffDict().lookupOrDefault("parcelsPerSplash", 2); + this->coeffDict().lookup("Adry") >> Adry_; + this->coeffDict().lookup("Awet") >> Awet_; + this->coeffDict().lookup("Cf") >> Cf_; + } +} + + +template +Foam::CloudFilmTransfer::CloudFilmTransfer +( + const CloudFilmTransfer& sfm +) +: + SurfaceFilmModel(sfm), + rndGen_(sfm.rndGen_), + UFilmPatch_(sfm.UFilmPatch_), + rhoFilmPatch_(sfm.rhoFilmPatch_), + TFilmPatch_(sfm.TFilmPatch_), + CpFilmPatch_(sfm.CpFilmPatch_), + interactionType_(sfm.interactionType_), + deltaWet_(sfm.deltaWet_), + splashParcelType_(sfm.splashParcelType_), + parcelsPerSplash_(sfm.parcelsPerSplash_), + Adry_(sfm.Adry_), + Awet_(sfm.Awet_), + Cf_(sfm.Cf_), + nParcelsSplashed_(sfm.nParcelsSplashed_) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::CloudFilmTransfer::~CloudFilmTransfer() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +bool Foam::CloudFilmTransfer::transferParcel +( + parcelType& p, + const polyPatch& pp, + bool& keepParticle +) +{ + const label patchi = pp.index(); + + forAll(this->filmTransferPtrs(), filmi) + { + if (patchi == filmPatches_[filmi]) + { + fv::filmCloudTransfer& filmCloudTransfer = + filmTransferPtrs()[filmi]; + + const label facei = pp.whichFace(p.face()); + + switch (interactionType_) + { + case interactionType::bounce: + { + bounceInteraction(p, pp, facei, keepParticle); + break; + } + case interactionType::absorb: + { + absorbInteraction + ( + filmCloudTransfer, + p, + pp, + facei, + p.nParticle()*p.mass(), + keepParticle + ); + break; + } + case interactionType::splashBai: + { + if (this->deltaFilmPatch_[facei] < deltaWet_) + { + drySplashInteraction + ( + filmCloudTransfer, + p, + pp, + facei, + keepParticle + ); + } + else + { + wetSplashInteraction + ( + filmCloudTransfer, + p, + pp, + facei, + keepParticle + ); + } + break; + } + default: + { + FatalErrorInFunction + << "Unknown interaction type enumeration" + << abort(FatalError); + } + } + + // Transfer parcel/parcel interactions complete + return true; + } + } + + // Parcel not interacting with film + return false; +} + + +template +void Foam::CloudFilmTransfer::info(Ostream& os) +{ + SurfaceFilmModel::info(os); + + label nSplash0 = this->template getModelProperty