TAB.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 "TAB.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, true),
38  SMDCalcMethod_(this->coeffDict().lookup("SMDCalculationMethod"))
39 {
40  // calculate the inverse function of the Rossin-Rammler Distribution
41  const scalar xx0 = 12.0;
42  const scalar rrd100 =
43  1.0/(1.0 - exp(-xx0)*(1.0 + xx0 + sqr(xx0)/2.0 + pow3(xx0)/6.0));
44 
45  forAll(rrd_, n)
46  {
47  scalar xx = 0.12*(n + 1);
48  rrd_[n] =
49  (1.0 - exp(-xx)*(1.0 + xx + sqr(xx)/2.0 + pow3(xx)/6.0))*rrd100;
50  }
51 
52  if (SMDCalcMethod_ == "method1")
53  {
54  SMDMethod_ = method1;
55  }
56  else if (SMDCalcMethod_ == "method2")
57  {
58  SMDMethod_ = method2;
59  }
60  else
61  {
62  SMDMethod_ = method2;
64  << "Unknown SMDCalculationMethod. Valid options are "
65  << "(method1 | method2). Using method2" << endl;
66  }
67 }
68 
69 
70 template<class CloudType>
72 :
73  BreakupModel<CloudType>(bum),
74  SMDCalcMethod_(bum.SMDCalcMethod_)
75 {}
76 
77 
78 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
79 
80 template<class CloudType>
82 {}
83 
84 
85 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
86 
87 template<class CloudType>
89 (
90  const scalar dt,
91  const vector& g,
92  scalar& d,
93  scalar& tc,
94  scalar& ms,
95  scalar& nParticle,
96  scalar& KHindex,
97  scalar& y,
98  scalar& yDot,
99  const scalar d0,
100  const scalar rho,
101  const scalar mu,
102  const scalar sigma,
103  const vector& U,
104  const scalar rhoc,
105  const scalar muc,
106  const vector& Urel,
107  const scalar Urmag,
108  const scalar tMom,
109  const label injectori,
110  scalar& dChild,
111  scalar& massChild
112 )
113 {
114  randomGenerator& rndGen = this->owner().rndGen();
115 
116  scalar r = 0.5*d;
117  scalar r2 = r*r;
118  scalar r3 = r*r2;
119 
120  scalar semiMass = nParticle*pow3(d);
121 
122  // inverse of characteristic viscous damping time
123  scalar rtd = 0.5*this->TABCmu_*mu/(rho*r2);
124 
125  // oscillation frequency (squared)
126  scalar omega2 = this->TABComega_*sigma/(rho*r3) - rtd*rtd;
127 
128  if (omega2 > 0)
129  {
130  scalar omega = sqrt(omega2);
131  scalar We = rhoc*sqr(Urmag)*r/sigma;
132  scalar Wetmp = We/this->TABtwoWeCrit_;
133 
134  scalar y1 = y - Wetmp;
135  scalar y2 = yDot/omega;
136 
137  scalar a = sqrt(y1*y1 + y2*y2);
138 
139  // scotty we may have break-up
140  if (a+Wetmp > 1.0)
141  {
142  scalar phic = y1/a;
143 
144  // constrain phic within -1 to 1
145  phic = max(min(phic, 1), -1);
146 
147  scalar phit = acos(phic);
148  scalar phi = phit;
149  scalar quad = -y2/a;
150  if (quad < 0)
151  {
152  phi = constant::mathematical::twoPi - phit;
153  }
154 
155  scalar tb = 0;
156 
157  if (mag(y) < 1.0)
158  {
159  scalar coste = 1.0;
160  if ((Wetmp - a < -1) && (yDot < 0))
161  {
162  coste = -1.0;
163  }
164 
165  scalar theta = acos((coste-Wetmp)/a);
166 
167  if (theta < phi)
168  {
169  if (constant::mathematical::twoPi - theta >= phi)
170  {
171  theta = -theta;
172  }
174  }
175  tb = (theta-phi)/omega;
176 
177  // breakup occurs
178  if (dt > tb)
179  {
180  y = 1.0;
181  yDot = -a*omega*sin(omega*tb + phi);
182  }
183 
184  }
185 
186  // update droplet size
187  if (dt > tb)
188  {
189  scalar rs =
190  r/(1.0 + (4.0/3.0)*sqr(y) + rho*r3/(8*sigma)*sqr(yDot));
191 
192  label n = 0;
193  scalar rNew = 0.0;
194  switch (SMDMethod_)
195  {
196  case method1:
197  {
198  #include "TABSMDCalcMethod1.H"
199  break;
200  }
201  case method2:
202  {
203  #include "TABSMDCalcMethod2.H"
204  break;
205  }
206  }
207 
208  if (rNew < r)
209  {
210  d = 2*rNew;
211  y = 0;
212  yDot = 0;
213  }
214  }
215  }
216  }
217  else
218  {
219  // reset droplet distortion parameters
220  y = 0;
221  yDot = 0;
222  }
223 
224  // update the nParticle count to conserve mass
225  nParticle = semiMass/pow3(d);
226 
227  // Do not add child parcel
228  return false;
229 }
230 
231 
232 // ************************************************************************* //
scalar y
label n
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Templated break-up model class.
Definition: BreakupModel.H:55
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:80
The TAB Method for Numerical Calculation of Spray Droplet Breakup.
Definition: TAB.H:63
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: TAB.C:89
virtual ~TAB()
Destructor.
Definition: TAB.C:81
TAB(const dictionary &dict, CloudType &owner)
Construct from dictionary.
Definition: TAB.C:32
@ method1
Definition: TAB.H:69
@ method2
Definition: TAB.H:70
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
Convenience class to handle the input of constant rotational speed. Reads an omega entry with default...
Definition: omega.H:54
Random number generator.
const scalar omega
U
Definition: pEqn.H:72
#define WarningInFunction
Report a warning using Foam::Warning.
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.
dimensionedScalar exp(const dimensionedScalar &ds)
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)
dimensionedScalar sin(const dimensionedScalar &ds)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar y1(const dimensionedScalar &ds)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar acos(const dimensionedScalar &ds)
dictionary dict
randomGenerator rndGen(653213)