VoF solvers: Updated PC-MULES

This commit is contained in:
Henry
2013-11-26 21:16:51 +00:00
parent 7a76df3074
commit 7429f59358
8 changed files with 77 additions and 26 deletions

View File

@ -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();

View File

@ -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())
{

View File

@ -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()

View File

@ -135,3 +135,6 @@
}
fv::IOoptionList fvOptions(mesh);
tmp<surfaceScalarField> tphiAlphaCorr0;

View File

@ -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())
{

View File

@ -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())
{

View File

@ -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())
{

View File

@ -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())
{