Merge branch 'film'

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
 

src/lagrangian/basic/Particle/ParticleI.H

@@ -1117,6 +1117,13 @@ inline Foam::label Foam::Particle<ParticleType>::origProc() const
 }
 
 
+template<class ParticleType>
+inline Foam::label& Foam::Particle<ParticleType>::origProc()
+{
+    return origProc_;
+}
+
+
 template<class ParticleType>
 inline Foam::label Foam::Particle<ParticleType>::origId() const
 {
@@ -1124,6 +1131,13 @@ inline Foam::label Foam::Particle<ParticleType>::origId() const
 }
 
 
+template<class ParticleType>
+inline Foam::label& Foam::Particle<ParticleType>::origId()
+{
+    return origId_;
+}
+
+
 template<class ParticleType>
 inline bool Foam::Particle<ParticleType>::softImpact() const
 {

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -426,11 +426,8 @@ bool Foam::KinematicParcel<ParcelType>::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;
     }

src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C

@@ -37,10 +37,12 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel(CloudType& owner)
 :
     SubModelBase<CloudType>(owner),
     g_(dimensionedVector("zero", dimAcceleration, vector::zero)),
+    ejectedParcelType_(0),
     massParcelPatch_(0),
     diameterParcelPatch_(0),
     UFilmPatch_(0),
     rhoFilmPatch_(0),
+    deltaFilmPatch_(0),
     nParcelsTransferred_(0),
     nParcelsInjected_(0)
 {}
@@ -57,10 +59,15 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel
 :
     SubModelBase<CloudType>(owner, dict, type),
     g_(g),
+    ejectedParcelType_
+    (
+        this->coeffDict().lookupOrDefault("ejectedParcelType", -1)
+    ),
     massParcelPatch_(0),
     diameterParcelPatch_(0),
     UFilmPatch_(0),
     rhoFilmPatch_(0),
+    deltaFilmPatch_(owner.mesh().boundary().size()),
     nParcelsTransferred_(0),
     nParcelsInjected_(0)
 {}
@@ -74,10 +81,12 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel
 :
     SubModelBase<CloudType>(sfm),
     g_(sfm.g_),
+    ejectedParcelType_(sfm.ejectedParcelType_),
     massParcelPatch_(sfm.massParcelPatch_),
     diameterParcelPatch_(sfm.diameterParcelPatch_),
     UFilmPatch_(sfm.UFilmPatch_),
     rhoFilmPatch_(sfm.rhoFilmPatch_),
+    deltaFilmPatch_(sfm.deltaFilmPatch_),
     nParcelsTransferred_(sfm.nParcelsTransferred_),
     nParcelsInjected_(sfm.nParcelsInjected_)
 {}
@@ -95,16 +104,18 @@ Foam::SurfaceFilmModel<CloudType>::~SurfaceFilmModel()
 template<class CloudType>
 bool Foam::SurfaceFilmModel<CloudType>::transferParcel
 (
-    const parcelType& p,
-    const label patchI
+    parcelType& p,
+    const polyPatch& pp,
+    bool& keepParticle
 )
 {
     notImplemented
     (
         "bool Foam::SurfaceFilmModel<CloudType>::transferParcel"
         "("
-            "const parcelType&, "
-            "const label"
+            "parcelType&, "
+            "const label, "
+            "const bool&"
         ")"
     );
 
@@ -145,7 +156,7 @@ void Foam::SurfaceFilmModel<CloudType>::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 +207,7 @@ template<class CloudType>
 void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
 (
     const label filmPatchI,
+    const label primaryPatchI,
     const mapDistribute& distMap,
     const surfaceFilmModels::surfaceFilmModel& filmModel
 )
@@ -212,6 +224,10 @@ void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
 
     rhoFilmPatch_ = filmModel.rho().boundaryField()[filmPatchI];
     distMap.distribute(rhoFilmPatch_);
+
+    deltaFilmPatch_[primaryPatchI] =
+        filmModel.delta().boundaryField()[filmPatchI];
+    distMap.distribute(deltaFilmPatch_[primaryPatchI]);
 }
 
 
@@ -229,6 +245,11 @@ void Foam::SurfaceFilmModel<CloudType>::setParcelProperties
     p.rho() = rhoFilmPatch_[filmFaceI];
 
     p.nParticle() = massParcelPatch_[filmFaceI]/p.rho()/vol;
+
+    if (ejectedParcelType_ >= 0)
+    {
+        p.typeId() = ejectedParcelType_;
+    }
 }
 
 

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
 
@@ -88,6 +93,9 @@ protected:
             //- Film density / patch face
             scalarList rhoFilmPatch_;
 
+            //- Film height of all film patches / patch face
+            scalarListList deltaFilmPatch_;
+
 
         // Counters
 
@@ -104,6 +112,7 @@ protected:
         virtual void cacheFilmFields
         (
             const label filmPatchI,
+            const label primaryPatchI,
             const mapDistribute& distMap,
             const surfaceFilmModels::surfaceFilmModel& filmModel
         );
@@ -206,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

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<class CloudType>
+Foam::wordList Foam::ThermoSurfaceFilm<CloudType>::interactionTypeNames_
+(
+    IStringStream
+    (
+        "(absorb bounce splashBai)"
+    )()
+);
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class CloudType>
+typename Foam::ThermoSurfaceFilm<CloudType>::interactionType
+Foam::ThermoSurfaceFilm<CloudType>::interactionTypeEnum(const word& it) const
+{
+    forAll(interactionTypeNames_, i)
+    {
+        if (interactionTypeNames_[i] == it)
+        {
+            return interactionType(i);
+        }
+    }
+
+    FatalErrorIn
+    (
+        "ThermoSurfaceFilm<CloudType>::interactionType "
+        "ThermoSurfaceFilm<CloudType>::interactionTypeEnum"
+        "("
+            "const word& it"
+        ") const"
+    )   << "Unknown interaction type " << it
+        << ". Valid interaction types include: " << interactionTypeNames_
+        << abort(FatalError);
+
+    return interactionType(0);
+}
+
+
+template<class CloudType>
+Foam::word Foam::ThermoSurfaceFilm<CloudType>::interactionTypeStr
+(
+    const interactionType& it
+) const
+{
+    if (it >= interactionTypeNames_.size())
+    {
+        FatalErrorIn
+        (
+            "ThermoSurfaceFilm<CloudType>::interactionType "
+            "ThermoSurfaceFilm<CloudType>::interactionTypeStr"
+            "("
+                "const interactionType& it"
+            ") const"
+        )   << "Unknown interaction type enumeration" << abort(FatalError);
+    }
+
+    return interactionTypeNames_[it];
+}
+
+
+template<class CloudType>
+Foam::vector Foam::ThermoSurfaceFilm<CloudType>::tangentVector
+(
+    const vector& v
+) const
+{
+    vector tangent = vector::zero;
+    scalar magTangent = 0.0;
+
+    while (magTangent < SMALL)
+    {
+        vector vTest = rndGen_.sample01<vector>();
+        tangent = vTest - (vTest & v)*v;
+        magTangent = mag(tangent);
+    }
+
+    return tangent/magTangent;
+}
+
+
+template<class CloudType>
+Foam::vector Foam::ThermoSurfaceFilm<CloudType>::splashDirection
+(
+    const vector& tanVec1,
+    const vector& tanVec2,
+    const vector& nf
+) const
+{
+    // azimuthal angle [rad]
+    const scalar phiSi = twoPi*rndGen_.sample01<scalar>();
+
+    // ejection angle [rad]
+    const scalar thetaSi = pi/180.0*(rndGen_.sample01<scalar>()*(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<class CloudType>
+void Foam::ThermoSurfaceFilm<CloudType>::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<class CloudType>
+void Foam::ThermoSurfaceFilm<CloudType>::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<class CloudType>
+void Foam::ThermoSurfaceFilm<CloudType>::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<scalar>();
+        splashInteraction
+            (filmModel, p, pp, faceI, mRatio, We, Wec, sigma, keepParticle);
+    }
+}
+
+
+template<class CloudType>
+void Foam::ThermoSurfaceFilm<CloudType>::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<scalar>();
+        splashInteraction
+            (filmModel, p, pp, faceI, mRatio, We, Wec, sigma, keepParticle);
+    }
+}
+
+
+template<class CloudType>
+void Foam::ThermoSurfaceFilm<CloudType>::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