FreeFOAM The Cross-Platform CFD Toolkit
kineticTheoryModel.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 "kineticTheoryModel.H"
29 #include <finiteVolume/fvCFD.H>
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 Foam::kineticTheoryModel::kineticTheoryModel
34 (
35  const Foam::phaseModel& phasea,
36  const Foam::volVectorField& Ub,
37  const Foam::volScalarField& alpha,
38  const Foam::dragModel& draga
39 )
40 :
41  phasea_(phasea),
42  Ua_(phasea.U()),
43  Ub_(Ub),
44  alpha_(alpha),
45  phia_(phasea.phi()),
46  draga_(draga),
47 
48  rhoa_(phasea.rho()),
49  da_(phasea.d()),
50  nua_(phasea.nu()),
51 
52  kineticTheoryProperties_
53  (
54  IOobject
55  (
56  "kineticTheoryProperties",
57  Ua_.time().constant(),
58  Ua_.mesh(),
59  IOobject::MUST_READ,
60  IOobject::NO_WRITE
61  )
62  ),
63  kineticTheory_(kineticTheoryProperties_.lookup("kineticTheory")),
64  equilibrium_(kineticTheoryProperties_.lookup("equilibrium")),
65 
66  viscosityModel_
67  (
68  kineticTheoryModels::viscosityModel::New
69  (
70  kineticTheoryProperties_
71  )
72  ),
73  conductivityModel_
74  (
76  (
77  kineticTheoryProperties_
78  )
79  ),
80  radialModel_
81  (
82  radialModel::New
83  (
84  kineticTheoryProperties_
85  )
86  ),
87  granularPressureModel_
88  (
90  (
91  kineticTheoryProperties_
92  )
93  ),
94  frictionalStressModel_
95  (
97  (
98  kineticTheoryProperties_
99  )
100  ),
101  e_(kineticTheoryProperties_.lookup("e")),
102  alphaMax_(kineticTheoryProperties_.lookup("alphaMax")),
103  alphaMinFriction_(kineticTheoryProperties_.lookup("alphaMinFriction")),
104  Fr_(kineticTheoryProperties_.lookup("Fr")),
105  eta_(kineticTheoryProperties_.lookup("eta")),
106  p_(kineticTheoryProperties_.lookup("p")),
107  phi_(dimensionedScalar(kineticTheoryProperties_.lookup("phi"))*M_PI/180.0),
108  Theta_
109  (
110  IOobject
111  (
112  "Theta",
113  Ua_.time().timeName(),
114  Ua_.mesh(),
115  IOobject::MUST_READ,
116  IOobject::AUTO_WRITE
117  ),
118  Ua_.mesh()
119  ),
120  mua_
121  (
122  IOobject
123  (
124  "mua",
125  Ua_.time().timeName(),
126  Ua_.mesh(),
127  IOobject::NO_READ,
128  IOobject::AUTO_WRITE
129  ),
130  Ua_.mesh(),
131  dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
132  ),
133  lambda_
134  (
135  IOobject
136  (
137  "lambda",
138  Ua_.time().timeName(),
139  Ua_.mesh(),
140  IOobject::NO_READ,
141  IOobject::NO_WRITE
142  ),
143  Ua_.mesh(),
144  dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
145  ),
146  pa_
147  (
148  IOobject
149  (
150  "pa",
151  Ua_.time().timeName(),
152  Ua_.mesh(),
153  IOobject::NO_READ,
154  IOobject::AUTO_WRITE
155  ),
156  Ua_.mesh(),
157  dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
158  ),
159  kappa_
160  (
161  IOobject
162  (
163  "kappa",
164  Ua_.time().timeName(),
165  Ua_.mesh(),
166  IOobject::NO_READ,
167  IOobject::NO_WRITE
168  ),
169  Ua_.mesh(),
170  dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
171  ),
172  gs0_
173  (
174  IOobject
175  (
176  "gs0",
177  Ua_.time().timeName(),
178  Ua_.mesh(),
179  IOobject::NO_READ,
180  IOobject::NO_WRITE
181  ),
182  Ua_.mesh(),
183  dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 1.0)
184  )
185 {}
186 
187 
188 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
189 
191 {}
192 
193 
194 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
195 
197 {
198  if (!kineticTheory_)
199  {
200  return;
201  }
202 
203  const scalar sqrtPi = sqrt(mathematicalConstant::pi);
204 
205  surfaceScalarField phi = 1.5*rhoa_*phia_*fvc::interpolate(alpha_);
206 
207  volTensorField dU = gradUat.T();//fvc::grad(Ua_);
208  volSymmTensorField D = symm(dU);
209 
210  // NB, drag = K*alpha*beta,
211  // (the alpha and beta has been extracted from the drag function for
212  // numerical reasons)
213  volScalarField Ur = mag(Ua_ - Ub_);
214  volScalarField betaPrim = alpha_*(1.0 - alpha_)*draga_.K(Ur);
215 
216  // Calculating the radial distribution function (solid volume fraction is
217  // limited close to the packing limit, but this needs improvements)
218  // The solution is higly unstable close to the packing limit.
219  gs0_ = radialModel_->g0
220  (
221  min(max(alpha_, scalar(1e-6)), alphaMax_ - 0.01),
222  alphaMax_
223  );
224 
225  // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
226  volScalarField PsCoeff = granularPressureModel_->granularPressureCoeff
227  (
228  alpha_,
229  gs0_,
230  rhoa_,
231  e_
232  );
233 
234  // 'thermal' conductivity (Table 3.3, p. 49)
235  kappa_ = conductivityModel_->kappa(alpha_, Theta_, gs0_, rhoa_, da_, e_);
236 
237  // particle viscosity (Table 3.2, p.47)
238  mua_ = viscosityModel_->mua(alpha_, Theta_, gs0_, rhoa_, da_, e_);
239 
240  dimensionedScalar Tsmall
241  (
242  "small",
243  dimensionSet(0 , 2 ,-2 ,0 , 0, 0, 0),
244  1.0e-6
245  );
246 
247  dimensionedScalar TsmallSqrt = sqrt(Tsmall);
248  volScalarField ThetaSqrt = sqrt(Theta_);
249 
250  // dissipation (Eq. 3.24, p.50)
251  volScalarField gammaCoeff =
252  12.0*(1.0 - sqr(e_))*sqr(alpha_)*rhoa_*gs0_*(1.0/da_)*ThetaSqrt/sqrtPi;
253 
254  // Eq. 3.25, p. 50 Js = J1 - J2
255  volScalarField J1 = 3.0*betaPrim;
256  volScalarField J2 =
257  0.25*sqr(betaPrim)*da_*sqr(Ur)
258  /(max(alpha_, scalar(1e-6))*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt));
259 
260  // bulk viscosity p. 45 (Lun et al. 1984).
261  lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi;
262 
263  // stress tensor, Definitions, Table 3.1, p. 43
264  volSymmTensorField tau = 2.0*mua_*D + (lambda_ - (2.0/3.0)*mua_)*tr(D)*I;
265 
266  if (!equilibrium_)
267  {
268  // construct the granular temperature equation (Eq. 3.20, p. 44)
269  // NB. note that there are two typos in Eq. 3.20
270  // no grad infront of Ps
271  // wrong sign infront of laplacian
272  fvScalarMatrix ThetaEqn
273  (
274  fvm::ddt(1.5*alpha_*rhoa_, Theta_)
275  + fvm::div(phi, Theta_, "div(phi,Theta)")
276  ==
277  fvm::SuSp(-((PsCoeff*I) && dU), Theta_)
278  + (tau && dU)
279  + fvm::laplacian(kappa_, Theta_, "laplacian(kappa,Theta)")
280  + fvm::Sp(-gammaCoeff, Theta_)
281  + fvm::Sp(-J1, Theta_)
282  + fvm::Sp(J2/(Theta_ + Tsmall), Theta_)
283  );
284 
285  ThetaEqn.relax();
286  ThetaEqn.solve();
287  }
288  else
289  {
290  // equilibrium => dissipation == production
291  // Eq. 4.14, p.82
292  volScalarField K1 = 2.0*(1.0 + e_)*rhoa_*gs0_;
293  volScalarField K3 = 0.5*da_*rhoa_*
294  (
295  (sqrtPi/(3.0*(3.0-e_)))
296  *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha_*gs0_)
297  +1.6*alpha_*gs0_*(1.0 + e_)/sqrtPi
298  );
299 
301  4.0*da_*rhoa_*(1.0 + e_)*alpha_*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0;
302 
303  volScalarField K4 = 12.0*(1.0 - sqr(e_))*rhoa_*gs0_/(da_*sqrtPi);
304 
305  volScalarField trD = tr(D);
306  volScalarField tr2D = sqr(trD);
307  volScalarField trD2 = tr(D & D);
308 
309  volScalarField t1 = K1*alpha_ + rhoa_;
310  volScalarField l1 = -t1*trD;
311  volScalarField l2 = sqr(t1)*tr2D;
312  volScalarField l3 = 4.0*K4*max(alpha_, scalar(1e-6))*(2.0*K3*trD2 + K2*tr2D);
313 
314  Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha_ + 1.0e-4)*K4));
315  }
316 
317  Theta_.max(1.0e-15);
318  Theta_.min(1.0e+3);
319 
320  volScalarField pf = frictionalStressModel_->frictionalPressure
321  (
322  alpha_,
323  alphaMinFriction_,
324  alphaMax_,
325  Fr_,
326  eta_,
327  p_
328  );
329 
330  PsCoeff += pf/(Theta_+Tsmall);
331 
332  PsCoeff.min(1.0e+10);
333  PsCoeff.max(-1.0e+10);
334 
335  // update particle pressure
336  pa_ = PsCoeff*Theta_;
337 
338  // frictional shear stress, Eq. 3.30, p. 52
339  volScalarField muf = frictionalStressModel_->muf
340  (
341  alpha_,
342  alphaMax_,
343  pf,
344  D,
345  phi_
346  );
347 
348  // add frictional stress
349  mua_ += muf;
350  mua_.min(1.0e+2);
351  mua_.max(0.0);
352 
353  Info<< "kinTheory: max(Theta) = " << max(Theta_).value() << endl;
354 
355  volScalarField ktn = mua_/rhoa_;
356 
357  Info<< "kinTheory: min(nua) = " << min(ktn).value()
358  << ", max(nua) = " << max(ktn).value() << endl;
359 
360  Info<< "kinTheory: min(pa) = " << min(pa_).value()
361  << ", max(pa) = " << max(pa_).value() << endl;
362 }
363 
364 
365 // ************************ vim: set sw=4 sts=4 et: ************************ //