From 9046d1110a0332b2b66b8f92789338cf8ddf6850 Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 11 Nov 2010 12:37:52 +0000 Subject: [PATCH 1/4] ENH: Added access to film height at top level --- .../surfaceFilmModel/noFilm/noFilm.C | 14 +++++++++++--- .../surfaceFilmModel/noFilm/noFilm.H | 3 +++ .../surfaceFilmModel/surfaceFilmModel.H | 3 +++ 3 files changed, 17 insertions(+), 3 deletions(-) diff --git a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C index 619226bb8f..63ea47880d 100644 --- a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C +++ b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C @@ -136,9 +136,17 @@ void Foam::surfaceFilmModels::noFilm::addSources } +const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::delta() const +{ + FatalErrorIn("const volScalarField& noFilm::delta() const") + << "delta field not available for " << type() << abort(FatalError); + + return volScalarField::null(); +} + const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::U() const { - FatalErrorIn("const volScalarField& noFilm::U() const") + FatalErrorIn("const volVectorField& noFilm::U() const") << "U field not available for " << type() << abort(FatalError); return volVectorField::null(); @@ -147,7 +155,7 @@ const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::U() const const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::Us() const { - FatalErrorIn("const volScalarField& noFilm::Us() const") + FatalErrorIn("const volVectorField& noFilm::Us() const") << "Us field not available for " << type() << abort(FatalError); return volVectorField::null(); @@ -156,7 +164,7 @@ const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::Us() const const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::Uw() const { - FatalErrorIn("const volScalarField& noFilm::Uw() const") + FatalErrorIn("const volVectorField& noFilm::Uw() const") << "Uw field not available for " << type() << abort(FatalError); return volVectorField::null(); diff --git a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.H b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.H index 7c3d7f89d9..40726066ac 100644 --- a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.H +++ b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.H @@ -126,6 +126,9 @@ public: // Fields + //- Return the film thickness [m] + virtual const volScalarField& delta() const; + //- Return the film velocity [m/s] virtual const volVectorField& U() const; diff --git a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.H b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.H index 511a6204e4..079f538c83 100644 --- a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.H +++ b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.H @@ -220,6 +220,9 @@ public: // Fields + //- Return the film thickness [m] + virtual const volScalarField& delta() const = 0; + //- Return the film velocity [m/s] virtual const volVectorField& U() const = 0; From 8b20d01a23f6d0310e87107346c34f696f4b7ec8 Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 11 Nov 2010 12:50:51 +0000 Subject: [PATCH 2/4] ENH: Cache film height field in cloud film models --- .../SurfaceFilmModel/SurfaceFilmModel.C | 10 +++++++++- .../SurfaceFilmModel/SurfaceFilmModel.H | 4 ++++ .../ThermoSurfaceFilm/ThermoSurfaceFilm.C | 2 ++ .../ThermoSurfaceFilm/ThermoSurfaceFilm.H | 1 + 4 files changed, 16 insertions(+), 1 deletion(-) diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C index d71ab36172..8b9294e57a 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C @@ -41,6 +41,7 @@ Foam::SurfaceFilmModel::SurfaceFilmModel(CloudType& owner) diameterParcelPatch_(0), UFilmPatch_(0), rhoFilmPatch_(0), + deltaFilmPatch_(0), nParcelsTransferred_(0), nParcelsInjected_(0) {} @@ -61,6 +62,7 @@ Foam::SurfaceFilmModel::SurfaceFilmModel diameterParcelPatch_(0), UFilmPatch_(0), rhoFilmPatch_(0), + deltaFilmPatch_(owner.mesh().boundary().size()), nParcelsTransferred_(0), nParcelsInjected_(0) {} @@ -78,6 +80,7 @@ Foam::SurfaceFilmModel::SurfaceFilmModel diameterParcelPatch_(sfm.diameterParcelPatch_), UFilmPatch_(sfm.UFilmPatch_), rhoFilmPatch_(sfm.rhoFilmPatch_), + deltaFilmPatch_(sfm.deltaFilmPatch_), nParcelsTransferred_(sfm.nParcelsTransferred_), nParcelsInjected_(sfm.nParcelsInjected_) {} @@ -145,7 +148,7 @@ void Foam::SurfaceFilmModel::inject(TrackData& td) const label filmPatchI = filmPatches[i]; const mapDistribute& distMap = wpp.map(); - cacheFilmFields(filmPatchI, distMap, filmModel); + cacheFilmFields(filmPatchI, primaryPatchI, distMap, filmModel); forAll(injectorCellsPatch, j) { @@ -196,6 +199,7 @@ template void Foam::SurfaceFilmModel::cacheFilmFields ( const label filmPatchI, + const label primaryPatchI, const mapDistribute& distMap, const surfaceFilmModels::surfaceFilmModel& filmModel ) @@ -212,6 +216,10 @@ void Foam::SurfaceFilmModel::cacheFilmFields rhoFilmPatch_ = filmModel.rho().boundaryField()[filmPatchI]; distMap.distribute(rhoFilmPatch_); + + deltaFilmPatch_[primaryPatchI] = + filmModel.delta().boundaryField()[filmPatchI]; + distMap.distribute(deltaFilmPatch_[primaryPatchI]); } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H index 25452ce1d7..54dad5abcf 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H @@ -88,6 +88,9 @@ protected: //- Film density / patch face scalarList rhoFilmPatch_; + //- Film height of all film patches / patch face + scalarListList deltaFilmPatch_; + // Counters @@ -104,6 +107,7 @@ protected: virtual void cacheFilmFields ( const label filmPatchI, + const label primaryPatchI, const mapDistribute& distMap, const surfaceFilmModels::surfaceFilmModel& filmModel ); diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C index efc60d46b6..2b28213c20 100644 --- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C +++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C @@ -131,6 +131,7 @@ template void Foam::ThermoSurfaceFilm::cacheFilmFields ( const label filmPatchI, + const label primaryPatchI, const mapDistribute& distMap, const surfaceFilmModels::surfaceFilmModel& filmModel ) @@ -138,6 +139,7 @@ void Foam::ThermoSurfaceFilm::cacheFilmFields SurfaceFilmModel::cacheFilmFields ( filmPatchI, + primaryPatchI, distMap, filmModel ); diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H index 74944fc58c..9f7e02e5e7 100644 --- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H +++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H @@ -75,6 +75,7 @@ protected: virtual void cacheFilmFields ( const label filmPatchI, + const label primaryPatchI, const mapDistribute& distMap, const surfaceFilmModels::surfaceFilmModel& filmModel ); From 0647c7df680fee273aa20880d12ae5947d96bdb2 Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 11 Nov 2010 13:29:06 +0000 Subject: [PATCH 3/4] ENH: Added write access to particle origId and origProc --- src/lagrangian/basic/Particle/Particle.H | 10 ++++++++-- src/lagrangian/basic/Particle/ParticleI.H | 14 ++++++++++++++ 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/src/lagrangian/basic/Particle/Particle.H b/src/lagrangian/basic/Particle/Particle.H index d17a9c995b..8e7604bef5 100644 --- a/src/lagrangian/basic/Particle/Particle.H +++ b/src/lagrangian/basic/Particle/Particle.H @@ -517,12 +517,18 @@ public: //- Return the fraction of time-step completed inline scalar stepFraction() const; - //- Return the originating processor id + //- Return const access to the originating processor id inline label origProc() const; - //- Return the particle id on originating processor + //- Return the originating processor id for manipulation + inline label& origProc(); + + //- Return const access to the particle id on originating processor inline label origId() const; + //- Return the particle id on originating processor for manipulation + inline label& origId(); + // Track diff --git a/src/lagrangian/basic/Particle/ParticleI.H b/src/lagrangian/basic/Particle/ParticleI.H index 50c727fe4a..0ef468ea2f 100644 --- a/src/lagrangian/basic/Particle/ParticleI.H +++ b/src/lagrangian/basic/Particle/ParticleI.H @@ -1117,6 +1117,13 @@ inline Foam::label Foam::Particle::origProc() const } +template +inline Foam::label& Foam::Particle::origProc() +{ + return origProc_; +} + + template inline Foam::label Foam::Particle::origId() const { @@ -1124,6 +1131,13 @@ inline Foam::label Foam::Particle::origId() const } +template +inline Foam::label& Foam::Particle::origId() +{ + return origId_; +} + + template inline bool Foam::Particle::softImpact() const { From 9a85453a5bd6c2d176dc3d01402bd7e1c9a76058 Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 11 Nov 2010 13:29:33 +0000 Subject: [PATCH 4/4] ENH: Added Bai particle/film splash modelling --- .../KinematicParcel/KinematicParcel.C | 5 +- .../SurfaceFilmModel/SurfaceFilmModel.C | 21 +- .../SurfaceFilmModel/SurfaceFilmModel.H | 10 +- .../ThermoSurfaceFilm/ThermoSurfaceFilm.C | 591 ++++++++++++++++-- .../ThermoSurfaceFilm/ThermoSurfaceFilm.H | 182 +++++- 5 files changed, 742 insertions(+), 67 deletions(-) diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index a5f1b44d84..cccb87fea1 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -426,11 +426,8 @@ bool Foam::KinematicParcel::hitPatch td.cloud().postProcessing().postPatch(p, patchI); // Invoke surface film model - if (td.cloud().surfaceFilm().transferParcel(p, patchI)) + if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle)) { - // Parcel transferred to the surface film - td.keepParticle = false; - // All interactions done return true; } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C index 8b9294e57a..f1ef48de9d 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C @@ -37,6 +37,7 @@ Foam::SurfaceFilmModel::SurfaceFilmModel(CloudType& owner) : SubModelBase(owner), g_(dimensionedVector("zero", dimAcceleration, vector::zero)), + ejectedParcelType_(0), massParcelPatch_(0), diameterParcelPatch_(0), UFilmPatch_(0), @@ -58,6 +59,10 @@ Foam::SurfaceFilmModel::SurfaceFilmModel : SubModelBase(owner, dict, type), g_(g), + ejectedParcelType_ + ( + this->coeffDict().lookupOrDefault("ejectedParcelType", -1) + ), massParcelPatch_(0), diameterParcelPatch_(0), UFilmPatch_(0), @@ -76,6 +81,7 @@ Foam::SurfaceFilmModel::SurfaceFilmModel : SubModelBase(sfm), g_(sfm.g_), + ejectedParcelType_(sfm.ejectedParcelType_), massParcelPatch_(sfm.massParcelPatch_), diameterParcelPatch_(sfm.diameterParcelPatch_), UFilmPatch_(sfm.UFilmPatch_), @@ -98,16 +104,18 @@ Foam::SurfaceFilmModel::~SurfaceFilmModel() template bool Foam::SurfaceFilmModel::transferParcel ( - const parcelType& p, - const label patchI + parcelType& p, + const polyPatch& pp, + bool& keepParticle ) { notImplemented ( "bool Foam::SurfaceFilmModel::transferParcel" "(" - "const parcelType&, " - "const label" + "parcelType&, " + "const label, " + "const bool&" ")" ); @@ -237,6 +245,11 @@ void Foam::SurfaceFilmModel::setParcelProperties p.rho() = rhoFilmPatch_[filmFaceI]; p.nParticle() = massParcelPatch_[filmFaceI]/p.rho()/vol; + + if (ejectedParcelType_ >= 0) + { + p.typeId() = ejectedParcelType_; + } } diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H index 54dad5abcf..35c4157512 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H +++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H @@ -73,6 +73,11 @@ protected: //- Gravitational acceleration constant const dimensionedVector& g_; + //- Ejected parcel type label - id assigned to identify parcel for + // post-processing. If not specified, defaults to originating cloud + // type + label ejectedParcelType_; + // Cached injector fields per film patch @@ -210,8 +215,9 @@ public: // Returns true if parcel is to be transferred virtual bool transferParcel ( - const parcelType& p, - const label patchI + parcelType& p, + const polyPatch& pp, + bool& keepParticle ); //- Inject parcels into the cloud diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C index 2b28213c20..1a037fdc77 100644 --- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C +++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C @@ -25,6 +25,460 @@ License #include "ThermoSurfaceFilm.H" #include "addToRunTimeSelectionTable.H" +#include "mathematicalConstants.H" +#include "Pstream.H" + +using namespace Foam::constant::mathematical; + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +template +Foam::wordList Foam::ThermoSurfaceFilm::interactionTypeNames_ +( + IStringStream + ( + "(absorb bounce splashBai)" + )() +); + + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +template +typename Foam::ThermoSurfaceFilm::interactionType +Foam::ThermoSurfaceFilm::interactionTypeEnum(const word& it) const +{ + forAll(interactionTypeNames_, i) + { + if (interactionTypeNames_[i] == it) + { + return interactionType(i); + } + } + + FatalErrorIn + ( + "ThermoSurfaceFilm::interactionType " + "ThermoSurfaceFilm::interactionTypeEnum" + "(" + "const word& it" + ") const" + ) << "Unknown interaction type " << it + << ". Valid interaction types include: " << interactionTypeNames_ + << abort(FatalError); + + return interactionType(0); +} + + +template +Foam::word Foam::ThermoSurfaceFilm::interactionTypeStr +( + const interactionType& it +) const +{ + if (it >= interactionTypeNames_.size()) + { + FatalErrorIn + ( + "ThermoSurfaceFilm::interactionType " + "ThermoSurfaceFilm::interactionTypeStr" + "(" + "const interactionType& it" + ") const" + ) << "Unknown interaction type enumeration" << abort(FatalError); + } + + return interactionTypeNames_[it]; +} + + +template +Foam::vector Foam::ThermoSurfaceFilm::tangentVector +( + const vector& v +) const +{ + vector tangent = vector::zero; + scalar magTangent = 0.0; + + while (magTangent < SMALL) + { + vector vTest = rndGen_.sample01(); + tangent = vTest - (vTest & v)*v; + magTangent = mag(tangent); + } + + return tangent/magTangent; +} + + +template +Foam::vector Foam::ThermoSurfaceFilm::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::ThermoSurfaceFilm::absorbInteraction +( + surfaceFilmModels::surfaceFilmModel& filmModel, + const parcelType& p, + const polyPatch& pp, + const label faceI, + const scalar mass, + bool& keepParticle +) +{ + if (debug) + { + Info<< "Parcel " << p.origId() << " absorbInteraction" << 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; + + // Parcel normal velocity + const vector Un = nf*(Urel & nf); + + // Parcel tangential velocity + const vector Ut = Urel - Un; + + filmModel.addSources + ( + pp.index(), + faceI, + mass, // mass + mass*Ut, // tangential momentum + mass*mag(Un), // impingement pressure + mass*p.hs() // energy + ); + + this->nParcelsTransferred()++; + + keepParticle = false; +} + + +template +void Foam::ThermoSurfaceFilm::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::ThermoSurfaceFilm::drySplashInteraction +( + surfaceFilmModels::surfaceFilmModel& filmModel, + const parcelType& p, + const polyPatch& pp, + const label faceI, + bool& keepParticle +) +{ + if (debug) + { + Info<< "Parcel " << p.origId() << " drySplashInteraction" << endl; + } + + const liquid& 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 = thermo_.thermo().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(filmModel, p, pp, faceI, m, keepParticle); + } + else // splash + { + // ratio of incident mass to splashing mass + const scalar mRatio = 0.2 + 0.6*rndGen_.sample01(); + splashInteraction + (filmModel, p, pp, faceI, mRatio, We, Wec, sigma, keepParticle); + } +} + + +template +void Foam::ThermoSurfaceFilm::wetSplashInteraction +( + surfaceFilmModels::surfaceFilmModel& filmModel, + parcelType& p, + const polyPatch& pp, + const label faceI, + bool& keepParticle +) +{ + if (debug) + { + Info<< "Parcel " << p.origId() << " wetSplashInteraction" << endl; + } + + const liquid& 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 = thermo_.thermo().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 < 1) // adhesion - assume absorb + { + absorbInteraction(filmModel, p, pp, faceI, m, keepParticle); + } + else if ((We >= 1) && (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/7*(Ut); + + keepParticle = true; + } + else if ((We >= 20) && (We < Wec)) // spread - assume absorb + { + absorbInteraction(filmModel, 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 + (filmModel, p, pp, faceI, mRatio, We, Wec, sigma, keepParticle); + } +} + + +template +void Foam::ThermoSurfaceFilm::splashInteraction +( + surfaceFilmModels::surfaceFilmModel& filmModel, + 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 = tangentVector(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 dMin = 0.01*d; + const scalar dMax = d; + const scalar K = exp(-dMin/dBarSplash) - exp(-dMax/dBarSplash); + + // surface energy of secondary parcels [J] + scalar ESigmaSec = 0; + + // sample splash distribution to detrmine secondary parcel diameters + scalarList dNew(parcelsPerSplash_); + forAll(dNew, i) + { + const scalar y = rndGen_.sample01(); + dNew[i] = -dBarSplash*log(exp(-dMin/dBarSplash) - y*K); + ESigmaSec += sigma*p.areaS(dNew[i]); + } + + // incident kinetic energy [J] + const scalar EKIn = 0.5*m*magSqr(Urel); + + // incident surface energy [J] + const scalar ESigmaIn = sigma*p.areaS(d); + + // dissipative energy + const scalar Ed = max(0.8*EKIn, 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(filmModel, 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 + 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() = this->owner().getNewParticleID(); + + pPtr->origProc() = Pstream::myProcNo(); + + if (splashParcelType_ >= 0) + { + pPtr->typeId() = splashParcelType_; + } + + // perturb new parcels towards the owner cell centre + pPtr->position() += 0.5*rndGen_.sample01()*(posC - posCf); + + pPtr->nParticle() = mRatio*np*pow3(d)/pow3(dNew[i])/parcelsPerSplash_; + + 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(filmModel, p, pp, faceI, mDash, keepParticle); +} + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -37,9 +491,37 @@ Foam::ThermoSurfaceFilm::ThermoSurfaceFilm ) : SurfaceFilmModel(dict, owner, g, typeName), + rndGen_(owner.rndGen()), + thermo_(owner.db().objectRegistry::lookupObject("SLGThermo")), TFilmPatch_(0), - CpFilmPatch_(0) -{} + CpFilmPatch_(0), + interactionType_ + ( + interactionTypeEnum(this->coeffDict().lookup("interactionType")) + ), + deltaWet_(0.0), + splashParcelType_(0), + parcelsPerSplash_(0), + Adry_(0.0), + Awet_(0.0), + Cf_(0.0), + nParcelsSplashed_(0) +{ + Info<< " Applying " << interactionTypeStr(interactionType_) + << " interaction model" << endl; + + if (interactionType_ == itSplashBai) + { + 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 @@ -49,8 +531,18 @@ Foam::ThermoSurfaceFilm::ThermoSurfaceFilm ) : SurfaceFilmModel(sfm), + rndGen_(sfm.rndGen_), + thermo_(sfm.thermo_), TFilmPatch_(sfm.TFilmPatch_), - CpFilmPatch_(sfm.CpFilmPatch_) + 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_) {} @@ -66,8 +558,9 @@ Foam::ThermoSurfaceFilm::~ThermoSurfaceFilm() template bool Foam::ThermoSurfaceFilm::transferParcel ( - const parcelType& p, - const label patchI + parcelType& p, + const polyPatch& pp, + bool& keepParticle ) { // Retrieve the film model from the owner database @@ -81,49 +574,61 @@ bool Foam::ThermoSurfaceFilm::transferParcel ) ); + const label patchI = pp.index(); + if (filmModel.isFilmPatch(patchI)) { - const polyPatch& pp = this->owner().mesh().boundaryMesh()[patchI]; - const label faceI = pp.whichFace(p.face()); - // Patch face normal - const vector& nf = pp.faceNormals()[faceI]; - - // Relative parcel velocity - const vector Urel = - p.U() - this->owner().U().boundaryField()[patchI][faceI]; - - // Parcel mass - const scalar m = p.nParticle()*p.mass(); - - // Add the particle properties as sources to the film model - filmModel.addSources - ( - patchI, - faceI, - m, // mass - m*(Urel - nf*(Urel & nf)), // tangential momentum - m*mag(Urel & nf), // impingement pressure - m*p.hs() // energy - ); - - if (debug) + switch (interactionType_) { - Info<< "ThermoSurfaceFilm::transferParcel:" << nl - << " Effective increase in film height = " - << p.nParticle()*p.volume()/mag(pp.faceAreas()[faceI]) << endl; + case itBounce: + { + bounceInteraction(p, pp, faceI, keepParticle); + + break; + } + case itAbsorb: + { + const scalar m = p.nParticle()*p.mass(); + absorbInteraction(filmModel, p, pp, faceI, m, keepParticle); + + break; + } + case itSplashBai: + { + bool dry = this->deltaFilmPatch_[patchI][faceI] < deltaWet_; + + if (dry) + { + drySplashInteraction(filmModel, p, pp, faceI, keepParticle); + } + else + { + wetSplashInteraction(filmModel, p, pp, faceI, keepParticle); + } + + break; + } + default: + { + FatalErrorIn + ( + "bool ThermoSurfaceFilm::transferParcel" + "(" + "const parcelType&, " + "const label" + ")" + ) << "Unknown interaction type enumeration" + << abort(FatalError); + } } - this->nParcelsTransferred()++; - - // Flag to remove parcel p from owner cloud return true; } - else - { - return false; - } + + // do not transfer parcel + return false; } @@ -170,10 +675,12 @@ void Foam::ThermoSurfaceFilm::setParcelProperties template void Foam::ThermoSurfaceFilm::info(Ostream& os) const { - os << " Parcels transferred to film = " + os << " Parcels absorbed into film = " << returnReduce(this->nParcelsTransferred(), sumOp