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-2020 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  scalar& dChild,
110  scalar& massChild
111 )
112 {
113  Random& rndGen = this->owner().rndGen();
114 
115  scalar r = 0.5*d;
116  scalar r2 = r*r;
117  scalar r3 = r*r2;
118 
119  scalar semiMass = nParticle*pow3(d);
120 
121  // inverse of characteristic viscous damping time
122  scalar rtd = 0.5*this->TABCmu_*mu/(rho*r2);
123 
124  // oscillation frequency (squared)
125  scalar omega2 = this->TABComega_*sigma/(rho*r3) - rtd*rtd;
126 
127  if (omega2 > 0)
128  {
129  scalar omega = sqrt(omega2);
130  scalar We = rhoc*sqr(Urmag)*r/sigma;
131  scalar Wetmp = We/this->TABtwoWeCrit_;
132 
133  scalar y1 = y - Wetmp;
134  scalar y2 = yDot/omega;
135 
136  scalar a = sqrt(y1*y1 + y2*y2);
137 
138  // scotty we may have break-up
139  if (a+Wetmp > 1.0)
140  {
141  scalar phic = y1/a;
142 
143  // constrain phic within -1 to 1
144  phic = max(min(phic, 1), -1);
145 
146  scalar phit = acos(phic);
147  scalar phi = phit;
148  scalar quad = -y2/a;
149  if (quad < 0)
150  {
151  phi = constant::mathematical::twoPi - phit;
152  }
153 
154  scalar tb = 0;
155 
156  if (mag(y) < 1.0)
157  {
158  scalar coste = 1.0;
159  if ((Wetmp - a < -1) && (yDot < 0))
160  {
161  coste = -1.0;
162  }
163 
164  scalar theta = acos((coste-Wetmp)/a);
165 
166  if (theta < phi)
167  {
168  if (constant::mathematical::twoPi - theta >= phi)
169  {
170  theta = -theta;
171  }
173  }
174  tb = (theta-phi)/omega;
175 
176  // breakup occurs
177  if (dt > tb)
178  {
179  y = 1.0;
180  yDot = -a*omega*sin(omega*tb + phi);
181  }
182 
183  }
184 
185  // update droplet size
186  if (dt > tb)
187  {
188  scalar rs =
189  r/(1.0 + (4.0/3.0)*sqr(y) + rho*r3/(8*sigma)*sqr(yDot));
190 
191  label n = 0;
192  scalar rNew = 0.0;
193  switch (SMDMethod_)
194  {
195  case method1:
196  {
197  #include "TABSMDCalcMethod1.H"
198  break;
199  }
200  case method2:
201  {
202  #include "TABSMDCalcMethod2.H"
203  break;
204  }
205  }
206 
207  if (rNew < r)
208  {
209  d = 2*rNew;
210  y = 0;
211  yDot = 0;
212  }
213  }
214  }
215  }
216  else
217  {
218  // reset droplet distortion parameters
219  y = 0;
220  yDot = 0;
221  }
222 
223  // update the nParticle count to conserve mass
224  nParticle = semiMass/pow3(d);
225 
226  // Do not add child parcel
227  return false;
228 }
229 
230 
231 // ************************************************************************* //
dictionary dict
dimensionedScalar acos(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const CloudType & owner() const
Return const access to the owner cloud.
Random rndGen(label(0))
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
dimensionedScalar exp(const dimensionedScalar &ds)
phic
Definition: correctPhic.H:2
Random & rndGen()
Return references to the random object.
Definition: DSMCCloudI.H:121
phi
Definition: correctPhi.H:3
virtual ~TAB()
Destructor.
Definition: TAB.C:81
dimensionedScalar y1(const dimensionedScalar &ds)
Random number generator.
Definition: Random.H:57
const scalar twoPi(2 *pi)
dimensionedScalar sin(const dimensionedScalar &ds)
The TAB Method for Numerical Calculation of Spray Droplet Breakup.
Definition: TAB.H:60
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
TAB(const dictionary &dict, CloudType &owner)
Construct from dictionary.
Definition: TAB.C:32
dimensionedScalar pow3(const dimensionedScalar &ds)
#define WarningInFunction
Report a warning using Foam::Warning.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Templated break-up model class.
Definition: SprayCloud.H:55
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
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, scalar &dChild, scalar &massChild)
Update the parcel diameter.
Definition: TAB.C:89