Skip to content
Snippets Groups Projects
Commit 3bab081c authored by Henry's avatar Henry
Browse files

multiphaseEulerFoam: Remove old file

parent 4351cb49
Branches
Tags
No related merge requests found
surfaceScalarField alphaPhi1("alphaPhi1", phi1);
surfaceScalarField alphaPhi2("alphaPhi2", phi2);
{
word scheme("div(phi,alpha1)");
word schemer("div(phir,alpha1)");
surfaceScalarField phic("phic", phi);
surfaceScalarField phir("phir", phi1 - phi2);
if (g0.value() > 0.0)
{
surfaceScalarField alpha1f = fvc::interpolate(alpha1);
surfaceScalarField phipp = ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
phir += phipp;
phic += fvc::interpolate(alpha1)*phipp;
}
for (int acorr=0; acorr<nAlphaCorr; acorr++)
{
volScalarField::DimensionedInternalField Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
);
volScalarField::DimensionedInternalField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
fvc::div(phi)*min(alpha1, scalar(1))
);
forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]*alpha1[celli];
Su[celli] += dgdt[celli]*alpha1[celli];
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
}
}
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1)
+ fvm::div(phic, alpha1, scheme)
+ fvm::div(-fvc::flux(-phir, alpha2, schemer), alpha1, schemer)
==
fvm::Sp(Sp, alpha1) + Su
);
if (g0.value() > 0.0)
{
ppMagf = rU1Af*fvc::interpolate
(
(1.0/(rho1*(alpha1 + scalar(0.0001))))
*g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
);
alpha1Eqn -= fvm::laplacian
(
(fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
alpha1,
"laplacian(alphaPpMag,alpha1)"
);
}
alpha1Eqn.relax();
alpha1Eqn.solve();
//***HGW temporary boundedness-fix pending the introduction of MULES
alpha1 = max(min(alpha1, 1.0), 0.0);
#include "packingLimiter.H"
alphaPhi1 = alpha1Eqn.flux();
alphaPhi2 = phi - alphaPhi1;
alpha2 = scalar(1) - alpha1;
Info<< "Dispersed phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}
}
rho = alpha1*rho1 + alpha2*rho2;
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment