164 lines
4.3 KiB
C
164 lines
4.3 KiB
C
{
|
|
word alphaScheme("div(phi,alpha)");
|
|
word alpharScheme("div(phirb,alpha)");
|
|
|
|
surfaceScalarField phir
|
|
(
|
|
IOobject
|
|
(
|
|
"phir",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf()
|
|
);
|
|
|
|
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
|
|
{
|
|
// Create the limiter to be used for all phase-fractions
|
|
scalarField allLambda(mesh.nFaces(), 1.0);
|
|
|
|
// Split the limiter into a surfaceScalarField
|
|
slicedSurfaceScalarField lambda
|
|
(
|
|
IOobject
|
|
(
|
|
"lambda",
|
|
mesh.time().timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
mesh,
|
|
dimless,
|
|
allLambda
|
|
);
|
|
|
|
|
|
// Create the complete convection flux for alpha1
|
|
surfaceScalarField phiAlpha1 =
|
|
fvc::flux
|
|
(
|
|
phi,
|
|
alpha1,
|
|
alphaScheme
|
|
)
|
|
+ fvc::flux
|
|
(
|
|
-fvc::flux(-phir, alpha2, alpharScheme),
|
|
alpha1,
|
|
alpharScheme
|
|
)
|
|
+ fvc::flux
|
|
(
|
|
-fvc::flux(-phir, alpha3, alpharScheme),
|
|
alpha1,
|
|
alpharScheme
|
|
);
|
|
|
|
// Create the bounded (upwind) flux for alpha1
|
|
surfaceScalarField phiAlpha1BD =
|
|
upwind<scalar>(mesh, phi).flux(alpha1);
|
|
|
|
// Calculate the flux correction for alpha1
|
|
phiAlpha1 -= phiAlpha1BD;
|
|
|
|
// Calculate the limiter for alpha1
|
|
MULES::limiter
|
|
(
|
|
allLambda,
|
|
geometricOneField(),
|
|
alpha1,
|
|
phiAlpha1BD,
|
|
phiAlpha1,
|
|
zeroField(),
|
|
zeroField(),
|
|
1,
|
|
0,
|
|
3
|
|
);
|
|
|
|
// Create the complete flux for alpha2
|
|
surfaceScalarField phiAlpha2 =
|
|
fvc::flux
|
|
(
|
|
phi,
|
|
alpha2,
|
|
alphaScheme
|
|
)
|
|
+ fvc::flux
|
|
(
|
|
-fvc::flux(phir, alpha1, alpharScheme),
|
|
alpha2,
|
|
alpharScheme
|
|
);
|
|
|
|
// Create the bounded (upwind) flux for alpha2
|
|
surfaceScalarField phiAlpha2BD =
|
|
upwind<scalar>(mesh, phi).flux(alpha2);
|
|
|
|
// Calculate the flux correction for alpha2
|
|
phiAlpha2 -= phiAlpha2BD;
|
|
|
|
// Further limit the limiter for alpha2
|
|
MULES::limiter
|
|
(
|
|
allLambda,
|
|
geometricOneField(),
|
|
alpha2,
|
|
phiAlpha2BD,
|
|
phiAlpha2,
|
|
zeroField(),
|
|
zeroField(),
|
|
1,
|
|
0,
|
|
3
|
|
);
|
|
|
|
// Construct the limited fluxes
|
|
phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
|
|
phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
|
|
|
|
// Solve for alpha1
|
|
solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
|
|
|
|
// Create the diffusion coefficients for alpha2<->alpha3
|
|
volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2);
|
|
volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3);
|
|
|
|
// Add the diffusive flux for alpha3->alpha2
|
|
phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
|
|
|
|
// Solve for alpha2
|
|
fvScalarMatrix alpha2Eqn
|
|
(
|
|
fvm::ddt(alpha2)
|
|
+ fvc::div(phiAlpha2)
|
|
- fvm::laplacian(Dc23 + Dc32, alpha2)
|
|
);
|
|
alpha2Eqn.solve();
|
|
|
|
// Construct the complete mass flux
|
|
rhoPhi =
|
|
phiAlpha1*(rho1 - rho3)
|
|
+ (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
|
|
+ phi*rho3;
|
|
|
|
alpha3 = 1.0 - alpha1 - alpha2;
|
|
}
|
|
|
|
Info<< "Air phase volume fraction = "
|
|
<< alpha1.weightedAverage(mesh.V()).value()
|
|
<< " Min(alpha1) = " << min(alpha1).value()
|
|
<< " Max(alpha1) = " << max(alpha1).value()
|
|
<< endl;
|
|
|
|
Info<< "Liquid phase volume fraction = "
|
|
<< alpha2.weightedAverage(mesh.V()).value()
|
|
<< " Min(alpha2) = " << min(alpha2).value()
|
|
<< " Max(alpha2) = " << max(alpha2).value()
|
|
<< endl;
|
|
}
|