FreeFOAM The Cross-Platform CFD Toolkit
TAB.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) 1991-2010 OpenCFD Ltd.
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"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
38 
40 (
41  breakupModel,
42  TAB,
43  dictionary
44 );
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
48 // Construct from components
50 (
51  const dictionary& dict,
52  spray& sm
53 )
54 :
55  breakupModel(dict, sm),
56  coeffsDict_(dict.subDict(typeName + "Coeffs")),
57  Cmu_(readScalar(coeffsDict_.lookup("Cmu"))),
58  Comega_(readScalar(coeffsDict_.lookup("Comega"))),
59  WeCrit_(readScalar(coeffsDict_.lookup("WeCrit")))
60 {
61 
62  // calculate the inverse function of the Rossin-Rammler Distribution
63  const scalar xx0 = 12.0;
64  const scalar rrd100 = 1.0/(1.0-exp(-xx0)*(1+xx0+pow(xx0, 2)/2+pow(xx0, 3)/6));
65 
66  for(label n=0; n<100; n++)
67  {
68  scalar xx = 0.12*(n+1);
69  rrd_[n] = (1-exp(-xx)*(1 + xx + pow(xx, 2)/2 + pow(xx, 3)/6))*rrd100;
70  }
71 }
72 
73 
74 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
75 
77 {}
78 
79 
80 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
81 
83 (
84  parcel& p,
85  const scalar deltaT,
86  const vector& Ug,
87  const liquidMixture& fuels
88 ) const
89 {
90 
91  scalar T = p.T();
92  scalar pc = spray_.p()[p.cell()];
93  scalar r = 0.5*p.d();
94  scalar r2 = r*r;
95  scalar r3 = r*r2;
96 
97  scalar rho = fuels.rho(pc, T, p.X());
98  scalar sigma = fuels.sigma(pc, T, p.X());
99  scalar mu = fuels.mu(pc, T, p.X());
100 
101  // inverse of characteristic viscous damping time
102  scalar rtd = 0.5*Cmu_*mu/(rho*r2);
103 
104  // oscillation frequency (squared)
105  scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
106 
107  if (omega2 > 0)
108  {
109  scalar omega = sqrt(omega2);
110  scalar rhog = spray_.rho()[p.cell()];
111  scalar We = p.We(Ug, rhog, sigma);
112  scalar Wetmp = We/WeCrit_;
113 
114  scalar y1 = p.dev() - Wetmp;
115  scalar y2 = p.ddev()/omega;
116 
117  scalar a = sqrt(y1*y1 + y2*y2);
118 
119  // scotty we may have break-up
120  if (a+Wetmp > 1.0)
121  {
122  scalar phic = y1/a;
123 
124  // constrain phic within -1 to 1
125  phic = max(min(phic, 1), -1);
126 
127  scalar phit = acos(phic);
128  scalar phi = phit;
129  scalar quad = -y2/a;
130  if (quad < 0)
131  {
132  phi = 2*mathematicalConstant::pi - phit;
133  }
134 
135  scalar tb = 0;
136 
137  if (mag(p.dev()) < 1.0)
138  {
139  scalar coste = 1.0;
140  if
141  (
142  (Wetmp - a < -1)
143  && (p.ddev() < 0)
144  )
145  {
146  coste = -1.0;
147  }
148 
149  scalar theta = acos((coste-Wetmp)/a);
150 
151  if (theta < phi)
152  {
153  if (2*mathematicalConstant::pi-theta >= phi)
154  {
155  theta = -theta;
156  }
157  theta += 2*mathematicalConstant::pi;
158  }
159  tb = (theta-phi)/omega;
160 
161  // breakup occurs
162  if (deltaT > tb)
163  {
164  p.dev() = 1.0;
165  p.ddev() = -a*omega*sin(omega*tb + phi);
166  }
167 
168  }
169 
170  // update droplet size
171  if (deltaT > tb)
172  {
173  scalar rs = r/
174  (
175  1
176  + (4.0/3.0)*pow(p.dev(), 2)
177  + rho*r3/(8*sigma)*pow(p.ddev(), 2)
178  );
179 
180  label n = 0;
181  bool found = false;
182  scalar random = rndGen_.scalar01();
183  while (!found && (n<99))
184  {
185  if (rrd_[n]>random)
186  {
187  found = true;
188  }
189  n++;
190 
191  }
192  scalar rNew = 0.04*n*rs;
193  if (rNew < r)
194  {
195  p.d() = 2*rNew;
196  p.dev() = 0;
197  p.ddev() = 0;
198  }
199 
200  }
201 
202  }
203 
204  }
205  else
206  {
207  // reset droplet distortion parameters
208  p.dev() = 0;
209  p.ddev() = 0;
210  }
211 
212 }
213 
214 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
215 
216 } // End namespace Foam
217 
218 // ************************ vim: set sw=4 sts=4 et: ************************ //