diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C index 1e870c42..597b5a49 100644 --- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C @@ -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" diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H index 8696f769..635e5287 100644 --- a/applications/solvers/multiphase/interFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interFoam/alphaEqn.H @@ -2,6 +2,47 @@ word alphaScheme("div(phi,alpha)"); word alpharScheme("div(phirb,alpha)"); + tmp > ddtAlpha + ( + fv::ddtScheme::New + ( + mesh, + mesh.ddtScheme("ddt(alpha)") + ) + ); + + // Set the off-centering coefficient according to ddt scheme + scalar ocCoeff = 0; + if + ( + isType >(ddtAlpha()) + || isType >(ddtAlpha()) + ) + { + ocCoeff = 0; + } + else if (isType >(ddtAlpha())) + { + if (nAlphaSubCycles > 1) + { + FatalErrorIn(args.executable()) + << "Sub-cycling is not supported " + "with the CrankNicolson ddt scheme" + << exit(FatalError); + } + + ocCoeff = + refCast >(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 tphiAlpha; + tmp 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 ( mesh, - phi, - upwind(mesh, phi) - ).fvmDiv(phi, alpha1) + phiCN, + upwind(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 tphiAlphaUD(alpha1Eqn.flux()); - tphiAlpha = tmp - ( - 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 0) + { + tphiAlphaUn = + cnCoeff*tphiAlphaUn + (1.0 - cnCoeff)*phiAlpha.oldTime(); + } + if (MULESCorr) { - tmp tphiAlphaCorr(tphiAlphaUn() - tphiAlpha()); - volScalarField alpha10(alpha1); + tmp 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::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; } diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H index a87673b0..2baf54bf 100644 --- a/applications/solvers/multiphase/interFoam/createFields.H +++ b/applications/solvers/multiphase/interFoam/createFields.H @@ -78,16 +78,6 @@ #include "readGravitationalAcceleration.H" - /* - dimensionedVector g0(g); - - // Read the data file and initialise the interpolation table - interpolationTable 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 tphiAlphaCorr0; diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C index 58d74a62..08f9b5e7 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C @@ -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" diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C index 6b6168bf..5b102360 100644 --- a/applications/solvers/multiphase/interFoam/interFoam.C +++ b/applications/solvers/multiphase/interFoam/interFoam.C @@ -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" diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C index c17847bd..da5a6198 100644 --- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C +++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C @@ -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" diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H index e94c8c2c..1568890d 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H +++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H @@ -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::mesh(); } + //- Return the off-centreing coefficient + scalar ocCoeff() const + { + return ocCoeff_; + } + tmp > fvcDdt ( const dimensioned& diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/dynamicMeshDict b/tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/dynamicMeshDict index e32ae227..68f291d7 100644 --- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/dynamicMeshDict +++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/dynamicMeshDict @@ -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 { diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes index c0a13530..e16dc819 100644 --- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes +++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes @@ -17,7 +17,7 @@ FoamFile ddtSchemes { - default Euler; + default CrankNicolson 0.9; } gradSchemes diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution index 2d0db409..38a7d042 100644 --- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution +++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution @@ -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;