All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SHF.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
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 
26 #include "SHF.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class CloudType>
32 (
33  const dictionary& dict,
34  CloudType& owner
35 )
36 :
37  BreakupModel<CloudType>(dict, owner, typeName),
38  weCorrCoeff_(this->coeffDict().template lookup<scalar>("weCorrCoeff")),
39  weBuCrit_(this->coeffDict().template lookup<scalar>("weBuCrit")),
40  weBuBag_(this->coeffDict().template lookup<scalar>("weBuBag")),
41  weBuMM_(this->coeffDict().template lookup<scalar>("weBuMM")),
42  ohnCoeffCrit_(this->coeffDict().template lookup<scalar>("ohnCoeffCrit")),
43  ohnCoeffBag_(this->coeffDict().template lookup<scalar>("ohnCoeffBag")),
44  ohnCoeffMM_(this->coeffDict().template lookup<scalar>("ohnCoeffMM")),
45  ohnExpCrit_(this->coeffDict().template lookup<scalar>("ohnExpCrit")),
46  ohnExpBag_(this->coeffDict().template lookup<scalar>("ohnExpBag")),
47  ohnExpMM_(this->coeffDict().template lookup<scalar>("ohnExpMM")),
48  cInit_(this->coeffDict().template lookup<scalar>("Cinit")),
49  c1_(this->coeffDict().template lookup<scalar>("C1")),
50  c2_(this->coeffDict().template lookup<scalar>("C2")),
51  c3_(this->coeffDict().template lookup<scalar>("C3")),
52  cExp1_(this->coeffDict().template lookup<scalar>("Cexp1")),
53  cExp2_(this->coeffDict().template lookup<scalar>("Cexp2")),
54  cExp3_(this->coeffDict().template lookup<scalar>("Cexp3")),
55  weConst_(this->coeffDict().template lookup<scalar>("Weconst")),
56  weCrit1_(this->coeffDict().template lookup<scalar>("Wecrit1")),
57  weCrit2_(this->coeffDict().template lookup<scalar>("Wecrit2")),
58  coeffD_(this->coeffDict().template lookup<scalar>("CoeffD")),
59  onExpD_(this->coeffDict().template lookup<scalar>("OnExpD")),
60  weExpD_(this->coeffDict().template lookup<scalar>("WeExpD")),
61  mu_(this->coeffDict().template lookup<scalar>("mu")),
62  sigma_(this->coeffDict().template lookup<scalar>("sigma")),
63  d32Coeff_(this->coeffDict().template lookup<scalar>("d32Coeff")),
64  cDmaxBM_(this->coeffDict().template lookup<scalar>("cDmaxBM")),
65  cDmaxS_(this->coeffDict().template lookup<scalar>("cDmaxS")),
66  corePerc_(this->coeffDict().template lookup<scalar>("corePerc"))
67 {}
68 
69 
70 template<class CloudType>
72 :
73  BreakupModel<CloudType>(bum),
74  weCorrCoeff_(bum.weCorrCoeff_),
75  weBuCrit_(bum.weBuCrit_),
76  weBuBag_(bum.weBuBag_),
77  weBuMM_(bum.weBuMM_),
78  ohnCoeffCrit_(bum.ohnCoeffCrit_),
79  ohnCoeffBag_(bum.ohnCoeffBag_),
80  ohnCoeffMM_(bum.ohnCoeffMM_),
81  ohnExpCrit_(bum.ohnExpCrit_),
82  ohnExpBag_(bum.ohnExpBag_),
83  ohnExpMM_(bum.ohnExpMM_),
84  cInit_(bum.cInit_),
85  c1_(bum.c1_),
86  c2_(bum.c2_),
87  c3_(bum.c3_),
88  cExp1_(bum.cExp1_),
89  cExp2_(bum.cExp2_),
90  cExp3_(bum.cExp3_),
91  weConst_(bum.weConst_),
92  weCrit1_(bum.weCrit1_),
93  weCrit2_(bum.weCrit2_),
94  coeffD_(bum.coeffD_),
95  onExpD_(bum.onExpD_),
96  weExpD_(bum.weExpD_),
97  mu_(bum.mu_),
98  sigma_(bum.sigma_),
99  d32Coeff_(bum.d32Coeff_),
100  cDmaxBM_(bum.cDmaxBM_),
101  cDmaxS_(bum.cDmaxS_),
102  corePerc_(bum.corePerc_)
103 {}
104 
105 
106 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
107 
108 template<class CloudType>
110 {}
111 
112 
113 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
114 
115 template<class CloudType>
117 (
118  const scalar dt,
119  const vector& g,
120  scalar& d,
121  scalar& tc,
122  scalar& ms,
123  scalar& nParticle,
124  scalar& KHindex,
125  scalar& y,
126  scalar& yDot,
127  const scalar d0,
128  const scalar rho,
129  const scalar mu,
130  const scalar sigma,
131  const vector& U,
132  const scalar rhoc,
133  const scalar muc,
134  const vector& Urel,
135  const scalar Urmag,
136  const scalar tMom,
137  scalar& dChild,
138  scalar& massChild
139 )
140 {
141  Random& rndGen = this->owner().rndGen();
142 
143  bool addChild = false;
144 
145  scalar d03 = pow3(d);
146  scalar rhopi6 = rho*constant::mathematical::pi/6.0;
147  scalar mass0 = nParticle*rhopi6*d03;
148  scalar mass = mass0;
149 
150  scalar weGas = 0.5*rhoc*sqr(Urmag)*d/sigma;
151  scalar weLiquid = 0.5*rho*sqr(Urmag)*d/sigma;
152 
153  // correct the Reynolds number. Reitz is using radius instead of diameter
154  scalar reLiquid = 0.5*Urmag*d/mu;
155  scalar ohnesorge = sqrt(weLiquid)/(reLiquid + vSmall);
156 
157  scalar weGasCorr = weGas/(1.0 + weCorrCoeff_*ohnesorge);
158 
159  // update the droplet characteristic time
160  tc += dt;
161 
162  // droplet deformation characteristic rate
163  scalar rChar = Urmag/d*sqrt(rhoc/rho);
164 
165  // return if the characteristic deformation rate is too low for the
166  // following modelling to be calculable
167  if (tc*rChar < small)
168  {
169  return false;
170  }
171 
172  // droplet deformation characteristic time
173  scalar tChar = 1/rChar;
174 
175  scalar tFirst = cInit_*tChar;
176 
177  scalar tSecond = 0;
178  scalar tCharSecond = 0;
179 
180  bool bag = false;
181  bool multimode = false;
182  bool shear = false;
183  bool success = false;
184 
185 
186  if (weGas > weConst_)
187  {
188  if (weGas < weCrit1_)
189  {
190  tCharSecond = c1_*pow((weGas - weConst_), cExp1_);
191  }
192  else if (weGas >= weCrit1_ && weGas <= weCrit2_)
193  {
194  tCharSecond = c2_*pow((weGas - weConst_), cExp2_);
195  }
196  else
197  {
198  tCharSecond = c3_*pow((weGas - weConst_), cExp3_);
199  }
200  }
201 
202  scalar weC = weBuCrit_*(1.0 + ohnCoeffCrit_*pow(ohnesorge, ohnExpCrit_));
203  scalar weB = weBuBag_*(1.0 + ohnCoeffBag_*pow(ohnesorge, ohnExpBag_));
204  scalar weMM = weBuMM_*(1.0 + ohnCoeffMM_*pow(ohnesorge, ohnExpMM_));
205 
206  if (weGas > weC && weGas < weB)
207  {
208  bag = true;
209  }
210 
211  if (weGas >= weB && weGas <= weMM)
212  {
213  multimode = true;
214  }
215 
216  if (weGas > weMM)
217  {
218  shear = true;
219  }
220 
221  tSecond = tCharSecond*tChar;
222 
223  scalar tBreakUP = tFirst + tSecond;
224  if (tc > tBreakUP)
225  {
226  scalar d32 = coeffD_*d*pow(ohnesorge, onExpD_)*pow(weGasCorr, weExpD_);
227 
228  if (bag || multimode)
229  {
230  scalar d05 = d32Coeff_ * d32;
231 
232  scalar x = 0.0;
233  scalar yGuess = 0.0;
234  scalar dGuess = 0.0;
235 
236  while(!success)
237  {
238  x = cDmaxBM_*rndGen.sample01<scalar>();
239  dGuess = sqr(x)*d05;
240  yGuess = rndGen.sample01<scalar>();
241 
242  scalar p =
243  x
244  /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
245  *exp(-0.5*sqr((x - mu_)/sigma_));
246 
247  if (yGuess < p)
248  {
249  success = true;
250  }
251  }
252 
253  d = dGuess;
254  tc = 0.0;
255  }
256 
257  if (shear)
258  {
259  scalar dC = weConst_*sigma/(rhoc*sqr(Urmag));
260  scalar d32Red = 4.0*(d32*dC)/(5.0*dC - d32);
261 
262  scalar d05 = d32Coeff_ * d32Red;
263 
264  scalar x = 0.0;
265  scalar yGuess = 0.0;
266  scalar dGuess = 0.0;
267 
268  while(!success)
269  {
270  x = cDmaxS_*rndGen.sample01<scalar>();
271  dGuess = sqr(x)*d05;
272  yGuess = rndGen.sample01<scalar>();
273 
274  scalar p =
275  x
276  /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
277  *exp(-0.5*sqr((x - mu_)/sigma_));
278 
279  if (yGuess<p)
280  {
281  success = true;
282  }
283  }
284 
285  d = dC;
286  dChild = dGuess;
287  massChild = corePerc_*mass0;
288  mass -= massChild;
289 
290  addChild = true;
291  // reset timer
292  tc = 0.0;
293  }
294 
295  // correct nParticle to conserve mass
296  nParticle = mass/(rhopi6*pow3(d));
297  }
298 
299  return addChild;
300 }
301 
302 
303 // ************************************************************************* //
virtual bool update(const scalar dt, const vector &g, scalar &d, scalar &tc, scalar &ms, scalar &nParticle, scalar &KHindex, scalar &y, scalar &yDot, const scalar d0, const scalar rho, const scalar mu, const scalar sigma, const vector &U, const scalar rhoc, const scalar muc, const vector &Urel, const scalar Urmag, const scalar tMom, scalar &dChild, scalar &massChild)
Update the parcel properties.
Definition: SHF.C:117
dictionary dict
SHF(const dictionary &, CloudType &)
Construct from dictionary.
Definition: SHF.C:32
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionedScalar & sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
dimensionedScalar sqrt(const dimensionedScalar &ds)
const CloudType & owner() const
Return const access to the owner cloud.
Secondary Breakup Model to take account of the different breakup regimes, bag, solutionmode, shear....
Definition: SHF.H:54
Random rndGen(label(0))
Type sample01()
Advance the state and return a sample of a given type from a.
Definition: RandomI.H:70
dimensionedScalar exp(const dimensionedScalar &ds)
Random & rndGen()
Return references to the random object.
Definition: DSMCCloudI.H:121
Random number generator.
Definition: Random.H:57
const scalar twoPi(2 *pi)
bool success
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
const dimensionedScalar & mu
Atomic mass unit.
Templated break-up model class.
Definition: SprayCloud.H:50
volScalarField & p
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69
virtual ~SHF()
Destructor.
Definition: SHF.C:109