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