All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MovingPhaseModel.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) 2015-2021 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 "MovingPhaseModel.H"
27 #include "phaseSystem.H"
29 #include "slipFvPatchFields.H"
31 
32 #include "fvmDdt.H"
33 #include "fvmDiv.H"
34 #include "fvmSup.H"
35 #include "fvcDdt.H"
36 #include "fvcDiv.H"
37 #include "fvcFlux.H"
38 
39 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
40 
41 template<class BasePhaseModel>
44 {
45  word phiName(IOobject::groupName("phi", this->name()));
46 
47  IOobject phiHeader
48  (
49  phiName,
50  U.mesh().time().timeName(),
51  U.mesh(),
52  IOobject::NO_READ
53  );
54 
55  if (phiHeader.typeHeaderOk<surfaceScalarField>(true))
56  {
57  Info<< "Reading face flux field " << phiName << endl;
58 
59  return tmp<surfaceScalarField>
60  (
62  (
63  IOobject
64  (
65  phiName,
66  U.mesh().time().timeName(),
67  U.mesh(),
68  IOobject::MUST_READ,
69  IOobject::AUTO_WRITE
70  ),
71  U.mesh()
72  )
73  );
74  }
75  else
76  {
77  Info<< "Calculating face flux field " << phiName << endl;
78 
79  wordList phiTypes
80  (
81  U.boundaryField().size(),
82  calculatedFvPatchScalarField::typeName
83  );
84 
85  forAll(U.boundaryField(), patchi)
86  {
87  if (!U.boundaryField()[patchi].assignable())
88  {
89  phiTypes[patchi] = fixedValueFvPatchScalarField::typeName;
90  }
91  }
92 
93  return tmp<surfaceScalarField>
94  (
96  (
97  IOobject
98  (
99  phiName,
100  U.mesh().time().timeName(),
101  U.mesh(),
102  IOobject::NO_READ,
103  IOobject::AUTO_WRITE
104  ),
105  fvc::flux(U),
106  phiTypes
107  )
108  );
109  }
110 }
111 
112 
113 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
114 
115 template<class BasePhaseModel>
117 (
118  const phaseSystem& fluid,
119  const word& phaseName,
120  const bool referencePhase,
121  const label index
122 )
123 :
124  BasePhaseModel(fluid, phaseName, referencePhase, index),
125  U_
126  (
127  IOobject
128  (
129  IOobject::groupName("U", this->name()),
130  fluid.mesh().time().timeName(),
131  fluid.mesh(),
132  IOobject::MUST_READ,
133  IOobject::AUTO_WRITE
134  ),
135  fluid.mesh()
136  ),
137  phi_(phi(U_)),
138  alphaPhi_
139  (
140  IOobject
141  (
142  IOobject::groupName("alphaPhi", this->name()),
143  fluid.mesh().time().timeName(),
144  fluid.mesh()
145  ),
146  fluid.mesh(),
147  dimensionedScalar(dimensionSet(0, 3, -1, 0, 0), 0)
148  ),
149  alphaRhoPhi_
150  (
151  IOobject
152  (
153  IOobject::groupName("alphaRhoPhi", this->name()),
154  fluid.mesh().time().timeName(),
155  fluid.mesh()
156  ),
157  fluid.mesh(),
158  dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), 0)
159  ),
160  Uf_(nullptr),
161  DUDt_(nullptr),
162  DUDtf_(nullptr),
163  divU_(nullptr),
164  turbulence_
165  (
166  phaseCompressible::momentumTransportModel::New
167  (
168  *this,
169  this->thermo().rho(),
170  U_,
171  alphaRhoPhi_,
172  phi_,
173  *this
174  )
175  ),
176  thermophysicalTransport_
177  (
178  PhaseThermophysicalTransportModel
179  <
180  phaseCompressible::momentumTransportModel,
181  transportThermoModel
182  >::New(turbulence_, this->thermo_)
183  ),
184  continuityError_
185  (
186  IOobject
187  (
188  IOobject::groupName("continuityError", this->name()),
189  fluid.mesh().time().timeName(),
190  fluid.mesh()
191  ),
192  fluid.mesh(),
194  ),
195  K_(nullptr)
196 {
197  phi_.writeOpt() = IOobject::AUTO_WRITE;
198 
199  if (fluid.mesh().dynamic())
200  {
201  Uf_ = new surfaceVectorField
202  (
203  IOobject
204  (
205  IOobject::groupName("Uf", this->name()),
206  fluid.mesh().time().timeName(),
207  fluid.mesh(),
210  ),
211  fvc::interpolate(U_)
212  );
213  }
214 
215  correctKinematics();
216 }
217 
218 
219 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
220 
221 template<class BasePhaseModel>
223 {}
224 
225 
226 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
227 
228 template<class BasePhaseModel>
230 (
231  const volScalarField& source
232 )
233 {
234  volScalarField& rho = this->thermoRef().rho();
235 
236  continuityError_ = fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_) - source;
237 }
238 
239 
240 template<class BasePhaseModel>
242 {
244  this->fluid().MRF().correctBoundaryVelocity(U_);
245 }
246 
247 
248 template<class BasePhaseModel>
250 {
251  BasePhaseModel::correctKinematics();
252 
253  if (DUDt_.valid())
254  {
255  DUDt_.clear();
256  DUDt();
257  }
258 
259  if (DUDtf_.valid())
260  {
261  DUDtf_.clear();
262  DUDtf();
263  }
264 
265  if (K_.valid())
266  {
267  K_.ref() = 0.5*magSqr(this->U());
268  }
269 }
270 
271 
272 template<class BasePhaseModel>
274 {
275  BasePhaseModel::correctTurbulence();
276 
277  turbulence_->correct();
278 }
279 
280 
281 template<class BasePhaseModel>
283 {
284  BasePhaseModel::correctEnergyTransport();
285  thermophysicalTransport_->correct();
286 }
287 
288 
289 template<class BasePhaseModel>
291 {
292  const fvMesh& mesh = this->fluid().mesh();
293 
294  if (mesh.dynamic())
295  {
296  Uf_.ref() = fvc::interpolate(U_);
297  surfaceVectorField n(mesh.Sf()/mesh.magSf());
298  Uf_.ref() += n*(fvc::absolute(phi_, U_)/mesh.magSf() - (n & Uf_()));
299 
300  surfaceVectorField::Boundary& UfBf = Uf_.ref().boundaryFieldRef();
301  const volVectorField::Boundary& UBf = U_.boundaryField();
302 
303  forAll(mesh.boundary(), patchi)
304  {
305  // Remove the flux correction on AMI patches to compensate for
306  // AMI non-conservation error
307  if (isA<cyclicAMIFvPatch>(mesh.boundary()[patchi]))
308  {
309  UfBf[patchi] = UBf[patchi];
310  }
311  }
312  }
313 }
314 
315 
316 template<class BasePhaseModel>
318 {
319  return false;
320 }
321 
322 
323 template<class BasePhaseModel>
326 {
327  const volScalarField& alpha = *this;
328  const volScalarField& rho = this->thermo().rho();
329 
330  return
331  (
332  fvm::ddt(alpha, rho, U_)
333  + fvm::div(alphaRhoPhi_, U_)
334  + fvm::SuSp(-this->continuityError(), U_)
335  + this->fluid().MRF().DDt(alpha*rho, U_)
336  + turbulence_->divDevTau(U_)
337  );
338 }
339 
340 
341 template<class BasePhaseModel>
344 {
345  // As the "normal" U-eqn but without the ddt terms
346 
347  const volScalarField& alpha = *this;
348  const volScalarField& rho = this->thermo().rho();
349 
350  return
351  (
352  fvm::div(alphaRhoPhi_, U_)
353  + fvm::SuSp(fvc::ddt(*this, rho) - this->continuityError(), U_)
354  + this->fluid().MRF().DDt(alpha*rho, U_)
355  + turbulence_->divDevTau(U_)
356  );
357 }
358 
359 
360 template<class BasePhaseModel>
363 {
364  return U_;
365 }
366 
367 
368 template<class BasePhaseModel>
371 {
372  return U_;
373 }
374 
375 
376 template<class BasePhaseModel>
379 {
380  return phi_;
381 }
382 
383 
384 template<class BasePhaseModel>
387 {
388  return phi_;
389 }
390 
391 
392 template<class BasePhaseModel>
395 {
396  return
397  Uf_.valid()
398  ? tmp<surfaceVectorField>(Uf_())
399  : tmp<surfaceVectorField>();
400 }
401 
402 
403 template<class BasePhaseModel>
406 {
407  if (Uf_.valid())
408  {
409  return Uf_.ref();
410  }
411  else
412  {
414  << "Uf has not been allocated."
415  << exit(FatalError);
416 
417  return const_cast<surfaceVectorField&>(surfaceVectorField::null());
418  }
419 }
420 
421 
422 template<class BasePhaseModel>
425 {
426  return alphaPhi_;
427 }
428 
429 
430 template<class BasePhaseModel>
433 {
434  return alphaPhi_;
435 }
436 
437 
438 template<class BasePhaseModel>
441 {
442  return alphaRhoPhi_;
443 }
444 
445 
446 template<class BasePhaseModel>
449 {
450  return alphaRhoPhi_;
451 }
452 
453 
454 template<class BasePhaseModel>
457 {
458  if (!DUDt_.valid())
459  {
460  const tmp<surfaceScalarField> taphi(fvc::absolute(phi_, U_));
461  const surfaceScalarField& aphi(taphi());
462  DUDt_ = fvc::ddt(U_) + fvc::div(aphi, U_) - fvc::div(aphi)*U_;
463  }
464 
465  return tmp<volVectorField>(DUDt_());
466 }
467 
468 
469 template<class BasePhaseModel>
472 {
473  if (!DUDtf_.valid())
474  {
475  DUDtf_ = byDt(phi_ - phi_.oldTime());
476  }
477 
478  return tmp<surfaceScalarField>(DUDtf_());
479 }
480 
481 
482 template<class BasePhaseModel>
485 {
486  return continuityError_;
487 }
488 
489 
490 template<class BasePhaseModel>
493 {
494  if (!K_.valid())
495  {
497  (
498  IOobject::groupName("K", this->name()),
499  0.5*magSqr(this->U())
500  );
501  }
502 
503  return tmp<volScalarField>(K_());
504 }
505 
506 
507 template<class BasePhaseModel>
510 {
511  return divU_.valid() ? tmp<volScalarField>(divU_()) : tmp<volScalarField>();
512 }
513 
514 
515 template<class BasePhaseModel>
516 void Foam::MovingPhaseModel<BasePhaseModel>::divU(tmp<volScalarField> divU)
517 {
518  divU_ = divU;
519 }
520 
521 
522 template<class BasePhaseModel>
525 {
526  return turbulence_->k();
527 }
528 
529 
530 template<class BasePhaseModel>
533 {
534  return turbulence_->pPrime();
535 }
536 
537 
538 template<class BasePhaseModel>
541 {
542  return thermophysicalTransport_->kappaEff(patchi);
543 }
544 
545 
546 template<class BasePhaseModel>
549 {
550  return thermophysicalTransport_->divq(he);
551 }
552 
553 
554 template<class BasePhaseModel>
557 {
558  return thermophysicalTransport_->divj(Yi);
559 }
560 
561 
562 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual tmp< volScalarField > divU() const
Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
tmp< GeometricField< Type, fvPatchField, volMesh > > DDt(const surfaceScalarField &phi, const GeometricField< Type, fvPatchField, volMesh > &psi)
Definition: fvcDDt.C:45
tmp< fvMatrix< Type > > SuSp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
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
fluidReactionThermo & thermo
Definition: createFields.H:28
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
virtual surfaceScalarField & alphaRhoPhiRef()
Access the mass flux of the phase.
IOMRFZoneList & MRF
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
virtual tmp< surfaceVectorField > Uf() const
Return the face velocity.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
virtual bool stationary() const
Return whether the phase is stationary.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual tmp< volVectorField > DUDt() const
Return the substantive acceleration.
virtual tmp< surfaceScalarField > alphaPhi() const
Return the volumetric flux of the phase.
virtual tmp< fvScalarMatrix > divj(volScalarField &Yi) const
Return the source term for the given specie mass-fraction.
virtual void correctTurbulence()
Correct the turbulence.
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:58
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
virtual tmp< volScalarField > k() const
Return the turbulent kinetic energy.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
Calculate the first temporal derivative.
CompressibleMomentumTransportModel< dynamicTransportModel > momentumTransportModel
const dimensionSet dimTime
MovingPhaseModel(const phaseSystem &fluid, const word &phaseName, const bool referencePhase, const label index)
virtual tmp< surfaceScalarField > alphaRhoPhi() const
Return the mass flux of the phase.
autoPtr< BasicCompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleMomentumTransportModel::transportModel &transport)
virtual tmp< fvVectorMatrix > UfEqn()
Return the momentum equation for the face-based algorithm.
virtual ~MovingPhaseModel()
Destructor.
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-9/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:68
Calculate the face-flux of the given field.
static word groupName(Name name, const word &group)
const dimensionSet dimDensity
virtual surfaceScalarField & phiRef()
Access the volumetric flux.
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
Calculate the matrix for the first temporal derivative.
phaseSystem & fluid
Definition: createFields.H:11
phi
Definition: correctPhi.H:3
word timeName
Definition: getTimeIndex.H:3
virtual surfaceVectorField & UfRef()
Access the face velocity.
virtual volVectorField & URef()
Access the velocity.
virtual tmp< volScalarField > pPrime() const
Return the phase-pressure&#39;.
Calculate the divergence of the given field.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
virtual void correctUf()
Correct the face velocity for moving meshes.
tmp< volScalarField > divU
Definition: alphaSuSp.H:9
virtual void correct()
Correct the phase properties other than the thermo and turbulence.
virtual tmp< volVectorField > U() const
Return the velocity.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
virtual void correctEnergyTransport()
Correct the energy transport e.g. alphat.
tmp< volScalarField > byDt(const volScalarField &vf)
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.
static const GeometricField< vector, fvsPatchField, surfaceMesh > & null()
Return a null geometric field.
virtual tmp< volScalarField > continuityError() const
Return the continuity error.
virtual tmp< fvScalarMatrix > divq(volScalarField &he) const
Return the source term for the energy equation.
Calculate the matrix for the divergence of the given field and flux.
List< word > wordList
A List of words.
Definition: fileName.H:54
virtual tmp< volScalarField > K() const
Return the phase kinetic energy.
label patchi
virtual void correctContinuityError(const volScalarField &source)
Correct the continuity error.
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:188
U
Definition: pEqn.H:72
fvModels source(alpha1, mixture.thermo1().rho())
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
virtual surfaceScalarField & alphaPhiRef()
Access the volumetric flux of the phase.
virtual tmp< fvVectorMatrix > UEqn()
Return the momentum equation.
virtual tmp< scalarField > kappaEff(const label patchi) const
Return the effective thermal conductivity on a patch.
messageStream Info
label n
virtual void correctKinematics()
Correct the kinematics.
virtual tmp< surfaceScalarField > phi() const
Return the volumetric flux.
A class for managing temporary objects.
Definition: PtrList.H:53
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
virtual tmp< surfaceScalarField > DUDtf() const
Return the substantive acceleration on the faces.
Calculate the matrix for implicit and explicit sources.