Commit f78d33b6 authored by Henry's avatar Henry
Browse files

interFoam family: Add support for MULES-bounded Crank-Nicolson 2nd-order ddt(alpha)

This is an experimental feature demonstrating the potential of MULES to
create bounded solution which are 2nd-order in time AND space.

Crank-Nicolson may be selected on U and/or alpha but will only be fully
2nd-order if used on both within the PIMPLE-loop to converge the
interaction between the flux and phase-fraction.  Note also that
Crank-Nicolson may not be used with sub-cycling but all the features of
semi-implicit MULES are available in particular MULESCorr and
alphaApplyPrevCorr.

Examples of ddt specification:

ddtSchemes
{
    default         Euler;
}

ddtSchemes
{
    default         CrankNicolson 0.9;
}

ddtSchemes
{
    default         none;
    ddt(alpha)      CrankNicolson 0.9;
    ddt(rho,U)      CrankNicolson 0.9;
}

ddtSchemes
{
    default         none;
    ddt(alpha)      Euler;
    ddt(rho,U)      CrankNicolson 0.9;
}

ddtSchemes
{
    default         none;
    ddt(alpha)      CrankNicolson 0.9;
    ddt(rho,U)      Euler;
}

In these examples a small amount of off-centering in used to stabilize
the Crank-Nicolson scheme.  Also the specification for alpha1 is via the
generic phase-fraction name to ensure in multiphase solvers (when
Crank-Nicolson support is added) the scheme is identical for all phase
fractions.
parent 514f9057
......@@ -39,6 +39,8 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulentTransportModel.H"
......
......@@ -2,6 +2,43 @@
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
tmp<fv::ddtScheme<scalar> > ddtAlpha
(
fv::ddtScheme<scalar>::New
(
mesh,
mesh.ddtScheme("ddt(alpha)")
)
);
// Set the off-centering coefficient according to ddt scheme
scalar ocCoeff = 0;
if (isType<fv::EulerDdtScheme<scalar> >(ddtAlpha()))
{
ocCoeff = 0;
}
else if (isType<fv::CrankNicolsonDdtScheme<scalar> >(ddtAlpha()))
{
if (nAlphaSubCycles > 1)
{
FatalErrorIn(args.executable())
<< "Sub-cycling is not supported "
"with the CrankNicolson ddt scheme"
<< exit(FatalError);
}
ocCoeff =
refCast<fv::CrankNicolsonDdtScheme<scalar> >(ddtAlpha()).ocCoeff();
}
else
{
FatalErrorIn(args.executable())
<< "Only Euler and CrankNicolson ddt schemes are supported"
<< exit(FatalError);
}
scalar cnCoeff = ocCoeff/2.0;
// Standard face-flux compression coefficient
surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
......@@ -24,10 +61,16 @@
}
}
tmp<surfaceScalarField> tphiAlpha;
if (MULESCorr)
{
tmp<surfaceScalarField> phiCN(phi);
// Calculate the Crank-Nicolson off-centred volumetric flux
if (ocCoeff > 0)
{
phiCN = (1.0 - cnCoeff)*phi + cnCoeff*phi.oldTime();
}
fvScalarMatrix alpha1Eqn
(
#ifdef LTSSOLVE
......@@ -38,9 +81,9 @@
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi,
upwind<scalar>(mesh, phi)
).fvmDiv(phi, alpha1)
phiCN,
upwind<scalar>(mesh, phiCN)
).fvmDiv(phiCN, alpha1)
);
alpha1Eqn.solve();
......@@ -52,21 +95,18 @@
<< endl;
tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
tphiAlpha = tmp<surfaceScalarField>
(
new surfaceScalarField(tphiAlphaUD())
);
phiAlpha = tphiAlphaUD();
if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
{
Info<< "Applying the previous iteration compression flux" << endl;
#ifdef LTSSOLVE
MULES::LTScorrect(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
MULES::LTScorrect(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
#else
MULES::correct(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
MULES::correct(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
#endif
tphiAlpha() += tphiAlphaCorr0();
phiAlpha += tphiAlphaCorr0();
}
// Cache the upwind-flux
......@@ -77,6 +117,7 @@
mixture.correct();
}
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
surfaceScalarField phir(phic*mixture.nHatf());
......@@ -97,9 +138,17 @@
)
);
// Calculate the Crank-Nicolson off-centred alpha flux
if (ocCoeff > 0)
{
tphiAlphaUn =
(1.0 - cnCoeff)*tphiAlphaUn
+ cnCoeff*phiAlpha.oldTime();
}
if (MULESCorr)
{
tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - tphiAlpha());
tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - phiAlpha);
volScalarField alpha10("alpha10", alpha1);
#ifdef LTSSOLVE
......@@ -111,22 +160,22 @@
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
tphiAlpha() += tphiAlphaCorr();
phiAlpha += tphiAlphaCorr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
tphiAlpha() += 0.5*tphiAlphaCorr();
phiAlpha += 0.5*tphiAlphaCorr();
}
}
else
{
tphiAlpha = tphiAlphaUn;
phiAlpha = tphiAlphaUn;
#ifdef LTSSOLVE
MULES::explicitLTSSolve(alpha1, phi, tphiAlpha(), 1, 0);
MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
#else
MULES::explicitSolve(alpha1, phi, tphiAlpha(), 1, 0);
MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
#endif
}
......@@ -135,11 +184,47 @@
mixture.correct();
}
rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
if (alphaApplyPrevCorr && MULESCorr)
{
tphiAlphaCorr0 = tphiAlpha() - tphiAlphaCorr0;
tphiAlphaCorr0 = phiAlpha - tphiAlphaCorr0;
}
if
(
word(mesh.ddtScheme("ddt(rho,U)"))
== fv::EulerDdtScheme<vector>::typeName
)
{
if (ocCoeff > 0)
{
// Calculate the Crank-Nicolson off-centred volumetric flux
surfaceScalarField phiCN
(
"phiCN",
(1.0 - cnCoeff)*phi + cnCoeff*phi.oldTime()
);
Info<< "Calculate the Crank-Nicolson off-centred mass flux" << endl;
// Calculate the Crank-Nicolson off-centred mass flux
rhoPhi = phiAlpha*(rho1 - rho2) + phiCN*rho2;
}
else
{
// Calculate the end-of-time-step mass flux
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
}
}
else
{
if (ocCoeff > 0)
{
Info<< "Calculate the end-of-time-step alpha flux" << endl;
// Calculate the end-of-time-step alpha flux
phiAlpha = (phiAlpha - cnCoeff*phiAlpha.oldTime())/(1.0 - cnCoeff);
}
// Calculate the end-of-time-step mass flux
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
}
Info<< "Phase-1 volume fraction = "
......
......@@ -78,16 +78,6 @@
#include "readGravitationalAcceleration.H"
/*
dimensionedVector g0(g);
// Read the data file and initialise the interpolation table
interpolationTable<vector> timeSeriesAcceleration
(
runTime.path()/runTime.caseConstant()/"acceleration.dat"
);
*/
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
......@@ -127,6 +117,19 @@
p_rgh = p - rho*gh;
}
// MULES flux from previous time-step
surfaceScalarField phiAlpha
(
IOobject
(
"phiAlpha",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
phi*fvc::interpolate(alpha1)
);
// MULES Correction
tmp<surfaceScalarField> tphiAlphaCorr0;
......@@ -35,6 +35,8 @@ Description
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulentTransportModel.H"
......
......@@ -39,6 +39,8 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulentTransportModel.H"
......
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