zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Multi-phase solvers: Improved handling of inflow/outflow BCs in MULES

Browse Source 
Avoids slight phase-fraction unboundedness at entertainment BCs and improved
robustness.

Additionally the phase-fractions in the multi-phase (rather than two-phase)
solvers are adjusted to avoid the slow growth of inconsistency ("drift") caused
by solving for all of the phase-fractions rather than deriving one from the
others.
This commit is contained in:
Henry Weller
2017-01-17 22:43:47 +00:00
parent 8b930836d3
commit ad92287afc
24 changed files with 214 additions and 606 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/finiteVolume/Make/files
View File

@ -235,7 +235,6 @@ 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
interpolation = interpolation/interpolation

50
src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -151,17 +151,17 @@ void Foam::MULES::limiterCorr
const dictionary& MULEScontrols = mesh.solverDict(psi.name());
label nLimiterIter
const label nLimiterIter
(
readLabel(MULEScontrols.lookup("nLimiterIter"))
);
scalar smoothLimiter
const scalar smoothLimiter
(
MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
);
scalar extremaCoeff
const scalar extremaCoeff
(
MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0)
);
@ -207,8 +207,8 @@ void Foam::MULES::limiterCorr
forAll(phiCorrIf, facei)
{
label own = owner[facei];
label nei = neighb[facei];
const label own = owner[facei];
const label nei = neighb[facei];
psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
psiMinn[own] = min(psiMinn[own], psiIf[nei]);
@ -216,9 +216,9 @@ void Foam::MULES::limiterCorr
psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
scalar phiCorrf = phiCorrIf[facei];
const scalar phiCorrf = phiCorrIf[facei];
if (phiCorrf > 0.0)
if (phiCorrf > 0)
{
sumPhip[own] += phiCorrf;
mSumPhim[nei] += phiCorrf;
@ -253,7 +253,7 @@ void Foam::MULES::limiterCorr
{
forAll(phiCorrPf, pFacei)
{
label pfCelli = pFaceCells[pFacei];
const label pfCelli = pFaceCells[pFacei];
psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
@ -262,11 +262,11 @@ void Foam::MULES::limiterCorr
forAll(phiCorrPf, pFacei)
{
label pfCelli = pFaceCells[pFacei];
const label pfCelli = pFaceCells[pFacei];
scalar phiCorrf = phiCorrPf[pFacei];
const scalar phiCorrf = phiCorrPf[pFacei];
if (phiCorrf > 0.0)
if (phiCorrf > 0)
{
sumPhip[pfCelli] += phiCorrf;
}
@ -309,17 +309,17 @@ void Foam::MULES::limiterCorr
for (int j=0; j<nLimiterIter; j++)
{
sumlPhip = 0.0;
mSumlPhim = 0.0;
sumlPhip = 0;
mSumlPhim = 0;
forAll(lambdaIf, facei)
{
label own = owner[facei];
label nei = neighb[facei];
const label own = owner[facei];
const label nei = neighb[facei];
scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
const scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
if (lambdaPhiCorrf > 0.0)
if (lambdaPhiCorrf > 0)
{
sumlPhip[own] += lambdaPhiCorrf;
mSumlPhim[nei] += lambdaPhiCorrf;
@ -344,7 +344,7 @@ void Foam::MULES::limiterCorr
scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
if (lambdaPhiCorrf > 0.0)
if (lambdaPhiCorrf > 0)
{
sumlPhip[pfCelli] += lambdaPhiCorrf;
}
@ -379,7 +379,7 @@ void Foam::MULES::limiterCorr
forAll(lambdaIf, facei)
{
if (phiCorrIf[facei] > 0.0)
if (phiCorrIf[facei] > 0)
{
lambdaIf[facei] = min
(
@ -415,9 +415,9 @@ void Foam::MULES::limiterCorr
forAll(lambdaPf, pFacei)
{
label pfCelli = pFaceCells[pFacei];
const label pfCelli = pFaceCells[pFacei];
if (phiCorrfPf[pFacei] > 0.0)
if (phiCorrfPf[pFacei] > 0)
{
lambdaPf[pFacei] =
min(lambdaPf[pFacei], lambdap[pfCelli]);
@ -438,11 +438,11 @@ void Foam::MULES::limiterCorr
forAll(lambdaPf, pFacei)
{
// Limit outlet faces only
if (phiPf[pFacei] > SMALL*SMALL)
if ((phiPf[pFacei] + phiCorrfPf[pFacei]) > SMALL*SMALL)
{
label pfCelli = pFaceCells[pFacei];
const label pfCelli = pFaceCells[pFacei];
if (phiCorrfPf[pFacei] > 0.0)
if (phiCorrfPf[pFacei] > 0)
{
lambdaPf[pFacei] =
min(lambdaPf[pFacei], lambdap[pfCelli]);

51
src/finiteVolume/fvMatrices/solvers/MULES/IMULES.C
View File

@ -1,51 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "IMULES.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::MULES::implicitSolve
(
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
const scalar psiMax,
const scalar psiMin
)
{
implicitSolve
(
geometricOneField(),
psi,
phi,
phiPsi,
zeroField(), zeroField(),
psiMax, psiMin
);
}
// ************************************************************************* //

94
src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H
View File

@ -1,94 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 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 <http://www.gnu.org/licenses/>.
Global
IMULES
Description
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
IMULESTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef IMULES_H
#define IMULES_H
#include "MULES.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace MULES
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class RhoType, class SpType, class SuType>
void implicitSolve
(
const RhoType& rho,
volScalarField& gamma,
const surfaceScalarField& phi,
surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
const scalar psiMin
);
void implicitSolve
(
volScalarField& gamma,
const surfaceScalarField& phi,
surfaceScalarField& phiCorr,
const scalar psiMax,
const scalar psiMin
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace MULES
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "IMULESTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

236
src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C
View File

@ -1,236 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "IMULES.H"
#include "gaussConvectionScheme.H"
#include "surfaceInterpolate.H"
#include "fvmDdt.H"
#include "fvmSup.H"
#include "fvcDiv.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace MULES
{
template<class RhoType>
inline tmp<surfaceScalarField> interpolate(const RhoType& rho)
{
NotImplemented;
return tmp<surfaceScalarField>(nullptr);
}
template<>
inline tmp<surfaceScalarField> interpolate(const volScalarField& rho)
{
return fvc::interpolate(rho);
}
}
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::implicitSolve
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
const scalar psiMin
)
{
const fvMesh& mesh = psi.mesh();
const dictionary& MULEScontrols = mesh.solverDict(psi.name());
label maxIter
(
readLabel(MULEScontrols.lookup("maxIter"))
);
scalar maxUnboundedness
(
readScalar(MULEScontrols.lookup("maxUnboundedness"))
);
scalar CoCoeff
(
readScalar(MULEScontrols.lookup("CoCoeff"))
);
scalarField allCoLambda(mesh.nFaces());
{
slicedSurfaceScalarField CoLambda
(
IOobject
(
"CoLambda",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimless,
allCoLambda,
false // Use slices for the couples
);
if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
{
tmp<surfaceScalarField> Cof =
mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
*mag(phi/interpolate(rho))/mesh.magSf();
CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
}
else
{
tmp<surfaceScalarField> Cof =
mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
*mag(phi)/mesh.magSf();
CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
}
}
scalarField allLambda(allCoLambda);
//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
);
linear<scalar> CDs(mesh);
upwind<scalar> UDs(mesh, phi);
//fv::uncorrectedSnGrad<scalar> snGrads(mesh);
fvScalarMatrix psiConvectionDiffusion
(
fvm::ddt(rho, psi)
+ fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
//- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
//.fvmLaplacian(Dpsif, psi)
- fvm::Sp(Sp, psi)
- Su
);
surfaceScalarField phiBD(psiConvectionDiffusion.flux());
surfaceScalarField& phiCorr = phiPsi;
phiCorr -= phiBD;
for (label i=0; i<maxIter; i++)
{
if (i != 0 && i < 4)
{
allLambda = allCoLambda;
}
limiter
(
allLambda,
1.0/mesh.time().deltaTValue(),
rho,
psi,
phiBD,
phiCorr,
Sp,
Su,
psiMax,
psiMin
);
solve
(
psiConvectionDiffusion + fvc::div(lambda*phiCorr),
MULEScontrols
);
scalar maxPsiM1 = gMax(psi.primitiveField()) - 1.0;
scalar minPsi = gMin(psi.primitiveField());
scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
if (unboundedness < maxUnboundedness)
{
break;
}
else
{
Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
<< " min(" << psi.name() << ") = " << minPsi << endl;
phiBD = psiConvectionDiffusion.flux();
/*
word gammaScheme("div(phi,gamma)");
word gammarScheme("div(phirb,gamma)");
const surfaceScalarField& phir =
mesh.lookupObject<surfaceScalarField>("phir");
phiCorr =
fvc::flux
(
phi,
psi,
gammaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, scalar(1) - psi, gammarScheme),
psi,
gammarScheme
)
- phiBD;
*/
}
}
phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
}
// ************************************************************************* //

92
src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
View File

@ -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
@ -185,12 +185,12 @@ void Foam::MULES::limiter
const dictionary& MULEScontrols = mesh.solverDict(psi.name());
label nLimiterIter
const label nLimiterIter
(
MULEScontrols.lookupOrDefault<label>("nLimiterIter", 3)
);
scalar smoothLimiter
const scalar smoothLimiter
(
MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
);
@ -228,8 +228,7 @@ void Foam::MULES::limiter
);
scalarField& lambdaIf = lambda;
surfaceScalarField::Boundary& lambdaBf =
lambda.boundaryFieldRef();
surfaceScalarField::Boundary& lambdaBf = lambda.boundaryFieldRef();
scalarField psiMaxn(psiIf.size(), psiMin);
scalarField psiMinn(psiIf.size(), psiMax);
@ -241,8 +240,8 @@ void Foam::MULES::limiter
forAll(phiCorrIf, facei)
{
label own = owner[facei];
label nei = neighb[facei];
const label own = owner[facei];
const label nei = neighb[facei];
psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
psiMinn[own] = min(psiMinn[own], psiIf[nei]);
@ -253,9 +252,9 @@ void Foam::MULES::limiter
sumPhiBD[own] += phiBDIf[facei];
sumPhiBD[nei] -= phiBDIf[facei];
scalar phiCorrf = phiCorrIf[facei];
const scalar phiCorrf = phiCorrIf[facei];
if (phiCorrf > 0.0)
if (phiCorrf > 0)
{
sumPhip[own] += phiCorrf;
mSumPhim[nei] += phiCorrf;
@ -281,7 +280,7 @@ void Foam::MULES::limiter
forAll(phiCorrPf, pFacei)
{
label pfCelli = pFaceCells[pFacei];
const label pfCelli = pFaceCells[pFacei];
psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
@ -291,7 +290,7 @@ void Foam::MULES::limiter
{
forAll(phiCorrPf, pFacei)
{
label pfCelli = pFaceCells[pFacei];
const label pfCelli = pFaceCells[pFacei];
psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
@ -300,13 +299,13 @@ void Foam::MULES::limiter
forAll(phiCorrPf, pFacei)
{
label pfCelli = pFaceCells[pFacei];
const label pfCelli = pFaceCells[pFacei];
sumPhiBD[pfCelli] += phiBDPf[pFacei];
scalar phiCorrf = phiCorrPf[pFacei];
const scalar phiCorrf = phiCorrPf[pFacei];
if (phiCorrf > 0.0)
if (phiCorrf > 0)
{
sumPhip[pfCelli] += phiCorrf;
}
@ -376,17 +375,17 @@ void Foam::MULES::limiter
for (int j=0; j<nLimiterIter; j++)
{
sumlPhip = 0.0;
mSumlPhim = 0.0;
sumlPhip = 0;
mSumlPhim = 0;
forAll(lambdaIf, facei)
{
label own = owner[facei];
label nei = neighb[facei];
const label own = owner[facei];
const label nei = neighb[facei];
scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
if (lambdaPhiCorrf > 0.0)
if (lambdaPhiCorrf > 0)
{
sumlPhip[own] += lambdaPhiCorrf;
mSumlPhim[nei] += lambdaPhiCorrf;
@ -407,11 +406,11 @@ void Foam::MULES::limiter
forAll(lambdaPf, pFacei)
{
label pfCelli = pFaceCells[pFacei];
const label pfCelli = pFaceCells[pFacei];
const scalar lambdaPhiCorrf =
lambdaPf[pFacei]*phiCorrfPf[pFacei];
scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
if (lambdaPhiCorrf > 0.0)
if (lambdaPhiCorrf > 0)
{
sumlPhip[pfCelli] += lambdaPhiCorrf;
}
@ -446,7 +445,7 @@ void Foam::MULES::limiter
forAll(lambdaIf, facei)
{
if (phiCorrIf[facei] > 0.0)
if (phiCorrIf[facei] > 0)
{
lambdaIf[facei] = min
(
@ -464,7 +463,6 @@ void Foam::MULES::limiter
}
}
forAll(lambdaBf, patchi)
{
fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
@ -482,9 +480,9 @@ void Foam::MULES::limiter
forAll(lambdaPf, pFacei)
{
label pfCelli = pFaceCells[pFacei];
const label pfCelli = pFaceCells[pFacei];
if (phiCorrfPf[pFacei] > 0.0)
if (phiCorrfPf[pFacei] > 0)
{
lambdaPf[pFacei] =
min(lambdaPf[pFacei], lambdap[pfCelli]);
@ -496,33 +494,6 @@ void Foam::MULES::limiter
}
}
}
else
{
const labelList& pFaceCells =
mesh.boundary()[patchi].faceCells();
const scalarField& phiBDPf = phiBDBf[patchi];
const scalarField& phiCorrPf = phiCorrBf[patchi];
forAll(lambdaPf, pFacei)
{
// Limit outlet faces only
if ((phiBDPf[pFacei] + 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<scalar>());
@ -549,6 +520,19 @@ void Foam::MULES::limit
surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
surfaceScalarField::Boundary& phiBDBf = phiBD.boundaryFieldRef();
const surfaceScalarField::Boundary& phiPsiBf = phiPsi.boundaryField();
forAll(phiBDBf, patchi)
{
fvsPatchScalarField& phiBDPf = phiBDBf[patchi];
if (!phiBDPf.coupled())
{
phiBDPf = phiPsiBf[patchi];
}
}
surfaceScalarField& phiCorr = phiPsi;
phiCorr -= phiBD;