From 571e2ce3c8a258e4e07622cf5bcd3323b8095d80 Mon Sep 17 00:00:00 2001
From: Thomas Lichtenegger <k3b1677@relise.jku.austriangrid.at>
Date: Thu, 8 Jul 2021 10:53:50 +0200
Subject: [PATCH] Turbulent dispersion model can read precalculated random
 displacements.

---
 .../turbulentDispersion/turbulentDispersion.C | 43 +++++++++++++------
 .../turbulentDispersion/turbulentDispersion.H |  6 +++
 2 files changed, 37 insertions(+), 12 deletions(-)

diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C
index a9d05327..a0084a1b 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C
@@ -66,6 +66,8 @@ turbulentDispersion::turbulentDispersion
     nutName_(propsDict_.lookupOrDefault<word>("nutName","nut")),
     usePreCalcK_(propsDict_.lookupOrDefault<bool>("usePreCalcK",false)),
     kName_(propsDict_.lookupOrDefault<word>("kName","k")),
+    usePreCalcDispField_(propsDict_.lookupOrDefault<bool>("usePreCalcDispField",false)),
+    dispFieldName_(propsDict_.lookupOrDefault<word>("dispFieldName","dispVarField")),
     nut_
     (   IOobject
         (
@@ -117,9 +119,9 @@ turbulentDispersion::turbulentDispersion
     }
     else existIgnoreCells_ = false;
 
-    if (usePreCalcNut_ && usePreCalcK_)
+    if ((usePreCalcNut_ && usePreCalcK_) || (usePreCalcNut_ && usePreCalcDispField_) || (usePreCalcDispField_ && usePreCalcK_))
     {
-        FatalError<< "Cannot use precalculated nut and k at the same time. Choose one." << abort(FatalError);
+        FatalError<< "Cannot use more than one precalculated nut, k and displacement fluctuations at the same time. Choose one." << abort(FatalError);
     }
 
     if (usePreCalcK_)
@@ -166,7 +168,7 @@ bool turbulentDispersion::ignoreCell(label cell) const
 
 void turbulentDispersion::setForce() const
 {
-    if (!usePreCalcNut_ && !usePreCalcK_)
+    if (!usePreCalcNut_ && !usePreCalcK_ && !usePreCalcDispField_)
     {
         nut_ = particleCloud_.turbulence().nut()();
     }
@@ -175,17 +177,23 @@ void turbulentDispersion::setForce() const
         volScalarField& nutRef (const_cast<volScalarField&>(mesh_.lookupObject<volScalarField> (nutName_)));
         nut_ = nutRef;
     }
-    else
+    else if (usePreCalcK_)
     {
         volScalarField& kRef (const_cast<volScalarField&>(mesh_.lookupObject<volScalarField> (kName_)));
         nut_ = Ck_ * delta_ * sqrt(kRef);
     }
+    else
+    {
+        volVectorField& dispFieldRef (const_cast<volVectorField&>(mesh_.lookupObject<volVectorField> (dispFieldName_)));
+        dispField_ = &dispFieldRef;
+    }
 
     label cellI = -1;
     label patchID = -1;
     label faceIGlobal = -1;
     scalar flucProjection = 0.0;
     scalar D = 0.0;
+    scalar randScalar = 0.0;
     vector faceINormal = vector::zero;
     vector flucU = vector::zero;
     vector position = vector::zero;
@@ -198,20 +206,31 @@ void turbulentDispersion::setForce() const
         cellI = particleCloud_.cellIDs()[index][0];
         if (cellI > -1 && !ignoreCell(cellI))
         {
-            if (interpolate_)
+            if (usePreCalcDispField_)
             {
-                position = particleCloud_.position(index);
-                D = nutInterpolator_.interpolate(position,cellI) / turbulentSchmidtNumber_;
+                for (label comp=0; comp<3; comp++)
+                {
+                    randScalar=ranGen_.scalarNormal();
+                    flucU.component(comp) = randScalar * (*dispField_)[cellI].component(comp);
+                }
             }
             else
             {
-                D = nut_[cellI] / turbulentSchmidtNumber_;
+                if (interpolate_)
+                {
+                    position = particleCloud_.position(index);
+                    D = nutInterpolator_.interpolate(position,cellI) / turbulentSchmidtNumber_;
+                }
+                else
+                {
+                    D = nut_[cellI] / turbulentSchmidtNumber_;
+                }
+
+                // include concentration dependence on the diffusivity at this point if necessary
+
+                flucU=unitFlucDir()*Foam::sqrt(6.0*D/dt_);
             }
 
-            // include concentration dependence on the diffusivity at this point if necessary
-
-            flucU=unitFlucDir()*Foam::sqrt(6.0*D/dt_);
-
             // prevent particles being pushed through walls by regulating velocity fluctuations
             // check if cell is adjacent to wall and remove corresponding components
             if (wallIndicatorField_[cellI] > 0.5)
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.H b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.H
index 609cb50a..c324b274 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.H
@@ -66,6 +66,12 @@ protected:
 
     word kName_;
 
+    bool usePreCalcDispField_;
+
+    word dispFieldName_;
+
+    mutable volVectorField* dispField_;
+
     mutable volScalarField nut_;
 
     mutable volScalarField wallIndicatorField_;