ETAB.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) 2011-2014 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 "ETAB.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  k1_(0.2),
39  k2_(0.2),
40  WeTransition_(100.0),
41  AWe_(0.0)
42 {
43  if (!this->defaultCoeffs(true))
44  {
45  this->coeffDict().lookup("k1") >> k1_;
46  this->coeffDict().lookup("k2") >> k2_;
47  this->coeffDict().lookup("WeTransition") >> WeTransition_;
48  }
49 
50  scalar k21 = k2_/k1_;
51  AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow4(WeTransition_);
52 }
53 
54 
55 template<class CloudType>
57 :
58  BreakupModel<CloudType>(bum),
59  k1_(bum.k1_),
60  k2_(bum.k2_),
61  WeTransition_(bum.WeTransition_),
62  AWe_(bum.AWe_)
63 {}
64 
65 
66 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
67 
68 template<class CloudType>
70 {}
71 
72 
73 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
74 
75 template<class CloudType>
77 (
78  const scalar dt,
79  const vector& g,
80  scalar& d,
81  scalar& tc,
82  scalar& ms,
83  scalar& nParticle,
84  scalar& KHindex,
85  scalar& y,
86  scalar& yDot,
87  const scalar d0,
88  const scalar rho,
89  const scalar mu,
90  const scalar sigma,
91  const vector& U,
92  const scalar rhoc,
93  const scalar muc,
94  const vector& Urel,
95  const scalar Urmag,
96  const scalar tMom,
97  scalar& dChild,
98  scalar& massChild
99 )
100 {
101  scalar r = 0.5*d;
102  scalar r2 = r*r;
103  scalar r3 = r*r2;
104 
105  scalar semiMass = nParticle*pow3(d);
106 
107  // inverse of characteristic viscous damping time
108  scalar rtd = 0.5*this->TABCmu_*mu/(rho*r2);
109 
110  // oscillation frequency (squared)
111  scalar omega2 = this->TABComega_*sigma/(rho*r3) - rtd*rtd;
112 
113  if (omega2 > 0)
114  {
115  scalar omega = sqrt(omega2);
116  scalar romega = 1.0/omega;
117 
118  scalar We = rhoc*sqr(Urmag)*r/sigma;
119  scalar Wetmp = We/this->TABtwoWeCrit_;
120 
121  scalar y1 = y - Wetmp;
122  scalar y2 = yDot*romega;
123 
124  scalar a = sqrt(y1*y1 + y2*y2);
125 
126  // scotty we may have break-up
127  if (a + Wetmp > 1.0)
128  {
129  scalar phic = y1/a;
130 
131  // constrain phic within -1 to 1
132  phic = max(min(phic, 1), -1);
133 
134  scalar phit = acos(phic);
135  scalar phi = phit;
136  scalar quad = -y2/a;
137  if (quad < 0)
138  {
139  phi = constant::mathematical::twoPi - phit;
140  }
141 
142  scalar tb = 0;
143 
144  if (mag(y) < 1.0)
145  {
146  scalar theta = acos((1.0 - Wetmp)/a);
147 
148  if (theta < phi)
149  {
150  if (constant::mathematical::twoPi - theta >= phi)
151  {
152  theta = -theta;
153  }
155  }
156  tb = (theta - phi)*romega;
157 
158  // breakup occurs
159  if (dt > tb)
160  {
161  y = 1.0;
162  yDot = -a*omega*sin(omega*tb + phi);
163  }
164  }
165 
166  // update droplet size
167  if (dt > tb)
168  {
169  scalar sqrtWe = AWe_*pow4(We) + 1.0;
170  scalar Kbr = k1_*omega*sqrtWe;
171 
172  if (We > WeTransition_)
173  {
174  sqrtWe = sqrt(We);
175  Kbr =k2_*omega*sqrtWe;
176  }
177 
178  scalar rWetmp = 1.0/Wetmp;
179  scalar cosdtbu = max(-1.0, min(1.0, 1.0 - rWetmp));
180  scalar dtbu = romega*acos(cosdtbu);
181  scalar decay = exp(-Kbr*dtbu);
182 
183  scalar rNew = decay*r;
184  if (rNew < r)
185  {
186  d = 2.0*rNew;
187  y = 0.0;
188  yDot = 0.0;
189  }
190  }
191  }
192  }
193  else
194  {
195  // reset droplet distortion parameters
196  y = 0;
197  yDot = 0;
198  }
199 
200  // update the nParticle count to conserve mass
201  nParticle = semiMass/pow3(d);
202 
203  // Do not add child parcel
204  return false;
205 }
206 
207 
208 // ************************************************************************* //
ETAB(const dictionary &, CloudType &)
Construct from dictionary.
Definition: ETAB.C:32
surfaceScalarField & phi
dictionary dict
dimensionedScalar acos(const dimensionedScalar &ds)
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
virtual ~ETAB()
Destructor.
Definition: ETAB.C:69
dimensionedScalar exp(const dimensionedScalar &ds)
dimensionedScalar y1(const dimensionedScalar &ds)
const scalar twoPi(2 *pi)
dimensionedScalar sin(const dimensionedScalar &ds)
The Enhanced TAB model.
Definition: ETAB.H:63
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar pow3(const dimensionedScalar &ds)
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 properties.
Definition: ETAB.C:77
surfaceScalarField phic(mixture.cAlpha() *mag(phi/mesh.magSf()))
dimensionedScalar pow4(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
Templated break-up model class.
Definition: SprayCloud.H:50
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69