Commit 3df71d18 authored by Henry Weller's avatar Henry Weller Committed by Andrew Heather
Browse files

compressibleInterFoam: Improved mass conservation

using the continuity error correction formulation developed for
twoPhaseEulerFoam and reactingEulerFoam.
parent 092d4f4f
......@@ -113,6 +113,8 @@
phiCN,
upwind<scalar>(mesh, phiCN)
).fvmDiv(phiCN, alpha1)
// - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh)
// + fvc::div(phiCN), alpha1)
==
Su + fvm::Sp(Sp + divU, alpha1)
);
......@@ -125,19 +127,19 @@
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
alphaPhi = talphaPhiUD();
tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
alphaPhi10 = talphaPhi1UD();
if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
{
Info<< "Applying the previous iteration compression flux" << endl;
MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0);
MULES::correct(alpha1, alphaPhi10, talphaPhi1Corr0.ref(), 1, 0);
alphaPhi += talphaPhiCorr0();
alphaPhi10 += talphaPhi1Corr0();
}
// Cache the upwind-flux
talphaPhiCorr0 = talphaPhiUD;
talphaPhi1Corr0 = talphaPhi1UD;
alpha2 = 1.0 - alpha1;
......@@ -151,7 +153,7 @@
surfaceScalarField phir(phic*mixture.nHatf());
tmp<surfaceScalarField> talphaPhiUn
tmp<surfaceScalarField> talphaPhi1Un
(
fvc::flux
(
......@@ -169,15 +171,15 @@
if (MULESCorr)
{
tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
volScalarField alpha10("alpha10", alpha1);
MULES::correct
(
geometricOneField(),
alpha1,
talphaPhiUn(),
talphaPhiCorr.ref(),
talphaPhi1Un(),
talphaPhi1Corr.ref(),
Sp,
(-Sp*alpha1)(),
1,
......@@ -187,24 +189,24 @@
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
alphaPhi += talphaPhiCorr();
alphaPhi10 += talphaPhi1Corr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
alphaPhi += 0.5*talphaPhiCorr();
alphaPhi10 += 0.5*talphaPhi1Corr();
}
}
else
{
alphaPhi = talphaPhiUn;
alphaPhi10 = talphaPhi1Un;
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phiCN,
alphaPhi,
alphaPhi10,
Sp,
(Su + divU*min(alpha1(), scalar(1)))(),
1,
......@@ -219,34 +221,37 @@
if (alphaApplyPrevCorr && MULESCorr)
{
talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
talphaPhiCorr0.ref().rename("alphaPhiCorr0");
talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0;
talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
}
else
{
talphaPhiCorr0.clear();
talphaPhi1Corr0.clear();
}
#include "rhofs.H"
if
(
word(mesh.ddtScheme("ddt(rho,U)"))
== fv::EulerDdtScheme<vector>::typeName
|| word(mesh.ddtScheme("ddt(rho,U)"))
== fv::localEulerDdtScheme<vector>::typeName
)
{
#include "rhofs.H"
rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f;
}
else
{
if (ocCoeff > 0)
{
// Calculate the end-of-time-step alpha flux
alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
alphaPhi10 =
(alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
}
// Calculate the end-of-time-step mass flux
#include "rhofs.H"
rhoPhi = alphaPhi*(rho1f - rho2f) + phi*rho2f;
rhoPhi = alphaPhi10*(rho1f - rho2f) + phi*rho2f;
}
Info<< "Phase-1 volume fraction = "
......
IOobject alphaPhiHeader
IOobject alphaPhi10Header
(
"alphaPhi",
"alphaPhi10",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
);
const bool alphaRestart = alphaPhiHeader.typeHeaderOk<surfaceScalarField>(true);
const bool alphaRestart =
alphaPhi10Header.typeHeaderOk<surfaceScalarField>(true);
// MULES flux from previous time-step
surfaceScalarField alphaPhi
surfaceScalarField alphaPhi10
(
alphaPhiHeader,
alphaPhi10Header,
phi*fvc::interpolate(alpha1)
);
// MULES Correction
tmp<surfaceScalarField> talphaPhiCorr0;
tmp<surfaceScalarField> talphaPhi1Corr0;
{
fvScalarMatrix TEqn
(
fvm::ddt(rho, T)
+ fvm::div(rhoPhi, T)
fvm::ddt(rho, T) + fvm::div(rhoPhi, T)
- fvm::Sp(contErr, T)
- fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
+ (
divU*p
......
fvVectorMatrix UEqn
(
fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
- fvm::Sp(contErr, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
......
......@@ -24,14 +24,14 @@ volScalarField::Internal Su
forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
if (dgdt[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]*alpha1[celli];
Su[celli] += dgdt[celli]*alpha1[celli];
Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
else if (dgdt[celli] < 0.0)
{
Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
}
}
......
tmp<surfaceScalarField> talphaPhi1(alphaPhi10);
if (nAlphaSubCycles > 1)
{
#include "alphaEqnSubCycle.H"
dimensionedScalar totalDeltaT = runTime.deltaT();
talphaPhi1 = new surfaceScalarField
(
IOobject
(
"alphaPhi1",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", alphaPhi10.dimensions(), 0)
);
surfaceScalarField rhoPhiSum
(
IOobject
(
"rhoPhiSum",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", rhoPhi.dimensions(), 0)
);
tmp<volScalarField> trSubDeltaT;
if (LTS)
{
trSubDeltaT =
fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
}
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
#include "alphaEqn.H"
talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi10;
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
#include "alphaEqn.H"
}
rho == alpha1*rho1 + alpha2*rho2;
const surfaceScalarField& alphaPhi1 = talphaPhi1();
surfaceScalarField alphaPhi2("alphaPhi2", phi - alphaPhi1);
volScalarField::Internal contErr
(
(
fvc::ddt(rho) + fvc::div(rhoPhi)
- (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
- (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
)()
);
......@@ -150,8 +150,6 @@ int main(int argc, char *argv[])
#include "alphaControls.H"
#include "compressibleAlphaEqnSubCycle.H"
solve(fvm::ddt(rho) + fvc::div(rhoPhi));
#include "UEqn.H"
#include "TEqn.H"
......
......@@ -56,13 +56,29 @@
}
else
{
#include "rhofs.H"
p_rghEqnComp1 =
fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+ fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
pos(alpha1)
*(
(
fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
- (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
)/rho1
- fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
);
p_rghEqnComp2 =
fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+ fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
pos(alpha2)
*(
(
fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
- (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
)/rho2
- fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh))
);
}
// Cache p_rgh prior to solve for density update
......@@ -78,16 +94,21 @@
solve
(
(
(max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
+ (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
)
+ p_rghEqnIncomp,
p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
mesh.solver(p_rgh.select(pimple.finalInnerIter()))
);
if (pimple.finalNonOrthogonalIter())
{
p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
dgdt =
(
alpha1*(p_rghEqnComp2 & p_rgh)
- alpha2*(p_rghEqnComp1 & p_rgh)
);
phi = phiHbyA + p_rghEqnIncomp.flux();
U = HbyA
......@@ -109,15 +130,9 @@
rho = alpha1*rho1 + alpha2*rho2;
p = max(p_rgh + rho*gh, pMin);
// Correct p_rgh for consistency with p and the updated densities
p_rgh = p - rho*gh;
p_rgh.correctBoundaryConditions();
dgdt =
(
pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
- pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
);
K = 0.5*magSqr(U);
}
......@@ -108,9 +108,7 @@ int main(int argc, char *argv[])
while (pimple.loop())
{
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
solve(fvm::ddt(rho) + fvc::div(rhoPhi));
#include "compressibleAlphaEqnSubCycle.H"
#include "UEqn.H"
volScalarField divU(fvc::div(fvc::absolute(phi, U)));
......
......@@ -89,7 +89,7 @@ surfaceScalarField rhoPhi
volScalarField dgdt
(
pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
pos0(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
);
// Construct compressible turbulence model
......
......@@ -53,13 +53,29 @@
}
else
{
#include "rhofs.H"
p_rghEqnComp1 =
fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+ fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
pos(alpha1)
*(
(
fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
- (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
)/rho1
- fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
);
p_rghEqnComp2 =
fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+ fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
pos(alpha2)
*(
(
fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
- (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
)/rho2
- fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh))
);
}
// Cache p_rgh prior to solve for density update
......@@ -75,11 +91,7 @@
solve
(
(
(max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
+ (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
)
+ p_rghEqnIncomp,
p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
mesh.solver(p_rgh.select(pimple.finalInnerIter()))
);
......@@ -90,8 +102,8 @@
dgdt =
(
pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
- pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
alpha1*(p_rghEqnComp2 & p_rgh)
- alpha2*(p_rghEqnComp1 & p_rgh)
);
phi = phiHbyA + p_rghEqnIncomp.flux();
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -950,7 +950,7 @@ Foam::multiphaseMixtureThermo::nearInterface() const
forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
{
tnearInt.ref() =
max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase()));
max(tnearInt(), pos0(phase() - 0.01)*pos0(0.99 - phase()));
}
return tnearInt;
......
......@@ -104,7 +104,7 @@
)
{
phase().dgdt() =
pos(phase())
pos0(phase())
*(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho();
}
......
......@@ -132,7 +132,7 @@ int main(int argc, char *argv[])
// if the mesh topology changed
if (mesh.topoChanging())
{
talphaPhiCorr0.clear();
talphaPhi1Corr0.clear();
}
gh = (g & mesh.C()) - ghRef;
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -434,7 +434,7 @@ const Foam::surfaceScalarField& Foam::fvMesh::phi() const
// Set zero current time
// mesh motion fluxes if the time has been incremented
if (phiPtr_->timeIndex() != time().timeIndex())
if (!time().subCycling() && phiPtr_->timeIndex() != time().timeIndex())
{
(*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
}
......
......@@ -31,8 +31,6 @@ divSchemes
div(phirb,alpha) Gauss linear;
div(rhoPhi,U) Gauss vanLeerV;
div(phi,thermo:rho.water) Gauss linear;
div(phi,thermo:rho.air) Gauss linear;
div(rhoPhi,T) Gauss linear;
div(rhoPhi,K) Gauss linear;
div((phi+meshPhi),p) Gauss linear;
......
......@@ -31,8 +31,6 @@ divSchemes
div(phirb,alpha) Gauss linear;
div(rhoPhi,U) Gauss upwind;
div(phi,thermo:rho.water) Gauss upwind;
div(phi,thermo:rho.air) Gauss upwind;
div(rhoPhi,T) Gauss upwind;
div(rhoPhi,K) Gauss upwind;
div(phi,p) Gauss upwind;
......
......@@ -31,8 +31,6 @@ divSchemes
div(phirb,alpha) Gauss linear;
div(rhoPhi,U) Gauss upwind;
div(phi,thermo:rho.water) Gauss upwind;
div(phi,thermo:rho.air) Gauss upwind;
div(rhoPhi,T) Gauss upwind;
div(rhoPhi,K) Gauss upwind;
div(phi,p) Gauss upwind;
......
Markdown is supported
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