Commit 8d21b380 authored by Henry Weller's avatar Henry Weller
Browse files

Renamed phiAlpha -P alphaPhi for consistency with Euler-Euler solvers

parent 85b27f67
......@@ -45,7 +45,7 @@
}
surfaceScalarField phiAlpha1
surfaceScalarField alphaPhi1
(
fvc::flux
(
......@@ -66,7 +66,7 @@
geometricOneField(),
alpha1,
phi,
phiAlpha1,
alphaPhi1,
Sp,
Su,
1,
......@@ -75,7 +75,7 @@
surfaceScalarField rho1f(fvc::interpolate(rho1));
surfaceScalarField rho2f(fvc::interpolate(rho2));
rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
rhoPhi = alphaPhi1*(rho1f - rho2f) + phi*rho2f;
alpha2 = scalar(1) - alpha1;
}
......
......@@ -942,14 +942,14 @@ void Foam::multiphaseMixtureThermo::solveAlphas
surfaceScalarField phic(mag(phi_/mesh_.magSf()));
phic = min(cAlpha*phic, max(phic));
PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
int phasei = 0;
forAllIter(PtrDictionary<phaseModel>, phases_, phase)
{
phaseModel& alpha = phase();
phiAlphaCorrs.set
alphaPhiCorrs.set
(
phasei,
new surfaceScalarField
......@@ -964,7 +964,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
)
);
surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
{
......@@ -974,7 +974,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
surfaceScalarField phir(phic*nHatf(alpha, alpha2));
phiAlphaCorr += fvc::flux
alphaPhiCorr += fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha,
......@@ -988,7 +988,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
geometricOneField(),
alpha,
phi_,
phiAlphaCorr,
alphaPhiCorr,
zeroField(),
zeroField(),
1,
......@@ -999,7 +999,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
phasei++;
}
MULES::limitSum(phiAlphaCorrs);
MULES::limitSum(alphaPhiCorrs);
rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
......@@ -1025,8 +1025,8 @@ void Foam::multiphaseMixtureThermo::solveAlphas
{
phaseModel& alpha = phase();
surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha);
surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
volScalarField::DimensionedInternalField Sp
(
......@@ -1096,12 +1096,12 @@ void Foam::multiphaseMixtureThermo::solveAlphas
(
geometricOneField(),
alpha,
phiAlpha,
alphaPhi,
Sp,
Su
);
rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*phiAlpha;
rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*alphaPhi;
Info<< alpha.name() << " volume fraction, min, max = "
<< alpha.weightedAverage(mesh_.V()).value()
......
......@@ -23,32 +23,32 @@
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
phiAlpha = tphiAlphaUD();
tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
alphaPhi = talphaPhiUD();
if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
{
Info<< "Applying the previous iteration correction flux" << endl;
MULES::correct
(
alpha1,
phiAlpha,
tphiAlphaCorr0(),
alphaPhi,
talphaPhiCorr0(),
mixture.alphaMax(),
0
);
phiAlpha += tphiAlphaCorr0();
alphaPhi += talphaPhiCorr0();
}
// Cache the upwind-flux
tphiAlphaCorr0 = tphiAlphaUD;
talphaPhiCorr0 = talphaPhiUD;
}
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
tmp<surfaceScalarField> tphiAlphaUn
tmp<surfaceScalarField> talphaPhiUn
(
fvc::flux
(
......@@ -66,14 +66,14 @@
if (MULESCorr)
{
tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - phiAlpha);
tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
volScalarField alpha10("alpha10", alpha1);
MULES::correct
(
alpha1,
tphiAlphaUn(),
tphiAlphaCorr(),
talphaPhiUn(),
talphaPhiCorr(),
mixture.alphaMax(),
0
);
......@@ -81,23 +81,23 @@
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
phiAlpha += tphiAlphaCorr();
alphaPhi += talphaPhiCorr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
phiAlpha += 0.5*tphiAlphaCorr();
alphaPhi += 0.5*talphaPhiCorr();
}
}
else
{
phiAlpha = tphiAlphaUn;
alphaPhi = talphaPhiUn;
MULES::explicitSolve
(
alpha1,
phi,
phiAlpha,
alphaPhi,
mixture.alphaMax(),
0
);
......@@ -106,7 +106,7 @@
if (alphaApplyPrevCorr && MULESCorr)
{
tphiAlphaCorr0 = phiAlpha - tphiAlphaCorr0;
talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
}
alpha2 = 1.0 - alpha1;
......
{
surfaceScalarField phiAlpha
surfaceScalarField alphaPhi
(
IOobject
(
"phiAlpha",
"alphaPhi",
runTime.timeName(),
mesh
),
......@@ -19,11 +19,11 @@
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField phiAlphaSum
surfaceScalarField alphaPhiSum
(
IOobject
(
"phiAlphaSum",
"alphaPhiSum",
runTime.timeName(),
mesh
),
......@@ -38,10 +38,10 @@
)
{
#include "alphaEqn.H"
phiAlphaSum += (runTime.deltaT()/totalDeltaT)*phiAlpha;
alphaPhiSum += (runTime.deltaT()/totalDeltaT)*alphaPhi;
}
phiAlpha = phiAlphaSum;
alphaPhi = alphaPhiSum;
}
else
{
......@@ -59,7 +59,7 @@
alpha1Eqn.solve(mesh.solver("alpha1Diffusion"));
phiAlpha += alpha1Eqn.flux();
alphaPhi += alpha1Eqn.flux();
alpha2 = 1.0 - alpha1;
Info<< "Phase-1 volume fraction = "
......@@ -69,6 +69,6 @@
<< endl;
}
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
rho = mixture.rho();
}
......@@ -135,4 +135,4 @@ mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
// MULES Correction
tmp<surfaceScalarField> tphiAlphaCorr0;
tmp<surfaceScalarField> talphaPhiCorr0;
......@@ -98,19 +98,19 @@
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
phiAlpha = tphiAlphaUD();
tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
alphaPhi = talphaPhiUD();
if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
{
Info<< "Applying the previous iteration compression flux" << endl;
MULES::correct(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
MULES::correct(alpha1, alphaPhi, talphaPhiCorr0(), 1, 0);
phiAlpha += tphiAlphaCorr0();
alphaPhi += talphaPhiCorr0();
}
// Cache the upwind-flux
tphiAlphaCorr0 = tphiAlphaUD;
talphaPhiCorr0 = talphaPhiUD;
alpha2 = 1.0 - alpha1;
......@@ -122,7 +122,7 @@
{
surfaceScalarField phir(phic*mixture.nHatf());
tmp<surfaceScalarField> tphiAlphaUn
tmp<surfaceScalarField> talphaPhiUn
(
fvc::flux
(
......@@ -141,33 +141,33 @@
// Calculate the Crank-Nicolson off-centred alpha flux
if (ocCoeff > 0)
{
tphiAlphaUn =
cnCoeff*tphiAlphaUn + (1.0 - cnCoeff)*phiAlpha.oldTime();
talphaPhiUn =
cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
}
if (MULESCorr)
{
tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - phiAlpha);
tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
volScalarField alpha10("alpha10", alpha1);
MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr(), 1, 0);
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
phiAlpha += tphiAlphaCorr();
alphaPhi += talphaPhiCorr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
phiAlpha += 0.5*tphiAlphaCorr();
alphaPhi += 0.5*talphaPhiCorr();
}
}
else
{
phiAlpha = tphiAlphaUn;
alphaPhi = talphaPhiUn;
MULES::explicitSolve(alpha1, phiCN, phiAlpha, 1, 0);
MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
}
alpha2 = 1.0 - alpha1;
......@@ -177,7 +177,7 @@
if (alphaApplyPrevCorr && MULESCorr)
{
tphiAlphaCorr0 = phiAlpha - tphiAlphaCorr0;
talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
}
if
......@@ -186,18 +186,18 @@
== fv::EulerDdtScheme<vector>::typeName
)
{
rhoPhi = phiAlpha*(rho1 - rho2) + phiCN*rho2;
rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
}
else
{
if (ocCoeff > 0)
{
// Calculate the end-of-time-step alpha flux
phiAlpha = (phiAlpha - (1.0 - cnCoeff)*phiAlpha.oldTime())/cnCoeff;
alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
}
// Calculate the end-of-time-step mass flux
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
}
Info<< "Phase-1 volume fraction = "
......
......@@ -121,11 +121,11 @@ mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
// MULES flux from previous time-step
surfaceScalarField phiAlpha
surfaceScalarField alphaPhi
(
IOobject
(
"phiAlpha",
"alphaPhi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
......@@ -135,4 +135,4 @@ surfaceScalarField phiAlpha
);
// MULES Correction
tmp<surfaceScalarField> tphiAlphaCorr0;
tmp<surfaceScalarField> talphaPhiCorr0;
......@@ -40,7 +40,7 @@
// Create the complete convection flux for alpha1
surfaceScalarField phiAlpha1
surfaceScalarField alphaPhi1
(
fvc::flux
(
......@@ -63,13 +63,13 @@
);
// Create the bounded (upwind) flux for alpha1
surfaceScalarField phiAlpha1BD
surfaceScalarField alphaPhi1BD
(
upwind<scalar>(mesh, phi).flux(alpha1)
);
// Calculate the flux correction for alpha1
phiAlpha1 -= phiAlpha1BD;
alphaPhi1 -= alphaPhi1BD;
// Calculate the limiter for alpha1
if (LTS)
......@@ -83,8 +83,8 @@
rDeltaT,
geometricOneField(),
alpha1,
phiAlpha1BD,
phiAlpha1,
alphaPhi1BD,
alphaPhi1,
zeroField(),
zeroField(),
1,
......@@ -99,8 +99,8 @@
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha1,
phiAlpha1BD,
phiAlpha1,
alphaPhi1BD,
alphaPhi1,
zeroField(),
zeroField(),
1,
......@@ -109,7 +109,7 @@
}
// Create the complete flux for alpha2
surfaceScalarField phiAlpha2
surfaceScalarField alphaPhi2
(
fvc::flux
(
......@@ -126,13 +126,13 @@
);
// Create the bounded (upwind) flux for alpha2
surfaceScalarField phiAlpha2BD
surfaceScalarField alphaPhi2BD
(
upwind<scalar>(mesh, phi).flux(alpha2)
);
// Calculate the flux correction for alpha2
phiAlpha2 -= phiAlpha2BD;
alphaPhi2 -= alphaPhi2BD;
// Further limit the limiter for alpha2
if (LTS)
......@@ -146,8 +146,8 @@
rDeltaT,
geometricOneField(),
alpha2,
phiAlpha2BD,
phiAlpha2,
alphaPhi2BD,
alphaPhi2,
zeroField(),
zeroField(),
1,
......@@ -162,8 +162,8 @@
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha2,
phiAlpha2BD,
phiAlpha2,
alphaPhi2BD,
alphaPhi2,
zeroField(),
zeroField(),
1,
......@@ -172,32 +172,32 @@
}
// Construct the limited fluxes
phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
alphaPhi1 = alphaPhi1BD + lambda*alphaPhi1;
alphaPhi2 = alphaPhi2BD + lambda*alphaPhi2;
// Solve for alpha1
solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
solve(fvm::ddt(alpha1) + fvc::div(alphaPhi1));
// 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);
alphaPhi2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
// Solve for alpha2
fvScalarMatrix alpha2Eqn
(
fvm::ddt(alpha2)
+ fvc::div(phiAlpha2)
+ fvc::div(alphaPhi2)
- fvm::laplacian(Dc23 + Dc32, alpha2)
);
alpha2Eqn.solve();
// Construct the complete mass flux
rhoPhi =
phiAlpha1*(rho1 - rho3)
+ (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
alphaPhi1*(rho1 - rho3)
+ (alphaPhi2 + alpha2Eqn.flux())*(rho2 - rho3)
+ phi*rho3;
alpha3 = 1.0 - alpha1 - alpha2;
......
......@@ -10,7 +10,7 @@
const volScalarField& vDotvAlphal = vDotAlphal[1]();
const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
tmp<surfaceScalarField> tphiAlpha;
tmp<surfaceScalarField> talphaPhi;
if (MULESCorr)
{
......@@ -37,14 +37,14 @@
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
tphiAlpha = alpha1Eqn.flux();
talphaPhi = alpha1Eqn.flux();
}
volScalarField alpha10("alpha10", alpha1);
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
tmp<surfaceScalarField> tphiAlphaCorr
tmp<surfaceScalarField> talphaPhiCorr
(
fvc::flux
(
......@@ -62,7 +62,7 @@
if (MULESCorr)
{
tphiAlphaCorr() -= tphiAlpha();
talphaPhiCorr() -= talphaPhi();
volScalarField alpha100("alpha100", alpha10);
alpha10 = alpha1;
......@@ -71,8 +71,8 @@
(
geometricOneField(),
alpha1,
tphiAlpha(),
tphiAlphaCorr(),
talphaPhi(),
talphaPhiCorr(),
vDotvmcAlphal,
(
divU*(alpha10 - alpha100)
......@@ -85,12 +85,12 @@
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
tphiAlpha() += tphiAlphaCorr();
talphaPhi() += talphaPhiCorr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
tphiAlpha() += 0.5*tphiAlphaCorr();
talphaPhi() += 0.5*talphaPhiCorr();
}
}
else
......@@ -100,20 +100,20 @@
geometricOneField(),
alpha1,
phi,
tphiAlphaCorr(),
talphaPhiCorr(),
vDotvmcAlphal,
(divU*alpha1 + vDotcAlphal)(),
1,
0
);
tphiAlpha = tphiAlphaCorr;
talphaPhi = talphaPhiCorr;
}
alpha2 = 1.0 - alpha1;
}
rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;