SHF.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2020 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "SHF.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class CloudType>
34 (
35  const dictionary& dict,
36  CloudType& owner
37 )
38 :
39  BreakupModel<CloudType>(dict, owner, typeName),
40  weCorrCoeff_(this->coeffDict().getScalar("weCorrCoeff")),
41  weBuCrit_(this->coeffDict().getScalar("weBuCrit")),
42  weBuBag_(this->coeffDict().getScalar("weBuBag")),
43  weBuMM_(this->coeffDict().getScalar("weBuMM")),
44  ohnCoeffCrit_(this->coeffDict().getScalar("ohnCoeffCrit")),
45  ohnCoeffBag_(this->coeffDict().getScalar("ohnCoeffBag")),
46  ohnCoeffMM_(this->coeffDict().getScalar("ohnCoeffMM")),
47  ohnExpCrit_(this->coeffDict().getScalar("ohnExpCrit")),
48  ohnExpBag_(this->coeffDict().getScalar("ohnExpBag")),
49  ohnExpMM_(this->coeffDict().getScalar("ohnExpMM")),
50  cInit_(this->coeffDict().getScalar("Cinit")),
51  c1_(this->coeffDict().getScalar("C1")),
52  c2_(this->coeffDict().getScalar("C2")),
53  c3_(this->coeffDict().getScalar("C3")),
54  cExp1_(this->coeffDict().getScalar("Cexp1")),
55  cExp2_(this->coeffDict().getScalar("Cexp2")),
56  cExp3_(this->coeffDict().getScalar("Cexp3")),
57  weConst_(this->coeffDict().getScalar("Weconst")),
58  weCrit1_(this->coeffDict().getScalar("Wecrit1")),
59  weCrit2_(this->coeffDict().getScalar("Wecrit2")),
60  coeffD_(this->coeffDict().getScalar("CoeffD")),
61  onExpD_(this->coeffDict().getScalar("OnExpD")),
62  weExpD_(this->coeffDict().getScalar("WeExpD")),
63  mu_(this->coeffDict().getScalar("mu")),
64  sigma_(this->coeffDict().getScalar("sigma")),
65  d32Coeff_(this->coeffDict().getScalar("d32Coeff")),
66  cDmaxBM_(this->coeffDict().getScalar("cDmaxBM")),
67  cDmaxS_(this->coeffDict().getScalar("cDmaxS")),
68  corePerc_(this->coeffDict().getScalar("corePerc"))
69 {}
70 
71 
72 template<class CloudType>
74 :
75  BreakupModel<CloudType>(bum),
76  weCorrCoeff_(bum.weCorrCoeff_),
77  weBuCrit_(bum.weBuCrit_),
78  weBuBag_(bum.weBuBag_),
79  weBuMM_(bum.weBuMM_),
80  ohnCoeffCrit_(bum.ohnCoeffCrit_),
81  ohnCoeffBag_(bum.ohnCoeffBag_),
82  ohnCoeffMM_(bum.ohnCoeffMM_),
83  ohnExpCrit_(bum.ohnExpCrit_),
84  ohnExpBag_(bum.ohnExpBag_),
85  ohnExpMM_(bum.ohnExpMM_),
86  cInit_(bum.cInit_),
87  c1_(bum.c1_),
88  c2_(bum.c2_),
89  c3_(bum.c3_),
90  cExp1_(bum.cExp1_),
91  cExp2_(bum.cExp2_),
92  cExp3_(bum.cExp3_),
93  weConst_(bum.weConst_),
94  weCrit1_(bum.weCrit1_),
95  weCrit2_(bum.weCrit2_),
96  coeffD_(bum.coeffD_),
97  onExpD_(bum.onExpD_),
98  weExpD_(bum.weExpD_),
99  mu_(bum.mu_),
100  sigma_(bum.sigma_),
101  d32Coeff_(bum.d32Coeff_),
102  cDmaxBM_(bum.cDmaxBM_),
103  cDmaxS_(bum.cDmaxS_),
104  corePerc_(bum.corePerc_)
105 {}
106 
107 
108 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
109 
110 template<class CloudType>
112 {}
113 
114 
115 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
116 
117 template<class CloudType>
119 (
120  const scalar dt,
121  const vector& g,
122  scalar& d,
123  scalar& tc,
124  scalar& ms,
125  scalar& nParticle,
126  scalar& KHindex,
127  scalar& y,
128  scalar& yDot,
129  const scalar d0,
130  const scalar rho,
131  const scalar mu,
132  const scalar sigma,
133  const vector& U,
134  const scalar rhoc,
135  const scalar muc,
136  const vector& Urel,
137  const scalar Urmag,
138  const scalar tMom,
139  scalar& dChild,
140  scalar& massChild
141 )
142 {
143  Random& rndGen = this->owner().rndGen();
144 
145  bool addChild = false;
146 
147  scalar d03 = pow3(d);
148  scalar rhopi6 = rho*constant::mathematical::pi/6.0;
149  scalar mass0 = nParticle*rhopi6*d03;
150  scalar mass = mass0;
151 
152  scalar weGas = 0.5*rhoc*sqr(Urmag)*d/sigma;
153  scalar weLiquid = 0.5*rho*sqr(Urmag)*d/sigma;
154 
155  // correct the Reynolds number. Reitz is using radius instead of diameter
156  scalar reLiquid = 0.5*Urmag*d/mu;
157  scalar ohnesorge = sqrt(weLiquid)/(reLiquid + VSMALL);
158 
159  scalar weGasCorr = weGas/(1.0 + weCorrCoeff_*ohnesorge);
160 
161  // update the droplet characteristic time
162  tc += dt;
163 
164  // droplet deformation characteristic rate
165  scalar rChar = Urmag/d*sqrt(rhoc/rho);
166 
167  // return if the characteristic deformation rate is too low for the
168  // following modelling to be calculable
169  if (tc*rChar < SMALL)
170  {
171  return false;
172  }
173 
174  // droplet deformation characteristic time
175  scalar tChar = 1/rChar;
176 
177  scalar tFirst = cInit_*tChar;
178 
179  scalar tSecond = 0;
180  scalar tCharSecond = 0;
181 
182  bool bag = false;
183  bool multimode = false;
184  bool shear = false;
185  bool success = false;
186 
187 
188  if (weGas > weConst_)
189  {
190  if (weGas < weCrit1_)
191  {
192  tCharSecond = c1_*pow((weGas - weConst_), cExp1_);
193  }
194  else if (weGas >= weCrit1_ && weGas <= weCrit2_)
195  {
196  tCharSecond = c2_*pow((weGas - weConst_), cExp2_);
197  }
198  else
199  {
200  tCharSecond = c3_*pow((weGas - weConst_), cExp3_);
201  }
202  }
203 
204  scalar weC = weBuCrit_*(1.0 + ohnCoeffCrit_*pow(ohnesorge, ohnExpCrit_));
205  scalar weB = weBuBag_*(1.0 + ohnCoeffBag_*pow(ohnesorge, ohnExpBag_));
206  scalar weMM = weBuMM_*(1.0 + ohnCoeffMM_*pow(ohnesorge, ohnExpMM_));
207 
208  if (weGas > weC && weGas < weB)
209  {
210  bag = true;
211  }
212 
213  if (weGas >= weB && weGas <= weMM)
214  {
215  multimode = true;
216  }
217 
218  if (weGas > weMM)
219  {
220  shear = true;
221  }
222 
223  tSecond = tCharSecond*tChar;
224 
225  scalar tBreakUP = tFirst + tSecond;
226  if (tc > tBreakUP)
227  {
228  scalar d32 = coeffD_*d*pow(ohnesorge, onExpD_)*pow(weGasCorr, weExpD_);
229 
230  if (bag || multimode)
231  {
232  scalar d05 = d32Coeff_ * d32;
233 
234  scalar x = 0.0;
235  scalar yGuess = 0.0;
236  scalar dGuess = 0.0;
237 
238  while(!success)
239  {
240  x = cDmaxBM_*rndGen.sample01<scalar>();
241  dGuess = sqr(x)*d05;
242  yGuess = rndGen.sample01<scalar>();
243 
244  scalar p =
245  x
246  /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
247  *exp(-0.5*sqr((x - mu_)/sigma_));
248 
249  if (yGuess < p)
250  {
251  success = true;
252  }
253  }
254 
255  d = dGuess;
256  tc = 0.0;
257  }
258 
259  if (shear)
260  {
261  scalar dC = weConst_*sigma/(rhoc*sqr(Urmag));
262  scalar d32Red = 4.0*(d32*dC)/(5.0*dC - d32);
263 
264  scalar d05 = d32Coeff_ * d32Red;
265 
266  scalar x = 0.0;
267  scalar yGuess = 0.0;
268  scalar dGuess = 0.0;
269 
270  while(!success)
271  {
272  x = cDmaxS_*rndGen.sample01<scalar>();
273  dGuess = sqr(x)*d05;
274  yGuess = rndGen.sample01<scalar>();
275 
276  scalar p =
277  x
278  /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
279  *exp(-0.5*sqr((x - mu_)/sigma_));
280 
281  if (yGuess<p)
282  {
283  success = true;
284  }
285  }
286 
287  d = dC;
288  dChild = dGuess;
289  massChild = corePerc_*mass0;
290  mass -= massChild;
291 
292  addChild = true;
293  // reset timer
294  tc = 0.0;
295  }
296 
297  // correct nParticle to conserve mass
298  nParticle = mass/(rhopi6*pow3(d));
299  }
300 
301  return addChild;
302 }
303 
304 
305 // ************************************************************************* //
Foam::Random
Random number generator.
Definition: Random.H:59
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::SHF::~SHF
virtual ~SHF()
Destructor.
Definition: SHF.C:111
Foam::constant::physicoChemical::mu
const dimensionedScalar mu
Atomic mass unit.
Definition: createFieldRefs.H:4
Urel
Urel
Definition: pEqn.H:56
success
bool success
Definition: LISASMDCalcMethod1.H:16
Foam::Random::sample01
Type sample01()
Return a sample whose components lie in the range [0,1].
Definition: RandomTemplates.C:36
Foam::exp
dimensionedScalar exp(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:261
Foam::SHF::update
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:119
rho
rho
Definition: readInitialConditions.H:88
Foam::SHF::SHF
SHF(const dictionary &, CloudType &)
Construct from dictionary.
Definition: SHF.C:34
Foam::DSMCCloud::rndGen
Random & rndGen()
Return reference to the random object.
Definition: DSMCCloudI.H:124
SHF.H
Foam::pow3
dimensionedScalar pow3(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:89
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
Foam::DSMCCloud
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:71
Foam::constant::mathematical::twoPi
constexpr scalar twoPi(2 *M_PI)
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
g
const uniformDimensionedVectorField & g
Definition: createFluidFields.H:26
Foam::SHF
Secondary Breakup Model to take account of the different breakup regimes, bag, molutimode,...
Definition: SHF.H:60
Foam::BreakupModel
Templated break-up model class.
Definition: SprayCloud.H:50
U
U
Definition: pEqn.H:72
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::Vector< scalar >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
x
x
Definition: LISASMDCalcMethod2.H:52
rndGen
Random rndGen
Definition: createFields.H:23
sigma
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
y
scalar y
Definition: LISASMDCalcMethod1.H:14