kineticTheoryModel.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-2018 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 "kineticTheoryModel.H"
27 #include "mathematicalConstants.H"
28 #include "twoPhaseSystem.H"
29 #include "fvOptions.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 Foam::RASModels::kineticTheoryModel::kineticTheoryModel
34 (
35  const volScalarField& alpha,
36  const volScalarField& rho,
37  const volVectorField& U,
38  const surfaceScalarField& alphaRhoPhi,
39  const surfaceScalarField& phi,
40  const transportModel& phase,
41  const word& propertiesName,
42  const word& type
43 )
44 :
45  eddyViscosity
46  <
48  >
49  (
50  type,
51  alpha,
52  rho,
53  U,
54  alphaRhoPhi,
55  phi,
56  phase,
57  propertiesName
58  ),
59 
60  phase_(phase),
61 
62  viscosityModel_
63  (
64  kineticTheoryModels::viscosityModel::New
65  (
66  coeffDict_
67  )
68  ),
69  conductivityModel_
70  (
71  kineticTheoryModels::conductivityModel::New
72  (
73  coeffDict_
74  )
75  ),
76  radialModel_
77  (
78  kineticTheoryModels::radialModel::New
79  (
80  coeffDict_
81  )
82  ),
83  granularPressureModel_
84  (
85  kineticTheoryModels::granularPressureModel::New
86  (
87  coeffDict_
88  )
89  ),
90  frictionalStressModel_
91  (
92  kineticTheoryModels::frictionalStressModel::New
93  (
94  coeffDict_
95  )
96  ),
97 
98  equilibrium_(coeffDict_.lookup("equilibrium")),
99  e_("e", dimless, coeffDict_),
100  alphaMax_("alphaMax", dimless, coeffDict_),
101  alphaMinFriction_
102  (
103  "alphaMinFriction",
104  dimless,
105  coeffDict_
106  ),
107  residualAlpha_
108  (
109  "residualAlpha",
110  dimless,
111  coeffDict_
112  ),
113 
114  maxNut_
115  (
116  "maxNut",
117  dimensionSet(0,2,-1,0,0),
118  coeffDict_.lookupOrDefault<scalar>("maxNut",1000)
119  ),
120 
121  Theta_
122  (
123  IOobject
124  (
125  IOobject::groupName("Theta", phase.name()),
126  U.time().timeName(),
127  U.mesh(),
128  IOobject::MUST_READ,
129  IOobject::AUTO_WRITE
130  ),
131  U.mesh()
132  ),
133 
134  lambda_
135  (
136  IOobject
137  (
138  IOobject::groupName("lambda", phase.name()),
139  U.time().timeName(),
140  U.mesh(),
141  IOobject::NO_READ,
142  IOobject::NO_WRITE
143  ),
144  U.mesh(),
145  dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
146  ),
147 
148  gs0_
149  (
150  IOobject
151  (
152  IOobject::groupName("gs0", phase.name()),
153  U.time().timeName(),
154  U.mesh(),
155  IOobject::NO_READ,
156  IOobject::NO_WRITE
157  ),
158  U.mesh(),
159  dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 0.0)
160  ),
161 
162  kappa_
163  (
164  IOobject
165  (
166  IOobject::groupName("kappa", phase.name()),
167  U.time().timeName(),
168  U.mesh(),
169  IOobject::NO_READ,
170  IOobject::NO_WRITE
171  ),
172  U.mesh(),
173  dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
174  ),
175 
176  nuFric_
177  (
178  IOobject
179  (
180  IOobject::groupName("nuFric", phase.name()),
181  U.time().timeName(),
182  U.mesh(),
183  IOobject::NO_READ,
184  IOobject::AUTO_WRITE
185  ),
186  U.mesh(),
187  dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
188  )
189 {
190  if (type == typeName)
191  {
192  printCoeffs(type);
193  }
194 }
195 
196 
197 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
198 
200 {}
201 
202 
203 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
204 
206 {
207  if
208  (
209  eddyViscosity
210  <
211  RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>
212  >::read()
213  )
214  {
215  coeffDict().lookup("equilibrium") >> equilibrium_;
216  e_.readIfPresent(coeffDict());
217  alphaMax_.readIfPresent(coeffDict());
218  alphaMinFriction_.readIfPresent(coeffDict());
219 
220  viscosityModel_->read();
221  conductivityModel_->read();
222  radialModel_->read();
223  granularPressureModel_->read();
224  frictionalStressModel_->read();
225 
226  return true;
227  }
228  else
229  {
230  return false;
231  }
232 }
233 
234 
237 {
239  return nut_;
240 }
241 
242 
245 {
247  return nut_;
248 }
249 
250 
253 {
254  return tmp<volSymmTensorField>
255  (
257  (
258  IOobject
259  (
260  IOobject::groupName("R", U_.group()),
261  runTime_.timeName(),
262  mesh_,
265  ),
266  - (nut_)*dev(twoSymm(fvc::grad(U_)))
267  - (lambda_*fvc::div(phi_))*symmTensor::I
268  )
269  );
270 }
271 
272 
275 {
276  const volScalarField& rho = phase_.rho();
277 
278  tmp<volScalarField> tpPrime
279  (
280  Theta_
281  *granularPressureModel_->granularPressureCoeffPrime
282  (
283  alpha_,
284  radialModel_->g0(alpha_, alphaMinFriction_, alphaMax_),
285  radialModel_->g0prime(alpha_, alphaMinFriction_, alphaMax_),
286  rho,
287  e_
288  )
289  + frictionalStressModel_->frictionalPressurePrime
290  (
291  phase_,
292  alphaMinFriction_,
293  alphaMax_
294  )
295  );
296 
297  volScalarField::Boundary& bpPrime =
298  tpPrime.ref().boundaryFieldRef();
299 
300  forAll(bpPrime, patchi)
301  {
302  if (!bpPrime[patchi].coupled())
303  {
304  bpPrime[patchi] == 0;
305  }
306  }
307 
308  return tpPrime;
309 }
310 
311 
314 {
315  return fvc::interpolate(pPrime());
316 }
317 
318 
321 {
322  return tmp<volSymmTensorField>
323  (
325  (
326  IOobject
327  (
328  IOobject::groupName("devRhoReff", U_.group()),
329  runTime_.timeName(),
330  mesh_,
333  ),
334  - (rho_*nut_)
335  *dev(twoSymm(fvc::grad(U_)))
336  - ((rho_*lambda_)*fvc::div(phi_))*symmTensor::I
337  )
338  );
339 }
340 
341 
344 (
345  volVectorField& U
346 ) const
347 {
348  return
349  (
350  - fvm::laplacian(rho_*nut_, U)
351  - fvc::div
352  (
353  (rho_*nut_)*dev2(T(fvc::grad(U)))
354  + ((rho_*lambda_)*fvc::div(phi_))
355  *dimensioned<symmTensor>("I", dimless, symmTensor::I)
356  )
357  );
358 }
359 
360 
362 {
363  // Local references
364  const twoPhaseSystem& fluid = refCast<const twoPhaseSystem>(phase_.fluid());
365  volScalarField alpha(max(alpha_, scalar(0)));
366  const volScalarField& rho = phase_.rho();
367  const surfaceScalarField& alphaRhoPhi = alphaRhoPhi_;
368  const volVectorField& U = U_;
369  const volVectorField& Uc_ = fluid.otherPhase(phase_).U();
370 
371  const scalar sqrtPi = sqrt(constant::mathematical::pi);
372  dimensionedScalar ThetaSmall("ThetaSmall", Theta_.dimensions(), 1.0e-6);
373  dimensionedScalar ThetaSmallSqrt(sqrt(ThetaSmall));
374 
375  tmp<volScalarField> tda(phase_.d());
376  const volScalarField& da = tda();
377 
378  tmp<volTensorField> tgradU(fvc::grad(U_));
379  const volTensorField& gradU(tgradU());
380  volSymmTensorField D(symm(gradU));
381 
382  // Calculating the radial distribution function
383  gs0_ = radialModel_->g0(alpha, alphaMinFriction_, alphaMax_);
384 
385  if (!equilibrium_)
386  {
387  // Particle viscosity (Table 3.2, p.47)
388  nut_ = viscosityModel_->nu(alpha, Theta_, gs0_, rho, da, e_);
389 
390  volScalarField ThetaSqrt("sqrtTheta", sqrt(Theta_));
391 
392  // Bulk viscosity p. 45 (Lun et al. 1984).
393  lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;
394 
395  // Stress tensor, Definitions, Table 3.1, p. 43
397  (
398  rho*(2.0*nut_*D + (lambda_ - (2.0/3.0)*nut_)*tr(D)*I)
399  );
400 
401  // Dissipation (Eq. 3.24, p.50)
402  volScalarField gammaCoeff
403  (
404  "gammaCoeff",
405  12.0*(1.0 - sqr(e_))
406  *max(sqr(alpha), residualAlpha_)
407  *rho*gs0_*(1.0/da)*ThetaSqrt/sqrtPi
408  );
409 
410  // Drag
412  (
413  fluid.lookupSubModel<dragModel>
414  (
415  phase_,
416  fluid.otherPhase(phase_)
417  ).K()
418  );
419 
420  // Eq. 3.25, p. 50 Js = J1 - J2
421  volScalarField J1("J1", 3.0*beta);
422  volScalarField J2
423  (
424  "J2",
425  0.25*sqr(beta)*da*magSqr(U - Uc_)
426  /(
427  max(alpha, residualAlpha_)*rho
428  *sqrtPi*(ThetaSqrt + ThetaSmallSqrt)
429  )
430  );
431 
432  // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
433  volScalarField PsCoeff
434  (
435  granularPressureModel_->granularPressureCoeff
436  (
437  alpha,
438  gs0_,
439  rho,
440  e_
441  )
442  );
443 
444  // 'thermal' conductivity (Table 3.3, p. 49)
445  kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
446 
447  fv::options& fvOptions(fv::options::New(mesh_));
448 
449  // Construct the granular temperature equation (Eq. 3.20, p. 44)
450  // NB. note that there are two typos in Eq. 3.20:
451  // Ps should be without grad
452  // the laplacian has the wrong sign
453  fvScalarMatrix ThetaEqn
454  (
455  1.5*
456  (
457  fvm::ddt(alpha, rho, Theta_)
458  + fvm::div(alphaRhoPhi, Theta_)
459  - fvc::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaRhoPhi), Theta_)
460  )
461  - fvm::laplacian(kappa_, Theta_, "laplacian(kappa,Theta)")
462  ==
463  - fvm::SuSp((PsCoeff*I) && gradU, Theta_)
464  + (tau && gradU)
465  + fvm::Sp(-gammaCoeff, Theta_)
466  + fvm::Sp(-J1, Theta_)
467  + fvm::Sp(J2/(Theta_ + ThetaSmall), Theta_)
468  + fvOptions(alpha, rho, Theta_)
469  );
470 
471  ThetaEqn.relax();
472  fvOptions.constrain(ThetaEqn);
473  ThetaEqn.solve();
474  fvOptions.correct(Theta_);
475  }
476  else
477  {
478  // Equilibrium => dissipation == production
479  // Eq. 4.14, p.82
480  volScalarField K1("K1", 2.0*(1.0 + e_)*rho*gs0_);
482  (
483  "K3",
484  0.5*da*rho*
485  (
486  (sqrtPi/(3.0*(3.0 - e_)))
487  *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha*gs0_)
488  +1.6*alpha*gs0_*(1.0 + e_)/sqrtPi
489  )
490  );
491 
493  (
494  "K2",
495  4.0*da*rho*(1.0 + e_)*alpha*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0
496  );
497 
498  volScalarField K4("K4", 12.0*(1.0 - sqr(e_))*rho*gs0_/(da*sqrtPi));
499 
500  volScalarField trD
501  (
502  "trD",
503  alpha/(alpha + residualAlpha_)
504  *fvc::div(phi_)
505  );
506  volScalarField tr2D("tr2D", sqr(trD));
507  volScalarField trD2("trD2", tr(D & D));
508 
509  volScalarField t1("t1", K1*alpha + rho);
510  volScalarField l1("l1", -t1*trD);
511  volScalarField l2("l2", sqr(t1)*tr2D);
512  volScalarField l3
513  (
514  "l3",
515  4.0
516  *K4
517  *alpha
518  *(2.0*K3*trD2 + K2*tr2D)
519  );
520 
521  Theta_ = sqr
522  (
523  (l1 + sqrt(l2 + l3))
524  /(2.0*max(alpha, residualAlpha_)*K4)
525  );
526 
527  kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
528  }
529 
530  Theta_.max(0);
531  Theta_.min(100);
532 
533  {
534  // particle viscosity (Table 3.2, p.47)
535  nut_ = viscosityModel_->nu(alpha, Theta_, gs0_, rho, da, e_);
536 
537  volScalarField ThetaSqrt("sqrtTheta", sqrt(Theta_));
538 
539  // Bulk viscosity p. 45 (Lun et al. 1984).
540  lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;
541 
542  // Frictional pressure
543  volScalarField pf
544  (
545  frictionalStressModel_->frictionalPressure
546  (
547  phase_,
548  alphaMinFriction_,
549  alphaMax_
550  )
551  );
552 
553  nuFric_ = frictionalStressModel_->nu
554  (
555  phase_,
556  alphaMinFriction_,
557  alphaMax_,
558  pf/rho,
559  D
560  );
561 
562  // Limit viscosity and add frictional viscosity
563  nut_.min(maxNut_);
564  nuFric_ = min(nuFric_, maxNut_ - nut_);
565  nut_ += nuFric_;
566  }
567 
568  if (debug)
569  {
570  Info<< typeName << ':' << nl
571  << " max(Theta) = " << max(Theta_).value() << nl
572  << " max(nut) = " << max(nut_).value() << endl;
573  }
574 }
575 
576 
577 // ************************************************************************* //
GeometricField< symmTensor, fvPatchField, volMesh > volSymmTensorField
Definition: volFieldsFwd.H:58
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
virtual tmp< volScalarField > epsilon() const
Return the turbulence kinetic energy dissipation rate.
ThermalDiffusivity< PhaseCompressibleTurbulenceModel< phaseModel > > phaseCompressibleTurbulenceModel
Typedef for phaseCompressibleTurbulenceModel.
tmp< fvMatrix< Type > > SuSp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
fv::options & fvOptions
#define K1
Definition: SHA1.C:167
#define K4
Definition: SHA1.C:170
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define K2
Definition: SHA1.C:168
virtual tmp< surfaceScalarField > pPrimef() const
Return the face-phase-pressure&#39;.
GeometricField< tensor, fvPatchField, volMesh > volTensorField
Definition: volFieldsFwd.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
#define K3
Definition: SHA1.C:169
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
CGAL::Exact_predicates_exact_constructions_kernel K
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
virtual tmp< volSymmTensorField > R() const
Return the Reynolds stress tensor.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
static const SymmTensor I
Definition: SymmTensor.H:72
virtual tmp< volSymmTensorField > devRhoReff() const
Return the effective stress tensor.
static word groupName(Name name, const word &group)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
virtual void correct()
Solve the kinetic theory equations and correct the viscosity.
word timeName
Definition: getTimeIndex.H:3
void min(const dimensioned< Type > &)
dimensionedScalar tr(const dimensionedSphericalTensor &dt)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
virtual ~kineticTheoryModel()
Destructor.
static const char nl
Definition: Ostream.H:265
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
RASModel< EddyDiffusivity< turbulenceModel > > RASModel
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
virtual tmp< fvVectorMatrix > divDevRhoReff(volVectorField &U) const
Return the source term for the momentum equation.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Internal & ref()
Return a reference to the dimensioned internal field.
label patchi
dimensionedSymmTensor dev2(const dimensionedSymmTensor &dt)
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
virtual tmp< volScalarField > pPrime() const
Return the phase-pressure&#39;.
virtual bool read()
Re-read model coefficients if they have changed.
messageStream Info
virtual tmp< volScalarField > k() const
Return the turbulence kinetic energy.
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
A class for managing temporary objects.
Definition: PtrList.H:53
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
static options & New(const fvMesh &mesh)
Construct fvOptions and register to datbase if not present.
Definition: fvOptions.C:101
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:366