reactingTwoPhaseEulerFoam: Update p_rgh following density changes

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=2073
This commit is contained in:
Henry Weller
2016-05-03 14:53:11 +01:00
parent 698f135426
commit 055b41e4da
2 changed files with 8 additions and 10 deletions

View File

@ -85,6 +85,10 @@ while (pimple.correct())
{ {
// Update continuity errors due to temperature changes // Update continuity errors due to temperature changes
fluid.correct(); fluid.correct();
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 // Correct fixed-flux BCs to be consistent with the velocity BCs
MRF.correctBoundaryFlux(U1, phi1); MRF.correctBoundaryFlux(U1, phi1);
@ -116,9 +120,6 @@ while (pimple.correct())
*rho2*U2.oldTime()/runTime.deltaT() *rho2*U2.oldTime()/runTime.deltaT()
); );
// Mean density for buoyancy force and p_rgh -> p
volScalarField rho("rho", fluid.rho());
surfaceScalarField ghSnGradRho surfaceScalarField ghSnGradRho
( (
"ghSnGradRho", "ghSnGradRho",

View File

@ -95,6 +95,10 @@ while (pimple.correct())
{ {
// Update continuity errors due to temperature changes // Update continuity errors due to temperature changes
fluid.correct(); fluid.correct();
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 rhof1(fvc::interpolate(rho1));
surfaceScalarField rhof2(fvc::interpolate(rho2)); surfaceScalarField rhof2(fvc::interpolate(rho2));
@ -115,7 +119,6 @@ while (pimple.correct())
max(alphaf2, phase2.residualAlpha())*rAUf2 max(alphaf2, phase2.residualAlpha())*rAUf2
); );
volScalarField rho("rho", fluid.rho());
surfaceScalarField ghSnGradRho surfaceScalarField ghSnGradRho
( (
"ghSnGradRho", "ghSnGradRho",
@ -389,12 +392,6 @@ while (pimple.correct())
} }
} }
Info<< "min(p) = " << min(p_rgh + rho*gh).value() << endl;
if (min(p_rgh + rho*gh) < pMin)
{
Info<< "Clipping p" << endl;
}
// Update and limit the static pressure // Update and limit the static pressure
p = max(p_rgh + rho*gh, pMin); p = max(p_rgh + rho*gh, pMin);