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-2024 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  const label injectori,
138  scalar& dChild,
139  scalar& massChild
140 )
141 {
142  randomGenerator& rndGen = this->owner().rndGen();
143 
144  bool addChild = false;
145 
146  scalar d03 = pow3(d);
147  scalar rhopi6 = rho*constant::mathematical::pi/6.0;
148  scalar mass0 = nParticle*rhopi6*d03;
149  scalar mass = mass0;
150 
151  scalar weGas = 0.5*rhoc*sqr(Urmag)*d/sigma;
152  scalar weLiquid = 0.5*rho*sqr(Urmag)*d/sigma;
153 
154  // correct the Reynolds number. Reitz is using radius instead of diameter
155  scalar reLiquid = 0.5*Urmag*d/mu;
156  scalar ohnesorge = sqrt(weLiquid)/(reLiquid + vSmall);
157 
158  scalar weGasCorr = weGas/(1.0 + weCorrCoeff_*ohnesorge);
159 
160  // update the droplet characteristic time
161  tc += dt;
162 
163  // droplet deformation characteristic rate
164  scalar rChar = Urmag/d*sqrt(rhoc/rho);
165 
166  // return if the characteristic deformation rate is too low for the
167  // following modelling to be calculable
168  if (tc*rChar < small)
169  {
170  return false;
171  }
172 
173  // droplet deformation characteristic time
174  scalar tChar = 1/rChar;
175 
176  scalar tFirst = cInit_*tChar;
177 
178  scalar tSecond = 0;
179  scalar tCharSecond = 0;
180 
181  bool bag = false;
182  bool multimode = false;
183  bool shear = false;
184  bool success = false;
185 
186 
187  if (weGas > weConst_)
188  {
189  if (weGas < weCrit1_)
190  {
191  tCharSecond = c1_*pow((weGas - weConst_), cExp1_);
192  }
193  else if (weGas >= weCrit1_ && weGas <= weCrit2_)
194  {
195  tCharSecond = c2_*pow((weGas - weConst_), cExp2_);
196  }
197  else
198  {
199  tCharSecond = c3_*pow((weGas - weConst_), cExp3_);
200  }
201  }
202 
203  scalar weC = weBuCrit_*(1.0 + ohnCoeffCrit_*pow(ohnesorge, ohnExpCrit_));
204  scalar weB = weBuBag_*(1.0 + ohnCoeffBag_*pow(ohnesorge, ohnExpBag_));
205  scalar weMM = weBuMM_*(1.0 + ohnCoeffMM_*pow(ohnesorge, ohnExpMM_));
206 
207  if (weGas > weC && weGas < weB)
208  {
209  bag = true;
210  }
211 
212  if (weGas >= weB && weGas <= weMM)
213  {
214  multimode = true;
215  }
216 
217  if (weGas > weMM)
218  {
219  shear = true;
220  }
221 
222  tSecond = tCharSecond*tChar;
223 
224  scalar tBreakUP = tFirst + tSecond;
225  if (tc > tBreakUP)
226  {
227  scalar d32 = coeffD_*d*pow(ohnesorge, onExpD_)*pow(weGasCorr, weExpD_);
228 
229  if (bag || multimode)
230  {
231  scalar d05 = d32Coeff_ * d32;
232 
233  scalar x = 0.0;
234  scalar yGuess = 0.0;
235  scalar dGuess = 0.0;
236 
237  while(!success)
238  {
239  x = cDmaxBM_*rndGen.sample01<scalar>();
240  dGuess = sqr(x)*d05;
241  yGuess = rndGen.sample01<scalar>();
242 
243  scalar p =
244  x
245  /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
246  *exp(-0.5*sqr((x - mu_)/sigma_));
247 
248  if (yGuess < p)
249  {
250  success = true;
251  }
252  }
253 
254  d = dGuess;
255  tc = 0.0;
256  }
257 
258  if (shear)
259  {
260  scalar dC = weConst_*sigma/(rhoc*sqr(Urmag));
261  scalar d32Red = 4.0*(d32*dC)/(5.0*dC - d32);
262 
263  scalar d05 = d32Coeff_ * d32Red;
264 
265  scalar x = 0.0;
266  scalar yGuess = 0.0;
267  scalar dGuess = 0.0;
268 
269  while(!success)
270  {
271  x = cDmaxS_*rndGen.sample01<scalar>();
272  dGuess = sqr(x)*d05;
273  yGuess = rndGen.sample01<scalar>();
274 
275  scalar p =
276  x
277  /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
278  *exp(-0.5*sqr((x - mu_)/sigma_));
279 
280  if (yGuess<p)
281  {
282  success = true;
283  }
284  }
285 
286  d = dC;
287  dChild = dGuess;
288  massChild = corePerc_*mass0;
289  mass -= massChild;
290 
291  addChild = true;
292  // reset timer
293  tc = 0.0;
294  }
295 
296  // correct nParticle to conserve mass
297  nParticle = mass/(rhopi6*pow3(d));
298  }
299 
300  return addChild;
301 }
302 
303 
304 // ************************************************************************* //
bool success
scalar y
Templated break-up model class.
Definition: BreakupModel.H:55
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:80
Secondary Breakup Model to take account of the different breakup regimes, bag, solutionmode,...
Definition: SHF.H:57
SHF(const dictionary &, CloudType &)
Construct from dictionary.
Definition: SHF.C:32
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, const label injectori, scalar &dChild, scalar &massChild)
Update the parcel properties.
Definition: SHF.C:117
virtual ~SHF()
Destructor.
Definition: SHF.C:109
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
Random number generator.
Type sample01()
Return a type with components uniformly distributed between.
U
Definition: pEqn.H:72
const scalar twoPi(2 *pi)
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
const dimensionedScalar mu
Atomic mass unit.
dimensionedScalar exp(const dimensionedScalar &ds)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dictionary dict
randomGenerator rndGen(653213)
volScalarField & p