ENH: Adding MPPIC sub-models to Kinematic cloud

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -38,6 +38,12 @@ License
 #include "SurfaceFilmModel.H"
 #include "profiling.H"
 
+#include "PackingModel.H"
+#include "ParticleStressModel.H"
+#include "DampingModel.H"
+#include "IsotropyModel.H"
+#include "TimeScaleModel.H"
+
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 template<class CloudType>
@@ -79,6 +85,33 @@ void Foam::KinematicCloud<CloudType>::setModels()
         ).ptr()
     );
 
+    packingModel_.reset
+    (
+        PackingModel<KinematicCloud<CloudType>>::New
+        (
+            subModelProperties_,
+            *this
+        ).ptr()
+    );
+
+    dampingModel_.reset
+    (
+        DampingModel<KinematicCloud<CloudType>>::New
+        (
+            subModelProperties_,
+            *this
+        ).ptr()
+    );
+
+    isotropyModel_.reset
+    (
+        IsotropyModel<KinematicCloud<CloudType>>::New
+        (
+            subModelProperties_,
+            *this
+        ).ptr()
+    );
+
     UIntegrator_.reset
     (
         integrationScheme::New
@@ -210,7 +243,6 @@ void Foam::KinematicCloud<CloudType>::evolveCloud
 
         injectors_.inject(cloud, td);
 
-
         // Assume that motion will update the cellOccupancy as necessary
         // before it is required.
         cloud.motion(cloud, td);
@@ -264,6 +296,15 @@ void Foam::KinematicCloud<CloudType>::postEvolve
             true
         );
     }
+
+    if (this->dampingModel().active())
+    {
+        this->dampingModel().cacheFields(false);
+    }
+    if (this->packingModel().active())
+    {
+        this->packingModel().cacheFields(false);
+    }
 }
 
 
@@ -285,6 +326,10 @@ void Foam::KinematicCloud<CloudType>::cloudReset(KinematicCloud<CloudType>& c)
     stochasticCollisionModel_.reset(c.stochasticCollisionModel_.ptr());
     surfaceFilmModel_.reset(c.surfaceFilmModel_.ptr());
 
+    packingModel_.reset(c.packingModel_.ptr());
+    dampingModel_.reset(c.dampingModel_.ptr());
+    isotropyModel_.reset(c.isotropyModel_.ptr());
+
     UIntegrator_.reset(c.UIntegrator_.ptr());
 }
 
@@ -375,6 +420,11 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     patchInteractionModel_(nullptr),
     stochasticCollisionModel_(nullptr),
     surfaceFilmModel_(nullptr),
+
+    packingModel_(nullptr),
+    dampingModel_(nullptr),
+    isotropyModel_(nullptr),
+
     UIntegrator_(nullptr),
     UTrans_
     (
@@ -458,6 +508,11 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     patchInteractionModel_(c.patchInteractionModel_->clone()),
     stochasticCollisionModel_(c.stochasticCollisionModel_->clone()),
     surfaceFilmModel_(c.surfaceFilmModel_->clone()),
+
+    packingModel_(c.packingModel_->clone()),
+    dampingModel_(c.dampingModel_->clone()),
+    isotropyModel_(c.isotropyModel_->clone()),
+
     UIntegrator_(c.UIntegrator_->clone()),
     UTrans_
     (
@@ -549,6 +604,11 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     patchInteractionModel_(nullptr),
     stochasticCollisionModel_(nullptr),
     surfaceFilmModel_(nullptr),
+
+    packingModel_(nullptr),
+    dampingModel_(nullptr),
+    isotropyModel_(nullptr),
+
     UIntegrator_(nullptr),
     UTrans_(nullptr),
     UCoeff_(nullptr)
@@ -683,15 +743,32 @@ void Foam::KinematicCloud<CloudType>::preEvolve
     // with topology change due to lazy evaluation of valid mesh dimensions
     label nGeometricD = mesh_.nGeometricD();
 
-    Info<< "\nSolving " << nGeometricD << "-D cloud " << this->name() << endl;
+    Info<< "\nSolving" << nGeometricD << "-D cloud " << this->name() << endl;
 
     this->dispersion().cacheFields(true);
     forces_.cacheFields(true);
-    updateCellOccupancy();
 
     pAmbient_ = constProps_.dict().template
         getOrDefault<scalar>("pAmbient", pAmbient_);
 
+    if (this->dampingModel().active() || this->packingModel().active())
+    {
+         const_cast<typename parcelType::trackingData&>(td).updateAverages(*this);
+    }
+
+    if (this->dampingModel().active())
+    {
+        DebugVar("dampingModel")
+        this->dampingModel().cacheFields(true);
+    }
+    if (this->packingModel().active())
+    {
+        DebugVar("packingModel")
+        this->packingModel().cacheFields(true);
+    }
+
+    updateCellOccupancy();
+
     functions_.preEvolve(td);
 }
 
@@ -702,7 +779,6 @@ void Foam::KinematicCloud<CloudType>::evolve()
     if (solution_.canEvolve())
     {
         typename parcelType::trackingData td(*this);
-
         solve(*this, td);
     }
 }
@@ -719,6 +795,12 @@ void Foam::KinematicCloud<CloudType>::motion
     td.part() = parcelType::trackingData::tpLinearTrack;
     CloudType::move(cloud, td, solution_.trackTime());
 
+    if (isotropyModel_->active())
+    {
+        td.updateAverages(cloud);
+        isotropyModel_->calculate();
+    }
+
     updateCellOccupancy();
 }
 
@@ -749,15 +831,15 @@ void Foam::KinematicCloud<CloudType>::patchData
         // just inside the domain rather than that of the wall itself.
         if (U_.boundaryField()[patchi].fixesValue())
         {
-            const vector Uw1 = U_.boundaryField()[patchi][patchFacei];
+            const vector Uw1(U_.boundaryField()[patchi][patchFacei]);
             const vector& Uw0 =
                 U_.oldTime().boundaryField()[patchi][patchFacei];
 
             const scalar f = p.currentTimeFraction();
 
-            const vector Uw = Uw0 + f*(Uw1 - Uw0);
+            const vector Uw(Uw0 + f*(Uw1 - Uw0));
 
-            const tensor nnw = nw*nw;
+            const tensor nnw(nw*nw);
 
             Up = (nnw & Up) + Uw - (nnw & Uw);
         }

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -95,6 +95,15 @@ class SurfaceFilmModel;
 template<class CloudType>
 class StochasticCollisionModel;
 
+template<class CloudType>
+class PackingModel;
+
+template<class CloudType>
+class DampingModel;
+
+template<class CloudType>
+class IsotropyModel;
+
 
 /*---------------------------------------------------------------------------*\
                        Class KinematicCloud Declaration
@@ -129,7 +138,7 @@ public:
 
 private:
 
-    // Private data
+    // Private Data
 
         //- Cloud copy pointer
         autoPtr<KinematicCloud<CloudType>> cloudCopyPtr_;
@@ -146,7 +155,7 @@ private:
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- References to the mesh and time databases
         const fvMesh& mesh_;
@@ -225,6 +234,18 @@ protected:
             autoPtr<SurfaceFilmModel<KinematicCloud<CloudType>>>
                 surfaceFilmModel_;
 
+            //- Packing model
+            autoPtr<PackingModel<KinematicCloud<CloudType>>>
+                packingModel_;
+
+            //- Damping model
+            autoPtr<DampingModel<KinematicCloud<CloudType>>>
+                dampingModel_;
+
+            //- Exchange model
+            autoPtr<IsotropyModel<KinematicCloud<CloudType>>>
+                isotropyModel_;
+
 
         // Reference to the particle integration schemes
 
@@ -464,6 +485,31 @@ public:
                     surfaceFilm();
 
 
+                //- Return const access to the packing model
+                inline const PackingModel<KinematicCloud<CloudType>>&
+                    packingModel() const;
+
+                //- Return a reference to the packing model
+                inline PackingModel<KinematicCloud<CloudType>>&
+                    packingModel();
+
+                //- Return const access to the damping model
+                inline const DampingModel<KinematicCloud<CloudType>>&
+                    dampingModel() const;
+
+                //- Return a reference to the damping model
+                inline DampingModel<KinematicCloud<CloudType>>&
+                    dampingModel();
+
+                //- Return const access to the isotropy model
+                inline const IsotropyModel<KinematicCloud<CloudType>>&
+                    isotropyModel() const;
+
+                //- Return a reference to the isotropy model
+                inline IsotropyModel<KinematicCloud<CloudType>>&
+                    isotropyModel();
+
+
             // Integration schemes
 
                 //-Return reference to velocity integration
@@ -590,7 +636,10 @@ public:
             void scaleSources();
 
             //- Pre-evolve
-            void preEvolve(const typename parcelType::trackingData& td);
+            void preEvolve
+            (
+                const typename parcelType::trackingData& td
+            );
 
             //- Evolve the cloud
             void evolve();

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -255,6 +255,54 @@ Foam::KinematicCloud<CloudType>::surfaceFilm()
 }
 
 
+template<class CloudType>
+inline const Foam::PackingModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::packingModel() const
+{
+    return *packingModel_;
+}
+
+
+template<class CloudType>
+inline Foam::PackingModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::packingModel()
+{
+    return *packingModel_;
+}
+
+
+template<class CloudType>
+inline const Foam::DampingModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::dampingModel() const
+{
+    return *dampingModel_;
+}
+
+
+template<class CloudType>
+inline Foam::DampingModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::dampingModel()
+{
+    return *dampingModel_;
+}
+
+
+template<class CloudType>
+inline const Foam::IsotropyModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::isotropyModel() const
+{
+    return *isotropyModel_;
+}
+
+
+template<class CloudType>
+inline Foam::IsotropyModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::isotropyModel()
+{
+    return *isotropyModel_;
+}
+
+
 template<class CloudType>
 inline const Foam::integrationScheme&
 Foam::KinematicCloud<CloudType>::UIntegrator() const

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

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -90,6 +91,28 @@ void Foam::KinematicParcel<ParcelType>::calcDispersion
 }
 
 
+template<class ParcelType>
+template<class TrackCloudType>
+void Foam::KinematicParcel<ParcelType>::calcUCorrection
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar dt
+)
+{
+    typename TrackCloudType::parcelType& p =
+        static_cast<typename TrackCloudType::parcelType&>(*this);
+
+    this->UCorrect_ = Zero;
+
+    this->UCorrect_ =
+        cloud.dampingModel().velocityCorrection(p, dt);
+
+    this->UCorrect_ +=
+        cloud.packingModel().velocityCorrection(p, dt);
+}
+
+
 template<class ParcelType>
 template<class TrackCloudType>
 void Foam::KinematicParcel<ParcelType>::cellValueSourceCorrection
@@ -141,6 +164,7 @@ void Foam::KinematicParcel<ParcelType>::calc
     this->U_ =
         calcVelocity(cloud, td, dt, Re, td.muc(), mass0, Su, dUTrans, Spu);
 
+    this->U_ += this->UCorrect_;
 
     // Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -249,7 +273,8 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(p.rho_),
     age_(p.age_),
     tTurb_(p.tTurb_),
-    UTurb_(p.UTurb_)
+    UTurb_(p.UTurb_),
+    UCorrect_(p.UCorrect_)
 {}
 
 
@@ -270,7 +295,8 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(p.rho_),
     age_(p.age_),
     tTurb_(p.tTurb_),
-    UTurb_(p.UTurb_)
+    UTurb_(p.UTurb_),
+    UCorrect_(p.UCorrect_)
 {}
 
 
@@ -335,6 +361,7 @@ bool Foam::KinematicParcel<ParcelType>::move
 
         const scalar dt = (p.stepFraction() - sfrac)*trackTime;
 
+
         // Avoid problems with extremely small timesteps
         if (dt > ROOTVSMALL)
         {
@@ -348,6 +375,8 @@ bool Foam::KinematicParcel<ParcelType>::move
                 p.cellValueSourceCorrection(cloud, ttd, dt);
             }
 
+            p.calcUCorrection(cloud, ttd, dt);
+
             p.calc(cloud, ttd, dt);
         }
 
@@ -362,7 +391,7 @@ bool Foam::KinematicParcel<ParcelType>::move
 
         if (p.active() && p.onFace() && ttd.keepParticle)
         {
-            p.hitFace(s, cloud, ttd);
+            p.hitFace(f*s - d, cloud, ttd);
         }
     }
 

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

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,8 +55,8 @@ SourceFiles
 #include "autoPtr.H"
 #include "interpolation.H"
 #include "demandDrivenEntry.H"
-
+#include "labelFieldIOField.H"
-// #include "ParticleForceList.H" // TODO
+#include "vectorFieldIOField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,6 +66,9 @@ namespace Foam
 template<class ParcelType>
 class KinematicParcel;
 
+template<class Type>
+class AveragingMethod;
+
 // Forward declaration of friend functions
 
 template<class ParcelType>
@@ -84,11 +87,12 @@ class KinematicParcel
 :
     public ParcelType
 {
-    // Private data
+    // Private Data
 
         //- Number of particle tracking attempts before we assume that it stalls
         static label maxTrackAttempts;
 
+
 public:
 
     //- Size in bytes of the fields
@@ -100,7 +104,7 @@ public:
     {
     protected:
 
-        // Protected data
+        // Protected Data
 
             //- Constant properties dictionary
             const dictionary dict_;
@@ -108,7 +112,7 @@ public:
 
     private:
 
-        // Private data
+        // Private Data
 
             //- Parcel type id - used for post-processing to flag the type
             //- of parcels issued by this cloud
@@ -138,7 +142,7 @@ public:
             constantProperties(const dictionary& parentDict);
 
 
-        // Member functions
+        // Member Functions
 
             //- Return const access to the constant properties dictionary
             inline const dictionary& dict() const;
@@ -173,7 +177,7 @@ public:
 
     private:
 
-        // Private data
+        // Private Data
 
             // Interpolators for continuous phase fields
 
@@ -199,6 +203,30 @@ public:
                 scalar muc_;
 
 
+            // MPPIC Averages
+
+                //- Volume average
+                autoPtr<AveragingMethod<scalar>> volumeAverage_;
+
+                //- Radius average [ volume^(1/3) ]
+                autoPtr<AveragingMethod<scalar>> radiusAverage_;
+
+                //- Density average
+                autoPtr<AveragingMethod<scalar>> rhoAverage_;
+
+                //- Velocity average
+                autoPtr<AveragingMethod<vector>> uAverage_;
+
+                //- Magnitude velocity squared average
+                autoPtr<AveragingMethod<scalar>> uSqrAverage_;
+
+                //- Frequency average
+                autoPtr<AveragingMethod<scalar>> frequencyAverage_;
+
+                //- Mass average
+                autoPtr<AveragingMethod<scalar>> massAverage_;
+
+
             //- Local gravitational or other body-force acceleration
             const vector& g_;
 
@@ -220,7 +248,7 @@ public:
             );
 
 
-        // Member functions
+        // Member Functions
 
             //- Return const access to the interpolator for continuous
             //- phase density field
@@ -260,12 +288,17 @@ public:
 
             //- Return access to the part of the tracking operation taking place
             inline trackPart& part();
+
+            //- Update the MPPIC averages
+            template<class TrackCloudType>
+            inline void updateAverages(const TrackCloudType& cloud);
+
     };
 
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         // Parcel properties
 
@@ -300,6 +333,9 @@ protected:
             //- Turbulent velocity fluctuation [m/s]
             vector UTurb_;
 
+            //- Velocity correction due to collisions MPPIC [m/s]
+            vector UCorrect_;
+
 
     // Protected Member Functions
 
@@ -340,6 +376,7 @@ public:
           + " age"
           + " tTurb"
           + " (UTurbx UTurby UTurbz)"
+          + " (UCorrectx UCorrecty UCorrectz)"
         );
 
 
@@ -465,6 +502,9 @@ public:
             //- Return const access to turbulent velocity fluctuation
             inline const vector& UTurb() const;
 
+            //- Return const access to correction velocity
+            inline const vector& UCorrect() const;
+
 
         // Edit
 
@@ -498,6 +538,9 @@ public:
             //- Return access to turbulent velocity fluctuation
             inline vector& UTurb();
 
+            //- Return access to correction velocity
+            inline vector& UCorrect();
+
 
         // Helper functions
 
@@ -611,6 +654,15 @@ public:
                 const scalar dt
             );
 
+            //- Correct U following MP-PIC sub-models
+            template<class TrackCloudType>
+            void calcUCorrection
+            (
+                TrackCloudType& cloud,
+                trackingData& td,
+                const scalar dt
+            );
+
 
         // Tracking
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -304,6 +304,20 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::U()
 }
 
 
+template<class ParcelType>
+inline const Foam::vector& Foam::KinematicParcel<ParcelType>::UCorrect() const
+{
+    return UCorrect_;
+}
+
+
+template<class ParcelType>
+inline Foam::vector& Foam::KinematicParcel<ParcelType>::UCorrect()
+{
+    return UCorrect_;
+}
+
+
 template<class ParcelType>
 inline Foam::scalar& Foam::KinematicParcel<ParcelType>::rho()
 {

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

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,7 +67,8 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(0.0),
     age_(0.0),
     tTurb_(0.0),
-    UTurb_(Zero)
+    UTurb_(Zero),
+    UCorrect_(Zero)
 {
     if (readFields)
     {
@@ -82,7 +83,8 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
                 >> rho_
                 >> age_
                 >> tTurb_
-                >> UTurb_;
+                >> UTurb_
+                >> UCorrect_;
         }
         else if (!is.checkLabelSize<>() || !is.checkScalarSize<>())
         {
@@ -100,6 +102,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
             readRawScalar(is, &age_);
             readRawScalar(is, &tTurb_);
             readRawScalar(is, UTurb_.data(), vector::nComponents);
+            readRawScalar(is, UCorrect_.data(), vector::nComponents);
 
             is.endRawRead();
         }
@@ -191,6 +194,13 @@ void Foam::KinematicParcel<ParcelType>::readFields(CloudType& c)
     );
     c.checkFieldIOobject(c, UTurb);
 
+    IOField<vector> UCorrect
+    (
+        c.fieldIOobject("UCorrect", IOobject::MUST_READ),
+        valid
+    );
+    c.checkFieldIOobject(c, UCorrect);
+
     label i = 0;
 
     for (KinematicParcel<ParcelType>& p : c)
@@ -205,6 +215,7 @@ void Foam::KinematicParcel<ParcelType>::readFields(CloudType& c)
         p.age_ = age[i];
         p.tTurb_ = tTurb[i];
         p.UTurb_ = UTurb[i];
+        p.UCorrect_ = UCorrect[i];
 
         ++i;
     }
@@ -234,6 +245,7 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const CloudType& c)
     IOField<scalar> age(c.fieldIOobject("age", IOobject::NO_READ), np);
     IOField<scalar> tTurb(c.fieldIOobject("tTurb", IOobject::NO_READ), np);
     IOField<vector> UTurb(c.fieldIOobject("UTurb", IOobject::NO_READ), np);
+    IOField<vector> UCorrect(c.fieldIOobject("UCorrect", IOobject::NO_READ), np);
 
     label i = 0;
 
@@ -249,6 +261,7 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const CloudType& c)
         age[i] = p.age();
         tTurb[i] = p.tTurb();
         UTurb[i] = p.UTurb();
+        UCorrect[i] = p.UCorrect();
 
         ++i;
     }
@@ -263,6 +276,7 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const CloudType& c)
     age.write(valid);
     tTurb.write(valid);
     UTurb.write(valid);
+    UCorrect.write(valid);
 }
 
 
@@ -291,6 +305,7 @@ void Foam::KinematicParcel<ParcelType>::writeProperties
     writeProp("age", age_);
     writeProp("tTurb", tTurb_);
     writeProp("UTurb", UTurb_);
+    writeProp("UCorrect", UCorrect_);
 
     #undef writeProp
 }
@@ -318,6 +333,7 @@ void Foam::KinematicParcel<ParcelType>::readObjects
     const auto& age = cloud::lookupIOField<scalar>("age", obr);
     const auto& tTurb = cloud::lookupIOField<scalar>("tTurb", obr);
     const auto& UTurb = cloud::lookupIOField<vector>("UTurb", obr);
+    const auto& UCorrect = cloud::lookupIOField<vector>("UCorrect", obr);
 
     label i = 0;
 
@@ -333,6 +349,7 @@ void Foam::KinematicParcel<ParcelType>::readObjects
         p.age_ = age[i];
         p.tTurb_ = tTurb[i];
         p.UTurb_ = UTurb[i];
+        p.UCorrect_ = UCorrect[i];
 
         ++i;
     }
@@ -361,6 +378,7 @@ void Foam::KinematicParcel<ParcelType>::writeObjects
     auto& age = cloud::createIOField<scalar>("age", np, obr);
     auto& tTurb = cloud::createIOField<scalar>("tTurb", np, obr);
     auto&& UTurb = cloud::createIOField<vector>("UTurb", np, obr);
+    auto&& UCorrect = cloud::createIOField<vector>("UCorrect", np, obr);
 
     label i = 0;
 
@@ -376,6 +394,7 @@ void Foam::KinematicParcel<ParcelType>::writeObjects
         age[i] = p.age();
         tTurb[i] = p.tTurb();
         UTurb[i] = p.UTurb();
+        UCorrect[i] = p.UCorrect();
 
         ++i;
     }
@@ -403,7 +422,8 @@ Foam::Ostream& Foam::operator<<
             << token::SPACE << p.rho()
             << token::SPACE << p.age()
             << token::SPACE << p.tTurb()
-            << token::SPACE << p.UTurb();
+            << token::SPACE << p.UTurb()
+            << token::SPACE << p.UCorrect();
     }
     else
     {

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelTrackingDataI.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +26,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "AveragingMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class ParcelType>
 template<class TrackCloudType>
 inline Foam::KinematicParcel<ParcelType>::trackingData::trackingData
@@ -61,6 +66,106 @@ inline Foam::KinematicParcel<ParcelType>::trackingData::trackingData
     rhoc_(Zero),
     Uc_(Zero),
     muc_(Zero),
+
+    volumeAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":volumeAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    radiusAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":radiusAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    rhoAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":rhoAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    uAverage_
+    (
+        AveragingMethod<vector>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":uAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    uSqrAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":uSqrAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    frequencyAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":frequencyAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    massAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":massAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+
     g_(cloud.g().value()),
     part_(part)
 {}
@@ -158,4 +263,107 @@ Foam::KinematicParcel<ParcelType>::trackingData::part()
 }
 
 
+template<class ParcelType>
+template<class TrackCloudType>
+inline void Foam::KinematicParcel<ParcelType>::trackingData::
+updateAverages
+(
+    const TrackCloudType& cloud
+)
+{
+    // zero the sums
+    volumeAverage_() = 0;
+    radiusAverage_() = 0;
+    rhoAverage_() = 0;
+    uAverage_() = Zero;
+    uSqrAverage_() = 0;
+    frequencyAverage_() = 0;
+    massAverage_() = 0;
+
+    // temporary weights
+    autoPtr<AveragingMethod<scalar>> weightAveragePtr
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":weightAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    );
+    AveragingMethod<scalar>& weightAverage = weightAveragePtr();
+
+    // averaging sums
+    for (const typename TrackCloudType::parcelType& p : cloud)
+    {
+        const tetIndices tetIs = p.currentTetIndices();
+
+        const scalar m = p.nParticle()*p.mass();
+
+        volumeAverage_->add(p.coordinates(), tetIs, p.nParticle()*p.volume());
+        rhoAverage_->add(p.coordinates(), tetIs, m*p.rho());
+        uAverage_->add(p.coordinates(), tetIs, m*p.U());
+        massAverage_->add(p.coordinates(), tetIs, m);
+    }
+    volumeAverage_->average();
+    massAverage_->average();
+    rhoAverage_->average(*massAverage_);
+    uAverage_->average(*massAverage_);
+
+    // squared velocity deviation
+    for (const typename TrackCloudType::parcelType& p : cloud)
+    {
+        const tetIndices tetIs = p.currentTetIndices();
+
+        const vector u = uAverage_->interpolate(p.coordinates(), tetIs);
+
+        uSqrAverage_->add
+        (
+            p.coordinates(),
+            tetIs,
+            p.nParticle()*p.mass()*magSqr(p.U() - u)
+        );
+    }
+    uSqrAverage_->average(*massAverage_);
+
+    // sauter mean radius
+    radiusAverage_() = volumeAverage_();
+    weightAverage = 0;
+    for (const typename TrackCloudType::parcelType& p : cloud)
+    {
+        const tetIndices tetIs = p.currentTetIndices();
+
+        weightAverage.add
+        (
+            p.coordinates(),
+            tetIs,
+            p.nParticle()*pow(p.volume(), 2.0/3.0)
+        );
+    }
+    weightAverage.average();
+    radiusAverage_->average(weightAverage);
+
+    // collision frequency
+    weightAverage = 0;
+    for (const typename TrackCloudType::parcelType& p : cloud)
+    {
+        const tetIndices tetIs = p.currentTetIndices();
+
+        const scalar a = volumeAverage_->interpolate(p.coordinates(), tetIs);
+        const scalar r = radiusAverage_->interpolate(p.coordinates(), tetIs);
+        const vector u = uAverage_->interpolate(p.coordinates(), tetIs);
+
+        const scalar f = 0.75*a/pow3(r)*sqr(0.5*p.d() + r)*mag(p.U() - u);
+
+        frequencyAverage_->add(p.coordinates(), tetIs, p.nParticle()*f*f);
+
+        weightAverage.add(p.coordinates(), tetIs, p.nParticle()*f);
+    }
+    frequencyAverage_->average(weightAverage);
+}
+
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousReactingParcelSubmodels.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,6 +45,11 @@ License
 #include "makeReactingMultiphaseParcelCompositionModels.H"
 #include "makeReactingParcelPhaseChangeModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicHeterogeneousReactingCloud);
@@ -68,4 +73,9 @@ makeHeterogeneousReactingParcelHeterogeneousReactingModels
     basicHeterogeneousReactingCloud
 );
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicHeterogeneousReactingCloud);
+makeMPPICParcelIsotropyModels(basicHeterogeneousReactingCloud);
+makeMPPICParcelPackingModels(basicHeterogeneousReactingCloud);
+
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicKinematicCollidingParcel/makeBasicKinematicCollidingParcelSubmodels.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,6 +39,11 @@ License
 #include "makeParcelStochasticCollisionModels.H"
 #include "makeParcelSurfaceFilmModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicKinematicCollidingCloud);
@@ -51,5 +57,9 @@ makeParcelPatchInteractionModels(basicKinematicCollidingCloud);
 makeParcelStochasticCollisionModels(basicKinematicCollidingCloud);
 makeParcelSurfaceFilmModels(basicKinematicCollidingCloud);
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicKinematicCollidingCloud);
+makeMPPICParcelIsotropyModels(basicKinematicCollidingCloud);
+makeMPPICParcelPackingModels(basicKinematicCollidingCloud);
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/makeBasicKinematicMPPICParcelSubmodels.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,9 +57,15 @@ makeParcelStochasticCollisionModels(basicKinematicMPPICCloud);
 makeParcelSurfaceFilmModels(basicKinematicMPPICCloud);
 
 // MPPIC sub-models
+// WIP: These models are defined in Kinematic and MPPIC clouds temporarily
 makeMPPICParcelDampingModels(basicKinematicMPPICCloud);
+makeMPPICCloudParcelDampingModels(basicKinematicMPPICCloud);
+
 makeMPPICParcelIsotropyModels(basicKinematicMPPICCloud);
+makeMPPICCloudParcelIsotropyModels(basicKinematicMPPICCloud);
+
 makeMPPICParcelPackingModels(basicKinematicMPPICCloud);
+makeMPPICCloudParcelPackingModels(basicKinematicMPPICCloud);
 
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelSubmodels.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,6 +38,11 @@ License
 #include "makeParcelStochasticCollisionModels.H"
 #include "makeParcelSurfaceFilmModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicKinematicCloud);
@@ -49,5 +55,10 @@ makeParcelPatchInteractionModels(basicKinematicCloud);
 makeParcelStochasticCollisionModels(basicKinematicCloud);
 makeParcelSurfaceFilmModels(basicKinematicCloud);
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicKinematicCloud);
+makeMPPICParcelIsotropyModels(basicKinematicCloud);
+makeMPPICParcelPackingModels(basicKinematicCloud);
+
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/makeBasicReactingMultiphaseParcelSubmodels.C

@@ -49,6 +49,11 @@ License
 #include "makeReactingMultiphaseParcelDevolatilisationModels.H"
 #include "makeReactingMultiphaseParcelSurfaceReactionModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicReactingMultiphaseCloud);
@@ -84,5 +89,9 @@ makeReactingMultiphaseParcelSurfaceReactionModels
     basicReactingMultiphaseCloud
 );
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicReactingMultiphaseCloud);
+makeMPPICParcelIsotropyModels(basicReactingMultiphaseCloud);
+makeMPPICParcelPackingModels(basicReactingMultiphaseCloud);
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicReactingParcel/makeBasicReactingParcelSubmodels.C

@@ -45,6 +45,11 @@ License
 #include "makeReactingParcelCompositionModels.H"
 #include "makeReactingParcelPhaseChangeModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicReactingCloud);
@@ -64,5 +69,10 @@ makeParcelHeatTransferModels(basicReactingCloud);
 makeReactingParcelCompositionModels(basicReactingCloud);
 makeReactingParcelPhaseChangeModels(basicReactingCloud);
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicReactingCloud);
+makeMPPICParcelIsotropyModels(basicReactingCloud);
+makeMPPICParcelPackingModels(basicReactingCloud);
+
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/makeBasicThermoParcelSubmodels.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +41,11 @@ License
 // Thermodynamic
 #include "makeParcelHeatTransferModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicThermoCloud);
@@ -56,4 +62,10 @@ makeParcelSurfaceFilmModels(basicThermoCloud);
 makeParcelHeatTransferModels(basicThermoCloud);
 
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicThermoCloud);
+makeMPPICParcelIsotropyModels(basicThermoCloud);
+makeMPPICParcelPackingModels(basicThermoCloud);
+
+
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/include/makeMPPICParcelDampingModels.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,6 +45,13 @@ License
     makeDampingModelType(Relaxation, CloudType);
 
 
+#define makeMPPICCloudParcelDampingModels(CloudType)                           \
+                                                                               \
+    makeDampingModelMPPIC(CloudType);                                          \
+                                                                               \
+    makeDampingModelTypeMPPIC(NoDamping, CloudType);                           \
+    makeDampingModelTypeMPPIC(Relaxation, CloudType);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/lagrangian/intermediate/parcels/include/makeMPPICParcelIsotropyModels.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,6 +44,13 @@ License
     makeIsotropyModelType(Stochastic, CloudType);
 
 
+#define makeMPPICCloudParcelIsotropyModels(CloudType)                          \
+                                                                               \
+    makeIsotropyModelMPPIC(CloudType);                                         \
+                                                                               \
+    makeIsotropyModelTypeMPPIC(NoIsotropy, CloudType);                         \
+    makeIsotropyModelTypeMPPIC(Stochastic, CloudType);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/lagrangian/intermediate/parcels/include/makeMPPICParcelPackingModels.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,6 +45,15 @@ License
     makePackingModelType(Explicit, CloudType);                                 \
     makePackingModelType(Implicit, CloudType);
 
+
+#define makeMPPICCloudParcelPackingModels(CloudType)                           \
+                                                                               \
+    makePackingModelMPPIC(CloudType);                                          \
+                                                                               \
+    makePackingModelTypeMPPIC(NoPacking, CloudType);                           \
+    makePackingModelTypeMPPIC(Explicit, CloudType);                            \
+    makePackingModelTypeMPPIC(Implicit, CloudType);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif