From 08baa6eda69a667a8560f612611386d9af83beed Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 11 Sep 2013 00:10:00 +0100
Subject: [PATCH] fixedFluxPressure BC: the snGrad is now pushed into the BC
from pEqn.H rather than being evaluated in the BC
---
applications/solvers/DNS/dnsFoam/dnsFoam.C | 6 +-
applications/solvers/combustion/XiFoam/pEqn.H | 10 +-
.../solvers/combustion/engineFoam/pEqn.H | 6 +-
.../solvers/combustion/fireFoam/fireFoam.C | 1 +
.../solvers/combustion/fireFoam/pEqn.H | 26 ++-
.../solvers/combustion/reactingFoam/pEqn.H | 6 +-
.../rhoReactingBuoyantFoam/pEqn.H | 21 +-
.../rhoReactingBuoyantFoam.C | 1 +
.../reactingFoam/rhoReactingFoam/pEqn.H | 10 +-
.../solvers/compressible/rhoPimpleFoam/pEqn.H | 12 +-
.../rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H | 10 +-
.../rhoPimpleFoam/rhoPimplecFoam/pEqn.H | 8 +-
.../rhoSimpleFoam/rhoSimplecFoam/pEqn.H | 8 +-
.../solvers/compressible/sonicFoam/pEqn.H | 6 +-
.../sonicFoam/sonicDyMFoam/pEqn.H | 6 +-
.../sonicLiquidFoam/sonicLiquidFoam.C | 12 +-
.../electromagnetics/mhdFoam/mhdFoam.C | 12 +-
.../buoyantBoussinesqPimpleFoam.C | 1 +
.../buoyantBoussinesqPimpleFoam/pEqn.H | 20 +-
.../buoyantBoussinesqSimpleFoam.C | 1 +
.../buoyantBoussinesqSimpleFoam/pEqn.H | 12 +-
.../buoyantPimpleFoam/buoyantPimpleFoam.C | 1 +
.../heatTransfer/buoyantPimpleFoam/pEqn.H | 21 +-
.../buoyantSimpleFoam/buoyantSimpleFoam.C | 1 +
.../heatTransfer/buoyantSimpleFoam/pEqn.H | 15 +-
.../chtMultiRegionFoam/chtMultiRegionFoam.C | 1 +
.../chtMultiRegionSimpleFoam.C | 1 +
.../chtMultiRegionSimpleFoam/fluid/pEqn.H | 15 +-
.../chtMultiRegionFoam/fluid/pEqn.H | 23 +-
.../solvers/incompressible/pimpleFoam/pEqn.H | 6 +-
.../pimpleFoam/pimpleDyMFoam/pEqn.H | 6 +-
.../lagrangian/coalChemistryFoam/pEqn.H | 10 +-
.../lagrangian/reactingParcelFilmFoam/pEqn.H | 21 +-
.../reactingParcelFilmFoam.C | 1 +
.../lagrangian/reactingParcelFoam/pEqn.H | 6 +-
.../solvers/lagrangian/sprayFoam/pEqn.H | 10 +-
.../sprayFoam/sprayEngineFoam/pEqn.H | 10 +-
.../cavitatingFoam/cavitatingDyMFoam/pEqn.H | 8 +-
.../solvers/multiphase/cavitatingFoam/pEqn.H | 8 +-
.../compressibleInterDyMFoam.C | 3 +-
.../compressibleInterFoam.C | 1 +
.../multiphase/compressibleInterFoam/pEqn.H | 19 +-
.../interFoam/LTSInterFoam/LTSInterFoam.C | 1 +
.../interFoam/MRFInterFoam/MRFInterFoam.C | 1 +
.../multiphase/interFoam/MRFInterFoam/pEqn.H | 19 +-
.../interFoam/interDyMFoam/interDyMFoam.C | 3 +-
.../multiphase/interFoam/interDyMFoam/pEqn.H | 20 +-
.../solvers/multiphase/interFoam/interFoam.C | 1 +
.../interMixingFoam/interMixingFoam.C | 1 +
.../solvers/multiphase/interFoam/pEqn.H | 19 +-
.../porousInterFoam/porousInterFoam.C | 1 +
.../interPhaseChangeDyMFoam.C | 5 +-
.../interPhaseChangeDyMFoam/pEqn.H | 20 +-
.../interPhaseChangeFoam.C | 1 +
.../multiphase/interPhaseChangeFoam/pEqn.H | 19 +-
.../multiphaseEulerFoam/multiphaseEulerFoam.C | 6 +-
.../multiphase/multiphaseEulerFoam/pEqn.H | 62 +++---
.../MRFMultiphaseInterFoam.C | 1 +
.../MRFMultiphaseInterFoam/pEqn.H | 19 +-
.../multiphaseInterFoam/multiphaseInterFoam.C | 1 +
.../multiphase/multiphaseInterFoam/pEqn.H | 19 +-
.../potentialFreeSurfaceFoam/pEqn.H | 11 +-
.../potentialFreeSurfaceFoam.C | 1 +
.../solvers/multiphase/settlingFoam/pEqn.H | 22 +-
.../multiphase/settlingFoam/settlingFoam.C | 1 +
.../multiphase/twoLiquidMixingFoam/pEqn.H | 19 +-
.../twoLiquidMixingFoam/twoLiquidMixingFoam.C | 1 +
.../multiphase/twoPhaseEulerFoam/pEqn.H | 36 ++--
.../twoPhaseEulerFoam/twoPhaseEulerFoam.C | 1 +
src/finiteVolume/Make/files | 1 -
.../fixedFluxPressureFvPatchScalarField.C | 105 ++-------
.../fixedFluxPressureFvPatchScalarField.H | 90 ++++----
...phaseFixedFluxPressureFvPatchScalarField.C | 197 -----------------
...phaseFixedFluxPressureFvPatchScalarField.H | 199 ------------------
.../laminar/depthCharge2D/0/p_rgh.org | 1 -
.../ras/damBreakWithObstacle/constant/p_rgh | 11 -
.../ras/sloshingTank2D3DoF/0/p_rgh | 2 -
77 files changed, 526 insertions(+), 773 deletions(-)
delete mode 100644 src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
delete mode 100644 src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H
diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index c29d4510915..865af56184b 100644
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -86,7 +86,7 @@ int main(int argc, char *argv[])
for (int corr=1; corr<=1; corr++)
{
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -94,12 +94,12 @@ int main(int argc, char *argv[])
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(U, phi)
+ + rAUf*fvc::ddtCorr(U, phi)
);
fvScalarMatrix pEqn
(
- fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
);
pEqn.solve();
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index 83a44ab7c22..995a59a737e 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -1,7 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
fvc::interpolate(psi)
*(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
);
@@ -26,7 +26,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
@@ -48,7 +48,7 @@ else
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
@@ -60,7 +60,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index a92e7285ab2..384ecbee598 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -1,7 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -15,7 +15,7 @@ if (pimple.transonic())
*(
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+ //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
- fvc::meshPhi(rho, U)
)
@@ -51,7 +51,7 @@ else
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+ //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
- fvc::interpolate(rho)*fvc::meshPhi(rho, U)
);
diff --git a/applications/solvers/combustion/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C
index 14e46bef661..b8c48549a68 100644
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@@ -41,6 +41,7 @@ Description
#include "psiCombustionModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index c98ac61f3b4..6652eab5649 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -1,27 +1,35 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
-phi.boundaryField() =
- fvc::interpolate(rho.boundaryField()*U.boundaryField())
- & mesh.Sf().boundaryField();
-surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
+ phig
);
-fvOptions.makeRelative(phiHbyA);
+fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the fixedFluxPressure BCs to ensure flux consistency
+setSnGrad<fixedFluxPressureFvPatchScalarField>
+(
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+ *rho.boundaryField()
+ )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+);
while (pimple.correctNonOrthogonal())
{
@@ -30,7 +38,7 @@ while (pimple.correctNonOrthogonal())
fvc::ddt(psi, rho)*gh
+ fvc::div(phiHbyA)
+ fvm::ddt(psi, p_rgh)
- - fvm::laplacian(Dp, p_rgh)
+ - fvm::laplacian(rhorAUf, p_rgh)
==
parcels.Srho()
+ surfaceFilm.Srho()
@@ -44,7 +52,7 @@ while (pimple.correctNonOrthogonal())
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/Dp);
+ U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index 0dc5d422db9..f7a3004ac21 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -1,7 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
fvc::interpolate(psi)
*(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
);
@@ -48,7 +48,7 @@ else
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
index 8602725bdb4..8894c7a9d76 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@@ -6,25 +6,36 @@
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
- surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+ surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
+ phig
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+ *rho.boundaryField()
+ )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+ );
+
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
@@ -38,7 +49,7 @@
fvScalarMatrix p_rghEqn
(
p_rghDDtEqn
- - fvm::laplacian(Dp, p_rgh)
+ - fvm::laplacian(rhorAUf, p_rgh)
);
fvOptions.constrain(p_rghEqn);
@@ -55,7 +66,7 @@
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
- U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
index 19e9f96f512..fd4b6f58eb5 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
@@ -36,6 +36,7 @@ Description
#include "multivariateScheme.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
index 3c92ce81197..f9b1bd65e41 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
@@ -6,7 +6,7 @@
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -18,7 +18,7 @@
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
);
@@ -39,7 +39,7 @@
fvScalarMatrix pEqn
(
pDDtEqn
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
@@ -61,7 +61,7 @@
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
@@ -80,7 +80,7 @@
fvScalarMatrix pEqn
(
pDDtEqn
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
);
fvOptions.constrain(pEqn);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 373bafcd3b8..a4e5fa20da7 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -4,7 +4,7 @@ rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn().A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -22,7 +22,7 @@ if (pimple.transonic())
fvc::interpolate(psi)
*(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
);
@@ -34,7 +34,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
@@ -56,14 +56,12 @@ else
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
- volScalarField Dp("Dp", rho*rAU);
-
while (pimple.correctNonOrthogonal())
{
// Pressure corrector
@@ -71,7 +69,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index d2681663e63..1d70a0f354b 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -4,7 +4,7 @@ rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn().A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -22,7 +22,7 @@ if (pimple.transonic())
fvc::interpolate(psi)
*(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phiAbs)
+ + rhorAUf*fvc::ddtCorr(rho, U, phiAbs)
)/fvc::interpolate(rho)
);
@@ -34,7 +34,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
@@ -55,7 +55,7 @@ else
(
"phiHbyA",
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phiAbs)
+ + rhorAUf*fvc::ddtCorr(rho, U, phiAbs)
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
@@ -67,7 +67,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
index f8837f0aa4d..0503db4207f 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
@@ -35,7 +35,7 @@ if (pimple.transonic())
HbyA -= (rAU - rAtU)*fvc::grad(p);
- volScalarField Dp("Dp", rho*rAtU);
+ volScalarField rhorAtU("rhorAtU", rho*rAtU);
while (pimple.correctNonOrthogonal())
{
@@ -44,7 +44,7 @@ if (pimple.transonic())
fvm::ddt(psi, p)
+ fvm::div(phid, p)
+ fvc::div(phic)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAtU, p)
==
fvOptions(psi, p, rho.name())
);
@@ -75,7 +75,7 @@ else
phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU)*fvc::grad(p);
- volScalarField Dp("Dp", rho*rAtU);
+ volScalarField rhorAtU("rhorAtU", rho*rAtU);
while (pimple.correctNonOrthogonal())
{
@@ -83,7 +83,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAtU, p)
==
fvOptions(psi, p, rho.name())
);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
index d61a80ffa81..c656b3c3d78 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
@@ -26,7 +26,7 @@ if (simple.transonic())
HbyA -= (rAU - rAtU)*fvc::grad(p);
- volScalarField Dp("Dp", rho*rAtU);
+ volScalarField rhorAtU("rhorAtU", rho*rAtU);
while (simple.correctNonOrthogonal())
{
@@ -34,7 +34,7 @@ if (simple.transonic())
(
fvm::div(phid, p)
+ fvc::div(phic)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAtU, p)
==
fvOptions(psi, p, rho.name())
);
@@ -67,14 +67,14 @@ else
phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU)*fvc::grad(p);
- volScalarField Dp("Dp", rho*rAtU);
+ volScalarField rhorAtU("rhorAtU", rho*rAtU);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAtU, p)
==
fvOptions(psi, p, rho.name())
);
diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H
index 2832bf421a9..edde4b198e5 100644
--- a/applications/solvers/compressible/sonicFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/pEqn.H
@@ -1,7 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -12,7 +12,7 @@ surfaceScalarField phid
fvc::interpolate(psi)*
(
(mesh.Sf() & fvc::interpolate(rho*HbyA))
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
);
@@ -23,7 +23,7 @@ while (pimple.correctNonOrthogonal())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
);
pEqn.solve();
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index 91becd8aa16..f8f99dafa54 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
@@ -1,7 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -13,7 +13,7 @@ surfaceScalarField phid
*(
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+ //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
- fvc::meshPhi(rho, U)
)
@@ -25,7 +25,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
);
pEqn.solve();
diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
index 2a0b08074a3..df3f3cdc61f 100644
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
@@ -75,8 +75,12 @@ int main(int argc, char *argv[])
// --- Pressure corrector loop
while (pimple.correct())
{
- volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+ volScalarField rAU("rAU", 1.0/UEqn.A());
+ surfaceScalarField rhorAUf
+ (
+ "rhorAUf",
+ fvc::interpolate(rho*rAU)
+ );
U = rAU*UEqn.H();
@@ -86,7 +90,7 @@ int main(int argc, char *argv[])
psi
*(
(fvc::interpolate(rho*U) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
);
@@ -97,7 +101,7 @@ int main(int argc, char *argv[])
fvm::ddt(psi, p)
+ fvc::div(phi)
+ fvm::div(phid, p)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
);
pEqn.solve();
diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
index 945aab6ba32..e9dedef5a93 100644
--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@@ -93,7 +93,7 @@ int main(int argc, char *argv[])
for (int corr=0; corr<nCorr; corr++)
{
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -102,14 +102,14 @@ int main(int argc, char *argv[])
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(U, phi)
+ + rAUf*fvc::ddtCorr(U, phi)
);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
- fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
@@ -143,14 +143,14 @@ int main(int argc, char *argv[])
BEqn.solve();
volScalarField rAB(1.0/BEqn.A());
- surfaceScalarField DpB("DpB", fvc::interpolate(rAB));
+ surfaceScalarField rABf("rABf", fvc::interpolate(rAB));
phiB = (fvc::interpolate(B) & mesh.Sf())
- + DpB*fvc::ddtCorr(B, phiB);
+ + rABf*fvc::ddtCorr(B, phiB);
fvScalarMatrix pBEqn
(
- fvm::laplacian(DpB, pB) == fvc::div(phiB)
+ fvm::laplacian(rABf, pB) == fvc::div(phiB)
);
pBEqn.solve();
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index 6168082e1de..43d8a4da534 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -51,6 +51,7 @@ Description
#include "radiationModel.H"
#include "fvIOoptionList.H"
#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
index 94299b84af1..d77eb12b0ee 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
@@ -1,25 +1,35 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
- surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rhok)*mesh.magSf());
+ surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(U, phi)
+ + rAUf*fvc::ddtCorr(U, phi)
+ phig
);
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -36,7 +46,7 @@
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
index b0d2f98aebf..483a723f7ce 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
@@ -50,6 +50,7 @@ Description
#include "RASModel.H"
#include "fvIOoptionList.H"
#include "simpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 53fec2e1b64..bae1bcef886 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn().A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -18,6 +18,16 @@
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
while (simple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index b1b6ceede26..9d62893f18a 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -39,6 +39,7 @@ Description
#include "radiationModel.H"
#include "fvIOoptionList.H"
#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index ad98545b407..5f6d70e054b 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -6,25 +6,36 @@
thermo.rho() -= psi*p_rgh;
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
- surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+ surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
(fvc::interpolate(rho*U) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rAUf*fvc::ddtCorr(rho, U, phi)
)
+ phig
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+ *rho.boundaryField()
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
@@ -38,7 +49,7 @@
fvScalarMatrix p_rghEqn
(
p_rghDDtEqn
- - fvm::laplacian(Dp, p_rgh)
+ - fvm::laplacian(rAUf, p_rgh)
);
fvOptions.constrain(p_rghEqn);
@@ -55,7 +66,7 @@
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
- U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index 67df5b61e05..2e52da8f531 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -36,6 +36,7 @@ Description
#include "radiationModel.H"
#include "simpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 41c8900dc08..60b3c1a0cfd 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -2,8 +2,8 @@
rho = thermo.rho();
rho.relax();
- volScalarField rAU(1.0/UEqn().A());
- surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+ volScalarField rAU("rAU", 1.0/UEqn().A());
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -23,6 +23,17 @@
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+ *rho.boundaryField()
+ )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+ );
+
while (simple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index b2cd9ac54e7..6aae0cf5b68 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -45,6 +45,7 @@ Description
#include "radiationModel.H"
#include "fvIOoptionList.H"
#include "coordinateSystem.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
index a3615c4309a..75720dd51c2 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
@@ -38,6 +38,7 @@ Description
#include "radiationModel.H"
#include "fvIOoptionList.H"
#include "coordinateSystem.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index f42ed7e39ce..fc854c40a50 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -4,8 +4,8 @@
rho = min(rho, rhoMax[i]);
rho.relax();
- volScalarField rAU(1.0/UEqn().A());
- surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+ volScalarField rAU("rAU", 1.0/UEqn().A());
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -25,6 +25,17 @@
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+ *rho.boundaryField()
+ )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+ );
+
dimensionedScalar compressibility = fvc::domainIntegrate(psi);
bool compressible = (compressibility.value() > SMALL);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index a712250f356..38ba49bd747 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -5,26 +5,37 @@
rho = thermo.rho();
- volScalarField rAU(1.0/UEqn().A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+ volScalarField rAU("rAU", 1.0/UEqn().A());
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
- surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+ surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
+ phig
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+ *rho.boundaryField()
+ )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+ );
+
{
fvScalarMatrix p_rghDDtEqn
(
@@ -41,7 +52,7 @@
fvScalarMatrix p_rghEqn
(
p_rghDDtEqn
- - fvm::laplacian(Dp, p_rgh)
+ - fvm::laplacian(rhorAUf, p_rgh)
);
p_rghEqn.solve
@@ -63,7 +74,7 @@
{
phi = phiHbyA + p_rghEqn.flux();
U = HbyA
- + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+ + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index 6824d8ee0c7..4799df09f72 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@@ -1,4 +1,4 @@
-surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -12,7 +12,7 @@ surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(U, phi)
+ + rAUf*fvc::ddtCorr(U, phi)
);
fvOptions.makeRelative(phiHbyA);
@@ -25,7 +25,7 @@ while (pimple.correctNonOrthogonal())
// Pressure corrector
fvScalarMatrix pEqn
(
- fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index 827daa4ce8c..7ecfe1f6623 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@@ -1,4 +1,4 @@
-surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -12,7 +12,7 @@ surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(U, Uf)
+ + rAUf*fvc::ddtCorr(U, Uf)
);
if (p.needReference())
@@ -26,7 +26,7 @@ while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
- fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index 200ea289d85..fce4b6e64da 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -1,7 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
fvc::interpolate(psi)
*(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
);
@@ -26,7 +26,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
coalParcels.Srho()
+ fvOptions(psi, p, rho.name())
@@ -49,7 +49,7 @@ else
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
@@ -61,7 +61,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
coalParcels.Srho()
+ fvOptions(psi, p, rho.name())
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 1e504462b6a..6652eab5649 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
@@ -1,25 +1,36 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
-surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
+ phig
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+// Update the fixedFluxPressure BCs to ensure flux consistency
+setSnGrad<fixedFluxPressureFvPatchScalarField>
+(
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+ *rho.boundaryField()
+ )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+);
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
@@ -27,7 +38,7 @@ while (pimple.correctNonOrthogonal())
fvc::ddt(psi, rho)*gh
+ fvc::div(phiHbyA)
+ fvm::ddt(psi, p_rgh)
- - fvm::laplacian(Dp, p_rgh)
+ - fvm::laplacian(rhorAUf, p_rgh)
==
parcels.Srho()
+ surfaceFilm.Srho()
@@ -41,7 +52,7 @@ while (pimple.correctNonOrthogonal())
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/Dp);
+ U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index 4067c125ebe..27f272100af 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -39,6 +39,7 @@ Description
#include "SLGThermo.H"
#include "fvIOoptionList.H"
#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 420bc1ef0cb..56b22dba9ca 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -6,7 +6,7 @@
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -16,7 +16,7 @@
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
@@ -36,7 +36,7 @@
fvScalarMatrix pEqn
(
pDDtEqn
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
);
fvOptions.constrain(pEqn);
diff --git a/applications/solvers/lagrangian/sprayFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/pEqn.H
index f62a17286c2..26511b001df 100644
--- a/applications/solvers/lagrangian/sprayFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/pEqn.H
@@ -1,7 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
fvc::interpolate(psi)
*(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
);
@@ -26,7 +26,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
@@ -49,7 +49,7 @@ else
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
@@ -61,7 +61,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
index b7de7e2f79e..62c3e0376e3 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
@@ -1,7 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -15,7 +15,7 @@ if (pimple.transonic())
*(
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+ //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
)/fvc::interpolate(rho)
- fvc::meshPhi(rho, U)
)
@@ -29,7 +29,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
@@ -52,7 +52,7 @@ else
"phiHbyA",
(
(fvc::interpolate(HbyA) & mesh.Sf())
- //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+ //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
- fvc::interpolate(rho)*fvc::meshPhi(rho, U)
);
@@ -65,7 +65,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index f869961a3c3..e7f078d1a65 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -12,16 +12,16 @@
surfaceScalarField rhof("rhof", fvc::interpolate(rho));
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
phiv = (fvc::interpolate(HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(U, phivAbs);
+ + rhorAUf*fvc::ddtCorr(U, phivAbs);
fvc::makeRelative(phiv, U);
- surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
+ surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
phiv -= phiGradp/rhof;
@@ -35,7 +35,7 @@
+ psi*correction(fvm::ddt(p))
+ fvc::div(phiv, rho)
+ fvc::div(phiGradp)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index e9d059aadbc..210a30fd143 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -12,15 +12,15 @@
surfaceScalarField rhof("rhof", fvc::interpolate(rho));
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
phiv = (fvc::interpolate(HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(U, phiv);
+ + rhorAUf*fvc::ddtCorr(U, phiv);
- surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
+ surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
phiv -= phiGradp/rhof;
@@ -32,7 +32,7 @@
- (rhol0 + (psil - psiv)*pSat)*fvc::ddt(alphav) - pSat*fvc::ddt(psi)
+ fvc::div(phiv, rho)
+ fvc::div(phiGradp)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rhorAUf, p)
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 7590e6ff45c..a9c3365bdf7 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -46,6 +46,7 @@ Description
#include "twoPhaseMixtureThermo.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -61,7 +62,7 @@ int main(int argc, char *argv[])
#include "readControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
- #include "createPcorrTypes.H"
+ #include "createPrghCorrTypes.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index ed5faa7c86a..030b997f486 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -44,6 +44,7 @@ Description
#include "twoPhaseMixtureThermo.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index adc380978ca..dbd859b954b 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -11,18 +11,27 @@
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
- phi = phiHbyA;
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*Dp*mesh.magSf()
+ )*rAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
tmp<fvScalarMatrix> p_rghEqnComp1;
tmp<fvScalarMatrix> p_rghEqnComp2;
@@ -70,7 +79,7 @@
fvScalarMatrix p_rghEqnIncomp
(
fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p_rgh)
+ - fvm::laplacian(rAUf, p_rgh)
);
solve
@@ -97,7 +106,7 @@
phi = phiHbyA + p_rghEqnIncomp.flux();
U = HbyA
- + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/Dp);
+ + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
index 3d1b93f429c..d1e164e3cb4 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@@ -46,6 +46,7 @@ Description
#include "fvcSmooth.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
index b819068ddd5..670ce727fda 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
@@ -45,6 +45,7 @@ Description
#include "IOMRFZoneList.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
index be0e7234938..6c315c8eb85 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -13,23 +13,32 @@
);
adjustPhi(phiHbyA, U, p_rgh);
mrfZones.makeRelative(phiHbyA);
- phi = phiHbyA;
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*Dp*mesh.magSf()
+ )*rAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - mrfZones.relative(mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -40,7 +49,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 75e4c2f8d46..4c986ef2d8a 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -41,6 +41,7 @@ Description
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -56,7 +57,7 @@ int main(int argc, char *argv[])
#include "createFields.H"
#include "createUf.H"
#include "readTimeControls.H"
- #include "createPcorrTypes.H"
+ #include "createPrghCorrTypes.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index 51190711c64..4fbd2666b4f 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -19,23 +19,31 @@
fvc::makeAbsolute(phiHbyA, U);
}
- surfaceScalarField phiAbs("phiAbs", phiHbyA);
-
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*Dp*mesh.magSf()
+ )*rAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -46,7 +54,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index dccf3f8b28d..fbcbd7c7e4a 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -45,6 +45,7 @@ Description
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
index 1fb80a4b2c9..1a10e45f3fd 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
@@ -38,6 +38,7 @@ Description
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index d8651bfa812..1194f970991 100644
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -13,23 +13,32 @@
);
adjustPhi(phiHbyA, U, p_rgh);
- phi = phiHbyA;
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*Dp*mesh.magSf()
+ )*rAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -40,7 +49,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
index 39c81c7a8dc..4e6177ef9ad 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@@ -47,6 +47,7 @@ Description
#include "IOporosityModelList.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
index 1246ee203d4..630b94534ee 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
@@ -50,6 +50,7 @@ Description
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -64,9 +65,9 @@ int main(int argc, char *argv[])
pimpleControl pimple(mesh);
#include "createFields.H"
- #include "../interFoam/interDyMFoam/createUf.H"
+ #include "createUf.H"
#include "readTimeControls.H"
- #include "../interFoam/interDyMFoam/createPcorrTypes.H"
+ #include "createPcorrTypes.H"
#include "../interFoam/interDyMFoam/correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
index 1169de5c54d..7863f5d7485 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -19,18 +19,26 @@
fvc::makeAbsolute(phiHbyA, U);
}
- surfaceScalarField phiAbs("phiAbs", phiHbyA);
-
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*Dp*mesh.magSf()
+ )*rAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
const volScalarField& vDotcP = vDotP[0]();
const volScalarField& vDotvP = vDotP[1]();
@@ -39,7 +47,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
+ fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
- (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
);
@@ -51,7 +59,7 @@
{
phi = phiHbyA + p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index bc6af055f0c..3dd2b50fed6 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -48,6 +48,7 @@ Description
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index 75f06420ee7..6714fce91c1 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -12,18 +12,27 @@
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
- phi = phiHbyA;
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*Dp*mesh.magSf()
+ )*rAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
const volScalarField& vDotcP = vDotP[0]();
const volScalarField& vDotvP = vDotP[1]();
@@ -32,7 +41,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
+ fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
- (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
);
@@ -44,7 +53,7 @@
{
phi = phiHbyA + p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
index b509b53562c..0c5c8a2a864 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
@@ -35,12 +35,12 @@ Description
#include "phaseModel.H"
#include "dragModel.H"
#include "heatTransferModel.H"
-#include "pimpleControl.H"
-
#include "singlePhaseTransportModel.H"
#include "LESModel.H"
-
+#include "pimpleControl.H"
#include "IOMRFZoneList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 00d41d3bb47..0f32807132b 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -97,6 +97,8 @@
+ rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
);
mrfZones.makeRelative(phiHbyAs[phasei]);
+ mrfZones.makeRelative(phase.phi().oldTime());
+ mrfZones.makeRelative(phase.phi());
phiHbyAs[phasei] +=
rAlphaAUfs[phasei]
@@ -157,54 +159,58 @@
phasei++;
}
- // Reset phi BCs
- phi.boundaryField() = 0;
-
- phasei = 0;
- forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
- {
- phaseModel& phase = iter();
-
- phase.phi().boundaryField() ==
- (mesh.Sf().boundaryField() & phase.U().boundaryField());
-
- mrfZones.makeRelative(phase.phi().oldTime());
- mrfZones.makeRelative(phase.phi());
-
- // Update phi BCs before pEqn
- phi.boundaryField() +=
- alphafs[phasei].boundaryField()*phase.phi().boundaryField();
-
- phasei++;
- }
-
- surfaceScalarField Dp
+ surfaceScalarField rAUf
(
IOobject
(
- "Dp",
+ "rAUf",
runTime.timeName(),
mesh
),
mesh,
- dimensionedScalar("Dp", dimensionSet(-1, 3, 1, 0, 0), 0)
+ dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
);
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
{
phaseModel& phase = iter();
- Dp += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();
+ rAUf += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();
phasei++;
}
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ {
+ surfaceScalarField::GeometricBoundaryField phib(phi.boundaryField());
+ phib = 0;
+ phasei = 0;
+ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+ {
+ phaseModel& phase = iter();
+
+ phib +=
+ alphafs[phasei].boundaryField()
+ *(mesh.Sf().boundaryField() & phase.U().boundaryField());
+
+ phasei++;
+ }
+
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p.boundaryField(),
+ (
+ phiHbyA.boundaryField() - mrfZones.relative(phib)
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+ }
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqnIncomp
(
fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rAUf, p)
);
pEqnIncomp.setReference(pRefCell, pRefValue);
@@ -221,7 +227,7 @@
if (pimple.finalNonOrthogonalIter())
{
- surfaceScalarField mSfGradp(pEqnIncomp.flux()/Dp);
+ surfaceScalarField mSfGradp(pEqnIncomp.flux()/rAUf);
phasei = 0;
phi = dimensionedScalar("phi", phi.dimensions(), 0);
@@ -248,7 +254,7 @@
// );
p.relax();
- mSfGradp = pEqnIncomp.flux()/Dp;
+ mSfGradp = pEqnIncomp.flux()/rAUf;
U = dimensionedVector("U", dimVelocity, vector::zero);
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
index 1f4b4cdb01c..bdcb8105810 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
@@ -37,6 +37,7 @@ Description
#include "turbulenceModel.H"
#include "IOMRFZoneList.H"
#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
index afa9d6058e7..12862266769 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -13,20 +13,29 @@
);
adjustPhi(phiHbyA, U, p_rgh);
mrfZones.makeRelative(phiHbyA);
- phi = phiHbyA;
surfaceScalarField phig
(
- - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
+ - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - mrfZones.relative(mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -37,7 +46,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
index 9f3fec616b3..9bfeef46745 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
@@ -36,6 +36,7 @@ Description
#include "multiphaseMixture.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index d092aab3e2e..a40c2be91b9 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -12,23 +12,32 @@
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
- phi = phiHbyA;
surfaceScalarField phig
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- )*Dp*mesh.magSf()
+ )*rAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -39,7 +48,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
index 9d5f21a94d6..fac21e8b1c8 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
@@ -20,8 +20,15 @@ adjustPhi(phiHbyA, U, p_gh);
fvOptions.makeRelative(phiHbyA);
-// Update the phi BCs from U before p BCs are updated
-phi.boundaryField() = mesh.Sf().boundaryField() & U.boundaryField();
+// Update the fixedFluxPressure BCs to ensure flux consistency
+setSnGrad<fixedFluxPressureFvPatchScalarField>
+(
+ p_gh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+);
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
index ca35e81116d..9b4caeed61e 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
@@ -39,6 +39,7 @@ Description
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index 6bc6e9a6c1a..13c3ddf7ff4 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+ surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -10,23 +10,33 @@
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
- + Dp*fvc::ddtCorr(rho, U, phi)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
- phi = phiHbyA;
surfaceScalarField phig
(
- - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
+ - ghf*fvc::snGrad(rho)*rhorAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ *rho.boundaryField()
+ )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+ );
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(Dp, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
+ fvm::laplacian(rhorAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -37,7 +47,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/settlingFoam/settlingFoam.C b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
index 57c531ba83d..3b2454c8bf5 100644
--- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C
+++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
@@ -38,6 +38,7 @@ Description
#include "plasticViscosity.H"
#include "yieldStress.H"
#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index 94495f4f30c..4a89f67d8fd 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+ surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -12,20 +12,29 @@
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
- phi = phiHbyA;
surfaceScalarField phig
(
- - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
+ - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
);
phiHbyA += phig;
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - (mesh.Sf().boundaryField() & U.boundaryField())
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -36,7 +45,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
index 8c51a402ee3..247c12a8b50 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@@ -37,6 +37,7 @@ Description
#include "incompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 6485ee03238..9e99c5cc61b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -80,22 +80,12 @@
surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
- // Update phi BCs before pEqn
- phi.boundaryField() =
- mrfZones.relative
- (
- alpha1f.boundaryField()
- *(mesh.Sf().boundaryField() & U1.boundaryField())
- + alpha2f.boundaryField()
- *(mesh.Sf().boundaryField() & U2.boundaryField())
- );
-
HbyA1 += (1.0/rho1)*rAU1*dragCoeff*U2;
HbyA2 += (1.0/rho2)*rAU2*dragCoeff*U1;
- surfaceScalarField Dp
+ surfaceScalarField rAUf
(
- "Dp",
+ "rAUf",
mag
(
alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
@@ -103,6 +93,22 @@
)
);
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p.boundaryField(),
+ (
+ phiHbyA.boundaryField()
+ - mrfZones.relative
+ (
+ alpha1f.boundaryField()
+ *(mesh.Sf().boundaryField() & U1.boundaryField())
+ + alpha2f.boundaryField()
+ *(mesh.Sf().boundaryField() & U2.boundaryField())
+ )
+ )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+ );
+
tmp<fvScalarMatrix> pEqnComp1;
tmp<fvScalarMatrix> pEqnComp2;
@@ -160,7 +166,7 @@
fvScalarMatrix pEqnIncomp
(
fvc::div(phiHbyA)
- - fvm::laplacian(Dp, p)
+ - fvm::laplacian(rAUf, p)
);
solve
@@ -175,7 +181,7 @@
if (pimple.finalNonOrthogonalIter())
{
- surfaceScalarField mSfGradp(pEqnIncomp.flux()/Dp);
+ surfaceScalarField mSfGradp(pEqnIncomp.flux()/rAUf);
phi1.boundaryField() ==
mrfZones.relative
@@ -200,7 +206,7 @@
);
p.relax();
- mSfGradp = pEqnIncomp.flux()/Dp;
+ mSfGradp = pEqnIncomp.flux()/rAUf;
U1 = HbyA1
+ fvc::reconstruct
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 7b8509337d2..1e7cbafea75 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -40,6 +40,7 @@ Description
#include "PhaseIncompressibleTurbulenceModel.H"
#include "pimpleControl.H"
#include "IOMRFZoneList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 2f0eaee518d..6d8dabcec83 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -154,7 +154,6 @@ $(derivedFvPatchFields)/mappedFixedValue/mappedFixedValueFvPatchFields.C
$(derivedFvPatchFields)/mappedFlowRate/mappedFlowRateFvPatchVectorField.C
$(derivedFvPatchFields)/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C
$(derivedFvPatchFields)/movingWallVelocity/movingWallVelocityFvPatchVectorField.C
-$(derivedFvPatchFields)/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
$(derivedFvPatchFields)/oscillatingFixedValue/oscillatingFixedValueFvPatchFields.C
$(derivedFvPatchFields)/outletInlet/outletInletFvPatchFields.C
$(derivedFvPatchFields)/outletMappedUniformInlet/outletMappedUniformInletFvPatchFields.C
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
index 12f59c1d11d..0f216379538 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
@@ -38,11 +38,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
)
:
fixedGradientFvPatchScalarField(p, iF),
- phiHbyAName_("phiHbyA"),
- phiName_("phi"),
- rhoName_("rho"),
- DpName_("Dp"),
- adjoint_(false)
+ curTimeIndex_(-1)
{}
@@ -55,11 +51,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
)
:
fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
- phiHbyAName_(ptf.phiHbyAName_),
- phiName_(ptf.phiName_),
- rhoName_(ptf.rhoName_),
- DpName_(ptf.DpName_),
- adjoint_(ptf.adjoint_)
+ curTimeIndex_(-1)
{}
@@ -71,11 +63,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
)
:
fixedGradientFvPatchScalarField(p, iF),
- phiHbyAName_(dict.lookupOrDefault<word>("phiHbyA", "phiHbyA")),
- phiName_(dict.lookupOrDefault<word>("phi", "phi")),
- rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
- DpName_(dict.lookupOrDefault<word>("Dp", "Dp")),
- adjoint_(dict.lookupOrDefault<Switch>("adjoint", false))
+ curTimeIndex_(-1)
{
if (dict.found("value") && dict.found("gradient"))
{
@@ -99,11 +87,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
)
:
fixedGradientFvPatchScalarField(wbppsf),
- phiHbyAName_(wbppsf.phiHbyAName_),
- phiName_(wbppsf.phiName_),
- rhoName_(wbppsf.rhoName_),
- DpName_(wbppsf.DpName_),
- adjoint_(wbppsf.adjoint_)
+ curTimeIndex_(-1)
{}
@@ -114,93 +98,49 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
)
:
fixedGradientFvPatchScalarField(wbppsf, iF),
- phiHbyAName_(wbppsf.phiHbyAName_),
- phiName_(wbppsf.phiName_),
- rhoName_(wbppsf.rhoName_),
- DpName_(wbppsf.DpName_),
- adjoint_(wbppsf.adjoint_)
+ curTimeIndex_(-1)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs()
+void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs
+(
+ const scalarField& snGradp
+)
{
if (updated())
{
return;
}
- const surfaceScalarField& phiHbyA =
- db().lookupObject<surfaceScalarField>(phiHbyAName_);
-
- const surfaceScalarField& phi =
- db().lookupObject<surfaceScalarField>(phiName_);
-
- fvsPatchField<scalar> phiHbyAp =
- patch().patchField<surfaceScalarField, scalar>(phiHbyA);
-
- fvsPatchField<scalar> phip =
- patch().patchField<surfaceScalarField, scalar>(phi);
-
- /*
- if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
- {
- const fvPatchField<scalar>& rhop =
- patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+ curTimeIndex_ = this->db().time().timeIndex();
- phip /= rhop;
- }
-
- if (phiHbyA.dimensions() == dimDensity*dimVelocity*dimArea)
- {
- const fvPatchField<scalar>& rhop =
- patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
- phiHbyAp /= rhop;
- }
- */
+ gradient() = snGradp;
+ fixedGradientFvPatchScalarField::updateCoeffs();
+}
- const scalarField *DppPtr = NULL;
- if (db().foundObject<volScalarField>(DpName_))
- {
- DppPtr =
- &patch().lookupPatchField<volScalarField, scalar>(DpName_);
- }
- else if (db().foundObject<surfaceScalarField>(DpName_))
+void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs()
+{
+ if (updated())
{
- const surfaceScalarField& Dp =
- db().lookupObject<surfaceScalarField>(DpName_);
-
- DppPtr =
- &patch().patchField<surfaceScalarField, scalar>(Dp);
+ return;
}
- if (adjoint_)
+ if (curTimeIndex_ != this->db().time().timeIndex())
{
- gradient() = (phip - phiHbyAp)/patch().magSf()/(*DppPtr);
+ FatalErrorIn("fixedFluxPressureFvPatchScalarField::updateCoeffs()")
+ << "updateCoeffs(const scalarField& snGradp) MUST be called before"
+ " updateCoeffs() or evaluate() to set the boundary gradient."
+ << exit(FatalError);
}
- else
- {
- gradient() = (phiHbyAp - phip)/patch().magSf()/(*DppPtr);
- }
-
- fixedGradientFvPatchScalarField::updateCoeffs();
}
void Foam::fixedFluxPressureFvPatchScalarField::write(Ostream& os) const
{
fixedGradientFvPatchScalarField::write(os);
- writeEntryIfDifferent<word>(os, "phiHbyA", "phiHbyA", phiHbyAName_);
- writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
- writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
- writeEntryIfDifferent<word>(os, "Dp", "Dp", DpName_);
- if (adjoint_)
- {
- os.writeKeyword("adjoint") << adjoint_ << token::END_STATEMENT << nl;
- }
writeEntry("value", os);
}
@@ -216,4 +156,5 @@ namespace Foam
);
}
+
// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H
index e18ef462eab..03691268f61 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H
@@ -28,44 +28,15 @@ Group
grpInletBoundaryConditions grpWallBoundaryConditions
Description
- This boundary condition adjusts the pressure gradient such that the flux
- on the boundary is that specified by the velocity boundary condition.
-
- The predicted flux to be compensated by the pressure gradient is evaluated
- as \f$(\phi - \phi_{H/A})\f$, both of which are looked-up from the database,
- as is the pressure diffusivity used to calculate the gradient using:
-
- \f[
- \nabla(p) = \frac{\phi_{H/A} - \phi}{|Sf| D_p}
- \f]
-
- where
-
- \vartable
- phi | flux
- D_p | pressure diffusivity
- Sf | patch face areas [m2]
- \endvartable
-
- \heading Patch usage
-
- \table
- Property | Description | Required | Default value
- phiHbyA | name of predicted flux field | no | phiHbyA
- phi | name of flux field | no | phi
- rho | name of density field | no | rho
- Dp | name of pressure diffusivity field | no | Dp
- \endtable
+ This boundary condition sets the pressure gradient to the provided value
+ such that the flux on the boundary is that specified by the velocity
+ boundary condition.
Example of the boundary condition specification:
\verbatim
myPatch
{
type fixedFluxPressure;
- phiHbyA phiHbyA;
- phi phi;
- rho rho;
- Dp Dp;
}
\endverbatim
@@ -99,21 +70,8 @@ class fixedFluxPressureFvPatchScalarField
{
// Private data
- //- Name of the predicted flux transporting the field
- word phiHbyAName_;
-
- //- Name of the flux transporting the field
- word phiName_;
-
- //- Name of the density field used to normalise the mass flux
- // if neccessary
- word rhoName_;
-
- //- Name of the pressure diffusivity field
- word DpName_;
-
- //- Is the pressure adjoint, i.e. has the opposite sign
- Switch adjoint_;
+ //- Current time index (used for updating)
+ label curTimeIndex_;
public:
@@ -186,7 +144,10 @@ public:
// Member functions
- //- Update the coefficients associated with the patch field
+ //- Update the patch pressure gradient field from the given snGradp
+ virtual void updateCoeffs(const scalarField& snGradp);
+
+ //- Update the patch pressure gradient field
virtual void updateCoeffs();
//- Write
@@ -200,6 +161,39 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include "volFields.H"
+
+namespace Foam
+{
+ template<class GradBC>
+ inline void setSnGrad
+ (
+ volScalarField::GeometricBoundaryField& bf,
+ const FieldField<fvsPatchField, double>& snGrad
+ )
+ {
+ forAll(bf, patchi)
+ {
+ if (isA<GradBC>(bf[patchi]))
+ {
+ refCast<GradBC>(bf[patchi]).updateCoeffs(snGrad[patchi]);
+ }
+ }
+ }
+
+ template<class GradBC>
+ inline void setSnGrad
+ (
+ volScalarField::GeometricBoundaryField& bf,
+ const tmp<FieldField<fvsPatchField, double> >& tsnGrad
+ )
+ {
+ setSnGrad<GradBC>(bf, tsnGrad());
+ }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
#endif
// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
deleted file mode 100644
index 89cf73e8ad7..00000000000
--- a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
+++ /dev/null
@@ -1,197 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "multiphaseFixedFluxPressureFvPatchScalarField.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
- const fvPatch& p,
- const DimensionedField<scalar, volMesh>& iF
-)
-:
- fixedGradientFvPatchScalarField(p, iF),
- phiHbyAName_("phiHbyA"),
- phiName_("phi"),
- rhoName_("rho"),
- DpName_("Dp")
-{}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
- const multiphaseFixedFluxPressureFvPatchScalarField& ptf,
- const fvPatch& p,
- const DimensionedField<scalar, volMesh>& iF,
- const fvPatchFieldMapper& mapper
-)
-:
- fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
- phiHbyAName_(ptf.phiHbyAName_),
- phiName_(ptf.phiName_),
- rhoName_(ptf.rhoName_),
- DpName_(ptf.DpName_)
-{}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
- const fvPatch& p,
- const DimensionedField<scalar, volMesh>& iF,
- const dictionary& dict
-)
-:
- fixedGradientFvPatchScalarField(p, iF),
- phiHbyAName_(dict.lookupOrDefault<word>("phiHbyA", "phiHbyA")),
- phiName_(dict.lookupOrDefault<word>("phi", "phi")),
- rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
- DpName_(dict.lookupOrDefault<word>("Dp", "Dp"))
-{
- if (dict.found("gradient"))
- {
- gradient() = scalarField("gradient", dict, p.size());
- fixedGradientFvPatchScalarField::updateCoeffs();
- fixedGradientFvPatchScalarField::evaluate();
- }
- else
- {
- fvPatchField<scalar>::operator=(patchInternalField());
- gradient() = 0.0;
- }
-}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
- const multiphaseFixedFluxPressureFvPatchScalarField& mfppsf
-)
-:
- fixedGradientFvPatchScalarField(mfppsf),
- phiHbyAName_(mfppsf.phiHbyAName_),
- phiName_(mfppsf.phiName_),
- rhoName_(mfppsf.rhoName_),
- DpName_(mfppsf.DpName_)
-{}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
- const multiphaseFixedFluxPressureFvPatchScalarField& mfppsf,
- const DimensionedField<scalar, volMesh>& iF
-)
-:
- fixedGradientFvPatchScalarField(mfppsf, iF),
- phiHbyAName_(mfppsf.phiHbyAName_),
- phiName_(mfppsf.phiName_),
- rhoName_(mfppsf.rhoName_),
- DpName_(mfppsf.DpName_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs()
-{
- if (updated())
- {
- return;
- }
-
- const surfaceScalarField& phiHbyA =
- db().lookupObject<surfaceScalarField>(phiHbyAName_);
-
- const surfaceScalarField& phi =
- db().lookupObject<surfaceScalarField>(phiName_);
-
- fvsPatchField<scalar> phiHbyAp =
- patch().patchField<surfaceScalarField, scalar>(phiHbyA);
-
- fvsPatchField<scalar> phip =
- patch().patchField<surfaceScalarField, scalar>(phi);
-
- /*
- if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
- {
- const fvPatchField<scalar>& rhop =
- patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
- phip /= rhop;
- }
-
- if (phiHbyA.dimensions() == dimDensity*dimVelocity*dimArea)
- {
- const fvPatchField<scalar>& rhop =
- patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
- phiHbyAp /= rhop;
- }
- */
-
- const fvsPatchField<scalar>& Dpp =
- patch().lookupPatchField<surfaceScalarField, scalar>(DpName_);
-
- gradient() = (phiHbyAp - phip)/patch().magSf()/Dpp;
-
- fixedGradientFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::multiphaseFixedFluxPressureFvPatchScalarField::write
-(
- Ostream& os
-) const
-{
- fvPatchScalarField::write(os);
- writeEntryIfDifferent<word>(os, "phiHbyA", "phiHbyA", phiHbyAName_);
- writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
- writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
- writeEntryIfDifferent<word>(os, "Dp", "Dp", DpName_);
- gradient().writeEntry("gradient", os);
- writeEntry("value", os);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
- makePatchTypeField
- (
- fvPatchScalarField,
- multiphaseFixedFluxPressureFvPatchScalarField
- );
-}
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H
deleted file mode 100644
index 18b9698c0e1..00000000000
--- a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H
+++ /dev/null
@@ -1,199 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
-
-Class
- Foam::multiphaseFixedFluxPressureFvPatchScalarField
-
-Group
- grpWallBoundaryConditions grpGenericBoundaryConditions
-
-Description
- This boundary condition adjusts the pressure gradient such that the flux
- on the boundary is that specified by the velocity boundary condition.
-
- The predicted flux to be compensated by the pressure gradient is evaluated
- as \f$(\phi - \phi_{H/A})\f$, both of which are looked-up from the database,
- as is the pressure diffusivity Dp used to calculate the gradient using:
-
- \f[
- \nabla(p) = \frac{\phi_{H/A} - \phi}{|Sf| Dp}
- \f]
-
- where
-
- \vartable
- \phi | flux
- Dp | pressure diffusivity
- Sf | patch face areas [m2]
- \endvartable
-
- \heading Patch usage
-
- \table
- Property | Description | Required | Default value
- phiHbyA | name of predicted flux field | no | phiHbyA
- phi | name of flux field | no | phi
- rho | name of density field | no | rho
- Dp | name of pressure diffusivity field | no | Dp
- \endtable
-
- Example of the boundary condition specification:
- \verbatim
- myPatch
- {
- type multiphaseFixedFluxPressure;
- phiHbyA phiHbyA;
- phi phi;
- rho rho;
- Dp Dp;
- }
- \endverbatim
-
-SourceFiles
- multiphaseFixedFluxPressureFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef multiphaseFixedFluxPressureFvPatchScalarFields_H
-#define multiphaseFixedFluxPressureFvPatchScalarFields_H
-
-#include "fvPatchFields.H"
-#include "fixedGradientFvPatchFields.H"
-#include "Switch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
- Class multiphaseFixedFluxPressureFvPatch Declaration
-\*---------------------------------------------------------------------------*/
-
-class multiphaseFixedFluxPressureFvPatchScalarField
-:
- public fixedGradientFvPatchScalarField
-{
- // Private data
-
- //- Name of the predicted flux transporting the field
- word phiHbyAName_;
-
- //- Name of the flux transporting the field
- word phiName_;
-
- //- Name of the density field used to normalise the mass flux
- // if neccessary
- word rhoName_;
-
- //- Name of the pressure diffusivity field
- word DpName_;
-
-
-public:
-
- //- Runtime type information
- TypeName("multiphaseFixedFluxPressure");
-
-
- // Constructors
-
- //- Construct from patch and internal field
- multiphaseFixedFluxPressureFvPatchScalarField
- (
- const fvPatch&,
- const DimensionedField<scalar, volMesh>&
- );
-
- //- Construct from patch, internal field and dictionary
- multiphaseFixedFluxPressureFvPatchScalarField
- (
- const fvPatch&,
- const DimensionedField<scalar, volMesh>&,
- const dictionary&
- );
-
- //- Construct by mapping given
- // multiphaseFixedFluxPressureFvPatchScalarField onto a new patch
- multiphaseFixedFluxPressureFvPatchScalarField
- (
- const multiphaseFixedFluxPressureFvPatchScalarField&,
- const fvPatch&,
- const DimensionedField<scalar, volMesh>&,
- const fvPatchFieldMapper&
- );
-
- //- Construct as copy
- multiphaseFixedFluxPressureFvPatchScalarField
- (
- const multiphaseFixedFluxPressureFvPatchScalarField&
- );
-
- //- Construct and return a clone
- virtual tmp<fvPatchScalarField> clone() const
- {
- return tmp<fvPatchScalarField>
- (
- new multiphaseFixedFluxPressureFvPatchScalarField(*this)
- );
- }
-
- //- Construct as copy setting internal field reference
- multiphaseFixedFluxPressureFvPatchScalarField
- (
- const multiphaseFixedFluxPressureFvPatchScalarField&,
- const DimensionedField<scalar, volMesh>&
- );
-
- //- Construct and return a clone setting internal field reference
- virtual tmp<fvPatchScalarField> clone
- (
- const DimensionedField<scalar, volMesh>& iF
- ) const
- {
- return tmp<fvPatchScalarField>
- (
- new multiphaseFixedFluxPressureFvPatchScalarField(*this, iF)
- );
- }
-
-
- // Member functions
-
- //- Update the coefficients associated with the patch field
- virtual void updateCoeffs();
-
- //- Write
- virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/0/p_rgh.org b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/0/p_rgh.org
index 396d380dadd..8bd9793d9af 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/0/p_rgh.org
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/0/p_rgh.org
@@ -23,7 +23,6 @@ boundaryField
walls
{
type fixedFluxPressure;
- value uniform 1e5;
}
defaultFaces
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh
index 85bfb4aef14..1b9491f31ea 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh
+++ b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh
@@ -30,14 +30,3 @@ boundaryField
p0 uniform 0;
value uniform 0;
}
- walls
- {
- type fixedFluxPressure;
- gradient uniform 0;
- phi phiAbs;
- value uniform 0;
- }
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/0/p_rgh b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/0/p_rgh
index dafcfc526a1..66122dcff0e 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/0/p_rgh
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/0/p_rgh
@@ -31,8 +31,6 @@ boundaryField
walls
{
type fixedFluxPressure;
- phi phiAbs;
- value uniform 0;
}
}
--
GitLab