Browse files

compressibleTwoPhaseEulerFoam: Solve alpha1 equation using MULES and …

…sub-cycling to guarantee boundedness
  • Loading branch information...
1 parent 757a4e7 commit 3e504a3746a33a0e2f65cae4fa207f874d773153 Henry committed Nov 14, 2012
View
4 applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
@@ -1,3 +1,7 @@
+mrfZones.correctBoundaryVelocity(U1);
+mrfZones.correctBoundaryVelocity(U2);
+mrfZones.correctBoundaryVelocity(U);
+
fvVectorMatrix U1Eqn(U1, U1.dimensions()*dimVol/dimTime);
fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
View
91 applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H
@@ -2,8 +2,8 @@ surfaceScalarField alphaPhi1("alphaPhi1", phi1);
surfaceScalarField alphaPhi2("alphaPhi2", phi2);
{
- word scheme("div(phi,alpha1)");
- word schemer("div(phir,alpha1)");
+ word alphaScheme("div(phi,alpha1)");
+ word alpharScheme("div(phir,alpha1)");
surfaceScalarField phic("phic", phi);
surfaceScalarField phir("phir", phi1 - phi2);
@@ -13,7 +13,7 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
surfaceScalarField alpha1f(fvc::interpolate(alpha1));
surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf());
phir += phipp;
- phic += fvc::interpolate(alpha1)*phipp;
+ phic += alpha1f*phipp;
}
for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -56,42 +56,83 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
}
}
+ dimensionedScalar totalDeltaT = runTime.deltaT();
+ if (nAlphaSubCycles > 1)
+ {
+ alphaPhi1 = dimensionedScalar("0", alphaPhi1.dimensions(), 0);
+ }
- fvScalarMatrix alpha1Eqn
+ for
(
- fvm::ddt(alpha1)
- + fvm::div(phic, alpha1, scheme)
- + fvm::div(-fvc::flux(-phir, alpha2, schemer), alpha1, schemer)
- ==
- fvm::Sp(Sp, alpha1) + Su
- );
+ subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+ !(++alphaSubCycle).end();
+ )
+ {
+ surfaceScalarField alphaPhic1
+ (
+ fvc::flux
+ (
+ phic,
+ alpha1,
+ alphaScheme
+ )
+ + fvc::flux
+ (
+ -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+ alpha1,
+ alpharScheme
+ )
+ );
+
+ MULES::explicitSolve
+ (
+ geometricOneField(),
+ alpha1,
+ phi,
+ alphaPhic1,
+ Sp,
+ Su,
+ (g0.value() > 0 ? alphaMax : 1),
+ 0
+ );
+
+ if (nAlphaSubCycles > 1)
+ {
+ alphaPhi1 += (runTime.deltaT()/totalDeltaT)*alphaPhic1;
+ }
+ else
+ {
+ alphaPhi1 = alphaPhic1;
+ }
+ }
if (g0.value() > 0.0)
{
+ surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+
ppMagf =
fvc::interpolate((1.0/rho1)*rAU1)
- *fvc::interpolate
- (
- g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
- );
+ *g0*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
- alpha1Eqn -= fvm::laplacian
+ fvScalarMatrix alpha1Eqn
(
- (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
- alpha1,
- "laplacian(alphaPpMag,alpha1)"
+ fvm::ddt(alpha1) - fvc::ddt(alpha1)
+ - fvm::laplacian
+ (
+ alpha1f*ppMagf,
+ alpha1,
+ "laplacian(alpha1PpMag,alpha1)"
+ )
);
- }
- alpha1Eqn.relax();
- alpha1Eqn.solve();
+ alpha1Eqn.relax();
+ alpha1Eqn.solve();
- //***HGW temporary boundedness-fix pending the introduction of MULES
- alpha1 = max(min(alpha1, scalar(1)), scalar(0));
+ #include "packingLimiter.H"
- #include "packingLimiter.H"
+ alphaPhi1 += alpha1Eqn.flux();
+ }
- alphaPhi1 = alpha1Eqn.flux();
alphaPhi2 = phi - alphaPhi1;
alpha2 = scalar(1) - alpha1;
View
2 ...ications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C
@@ -31,6 +31,8 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
#include "nearWallDist.H"
#include "wallFvPatch.H"
#include "fixedValueFvsPatchFields.H"
View
1 ...ications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
@@ -1,3 +1,4 @@
#include "readTimeControls.H"
int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
+ int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));
View
4 tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSolution
@@ -93,8 +93,8 @@ PIMPLE
nOuterCorrectors 1;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
- nAlphaCorr 2;
- correctAlpha no;
+ nAlphaCorr 1;
+ nAlphaSubCycles 2;
}
relaxationFactors
View
4 tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSolution
@@ -93,8 +93,8 @@ PIMPLE
nOuterCorrectors 1;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
- nAlphaCorr 2;
- correctAlpha yes;
+ nAlphaCorr 1;
+ nAlphaSubCycles 2;
}
relaxationFactors
View
2 tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/fvSolution
@@ -76,7 +76,7 @@ PIMPLE
nCorrectors 3;
nNonOrthogonalCorrectors 0;
nAlphaCorr 1;
- correctAlpha yes;
+ nAlphaSubCycles 2;
pRefCell 0;
pRefValue 0;
}

0 comments on commit 3e504a3

Please sign in to comment.