diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index c29d45109150fbcfb45728cdf3bd6ce93f619157..865af56184bf5f5baff6ce1732dad1a61d967de0 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 83a44ab7c2284f5c60dc66a6bd6714d16464dd09..995a59a737e66fca4ba33ac70b917af8f3d69d4d 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 a92e7285ab2cef90d59b7f9aed3f165b5b5ecfd4..384ecbee598a7f03806862c4c74645f1ba142c68 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 14e46bef66199cc2a28d6e1ab8f6e97a988952e4..b8c48549a68bf3b652e49fb3f8f3716f50b130a8 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 c98ac61f3b491d48df16df2c3a015e732d80c550..6652eab56495d5e1ea25a4544d944a85ed61206b 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 0dc5d422db91a34d183d06e3a2b63afe80bcb8b8..f7a3004ac212d83ab7ebfd4ce4f877e4d4d5d6f6 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 8602725bdb47d3885b91804bf4b8768e0792202c..8894c7a9d766b783292a22de3ff602158d7be028 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 19e9f96f512a871c54ac36ad3d37bc83e52bcf72..fd4b6f58eb5410fc0cc042e2b71011c9099dae2a 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 3c92ce8119733f0aa2c6de45e86b147d2894b900..f9b1bd65e41ea05198f2ee77603b2eaaaae6155f 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 373bafcd3b85bedfadac2f5b0adf31dd3da761b5..a4e5fa20da7f13b19caac7322160a0364248d29d 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 d2681663e634b33de6d1c6da211417b28f9a9ffc..1d70a0f354b1ab8a11fdd5ecd6978e0a3effb5d6 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 f8837f0aa4d532d742dfa5c4f211235565a76c70..0503db4207f7d64ef8dfb12b017c30a1076b201b 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 d61a80ffa8114e199e731494f98a1f09cb74da3a..c656b3c3d78e0273e6605733cc5648e98ef6c316 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 2832bf421a989214922e0a9aacc2a3984a98abed..edde4b198e579f7783e781b14b30678ca5e7d906 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 91becd8aa16fdf5ed7203c3e6a8eb80a83ff2f7c..f8f99dafa545fa0758ce83ece4778e51083b0a39 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 2a0b08074a31260b5e4a81e3c1b9caf8fd36a4be..df3f3cdc61f2747dc461022a82b71cd37f273809 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 945aab6ba32a60642f1c414f0e28a3df35167c67..e9dedef5a93078456b8725cb62cf47bb488851f5 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 6168082e1de26c7b9f055685ebac0d9348036ec7..43d8a4da534b4978ca6571bc736a235e1c52bb4f 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 94299b84af182f22d17fc30b11de8581cfa40ea3..d77eb12b0ee1d1fa19fc4a695fc3b370c9359194 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 b0d2f98aebf66a1d69d0f55779b09c0560f76dc6..483a723f7cedad96bb4865a341b9fa80cf7d6ac7 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 53fec2e1b64e9944244e7bbf9d61ee37ca036625..bae1bcef886ab4eca5d03b6c9035c048bff935f5 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 b1b6ceede2604aff349e25f489705b1b02e4a27f..9d62893f18ae4a59a54596313c14dc0321bff2c9 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 ad98545b4076f89b412a0c33c992523a28c24c2c..5f6d70e054bd44eea65aaf089072376148349d99 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 67df5b61e054ef46b32c8daefe3e367023496248..2e52da8f5310ee2edcfe1a56c5979036bdf24727 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 41c8900dc086833f35f3344ea10ea356dde21caa..60b3c1a0cfd524a657cae49185a179b091657fa3 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 b2cd9ac54e7f208b666fbf528b17e319dcee5bfc..6aae0cf5b68eee3598da96ede3b920fc07d37926 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 a3615c4309ad487272c631ce96d2c81a9834e0f8..75720dd51c28d1f8a9f1cf9ca0c21bd0b7d3286e 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 f42ed7e39ceaf313a3decd035ac5e1ac60dddf78..fc854c40a5019eb34709d5997814f19433c0de01 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 a712250f35646fd393a0a33d02087fab66a1014f..38ba49bd74799aa78a81265dc03017326f3c8852 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 6824d8ee0c7564f543e10f8f741106e4595d87c8..4799df09f72a7b5b1052c6ce5b7fb217f2c0b40e 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 827daa4ce8cf5b242205b1470b5e5d6637daafbb..7ecfe1f66236e7458835eb35dde5a76a7038b573 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 200ea289d857b60c111f6d7b503c66b3fe3cf59a..fce4b6e64da208836fbe0d09d78d05797a22a287 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 1e504462b6a7a786d041dc2d692a733d44ff8a77..6652eab56495d5e1ea25a4544d944a85ed61206b 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 4067c125ebecbe9ce5eda92ec1916ecd534e2478..27f272100af795893e3c5f3740c140565fc39462 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 420bc1ef0cb12b2ffc51ed2f51f1d83bce190b75..56b22dba9caab96e4769eae137c6100d29a49466 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 f62a17286c2c4d52b1f6941c0fe669474e92b999..26511b001df3e52986ce6bd2e5244ce582c56aca 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 b7de7e2f79ec2f5eb9d5b468fab01b7972fd9c1c..62c3e0376e363ef3f9a8d9a459d33bf10773dc33 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 f869961a3c3161da03ec9b1915a90b9c87882948..e7f078d1a6578be726065c128ede11fdacf66a44 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 e9d059aadbc47f951f494e7774c73f1dacd0e0f6..210a30fd14371e7983722106af70c3697eed66fa 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 7590e6ff45c493e3d0a3be2da8decbe9fd62cd37..a9c3365bdf73fedfa445979e9cf7e5b2f2ebf937 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 ed5faa7c86a515ad69565b24d1269267b1e2ade1..030b997f486bf0c15ccf4c86531f16407d24cdf6 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 adc380978ca94f57c7176be91356662deb998504..dbd859b954b7da6f4ed8fb154a39f852e449f47c 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 3d1b93f429c460c6987f6260dc3a492d8b8f0926..d1e164e3cb4f5c869a95b10fe32d68df6755ec6b 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 b819068ddd53f56d6af4456906466158eace7f51..670ce727fda8e3457badeb753564cb33eafd2c97 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 be0e7234938396f9369419fa64c7a258bdc5c0bd..6c315c8eb85421507cf8899fa7b6ca5e77225f73 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 75e4c2f8d46f69d15b73236f1d677d484c9983ce..4c986ef2d8ab48263f7f520814aabded1a11df2d 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 51190711c647047450cdfe0df0f181107953a6f9..4fbd2666b4ffd108ba263561c6076e4884fb69fa 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 dccf3f8b28dbcc1ed1c2a8d3668619818ce5a24c..fbcbd7c7e4a34864a96b1931a35ceb1ccfcce9e4 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 1fb80a4b2c9dd40a8a2349dd9be4e973eeab02af..1a10e45f3fda491e7fb8495a2fda79e389aee125 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 d8651bfa8120b0f368035f89570d942fe7d8a86c..1194f9709916f358b912b4874f124d573ccf5252 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 39c81c7a8dc9fb58cc2bda021b6ed73177803688..4e6177ef9ad4680fef0369f685416f7d581af2ea 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 1246ee203d4b91ccd0330895c7c24748eab28cc8..630b94534eec7a95298de13cda999fa07b8b4ebb 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 1169de5c54d90a403c51be15f9a37574a7d95fd4..7863f5d74851ebba806897821174ff43b6484bea 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 bc6af055f0c7e39e904101bf8695299e4dbd311b..3dd2b50fed60412d4284d3e542ff3ce2d07023fc 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 75f06420ee7df104c07e546d74bdf6273b47015a..6714fce91c18c54c453c0ad43b636255278d3a79 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 b509b53562cf63daee243cf275856c79e549c19d..0c5c8a2a86440d8c6afd828fd1bd95b40919adf9 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 00d41d3bb476672177ccb86896788bf5bcdeb6b6..0f32807132bf091891142a8bb4ac514495b65cc2 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 1f4b4cdb01c8b2cb41a85e092bc3a8772c0269d8..bdcb8105810cca9ca8e9046ef443fb3cf6503264 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 afa9d6058e7a2156535d9b5398a44ff9303c0440..128622667691859179b0cc34c34eaa357580d2fc 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 9f3fec616b3fcea27a13b06b7b45e2d695705aa4..9bfeef4674584a5c3fc6f3a0d30b54a524a31e8e 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 d092aab3e2ee8cb0b6813d7dbec44b64ee2b842f..a40c2be91b917517513202a15c2bf2eb5ad843ad 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 9d5f21a94d609809f88436eab024079983d6a059..fac21e8b1c84a0732a8ef0480e0568edfb3fab73 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 ca35e81116db2bae461b4029efc65037cebc6691..9b4caeed61ed0149dc13d98d02e9b493dc7489c9 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 6bc6e9a6c1a7a7cfb5c9b6ccafaccbbb5afbbed9..13c3ddf7ff4a306ef4a5159108965fb19bc1d0f2 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 57c531ba83d8f3fc2ecf9761a92c738343dd5e58..3b2454c8bf5f0c4d57f3e35349d715225ad8e1b3 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 94495f4f30cca92a1b5b5ef65febb7ed467275cb..4a89f67d8fde38188622defe6cdb80f9f31f08cf 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 8c51a402ee377d49a61b1429fbcfc4c7ddb34afc..247c12a8b508e14c431437da31851532c824ead2 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 6485ee032387e9995b741f69954c68bbfbe92b0b..9e99c5cc61bfd9daec53705af3f31add2ec095ab 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 7b8509337d29d592c1cc961bab971ffd111b9d0a..1e7cbafea75b13045ce2df010d2a4f164bb7f3ad 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 2f0eaee518d297e869217d3862052a874fd87ba8..6d8dabcec837cae6191202df20d133a7102bf14a 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 12f59c1d11dc217c4bf9eda0f29e20f027904f46..0f21637953813fc8188241451707269e69d51797 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 e18ef462eab6df881717561c60a127fc2ff16f3b..03691268f6141703d467da09c2843845ced53d89 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 89cf73e8ad7d4ea8de62c4958bc78d43b88ba647..0000000000000000000000000000000000000000
--- 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 18b9698c0e1e70092fe0c0aedd116dad44e1b046..0000000000000000000000000000000000000000
--- 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 396d380daddd7bc4e4a3c9d427a5226c6c3f44af..8bd9793d9af294b49a3ed90905a2c751c769770f 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 85bfb4aef1475eac599bd35b72b09081d918d929..1b9491f31eaf3e549e37e506d6f04eb8eedb2b5c 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 dafcfc526a1009b4eb35db6050301abce604f11a..66122dcff0e940326df3740c0bf86c718991a247 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;
}
}