FreeFOAM The Cross-Platform CFD Toolkit
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 | Copyright (C) 1991-2010 OpenCFD Ltd.
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"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
38 
40 (
41  breakupModel,
42  SHF,
43  dictionary
44 );
45 
46 
47 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48 
49 // Construct from components
51 (
52  const dictionary& dict,
53  spray& sm
54 )
55 :
56  breakupModel(dict, sm),
57  coeffsDict_(dict.subDict(typeName + "Coeffs")),
58  g_(sm.g()),
59  rndGen_(sm.rndGen()),
60  weCorrCoeff_(readScalar(coeffsDict_.lookup("weCorrCoeff"))),
61  weBuCrit_(readScalar(coeffsDict_.lookup("weBuCrit"))),
62  weBuBag_(readScalar(coeffsDict_.lookup("weBuBag"))),
63  weBuMM_(readScalar(coeffsDict_.lookup("weBuMM"))),
64  ohnCoeffCrit_(readScalar(coeffsDict_.lookup("ohnCoeffCrit"))),
65  ohnCoeffBag_(readScalar(coeffsDict_.lookup("ohnCoeffBag"))),
66  ohnCoeffMM_(readScalar(coeffsDict_.lookup("ohnCoeffMM"))),
67  ohnExpCrit_(readScalar(coeffsDict_.lookup("ohnExpCrit"))),
68  ohnExpBag_(readScalar(coeffsDict_.lookup("ohnExpBag"))),
69  ohnExpMM_(readScalar(coeffsDict_.lookup("ohnExpMM"))),
70  cInit_(readScalar(coeffsDict_.lookup("Cinit"))),
71  c1_(readScalar(coeffsDict_.lookup("C1"))),
72  c2_(readScalar(coeffsDict_.lookup("C2"))),
73  c3_(readScalar(coeffsDict_.lookup("C3"))),
74  cExp1_(readScalar(coeffsDict_.lookup("Cexp1"))),
75  cExp2_(readScalar(coeffsDict_.lookup("Cexp2"))),
76  cExp3_(readScalar(coeffsDict_.lookup("Cexp3"))),
77  weConst_(readScalar(coeffsDict_.lookup("Weconst"))),
78  weCrit1_(readScalar(coeffsDict_.lookup("Wecrit1"))),
79  weCrit2_(readScalar(coeffsDict_.lookup("Wecrit2"))),
80  coeffD_(readScalar(coeffsDict_.lookup("CoeffD"))),
81  onExpD_(readScalar(coeffsDict_.lookup("OnExpD"))),
82  weExpD_(readScalar(coeffsDict_.lookup("WeExpD"))),
83  mu_(readScalar(coeffsDict_.lookup("mu"))),
84  sigma_(readScalar(coeffsDict_.lookup("sigma"))),
85  d32Coeff_(readScalar(coeffsDict_.lookup("d32Coeff"))),
86  cDmaxBM_(readScalar(coeffsDict_.lookup("cDmaxBM"))),
87  cDmaxS_(readScalar(coeffsDict_.lookup("cDmaxS"))),
88  corePerc_(readScalar(coeffsDict_.lookup("corePerc")))
89 {}
90 
91 
92 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
93 
95 {}
96 
97 
98 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
99 
101 (
102  parcel& p,
103  const scalar deltaT,
104  const vector& vel,
105  const liquidMixture& fuels
106 ) const
107 {
108 
109  label celli = p.cell();
110  scalar T = p.T();
111  scalar pc = spray_.p()[celli];
112 
113  scalar sigma = fuels.sigma(pc, T, p.X());
114  scalar rhoLiquid = fuels.rho(pc, T, p.X());
115  scalar muLiquid = fuels.mu(pc, T, p.X());
116  scalar rhoGas = spray_.rho()[celli];
117 
118  scalar weGas = p.We(vel, rhoGas, sigma);
119  scalar weLiquid = p.We(vel, rhoLiquid, sigma);
120 
121  // correct the Reynolds number. Reitz is using radius instead of diameter
122 
123  scalar reLiquid = p.Re(rhoLiquid, vel, muLiquid);
124  scalar ohnesorge = sqrt(weLiquid)/(reLiquid + VSMALL);
125 
126  vector acceleration = p.Urel(vel)/p.tMom();
127  vector trajectory = p.U()/mag(p.U());
128 
129  vector vRel = p.Urel(vel);
130 
131  scalar weGasCorr = weGas/(1.0 + weCorrCoeff_ * ohnesorge);
132 
133  // droplet deformation characteristic time
134 
135  scalar tChar = p.d()/mag(vRel)*sqrt(rhoLiquid/rhoGas);
136 
137  scalar tFirst = cInit_ * tChar;
138 
139  scalar tSecond = 0;
140  scalar tCharSecond = 0;
141 
142  bool bag = false;
143  bool multimode = false;
144  bool shear = false;
145  bool success = false;
146 
147 
148  // updating the droplet characteristic time
149  p.ct() += deltaT;
150 
151  if(weGas > weConst_)
152  {
153  if(weGas < weCrit1_)
154  {
155  tCharSecond = c1_*pow((weGas - weConst_),cExp1_);
156  }
157  else if(weGas >= weCrit1_ && weGas <= weCrit2_)
158  {
159  tCharSecond = c2_*pow((weGas - weConst_),cExp2_);
160  }
161  else
162  {
163  tCharSecond = c3_*pow((weGas - weConst_),cExp3_);
164  }
165  }
166 
167  scalar weC = weBuCrit_*(1.0+ohnCoeffCrit_*pow(ohnesorge,ohnExpCrit_));
168  scalar weB = weBuBag_*(1.0+ohnCoeffBag_*pow(ohnesorge, ohnExpBag_));
169  scalar weMM = weBuMM_*(1.0+ohnCoeffMM_*pow(ohnesorge, ohnExpMM_));
170 
171  if(weGas > weC && weGas < weB)
172  {
173  bag = true;
174  }
175 
176  if(weGas >= weB && weGas <= weMM)
177  {
178  multimode = true;
179  }
180 
181  if(weGas > weMM)
182  {
183  shear = true;
184  }
185 
186  tSecond = tCharSecond * tChar;
187 
188  scalar tBreakUP = tFirst + tSecond;
189  if(p.ct() > tBreakUP)
190  {
191 
192  scalar d32 = coeffD_*p.d()*pow(ohnesorge,onExpD_)*pow(weGasCorr,weExpD_);
193 
194  if(bag || multimode)
195  {
196 
197  scalar d05 = d32Coeff_ * d32;
198 
199  scalar x = 0.0;
200  scalar y = 0.0;
201  scalar d = 0.0;
202 
203  while(!success)
204  {
205  x = cDmaxBM_*rndGen_.scalar01();
206  d = sqr(x)*d05;
207  y = rndGen_.scalar01();
208 
209  scalar p = x/(2.0*sqrt(2.0*mathematicalConstant::pi)*sigma_)*exp(-0.5*sqr((x-mu_)/sigma_));
210 
211  if (y<p)
212  {
213  success = true;
214  }
215  }
216 
217  p.d() = d;
218  p.ct() = 0.0;
219  }
220 
221  if(shear)
222  {
223  scalar dC = weConst_*sigma/(rhoGas*sqr(mag(vRel)));
224  scalar d32Red = 4.0*(d32 * dC)/(5.0 * dC - d32);
225  scalar initMass = p.m();
226 
227  scalar d05 = d32Coeff_ * d32Red;
228 
229  scalar x = 0.0;
230  scalar y = 0.0;
231  scalar d = 0.0;
232 
233  while(!success)
234  {
235 
236  x = cDmaxS_*rndGen_.scalar01();
237  d = sqr(x)*d05;
238  y = rndGen_.scalar01();
239 
240  scalar p = x/(2.0*sqrt(2.0*mathematicalConstant::pi)*sigma_)*exp(-0.5*sqr((x-mu_)/sigma_));
241 
242  if (y<p)
243  {
244  success = true;
245  }
246  }
247 
248  p.d() = dC;
249  p.m() = corePerc_ * initMass;
250 
251  spray_.addParticle
252  (
253  new parcel
254  (
255  spray_,
256  p.position(),
257  p.cell(),
258  p.n(),
259  d,
260  p.T(),
261  (1.0 - corePerc_)* initMass,
262  0.0,
263  0.0,
264  0.0,
265  -GREAT,
266  p.tTurb(),
267  0.0,
268  scalar(p.injector()),
269  p.U(),
270  p.Uturb(),
271  p.X(),
272  p.fuelNames()
273  )
274  );
275 
276  p.ct() = 0.0;
277  }
278  }
279 }
280 
281 
282 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
283 
284 } // End namespace Foam
285 
286 // ************************ vim: set sw=4 sts=4 et: ************************ //