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