From ee9e7cf0375147af1e9d57a039ad4f241bad8e4f Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sat, 23 Jul 2016 11:53:48 +0100
Subject: [PATCH] lagrangian::BrownianMotionForce: Changed from a cubic to a
 spherical distribution See also: StochasticDispersionRAS Resolves bug-report
 http://bugs.openfoam.org/view.php?id=2153

---
 .../Templates/KinematicCloud/KinematicCloud.H |  4 +--
 .../KinematicCloud/KinematicCloudI.H          |  2 +-
 .../BrownianMotion/BrownianMotionForce.C      | 29 ++++++++++++++-----
 3 files changed, 25 insertions(+), 10 deletions(-)

diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index 20916540674..42c3aee5f15 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -160,7 +160,7 @@ protected:
         const dictionary subModelProperties_;
 
         //- Random number generator - used by some injection routines
-        cachedRandom rndGen_;
+        mutable cachedRandom rndGen_;
 
         //- Cell occupancy information for each parcel, (demand driven)
         autoPtr<List<DynamicList<parcelType*>>> cellOccupancyPtr_;
@@ -363,7 +363,7 @@ public:
             // Cloud data
 
                 //- Return reference to the random object
-                inline cachedRandom& rndGen();
+                inline cachedRandom& rndGen() const;
 
                 //- Return the cell occupancy information for each
                 //  parcel, non-const access, the caller is
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index 37a1a43d176..53b6c713160 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -356,7 +356,7 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
 
 
 template<class CloudType>
-inline Foam::cachedRandom& Foam::KinematicCloud<CloudType>::rndGen()
+inline Foam::cachedRandom& Foam::KinematicCloud<CloudType>::rndGen() const
 {
     return rndGen_;
 }
diff --git a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
index af762752cf7..313cbf0b171 100644
--- a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
+++ b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
@@ -193,13 +193,28 @@ Foam::forceSuSp Foam::BrownianMotionForce<CloudType>::calcCoupled
         f = mass*sqrt(mathematical::pi*s0/dt);
     }
 
-    const scalar sqrt2 = sqrt(2.0);
-    for (direction dir = 0; dir < vector::nComponents; dir++)
-    {
-        const scalar x = rndGen_.sample01<scalar>();
-        const scalar eta = sqrt2*erfInv(2*x - 1.0);
-        value.Su()[dir] = f*eta;
-    }
+
+    // To generate a cubic distribution (3 independent directions) :
+    // const scalar sqrt2 = sqrt(2.0);
+    // for (direction dir = 0; dir < vector::nComponents; dir++)
+    // {
+    //     const scalar x = rndGen_.sample01<scalar>();
+    //     const scalar eta = sqrt2*erfInv(2*x - 1.0);
+    //     value.Su()[dir] = f*eta;
+    // }
+
+
+    // To generate a spherical distribution:
+
+    cachedRandom& rnd = this->owner().rndGen();
+
+    const scalar theta = rnd.sample01<scalar>()*twoPi;
+    const scalar u = 2*rnd.sample01<scalar>() - 1;
+
+    const scalar a = sqrt(1 - sqr(u));
+    const vector dir(a*cos(theta), a*sin(theta), u);
+
+    value.Su() = f*mag(rnd.GaussNormal<scalar>())*dir;
 
     return value;
 }
-- 
GitLab