All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
correctBuoyantPressure.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-2023 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 "isothermalFluid.H"
27 #include "constrainHbyA.H"
28 #include "constrainPressure.H"
29 #include "adjustPhi.H"
30 #include "fvcMeshPhi.H"
31 #include "fvcFlux.H"
32 #include "fvcDdt.H"
33 #include "fvcSnGrad.H"
34 #include "fvcReconstruct.H"
35 #include "fvcVolumeIntegrate.H"
36 #include "fvmDiv.H"
37 #include "fvmLaplacian.H"
38 
39 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
40 
41 void Foam::solvers::isothermalFluid::correctBuoyantPressure()
42 {
47 
48  // Local references to the buoyancy parameters
49  const volScalarField& gh = buoyancy->gh;
50  const surfaceScalarField& ghf = buoyancy->ghf;
51  const uniformDimensionedScalarField pRef = buoyancy->pRef;
52 
53  const volScalarField& psi = thermo.psi();
54  rho = thermo.rho();
55  rho.relax();
56 
57  fvVectorMatrix& UEqn = tUEqn.ref();
58 
59  // Thermodynamic density needs to be updated by psi*d(p) after the
60  // pressure solution
61  const volScalarField psip0(psi*p);
62 
64 
65  const volScalarField rAU("rAU", 1.0/UEqn.A());
66  const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
67 
68  tmp<volScalarField> rAtU
69  (
71  ? volScalarField::New("rAtU", 1.0/(1.0/rAU - UEqn.H1()))
72  : tmp<volScalarField>(nullptr)
73  );
74 
75  tmp<surfaceScalarField> rhorAtUf
76  (
78  ? surfaceScalarField::New("rhoRAtUf", fvc::interpolate(rho*rAtU()))
79  : tmp<surfaceScalarField>(nullptr)
80  );
81 
82  const volScalarField& rAAtU = pimple.consistent() ? rAtU() : rAU;
83  const surfaceScalarField& rhorAAtUf =
84  pimple.consistent() ? rhorAtUf() : rhorAUf;
85 
87 
88  if (pimple.nCorrPiso() <= 1)
89  {
90  tUEqn.clear();
91  }
92 
94  (
95  "phiHbyA",
96  rhof*fvc::flux(HbyA)
97  + rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf)
98  );
99 
100  MRF.makeRelative(rhof, phiHbyA);
101 
102  const bool adjustMass =
104 
105  const surfaceScalarField ghGradRhof(-ghf*fvc::snGrad(rho)*mesh.magSf());
106 
107  phiHbyA += rhorAUf*ghGradRhof;
108 
109  tmp<fvScalarMatrix> tp_rghEqn;
110 
111  if (pimple.transonic())
112  {
113  const surfaceScalarField phidByPsi
114  (
116  (
117  fvc::relative(phiHbyA, rho, U)/rhof,
118  p_rgh
119  )
120  );
121 
122  const surfaceScalarField phid("phid", fvc::interpolate(psi)*phidByPsi);
123 
124  // Subtract the compressible part
125  // The resulting flux will be zero for a perfect gas
126  phiHbyA -= fvc::interpolate(psi*p_rgh)*phidByPsi;
127 
128  if (pimple.consistent())
129  {
130  const surfaceScalarField gradpf(fvc::snGrad(p_rgh)*mesh.magSf());
131  phiHbyA += (rhorAAtUf - rhorAUf)*gradpf;
132  HbyA += (rAAtU - rAU)*fvc::reconstruct(gradpf - ghGradRhof);
133  }
134 
135  // Update the pressure BCs to ensure flux consistency
136  constrainPressure(p_rgh, rho, U, phiHbyA, rhorAAtUf, MRF);
137 
139 
140  fvScalarMatrix p_rghDDtEqn
141  (
143  + fvc::div(phiHbyA) + fvm::div(phid, p_rgh)
144  ==
145  fvModels().sourceProxy(rho, p_rgh)
146  );
147 
148  while (pimple.correctNonOrthogonal())
149  {
150  tp_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAAtUf, p);
151  fvScalarMatrix& p_rghEqn = tp_rghEqn.ref();
152 
153  // Relax the pressure equation to ensure diagonal-dominance
154  p_rghEqn.relax();
155 
156  p_rghEqn.setReference
157  (
160  );
161 
162  fvConstraints().constrain(p_rghEqn);
163 
164  p_rghEqn.solve();
165  }
166  }
167  else
168  {
169  if (pimple.consistent())
170  {
171  const surfaceScalarField gradpf(fvc::snGrad(p_rgh)*mesh.magSf());
172  phiHbyA += (rhorAAtUf - rhorAUf)*gradpf;
173  HbyA += (rAAtU - rAU)*fvc::reconstruct(gradpf - ghGradRhof);
174  }
175 
176  // Update the pressure BCs to ensure flux consistency
177  constrainPressure(p_rgh, rho, U, phiHbyA, rhorAAtUf, MRF);
178 
180 
181  fvScalarMatrix p_rghDDtEqn
182  (
184  + fvc::div(phiHbyA)
185  ==
186  fvModels().sourceProxy(rho, p_rgh)
187  );
188 
189  while (pimple.correctNonOrthogonal())
190  {
191  tp_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAAtUf, p_rgh);
192  fvScalarMatrix& p_rghEqn = tp_rghEqn.ref();
193 
194  p_rghEqn.setReference
195  (
198  );
199 
200  fvConstraints().constrain(p_rghEqn);
201 
202  p_rghEqn.solve();
203  }
204  }
205 
206  const fvScalarMatrix& p_rghEqn = tp_rghEqn();
207 
208  phi = phiHbyA + p_rghEqn.flux();
209 
210  // Calculate and relax the net pressure-buoyancy force
211  netForce.ref().relax
212  (
213  fvc::reconstruct((ghGradRhof + p_rghEqn.flux()/rhorAAtUf)),
215  );
216 
217  // Correct the momentum source with the pressure gradient flux
218  // calculated from the relaxed pressure
219  U = HbyA + rAAtU*netForce();
222 
223  K = 0.5*magSqr(U);
224 
225  if (!mesh.schemes().steady())
226  {
227  p = p_rgh + rho*gh + pRef;
228 
229  const bool constrained = fvConstraints().constrain(p);
230 
231  // Thermodynamic density update
232  thermo_.correctRho(psi*p - psip0);
233 
234  if (constrained)
235  {
236  rho = thermo.rho();
237  }
238 
239  correctDensity();
240  }
241 
242  continuityErrors();
243 
244  p = p_rgh + rho*gh + pRef;
245 
246  if (mesh.schemes().steady())
247  {
248  if (fvConstraints().constrain(p))
249  {
250  p_rgh = p - rho*gh - pRef;
252  }
253  }
254 
255  // For steady compressible closed-volume cases adjust the pressure level
256  // to obey overall mass continuity
257  if (adjustMass && !thermo.incompressible())
258  {
261  p_rgh = p - rho*gh - pRef;
263  }
264 
265  // Optionally relax pressure for the thermophysics
266  p.relax();
267 
268  if (mesh.schemes().steady() || pimple.simpleRho() || adjustMass)
269  {
270  rho = thermo.rho();
271  }
272 
273  if (mesh.schemes().steady() || pimple.simpleRho())
274  {
275  rho.relax();
276  }
277 
278  // Correct rhoUf if the mesh is moving
280 
281  if (thermo.dpdt())
282  {
283  dpdt = fvc::ddt(p);
284  }
285 }
286 
287 
288 // ************************************************************************* //
For cases which do no have a pressure boundary adjust the balance of fluxes to obey continuity....
void relax(const scalar alpha)
Relax field (for steady-state solution).
scalar relaxationFactor() const
Return the field relaxation factor read from fvSolution.
void correctBoundaryConditions()
Correct boundary field.
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &, const HashPtrTable< Source > &=HashPtrTable< Source >())
Return a temporary field constructed from name,.
void makeRelative(volVectorField &U) const
Make the given absolute velocity relative within the MRF region.
Definition: MRFZoneList.C:157
virtual Switch dpdt() const =0
Should the dpdt term be included in the enthalpy equation.
virtual bool incompressible() const =0
Return true if the equation of state is incompressible.
virtual tmp< volScalarField > rho() const =0
Density [kg/m^3].
bool transonic() const
Flag to indicate to solve using transonic algorithm.
bool consistent() const
Flag to indicate to relax pressure using the "consistent".
virtual void correctRho(const volScalarField &deltaRho)=0
Add the given density correction to the density field.
virtual const volScalarField & psi() const =0
Compressibility [s^2/m^2].
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
const fvSchemes & schemes() const
Return the fvSchemes.
Definition: fvMesh.C:1751
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
bool steady() const
Return true if the default ddtScheme is steadyState.
Definition: fvSchemes.H:147
bool simpleRho() const
Switch to indicate whether to update the density in simple mode.
bool correctNonOrthogonal()
Non-orthogonal corrector loop.
label nCorrPiso() const
Maximum number of piso correctors.
Definition: pisoControlI.H:28
scalar refValue() const
Return the pressure reference level.
label refCell() const
Return the cell in which the reference pressure is set.
Foam::fvModels & fvModels() const
Return the fvModels that are created on demand.
Definition: solver.C:96
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 fvMesh & mesh
Region mesh.
Definition: solver.H:101
volScalarField rho_
The continuity density field.
volScalarField & p_rgh
Reference to the buoyant pressure for buoyant cases.
const surfaceScalarField & phi
Mass-flux field.
volVectorField U_
Velocity field.
tmp< fvVectorMatrix > tUEqn
Cached momentum matrix.
tmp< volVectorField > netForce
Momentum equation net force source term.
autoPtr< surfaceVectorField > rhoUf
Pointer to the surface momentum field.
fluidThermo & thermo_
Reference to the fluid thermophysical properties.
volScalarField K
Kinetic energy field.
autoPtr< solvers::buoyancy > buoyancy
Pointer to the optional buoyancy force.
const volVectorField & U
Velocity field.
IOMRFZoneList MRF
MRF zone list.
volScalarField & p_
Reference to the pressure field.
surfaceScalarField phi_
Mass-flux field.
const volScalarField & rho
Reference to the continuity density field.
volScalarField::Internal dpdt
Rate of change of the pressure.
Foam::pressureReference pressureReference
Pressure reference.
const fluidThermo & thermo
Reference to the fluid thermophysical properties.
const volScalarField & p
Reference to the pressure field.
dimensionedScalar initialMass
Initial mass in the region.
fvVectorMatrix & UEqn
Definition: UEqn.H:11
Calculate the first temporal derivative.
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.
Volume integrate volField creating a volField.
Calculate the matrix for the divergence of the given field and flux.
Calculate the matrix for the laplacian of the field.
const volScalarField & psi
volScalarField rAU(1.0/UEqn.A())
volVectorField & HbyA
Definition: pEqn.H:13
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho) *fvc::flux(HbyA))
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< Type > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
void correctRhoUf(autoPtr< surfaceVectorField > &rhoUf, const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const MRFType &MRF)
tmp< surfaceScalarField > relative(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given absolute flux in relative form.
Definition: fvcMeshPhi.C:167
tmp< VolField< typename outerProduct< vector, Type >::type > > reconstruct(const SurfaceField< Type > &ssf)
dimensioned< Type > domainIntegrate(const VolField< Type > &vf)
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< SurfaceField< typename Foam::flux< Type >::type > > ddtCorr(const VolField< Type > &U, const SurfaceField< Type > &Uf)
Definition: fvcDdt.C:196
tmp< fvMatrix< Type > > laplacian(const VolField< Type > &vf, const word &name)
Definition: fvmLaplacian.C:47
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const VolField< Type > &vf, const word &name)
Definition: fvmDiv.C:48
tmp< fvMatrix< Type > > ddt(const VolField< Type > &vf)
Definition: fvmDdt.C:46
void constrainPressure(volScalarField &p, const RhoType &rho, const volVectorField &U, const surfaceScalarField &phiHbyA, const RAUType &rhorAU, const MRFType &MRF)
bool adjustPhi(surfaceScalarField &phi, const volVectorField &U, volScalarField &p)
Adjust the balance of fluxes to obey continuity.
Definition: adjustPhi.C:35
VolField< vector > volVectorField
Definition: volFieldsFwd.H:65
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
UniformDimensionedField< scalar > uniformDimensionedScalarField
SurfaceField< scalar > surfaceScalarField
void constrainHbyA(volVectorField &HbyA, const volVectorField &U, const volScalarField &p)
Definition: constrainHbyA.C:34
tmp< surfaceScalarField > constrainPhid(const tmp< surfaceScalarField > &tphiHbyA, const volScalarField &p)
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:64
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
fvConstraints constrain(rhoEqn)