research/OpenFOAM-2.4.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

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

Browse source 
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.

This development is back-ported from OpenFOAM-dev:
f78d33b634
This commit is contained in:
Henry 2015-02-25 16:32:06 +00:00
parent 59603afc8f
commit 1cacdf0d89
10 changed files with 188 additions and 60 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

9
applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -39,6 +39,9 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"
@ -54,11 +57,11 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "initContinuityErrs.H"
#include "createFields.H"
pimpleControl pimple(mesh);
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createPrghCorrTypes.H"
#include "correctPhi.H"
#include "CourantNo.H"

120
applications/solvers/multiphase/interFoam/alphaEqn.H
View file

@ -2,6 +2,47 @@
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())
|| isType<fv::localEulerDdtScheme<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 = 1.0/(1.0 + ocCoeff);
// Standard face-flux compression coefficient
surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
@ -24,7 +65,13 @@
}
}
tmp<surfaceScalarField> tphiAlpha;
tmp<surfaceScalarField> phiCN(phi);
// Calculate the Crank-Nicolson off-centred volumetric flux
if (ocCoeff > 0)
{
phiCN = cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime();
}
if (MULESCorr)
{
@ -38,35 +85,32 @@
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi,
upwind<scalar>(mesh, phi)
).fvmDiv(phi, alpha1)
phiCN,
upwind<scalar>(mesh, phiCN)
).fvmDiv(phiCN, alpha1)
);
alpha1Eqn.solve();
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() <<") = " << min(alpha1).value()
<< " Max(" << alpha1.name() <<") = " << max(alpha1).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< 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 +121,7 @@
mixture.correct();
}
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
surfaceScalarField phir(phic*mixture.nHatf());
@ -97,10 +142,17 @@
)
);
// Calculate the Crank-Nicolson off-centred alpha flux
if (ocCoeff > 0)
{
tphiAlphaUn =
cnCoeff*tphiAlphaUn + (1.0 - cnCoeff)*phiAlpha.oldTime();
}
if (MULESCorr)
{
tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - tphiAlpha());
volScalarField alpha10(alpha1);
tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - phiAlpha);
volScalarField alpha10("alpha10", alpha1);
#ifdef LTSSOLVE
MULES::LTScorrect(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
@ -111,22 +163,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, phiCN, phiAlpha, 1, 0);
#endif
}
@ -135,16 +187,34 @@
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
)
{
rhoPhi = phiAlpha*(rho1 - rho2) + phiCN*rho2;
}
else
{
if (ocCoeff > 0)
{
// Calculate the end-of-time-step alpha flux
phiAlpha = (phiAlpha - (1.0 - cnCoeff)*phiAlpha.oldTime())/cnCoeff;
}
// Calculate the end-of-time-step mass flux
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
}
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() <<") = " << min(alpha1).value()
<< " Max(" << alpha1.name() <<") = " << max(alpha1).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
}

26
applications/solvers/multiphase/interFoam/createFields.H
View file

@ -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());
@ -111,7 +101,7 @@
(
p,
p_rgh,
mesh.solutionDict().subDict("PIMPLE"),
pimple.dict(),
pRefCell,
pRefValue
);
@ -131,5 +121,19 @@
fv::IOoptionList fvOptions(mesh);
// 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;

7
applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,6 +35,9 @@ Description
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"
@ -49,10 +52,10 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
#include "initContinuityErrs.H"
#include "createFields.H"
#include "readTimeControls.H"
#include "createPrghCorrTypes.H"

5
applications/solvers/multiphase/interFoam/interFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -39,6 +39,9 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"

11
applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -40,6 +40,9 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"
@ -55,13 +58,13 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
pimpleControl pimple(mesh);
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createPorousZones.H"
#include "readTimeControls.H"
pimpleControl pimple(mesh);
#include "createPrghCorrTypes.H"
#include "correctPhi.H"
#include "CourantNo.H"

51
src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,9 +25,48 @@ Class
Foam::fv::CrankNicolsonDdtScheme
Description
Second-oder CrankNicolson implicit ddt using the current and
Second-oder Crank-Nicolson implicit ddt using the current and
previous time-step fields as well as the previous time-step ddt.
The Crank-Nicolson scheme is often unstable for complex flows in complex
geometries and it is necessary to "off-centre" the scheme to stabilize it
while retaining greater temporal accuracy than the first-order
Euler-implicit scheme. Off-centering is specified via the mandatory
coefficient in the range [0,1] following the scheme name e.g.
\verbatim
ddtSchemes
{
default CrankNicolson 0.9;
}
\endverbatim
With a coefficient of 1 the scheme is fully centred and second-order,
with a coefficient of 0 the scheme is equivalent to Euler-implicit.
A coefficient of 0.9 has been found to be suitable for a range of cases for
which higher-order accuracy in time is needed and provides similar accuracy
and stability to the backward scheme. However, the backward scheme has
been found to be more robust and provides formal second-order accuracy in
time.
The advantage of the Crank-Nicolson scheme over backward is that only the
new and old-time values are needed, the additional terms relating to the
fluxes and sources are evaluated at the mid-point of the time-step which
provides the opportunity to limit the fluxes in such a way as to ensure
boundedness while maintaining greater accuracy in time compared to the
Euler-implicit scheme. This approach is now used with MULES in the
interFoam family of solvers. Boundedness cannot be guaranteed with the
backward scheme.
Note
The Crank-Nicolson coefficient for the implicit part of the RHS is related
to the off-centering coefficient by
\verbatim
cnCoeff = 1.0/(1.0 + ocCoeff);
\endverbatim
See Also
Foam::Euler
Foam::backward
SourceFiles
CrankNicolsonDdtScheme.C
@ -173,7 +212,7 @@ public:
(
"CrankNicolsonDdtScheme(const fvMesh& mesh, Istream& is)",
is
) << "coefficient = " << ocCoeff_
) << "Off-centreing coefficient = " << ocCoeff_
<< " should be >= 0 and <= 1"
<< exit(FatalIOError);
}
@ -188,6 +227,12 @@ public:
return fv::ddtScheme<Type>::mesh();
}
//- Return the off-centreing coefficient
scalar ocCoeff() const
{
return ocCoeff_;
}
tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDdt
(
const dimensioned<Type>&

4
tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/dynamicMeshDict
View file

@ -1,7 +1,7 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.x |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
@ -59,7 +59,7 @@ sixDoFRigidBodyMotionCoeffs
};
report on;
accelerationRelaxation 0.5;
accelerationRelaxation 0.3;
constraints
{

2
tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
View file

@ -17,7 +17,7 @@ FoamFile
ddtSchemes
{
default Euler;
default CrankNicolson 0.9;
}
gradSchemes

13
tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
View file

@ -19,7 +19,7 @@ solvers
{
"alpha.water.*"
{
nAlphaCorr 1;
nAlphaCorr 2;
nAlphaSubCycles 1;
cAlpha 1;
@ -116,19 +116,16 @@ solvers
PIMPLE
{
momentumPredictor no;
nOuterCorrectors 5;
nCorrectors 1;
momentumPredictor no;
nOuterCorrectors 5;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
correctPhi yes;
correctPhi yes;
moveMeshOuterCorrectors yes;
}
relaxationFactors
{
fields
{
}
equations
{
".*" 1;