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  <
47  >
48  (
49  type,
50  alpha,
51  rho,
52  U,
53  alphaRhoPhi,
54  phi,
55  phase,
56  propertiesName
57  ),
58 
59  phase_(phase),
60 
61  viscosityModel_
62  (
63  kineticTheoryModels::viscosityModel::New
64  (
65  coeffDict_
66  )
67  ),
68  conductivityModel_
69  (
70  kineticTheoryModels::conductivityModel::New
71  (
72  coeffDict_
73  )
74  ),
75  radialModel_
76  (
77  kineticTheoryModels::radialModel::New
78  (
79  coeffDict_
80  )
81  ),
82  granularPressureModel_
83  (
84  kineticTheoryModels::granularPressureModel::New
85  (
86  coeffDict_
87  )
88  ),
89  frictionalStressModel_
90  (
91  kineticTheoryModels::frictionalStressModel::New
92  (
93  coeffDict_
94  )
95  ),
96 
97  equilibrium_(coeffDict_.lookup("equilibrium")),
98  e_("e", dimless, coeffDict_),
99  alphaMax_("alphaMax", dimless, coeffDict_),
100  alphaMinFriction_
101  (
102  "alphaMinFriction",
103  dimless,
104  coeffDict_
105  ),
106  residualAlpha_
107  (
108  "residualAlpha",
109  dimless,
110  coeffDict_
111  ),
112 
113  Theta_
114  (
115  IOobject
116  (
117  IOobject::groupName("Theta", phase.name()),
118  U.time().timeName(),
119  U.mesh(),
120  IOobject::MUST_READ,
121  IOobject::AUTO_WRITE
122  ),
123  U.mesh()
124  ),
125 
126  lambda_
127  (
128  IOobject
129  (
130  IOobject::groupName("lambda", phase.name()),
131  U.time().timeName(),
132  U.mesh(),
133  IOobject::NO_READ,
134  IOobject::NO_WRITE
135  ),
136  U.mesh(),
137  dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
138  ),
139 
140  gs0_
141  (
142  IOobject
143  (
144  IOobject::groupName("gs0", phase.name()),
145  U.time().timeName(),
146  U.mesh(),
147  IOobject::NO_READ,
148  IOobject::NO_WRITE
149  ),
150  U.mesh(),
151  dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 0.0)
152  ),
153 
154  kappa_
155  (
156  IOobject
157  (
158  IOobject::groupName("kappa", phase.name()),
159  U.time().timeName(),
160  U.mesh(),
161  IOobject::NO_READ,
162  IOobject::NO_WRITE
163  ),
164  U.mesh(),
165  dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
166  )
167 {
168  if (type == typeName)
169  {
170  printCoeffs(type);
171  }
172 }
173 
174 
175 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
176 
178 {}
179 
180 
181 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
182 
184 {
185  if
186  (
187  eddyViscosity
188  <
189  RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel> >
190  >::read()
191  )
192  {
193  coeffDict().lookup("equilibrium") >> equilibrium_;
194  e_.readIfPresent(coeffDict());
195  alphaMax_.readIfPresent(coeffDict());
196  alphaMinFriction_.readIfPresent(coeffDict());
197 
198  viscosityModel_->read();
199  conductivityModel_->read();
200  radialModel_->read();
201  granularPressureModel_->read();
202  frictionalStressModel_->read();
203 
204  return true;
205  }
206  else
207  {
208  return false;
209  }
210 }
211 
212 
215 {
216  notImplemented("kineticTheoryModel::k()");
217  return nut_;
218 }
219 
220 
223 {
224  notImplemented("kineticTheoryModel::epsilon()");
225  return nut_;
226 }
227 
228 
231 {
232  return tmp<volSymmTensorField>
233  (
235  (
236  IOobject
237  (
238  IOobject::groupName("R", U_.group()),
239  runTime_.timeName(),
240  mesh_,
243  ),
244  - (nut_)*dev(twoSymm(fvc::grad(U_)))
245  - (lambda_*fvc::div(phi_))*symmTensor::I
246  )
247  );
248 }
249 
250 
253 {
254  const volScalarField& rho = phase_.rho();
255 
256  tmp<volScalarField> tpPrime
257  (
258  Theta_
259  *granularPressureModel_->granularPressureCoeffPrime
260  (
261  alpha_,
262  radialModel_->g0(alpha_, alphaMinFriction_, alphaMax_),
263  radialModel_->g0prime(alpha_, alphaMinFriction_, alphaMax_),
264  rho,
265  e_
266  )
267  + frictionalStressModel_->frictionalPressurePrime
268  (
269  alpha_,
270  alphaMinFriction_,
271  alphaMax_
272  )
273  );
274 
275  volScalarField::GeometricBoundaryField& bpPrime = tpPrime().boundaryField();
276 
277  forAll(bpPrime, patchi)
278  {
279  if (!bpPrime[patchi].coupled())
280  {
281  bpPrime[patchi] == 0;
282  }
283  }
284 
285  return tpPrime;
286 }
287 
288 
291 {
292  return fvc::interpolate(pPrime());
293 }
294 
295 
298 {
299  return tmp<volSymmTensorField>
300  (
302  (
303  IOobject
304  (
305  IOobject::groupName("devRhoReff", U_.group()),
306  runTime_.timeName(),
307  mesh_,
310  ),
311  - (rho_*nut_)
312  *dev(twoSymm(fvc::grad(U_)))
313  - ((rho_*lambda_)*fvc::div(phi_))*symmTensor::I
314  )
315  );
316 }
317 
318 
321 (
322  volVectorField& U
323 ) const
324 {
325  return
326  (
327  - fvm::laplacian(rho_*nut_, U)
328  - fvc::div
329  (
330  (rho_*nut_)*dev2(T(fvc::grad(U)))
331  + ((rho_*lambda_)*fvc::div(phi_))
332  *dimensioned<symmTensor>("I", dimless, symmTensor::I)
333  )
334  );
335 }
336 
337 
339 {
340  // Local references
341  volScalarField alpha(max(alpha_, scalar(0)));
342  const volScalarField& rho = phase_.rho();
343  const surfaceScalarField& alphaRhoPhi = alphaRhoPhi_;
344  const volVectorField& U = U_;
345  const volVectorField& Uc_ =
346  refCast<const twoPhaseSystem>(phase_.fluid()).otherPhase(phase_).U();
347 
348  const scalar sqrtPi = sqrt(constant::mathematical::pi);
349  dimensionedScalar ThetaSmall("ThetaSmall", Theta_.dimensions(), 1.0e-6);
350  dimensionedScalar ThetaSmallSqrt(sqrt(ThetaSmall));
351 
352  tmp<volScalarField> tda(phase_.d());
353  const volScalarField& da = tda();
354 
355  tmp<volTensorField> tgradU(fvc::grad(U_));
356  const volTensorField& gradU(tgradU());
357  volSymmTensorField D(symm(gradU));
358 
359  // Calculating the radial distribution function
360  gs0_ = radialModel_->g0(alpha, alphaMinFriction_, alphaMax_);
361 
362  if (!equilibrium_)
363  {
364  // Particle viscosity (Table 3.2, p.47)
365  nut_ = viscosityModel_->nu(alpha, Theta_, gs0_, rho, da, e_);
366 
367  volScalarField ThetaSqrt("sqrtTheta", sqrt(Theta_));
368 
369  // Bulk viscosity p. 45 (Lun et al. 1984).
370  lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;
371 
372  // Stress tensor, Definitions, Table 3.1, p. 43
374  (
375  rho*(2.0*nut_*D + (lambda_ - (2.0/3.0)*nut_)*tr(D)*I)
376  );
377 
378  // Dissipation (Eq. 3.24, p.50)
379  volScalarField gammaCoeff
380  (
381  "gammaCoeff",
382  12.0*(1.0 - sqr(e_))
383  *max(sqr(alpha), residualAlpha_)
384  *rho*gs0_*(1.0/da)*ThetaSqrt/sqrtPi
385  );
386 
387  // Drag
389  (
390  refCast<const twoPhaseSystem>(phase_.fluid()).drag(phase_).K()
391  );
392 
393  // Eq. 3.25, p. 50 Js = J1 - J2
394  volScalarField J1("J1", 3.0*beta);
395  volScalarField J2
396  (
397  "J2",
398  0.25*sqr(beta)*da*magSqr(U - Uc_)
399  /(
400  max(alpha, residualAlpha_)*rho
401  *sqrtPi*(ThetaSqrt + ThetaSmallSqrt)
402  )
403  );
404 
405  // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
406  volScalarField PsCoeff
407  (
408  granularPressureModel_->granularPressureCoeff
409  (
410  alpha,
411  gs0_,
412  rho,
413  e_
414  )
415  );
416 
417  // 'thermal' conductivity (Table 3.3, p. 49)
418  kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
419 
420  // Construct the granular temperature equation (Eq. 3.20, p. 44)
421  // NB. note that there are two typos in Eq. 3.20:
422  // Ps should be without grad
423  // the laplacian has the wrong sign
424  fvScalarMatrix ThetaEqn
425  (
426  1.5*
427  (
428  fvm::ddt(alpha, rho, Theta_)
429  + fvm::div(alphaRhoPhi, Theta_)
430  - fvc::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaRhoPhi), Theta_)
431  )
432  - fvm::laplacian(kappa_, Theta_, "laplacian(kappa,Theta)")
433  ==
434  fvm::SuSp(-((PsCoeff*I) && gradU), Theta_)
435  + (tau && gradU)
436  + fvm::Sp(-gammaCoeff, Theta_)
437  + fvm::Sp(-J1, Theta_)
438  + fvm::Sp(J2/(Theta_ + ThetaSmall), Theta_)
439  );
440 
441  ThetaEqn.relax();
442  ThetaEqn.solve();
443  }
444  else
445  {
446  // Equilibrium => dissipation == production
447  // Eq. 4.14, p.82
448  volScalarField K1("K1", 2.0*(1.0 + e_)*rho*gs0_);
450  (
451  "K3",
452  0.5*da*rho*
453  (
454  (sqrtPi/(3.0*(3.0 - e_)))
455  *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha*gs0_)
456  +1.6*alpha*gs0_*(1.0 + e_)/sqrtPi
457  )
458  );
459 
461  (
462  "K2",
463  4.0*da*rho*(1.0 + e_)*alpha*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0
464  );
465 
466  volScalarField K4("K4", 12.0*(1.0 - sqr(e_))*rho*gs0_/(da*sqrtPi));
467 
468  volScalarField trD
469  (
470  "trD",
471  alpha/(alpha + residualAlpha_)
472  *fvc::div(phi_)
473  );
474  volScalarField tr2D("tr2D", sqr(trD));
475  volScalarField trD2("trD2", tr(D & D));
476 
477  volScalarField t1("t1", K1*alpha + rho);
478  volScalarField l1("l1", -t1*trD);
479  volScalarField l2("l2", sqr(t1)*tr2D);
480  volScalarField l3
481  (
482  "l3",
483  4.0
484  *K4
485  *alpha
486  *(2.0*K3*trD2 + K2*tr2D)
487  );
488 
489  Theta_ = sqr
490  (
491  (l1 + sqrt(l2 + l3))
492  /(2.0*max(alpha, residualAlpha_)*K4)
493  );
494 
495  kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
496  }
497 
498  Theta_.max(0);
499  Theta_.min(100);
500 
501  {
502  // particle viscosity (Table 3.2, p.47)
503  nut_ = viscosityModel_->nu(alpha, Theta_, gs0_, rho, da, e_);
504 
505  volScalarField ThetaSqrt("sqrtTheta", sqrt(Theta_));
506 
507  // Bulk viscosity p. 45 (Lun et al. 1984).
508  lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;
509 
510  // Frictional pressure
511  volScalarField pf
512  (
513  frictionalStressModel_->frictionalPressure
514  (
515  alpha,
516  alphaMinFriction_,
517  alphaMax_
518  )
519  );
520 
521  // Add frictional shear viscosity, Eq. 3.30, p. 52
522  nut_ += frictionalStressModel_->nu
523  (
524  alpha,
525  alphaMinFriction_,
526  alphaMax_,
527  pf/rho,
528  D
529  );
530 
531  // Limit viscosity
532  nut_.min(100);
533  }
534 
535  if (debug)
536  {
537  Info<< typeName << ':' << nl
538  << " max(Theta) = " << max(Theta_).value() << nl
539  << " max(nut) = " << max(nut_).value() << endl;
540  }
541 }
542 
543 
544 // ************************************************************************* //
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.
ThermalDiffusivity< PhaseCompressibleTurbulenceModel< phaseModel > > phaseCompressibleTurbulenceModel
Typedef for phaseCompressibleTurbulenceModel.
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