Given that the type of the dimensioned internal field is encapsulated in the GeometricField class the name need not include "Field"; the type name is "Internal" so volScalarField::DimensionedInternalField -> volScalarField::Internal In addition to the ".dimensionedInternalField()" access function the simpler "()" de-reference operator is also provided to greatly simplify FV equation source term expressions which need not evaluate boundary conditions. To demonstrate this kEpsilon.C has been updated to use dimensioned internal field expressions in the k and epsilon equation source terms.
88 lines
2.2 KiB
C
88 lines
2.2 KiB
C
{
|
|
word alphaScheme("div(phi,alpha)");
|
|
word alpharScheme("div(phirb,alpha)");
|
|
|
|
surfaceScalarField phir(phic*interface.nHatf());
|
|
|
|
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
|
|
{
|
|
volScalarField::Internal Sp
|
|
(
|
|
IOobject
|
|
(
|
|
"Sp",
|
|
runTime.timeName(),
|
|
mesh
|
|
),
|
|
mesh,
|
|
dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
|
|
);
|
|
|
|
volScalarField::Internal Su
|
|
(
|
|
IOobject
|
|
(
|
|
"Su",
|
|
runTime.timeName(),
|
|
mesh
|
|
),
|
|
// Divergence term is handled explicitly to be
|
|
// consistent with the explicit transport solution
|
|
divU*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]);
|
|
}
|
|
}
|
|
|
|
|
|
surfaceScalarField alphaPhi1
|
|
(
|
|
fvc::flux
|
|
(
|
|
phi,
|
|
alpha1,
|
|
alphaScheme
|
|
)
|
|
+ fvc::flux
|
|
(
|
|
-fvc::flux(-phir, alpha2, alpharScheme),
|
|
alpha1,
|
|
alpharScheme
|
|
)
|
|
);
|
|
|
|
MULES::explicitSolve
|
|
(
|
|
geometricOneField(),
|
|
alpha1,
|
|
phi,
|
|
alphaPhi1,
|
|
Sp,
|
|
Su,
|
|
1,
|
|
0
|
|
);
|
|
|
|
surfaceScalarField rho1f(fvc::interpolate(rho1));
|
|
surfaceScalarField rho2f(fvc::interpolate(rho2));
|
|
rhoPhi = alphaPhi1*(rho1f - rho2f) + phi*rho2f;
|
|
|
|
alpha2 = scalar(1) - alpha1;
|
|
}
|
|
|
|
Info<< "Liquid phase volume fraction = "
|
|
<< alpha1.weightedAverage(mesh.V()).value()
|
|
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
|
|
<< " Min(" << alpha2.name() << ") = " << min(alpha2).value()
|
|
<< endl;
|
|
}
|