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-2022 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 =
108  this->owner().injectors().averageParcelMass();
109 
110  scalar r = 0.5*d;
111  scalar d3 = pow3(d);
112  scalar d03 = pow3(d0);
113 
114  scalar rhopi6 = rho*constant::mathematical::pi/6.0;
115  scalar mass = nParticle*d3*rhopi6;
116  scalar mass0 = nParticle*d03*rhopi6;
117 
118  scalar weGas = 0.5*rhoc*sqr(Urmag)*d/sigma;
119  scalar weLiquid = 0.5*rho*sqr(Urmag)*d/sigma;
120 
121  // Note: Reitz is using radius instead of diameter for Re-number
122  scalar reLiquid = rho*Urmag*r/mu;
123  scalar ohnesorge = sqrt(weLiquid)/(reLiquid + vSmall);
124  scalar taylor = ohnesorge*sqrt(weGas);
125 
126  vector acceleration = Urel/tMom;
127  vector trajectory = U/mag(U);
128  scalar gt = (g + acceleration) & trajectory;
129 
130  // frequency of the fastest growing KH-wave
131  scalar omegaKH =
132  (0.34 + 0.38*pow(weGas, 1.5))
133  /((1.0 + ohnesorge)*(1.0 + 1.4*pow(taylor, 0.6)))
134  *sqrt(sigma/(rho*pow3(r)));
135 
136  // corresponding KH wave-length.
137  scalar lambdaKH =
138  9.02
139  *r
140  *(1.0 + 0.45*sqrt(ohnesorge))
141  *(1.0 + 0.4*pow(taylor, 0.7))
142  /pow(1.0 + 0.865*pow(weGas, 1.67), 0.6);
143 
144  // characteristic Kelvin-Helmholtz breakup time
145  scalar tauKH = 3.726*b1_*r/(omegaKH*lambdaKH);
146 
147  // stable KH diameter
148  scalar dc = 2.0*b0_*lambdaKH;
149 
150  // the frequency of the fastest growing RT wavelength.
151  scalar helpVariable = mag(gt*(rho - rhoc));
152  scalar omegaRT = sqrt
153  (
154  2.0*pow(helpVariable, 1.5)
155  /(3.0*sqrt(3.0*sigma)*(rhoc + rho))
156  );
157 
158  // RT wave number
159  scalar KRT = sqrt(helpVariable/(3.0*sigma + vSmall));
160 
161  // wavelength of the fastest growing RT frequency
162  scalar lambdaRT = constant::mathematical::twoPi*cRT_/(KRT + vSmall);
163 
164  // if lambdaRT < diameter, then RT waves are growing on the surface
165  // and we start to keep track of how long they have been growing
166  if ((tc > 0) || (lambdaRT < d) )
167  {
168  tc += dt;
169  }
170 
171  // characteristic RT breakup time
172  scalar tauRT = cTau_/(omegaRT + vSmall);
173 
174  // check if we have RT breakup
175  if ((tc > tauRT) && (lambdaRT < d))
176  {
177  // the RT breakup creates diameter/lambdaRT new droplets
178  tc = -great;
179  scalar nDrops = d/lambdaRT;
180  d = cbrt(d3/nDrops);
181  }
182  // otherwise check for KH breakup
183  else if (dc < d)
184  {
185  // no breakup below Weber = 12
186  if (weGas > weberLimit_)
187  {
188  scalar fraction = dt/tauKH;
189 
190  // reduce the diameter according to the rate-equation
191  d = (fraction*dc + d)/(1.0 + fraction);
192 
193  // scalar ms0 = rho*pow3(dc)*mathematicalConstant::pi/6.0;
194  scalar ms0 = mass0*(1.0 - pow3(d/d0));
195  ms += ms0;
196 
197  if (ms/averageParcelMass > msLimit_)
198  {
199  // Correct evaluation of the number of child droplets and the
200  // diameter of parcel droplets after breaukp
201  // Solution of cubic equation for the diameter of the parent
202  // drops after breakup, see Eq. 18 in
203  // Patterson & Reitz, SAE 980131
204  bool br3 = true;
205  scalar ae3 = 1.0;
206  scalar be3 = -dc;
207  scalar ce3 = 0.0;
208  scalar de3 = d*d*(dc - d);
209  scalar qe3 =
210  pow3(be3/(3.0*ae3)) - be3*ce3/(6.0*ae3*ae3) + de3/(2.0*ae3);
211  scalar pe3 = (3.0*ae3*ce3 - be3*be3)/(9.0*ae3*ae3);
212  scalar D3 = qe3*qe3 + pe3*pe3*pe3;
213 
214  if (D3 < 0) br3 = false;
215 
216  if (br3)
217  {
218  D3 = sqrt(D3);
219  scalar ue3 = cbrt(-qe3 + D3);
220  scalar ve3 = cbrt(-qe3 - D3);
221  scalar dParenDrops = ue3 + ve3 - be3/3.;
222  scalar mc = nParticle*(pow3(d) - pow3(dParenDrops));
223  scalar nChildDrops = mc/pow3(dc);
224 
225  if (nChildDrops >= nParticle)
226  {
227  addParcel = true;
228  d = dParenDrops;
229  ms = 0.0;
230  dChild = dc;
231  massChild = mc*rhopi6;
232 
233  // reduce the parent mass by reducing nParticle
234  mass -= massChild;
235  }
236  }
237  }
238  }
239  }
240  else if (KHindex < 0.5)
241  {
242  // Case of larger drops after breakup (Reitz, Atomisation & Spray
243  // Technology 3 (1987) 309-337, p.322) pIndKH() should be introduced
244 
245  scalar lengthScale =
246  min(lambdaKH, constant::mathematical::twoPi*Urmag/omegaKH);
247  scalar diameterLargerDrop = cbrt(1.5*d*d*lengthScale);
248  d = diameterLargerDrop;
249  ms = 0.0;
250  KHindex = 1.0;
251  }
252 
253  // correct the number of parcels in parent
254  scalar massDrop = pow3(d)*rhopi6;
255  nParticle = mass/massDrop;
256 
257  return addParcel;
258 }
259 
260 
261 // ************************************************************************* //
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:156
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 CloudType & owner() const
Return const access to the owner cloud.
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
dimensionedScalar cbrt(const dimensionedScalar &ds)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
const scalar twoPi(2 *pi)
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:55
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
virtual ~ReitzKHRT()
Destructor.
Definition: ReitzKHRT.C:73