From 92c4c12a08efd2f49e338e806731efde4a4ef671 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Tue, 3 May 2016 15:51:15 +0100
Subject: [PATCH] Multiphase solvers: Update p_rgh following density changes
---
.../compressibleInterDyMFoam/compressibleInterDyMFoam.C | 4 ++++
.../compressibleInterFoam/compressibleInterDyMFoam/pEqn.H | 4 ++++
.../solvers/multiphase/compressibleInterFoam/pEqn.H | 6 ++++--
.../multiphase/compressibleMultiphaseInterFoam/pEqn.H | 7 ++++---
.../reactingMultiphaseEulerFoam/pU/pEqn.H | 8 +++++---
.../reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H | 1 +
.../reactingTwoPhaseEulerFoam/pUf/pEqn.H | 1 +
.../solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H | 8 +++++---
.../solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H | 6 +++++-
9 files changed, 33 insertions(+), 12 deletions(-)
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 65cc137acb..521628fa10 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -146,6 +146,10 @@ int main(int argc, char *argv[])
rho = alpha1*rho1 + alpha2*rho2;
+ // Correct p_rgh for consistency with p and the updated densities
+ p_rgh = p - rho*gh;
+ p_rgh.correctBoundaryConditions();
+
runTime.write();
Info<< "ExecutionTime = "
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index 859fc9cc4e..fef73aa365 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -116,6 +116,10 @@
rho = alpha1*rho1 + alpha2*rho2;
+ // Correct p_rgh for consistency with p and the updated densities
+ p_rgh = p - rho*gh;
+ p_rgh.correctBoundaryConditions();
+
K = 0.5*magSqr(U);
Info<< "max(U) " << max(mag(U)).value() << endl;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index a4a3fb2e29..78e46f4664 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -101,14 +101,16 @@
}
}
- // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-
// Update densities from change in p_rgh
rho1 += psi1*(p_rgh - p_rgh_0);
rho2 += psi2*(p_rgh - p_rgh_0);
rho = alpha1*rho1 + alpha2*rho2;
+ // Correct p_rgh for consistency with p and the updated densities
+ p_rgh = p - rho*gh;
+ p_rgh.correctBoundaryConditions();
+
K = 0.5*magSqr(U);
Info<< "max(U) " << max(mag(U)).value() << endl;
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
index 70ac37c1b8..6526dacb58 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
@@ -95,9 +95,6 @@
if (pimple.finalNonOrthogonalIter())
{
- // p = max(p_rgh + mixture.rho()*gh, pMin);
- // p_rgh = p - mixture.rho()*gh;
-
phasei = 0;
forAllIter
(
@@ -125,6 +122,10 @@
mixture.correctRho(p_rgh - p_rgh_0);
rho = mixture.rho();
+ // Correct p_rgh for consistency with p and the updated densities
+ p_rgh = p - rho*gh;
+ p_rgh.correctBoundaryConditions();
+
K = 0.5*magSqr(U);
Info<< "max(U) " << max(mag(U)).value() << endl;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index 3431a1db24..7e43e67749 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -92,6 +92,11 @@ while (pimple.correct())
// Update continuity errors due to temperature changes
fluid.correct();
+ volScalarField rho("rho", fluid.rho());
+
+ // Correct p_rgh for consistency with p and the updated densities
+ p_rgh = p - rho*gh;
+
PtrList<volVectorField> HbyAs(phases.size());
forAll(phases, phasei)
@@ -121,9 +126,6 @@ while (pimple.correct())
);
}
- // Mean density for buoyancy force and p_rgh -> p
- volScalarField rho("rho", fluid.rho());
-
surfaceScalarField ghSnGradRho
(
"ghSnGradRho",
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 6922869aa2..9020cd6d5d 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -85,6 +85,7 @@ while (pimple.correct())
{
// Update continuity errors due to temperature changes
fluid.correct();
+
volScalarField rho("rho", fluid.rho());
// Correct p_rgh for consistency with p and the updated densities
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index ddc262ba78..b27c9f4646 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -95,6 +95,7 @@ while (pimple.correct())
{
// Update continuity errors due to temperature changes
fluid.correct();
+
volScalarField rho("rho", fluid.rho());
// Correct p_rgh for consistency with p and the updated densities
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
index a64aeea388..2010e65ebc 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
@@ -83,6 +83,11 @@ while (pimple.correct())
// Update continuity errors due to temperature changes
#include "correctContErrs.H"
+ volScalarField rho("rho", fluid.rho());
+
+ // Correct p_rgh for consistency with p and the updated densities
+ p_rgh = p - rho*gh;
+
// Correct fixed-flux BCs to be consistent with the velocity BCs
MRF.correctBoundaryFlux(U1, phi1);
MRF.correctBoundaryFlux(U2, phi2);
@@ -113,9 +118,6 @@ while (pimple.correct())
*rho2*U2.oldTime()/runTime.deltaT()
);
- // Mean density for buoyancy force and p_rgh -> p
- volScalarField rho("rho", fluid.rho());
-
surfaceScalarField ghSnGradRho
(
"ghSnGradRho",
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
index bb78b2799d..f8388b2d45 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
@@ -95,6 +95,11 @@ while (pimple.correct())
// Update continuity errors due to temperature changes
#include "correctContErrs.H"
+ volScalarField rho("rho", fluid.rho());
+
+ // Correct p_rgh for consistency with p and the updated densities
+ p_rgh = p - rho*gh;
+
surfaceScalarField rhof1(fvc::interpolate(rho1));
surfaceScalarField rhof2(fvc::interpolate(rho2));
@@ -114,7 +119,6 @@ while (pimple.correct())
max(alphaf2, phase2.residualAlpha())*rAUf2
);
- volScalarField rho("rho", fluid.rho());
surfaceScalarField ghSnGradRho
(
"ghSnGradRho",
--
GitLab