All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cellPressureCorrector.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) 2022-2024 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 "multiphaseEuler.H"
27 #include "constrainHbyA.H"
28 #include "constrainPressure.H"
29 #include "findRefCell.H"
30 #include "fvcDdt.H"
31 #include "fvcDiv.H"
32 #include "fvcSup.H"
33 #include "fvcSnGrad.H"
34 #include "fvmDdt.H"
35 #include "fvmDiv.H"
36 #include "fvmLaplacian.H"
37 #include "fvmSup.H"
38 #include "fvcFlux.H"
39 #include "fvcMeshPhi.H"
40 #include "fvcReconstruct.H"
41 
42 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
43 
44 void Foam::solvers::multiphaseEuler::cellPressureCorrector()
45 {
47 
48  volScalarField rho("rho", fluid.rho());
49 
50  // Correct p_rgh for consistency with the current density
52 
53  // Face volume fractions
54  PtrList<surfaceScalarField> alphafs(phases.size());
55  forAll(phases, phasei)
56  {
57  const phaseModel& phase = phases[phasei];
58  const volScalarField& alpha = phase;
59 
60  alphafs.set(phasei, fvc::interpolate(max(alpha, scalar(0))).ptr());
61  alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
62  }
63 
64  // Diagonal coefficients
65  rAs.clear();
67  {
68  rAs.setSize(phases.size());
69  }
70 
71  PtrList<volVectorField> HVms(movingPhases.size());
72  PtrList<PtrList<volScalarField>> invADVs;
73  PtrList<PtrList<surfaceScalarField>> invADVfs;
74  {
75  PtrList<volScalarField> As(movingPhases.size());
76  forAll(movingPhases, movingPhasei)
77  {
78  const phaseModel& phase = movingPhases[movingPhasei];
79  const volScalarField& alpha = phase;
80 
81  As.set
82  (
83  movingPhasei,
84  UEqns[phase.index()].A()
85  + byDt
86  (
87  max(phase.residualAlpha() - alpha, scalar(0))
88  *phase.rho()
89  )
90  );
91 
93  {
94  rAs.set
95  (
96  phase.index(),
97  new volScalarField
98  (
99  IOobject::groupName("rA", phase.name()),
100  1/As[movingPhasei]
101  )
102  );
103  }
104  }
105 
106  fluid.invADVs(As, HVms, invADVs, invADVfs);
107  }
108 
109  // Explicit force fluxes
110  PtrList<surfaceScalarField> alphaByADfs;
111  PtrList<surfaceScalarField> FgByADfs;
112  {
113  PtrList<surfaceScalarField> Ffs(fluid.Fs());
114 
115  const surfaceScalarField ghSnGradRho
116  (
117  "ghSnGradRho",
119  );
120 
121  PtrList<surfaceScalarField> lalphafs(movingPhases.size());
122  PtrList<surfaceScalarField> Fgfs(movingPhases.size());
123 
124  forAll(movingPhases, movingPhasei)
125  {
126  const phaseModel& phase = movingPhases[movingPhasei];
127  const volScalarField& alpha = phase;
128 
129  lalphafs.set
130  (
131  movingPhasei,
132  fvc::interpolate(max(alpha, phase.residualAlpha()))
133  );
134 
135  Fgfs.set
136  (
137  movingPhasei,
138  (
139  Ffs[phase.index()]
140  + lalphafs[movingPhasei]
141  *(
142  ghSnGradRho
143  - (fvc::interpolate(phase.rho() - rho))
144  *(buoyancy.g & mesh.Sf())
145  - fluid.surfaceTension(phase)*mesh.magSf()
146  )
147  ).ptr()
148  );
149  }
150 
151  alphaByADfs = invADVfs & lalphafs;
152  FgByADfs = invADVfs & Fgfs;
153  }
154 
155 
156  // Mass transfer rates
157  PtrList<volScalarField> dmdts(fluid.dmdts());
158  PtrList<volScalarField> d2mdtdps(fluid.d2mdtdps());
159 
160  // --- Optional momentum predictor
161  if (predictMomentum)
162  {
163  PtrList<volVectorField> HbyADs;
164  {
165  PtrList<volVectorField> Hs(movingPhases.size());
166 
167  forAll(movingPhases, movingPhasei)
168  {
169  const phaseModel& phase = movingPhases[movingPhasei];
170  const volScalarField& alpha = phase;
171 
172  Hs.set
173  (
174  movingPhasei,
175  UEqns[phase.index()].H()
176  + byDt
177  (
178  max(phase.residualAlpha() - alpha, scalar(0))
179  *phase.rho()
180  )
181  *phase.U()().oldTime()
182  );
183 
184  if (HVms.set(movingPhasei))
185  {
186  Hs[movingPhasei] += HVms[movingPhasei];
187  }
188  }
189 
190  HbyADs = invADVs & Hs;
191  }
192 
193  forAll(movingPhases, movingPhasei)
194  {
195  const phaseModel& phase = movingPhases[movingPhasei];
196  constrainHbyA(HbyADs[movingPhasei], phase.U(), p_rgh);
197  }
198 
199  const surfaceScalarField mSfGradp(-mesh.magSf()*fvc::snGrad(p_rgh));
200 
201  forAll(movingPhases, movingPhasei)
202  {
203  phaseModel& phase = movingPhases_[movingPhasei];
204 
205  phase.URef() =
206  HbyADs[movingPhasei]
208  (
209  alphaByADfs[movingPhasei]*mSfGradp
210  - FgByADfs[movingPhasei]
211  );
212 
213  phase.URef().correctBoundaryConditions();
214  fvConstraints().constrain(phase.URef());
215  }
216  }
217 
218  // --- Pressure corrector loop
219  while (pimple.correct())
220  {
221  // Correct fixed-flux BCs to be consistent with the velocity BCs
223 
224  PtrList<volVectorField> HbyADs;
225  PtrList<surfaceScalarField> phiHbyADs;
226  {
227  // Predicted velocities and fluxes for each phase
228  PtrList<volVectorField> Hs(movingPhases.size());
229  PtrList<surfaceScalarField> phiHs(movingPhases.size());
230 
231  // Correction force fluxes
232  PtrList<surfaceScalarField> ddtCorrs(fluid.ddtCorrs());
233 
234  forAll(movingPhases, movingPhasei)
235  {
236  const phaseModel& phase = movingPhases[movingPhasei];
237  const volScalarField& alpha = phase;
238 
239  Hs.set
240  (
241  movingPhasei,
242  UEqns[phase.index()].H()
243  + byDt
244  (
245  max(phase.residualAlpha() - alpha, scalar(0))
246  *phase.rho()
247  )
248  *phase.U()().oldTime()
249  );
250 
251  if (HVms.set(movingPhasei))
252  {
253  Hs[movingPhasei] += HVms[movingPhasei];
254  }
255 
256  phiHs.set
257  (
258  movingPhasei,
259  fvc::flux(Hs[movingPhasei]) + ddtCorrs[phase.index()]
260  );
261  }
262 
263  HbyADs = invADVs & Hs;
264  phiHbyADs = invADVfs & phiHs;
265  }
266 
267  forAll(movingPhases, movingPhasei)
268  {
269  const phaseModel& phase = movingPhases[movingPhasei];
270 
271  constrainHbyA(HbyADs[movingPhasei], phase.U(), p_rgh);
272  constrainPhiHbyA(phiHbyADs[movingPhasei], phase.U(), p_rgh);
273 
274  phiHbyADs[movingPhasei] -= FgByADfs[movingPhasei];
275  }
276 
277  // Total predicted flux
279  (
280  IOobject
281  (
282  "phiHbyA",
283  runTime.name(),
284  mesh
285  ),
286  mesh,
288  );
289 
290  forAll(movingPhases, movingPhasei)
291  {
292  const phaseModel& phase = movingPhases[movingPhasei];
293 
294  phiHbyA += alphafs[phase.index()]*phiHbyADs[movingPhasei];
295  }
296 
299 
300  // Pressure "diffusivity"
302  (
303  IOobject
304  (
305  "rAf",
306  runTime.name(),
307  mesh
308  ),
309  mesh,
311  );
312 
313  forAll(movingPhases, movingPhasei)
314  {
315  const phaseModel& phase = movingPhases[movingPhasei];
316 
317  rAf += alphafs[phase.index()]*alphaByADfs[movingPhasei];
318  }
319 
320  // Update the fixedFluxPressure BCs to ensure flux consistency
321  {
323  (
326  );
327  phib = 0;
328 
329  forAll(movingPhases, movingPhasei)
330  {
331  phaseModel& phase = movingPhases_[movingPhasei];
332 
333  phib +=
334  alphafs[phase.index()].boundaryField()
335  *phase.phi()().boundaryField();
336  }
337 
338  setSnGrad<fixedFluxPressureFvPatchScalarField>
339  (
341  (
342  phiHbyA.boundaryField() - phib
343  )/(mesh.magSf().boundaryField()*rAf.boundaryField())
344  );
345  }
346 
347  // Compressible pressure equations
348  PtrList<fvScalarMatrix> pEqnComps(compressibilityEqns(dmdts, d2mdtdps));
349 
350  // Cache p prior to solve for density update
351  volScalarField p_rgh_0(p_rgh);
352 
353  // Iterate over the pressure equation to correct for non-orthogonality
354  while (pimple.correctNonOrthogonal())
355  {
356  // Construct the transport part of the pressure equation
357  fvScalarMatrix pEqnIncomp
358  (
360  - fvm::laplacian(rAf, p_rgh)
361  );
362 
363  // Solve
364  {
365  fvScalarMatrix pEqn(pEqnIncomp);
366 
367  forAll(phases, phasei)
368  {
369  pEqn += pEqnComps[phasei];
370  }
371 
372  if (fluid.incompressible())
373  {
374  pEqn.setReference
375  (
378  );
379  }
380 
381  fvConstraints().constrain(pEqn);
382 
383  pEqn.solve();
384  }
385 
386  // Correct fluxes and velocities on last non-orthogonal iteration
388  {
389  phi_ = phiHbyA + pEqnIncomp.flux();
390 
391  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAf);
392 
393  forAll(movingPhases, movingPhasei)
394  {
395  phaseModel& phase = movingPhases_[movingPhasei];
396 
397  phase.phiRef() =
398  phiHbyADs[movingPhasei]
399  + alphaByADfs[movingPhasei]*mSfGradp;
400 
401  // Set the phase dilatation rate
402  phase.divU(-pEqnComps[phase.index()] & p_rgh);
403 
404  MRF.makeRelative(phase.phiRef());
405  fvc::makeRelative(phase.phiRef(), phase.U());
406  }
407 
408  // Optionally relax pressure for velocity correction
409  p_rgh.relax();
410 
411  mSfGradp = pEqnIncomp.flux()/rAf;
412 
413  if (dragCorrection)
414  {
415  PtrList<volVectorField> dragCorrs(movingPhases.size());
416  PtrList<surfaceScalarField> dragCorrfs(movingPhases.size());
417  fluid.dragCorrs(dragCorrs, dragCorrfs);
418 
419  PtrList<volVectorField> dragCorrByADs
420  (
421  invADVs & dragCorrs
422  );
423 
424  PtrList<surfaceScalarField> dragCorrByADfs
425  (
426  invADVfs & dragCorrfs
427  );
428 
429  forAll(movingPhases, movingPhasei)
430  {
431  phaseModel& phase = movingPhases_[movingPhasei];
432 
433  phase.URef() =
434  HbyADs[movingPhasei]
436  (
437  alphaByADfs[movingPhasei]*mSfGradp
438  - FgByADfs[movingPhasei]
439  + dragCorrByADfs[movingPhasei]
440  )
441  - dragCorrByADs[movingPhasei];
442  }
443  }
444  else
445  {
446  forAll(movingPhases, movingPhasei)
447  {
448  phaseModel& phase = movingPhases_[movingPhasei];
449 
450  phase.URef() =
451  HbyADs[movingPhasei]
453  (
454  alphaByADfs[movingPhasei]*mSfGradp
455  - FgByADfs[movingPhasei]
456  );
457  }
458  }
459 
460  forAll(movingPhases, movingPhasei)
461  {
462  phaseModel& phase = movingPhases_[movingPhasei];
463 
464  phase.URef().correctBoundaryConditions();
465  phase.correctUf();
466  fvConstraints().constrain(phase.URef());
467  }
468  }
469  }
470 
471  // Update and limit the static pressure
474 
475  // Account for static pressure reference
477  {
479  (
480  "p",
481  p.dimensions(),
484  );
485  }
486 
487  // Limit p_rgh
489 
490  // Update densities from change in p_rgh
491  forAll(phases, phasei)
492  {
493  phaseModel& phase = phases_[phasei];
494  if (!phase.incompressible())
495  {
496  phase.rho() += phase.fluidThermo().psi()*(p_rgh - p_rgh_0);
497  }
498  }
499 
500  // Update mass transfer rates for change in p_rgh
501  forAll(phases, phasei)
502  {
503  if (dmdts.set(phasei) && d2mdtdps.set(phasei))
504  {
505  dmdts[phasei] += d2mdtdps[phasei]*(p_rgh - p_rgh_0);
506  }
507  }
508 
509  // Correct p_rgh for consistency with p and the updated densities
510  rho = fluid.rho();
512  }
513 
514  UEqns.clear();
515 }
516 
517 
518 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const dimensionSet & dimensions() const
Return dimensions.
static const DimensionedField< Type, GeoMesh > & null()
Return a null DimensionedField.
void relax(const scalar alpha)
Relax field (for steady-state solution).
GeometricBoundaryField< Type, PatchField, GeoMesh > Boundary
Type of the boundary field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
bool needReference() const
Does the field need a reference level for solution.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
IOobject(const word &name, const fileName &instance, const objectRegistry &registry, readOption r=NO_READ, writeOption w=NO_WRITE, bool registerObject=true)
Construct from name, instance, registry, io options.
Definition: IOobject.C:167
const word & name() const
Return name.
Definition: IOobject.H:310
static word groupName(Name name, const word &group)
void makeRelative(volVectorField &U) const
Make the given absolute velocity relative within the MRF region.
Definition: MRFZoneList.C:157
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
const word & name() const
Return const reference to name.
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
const surfaceVectorField & Sf() const
Return cell face area vectors.
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
bool finalNonOrthogonalIter() const
Flag to indicate the last non-orthogonal iteration.
virtual PtrList< surfaceScalarField > ddtCorrs() const =0
Return the flux corrections for the cell-based algorithm.
void correctBoundaryFlux()
Correct fixed-flux BCs to be consistent with the velocity BCs.
Definition: phaseSystem.C:690
virtual bool implicitPhasePressure(const phaseModel &phase) const
Returns true if the phase pressure is treated implicitly.
Definition: phaseSystem.C:500
virtual void invADVs(const PtrList< volScalarField > &As, PtrList< volVectorField > &HVms, PtrList< PtrList< volScalarField >> &invADVs, PtrList< PtrList< surfaceScalarField >> &invADVfs) const =0
Return the inverse of the central + drag + virtual mass.
virtual void dragCorrs(PtrList< volVectorField > &dragCorrs, PtrList< surfaceScalarField > &dragCorrf) const =0
Set the cell and faces drag correction fields.
virtual PtrList< volScalarField > dmdts() const
Return the mass transfer rates for each phase.
Definition: phaseSystem.C:474
virtual PtrList< surfaceScalarField > Fs() const =0
Return the force fluxes for the cell-based algorithm.
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
Return the surface tension force.
Definition: phaseSystem.C:513
tmp< volScalarField > rho() const
Return the mixture density.
Definition: phaseSystem.C:344
virtual PtrList< volScalarField > d2mdtdps() const
Return the mass transfer pressure implicit coefficients.
Definition: phaseSystem.C:480
bool incompressible() const
Return incompressibility.
Definition: phaseSystem.C:486
bool correct()
Piso loop within outer loop.
bool correctNonOrthogonal()
Non-orthogonal corrector loop.
scalar refValue() const
Return the pressure reference level.
label refCell() const
Return the cell in which the reference pressure is set.
pimpleNoLoopControl pimple
PIMPLE inner-loop controls.
Definition: solver.H:107
Foam::fvConstraints & fvConstraints() const
Return the fvConstraints that are created on demand.
Definition: solver.C:107
const Time & runTime
Time.
Definition: solver.H:104
const fvMesh & mesh
Region mesh.
Definition: solver.H:101
uniformDimensionedVectorField g
Gravitational acceleration.
Definition: buoyancy.H:83
volScalarField gh
(g & h) - ghRef
Definition: buoyancy.H:95
uniformDimensionedScalarField pRef
Reference pressure.
Definition: buoyancy.H:89
surfaceScalarField ghf
(g & hf) - ghRef
Definition: buoyancy.H:98
volScalarField & p_rgh
Reference to the buoyant pressure for buoyant cases.
const surfaceScalarField & phi
Reference to the mass-flux field.
solvers::buoyancy buoyancy
Buoyancy force.
phaseSystem::phaseModelList & phases_
PtrList< fvVectorMatrix > UEqns
Temporary phase momentum matrices.
const phaseSystem::phaseModelPartialList & movingPhases
Reference to the moving phases.
Switch predictMomentum
Momentum equation predictor switch.
Switch dragCorrection
Cell/face drag correction for cell momentum corrector.
Foam::pressureReference pressureReference
Pressure reference.
PtrList< volScalarField > rAs
Temporary storage for the reciprocal momentum equation diagonal.
const phaseSystem::phaseModelList & phases
Reference to the phases.
phaseSystem::phaseModelPartialList & movingPhases_
const volScalarField & p
Reference to the pressure field.
const phaseSystem & fluid
Reference to the multiphase fluid.
Find the reference cell nearest (in index) to the given cell but which is not on a cyclic,...
Calculate the first temporal derivative.
Calculate the divergence of the given field.
Calculate the face-flux of the given field.
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
Reconstruct volField from a face flux field.
Calculate the snGrad of the given volField.
Calculate the field for explicit evaluation of implicit and explicit sources.
Calculate the matrix for the first temporal derivative.
Calculate the matrix for the divergence of the given field and flux.
Calculate the matrix for the laplacian of the field.
Calculate the matrix for implicit and explicit sources.
U
Definition: pEqn.H:72
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho) *fvc::flux(HbyA))
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
tmp< SurfaceField< typename innerProduct< vector, Type >::type > > flux(const VolField< Type > &vf)
Return the face-flux field obtained from the given volVectorField.
static tmp< SurfaceField< Type > > interpolate(const VolField< Type > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
tmp< VolField< typename outerProduct< vector, Type >::type > > reconstruct(const SurfaceField< Type > &ssf)
tmp< VolField< Type > > div(const SurfaceField< Type > &ssf)
Definition: fvcDiv.C:47
void makeRelative(surfaceScalarField &phi, const volVectorField &U)
Make the given flux relative.
Definition: fvcMeshPhi.C:89
tmp< SurfaceField< Type > > snGrad(const VolField< Type > &vf, const word &name)
Definition: fvcSnGrad.C:45
tmp< fvMatrix< Type > > laplacian(const VolField< Type > &vf, const word &name)
Definition: fvmLaplacian.C:47
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
const dimensionSet dimVolumetricFlux
SurfaceField< scalar > surfaceScalarField
tmp< volScalarField > byDt(const volScalarField &vf)
Definition: phaseSystem.C:855
void constrainHbyA(volVectorField &HbyA, const volVectorField &U, const volScalarField &p)
Definition: constrainHbyA.C:34
const dimensionSet dimTime
const dimensionSet dimDensity
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:64
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
scalar getRefCellValue(const volScalarField &field, const label refCelli)
Return the current value of field in the reference cell.
Definition: findRefCell.C:136
void constrainPhiHbyA(surfaceScalarField &phiHbyA, const volVectorField &U, const volScalarField &p)
Definition: constrainHbyA.C:60
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.