facePressureCorrector.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::facePressureCorrector()
45 {
47 
48  // Face volume fractions
49  PtrList<surfaceScalarField> alphafs(phases.size());
50  forAll(phases, phasei)
51  {
52  const phaseModel& phase = phases[phasei];
53  const volScalarField& alpha = phase;
54 
55  alphafs.set(phasei, fvc::interpolate(alpha).ptr());
56  alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
57  }
58 
59  // Diagonal coefficients
60  rAs.clear();
62  {
63  rAs.setSize(phases.size());
64  }
65 
66  PtrList<surfaceScalarField> HVmfs(movingPhases.size());
67  PtrList<PtrList<surfaceScalarField>> invADVfs;
68  {
69  PtrList<surfaceScalarField> Afs(movingPhases.size());
70 
71  forAll(movingPhases, movingPhasei)
72  {
73  const phaseModel& phase = movingPhases[movingPhasei];
74  const volScalarField& alpha = phase;
75 
76  const volScalarField A
77  (
78  byDt
79  (
80  max(alpha.oldTime(), phase.residualAlpha())
81  *phase.rho().oldTime()
82  )
83  + UEqns[phase.index()].A()
84  );
85 
87  {
88  rAs.set
89  (
90  phase.index(),
91  new volScalarField
92  (
93  IOobject::groupName("rA", phase.name()),
94  1/A
95  )
96  );
97  }
98 
99  Afs.set
100  (
101  movingPhasei,
103  (
104  IOobject::groupName("rAf", phase.name()),
106  )
107  );
108  }
109 
110  invADVfs = fluid.invADVfs(Afs, HVmfs);
111  }
112 
113  volScalarField rho("rho", fluid.rho());
114 
115  // Phase diagonal coefficients
116  PtrList<surfaceScalarField> alphaByADfs;
117  PtrList<surfaceScalarField> FgByADfs;
118  {
119  // Explicit force fluxes
120  PtrList<surfaceScalarField> Ffs(fluid.Ffs());
121 
122  const surfaceScalarField ghSnGradRho
123  (
124  "ghSnGradRho",
126  );
127 
128  PtrList<surfaceScalarField> lalphafs(movingPhases.size());
129  PtrList<surfaceScalarField> Fgfs(movingPhases.size());
130 
131  forAll(movingPhases, movingPhasei)
132  {
133  const phaseModel& phase = movingPhases[movingPhasei];
134 
135  lalphafs.set
136  (
137  movingPhasei,
138  max(alphafs[phase.index()], phase.residualAlpha())
139  );
140 
141  Fgfs.set
142  (
143  movingPhasei,
144  Ffs[phase.index()]
145  + lalphafs[movingPhasei]
146  *(
147  ghSnGradRho
148  - (fvc::interpolate(phase.rho() - rho))
149  *(buoyancy.g & mesh.Sf())
150  - fluid.surfaceTension(phase)*mesh.magSf()
151  )
152  );
153  }
154 
155  alphaByADfs = invADVfs & lalphafs;
156  FgByADfs = invADVfs & Fgfs;
157  }
158 
159 
160  // Mass transfer rates
161  PtrList<volScalarField> dmdts(fluid.dmdts());
162  PtrList<volScalarField> d2mdtdps(fluid.d2mdtdps());
163 
164  // --- Pressure corrector loop
165  while (pimple.correct())
166  {
167  // Correct fixed-flux BCs to be consistent with the velocity BCs
169 
170  // Predicted fluxes for each phase
171  PtrList<surfaceScalarField> phiHbyADs;
172  {
173  PtrList<surfaceScalarField> phiHs(movingPhases.size());
174 
175  forAll(movingPhases, movingPhasei)
176  {
177  const phaseModel& phase = movingPhases[movingPhasei];
178  const volScalarField& alpha = phase;
179 
180  phiHs.set
181  (
182  movingPhasei,
183  (
185  (
186  max(alpha.oldTime(), phase.residualAlpha())
187  *phase.rho().oldTime()
188  )
189  *byDt
190  (
191  phase.Uf().valid()
192  ? (mesh.Sf() & phase.Uf()().oldTime())
193  : MRF.absolute(phase.phi()().oldTime())
194  )
195  + fvc::flux(UEqns[phase.index()].H())
196  )
197  );
198 
199  if (HVmfs.set(movingPhasei))
200  {
201  phiHs[movingPhasei] += HVmfs[movingPhasei];
202  }
203  }
204 
205  phiHbyADs = invADVfs & phiHs;
206  }
207 
208  forAll(movingPhases, movingPhasei)
209  {
210  const phaseModel& phase = movingPhases[movingPhasei];
211 
212  constrainPhiHbyA(phiHbyADs[movingPhasei], phase.U(), p_rgh);
213 
214  phiHbyADs[movingPhasei] -= FgByADfs[movingPhasei];
215  }
216 
217  // Total predicted flux
219  (
220  IOobject
221  (
222  "phiHbyA",
223  runTime.name(),
224  mesh,
227  ),
228  mesh,
230  );
231 
232  forAll(movingPhases, movingPhasei)
233  {
234  const phaseModel& phase = movingPhases[movingPhasei];
235 
236  phiHbyA += alphafs[phase.index()]*phiHbyADs[movingPhasei];
237  }
238 
241 
242  // Construct pressure "diffusivity"
244  (
245  IOobject
246  (
247  "rAf",
248  runTime.name(),
249  mesh
250  ),
251  mesh,
252  dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), 0)
253  );
254 
255  forAll(movingPhases, movingPhasei)
256  {
257  const phaseModel& phase = movingPhases[movingPhasei];
258 
259  rAf += alphafs[phase.index()]*alphaByADfs[movingPhasei];
260  }
261 
262  // Update the fixedFluxPressure BCs to ensure flux consistency
263  {
265  (
268  );
269  phib = 0;
270 
271  forAll(movingPhases, movingPhasei)
272  {
273  phaseModel& phase = movingPhases_[movingPhasei];
274 
275  phib +=
276  alphafs[phase.index()].boundaryField()
277  *phase.phi()().boundaryField();
278  }
279 
280  setSnGrad<fixedFluxPressureFvPatchScalarField>
281  (
283  (
284  phiHbyA.boundaryField() - phib
285  )/(mesh.magSf().boundaryField()*rAf.boundaryField())
286  );
287  }
288 
289  // Compressible pressure equations
290  PtrList<fvScalarMatrix> pEqnComps(compressibilityEqns(dmdts, d2mdtdps));
291 
292  // Cache p prior to solve for density update
293  volScalarField p_rgh_0(p_rgh);
294 
295  // Iterate over the pressure equation to correct for non-orthogonality
296  while (pimple.correctNonOrthogonal())
297  {
298  // Construct the transport part of the pressure equation
299  fvScalarMatrix pEqnIncomp
300  (
302  - fvm::laplacian(rAf, p_rgh)
303  );
304 
305  // Solve
306  {
307  fvScalarMatrix pEqn(pEqnIncomp);
308 
309  forAll(phases, phasei)
310  {
311  pEqn += pEqnComps[phasei];
312  }
313 
314  if (fluid.incompressible())
315  {
316  pEqn.setReference
317  (
320  );
321  }
322 
323  fvConstraints().constrain(pEqn);
324 
325  pEqn.solve();
326  }
327 
328  // Correct fluxes and velocities on last non-orthogonal iteration
330  {
331  phi_ = phiHbyA + pEqnIncomp.flux();
332 
333  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAf);
334 
335  forAll(movingPhases, movingPhasei)
336  {
337  phaseModel& phase = movingPhases_[movingPhasei];
338  const label phasei = phase.index();
339 
340  phase.phiRef() =
341  phiHbyADs[movingPhasei]
342  + alphaByADfs[movingPhasei]*mSfGradp;
343 
344  // Set the phase dilatation rate
345  phase.divU(-pEqnComps[phasei] & p_rgh);
346  }
347 
348  forAll(movingPhases, movingPhasei)
349  {
350  phaseModel& phase = movingPhases_[movingPhasei];
351 
352  MRF.makeRelative(phase.phiRef());
353  fvc::makeRelative(phase.phiRef(), phase.U());
354 
355  phase.URef() = fvc::reconstruct
356  (
357  fvc::absolute(MRF.absolute(phase.phi()), phase.U())
358  );
359 
360  phase.URef().correctBoundaryConditions();
361  phase.correctUf();
362  fvConstraints().constrain(phase.URef());
363  }
364  }
365  }
366 
367  // Update and limit the static pressure
370 
371  // Account for static pressure reference
373  {
375  (
376  "p",
377  p.dimensions(),
380  );
381  }
382 
383  // Limit p_rgh
385 
386  // Update densities from change in p_rgh
387  forAll(phases, phasei)
388  {
389  phaseModel& phase = phases_[phasei];
390  if (!phase.incompressible())
391  {
392  phase.rho() += phase.fluidThermo().psi()*(p_rgh - p_rgh_0);
393  }
394  }
395 
396  // Update mass transfer rates for change in p_rgh
397  forAll(phases, phasei)
398  {
399  if (dmdts.set(phasei) && d2mdtdps.set(phasei))
400  {
401  dmdts[phasei] += d2mdtdps[phasei]*(p_rgh - p_rgh_0);
402  }
403  }
404 
405  // Correct p_rgh for consistency with p and the updated densities
406  rho = fluid.rho();
409  }
410 
411  UEqns.clear();
412 }
413 
414 
415 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
#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.
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.
void correctBoundaryConditions()
Correct boundary field.
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)
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &phi) const
Return the given relative flux absolute within the MRF region.
Definition: MRFZoneList.C:291
void makeRelative(volVectorField &U) const
Make the given absolute velocity relative within the MRF region.
Definition: MRFZoneList.C:157
const FieldType & oldTime() const
Return the old time field.
Definition: OldTimeField.C:315
bool set(const label) const
Is element set.
Definition: UPtrListI.H:87
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 > Ffs() const =0
Return the force fluxes for the face-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 PtrList< volScalarField > dmdts() const
Return the mass transfer rates for each phase.
Definition: phaseSystem.C:474
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< PtrList< surfaceScalarField > > invADVfs(const PtrList< surfaceScalarField > &Afs, PtrList< surfaceScalarField > &HVmfs) const =0
Return the inverse of the face central + drag + virtual mass.
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.
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
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:202
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
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
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
const dimensionSet dimVolumetricFlux
SurfaceField< scalar > surfaceScalarField
tmp< volScalarField > byDt(const volScalarField &vf)
Definition: phaseSystem.C:855
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.