ETAB.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-2023 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  const label injectori,
98  scalar& dChild,
99  scalar& massChild
100 )
101 {
102  scalar r = 0.5*d;
103  scalar r2 = r*r;
104  scalar r3 = r*r2;
105 
106  scalar semiMass = nParticle*pow3(d);
107 
108  // inverse of characteristic viscous damping time
109  scalar rtd = 0.5*this->TABCmu_*mu/(rho*r2);
110 
111  // oscillation frequency (squared)
112  scalar omega2 = this->TABComega_*sigma/(rho*r3) - rtd*rtd;
113 
114  if (omega2 > 0)
115  {
116  scalar omega = sqrt(omega2);
117  scalar romega = 1.0/omega;
118 
119  scalar We = rhoc*sqr(Urmag)*r/sigma;
120  scalar Wetmp = We/this->TABtwoWeCrit_;
121 
122  scalar y1 = y - Wetmp;
123  scalar y2 = yDot*romega;
124 
125  scalar a = sqrt(y1*y1 + y2*y2);
126 
127  // scotty we may have break-up
128  if (a + Wetmp > 1.0)
129  {
130  scalar phic = y1/a;
131 
132  // constrain phic within -1 to 1
133  phic = max(min(phic, 1), -1);
134 
135  scalar phit = acos(phic);
136  scalar phi = phit;
137  scalar quad = -y2/a;
138  if (quad < 0)
139  {
140  phi = constant::mathematical::twoPi - phit;
141  }
142 
143  scalar tb = 0;
144 
145  if (mag(y) < 1.0)
146  {
147  scalar theta = acos((1.0 - Wetmp)/a);
148 
149  if (theta < phi)
150  {
151  if (constant::mathematical::twoPi - theta >= phi)
152  {
153  theta = -theta;
154  }
156  }
157  tb = (theta - phi)*romega;
158 
159  // breakup occurs
160  if (dt > tb)
161  {
162  y = 1.0;
163  yDot = -a*omega*sin(omega*tb + phi);
164  }
165  }
166 
167  // update droplet size
168  if (dt > tb)
169  {
170  scalar sqrtWe = AWe_*pow4(We) + 1.0;
171  scalar Kbr = k1_*omega*sqrtWe;
172 
173  if (We > WeTransition_)
174  {
175  sqrtWe = sqrt(We);
176  Kbr =k2_*omega*sqrtWe;
177  }
178 
179  scalar rWetmp = 1.0/Wetmp;
180  scalar cosdtbu = max(-1.0, min(1.0, 1.0 - rWetmp));
181  scalar dtbu = romega*acos(cosdtbu);
182  scalar decay = exp(-Kbr*dtbu);
183 
184  scalar rNew = decay*r;
185  if (rNew < r)
186  {
187  d = 2.0*rNew;
188  y = 0.0;
189  yDot = 0.0;
190  }
191  }
192  }
193  }
194  else
195  {
196  // reset droplet distortion parameters
197  y = 0;
198  yDot = 0;
199  }
200 
201  // update the nParticle count to conserve mass
202  nParticle = semiMass/pow3(d);
203 
204  // Do not add child parcel
205  return false;
206 }
207 
208 
209 // ************************************************************************* //
scalar y
Templated break-up model class.
Definition: BreakupModel.H:55
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:80
The Enhanced TAB model.
Definition: ETAB.H:66
virtual ~ETAB()
Destructor.
Definition: ETAB.C:69
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 properties.
Definition: ETAB.C:77
ETAB(const dictionary &, CloudType &)
Construct from dictionary.
Definition: ETAB.C:32
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:710
Convenience class to handle the input of constant rotational speed. Reads an omega entry with default...
Definition: omega.H:54
const dictionary & coeffDict() const
Return const access to the coefficients dictionary.
Definition: subModelBase.C:128
virtual bool defaultCoeffs(const bool printMsg) const
Returns true if defaultCoeffs is true and outputs on printMsg.
Definition: subModelBase.C:140
const scalar omega
U
Definition: pEqn.H:72
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)
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 > &)
dimensionedScalar pow4(const dimensionedScalar &ds)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar acos(const dimensionedScalar &ds)
dictionary dict