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-2018 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_(readScalar(this->coeffDict().lookup("weCorrCoeff"))),
39  weBuCrit_(readScalar(this->coeffDict().lookup("weBuCrit"))),
40  weBuBag_(readScalar(this->coeffDict().lookup("weBuBag"))),
41  weBuMM_(readScalar(this->coeffDict().lookup("weBuMM"))),
42  ohnCoeffCrit_(readScalar(this->coeffDict().lookup("ohnCoeffCrit"))),
43  ohnCoeffBag_(readScalar(this->coeffDict().lookup("ohnCoeffBag"))),
44  ohnCoeffMM_(readScalar(this->coeffDict().lookup("ohnCoeffMM"))),
45  ohnExpCrit_(readScalar(this->coeffDict().lookup("ohnExpCrit"))),
46  ohnExpBag_(readScalar(this->coeffDict().lookup("ohnExpBag"))),
47  ohnExpMM_(readScalar(this->coeffDict().lookup("ohnExpMM"))),
48  cInit_(readScalar(this->coeffDict().lookup("Cinit"))),
49  c1_(readScalar(this->coeffDict().lookup("C1"))),
50  c2_(readScalar(this->coeffDict().lookup("C2"))),
51  c3_(readScalar(this->coeffDict().lookup("C3"))),
52  cExp1_(readScalar(this->coeffDict().lookup("Cexp1"))),
53  cExp2_(readScalar(this->coeffDict().lookup("Cexp2"))),
54  cExp3_(readScalar(this->coeffDict().lookup("Cexp3"))),
55  weConst_(readScalar(this->coeffDict().lookup("Weconst"))),
56  weCrit1_(readScalar(this->coeffDict().lookup("Wecrit1"))),
57  weCrit2_(readScalar(this->coeffDict().lookup("Wecrit2"))),
58  coeffD_(readScalar(this->coeffDict().lookup("CoeffD"))),
59  onExpD_(readScalar(this->coeffDict().lookup("OnExpD"))),
60  weExpD_(readScalar(this->coeffDict().lookup("WeExpD"))),
61  mu_(readScalar(this->coeffDict().lookup("mu"))),
62  sigma_(readScalar(this->coeffDict().lookup("sigma"))),
63  d32Coeff_(readScalar(this->coeffDict().lookup("d32Coeff"))),
64  cDmaxBM_(readScalar(this->coeffDict().lookup("cDmaxBM"))),
65  cDmaxS_(readScalar(this->coeffDict().lookup("cDmaxS"))),
66  corePerc_(readScalar(this->coeffDict().lookup("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  // droplet deformation characteristic time
160 
161  scalar tChar = d/Urmag*sqrt(rho/rhoc);
162 
163  scalar tFirst = cInit_*tChar;
164 
165  scalar tSecond = 0;
166  scalar tCharSecond = 0;
167 
168  bool bag = false;
169  bool multimode = false;
170  bool shear = false;
171  bool success = false;
172 
173 
174  // update the droplet characteristic time
175  tc += dt;
176 
177  if (weGas > weConst_)
178  {
179  if (weGas < weCrit1_)
180  {
181  tCharSecond = c1_*pow((weGas - weConst_), cExp1_);
182  }
183  else if (weGas >= weCrit1_ && weGas <= weCrit2_)
184  {
185  tCharSecond = c2_*pow((weGas - weConst_), cExp2_);
186  }
187  else
188  {
189  tCharSecond = c3_*pow((weGas - weConst_), cExp3_);
190  }
191  }
192 
193  scalar weC = weBuCrit_*(1.0 + ohnCoeffCrit_*pow(ohnesorge, ohnExpCrit_));
194  scalar weB = weBuBag_*(1.0 + ohnCoeffBag_*pow(ohnesorge, ohnExpBag_));
195  scalar weMM = weBuMM_*(1.0 + ohnCoeffMM_*pow(ohnesorge, ohnExpMM_));
196 
197  if (weGas > weC && weGas < weB)
198  {
199  bag = true;
200  }
201 
202  if (weGas >= weB && weGas <= weMM)
203  {
204  multimode = true;
205  }
206 
207  if (weGas > weMM)
208  {
209  shear = true;
210  }
211 
212  tSecond = tCharSecond*tChar;
213 
214  scalar tBreakUP = tFirst + tSecond;
215  if (tc > tBreakUP)
216  {
217  scalar d32 = coeffD_*d*pow(ohnesorge, onExpD_)*pow(weGasCorr, weExpD_);
218 
219  if (bag || multimode)
220  {
221  scalar d05 = d32Coeff_ * d32;
222 
223  scalar x = 0.0;
224  scalar yGuess = 0.0;
225  scalar dGuess = 0.0;
226 
227  while(!success)
228  {
229  x = cDmaxBM_*rndGen.sample01<scalar>();
230  dGuess = sqr(x)*d05;
231  yGuess = rndGen.sample01<scalar>();
232 
233  scalar p =
234  x
235  /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
236  *exp(-0.5*sqr((x - mu_)/sigma_));
237 
238  if (yGuess < p)
239  {
240  success = true;
241  }
242  }
243 
244  d = dGuess;
245  tc = 0.0;
246  }
247 
248  if (shear)
249  {
250  scalar dC = weConst_*sigma/(rhoc*sqr(Urmag));
251  scalar d32Red = 4.0*(d32*dC)/(5.0*dC - d32);
252 
253  scalar d05 = d32Coeff_ * d32Red;
254 
255  scalar x = 0.0;
256  scalar yGuess = 0.0;
257  scalar dGuess = 0.0;
258 
259  while(!success)
260  {
261  x = cDmaxS_*rndGen.sample01<scalar>();
262  dGuess = sqr(x)*d05;
263  yGuess = rndGen.sample01<scalar>();
264 
265  scalar p =
266  x
267  /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
268  *exp(-0.5*sqr((x - mu_)/sigma_));
269 
270  if (yGuess<p)
271  {
272  success = true;
273  }
274  }
275 
276  d = dC;
277  dChild = dGuess;
278  massChild = corePerc_*mass0;
279  mass -= massChild;
280 
281  addChild = true;
282  // reset timer
283  tc = 0.0;
284  }
285 
286  // correct nParticle to conserve mass
287  nParticle = mass/(rhopi6*pow3(d));
288  }
289 
290  return addChild;
291 }
292 
293 
294 // ************************************************************************* //
#define readScalar
Definition: doubleScalar.C:38
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:137
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
const CloudType & owner() const
Return const access to the owner cloud.
Secondary Breakup Model to take account of the different breakup regimes, bag, molutimode, 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
stressControl lookup("compactNormalStress") >> compactNormalStress
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
const dimensionedScalar mu
Atomic mass unit.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
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