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-2023 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  const label injectori,
102  scalar& dChild,
103  scalar& massChild
104 )
105 {
106  bool addParcel = false;
107 
108  const scalar averageParcelMass =
109  this->owner().injectors()[injectori].averageParcelMass();
110 
111  scalar r = 0.5*d;
112  scalar d3 = pow3(d);
113  scalar d03 = pow3(d0);
114 
115  scalar rhopi6 = rho*constant::mathematical::pi/6.0;
116  scalar mass = nParticle*d3*rhopi6;
117  scalar mass0 = nParticle*d03*rhopi6;
118 
119  scalar weGas = 0.5*rhoc*sqr(Urmag)*d/sigma;
120  scalar weLiquid = 0.5*rho*sqr(Urmag)*d/sigma;
121 
122  // Note: Reitz is using radius instead of diameter for Re-number
123  scalar reLiquid = rho*Urmag*r/mu;
124  scalar ohnesorge = sqrt(weLiquid)/(reLiquid + vSmall);
125  scalar taylor = ohnesorge*sqrt(weGas);
126 
127  vector acceleration = Urel/tMom;
128  vector trajectory = U/mag(U);
129  scalar gt = (g + acceleration) & trajectory;
130 
131  // frequency of the fastest growing KH-wave
132  scalar omegaKH =
133  (0.34 + 0.38*pow(weGas, 1.5))
134  /((1.0 + ohnesorge)*(1.0 + 1.4*pow(taylor, 0.6)))
135  *sqrt(sigma/(rho*pow3(r)));
136 
137  // corresponding KH wave-length.
138  scalar lambdaKH =
139  9.02
140  *r
141  *(1.0 + 0.45*sqrt(ohnesorge))
142  *(1.0 + 0.4*pow(taylor, 0.7))
143  /pow(1.0 + 0.865*pow(weGas, 1.67), 0.6);
144 
145  // characteristic Kelvin-Helmholtz breakup time
146  scalar tauKH = 3.726*b1_*r/(omegaKH*lambdaKH);
147 
148  // stable KH diameter
149  scalar dc = 2.0*b0_*lambdaKH;
150 
151  // the frequency of the fastest growing RT wavelength.
152  scalar helpVariable = mag(gt*(rho - rhoc));
153  scalar omegaRT = sqrt
154  (
155  2.0*pow(helpVariable, 1.5)
156  /(3.0*sqrt(3.0*sigma)*(rhoc + rho))
157  );
158 
159  // RT wave number
160  scalar KRT = sqrt(helpVariable/(3.0*sigma + vSmall));
161 
162  // wavelength of the fastest growing RT frequency
163  scalar lambdaRT = constant::mathematical::twoPi*cRT_/(KRT + vSmall);
164 
165  // if lambdaRT < diameter, then RT waves are growing on the surface
166  // and we start to keep track of how long they have been growing
167  if ((tc > 0) || (lambdaRT < d) )
168  {
169  tc += dt;
170  }
171 
172  // characteristic RT breakup time
173  scalar tauRT = cTau_/(omegaRT + vSmall);
174 
175  // check if we have RT breakup
176  if ((tc > tauRT) && (lambdaRT < d))
177  {
178  // the RT breakup creates diameter/lambdaRT new droplets
179  tc = -great;
180  scalar nDrops = d/lambdaRT;
181  d = cbrt(d3/nDrops);
182  }
183  // otherwise check for KH breakup
184  else if (dc < d)
185  {
186  // no breakup below Weber = 12
187  if (weGas > weberLimit_)
188  {
189  scalar fraction = dt/tauKH;
190 
191  // reduce the diameter according to the rate-equation
192  d = (fraction*dc + d)/(1.0 + fraction);
193 
194  // calculate the stripped mass
195  ms += mass0*(1.0 - pow3(d/d0));
196 
197  if (ms/averageParcelMass > msLimit_)
198  {
199  // Correct evaluation of the number of child droplets and the
200  // diameter of parcel droplets after breakup
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 // ************************************************************************* //
scalar y
Templated break-up model class.
Definition: BreakupModel.H:55
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:80
secondary breakup model which uses the Kelvin-Helmholtz instability theory to predict the 'stripped' ...
Definition: ReitzKHRT.H:50
virtual ~ReitzKHRT()
Destructor.
Definition: ReitzKHRT.C:73
ReitzKHRT(const dictionary &, CloudType &)
Construct from dictionary.
Definition: ReitzKHRT.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 diameter.
Definition: ReitzKHRT.C:81
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:710
const dictionary & coeffDict() const
Return const access to the coefficients dictionary.
Definition: subModelBase.C:128
virtual bool defaultCoeffs(const bool printMsg) const
Returns true if defaultCoeffs is true and outputs on printMsg.
Definition: subModelBase.C:140
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.
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)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
dimensionedScalar cbrt(const dimensionedScalar &ds)
dictionary dict