diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
index e9aa9fb9b9efa54295d25c629fa901ca72bc9073..494acaf05fd83d7cc934cb0dd737d1d90af1b4af 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@@ -77,8 +77,6 @@ int main(int argc, char *argv[])
#include "setrDeltaT.H"
- tmp<surfaceScalarField> tphiAlpha;
-
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
@@ -91,8 +89,6 @@ int main(int argc, char *argv[])
#define LTSSOLVE
#include "alphaEqnSubCycle.H"
#undef LTSSOLVE
-
- interface.correct();
}
turbulence->correct();
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
index 56527781a25b4f3b5093f04ccd5c370438e48635..ccf89e5467dc309e2162679dd8d1445a7fa46329 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
@@ -81,8 +81,6 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
- tmp<surfaceScalarField> tphiAlpha;
-
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index f0e9dab4036dd0cde2e72f377e37290158bd008f..195ed09fc903d15ab4a23b357926364e44fceebe 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -2,11 +2,29 @@
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
- surfaceScalarField phic(mag(phi/mesh.magSf()));
- phic = min(interface.cAlpha()*phic, max(phic));
- surfaceScalarField phir(phic*interface.nHatf());
+ //surfaceScalarField phic(mag(phi/mesh.magSf()));
+ //phic = min(interface.cAlpha()*phic, max(phic));
+ surfaceScalarField phic
+ (
+ (0.5*interface.cAlpha())
+ *fvc::interpolate(volScalarField("tauc", fvc::surfaceSum(mag(phi))))
+ /mesh.magSf()
+ );
+
+ // Do not compress interface at non-coupled boundary faces
+ // (inlets, outlets etc.)
+ forAll(phic.boundaryField(), patchi)
+ {
+ fvsPatchScalarField& phicp = phic.boundaryField()[patchi];
+
+ if (!phicp.coupled())
+ {
+ phicp == 0;
+ }
+ }
+
+ tmp<surfaceScalarField> tphiAlpha;
- //***HGW if (pimple.firstIter() && MULESCorr)
if (MULESCorr)
{
fvScalarMatrix alpha1Eqn
@@ -32,12 +50,37 @@
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
- tphiAlpha = alpha1Eqn.flux();
+ tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
+ tphiAlpha = tmp<surfaceScalarField>
+ (
+ new surfaceScalarField(tphiAlphaUD())
+ );
+
+ if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
+ {
+ Info<< "Applying the previous iteration compression flux" << endl;
+ #ifdef LTSSOLVE
+ MULES::LTScorrect(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
+ #else
+ MULES::correct(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
+ #endif
+
+ tphiAlpha() += tphiAlphaCorr0();
+ }
+
+ // Cache the upwind-flux
+ tphiAlphaCorr0 = tphiAlphaUD;
+
+ alpha2 = 1.0 - alpha1;
+
+ interface.correct();
}
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
- tmp<surfaceScalarField> tphiAlpha0
+ surfaceScalarField phir(phic*interface.nHatf());
+
+ tmp<surfaceScalarField> tphiAlphaUn
(
fvc::flux
(
@@ -47,7 +90,7 @@
)
+ fvc::flux
(
- -fvc::flux(-phir, alpha2, alpharScheme),
+ -fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
@@ -55,17 +98,29 @@
if (MULESCorr)
{
- tphiAlpha0() -= tphiAlpha();
+ tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - tphiAlpha());
+ volScalarField alpha10(alpha1);
+
#ifdef LTSSOLVE
- MULES::LTScorrect(alpha1, tphiAlpha0(), 1, 0);
+ MULES::LTScorrect(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
#else
- MULES::correct(alpha1, tphiAlpha0(), 1, 0);
+ MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
#endif
- tphiAlpha() += tphiAlpha0();
+
+ // Under-relax the correction for more than 3 correctors
+ if (aCorr < 3)
+ {
+ tphiAlpha() += tphiAlphaCorr();
+ }
+ else
+ {
+ alpha1 = 0.5*alpha1 + 0.5*alpha10;
+ tphiAlpha() += 0.5*tphiAlphaCorr();
+ }
}
else
{
- tphiAlpha = tphiAlpha0;
+ tphiAlpha = tphiAlphaUn;
#ifdef LTSSOLVE
MULES::explicitLTSSolve(alpha1, phi, tphiAlpha(), 1, 0);
@@ -75,10 +130,17 @@
}
alpha2 = 1.0 - alpha1;
+
+ interface.correct();
}
rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
+ if (alphaApplyPrevCorr && MULESCorr)
+ {
+ tphiAlphaCorr0 = tphiAlpha() - tphiAlphaCorr0;
+ }
+
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(alpha1) = " << min(alpha1).value()
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 3039171f77a74c38d46bb5770ed85d8357f3d597..00a7d12db6ea3167de00e508756cce472626ffe6 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -135,3 +135,6 @@
}
fv::IOoptionList fvOptions(mesh);
+
+
+ tmp<surfaceScalarField> tphiAlphaCorr0;
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index d56e300581d9d9e4cbd46dc79270392f4d51f887..0a3d7ed11e97ed60f03d44a5b8e61705d303c4f5 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -77,8 +77,6 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
- tmp<surfaceScalarField> tphiAlpha;
-
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 0ddfed50396eefbe37953f16d608fa60b69d83f1..f221ff9ad392adecef1637c88d4cf96a2acb75b4 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -80,8 +80,6 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
- tmp<surfaceScalarField> tphiAlpha;
-
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
index 3321da57f9dffa699ab0a136187f94c34810f49f..fad22244619414c05db570ecc9ebc26b5d81c500 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
@@ -74,8 +74,6 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
- tmp<surfaceScalarField> tphiAlpha;
-
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
index 15b4439220846d0079623443d49ef38b0530ac86..d81f5c93d383300994b3827f11559d6d55935a79 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@@ -83,8 +83,6 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
- tmp<surfaceScalarField> tphiAlpha;
-
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{