diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C index 980d1be630..317a18120c 100644 --- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C @@ -38,7 +38,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" @@ -79,7 +79,10 @@ int main(int argc, char *argv[]) twoPhaseProperties.correct(); + #define LTSSOLVE #include "alphaEqnSubCycle.H" + #undef LTSSOLVE + interface.correct(); turbulence->correct(); diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H deleted file mode 100644 index faae197670..0000000000 --- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H +++ /dev/null @@ -1,39 +0,0 @@ -{ - word alphaScheme("div(phi,alpha)"); - word alpharScheme("div(phirb,alpha)"); - - surfaceScalarField phic(mag(phi/mesh.magSf())); - phic = min(interface.cAlpha()*phic, max(phic)); - surfaceScalarField phir(phic*interface.nHatf()); - - for (int aCorr=0; aCorr 1) -{ - dimensionedScalar totalDeltaT = runTime.deltaT(); - surfaceScalarField rhoPhiSum - ( - IOobject - ( - "rhoPhiSum", - runTime.timeName(), - mesh - ), - mesh, - dimensionedScalar("0", rhoPhi.dimensions(), 0) - ); - - for - ( - subCycle alphaSubCycle(alpha1, nAlphaSubCycles); - !(++alphaSubCycle).end(); - ) - { - #include "alphaEqn.H" - rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; - } - - rhoPhi = rhoPhiSum; -} -else -{ - #include "alphaEqn.H" -} - -rho == alpha1*rho1 + alpha2*rho2; diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C index 887ef9e34b..0fa1a50ee2 100644 --- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C +++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C @@ -37,7 +37,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H index a2720e20ee..b0dd8ebef2 100644 --- a/applications/solvers/multiphase/interFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interFoam/alphaEqn.H @@ -6,9 +6,39 @@ phic = min(interface.cAlpha()*phic, max(phic)); surfaceScalarField phir(phic*interface.nHatf()); + tmp tphiAlpha; + + if (MULESCorr) + { + fvScalarMatrix alpha1Eqn + ( + #ifdef LTSSOLVE + fv::localEulerDdtScheme(mesh, rDeltaT.name()).fvmDdt(alpha1) + #else + fv::EulerDdtScheme(mesh).fvmDdt(alpha1) + #endif + + fv::gaussConvectionScheme + ( + mesh, + phi, + upwind(mesh, phi) + ).fvmDiv(phi, alpha1) + ); + + alpha1Eqn.solve(); + + Info<< "Phase-1 volume fraction = " + << alpha1.weightedAverage(mesh.Vsc()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() + << endl; + + tphiAlpha = alpha1Eqn.flux(); + } + for (int aCorr=0; aCorr tphiAlpha0 ( fvc::flux ( @@ -24,12 +54,32 @@ ) ); - MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0); + if (MULESCorr) + { + tphiAlpha0() -= tphiAlpha(); + #ifdef LTSSOLVE + MULES::LTScorrect(alpha1, tphiAlpha0(), 1, 0); + #else + MULES::correct(alpha1, tphiAlpha0(), 1, 0); + #endif + tphiAlpha() += tphiAlpha0(); + } + else + { + tphiAlpha = tphiAlpha0; + + #ifdef LTSSOLVE + MULES::explicitLTSSolve(alpha1, phi, tphiAlpha(), 1, 0); + #else + MULES::explicitSolve(alpha1, phi, tphiAlpha(), 1, 0); + #endif + } alpha2 = 1.0 - alpha1; - rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2; } + rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2; + Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(alpha1) = " << min(alpha1).value() diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C index 032f2db730..4a673c6a3f 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C @@ -34,7 +34,7 @@ Description #include "fvCFD.H" #include "dynamicFvMesh.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C index fa9b4fffe4..4ce08d27db 100644 --- a/applications/solvers/multiphase/interFoam/interFoam.C +++ b/applications/solvers/multiphase/interFoam/interFoam.C @@ -38,7 +38,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C index 93cb8a0c3d..12bc57d785 100644 --- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C +++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C @@ -31,7 +31,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "threePhaseMixture.H" #include "threePhaseInterfaceProperties.H" diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C index 64938174dc..3423a4654a 100644 --- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C +++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C @@ -39,7 +39,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "incompressibleTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H index b6e70fc295..45cba01d5b 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H @@ -16,8 +16,14 @@ { fvScalarMatrix alpha1Eqn ( - fvm::ddt(alpha1) - + fvm::div(phi, alpha1, "UD") - fvm::Sp(divU, alpha1) + fv::EulerDdtScheme(mesh).fvmDdt(alpha1) + + fv::gaussConvectionScheme + ( + mesh, + phi, + upwind(mesh, phi) + ).fvmDiv(phi, alpha1) + - fvm::Sp(divU, alpha1) == fvm::Sp(vDotvmcAlphal, alpha1) + vDotcAlphal diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C index 630b94534e..31fa619754 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C @@ -43,7 +43,7 @@ Description #include "fvCFD.H" #include "dynamicFvMesh.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "phaseChangeTwoPhaseMixture.H" diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C index c1d13ca926..28ef7f0a69 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C @@ -41,7 +41,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "MULES.H" +#include "CMULES.H" #include "subCycle.H" #include "interfaceProperties.H" #include "phaseChangeTwoPhaseMixture.H" diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index a99d3e3a31..212a22baf5 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -220,6 +220,7 @@ fields/surfaceFields/surfaceFields.C fvMatrices/fvMatrices.C fvMatrices/fvScalarMatrix/fvScalarMatrix.C fvMatrices/solvers/MULES/MULES.C +fvMatrices/solvers/MULES/CMULES.C fvMatrices/solvers/MULES/IMULES.C fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C new file mode 100644 index 0000000000..a954dccf56 --- /dev/null +++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C @@ -0,0 +1,68 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "CMULES.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void Foam::MULES::correct +( + volScalarField& psi, + surfaceScalarField& phiPsiCorr, + const scalar psiMax, + const scalar psiMin +) +{ + correct + ( + geometricOneField(), + psi, + phiPsiCorr, + zeroField(), zeroField(), + psiMax, psiMin + ); +} + + +void Foam::MULES::LTScorrect +( + volScalarField& psi, + surfaceScalarField& phiPsiCorr, + const scalar psiMax, + const scalar psiMin +) +{ + LTScorrect + ( + geometricOneField(), + psi, + phiPsiCorr, + zeroField(), zeroField(), + psiMax, psiMin + ); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H new file mode 100644 index 0000000000..2786a86aae --- /dev/null +++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H @@ -0,0 +1,160 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Global + MULES + +Description + CMULES: Multidimensional universal limiter for + explicit corrected implicit solution. + + Solve a convective-only transport equation using an explicit universal + multi-dimensional limiter to correct an implicit conservative bounded + obtained using rigorously bounded schemes such as Euler-implicit in time + upwind in space. + + Parameters are the variable to solve, the normal convective flux and the + actual explicit flux of the variable which is also used to return limited + flux used in the bounded-solution. + +SourceFiles + CMULES.C + CMULESTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef CMULES_H +#define CMULES_H + +#include "MULES.H" +#include "EulerDdtScheme.H" +#include "localEulerDdtScheme.H" +#include "gaussConvectionScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace MULES +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +void correct +( + const RdeltaTType& rDeltaT, + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su +); + +template +void correct +( + const RhoType& rho, + volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +); + +void correct +( + volScalarField& psi, + surfaceScalarField& phiCorr, + const scalar psiMax, + const scalar psiMin +); + +template +void LTScorrect +( + const RhoType& rho, + volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +); + +void LTScorrect +( + volScalarField& psi, + surfaceScalarField& phiCorr, + const scalar psiMax, + const scalar psiMin +); + + +template +void limiterCorr +( + scalarField& allLambda, + const RdeltaTType& rDeltaT, + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +); + +template +void limitCorr +( + const RdeltaTType& rDeltaT, + const RhoType& rho, + const volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace MULES +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "CMULESTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C new file mode 100644 index 0000000000..2ff99745da --- /dev/null +++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C @@ -0,0 +1,493 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "CMULES.H" +#include "fvcSurfaceIntegrate.H" +#include "slicedSurfaceFields.H" +#include "wedgeFvPatch.H" +#include "syncTools.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +void Foam::MULES::correct +( + const RdeltaTType& rDeltaT, + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su +) +{ + Info<< "MULES: Correcting " << psi.name() << endl; + + const fvMesh& mesh = psi.mesh(); + + scalarField psiIf(psi.size(), 0); + fvc::surfaceIntegrate(psiIf, phiCorr); + + if (mesh.moving()) + { + psi.internalField() = + ( + rho.field()*psi.internalField()*rDeltaT + + Su.field() + - psiIf + )/(rho.field()*rDeltaT - Sp.field()); + } + else + { + psi.internalField() = + ( + rho.field()*psi.internalField()*rDeltaT + + Su.field() + - psiIf + )/(rho.field()*rDeltaT - Sp.field()); + } + + psi.correctBoundaryConditions(); +} + + +template +void Foam::MULES::correct +( + const RhoType& rho, + volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +) +{ + const fvMesh& mesh = psi.mesh(); + const scalar rDeltaT = 1.0/mesh.time().deltaTValue(); + + const dictionary& MULEScontrols = mesh.solverDict(psi.name()); + + label nLimiterIter + ( + readLabel(MULEScontrols.lookup("nLimiterIter")) + ); + + limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter); + correct(rDeltaT, rho, psi, phiCorr, Sp, Su); +} + + +template +void Foam::MULES::LTScorrect +( + const RhoType& rho, + volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +) +{ + const fvMesh& mesh = psi.mesh(); + + const volScalarField& rDeltaT = + mesh.objectRegistry::lookupObject("rSubDeltaT"); + + const dictionary& MULEScontrols = mesh.solverDict(psi.name()); + + label nLimiterIter + ( + readLabel(MULEScontrols.lookup("nLimiterIter")) + ); + + limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter); + correct(rDeltaT, rho, psi, phiCorr, Sp, Su); +} + + +template +void Foam::MULES::limiterCorr +( + scalarField& allLambda, + const RdeltaTType& rDeltaT, + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +) +{ + const scalarField& psiIf = psi; + const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField(); + + const fvMesh& mesh = psi.mesh(); + + const labelUList& owner = mesh.owner(); + const labelUList& neighb = mesh.neighbour(); + tmp tVsc = mesh.Vsc(); + const scalarField& V = tVsc(); + + const scalarField& phiCorrIf = phiCorr; + const surfaceScalarField::GeometricBoundaryField& phiCorrBf = + phiCorr.boundaryField(); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + scalarField& lambdaIf = lambda; + surfaceScalarField::GeometricBoundaryField& lambdaBf = + lambda.boundaryField(); + + scalarField psiMaxn(psiIf.size(), psiMin); + scalarField psiMinn(psiIf.size(), psiMax); + + scalarField sumPhip(psiIf.size(), VSMALL); + scalarField mSumPhim(psiIf.size(), VSMALL); + + forAll(phiCorrIf, facei) + { + label own = owner[facei]; + label nei = neighb[facei]; + + psiMaxn[own] = max(psiMaxn[own], psiIf[nei]); + psiMinn[own] = min(psiMinn[own], psiIf[nei]); + + psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]); + psiMinn[nei] = min(psiMinn[nei], psiIf[own]); + + scalar phiCorrf = phiCorrIf[facei]; + + if (phiCorrf > 0.0) + { + sumPhip[own] += phiCorrf; + mSumPhim[nei] += phiCorrf; + } + else + { + mSumPhim[own] -= phiCorrf; + sumPhip[nei] -= phiCorrf; + } + } + + forAll(phiCorrBf, patchi) + { + const fvPatchScalarField& psiPf = psiBf[patchi]; + const scalarField& phiCorrPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + if (psiPf.coupled()) + { + const scalarField psiPNf(psiPf.patchNeighbourField()); + + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]); + } + } + else + { + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]); + } + } + + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + scalar phiCorrf = phiCorrPf[pFacei]; + + if (phiCorrf > 0.0) + { + sumPhip[pfCelli] += phiCorrf; + } + else + { + mSumPhim[pfCelli] -= phiCorrf; + } + } + } + + psiMaxn = min(psiMaxn, psiMax); + psiMinn = max(psiMinn, psiMin); + + // scalar smooth = 0.5; + // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax); + // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin); + + psiMaxn = + V + *( + (rho.field()*rDeltaT - Sp.field())*psiMaxn + - Su.field() + - rho.field()*psi.internalField()*rDeltaT + ); + + psiMinn = + V + *( + Su.field() + - (rho.field()*rDeltaT - Sp.field())*psiMinn + + rho.field()*psi.internalField()*rDeltaT + ); + + scalarField sumlPhip(psiIf.size()); + scalarField mSumlPhim(psiIf.size()); + + for (int j=0; j 0.0) + { + sumlPhip[own] += lambdaPhiCorrf; + mSumlPhim[nei] += lambdaPhiCorrf; + } + else + { + mSumlPhim[own] -= lambdaPhiCorrf; + sumlPhip[nei] -= lambdaPhiCorrf; + } + } + + forAll(lambdaBf, patchi) + { + scalarField& lambdaPf = lambdaBf[patchi]; + const scalarField& phiCorrfPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + forAll(lambdaPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei]; + + if (lambdaPhiCorrf > 0.0) + { + sumlPhip[pfCelli] += lambdaPhiCorrf; + } + else + { + mSumlPhim[pfCelli] -= lambdaPhiCorrf; + } + } + } + + forAll(sumlPhip, celli) + { + sumlPhip[celli] = + max(min + ( + (sumlPhip[celli] + psiMaxn[celli]) + /(mSumPhim[celli] - SMALL), + 1.0), 0.0 + ); + + mSumlPhim[celli] = + max(min + ( + (mSumlPhim[celli] + psiMinn[celli]) + /(sumPhip[celli] + SMALL), + 1.0), 0.0 + ); + } + + const scalarField& lambdam = sumlPhip; + const scalarField& lambdap = mSumlPhim; + + forAll(lambdaIf, facei) + { + if (phiCorrIf[facei] > 0.0) + { + lambdaIf[facei] = min + ( + lambdaIf[facei], + min(lambdap[owner[facei]], lambdam[neighb[facei]]) + ); + } + else + { + lambdaIf[facei] = min + ( + lambdaIf[facei], + min(lambdam[owner[facei]], lambdap[neighb[facei]]) + ); + } + } + + + forAll(lambdaBf, patchi) + { + fvsPatchScalarField& lambdaPf = lambdaBf[patchi]; + const scalarField& phiCorrfPf = phiCorrBf[patchi]; + const fvPatchScalarField& psiPf = psiBf[patchi]; + + if (isA(mesh.boundary()[patchi])) + { + lambdaPf = 0; + } + else if (psiPf.coupled()) + { + const labelList& pFaceCells = + mesh.boundary()[patchi].faceCells(); + + forAll(lambdaPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + if (phiCorrfPf[pFacei] > 0.0) + { + lambdaPf[pFacei] = + min(lambdaPf[pFacei], lambdap[pfCelli]); + } + else + { + lambdaPf[pFacei] = + min(lambdaPf[pFacei], lambdam[pfCelli]); + } + } + } + else + { + const labelList& pFaceCells = + mesh.boundary()[patchi].faceCells(); + const scalarField& phiCorrPf = phiCorrBf[patchi]; + + forAll(lambdaPf, pFacei) + { + // Limit outlet faces only + if (phiCorrPf[pFacei] > SMALL*SMALL) + { + label pfCelli = pFaceCells[pFacei]; + + if (phiCorrfPf[pFacei] > 0.0) + { + lambdaPf[pFacei] = + min(lambdaPf[pFacei], lambdap[pfCelli]); + } + else + { + lambdaPf[pFacei] = + min(lambdaPf[pFacei], lambdam[pfCelli]); + } + } + } + } + } + + syncTools::syncFaceList(mesh, allLambda, minEqOp()); + } +} + + +template +void Foam::MULES::limitCorr +( + const RdeltaTType& rDeltaT, + const RhoType& rho, + const volScalarField& psi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +) +{ + const fvMesh& mesh = psi.mesh(); + + scalarField allLambda(mesh.nFaces(), 1.0); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + limiterCorr + ( + allLambda, + rDeltaT, + rho, + psi, + phiCorr, + Sp, + Su, + psiMax, + psiMin, + nLimiterIter + ); + + phiCorr *= lambda; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H index 1866e8022f..2d6e4a26a0 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H +++ b/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H @@ -25,7 +25,14 @@ Global IMULES Description - IMULES: Multidimensional universal limiter with implicit solution. + IMULES: Multidimensional universal limiter for implicit solution. + + Solve a convective-only transport equation using an explicit universal + multi-dimensional limiter applied to an implicit formulation requiring + iteration to guarantee boundedness. The number of iterations required + to obtain boundedness increases with the Courant number of the simulation. + + It may be more efficient to use CMULES. SourceFiles IMULES.C diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C index f2898b82b8..cb0970b5c4 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C @@ -69,25 +69,6 @@ void Foam::MULES::explicitLTSSolve } -void Foam::MULES::correct -( - volScalarField& psi, - surfaceScalarField& phiPsiCorr, - const scalar psiMax, - const scalar psiMin -) -{ - correct - ( - geometricOneField(), - psi, - phiPsiCorr, - zeroField(), zeroField(), - psiMax, psiMin - ); -} - - void Foam::MULES::limitSum(UPtrList& phiPsiCorrs) { forAll(phiPsiCorrs[0], facei) diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H index b4253f6ee8..f0525ddc9e 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H @@ -25,7 +25,7 @@ Global MULES Description - MULES: Multidimensional universal limiter with explicit solution. + MULES: Multidimensional universal limiter for explicit solution. Solve a convective-only transport equation using an explicit universal multi-dimensional limiter. @@ -126,38 +126,6 @@ void explicitLTSSolve ); -template -void correct -( - const RdeltaTType& rDeltaT, - const RhoType& rho, - volScalarField& psi, - const surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su -); - -template -void correct -( - const RhoType& rho, - volScalarField& psi, - surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin -); - -void correct -( - volScalarField& psi, - surfaceScalarField& phiCorr, - const scalar psiMax, - const scalar psiMin -); - - template void limiter ( @@ -191,41 +159,12 @@ void limit ); -template -void limiterCorr -( - scalarField& allLambda, - const RdeltaTType& rDeltaT, - const RhoType& rho, - const volScalarField& psi, - const surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin, - const label nLimiterIter -); - -template -void limitCorr -( - const RdeltaTType& rDeltaT, - const RhoType& rho, - const volScalarField& psi, - surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin, - const label nLimiterIter -); - - void limitSum(UPtrList& phiPsiCorrs); template void limitSum(SurfaceScalarFieldList& phiPsiCorrs); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace MULES diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C index a0d2e94daa..a7a5282818 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C +++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C @@ -138,74 +138,6 @@ void Foam::MULES::explicitLTSSolve } -template -void Foam::MULES::correct -( - const RdeltaTType& rDeltaT, - const RhoType& rho, - volScalarField& psi, - const surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su -) -{ - Info<< "MULES: Correcting " << psi.name() << endl; - - const fvMesh& mesh = psi.mesh(); - - scalarField psiIf(psi.size(), 0); - fvc::surfaceIntegrate(psiIf, phiCorr); - - if (mesh.moving()) - { - psi.internalField() = - ( - rho.field()*psi.internalField()*rDeltaT - + Su.field() - - psiIf - )/(rho.field()*rDeltaT - Sp.field()); - } - else - { - psi.internalField() = - ( - rho.field()*psi.internalField()*rDeltaT - + Su.field() - - psiIf - )/(rho.field()*rDeltaT - Sp.field()); - } - - psi.correctBoundaryConditions(); -} - - -template -void Foam::MULES::correct -( - const RhoType& rho, - volScalarField& psi, - surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin -) -{ - const fvMesh& mesh = psi.mesh(); - const scalar rDeltaT = 1.0/mesh.time().deltaTValue(); - - const dictionary& MULEScontrols = mesh.solverDict(psi.name()); - - label nLimiterIter - ( - readLabel(MULEScontrols.lookup("nLimiterIter")) - ); - - limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter); - correct(rDeltaT, rho, psi, phiCorr, Sp, Su); -} - - template void Foam::MULES::limiter ( @@ -626,368 +558,6 @@ void Foam::MULES::limit } -template -void Foam::MULES::limiterCorr -( - scalarField& allLambda, - const RdeltaTType& rDeltaT, - const RhoType& rho, - const volScalarField& psi, - const surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin, - const label nLimiterIter -) -{ - const scalarField& psiIf = psi; - const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField(); - - const fvMesh& mesh = psi.mesh(); - - const labelUList& owner = mesh.owner(); - const labelUList& neighb = mesh.neighbour(); - tmp tVsc = mesh.Vsc(); - const scalarField& V = tVsc(); - - const scalarField& phiCorrIf = phiCorr; - const surfaceScalarField::GeometricBoundaryField& phiCorrBf = - phiCorr.boundaryField(); - - slicedSurfaceScalarField lambda - ( - IOobject - ( - "lambda", - mesh.time().timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh, - dimless, - allLambda, - false // Use slices for the couples - ); - - scalarField& lambdaIf = lambda; - surfaceScalarField::GeometricBoundaryField& lambdaBf = - lambda.boundaryField(); - - scalarField psiMaxn(psiIf.size(), psiMin); - scalarField psiMinn(psiIf.size(), psiMax); - - scalarField sumPhip(psiIf.size(), VSMALL); - scalarField mSumPhim(psiIf.size(), VSMALL); - - forAll(phiCorrIf, facei) - { - label own = owner[facei]; - label nei = neighb[facei]; - - psiMaxn[own] = max(psiMaxn[own], psiIf[nei]); - psiMinn[own] = min(psiMinn[own], psiIf[nei]); - - psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]); - psiMinn[nei] = min(psiMinn[nei], psiIf[own]); - - scalar phiCorrf = phiCorrIf[facei]; - - if (phiCorrf > 0.0) - { - sumPhip[own] += phiCorrf; - mSumPhim[nei] += phiCorrf; - } - else - { - mSumPhim[own] -= phiCorrf; - sumPhip[nei] -= phiCorrf; - } - } - - forAll(phiCorrBf, patchi) - { - const fvPatchScalarField& psiPf = psiBf[patchi]; - const scalarField& phiCorrPf = phiCorrBf[patchi]; - - const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); - - if (psiPf.coupled()) - { - const scalarField psiPNf(psiPf.patchNeighbourField()); - - forAll(phiCorrPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]); - psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]); - } - } - else - { - forAll(phiCorrPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]); - psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]); - } - } - - forAll(phiCorrPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - scalar phiCorrf = phiCorrPf[pFacei]; - - if (phiCorrf > 0.0) - { - sumPhip[pfCelli] += phiCorrf; - } - else - { - mSumPhim[pfCelli] -= phiCorrf; - } - } - } - - psiMaxn = min(psiMaxn, psiMax); - psiMinn = max(psiMinn, psiMin); - - // scalar smooth = 0.5; - // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax); - // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin); - - psiMaxn = - V - *( - (rho.field()*rDeltaT - Sp.field())*psiMaxn - - Su.field() - - rho.field()*psi.internalField()*rDeltaT - ); - - psiMinn = - V - *( - Su.field() - - (rho.field()*rDeltaT - Sp.field())*psiMinn - + rho.field()*psi.internalField()*rDeltaT - ); - - scalarField sumlPhip(psiIf.size()); - scalarField mSumlPhim(psiIf.size()); - - for (int j=0; j 0.0) - { - sumlPhip[own] += lambdaPhiCorrf; - mSumlPhim[nei] += lambdaPhiCorrf; - } - else - { - mSumlPhim[own] -= lambdaPhiCorrf; - sumlPhip[nei] -= lambdaPhiCorrf; - } - } - - forAll(lambdaBf, patchi) - { - scalarField& lambdaPf = lambdaBf[patchi]; - const scalarField& phiCorrfPf = phiCorrBf[patchi]; - - const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); - - forAll(lambdaPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei]; - - if (lambdaPhiCorrf > 0.0) - { - sumlPhip[pfCelli] += lambdaPhiCorrf; - } - else - { - mSumlPhim[pfCelli] -= lambdaPhiCorrf; - } - } - } - - forAll(sumlPhip, celli) - { - sumlPhip[celli] = - max(min - ( - (sumlPhip[celli] + psiMaxn[celli]) - /(mSumPhim[celli] - SMALL), - 1.0), 0.0 - ); - - mSumlPhim[celli] = - max(min - ( - (mSumlPhim[celli] + psiMinn[celli]) - /(sumPhip[celli] + SMALL), - 1.0), 0.0 - ); - } - - const scalarField& lambdam = sumlPhip; - const scalarField& lambdap = mSumlPhim; - - forAll(lambdaIf, facei) - { - if (phiCorrIf[facei] > 0.0) - { - lambdaIf[facei] = min - ( - lambdaIf[facei], - min(lambdap[owner[facei]], lambdam[neighb[facei]]) - ); - } - else - { - lambdaIf[facei] = min - ( - lambdaIf[facei], - min(lambdam[owner[facei]], lambdap[neighb[facei]]) - ); - } - } - - - forAll(lambdaBf, patchi) - { - fvsPatchScalarField& lambdaPf = lambdaBf[patchi]; - const scalarField& phiCorrfPf = phiCorrBf[patchi]; - const fvPatchScalarField& psiPf = psiBf[patchi]; - - if (isA(mesh.boundary()[patchi])) - { - lambdaPf = 0; - } - else if (psiPf.coupled()) - { - const labelList& pFaceCells = - mesh.boundary()[patchi].faceCells(); - - forAll(lambdaPf, pFacei) - { - label pfCelli = pFaceCells[pFacei]; - - if (phiCorrfPf[pFacei] > 0.0) - { - lambdaPf[pFacei] = - min(lambdaPf[pFacei], lambdap[pfCelli]); - } - else - { - lambdaPf[pFacei] = - min(lambdaPf[pFacei], lambdam[pfCelli]); - } - } - } - else - { - const labelList& pFaceCells = - mesh.boundary()[patchi].faceCells(); - const scalarField& phiCorrPf = phiCorrBf[patchi]; - - forAll(lambdaPf, pFacei) - { - // Limit outlet faces only - if (phiCorrPf[pFacei] > SMALL*SMALL) - { - label pfCelli = pFaceCells[pFacei]; - - if (phiCorrfPf[pFacei] > 0.0) - { - lambdaPf[pFacei] = - min(lambdaPf[pFacei], lambdap[pfCelli]); - } - else - { - lambdaPf[pFacei] = - min(lambdaPf[pFacei], lambdam[pfCelli]); - } - } - } - } - } - - syncTools::syncFaceList(mesh, allLambda, minEqOp()); - } -} - - -template -void Foam::MULES::limitCorr -( - const RdeltaTType& rDeltaT, - const RhoType& rho, - const volScalarField& psi, - surfaceScalarField& phiCorr, - const SpType& Sp, - const SuType& Su, - const scalar psiMax, - const scalar psiMin, - const label nLimiterIter -) -{ - const fvMesh& mesh = psi.mesh(); - - scalarField allLambda(mesh.nFaces(), 1.0); - - slicedSurfaceScalarField lambda - ( - IOobject - ( - "lambda", - mesh.time().timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh, - dimless, - allLambda, - false // Use slices for the couples - ); - - limiterCorr - ( - allLambda, - rDeltaT, - rho, - psi, - phiCorr, - Sp, - Su, - psiMax, - psiMin, - nLimiterIter - ); - - phiCorr *= lambda; -} - - template void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs) { diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict b/tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict index be315b40ac..05b109778c 100644 --- a/tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict +++ b/tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict @@ -17,7 +17,7 @@ FoamFile application interFoam; -startFrom latestTime; +startFrom startTime; startTime 0; @@ -25,7 +25,7 @@ stopAt endTime; endTime 100; -deltaT 0.01; +deltaT 0.1; writeControl adjustableRunTime; @@ -46,8 +46,9 @@ timePrecision 6; runTimeModifiable yes; adjustTimeStep yes; -maxCo 0.5; -maxAlphaCo 0.5; + +maxCo 6; +maxAlphaCo 6; maxDeltaT 1; functions diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes index 40a4e9640d..11fc30f272 100644 --- a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes +++ b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes @@ -27,12 +27,13 @@ gradSchemes divSchemes { - div(rho*phi,U) Gauss linearUpwind grad(U); - div(phi,alpha) Gauss vanLeer; - div(phirb,alpha) Gauss interfaceCompression; + default none; - div(phi,k) Gauss upwind; - div(phi,omega) $div(phi,k); + div(rho*phi,U) Gauss linearUpwind grad(U); + div(phi,alpha) Gauss vanLeer; + div(phirb,alpha) Gauss interfaceCompression; + + "div\(phi,(k|omega)\)" Gauss upwind; div((muEff*dev(T(grad(U))))) Gauss linear; } diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution index d2428017c4..f6208ef3d2 100644 --- a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution +++ b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution @@ -20,8 +20,16 @@ solvers alpha.water { nAlphaCorr 1; - nAlphaSubCycles 3; + nAlphaSubCycles 1; cAlpha 1; + + MULESCorr yes; + nLimiterIter 3; + + solver PBiCG; + preconditioner DILU; + tolerance 1e-8; + relTol 0; } pcorr @@ -30,9 +38,9 @@ solvers preconditioner { preconditioner GAMG; - tolerance 1e-05; + tolerance 1e-5; relTol 0; - smoother DICGaussSeidel; + smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; @@ -41,31 +49,43 @@ solvers agglomerator faceAreaPair; mergeLevels 1; } - tolerance 1e-10; + tolerance 1e-5; relTol 0; - maxIter 100; + maxIter 50; } p_rgh { - $pcorr; - tolerance 1e-6; - relTol 0.01; + solver GAMG; + tolerance 5e-9; + relTol 0.01; + + smoother GaussSeidel; + nPreSweeps 0; + nPostSweeps 2; + + cacheAgglomeration true; + + nCellsInCoarsestLevel 10; + agglomerator faceAreaPair; + mergeLevels 1; + + maxIter 50; }; p_rghFinal { $p_rgh; - tolerance 1e-6; + tolerance 5e-9; relTol 0; } "(U|k|omega).*" { solver smoothSolver; - smoother GaussSeidel; + smoother symGaussSeidel; nSweeps 1; - tolerance 1e-7; + tolerance 1e-6; relTol 0.1; }; } @@ -84,6 +104,7 @@ relaxationFactors } equations { + ".*" 1; } } diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes index ded24b94b6..bc5571af9e 100644 --- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes +++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes @@ -32,7 +32,6 @@ divSchemes div(phi,k) Gauss linearUpwind grad(k); div(phi,alpha) Gauss vanLeer; div(phirb,alpha) Gauss interfaceCompression; - UD Gauss upwind; div((muEff*dev(T(grad(U))))) Gauss linear; }