ReitzKHRT.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 "ReitzKHRT.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  b0_(0.61),
39  b1_(40.0),
40  cTau_(1.0),
41  cRT_(0.1),
42  msLimit_(0.03),
43  weberLimit_(6.0)
44 {
45  if (!this->defaultCoeffs(true))
46  {
47  this->coeffDict().lookup("B0") >> b0_;
48  this->coeffDict().lookup("B1") >> b1_;
49  this->coeffDict().lookup("Ctau") >> cTau_;
50  this->coeffDict().lookup("CRT") >> cRT_;
51  this->coeffDict().lookup("msLimit") >> msLimit_;
52  this->coeffDict().lookup("WeberLimit") >> weberLimit_;
53  }
54 }
55 
56 
57 template<class CloudType>
59 :
60  BreakupModel<CloudType>(bum),
61  b0_(bum.b0_),
62  b1_(bum.b1_),
63  cTau_(bum.cTau_),
64  cRT_(bum.cRT_),
65  msLimit_(bum.msLimit_),
66  weberLimit_(bum.weberLimit_)
67 {}
68 
69 
70 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
71 
72 template<class CloudType>
74 {}
75 
76 
77 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
78 
79 template<class CloudType>
81 (
82  const scalar dt,
83  const vector& g,
84  scalar& d,
85  scalar& tc,
86  scalar& ms,
87  scalar& nParticle,
88  scalar& KHindex,
89  scalar& y,
90  scalar& yDot,
91  const scalar d0,
92  const scalar rho,
93  const scalar mu,
94  const scalar sigma,
95  const vector& U,
96  const scalar rhoc,
97  const scalar muc,
98  const vector& Urel,
99  const scalar Urmag,
100  const scalar tMom,
101  scalar& dChild,
102  scalar& massChild
103 )
104 {
105  bool addParcel = false;
106 
107  const scalar averageParcelMass = this->owner().averageParcelMass();
108 
109  scalar r = 0.5*d;
110  scalar d3 = pow3(d);
111  scalar d03 = pow3(d0);
112 
113  scalar rhopi6 = rho*constant::mathematical::pi/6.0;
114  scalar mass = nParticle*d3*rhopi6;
115  scalar mass0 = nParticle*d03*rhopi6;
116 
117  scalar weGas = 0.5*rhoc*sqr(Urmag)*d/sigma;
118  scalar weLiquid = 0.5*rho*sqr(Urmag)*d/sigma;
119 
120  // Note: Reitz is using radius instead of diameter for Re-number
121  scalar reLiquid = rho*Urmag*r/mu;
122  scalar ohnesorge = sqrt(weLiquid)/(reLiquid + vSmall);
123  scalar taylor = ohnesorge*sqrt(weGas);
124 
125  vector acceleration = Urel/tMom;
126  vector trajectory = U/mag(U);
127  scalar gt = (g + acceleration) & trajectory;
128 
129  // frequency of the fastest growing KH-wave
130  scalar omegaKH =
131  (0.34 + 0.38*pow(weGas, 1.5))
132  /((1.0 + ohnesorge)*(1.0 + 1.4*pow(taylor, 0.6)))
133  *sqrt(sigma/(rho*pow3(r)));
134 
135  // corresponding KH wave-length.
136  scalar lambdaKH =
137  9.02
138  *r
139  *(1.0 + 0.45*sqrt(ohnesorge))
140  *(1.0 + 0.4*pow(taylor, 0.7))
141  /pow(1.0 + 0.865*pow(weGas, 1.67), 0.6);
142 
143  // characteristic Kelvin-Helmholtz breakup time
144  scalar tauKH = 3.726*b1_*r/(omegaKH*lambdaKH);
145 
146  // stable KH diameter
147  scalar dc = 2.0*b0_*lambdaKH;
148 
149  // the frequency of the fastest growing RT wavelength.
150  scalar helpVariable = mag(gt*(rho - rhoc));
151  scalar omegaRT = sqrt
152  (
153  2.0*pow(helpVariable, 1.5)
154  /(3.0*sqrt(3.0*sigma)*(rhoc + rho))
155  );
156 
157  // RT wave number
158  scalar KRT = sqrt(helpVariable/(3.0*sigma + vSmall));
159 
160  // wavelength of the fastest growing RT frequency
161  scalar lambdaRT = constant::mathematical::twoPi*cRT_/(KRT + vSmall);
162 
163  // if lambdaRT < diameter, then RT waves are growing on the surface
164  // and we start to keep track of how long they have been growing
165  if ((tc > 0) || (lambdaRT < d) )
166  {
167  tc += dt;
168  }
169 
170  // characteristic RT breakup time
171  scalar tauRT = cTau_/(omegaRT + vSmall);
172 
173  // check if we have RT breakup
174  if ((tc > tauRT) && (lambdaRT < d))
175  {
176  // the RT breakup creates diameter/lambdaRT new droplets
177  tc = -great;
178  scalar nDrops = d/lambdaRT;
179  d = cbrt(d3/nDrops);
180  }
181  // otherwise check for KH breakup
182  else if (dc < d)
183  {
184  // no breakup below Weber = 12
185  if (weGas > weberLimit_)
186  {
187  scalar fraction = dt/tauKH;
188 
189  // reduce the diameter according to the rate-equation
190  d = (fraction*dc + d)/(1.0 + fraction);
191 
192  // scalar ms0 = rho*pow3(dc)*mathematicalConstant::pi/6.0;
193  scalar ms0 = mass0*(1.0 - pow3(d/d0));
194  ms += ms0;
195 
196  if (ms/averageParcelMass > msLimit_)
197  {
198  // Correct evaluation of the number of child droplets and the
199  // diameter of parcel droplets after breaukp
200  // Solution of cubic equation for the diameter of the parent
201  // drops after breakup, see Eq. 18 in
202  // Patterson & Reitz, SAE 980131
203  bool br3 = true;
204  scalar ae3 = 1.0;
205  scalar be3 = -dc;
206  scalar ce3 = 0.0;
207  scalar de3 = d*d*(dc - d);
208  scalar qe3 =
209  pow3(be3/(3.0*ae3)) - be3*ce3/(6.0*ae3*ae3) + de3/(2.0*ae3);
210  scalar pe3 = (3.0*ae3*ce3 - be3*be3)/(9.0*ae3*ae3);
211  scalar D3 = qe3*qe3 + pe3*pe3*pe3;
212 
213  if (D3 < 0) br3 = false;
214 
215  if (br3)
216  {
217  D3 = sqrt(D3);
218  scalar ue3 = cbrt(-qe3 + D3);
219  scalar ve3 = cbrt(-qe3 - D3);
220  scalar dParenDrops = ue3 + ve3 - be3/3.;
221  scalar mc = nParticle*(pow3(d) - pow3(dParenDrops));
222  scalar nChildDrops = mc/pow3(dc);
223 
224  if (nChildDrops >= nParticle)
225  {
226  addParcel = true;
227  d = dParenDrops;
228  ms = 0.0;
229  dChild = dc;
230  massChild = mc*rhopi6;
231 
232  // reduce the parent mass by reducing nParticle
233  mass -= massChild;
234  }
235  }
236  }
237  }
238  }
239  else if (KHindex < 0.5)
240  {
241  // Case of larger drops after breakup (Reitz, Atomization & Spray
242  // Technology 3 (1987) 309-337, p.322) pIndKH() should be introduced
243 
244  scalar lengthScale =
245  min(lambdaKH, constant::mathematical::twoPi*Urmag/omegaKH);
246  scalar diameterLargerDrop = cbrt(1.5*d*d*lengthScale);
247  d = diameterLargerDrop;
248  ms = 0.0;
249  KHindex = 1.0;
250  }
251 
252  // correct the number of parcels in parent
253  scalar massDrop = pow3(d)*rhopi6;
254  nParticle = mass/massDrop;
255 
256  return addParcel;
257 }
258 
259 
260 // ************************************************************************* //
dictionary dict
secondary breakup model which uses the Kelvin-Helmholtz instability theory to predict the &#39;stripped&#39; ...
Definition: ReitzKHRT.H:47
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
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 diameter.
Definition: ReitzKHRT.C:81
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
const CloudType & owner() const
Return const access to the owner cloud.
dimensionedScalar cbrt(const dimensionedScalar &ds)
const scalar twoPi(2 *pi)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const dimensionedScalar mu
Atomic mass unit.
ReitzKHRT(const dictionary &, CloudType &)
Construct from dictionary.
Definition: ReitzKHRT.C:32
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
Templated break-up model class.
Definition: SprayCloud.H:50
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69
virtual ~ReitzKHRT()
Destructor.
Definition: ReitzKHRT.C:73