phaseScalarTransport.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) 2019 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 
27 #include "fixedValueFvPatchField.H"
28 #include "fvcDdt.H"
29 #include "fvcDiv.H"
30 #include "fvmDdt.H"
31 #include "fvmDiv.H"
32 #include "fvmLaplacian.H"
34 #include "phaseScalarTransport.H"
35 #include "surfaceFields.H"
36 #include "turbulenceModel.H"
37 #include "wallFvPatch.H"
39 
40 #define PhiDimensionErrorInFunction(phi) \
41  FatalErrorInFunction \
42  << "Incompatible dimensions for " << phi.name() << ": " \
43  << phi.dimensions() << nl \
44  << "Dimensions should be " << dimMass/dimTime << " or " \
45  << dimVolume/dimTime << exit(FatalError)
46 
47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
48 
49 namespace Foam
50 {
51 namespace functionObjects
52 {
53  defineTypeNameAndDebug(phaseScalarTransport, 0);
54 
56  (
57  functionObject,
58  phaseScalarTransport,
59  dictionary
60  );
61 }
62 }
63 
64 
65 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
66 
67 Foam::volScalarField& Foam::functionObjects::phaseScalarTransport::Phi()
68 {
69  if (!PhiPtr_.valid())
70  {
71  const surfaceScalarField& phi =
72  mesh_.lookupObject<surfaceScalarField>(phiName_);
73  const volScalarField& p =
74  mesh_.lookupObject<volScalarField>(pName_);
75 
76  wordList PhiPatchFieldTypes(mesh_.boundaryMesh().size());
77  forAll(p.boundaryField(), patchi)
78  {
79  PhiPatchFieldTypes[patchi] =
80  p.boundaryField()[patchi].fixesValue()
83  }
84 
85  PhiPtr_.set
86  (
87  new volScalarField
88  (
89  IOobject
90  (
91  "Phi" + s_.name(),
92  mesh_.time().timeName(),
93  mesh_,
96  ),
97  mesh_,
98  dimensionedScalar(phi.dimensions()/dimLength, Zero),
99  PhiPatchFieldTypes
100  )
101  );
102 
103  mesh_.setFluxRequired(PhiPtr_->name());
104  }
105 
106  return PhiPtr_();
107 }
108 
109 
111 Foam::functionObjects::phaseScalarTransport::alphaPhi()
112 {
113  // If alphaPhi exists then return it
114  if (mesh_.foundObject<surfaceScalarField>(alphaPhiName_))
115  {
116  return mesh_.lookupObject<surfaceScalarField>(alphaPhiName_);
117  }
118 
119  // Otherwise generate it ...
120  Info<< type() << ": " << surfaceScalarField::typeName << " "
121  << alphaPhiName_ << " was not found, so generating it" << endl;
122 
123  const volScalarField& alpha =
124  mesh_.lookupObject<volScalarField>(alphaName_);
125  const surfaceScalarField& phi =
126  mesh_.lookupObject<surfaceScalarField>(phiName_);
127 
128  // Make a crude guess of the phase flux using default interpolation
129  tmp<surfaceScalarField> tAlphaPhi
130  (
132  (
133  alphaPhiName_,
134  phi*fvc::interpolate(alpha)
135  )
136  );
137  surfaceScalarField& alphaPhi = tAlphaPhi.ref();
138 
139  // Get the potential field
140  volScalarField& Phi(this->Phi());
141 
142  // Construct the scheme names
143  const word laplacianScheme = "laplacian(" + pName_ + ")";
144 
145  // Debug writing. Write the material derivative of alpha, before and after
146  // the solution of the potential and the correction of alphaPhi. Before
147  // correction the field should be non-zero, and after it should be
148  // comparable to the solution tolerance.
149  auto writeDDt = [&](const label i)
150  {
151  const volScalarField DDtAlpha
152  (
153  "DDt("
155  (
156  IOobject::member(alpha.name()) + Foam::name(i),
157  IOobject::group(alpha.name())
158  )
159  + ")",
160  fvc::ddt(alpha) + fvc::div(alphaPhi)
161  );
162  Info<< type() << ": Writing " << DDtAlpha.name() << endl;
163  DDtAlpha.write();
164  };
165  if (debug && mesh_.time().writeTime())
166  {
167  writeDDt(0);
168  }
169 
170  // Lookup the non-orthogonal solution control
171  nonOrthogonalSolutionControl& control =
172  mesh_.lookupObjectRef<nonOrthogonalSolutionControl>
173  (
174  solutionControl::typeName
175  );
176 
177  // Solve for the potential and correct alphaPhi with the resulting flux
178  if (phi.dimensions() == dimVolume/dimTime)
179  {
180  while (control.correctNonOrthogonal())
181  {
182  fvScalarMatrix PhiEqn
183  (
184  fvm::laplacian(Phi, laplacianScheme)
185  + fvc::ddt(alpha)
186  + fvc::div(alphaPhi)
187  );
188 
189  PhiEqn.solve(pName_);
190 
191  if (control.finalNonOrthogonalIter())
192  {
193  alphaPhi += PhiEqn.flux();
194  }
195  }
196  }
197  else if (phi.dimensions() == dimMass/dimTime)
198  {
199  const volScalarField& rho =
200  mesh_.lookupObject<volScalarField>(rhoName_);
201 
202  while (control.correctNonOrthogonal())
203  {
204  fvScalarMatrix PhiEqn
205  (
206  fvm::laplacian(Phi, laplacianScheme)
207  + fvc::ddt(rho, alpha)
208  + fvc::div(alphaPhi)
209  );
210 
211  PhiEqn.solve(pName_);
212 
213  if (control.finalNonOrthogonalIter())
214  {
215  alphaPhi += PhiEqn.flux();
216  }
217  }
218  }
219  else
220  {
222  }
223 
224  // Debug writing
225  if (debug && mesh_.time().writeTime())
226  {
227  writeDDt(1);
228  }
229 
230  return tAlphaPhi;
231 }
232 
233 
235 Foam::functionObjects::phaseScalarTransport::D
236 (
238 ) const
239 {
240  if (constantD_)
241  {
242  return volScalarField::New
243  (
244  "D" + s_.name(),
245  mesh_,
246  dimensionedScalar(alphaPhi.dimensions()/dimLength, D_)
247  );
248  }
249 
250  const word& nameNoPhase = turbulenceModel::propertiesName;
251  const word namePhase = IOobject::groupName(nameNoPhase, phaseName_);
252 
253  const word& name =
254  mesh_.foundObject<turbulenceModel>(namePhase)
255  ? namePhase
256  : mesh_.foundObject<turbulenceModel>(nameNoPhase)
257  ? nameNoPhase
258  : word::null;
259 
260  if (name == word::null)
261  {
262  return volScalarField::New
263  (
264  "D" + s_.name(),
265  mesh_,
266  dimensionedScalar(alphaPhi.dimensions()/dimLength, 0)
267  );
268  }
269 
270  const turbulenceModel& turbulence =
271  mesh_.lookupObject<turbulenceModel>(name);
272 
273  if (alphaPhi.dimensions() == dimVolume/dimTime)
274  {
275  return alphaD_*turbulence.nu() + alphaDt_*turbulence.nut();
276  }
277  else if (alphaPhi.dimensions() == dimMass/dimTime)
278  {
279  return alphaD_*turbulence.mu() + alphaDt_*turbulence.mut();
280  }
281  else
282  {
283  PhiDimensionErrorInFunction(alphaPhi);
284  return tmp<volScalarField>(nullptr);
285  }
286 }
287 
288 
289 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
290 
292 (
293  const word& name,
294  const Time& runTime,
295  const dictionary& dict
296 )
297 :
298  fvMeshFunctionObject(name, runTime, dict),
299  fieldName_(dict.lookup("field")),
300  phaseName_(IOobject::group(fieldName_)),
301  nCorr_(0),
302  residualAlpha_(rootSmall),
303  fvOptions_(mesh_),
304  s_
305  (
306  IOobject
307  (
308  fieldName_,
309  mesh_.time().timeName(),
310  mesh_,
313  ),
314  mesh_
315  ),
316  alphaSPtr_(nullptr),
317  PhiPtr_(nullptr)
318 {
319  if (phaseName_ == word::null)
320  {
322  << "Field \"" << fieldName_ << "\" does not have a phase extension "
323  << "in its name. If it is associated with \"phaseA\" then it "
324  << "should be named \"" << fieldName_ << ".phaseA\"."
325  << exit(FatalError);
326  }
327 
328  read(dict);
329 }
330 
331 
332 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
333 
335 {}
336 
337 
338 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
339 
341 {
343 
344  alphaName_ =
345  dict.lookupOrDefault<word>
346  (
347  "alpha",
348  IOobject::groupName("alpha", phaseName_)
349  );
350  alphaPhiName_ =
351  dict.lookupOrDefault<word>
352  (
353  "alphaPhi",
354  IOobject::groupName("alphaPhi", phaseName_)
355  );
356  phiName_ = dict.lookupOrDefault<word>("phi", "phi");
357  rhoName_ =
358  dict.lookupOrDefault<word>
359  (
360  "rho",
361  IOobject::groupName("rho", phaseName_)
362  );
363  pName_ = dict.lookupOrDefault<word>("p", "p");
364  schemesField_ = dict.lookupOrDefault<word>("schemesField", fieldName_);
365 
366  constantD_ = dict.readIfPresent("D", D_);
367  alphaD_ = dict.lookupOrDefault<scalar>("alphaD", 1);
368  alphaDt_ = dict.lookupOrDefault<scalar>("alphaDt", 1);
369 
370  dict.readIfPresent("nCorr", nCorr_);
371  dict.readIfPresent("residualAlpha", residualAlpha_);
372  writeAlphaField_ = dict.lookupOrDefault<bool>("writeAlphaField", true);
373 
374  if (dict.found("fvOptions"))
375  {
376  fvOptions_.reset(dict.subDict("fvOptions"));
377  }
378 
379  return true;
380 }
381 
382 
384 {
385  Info<< type() << ": Executing" << endl;
386 
387  const volScalarField& alpha =
388  mesh_.lookupObject<volScalarField>(alphaName_);
389 
390  // Get the phase flux
391  tmp<surfaceScalarField> tAlphaPhi(this->alphaPhi());
392  const surfaceScalarField& alphaPhi = tAlphaPhi();
393 
394  // Get the diffusivity
395  const volScalarField D(this->D(alphaPhi));
396 
397  // Construct the scheme names
398  const word divScheme =
399  "div(" + alphaPhi.name() + "," + schemesField_ + ")";
400  const word laplacianScheme =
401  "laplacian(" + D.name() + "," + schemesField_ + ")";
402 
403  // Get the relaxation coefficient
404  const scalar relaxCoeff =
405  mesh_.relaxEquation(schemesField_)
406  ? mesh_.equationRelaxationFactor(schemesField_)
407  : 0;
408 
409  // Solve
410  if (alphaPhi.dimensions() == dimVolume/dimTime)
411  {
412  for (label i = 0; i <= nCorr_; ++ i)
413  {
414  fvScalarMatrix fieldEqn
415  (
416  fvm::ddt(alpha, s_)
417  + fvm::div(alphaPhi, s_, divScheme)
419  (
421  s_,
422  laplacianScheme
423  )
424  ==
425  fvOptions_(alpha, s_)
426  - fvm::ddt(residualAlpha_, s_)
427  + fvc::ddt(residualAlpha_, s_)
428  );
429 
430  fieldEqn.relax(relaxCoeff);
431  fvOptions_.constrain(fieldEqn);
432  fieldEqn.solve(schemesField_);
433  }
434  }
435  else if (alphaPhi.dimensions() == dimMass/dimTime)
436  {
437  const volScalarField& rho =
438  mesh_.lookupObject<volScalarField>(rhoName_);
439 
440  for (label i = 0; i <= nCorr_; ++ i)
441  {
442  fvScalarMatrix fieldEqn
443  (
444  fvm::ddt(alpha, rho, s_)
445  + fvm::div(alphaPhi, s_, divScheme)
447  (
449  s_,
450  laplacianScheme
451  )
452  ==
453  fvOptions_(alpha, rho, s_)
454  - fvm::ddt(residualAlpha_*rho, s_)
455  + fvc::ddt(residualAlpha_*rho, s_)
456  );
457 
458  fieldEqn.relax(relaxCoeff);
459  fvOptions_.constrain(fieldEqn);
460  fieldEqn.solve(schemesField_);
461  }
462  }
463  else
464  {
465  PhiDimensionErrorInFunction(alphaPhi);
466  }
467 
468  // Update
469  if (writeAlphaField_)
470  {
471  if (!alphaSPtr_.valid())
472  {
473  alphaSPtr_.set
474  (
475  new volScalarField
476  (
477  IOobject
478  (
479  "alpha"
480  + word(toupper(fieldName_[0]))
481  + fieldName_(1, fieldName_.size() - 1),
482  mesh_.time().timeName(),
483  mesh_,
486  ),
487  mesh_,
488  dimensionedScalar(s_.dimensions(), Zero)
489  )
490  );
491  }
492 
493  alphaSPtr_() = alpha*s_;
494  }
495  else
496  {
497  if (alphaSPtr_.valid())
498  {
499  alphaSPtr_().clear();
500  }
501  }
502 
503  Info<< endl;
504 
505  return true;
506 }
507 
508 
510 {
511  return true;
512 }
513 
514 
515 // ************************************************************************* //
Foam::surfaceFields.
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:438
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const word & name() const
Return name.
Definition: IOobject.H:295
surfaceScalarField & phi
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
static const char *const typeName
Definition: Field.H:105
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
addToRunTimeSelectionTable(functionObject, Qdot, dictionary)
ThermalDiffusivity< CompressibleTurbulenceModel< fluidThermo > > turbulenceModel
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
Calculate the matrix for the laplacian of the field.
word member() const
Return member (name without the extension)
Definition: IOobject.C:378
Info<< "Reading field U\"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);volScalarField rho(IOobject("rho", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), thermo.rho());volVectorField rhoU(IOobject("rhoU", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *U);volScalarField rhoE(IOobject("rhoE", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), rho *(e+0.5 *magSqr(U)));surfaceScalarField pos(IOobject("pos", runTime.timeName(), mesh), mesh, dimensionedScalar(dimless, 1.0));surfaceScalarField neg(IOobject("neg", runTime.timeName(), mesh), mesh, dimensionedScalar(dimless, -1.0));surfaceScalarField phi("phi", fvc::flux(rhoU));Info<< "Creating turbulence model\"<< endl;autoPtr< compressible::turbulenceModel > turbulence(compressible::turbulenceModel::New(rho, U, phi, thermo))
Definition: createFields.H:94
word group() const
Return group (extension part of name)
Definition: IOobject.C:372
virtual bool read(const dictionary &)
Read the settings from the given dictionary.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
Macros for easy insertion into run-time selection tables.
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Mesh &, const dimensionSet &, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return a temporary field constructed from name, mesh, dimensionSet.
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:699
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
Calculate the first temporal derivative.
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
const dimensionSet & dimensions() const
Return dimensions.
virtual bool read(const dictionary &)
Read optional controls.
bool read(const char *, int32_t &)
Definition: int32IO.C:85
virtual bool write()
Do nothing. The field is registered and written automatically.
static const word propertiesName
Default name of the turbulence properties dictionary.
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:124
A class for handling words, derived from string.
Definition: word.H:59
static word groupName(Name name, const word &group)
#define PhiDimensionErrorInFunction(phi)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
Calculate the matrix for the first temporal derivative.
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
static const word null
An empty word.
Definition: word.H:77
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
static const zero Zero
Definition: zero.H:97
Calculate the divergence of the given field.
SolverPerformance< Type > solve(const dictionary &)
Solve segregated or coupled returning the solution statistics.
Definition: fvMatrixSolve.C:58
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:521
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
defineTypeNameAndDebug(Qdot, 0)
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.
Calculate the matrix for the divergence of the given field and flux.
List< word > wordList
A List of words.
Definition: fileName.H:54
label patchi
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
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
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
messageStream Info
virtual Ostream & write(const token &)=0
Write next token to stream.
Specialization of Foam::functionObject for an Foam::fvMesh, providing a reference to the Foam::fvMesh...
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
volScalarField & p
A class for managing temporary objects.
Definition: PtrList.H:53
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
virtual bool execute()
Solve for the evolution of the field.
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
phaseScalarTransport(const word &name, const Time &runTime, const dictionary &dict)
Construct from Time and dictionary.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:583