StochasticDispersionRAS: Corrected spherical distribution of UTurb direction

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1650
This commit is contained in:
Henry 2015-04-08 13:08:08 +01:00
parent d7d91261a9
commit a869754bc4
2 changed files with 26 additions and 24 deletions

View file

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -24,6 +24,9 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "StochasticDispersionRAS.H" #include "StochasticDispersionRAS.H"
#include "constants.H"
using namespace Foam::constant::mathematical;
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -89,29 +92,34 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
if (tTurb > tTurbLoc) if (tTurb > tTurbLoc)
{ {
tTurb = 0.0; tTurb = 0;
scalar sigma = sqrt(2.0*k/3.0); const scalar sigma = sqrt(2*k/3.0);
vector dir = 2.0*rnd.sample01<vector>() - vector::one;
dir /= mag(dir) + SMALL; // Calculate a random direction dir distributed uniformly
// in spherical coordinates
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);
// Numerical Recipes... Ch. 7. Random Numbers... // Numerical Recipes... Ch. 7. Random Numbers...
scalar x1 = 0.0; scalar x1 = 0;
scalar x2 = 0.0; scalar x2 = 0;
scalar rsq = 10.0; scalar rsq = 10;
while ((rsq > 1.0) || (rsq == 0.0)) while ((rsq > 1) || (rsq == 0))
{ {
x1 = 2.0*rnd.sample01<scalar>() - 1.0; x1 = 2*rnd.sample01<scalar>() - 1;
x2 = 2.0*rnd.sample01<scalar>() - 1.0; x2 = 2*rnd.sample01<scalar>() - 1;
rsq = x1*x1 + x2*x2; rsq = x1*x1 + x2*x2;
} }
scalar fac = sqrt(-2.0*log(rsq)/rsq); scalar fac = sqrt(-2*log(rsq)/rsq);
fac *= mag(x1); fac *= mag(x1);
UTurb = sigma*fac*dir; UTurb = sigma*fac*dir;
} }
} }
else else

View file

@ -99,17 +99,11 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
// Calculate a random direction dir distributed uniformly // Calculate a random direction dir distributed uniformly
// in spherical coordinates // in spherical coordinates
const scalar theta = rnd.sample01<scalar>()*pi; const scalar theta = rnd.sample01<scalar>()*twoPi;
const scalar phi = rnd.sample01<scalar>()*twoPi; const scalar u = 2*rnd.sample01<scalar>() - 1;
// Optimising compilers will use the sincos function const scalar a = sqrt(1 - sqr(u));
const scalar sinTheta = sin(theta); const vector dir(a*cos(theta), a*sin(theta), u);
const scalar cosTheta = cos(theta);
const scalar sinPhi = sin(phi);
const scalar cosPhi = cos(phi);
const vector dir(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta);
// Numerical Recipes... Ch. 7. Random Numbers... // Numerical Recipes... Ch. 7. Random Numbers...
scalar x1 = 0; scalar x1 = 0;