diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 15b97633166a84306e64b3198952ae08cea0db1a..9e015406d187ba134a2c5f709bdee986ef91ae18 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -1,3 +1,5 @@
+// --- Pressure corrector loop
+while (pimple.correct())
{
surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1));
surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqnElim.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqnElim.H
new file mode 100644
index 0000000000000000000000000000000000000000..638dc11b8686514086e8107843ce957cfc09a67e
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqnElim.H
@@ -0,0 +1,339 @@
+surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1));
+surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f);
+
+rAU1 = 1.0/U1Eqn.A();
+rAU2 = 1.0/U2Eqn.A();
+
+surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rho1*rAU1));
+surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rho2*rAU2));
+
+// --- Pressure corrector loop
+while (pimple.correct())
+{
+ volVectorField HbyA1
+ (
+ IOobject::groupName("HbyA", phase1.name()),
+ U1
+ );
+ HbyA1 = rAU1*U1Eqn.H();
+
+ volVectorField HbyA2
+ (
+ IOobject::groupName("HbyA", phase2.name()),
+ U2
+ );
+ HbyA2 = rAU2*U2Eqn.H();
+
+ // Phase pressure flux (e.g. due to particle-particle pressure)
+ surfaceScalarField phiF1
+ (
+ "phiF1",
+ fvc::interpolate(rAU1*phase1.turbulence().pPrime())
+ *fvc::snGrad(alpha1)*mesh.magSf()
+ );
+
+ // Phase-2 pressure flux (e.g. due to particle-particle pressure)
+ surfaceScalarField phiF2
+ (
+ "phiF2",
+ fvc::interpolate(rAU2*phase2.turbulence().pPrime())
+ *fvc::snGrad(alpha2)*mesh.magSf()
+ );
+
+ // Mean density for buoyancy force and p_rgh -> p
+ volScalarField rho("rho", fluid.rho());
+
+ // Add Buoyancy force
+ {
+ surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
+
+ phiF1 +=
+ rAlphaAU1f
+ *(
+ ghSnGradRho
+ - alpha2f*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
+ )/fvc::interpolate(rho1);
+
+ phiF2 +=
+ rAlphaAU2f
+ *(
+ ghSnGradRho
+ - alpha1f*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
+ )/fvc::interpolate(rho2);
+ }
+
+ // ddtPhiCorr filter -- only apply in pure(ish) phases
+ surfaceScalarField alpha1fBar(fvc::interpolate(fvc::average(alpha1f)));
+ surfaceScalarField phiCorrCoeff1(pos(alpha1fBar - 0.99));
+ surfaceScalarField phiCorrCoeff2(pos(0.01 - alpha1fBar));
+
+ // Set ddtPhiCorr to 0 on non-coupled boundaries
+ forAll(mesh.boundary(), patchi)
+ {
+ if
+ (
+ !mesh.boundary()[patchi].coupled()
+ || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
+ )
+ {
+ phiCorrCoeff1.boundaryField()[patchi] = 0;
+ phiCorrCoeff2.boundaryField()[patchi] = 0;
+ }
+ }
+
+ // Phase-1 predicted flux
+ surfaceScalarField phiHbyA1
+ (
+ IOobject::groupName("phiHbyA", phase1.name()),
+ (fvc::interpolate(HbyA1) & mesh.Sf())
+ + phiCorrCoeff1*rAlphaAU1f
+ *(
+ mrfZones.absolute(phi1.oldTime())
+ - (fvc::interpolate(U1.oldTime()) & mesh.Sf())
+ )/runTime.deltaT()
+ - phiF1
+ );
+
+ // Phase-2 predicted flux
+ surfaceScalarField phiHbyA2
+ (
+ IOobject::groupName("phiHbyA", phase2.name()),
+ (fvc::interpolate(HbyA2) & mesh.Sf())
+ + phiCorrCoeff2*rAlphaAU2f
+ *(
+ mrfZones.absolute(phi2.oldTime())
+ - (fvc::interpolate(U2.oldTime()) & mesh.Sf())
+ )/runTime.deltaT()
+ - phiF2
+ );
+
+ // Face-drag coefficients
+ surfaceScalarField D1f(fvc::interpolate(rAU1*dragCoeff));
+ surfaceScalarField D2f(fvc::interpolate(rAU2*dragCoeff));
+
+ // Construct the mean predicted flux
+ // including explicit drag contributions based on absolute fluxes
+ surfaceScalarField phiHbyA
+ (
+ "phiHbyA",
+ alpha1f*(phiHbyA1 + D1f*mrfZones.absolute(phi2))
+ + alpha2f*(phiHbyA2 + D2f*mrfZones.absolute(phi1))
+ );
+ mrfZones.makeRelative(phiHbyA);
+
+ // Construct pressure "diffusivity"
+ surfaceScalarField rAUf
+ (
+ "rAUf",
+ mag
+ (
+ alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
+ + alpha2f*rAlphaAU2f/fvc::interpolate(rho2)
+ )
+ );
+
+ // Update the fixedFluxPressure BCs to ensure flux consistency
+ setSnGrad<fixedFluxPressureFvPatchScalarField>
+ (
+ p_rgh.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;
+
+ // Construct the compressibility parts of the pressure equation
+ if (pimple.transonic())
+ {
+ surfaceScalarField phid1
+ (
+ IOobject::groupName("phid", phase1.name()),
+ fvc::interpolate(psi1)*phi1
+ );
+ surfaceScalarField phid2
+ (
+ IOobject::groupName("phid", phase2.name()),
+ fvc::interpolate(psi2)*phi2
+ );
+
+ pEqnComp1 =
+ (
+ contErr1
+ - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+ )/rho1
+ + (alpha1/rho1)*correction
+ (
+ psi1*fvm::ddt(p_rgh)
+ + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+ );
+ deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
+ pEqnComp1().relax();
+
+ pEqnComp2 =
+ (
+ contErr2
+ - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+ )/rho2
+ + (alpha2/rho2)*correction
+ (
+ psi2*fvm::ddt(p_rgh)
+ + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+ );
+ deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
+ pEqnComp2().relax();
+ }
+ else
+ {
+ pEqnComp1 =
+ (
+ contErr1
+ - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+ )/rho1
+ + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
+
+ pEqnComp2 =
+ (
+ contErr2
+ - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+ )/rho2
+ + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
+ }
+
+ // Cache p prior to solve for density update
+ volScalarField p_rgh_0(p_rgh);
+
+ // Iterate over the pressure equation to correct for non-orthogonality
+ while (pimple.correctNonOrthogonal())
+ {
+ // Construct the transport part of the pressure equation
+ fvScalarMatrix pEqnIncomp
+ (
+ fvc::div(phiHbyA)
+ - fvm::laplacian(rAUf, p_rgh)
+ );
+
+ solve
+ (
+ pEqnComp1() + pEqnComp2() + pEqnIncomp,
+ mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+ );
+
+ // Correct fluxes and velocities on last non-orthogonal iteration
+ if (pimple.finalNonOrthogonalIter())
+ {
+ phi = phiHbyA + pEqnIncomp.flux();
+
+ surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
+
+ // Partial-elimination phase-flux corrector
+ {
+ surfaceScalarField phi1s
+ (
+ phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1)
+ );
+
+ surfaceScalarField phi2s
+ (
+ phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2)
+ );
+
+ surfaceScalarField phir
+ (
+ ((phi1s + D1f*phi2s) - (phi2s + D2f*phi1s))/(1 - D1f*D2f)
+ );
+
+ phi1.boundaryField() ==
+ mrfZones.relative
+ (
+ mesh.Sf().boundaryField() & U1.boundaryField()
+ );
+ phi1 = phi + alpha2f*phir;
+
+ phi2.boundaryField() ==
+ mrfZones.relative
+ (
+ mesh.Sf().boundaryField() & U2.boundaryField()
+ );
+ phi2 = phi - alpha1f*phir;
+ }
+
+ // Compressibility correction for phase-fraction equations
+ fluid.dgdt() =
+ (
+ alpha1*(pEqnComp2 & p_rgh)
+ - alpha2*(pEqnComp1 & p_rgh)
+ );
+
+ // Optionally relax pressure for velocity correction
+ p_rgh.relax();
+
+ // Update the static pressure
+ p = max(p_rgh + rho*gh, pMin);
+ p_rgh = p - rho*gh;
+
+ mSfGradp = pEqnIncomp.flux()/rAUf;
+
+ // Partial-elimination phase-velocity corrector
+ {
+ volVectorField U1s
+ (
+ HbyA1
+ + fvc::reconstruct
+ (
+ rAlphaAU1f*mSfGradp/fvc::interpolate(rho1)
+ - phiF1
+ )
+ );
+
+ volVectorField U2s
+ (
+ HbyA2
+ + fvc::reconstruct
+ (
+ rAlphaAU2f*mSfGradp/fvc::interpolate(rho2)
+ - phiF2
+ )
+ );
+
+ volScalarField D1(rAU1*dragCoeff);
+ volScalarField D2(rAU2*dragCoeff);
+
+ U = alpha1*(U1s + D1*U2) + alpha2*(U2s + D2*U1);
+ volVectorField Ur(((1 - D2)*U1s - (1 - D1)*U2s)/(1 - D1*D2));
+
+ U1 = U + alpha2*Ur;
+ U1.correctBoundaryConditions();
+ fvOptions.correct(U1);
+
+ U2 = U - alpha1*Ur;
+ U2.correctBoundaryConditions();
+ fvOptions.correct(U2);
+
+ U = fluid.U();
+ }
+ }
+ }
+
+ // Update densities from change in p
+ rho1 += psi1*(p_rgh - p_rgh_0);
+ rho2 += psi2*(p_rgh - p_rgh_0);
+
+ // Update the phase kinetic energies
+ K1 = 0.5*magSqr(U1);
+ K2 = 0.5*magSqr(U2);
+
+ // Update the pressure time-derivative if required
+ if (thermo1.dpdt() || thermo2.dpdt())
+ {
+ dpdt = fvc::ddt(p);
+ }
+}
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 85eb40546558c1544c2c8b54da05e22280fbb02c..766536952a9a31d5bec1e399e1d9d2f05cfb45c2 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -90,16 +90,10 @@ int main(int argc, char *argv[])
- (fvOptions(alpha2, rho2)&rho2)
);
-
#include "UEqns.H"
#include "EEqns.H"
-
- // --- Pressure corrector loop
- while (pimple.correct())
- {
- #include "pEqn.H"
- }
-
+ //#include "pEqn.H"
+ #include "pEqnElim.H"
#include "DDtU.H"
if (pimple.turbCorr())
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
index 1d35530ae0d34292c529318603372e35d472d3a1..a9197ce5a5c221bb884b89325da8cf14f0bf49d2 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -218,6 +218,17 @@ void Foam::MRFZoneList::makeRelative(surfaceScalarField& phi) const
}
+Foam::tmp<Foam::surfaceScalarField> Foam::MRFZoneList::relative
+(
+ const tmp<surfaceScalarField>& phi
+) const
+{
+ tmp<surfaceScalarField> rphi(phi.ptr());
+ makeRelative(rphi());
+ return rphi;
+}
+
+
Foam::tmp<Foam::FieldField<Foam::fvsPatchField, Foam::scalar> >
Foam::MRFZoneList::relative
(
@@ -266,6 +277,17 @@ void Foam::MRFZoneList::makeAbsolute(surfaceScalarField& phi) const
}
+Foam::tmp<Foam::surfaceScalarField> Foam::MRFZoneList::absolute
+(
+ const tmp<surfaceScalarField>& phi
+) const
+{
+ tmp<surfaceScalarField> rphi(phi.ptr());
+ makeAbsolute(rphi());
+ return rphi;
+}
+
+
void Foam::MRFZoneList::makeAbsolute
(
const surfaceScalarField& rho,
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H
index 2cf759c3e94b1353efc086dbac4e5341cb38f2c7..52d992d1ce6a26c2918aaa2c8665aa317db4b7cc 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -118,12 +118,15 @@ public:
//- Make the given absolute velocity relative within the MRF region
void makeRelative(volVectorField& U) const;
- //- Make the given relative velocity absolute within the MRF region
- void makeAbsolute(volVectorField& U) const;
-
//- Make the given absolute flux relative within the MRF region
void makeRelative(surfaceScalarField& phi) const;
+ //- Return the given absolute flux relative within the MRF region
+ tmp<surfaceScalarField> relative
+ (
+ const tmp<surfaceScalarField>& phi
+ ) const;
+
//- Return the given absolute boundary flux relative within
// the MRF region
tmp<FieldField<fvsPatchField, scalar> > relative
@@ -138,9 +141,18 @@ public:
surfaceScalarField& phi
) const;
+ //- Make the given relative velocity absolute within the MRF region
+ void makeAbsolute(volVectorField& U) const;
+
//- Make the given relative flux absolute within the MRF region
void makeAbsolute(surfaceScalarField& phi) const;
+ //- Return the given relative flux absolute within the MRF region
+ tmp<surfaceScalarField> absolute
+ (
+ const tmp<surfaceScalarField>& phi
+ ) const;
+
//- Make the given relative mass-flux absolute within the MRF region
void makeAbsolute
(