LISAAtomisation.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-2021 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 "LISAAtomisation.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class CloudType>
32 (
33  const dictionary& dict,
34  CloudType& owner
35 )
36 :
37  AtomisationModel<CloudType>(dict, owner, typeName),
38  Cl_(this->coeffDict().template lookup<scalar>("Cl")),
39  cTau_(this->coeffDict().template lookup<scalar>("cTau")),
40  Q_(this->coeffDict().template lookup<scalar>("Q")),
41  lisaExp_(this->coeffDict().template lookup<scalar>("lisaExp")),
42  injectorDirection_(this->coeffDict().lookup("injectorDirection")),
43  SMDCalcMethod_(this->coeffDict().lookup("SMDCalculationMethod"))
44 {
45  // Note: Would be good if this could be picked up from the injector
46  injectorDirection_ /= mag(injectorDirection_);
47 
48  if (SMDCalcMethod_ == "method1")
49  {
50  SMDMethod_ = method1;
51  }
52  else if (SMDCalcMethod_ == "method2")
53  {
54  SMDMethod_ = method2;
55  }
56  else
57  {
58  SMDMethod_ = method2;
59  Info<< "Warning: SMDCalculationMethod " << SMDCalcMethod_
60  << " unknown. Options are (method1 | method2). Using method2"
61  << endl;
62  }
63 }
64 
65 template<class CloudType>
67 (
69 )
70 :
72  Cl_(am.Cl_),
73  cTau_(am.cTau_),
74  Q_(am.Q_),
75  lisaExp_(am.lisaExp_),
76  injectorDirection_(am.injectorDirection_),
77  SMDCalcMethod_(am.SMDCalcMethod_)
78 {}
79 
80 
81 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
82 
83 template<class CloudType>
85 {}
86 
87 
88 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
89 
90 template<class CloudType>
92 {
93  return 1.0;
94 }
95 
96 
97 template<class CloudType>
99 {
100  return true;
101 }
102 
103 
104 template<class CloudType>
106 (
107  const scalar dt,
108  scalar& d,
109  scalar& liquidCore,
110  scalar& tc,
111  const scalar rho,
112  const scalar mu,
113  const scalar sigma,
114  const scalar volFlowRate,
115  const scalar rhoAv,
116  const scalar Urel,
117  const vector& pos,
118  const vector& injectionPos,
119  const scalar pAmbient,
120  const scalar chi,
121  Random& rndGen
122 ) const
123 {
124  if (volFlowRate < small)
125  {
126  return;
127  }
128 
129  scalar tau = 0.0;
130  scalar dL = 0.0;
131  scalar k = 0.0;
132 
133  // update atomisation characteristic time
134  tc += dt;
135 
136  scalar We = 0.5*rhoAv*sqr(Urel)*d/sigma;
137  scalar nu = mu/rho;
138 
139  scalar Q = rhoAv/rho;
140 
141  vector diff = pos - injectionPos;
142  scalar pWalk = mag(diff);
143  scalar traveledTime = pWalk/Urel;
144 
145  scalar h = diff & injectorDirection_;
146  scalar delta = sqrt(sqr(pWalk) - sqr(h));
147 
148  scalar hSheet = volFlowRate/(constant::mathematical::pi*delta*Urel);
149 
150  // update drop diameter
151  d = min(d, hSheet);
152 
153  if (We > 27.0/16.0)
154  {
155  scalar kPos = 0.0;
156  scalar kNeg = Q*sqr(Urel)*rho/sigma;
157 
158  scalar derivPos = sqrt(Q*sqr(Urel));
159 
160  scalar derivNeg =
161  (
162  8.0*sqr(nu)*pow3(kNeg)
163  + Q*sqr(Urel)*kNeg
164  - 3.0*sigma/2.0/rho*sqr(kNeg)
165  )
166  / sqrt
167  (
168  4.0*sqr(nu)*pow4(kNeg)
169  + Q*sqr(Urel)*sqr(kNeg)
170  - sigma*pow3(kNeg)/rho
171  )
172  - 4.0*nu*kNeg;
173 
174  scalar kOld = 0.0;
175 
176  for (label i=0; i<40; i++)
177  {
178  k = kPos - (derivPos/((derivNeg - derivPos)/(kNeg - kPos)));
179 
180  scalar derivk =
181  (
182  8.0*sqr(nu)*pow3(k)
183  + Q*sqr(Urel)*k
184  - 3.0*sigma/2.0/rho*sqr(k)
185  )
186  / sqrt
187  (
188  4.0*sqr(nu)*pow4(k)
189  + Q*sqr(Urel)*sqr(k)
190  - sigma*pow3(k)/rho
191  )
192  - 4.0*nu*k;
193 
194  if (derivk > 0)
195  {
196  derivPos = derivk;
197  kPos = k;
198  }
199  else
200  {
201  derivNeg = derivk;
202  kNeg = k;
203  }
204 
205  if (mag(k - kOld)/k < 1e-4)
206  {
207  break;
208  }
209 
210  kOld = k;
211  }
212 
213  scalar omegaS =
214  - 2.0*nu*sqr(k)
215  + sqrt
216  (
217  4.0*sqr(nu)*pow4(k)
218  + Q*sqr(Urel)*sqr(k)
219  - sigma*pow3(k)/rho
220  );
221 
222  tau = cTau_/omegaS;
223 
224  dL = sqrt(8.0*d/k);
225  }
226  else
227  {
228  k = rhoAv*sqr(Urel)/(2.0*sigma);
229 
230  scalar J = 0.5*traveledTime*hSheet;
231 
232  tau = pow(3.0*cTau_, 2.0/3.0)*cbrt(J*sigma/(sqr(Q)*pow4(Urel)*rho));
233 
234  dL = sqrt(4.0*d/k);
235  }
236 
237  scalar kL = 1.0/(dL*sqrt(0.5 + 1.5 * mu/sqrt(rho*sigma*dL)));
238 
239  scalar dD = cbrt(3.0*constant::mathematical::pi*sqr(dL)/kL);
240 
241  scalar atmPressure = 1.0e+5;
242 
243  scalar pRatio = pAmbient/atmPressure;
244 
245  dD = dD*pow(pRatio, lisaExp_);
246 
247  scalar pExp = 0.135;
248 
249  // modifying dD to take account of flash boiling
250  dD = dD*(1.0 - chi*pow(pRatio, -pExp));
251  scalar lBU = Cl_ * mag(Urel)*tau;
252 
253  if (pWalk > lBU)
254  {
255  scalar x = 0;
256 
257  switch (SMDMethod_)
258  {
259  case method1:
260  {
261  #include "LISASMDCalcMethod1.H"
262  break;
263  }
264  case method2:
265  {
266  #include "LISASMDCalcMethod2.H"
267  break;
268  }
269  }
270 
271  // New droplet properties
272  liquidCore = 0.0;
273  d = x;
274  tc = 0.0;
275  }
276 }
277 
278 
279 // ************************************************************************* //
scalar delta
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:403
dictionary dict
LISAAtomisation(const dictionary &, CloudType &)
Construct from dictionary.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Templated atomisation model class.
Definition: SprayCloud.H:52
label k
Boltzmann constant.
const dimensionedScalar h
Planck constant.
virtual bool calcChi() const
Flag to indicate if chi needs to be calculated.
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
stressControl lookup("compactNormalStress") >> compactNormalStress
Primary Breakup Model for pressure swirl atomisers.
dimensionedScalar cbrt(const dimensionedScalar &ds)
virtual scalar initLiquidCore() const
Initial value of liquidCore.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
Urel
Definition: pEqn.H:60
Random number generator.
Definition: Random.H:57
virtual void update(const scalar dt, scalar &d, scalar &liquidCore, scalar &tc, const scalar rho, const scalar mu, const scalar sigma, const scalar volFlowRate, const scalar rhoAv, const scalar Urel, const vector &pos, const vector &injectionPos, const scalar pAmbient, const scalar chi, Random &rndGen) const
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar pow4(const dimensionedScalar &ds)
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual ~LISAAtomisation()
Destructor.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105