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) 2011-2016 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  volScalarField alpha(max(alpha_, scalar(0)));
365  const volScalarField& rho = phase_.rho();
366  const surfaceScalarField& alphaRhoPhi = alphaRhoPhi_;
367  const volVectorField& U = U_;
368  const volVectorField& Uc_ =
369  refCast<const twoPhaseSystem>(phase_.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  refCast<const twoPhaseSystem>(phase_.fluid()).drag(phase_).K()
414  );
415 
416  // Eq. 3.25, p. 50 Js = J1 - J2
417  volScalarField J1("J1", 3.0*beta);
418  volScalarField J2
419  (
420  "J2",
421  0.25*sqr(beta)*da*magSqr(U - Uc_)
422  /(
423  max(alpha, residualAlpha_)*rho
424  *sqrtPi*(ThetaSqrt + ThetaSmallSqrt)
425  )
426  );
427 
428  // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
429  volScalarField PsCoeff
430  (
431  granularPressureModel_->granularPressureCoeff
432  (
433  alpha,
434  gs0_,
435  rho,
436  e_
437  )
438  );
439 
440  // 'thermal' conductivity (Table 3.3, p. 49)
441  kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
442 
443  fv::options& fvOptions(fv::options::New(mesh_));
444 
445  // Construct the granular temperature equation (Eq. 3.20, p. 44)
446  // NB. note that there are two typos in Eq. 3.20:
447  // Ps should be without grad
448  // the laplacian has the wrong sign
449  fvScalarMatrix ThetaEqn
450  (
451  1.5*
452  (
453  fvm::ddt(alpha, rho, Theta_)
454  + fvm::div(alphaRhoPhi, Theta_)
455  - fvc::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaRhoPhi), Theta_)
456  )
457  - fvm::laplacian(kappa_, Theta_, "laplacian(kappa,Theta)")
458  ==
459  - fvm::SuSp((PsCoeff*I) && gradU, Theta_)
460  + (tau && gradU)
461  + fvm::Sp(-gammaCoeff, Theta_)
462  + fvm::Sp(-J1, Theta_)
463  + fvm::Sp(J2/(Theta_ + ThetaSmall), Theta_)
464  + fvOptions(alpha, rho, Theta_)
465  );
466 
467  ThetaEqn.relax();
468  fvOptions.constrain(ThetaEqn);
469  ThetaEqn.solve();
470  fvOptions.correct(Theta_);
471  }
472  else
473  {
474  // Equilibrium => dissipation == production
475  // Eq. 4.14, p.82
476  volScalarField K1("K1", 2.0*(1.0 + e_)*rho*gs0_);
478  (
479  "K3",
480  0.5*da*rho*
481  (
482  (sqrtPi/(3.0*(3.0 - e_)))
483  *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha*gs0_)
484  +1.6*alpha*gs0_*(1.0 + e_)/sqrtPi
485  )
486  );
487 
489  (
490  "K2",
491  4.0*da*rho*(1.0 + e_)*alpha*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0
492  );
493 
494  volScalarField K4("K4", 12.0*(1.0 - sqr(e_))*rho*gs0_/(da*sqrtPi));
495 
496  volScalarField trD
497  (
498  "trD",
499  alpha/(alpha + residualAlpha_)
500  *fvc::div(phi_)
501  );
502  volScalarField tr2D("tr2D", sqr(trD));
503  volScalarField trD2("trD2", tr(D & D));
504 
505  volScalarField t1("t1", K1*alpha + rho);
506  volScalarField l1("l1", -t1*trD);
507  volScalarField l2("l2", sqr(t1)*tr2D);
508  volScalarField l3
509  (
510  "l3",
511  4.0
512  *K4
513  *alpha
514  *(2.0*K3*trD2 + K2*tr2D)
515  );
516 
517  Theta_ = sqr
518  (
519  (l1 + sqrt(l2 + l3))
520  /(2.0*max(alpha, residualAlpha_)*K4)
521  );
522 
523  kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
524  }
525 
526  Theta_.max(0);
527  Theta_.min(100);
528 
529  {
530  // particle viscosity (Table 3.2, p.47)
531  nut_ = viscosityModel_->nu(alpha, Theta_, gs0_, rho, da, e_);
532 
533  volScalarField ThetaSqrt("sqrtTheta", sqrt(Theta_));
534 
535  // Bulk viscosity p. 45 (Lun et al. 1984).
536  lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;
537 
538  // Frictional pressure
539  volScalarField pf
540  (
541  frictionalStressModel_->frictionalPressure
542  (
543  phase_,
544  alphaMinFriction_,
545  alphaMax_
546  )
547  );
548 
549  nuFric_ = frictionalStressModel_->nu
550  (
551  phase_,
552  alphaMinFriction_,
553  alphaMax_,
554  pf/rho,
555  D
556  );
557 
558  // Limit viscosity and add frictional viscosity
559  nut_.min(maxNut_);
560  nuFric_ = min(nuFric_, maxNut_ - nut_);
561  nut_ += nuFric_;
562  }
563 
564  if (debug)
565  {
566  Info<< typeName << ':' << nl
567  << " max(Theta) = " << max(Theta_).value() << nl
568  << " max(nut) = " << max(nut_).value() << endl;
569  }
570 }
571 
572 
573 // ************************************************************************* //
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
ThermalDiffusivity< PhaseCompressibleTurbulenceModel< phaseModel > > phaseCompressibleTurbulenceModel
Typedef for phaseCompressibleTurbulenceModel.
#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
tmp< fvMatrix< Type > > Sp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
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:253
virtual tmp< volScalarField > pPrime() const
Return the phase-pressure&#39;.
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
virtual tmp< volScalarField > k() const
Return the turbulence kinetic energy.
Info<< "Reading strained laminar flame speed field Su\n"<< endl;volScalarField Su(IOobject("Su", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Reading field betav\n"<< endl;volScalarField betav(IOobject("betav", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);Info<< "Reading field Lobs\n"<< endl;volScalarField Lobs(IOobject("Lobs", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);Info<< "Reading field CT\n"<< endl;volSymmTensorField CT(IOobject("CT", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);Info<< "Reading field Nv\n"<< endl;volScalarField Nv(IOobject("Nv", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);Info<< "Reading field nsv\n"<< endl;volSymmTensorField nsv(IOobject("nsv", mesh.facesInstance(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE), mesh);IOdictionary PDRProperties(IOobject("PDRProperties", runTime.constant(), mesh, IOobject::MUST_READ_IF_MODIFIED, IOobject::NO_WRITE));autoPtr< PDRDragModel > drag
Definition: createFields.H:182
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
fv::options & fvOptions
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
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
tmp< fvMatrix< Type > > SuSp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
static word groupName(Name name, const word &group)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
virtual tmp< volSymmTensorField > devRhoReff() const
Return the effective stress tensor.
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)
virtual tmp< surfaceScalarField > pPrimef() const
Return the face-phase-pressure&#39;.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
virtual ~kineticTheoryModel()
Destructor.
static const char nl
Definition: Ostream.H:262
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
virtual tmp< volSymmTensorField > R() const
Return the Reynolds stress tensor.
virtual tmp< fvVectorMatrix > divDevRhoReff(volVectorField &U) const
Return the source term for the momentum equation.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
RASModel< EddyDiffusivity< turbulenceModel > > RASModel
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
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 bool read()
Re-read model coefficients if they have changed.
messageStream Info
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
virtual tmp< volScalarField > epsilon() const
Return the turbulence kinetic energy dissipation rate.
A class for managing temporary objects.
Definition: PtrList.H:54
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
static options & New(const fvMesh &mesh)
Construct fvOptions and register to datbase if not present.
Definition: fvOptions.C:103
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:366