FreeFOAM The Cross-Platform CFD Toolkit
StochasticDispersionRAS.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class CloudType>
32 (
33  const dictionary& dict,
34  CloudType& owner
35 )
36 :
38 {}
39 
40 
41 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
42 
43 template<class CloudType>
45 {}
46 
47 
48 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
49 
50 template<class CloudType>
52 {
53  return true;
54 }
55 
56 
57 template<class CloudType>
59 (
60  const scalar dt,
61  const label celli,
62  const vector& U,
63  const vector& Uc,
64  vector& UTurb,
65  scalar& tTurb
66 )
67 {
68  const scalar cps = 0.16432;
69 
70  const volScalarField& k = *this->kPtr_;
71  const volScalarField& epsilon = *this->epsilonPtr_;
72 
73  const scalar UrelMag = mag(U - Uc - UTurb);
74 
75  const scalar tTurbLoc = min
76  (
77  k[celli]/epsilon[celli],
78  cps*pow(k[celli], 1.5)/epsilon[celli]/(UrelMag + SMALL)
79  );
80 
81  // Parcel is perturbed by the turbulence
82  if (dt < tTurbLoc)
83  {
84  tTurb += dt;
85 
86  if (tTurb > tTurbLoc)
87  {
88  tTurb = 0.0;
89 
90  scalar sigma = sqrt(2.0*k[celli]/3.0);
91  vector dir = 2.0*this->owner().rndGen().vector01() - vector::one;
92  dir /= mag(dir) + SMALL;
93 
94  // Numerical Recipes... Ch. 7. Random Numbers...
95  scalar x1 = 0.0;
96  scalar x2 = 0.0;
97  scalar rsq = 10.0;
98  while ((rsq > 1.0) || (rsq == 0.0))
99  {
100  x1 = 2.0*this->owner().rndGen().scalar01() - 1.0;
101  x2 = 2.0*this->owner().rndGen().scalar01() - 1.0;
102  rsq = x1*x1 + x2*x2;
103  }
104 
105  scalar fac = sqrt(-2.0*log(rsq)/rsq);
106 
107  fac *= mag(x1);
108 
109  UTurb = sigma*fac*dir;
110 
111  }
112  }
113  else
114  {
115  tTurb = GREAT;
116  UTurb = vector::zero;
117  }
118 
119  return Uc + UTurb;
120 }
121 
122 
123 // ************************ vim: set sw=4 sts=4 et: ************************ //