Commit b167c95f authored by Henry Weller's avatar Henry Weller
Browse files

compressibleInterFoam: Completed LTS and semi-implicit MULES support

Now the interFoam and compressibleInterFoam families of solvers use the same
alphaEqn formulation and supporting all of the MULES options without
code-duplication.

The semi-implicit MULES support allows running with significantly larger
time-steps but this does reduce the interface sharpness.
parent 6fc2a3dc
......@@ -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
......
......@@ -79,6 +79,8 @@
if (MULESCorr)
{
#include "alphaSuSp.H"
fvScalarMatrix alpha1Eqn
(
(
......@@ -92,6 +94,8 @@
phiCN,
upwind<scalar>(mesh, phiCN)
).fvmDiv(phiCN, alpha1)
==
Su + fvm::Sp(Sp + divU, alpha1)
);
alpha1Eqn.solve();
......@@ -124,6 +128,8 @@
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
#include "alphaSuSp.H"
surfaceScalarField phir(phic*mixture.nHatf());
tmp<surfaceScalarField> talphaPhiUn
......@@ -154,7 +160,17 @@
tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
volScalarField alpha10("alpha10", alpha1);
MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr.ref(), 1, 0);
MULES::correct
(
geometricOneField(),
alpha1,
talphaPhiUn(),
talphaPhiCorr.ref(),
Sp,
(-Sp*alpha1)(),
1,
0
);
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
......@@ -171,7 +187,17 @@
{
alphaPhi = talphaPhiUn;
MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phiCN,
alphaPhi,
Sp,
(Su + divU*min(alpha1(), scalar(1)))(),
1,
0
);
}
alpha2 = 1.0 - alpha1;
......@@ -195,7 +221,8 @@
== fv::EulerDdtScheme<vector>::typeName
)
{
rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
#include "rhofs.H"
rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
}
else
{
......@@ -206,7 +233,8 @@
}
// Calculate the end-of-time-step mass flux
rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
#include "rhofs.H"
rhoPhi = alphaPhi*(rho1f - rho2f) + phi*rho2f;
}
Info<< "Phase-1 volume fraction = "
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
EXE_INC = \
-I. \
-I../interFoam \
-I../VoF \
-ItwoPhaseMixtureThermo \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
......
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
fvOptions(rho, U)
......
const dictionary& alphaControls = mesh.solverDict(alpha1.name());
label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir(phic*mixture.nHatf());
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
#include "alphaSuSp.H"
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;
}
{
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(mixture.cAlpha()*phic, max(phic));
volScalarField divU(fvc::div(fvc::absolute(phi, U)));
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum
(
IOobject
(
"rhoPhiSum",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", rhoPhi.dimensions(), 0)
);
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
#include "alphaEqn.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
#include "alphaEqn.H"
}
}
volScalarField::Internal Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
);
volScalarField::Internal Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Sp", dgdt.dimensions(), 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))
);
volScalarField::Internal Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Su", dgdt.dimensions(), 0)
);
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]);
}
}
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]);
}
}
volScalarField::Internal divU
(
mesh.moving()
? fvc::div(phiCN() + mesh.phi())
: fvc::div(phiCN())
);
EXE_INC = \
-I. \
-I.. \
-I../../interFoam \
-I../../VoF \
-I../twoPhaseMixtureThermo \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
......
......@@ -39,7 +39,10 @@ Description
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "MULES.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "twoPhaseMixture.H"
#include "twoPhaseMixtureThermo.H"
......@@ -47,7 +50,6 @@ Description
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -141,7 +143,7 @@ int main(int argc, char *argv[])
}
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
#include "compressibleAlphaEqnSubCycle.H"
solve(fvm::ddt(rho) + fvc::div(rhoPhi));
......
......@@ -8,6 +8,7 @@
fvc::flux(HbyA)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
);
MRF.makeRelative(phiHbyA);
surfaceScalarField phig
(
......@@ -20,7 +21,7 @@
phiHbyA += phig;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, U, phiHbyA, rAUf);
constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phiHbyA, U);
......
......@@ -36,7 +36,10 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "rhoThermo.H"
#include "twoPhaseMixture.H"
......@@ -44,7 +47,6 @@ Description
#include "turbulentFluidThermoModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -63,6 +63,7 @@ dimensionedScalar pMin
);
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
#include "readGravitationalAcceleration.H"
......@@ -99,3 +100,22 @@ autoPtr<compressible::turbulenceModel> turbulence
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
// MULES flux from previous time-step
surfaceScalarField alphaPhi
(
IOobject
(
"alphaPhi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
phi*fvc::interpolate(alpha1)
);
// MULES Correction
tmp<surfaceScalarField> talphaPhiCorr0;
#include "createMRF.H"
......@@ -8,6 +8,7 @@
fvc::flux(HbyA)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
MRF.makeRelative(phiHbyA);
surfaceScalarField phig
(
......@@ -20,7 +21,7 @@
phiHbyA += phig;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, U, phiHbyA, rAUf);
constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
tmp<fvScalarMatrix> p_rghEqnComp1;
tmp<fvScalarMatrix> p_rghEqnComp2;
......
surfaceScalarField rho1f(fvc::interpolate(rho1));
surfaceScalarField rho2f(fvc::interpolate(rho2));
EXE_INC = \
-I. \
-I../VoF \
-I../interFoam \
-ImultiphaseMixtureThermo/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
......
Supports Markdown
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