FreeFOAM The Cross-Platform CFD Toolkit
GradientDispersionRAS.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 :
37  DispersionRASModel<CloudType>(dict, owner),
38  gradkPtr_(NULL)
39 {}
40 
41 
42 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
43 
44 template<class CloudType>
46 {
47  cacheFields(false);
48 }
49 
50 
51 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
52 
53 template<class CloudType>
55 {
56  return true;
57 }
58 
59 
60 template<class CloudType>
62 {
64 
65  if (store)
66  {
67  gradkPtr_ = fvc::grad(*this->kPtr_).ptr();
68  }
69  else
70  {
71  if (gradkPtr_)
72  {
73  delete gradkPtr_;
74  gradkPtr_ = NULL;
75  }
76  }
77 }
78 
79 
80 template<class CloudType>
82 (
83  const scalar dt,
84  const label celli,
85  const vector& U,
86  const vector& Uc,
87  vector& UTurb,
88  scalar& tTurb
89 )
90 {
91  const scalar cps = 0.16432;
92 
93  const volScalarField& k = *this->kPtr_;
94  const volScalarField& epsilon = *this->epsilonPtr_;
95  const volVectorField& gradk = *this->gradkPtr_;
96 
97  const scalar UrelMag = mag(U - Uc - UTurb);
98 
99  const scalar tTurbLoc = min
100  (
101  k[celli]/epsilon[celli],
102  cps*pow(k[celli], 1.5)/epsilon[celli]/(UrelMag + SMALL)
103  );
104 
105  // Parcel is perturbed by the turbulence
106  if (dt < tTurbLoc)
107  {
108  tTurb += dt;
109 
110  if (tTurb > tTurbLoc)
111  {
112  tTurb = 0.0;
113 
114  scalar sigma = sqrt(2.0*k[celli]/3.0);
115  vector dir = -gradk[celli]/(mag(gradk[celli]) + SMALL);
116 
117  // Numerical Recipes... Ch. 7. Random Numbers...
118  scalar x1 = 0.0;
119  scalar x2 = 0.0;
120  scalar rsq = 10.0;
121  while ((rsq > 1.0) || (rsq == 0.0))
122  {
123  x1 = 2.0*this->owner().rndGen().scalar01() - 1.0;
124  x2 = 2.0*this->owner().rndGen().scalar01() - 1.0;
125  rsq = x1*x1 + x2*x2;
126  }
127 
128  scalar fac = sqrt(-2.0*log(rsq)/rsq);
129 
130  // In 2D calculations the -grad(k) is always
131  // away from the axis of symmetry
132  // This creates a 'hole' in the spray and to
133  // prevent this we let x1 be both negative/positive
134  if (this->owner().mesh().nSolutionD() == 2)
135  {
136  fac *= x1;
137  }
138  else
139  {
140  fac *= mag(x1);
141  }
142 
143  UTurb = sigma*fac*dir;
144  }
145  }
146  else
147  {
148  tTurb = GREAT;
149  UTurb = vector::zero;
150  }
151 
152  return Uc + UTurb;
153 }
154 
155 
156 // ************************ vim: set sw=4 sts=4 et: ************************ //