Commit 46d69e1d authored by Henry Weller's avatar Henry Weller
Browse files

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

See also: StochasticDispersionRAS
Resolves bug-report http://bugs.openfoam.org/view.php?id=2153
parent 96645225
......@@ -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
......
......@@ -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_;
}
......
......@@ -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;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment