lagrangian::BrownianMotionForce: Changed from a cubic to a spherical distribution

See also: StochasticDispersionRAS
Resolves bug-report http://bugs.openfoam.org/view.php?id=2153

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

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_;
 }

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;
 }