zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-12-28 03:37:59 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
mattijs
2013-10-29 09:18:25 +00:00
parent f5d9d59436 b994d6a11f
commit 95c06af87e
40 changed files with 573 additions and 958 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
View File

@ -986,6 +986,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
MULES::limit
(
1.0/mesh_.time().deltaT().value(),
geometricOneField(),
alpha,
phi_,

5
applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
View File

@ -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();

1
applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
View File

@ -1,4 +1,3 @@
LTSInterFoam.C
MULES.C
EXE = $(FOAM_APPBIN)/LTSInterFoam

39
applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
View File

@ -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<nAlphaCorr; aCorr++)
{
surfaceScalarField phiAlpha
(
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
);
MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
//MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
alpha2 = 1.0 - alpha1;
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
}
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}

35
applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
View File

@ -1,35 +0,0 @@
#include "alphaControls.H"
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum
(
IOobject
(
"rhoPhiSum",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", rhoPhi.dimensions(), 0)
);
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
#include "alphaEqn.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
#include "alphaEqn.H"
}
rho == alpha1*rho1 + alpha2*rho2;

2
applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
View File

@ -37,7 +37,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "incompressibleTwoPhaseMixture.H"

56
applications/solvers/multiphase/interFoam/alphaEqn.H
View File

@ -6,9 +6,39 @@
phic = min(interface.cAlpha()*phic, max(phic));
surfaceScalarField phir(phic*interface.nHatf());
tmp<surfaceScalarField> tphiAlpha;
if (MULESCorr)
{
fvScalarMatrix alpha1Eqn
(
#ifdef LTSSOLVE
fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
#else
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
#endif
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi,
upwind<scalar>(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<nAlphaCorr; aCorr++)
{
surfaceScalarField phiAlpha
tmp<surfaceScalarField> 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()

2
applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
View File

@ -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"

2
applications/solvers/multiphase/interFoam/interFoam.C
View File

@ -38,7 +38,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "incompressibleTwoPhaseMixture.H"

2
applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H
View File

@ -75,6 +75,7 @@
MULES::limiter
(
allLambda,
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha1,
phiAlpha1BD,
@ -116,6 +117,7 @@
MULES::limiter
(
allLambda,
1.0/runTime.deltaT().value(),
geometricOneField(),
alpha2,
phiAlpha2BD,

2
applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
View File

@ -31,7 +31,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "CMULES.H"
#include "subCycle.H"
#include "threePhaseMixture.H"
#include "threePhaseInterfaceProperties.H"

2
applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
View File

@ -39,7 +39,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "incompressibleTwoPhaseMixture.H"

10
applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
View File

@ -16,8 +16,14 @@
{
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1)
+ fvm::div(phi, alpha1, "UD") - fvm::Sp(divU, alpha1)
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi,
upwind<scalar>(mesh, phi)
).fvmDiv(phi, alpha1)
- fvm::Sp(divU, alpha1)
==
fvm::Sp(vDotvmcAlphal, alpha1)
+ vDotcAlphal

2
applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
View File

@ -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"

2
applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
View File

@ -41,7 +41,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "phaseChangeTwoPhaseMixture.H"

1
applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
View File

@ -148,6 +148,7 @@ void Foam::multiphaseSystem::solveAlphas()
MULES::limit
(
1.0/mesh_.time().deltaT().value(),
geometricOneField(),
phase1,
phi_,

1
applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
View File

@ -598,6 +598,7 @@ void Foam::multiphaseMixture::solveAlphas
MULES::limit
(
1.0/mesh_.time().deltaT().value(),
geometricOneField(),
alpha,
phi_,

1
src/finiteVolume/Make/files
View File

@ -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

39
src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C
View File

@ -38,7 +38,7 @@ interstitialInletVelocityFvPatchVectorField
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchField<vector>(p, iF),
fixedValueFvPatchVectorField(p, iF),
inletVelocity_(p.size(), vector::zero),
alphaName_("alpha")
{}
@ -53,8 +53,8 @@ interstitialInletVelocityFvPatchVectorField
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
inletVelocity_(ptf.inletVelocity_),
fixedValueFvPatchVectorField(ptf, p, iF, mapper),
inletVelocity_(ptf.inletVelocity_, mapper),
alphaName_(ptf.alphaName_)
{}
@ -67,7 +67,7 @@ interstitialInletVelocityFvPatchVectorField
const dictionary& dict
)
:
fixedValueFvPatchField<vector>(p, iF, dict),
fixedValueFvPatchVectorField(p, iF, dict),
inletVelocity_("inletVelocity", dict, p.size()),
alphaName_(dict.lookupOrDefault<word>("alpha", "alpha"))
{}
@ -79,7 +79,7 @@ interstitialInletVelocityFvPatchVectorField
const interstitialInletVelocityFvPatchVectorField& ptf
)
:
fixedValueFvPatchField<vector>(ptf),
fixedValueFvPatchVectorField(ptf),
inletVelocity_(ptf.inletVelocity_),
alphaName_(ptf.alphaName_)
{}
@ -92,7 +92,7 @@ interstitialInletVelocityFvPatchVectorField
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchField<vector>(ptf, iF),
fixedValueFvPatchVectorField(ptf, iF),
inletVelocity_(ptf.inletVelocity_),
alphaName_(ptf.alphaName_)
{}
@ -100,6 +100,31 @@ interstitialInletVelocityFvPatchVectorField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::interstitialInletVelocityFvPatchVectorField::autoMap
(
const fvPatchFieldMapper& m
)
{
fixedValueFvPatchVectorField::autoMap(m);
inletVelocity_.autoMap(m);
}
void Foam::interstitialInletVelocityFvPatchVectorField::rmap
(
const fvPatchVectorField& ptf,
const labelList& addr
)
{
fixedValueFvPatchVectorField::rmap(ptf, addr);
const interstitialInletVelocityFvPatchVectorField& tiptf =
refCast<const interstitialInletVelocityFvPatchVectorField>(ptf);
inletVelocity_.rmap(tiptf.inletVelocity_, addr);
}
void Foam::interstitialInletVelocityFvPatchVectorField::updateCoeffs()
{
if (updated())
@ -111,7 +136,7 @@ void Foam::interstitialInletVelocityFvPatchVectorField::updateCoeffs()
patch().lookupPatchField<volScalarField, scalar>(alphaName_);
operator==(inletVelocity_/alphap);
fixedValueFvPatchField<vector>::updateCoeffs();
fixedValueFvPatchVectorField::updateCoeffs();
}

22
src/finiteVolume/fields/fvPatchFields/derived/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.H
View File

@ -141,8 +141,26 @@ public:
// Member functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
// Mapping functions
//- Map (and resize as needed) from self given a mapping object
virtual void autoMap
(
const fvPatchFieldMapper&
);
//- Reverse map the given fvPatchField onto this fvPatchField
virtual void rmap
(
const fvPatchVectorField&,
const labelList&
);
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;

8
src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
View File

@ -198,7 +198,6 @@ tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
template<>
tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtPhiCorr
(
@ -206,6 +205,13 @@ tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
template<>
tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtUfCorr
(
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& Uf
);
template<>
tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtPhiCorr

8
src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
View File

@ -270,7 +270,6 @@ tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
template<>
tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
(
@ -278,6 +277,13 @@ tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
template<>
tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtUfCorr
(
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& Uf
);
template<>
tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr

8
src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
View File

@ -176,7 +176,6 @@ tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
template<>
tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtPhiCorr
(
@ -184,6 +183,13 @@ tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
template<>
tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtUfCorr
(
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& Uf
);
template<>
tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtPhiCorr

8
src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
View File

@ -199,7 +199,6 @@ tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
template<>
tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtPhiCorr
(
@ -207,6 +206,13 @@ tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
template<>
tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtUfCorr
(
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& Uf
);
template<>
tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtPhiCorr

8
src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H
View File

@ -187,7 +187,6 @@ tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
template<>
tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtPhiCorr
(
@ -195,6 +194,13 @@ tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
template<>
tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtUfCorr
(
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& Uf
);
template<>
tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtPhiCorr

8
src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H
View File

@ -185,7 +185,6 @@ tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
template<>
tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtPhiCorr
(
@ -193,6 +192,13 @@ tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
template<>
tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtUfCorr
(
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& Uf
);
template<>
tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtPhiCorr

15
src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H
View File

@ -91,6 +91,13 @@ public:
// Constructors
//- Construct from mesh and name of the rDeltaT field
localEulerDdtScheme(const fvMesh& mesh, const word& rDeltaTName)
:
ddtScheme<Type>(mesh),
rDeltaTName_(rDeltaTName)
{}
//- Construct from mesh and Istream
localEulerDdtScheme(const fvMesh& mesh, Istream& is)
:
@ -188,7 +195,6 @@ tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
template<>
tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
(
@ -196,6 +202,13 @@ tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
template<>
tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtUfCorr
(
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& Uf
);
template<>
tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr

8
src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H
View File

@ -175,7 +175,6 @@ tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtUfCorr
const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
);
template<>
tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtPhiCorr
(
@ -183,6 +182,13 @@ tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtPhiCorr
const surfaceScalarField& phi
);
template<>
tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtUfCorr
(
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& Uf
);
template<>
tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtPhiCorr

35
applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C → src/finiteVolume/fvMatrices/solvers/MULES/CMULES.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,57 +23,42 @@ License
\*---------------------------------------------------------------------------*/
#include "MULES.H"
#include "upwind.H"
#include "uncorrectedSnGrad.H"
#include "gaussConvectionScheme.H"
#include "gaussLaplacianScheme.H"
#include "uncorrectedSnGrad.H"
#include "surfaceInterpolate.H"
#include "fvcSurfaceIntegrate.H"
#include "slicedSurfaceFields.H"
#include "syncTools.H"
#include "fvCFD.H"
#include "CMULES.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::MULES::explicitLTSSolve
void Foam::MULES::correct
(
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
surfaceScalarField& phiPsiCorr,
const scalar psiMax,
const scalar psiMin
)
{
explicitLTSSolve
correct
(
geometricOneField(),
psi,
phi,
phiPsi,
phiPsiCorr,
zeroField(), zeroField(),
psiMax, psiMin
);
}
void Foam::MULES::implicitSolve
void Foam::MULES::LTScorrect
(
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
surfaceScalarField& phiPsiCorr,
const scalar psiMax,
const scalar psiMin
)
{
implicitSolve
LTScorrect
(
geometricOneField(),
psi,
phi,
phiPsi,
phiPsiCorr,
zeroField(), zeroField(),
psiMax, psiMin
);

104
applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H → src/finiteVolume/fvMatrices/solvers/MULES/CMULES.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,67 +25,57 @@ Global
MULES
Description
MULES: Multidimensional universal limiter with explicit solution.
CMULES: Multidimensional universal limiter for
explicit corrected implicit solution.
Solve a convective-only transport equation using an explicit universal
multi-dimensional limiter.
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
MULES.C
CMULES.C
CMULESTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef MULES_H
#define MULES_H
#ifndef CMULES_H
#define CMULES_H
#include "volFields.H"
#include "surfaceFieldsFwd.H"
#include "primitiveFieldsFwd.H"
#include "zeroField.H"
#include "geometricOneField.H"
#include "MULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "gaussConvectionScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace MULES
{
template<class RhoType, class SpType, class SuType>
void explicitLTSSolve
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void correct
(
const RdeltaTType& rDeltaT,
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiBD,
surfaceScalarField& phiPsi,
const surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
const scalar psiMin
);
void explicitLTSSolve
(
volScalarField& psi,
const surfaceScalarField& phiBD,
surfaceScalarField& phiPsi,
const scalar psiMax,
const scalar psiMin
const SuType& Su
);
template<class RhoType, class SpType, class SuType>
void implicitSolve
void correct
(
const RhoType& rho,
volScalarField& gamma,
const surfaceScalarField& phi,
volScalarField& psi,
surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
@ -93,22 +83,42 @@ void implicitSolve
const scalar psiMin
);
void implicitSolve
void correct
(
volScalarField& gamma,
const surfaceScalarField& phi,
volScalarField& psi,
surfaceScalarField& phiCorr,
const scalar psiMax,
const scalar psiMin
);
template<class RhoType, class SpType, class SuType>
void limiter
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<class RdeltaTType, class RhoType, class SpType, class SuType>
void limiterCorr
(
scalarField& allLambda,
const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phiBD,
const surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
@ -117,16 +127,30 @@ void limiter
const label nLimiterIter
);
} // End namespace MULES
template<class RdeltaTType, class RhoType, class SpType, class SuType>
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 "MULESTemplates.C"
# include "CMULESTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

461
applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C → 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) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,105 +23,49 @@ License
\*---------------------------------------------------------------------------*/
#include "MULES.H"
#include "upwind.H"
#include "uncorrectedSnGrad.H"
#include "gaussConvectionScheme.H"
#include "gaussLaplacianScheme.H"
#include "uncorrectedSnGrad.H"
#include "surfaceInterpolate.H"
#include "CMULES.H"
#include "fvcSurfaceIntegrate.H"
#include "slicedSurfaceFields.H"
#include "wedgeFvPatch.H"
#include "syncTools.H"
#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitLTSSolve
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::correct
(
const RdeltaTType& rDeltaT,
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
const surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
const scalar psiMin
const SuType& Su
)
{
Info<< "MULES: Solving for " << psi.name() << endl;
Info<< "MULES: Correcting " << psi.name() << endl;
const fvMesh& mesh = psi.mesh();
psi.correctBoundaryConditions();
surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
surfaceScalarField& phiCorr = phiPsi;
phiCorr -= phiBD;
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
);
limiter
(
allLambda,
rho,
psi,
phiBD,
phiCorr,
Sp.field(),
Su.field(),
psiMax,
psiMin,
3
);
phiPsi = phiBD + lambda*phiCorr;
scalarField& psiIf = psi;
const scalarField& psi0 = psi.oldTime();
const volScalarField& rDeltaT =
mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
psiIf = 0.0;
fvc::surfaceIntegrate(psiIf, phiPsi);
scalarField psiIf(psi.size(), 0);
fvc::surfaceIntegrate(psiIf, phiCorr);
if (mesh.moving())
{
psiIf =
psi.internalField() =
(
mesh.Vsc0()*rho.oldTime()*psi0*rDeltaT/mesh.Vsc()
rho.field()*psi.internalField()*rDeltaT
+ Su.field()
- psiIf
)/(rho*rDeltaT - Sp.field());
)/(rho.field()*rDeltaT - Sp.field());
}
else
{
psiIf =
psi.internalField() =
(
rho.oldTime()*psi0*rDeltaT
rho.field()*psi.internalField()*rDeltaT
+ Su.field()
- psiIf
)/(rho*rDeltaT - Sp.field());
)/(rho.field()*rDeltaT - Sp.field());
}
psi.correctBoundaryConditions();
@ -129,12 +73,11 @@ void Foam::MULES::explicitLTSSolve
template<class RhoType, class SpType, class SuType>
void Foam::MULES::implicitSolve
void Foam::MULES::correct
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
@ -142,176 +85,56 @@ void Foam::MULES::implicitSolve
)
{
const fvMesh& mesh = psi.mesh();
const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
const dictionary& MULEScontrols = mesh.solverDict(psi.name());
label maxIter
(
readLabel(MULEScontrols.lookup("maxIter"))
);
label nLimiterIter
(
readLabel(MULEScontrols.lookup("nLimiterIter"))
);
scalar maxUnboundedness
(
readScalar(MULEScontrols.lookup("maxUnboundedness"))
);
scalar CoCoeff
(
readScalar(MULEScontrols.lookup("CoCoeff"))
);
scalarField allCoLambda(mesh.nFaces());
{
surfaceScalarField Cof
(
mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
*mag(phi)/mesh.magSf()
);
slicedSurfaceScalarField CoLambda
(
IOobject
(
"CoLambda",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimless,
allCoLambda,
false // Use slices for the couples
);
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,
rho,
psi,
phiBD,
phiCorr,
Sp.field(),
Su.field(),
psiMax,
psiMin,
nLimiterIter
);
solve
(
psiConvectionDiffusion + fvc::div(lambda*phiCorr),
MULEScontrols
);
scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
scalar minPsi = gMin(psi.internalField());
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;
limitCorr(rDeltaT, rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
correct(rDeltaT, rho, psi, phiCorr, Sp, Su);
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::limiter
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<volScalarField>("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<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::limiterCorr
(
scalarField& allLambda,
const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phiBD,
const surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su,
@ -323,8 +146,6 @@ void Foam::MULES::limiter
const scalarField& psiIf = psi;
const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
const scalarField& psi0 = psi.oldTime();
const fvMesh& mesh = psi.mesh();
const labelUList& owner = mesh.owner();
@ -332,13 +153,6 @@ void Foam::MULES::limiter
tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
const scalarField& V = tVsc();
const volScalarField& rDeltaT =
mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
const scalarField& phiBDIf = phiBD;
const surfaceScalarField::GeometricBoundaryField& phiBDBf =
phiBD.boundaryField();
const scalarField& phiCorrIf = phiCorr;
const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
phiCorr.boundaryField();
@ -367,8 +181,6 @@ void Foam::MULES::limiter
scalarField psiMaxn(psiIf.size(), psiMin);
scalarField psiMinn(psiIf.size(), psiMax);
scalarField sumPhiBD(psiIf.size(), 0.0);
scalarField sumPhip(psiIf.size(), VSMALL);
scalarField mSumPhim(psiIf.size(), VSMALL);
@ -383,9 +195,6 @@ void Foam::MULES::limiter
psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
sumPhiBD[own] += phiBDIf[facei];
sumPhiBD[nei] -= phiBDIf[facei];
scalar phiCorrf = phiCorrIf[facei];
if (phiCorrf > 0.0)
@ -403,14 +212,13 @@ void Foam::MULES::limiter
forAll(phiCorrBf, patchi)
{
const fvPatchScalarField& psiPf = psiBf[patchi];
const scalarField& phiBDPf = phiBDBf[patchi];
const scalarField& phiCorrPf = phiCorrBf[patchi];
const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
if (psiPf.coupled())
{
scalarField psiPNf(psiPf.patchNeighbourField());
const scalarField psiPNf(psiPf.patchNeighbourField());
forAll(phiCorrPf, pFacei)
{
@ -435,8 +243,6 @@ void Foam::MULES::limiter
{
label pfCelli = pFaceCells[pFacei];
sumPhiBD[pfCelli] += phiBDPf[pFacei];
scalar phiCorrf = phiCorrPf[pFacei];
if (phiCorrf > 0.0)
@ -453,39 +259,30 @@ void Foam::MULES::limiter
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);
// scalar smooth = 0.5;
// psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
// psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
if (mesh.moving())
{
tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
psiMaxn =
V
*(
(rho.field()*rDeltaT - Sp.field())*psiMaxn
- Su.field()
- rho.field()*psi.internalField()*rDeltaT
);
psiMaxn =
V*((rho*rDeltaT - Sp)*psiMaxn - Su)
- (V0()*rDeltaT)*rho.oldTime()*psi0
+ sumPhiBD;
psiMinn =
V*(Su - (rho*rDeltaT - Sp)*psiMinn)
+ (V0*rDeltaT)*rho.oldTime()*psi0
- sumPhiBD;
}
else
{
psiMaxn =
V*((rho*rDeltaT - Sp)*psiMaxn - (rho.oldTime()*rDeltaT)*psi0 - Su)
+ sumPhiBD;
psiMinn =
V*((rho*rDeltaT)*psi0 - (rho.oldTime()*rDeltaT - Sp)*psiMinn + Su)
- sumPhiBD;
}
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<nLimiterIter; j++)
for (int j=0; j<nLimiterIter; j++)
{
sumlPhip = 0.0;
mSumlPhim = 0.0;
@ -533,19 +330,21 @@ void Foam::MULES::limiter
}
}
forAll (sumlPhip, celli)
forAll(sumlPhip, celli)
{
sumlPhip[celli] =
max(min
(
(sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
(sumlPhip[celli] + psiMaxn[celli])
/(mSumPhim[celli] - SMALL),
1.0), 0.0
);
mSumlPhim[celli] =
max(min
(
(mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
(mSumlPhim[celli] + psiMinn[celli])
/(sumPhip[celli] + SMALL),
1.0), 0.0
);
}
@ -578,20 +377,57 @@ void Foam::MULES::limiter
{
fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
const scalarField& phiCorrfPf = phiCorrBf[patchi];
const fvPatchScalarField& psiPf = psiBf[patchi];
const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
forAll(lambdaPf, pFacei)
if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
{
label pfCelli = pFaceCells[pFacei];
lambdaPf = 0;
}
else if (psiPf.coupled())
{
const labelList& pFaceCells =
mesh.boundary()[patchi].faceCells();
if (phiCorrfPf[pFacei] > 0.0)
forAll(lambdaPf, pFacei)
{
lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]);
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
}
else
{
const labelList& pFaceCells =
mesh.boundary()[patchi].faceCells();
const scalarField& phiCorrPf = phiCorrBf[patchi];
forAll(lambdaPf, pFacei)
{
lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]);
// 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]);
}
}
}
}
}
@ -601,4 +437,57 @@ void Foam::MULES::limiter
}
template<class RdeltaTType, class RhoType, class SpType, class SuType>
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;
}
// ************************************************************************* //

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

@ -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

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

@ -179,6 +179,7 @@ void Foam::MULES::implicitSolve
limiter
(
allLambda,
1.0/mesh.time().deltaTValue(),
rho,
psi,
phiBD,

10
src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
View File

@ -48,19 +48,21 @@ void Foam::MULES::explicitSolve
}
void Foam::MULES::correct
void Foam::MULES::explicitLTSSolve
(
volScalarField& psi,
surfaceScalarField& phiPsiCorr,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
const scalar psiMax,
const scalar psiMin
)
{
correct
explicitLTSSolve
(
geometricOneField(),
psi,
phiPsiCorr,
phi,
phiPsi,
zeroField(), zeroField(),
psiMax, psiMin
);

69
src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
View File

@ -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.
@ -60,6 +60,17 @@ namespace MULES
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void explicitSolve
(
const RdeltaTType& rDeltaT,
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiPsi,
const SpType& Sp,
const SuType& Su
);
template<class RhoType, class SpType, class SuType>
void explicitSolve
(
@ -92,42 +103,34 @@ void explicitSolve
const scalar psiMin
);
template<class RhoType, class SpType, class SuType>
void correct
void explicitLTSSolve
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiCorr,
const SpType& Sp,
const SuType& Su
);
template<class RhoType, class SpType, class SuType>
void correct
(
const RhoType& rho,
volScalarField& psi,
surfaceScalarField& phiCorr,
const surfaceScalarField& phiBD,
surfaceScalarField& phiPsi,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
const scalar psiMin
);
void correct
void explicitLTSSolve
(
volScalarField& psi,
surfaceScalarField& phiCorr,
const surfaceScalarField& phiBD,
surfaceScalarField& phiPsi,
const scalar psiMax,
const scalar psiMin
);
template<class RhoType, class SpType, class SuType>
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void limiter
(
scalarField& allLambda,
const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phiBD,
@ -139,9 +142,10 @@ void limiter
const label nLimiterIter
);
template<class RhoType, class SpType, class SuType>
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void limit
(
const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phi,
@ -155,39 +159,12 @@ void limit
);
template<class RhoType, class SpType, class SuType>
void limiterCorr
(
scalarField& allLambda,
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<class RhoType, class SpType, class SuType>
void limitCorr
(
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<scalarField>& phiPsiCorrs);
template<class SurfaceScalarFieldList>
void limitSum(SurfaceScalarFieldList& phiPsiCorrs);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace MULES

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

@ -32,9 +32,10 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class RhoType, class SpType, class SuType>
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
const RdeltaTType& rDeltaT,
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiPsi,
@ -48,7 +49,6 @@ void Foam::MULES::explicitSolve
scalarField& psiIf = psi;
const scalarField& psi0 = psi.oldTime();
const scalar deltaT = mesh.time().deltaTValue();
psiIf = 0.0;
fvc::surfaceIntegrate(psiIf, phiPsi);
@ -58,25 +58,41 @@ void Foam::MULES::explicitSolve
psiIf =
(
mesh.Vsc0()().field()*rho.oldTime().field()
*psi0/(deltaT*mesh.Vsc()().field())
*psi0*rDeltaT/mesh.Vsc()().field()
+ Su.field()
- psiIf
)/(rho.field()/deltaT - Sp.field());
)/(rho.field()*rDeltaT - Sp.field());
}
else
{
psiIf =
(
rho.oldTime().field()*psi0/deltaT
rho.oldTime().field()*psi0*rDeltaT
+ Su.field()
- psiIf
)/(rho.field()/deltaT - Sp.field());
)/(rho.field()*rDeltaT - Sp.field());
}
psi.correctBoundaryConditions();
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiPsi,
const SpType& Sp,
const SuType& Su
)
{
const fvMesh& mesh = psi.mesh();
const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
@ -90,60 +106,21 @@ void Foam::MULES::explicitSolve
const scalar psiMin
)
{
psi.correctBoundaryConditions();
limit(rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
explicitSolve(rho, psi, phiPsi, Sp, Su);
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::correct
(
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();
const scalar deltaT = mesh.time().deltaTValue();
scalarField psiIf(psi.size(), 0);
fvc::surfaceIntegrate(psiIf, phiCorr);
if (mesh.moving())
{
psi.internalField() =
(
rho.field()*psi.internalField()/deltaT
+ Su.field()
- psiIf
)/(rho.field()/deltaT - Sp.field());
}
else
{
psi.internalField() =
(
rho.field()*psi.internalField()/deltaT
+ Su.field()
- psiIf
)/(rho.field()/deltaT - Sp.field());
}
const scalar rDeltaT = 1.0/mesh.time().deltaTValue();
psi.correctBoundaryConditions();
limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::correct
void Foam::MULES::explicitLTSSolve
(
const RhoType& rho,
volScalarField& psi,
surfaceScalarField& phiCorr,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
@ -152,22 +129,20 @@ void Foam::MULES::correct
{
const fvMesh& mesh = psi.mesh();
const dictionary& MULEScontrols = mesh.solverDict(psi.name());
const volScalarField& rDeltaT =
mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
label nLimiterIter
(
readLabel(MULEScontrols.lookup("nLimiterIter"))
);
limitCorr(rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
correct(rho, psi, phiCorr, Sp, Su);
psi.correctBoundaryConditions();
limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);
explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);
}
template<class RhoType, class SpType, class SuType>
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::limiter
(
scalarField& allLambda,
const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phiBD,
@ -190,7 +165,6 @@ void Foam::MULES::limiter
const labelUList& neighb = mesh.neighbour();
tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
const scalarField& V = tVsc();
const scalar deltaT = mesh.time().deltaTValue();
const scalarField& phiBDIf = phiBD;
const surfaceScalarField::GeometricBoundaryField& phiBDBf =
@ -321,19 +295,19 @@ void Foam::MULES::limiter
psiMaxn =
V
*(
(rho.field()/deltaT - Sp.field())*psiMaxn
(rho.field()*rDeltaT - Sp.field())*psiMaxn
- Su.field()
)
- (V0().field()/deltaT)*rho.oldTime().field()*psi0
- (V0().field()*rDeltaT)*rho.oldTime().field()*psi0
+ sumPhiBD;
psiMinn =
V
*(
Su.field()
- (rho.field()/deltaT - Sp.field())*psiMinn
- (rho.field()*rDeltaT - Sp.field())*psiMinn
)
+ (V0().field()/deltaT)*rho.oldTime().field()*psi0
+ (V0().field()*rDeltaT)*rho.oldTime().field()*psi0
- sumPhiBD;
}
else
@ -341,9 +315,9 @@ void Foam::MULES::limiter
psiMaxn =
V
*(
(rho.field()/deltaT - Sp.field())*psiMaxn
(rho.field()*rDeltaT - Sp.field())*psiMaxn
- Su.field()
- (rho.oldTime().field()/deltaT)*psi0
- (rho.oldTime().field()*rDeltaT)*psi0
)
+ sumPhiBD;
@ -351,8 +325,8 @@ void Foam::MULES::limiter
V
*(
Su.field()
- (rho.field()/deltaT - Sp.field())*psiMinn
+ (rho.oldTime().field()/deltaT)*psi0
- (rho.field()*rDeltaT - Sp.field())*psiMinn
+ (rho.oldTime().field()*rDeltaT)*psi0
)
- sumPhiBD;
}
@ -413,14 +387,16 @@ void Foam::MULES::limiter
sumlPhip[celli] =
max(min
(
(sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
(sumlPhip[celli] + psiMaxn[celli])
/(mSumPhim[celli] - SMALL),
1.0), 0.0
);
mSumlPhim[celli] =
max(min
(
(mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
(mSumlPhim[celli] + psiMinn[celli])
/(sumPhip[celli] + SMALL),
1.0), 0.0
);
}
@ -514,9 +490,10 @@ void Foam::MULES::limiter
}
template<class RhoType, class SpType, class SuType>
template<class RdeltaTType, class RhoType, class SpType, class SuType>
void Foam::MULES::limit
(
const RdeltaTType& rDeltaT,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phi,
@ -558,6 +535,7 @@ void Foam::MULES::limit
limiter
(
allLambda,
rDeltaT,
rho,
psi,
phiBD,
@ -580,364 +558,6 @@ void Foam::MULES::limit
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::limiterCorr
(
scalarField& allLambda,
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<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
const scalarField& V = tVsc();
const scalar deltaT = mesh.time().deltaTValue();
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()/deltaT - Sp.field())*psiMaxn
- Su.field()
- rho.field()*psi.internalField()/deltaT
);
psiMinn =
V
*(
Su.field()
- (rho.field()/deltaT - Sp.field())*psiMinn
+ rho.field()*psi.internalField()/deltaT
);
scalarField sumlPhip(psiIf.size());
scalarField mSumlPhim(psiIf.size());
for (int j=0; j<nLimiterIter; j++)
{
sumlPhip = 0.0;
mSumlPhim = 0.0;
forAll(lambdaIf, facei)
{
label own = owner[facei];
label nei = neighb[facei];
scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
if (lambdaPhiCorrf > 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],
1.0), 0.0
);
mSumlPhim[celli] =
max(min
(
(mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
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<wedgeFvPatch>(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<scalar>());
}
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::limitCorr
(
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,
rho,
psi,
phiCorr,
Sp,
Su,
psiMax,
psiMin,
nLimiterIter
);
phiCorr *= lambda;
}
template<class SurfaceScalarFieldList>
void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
{

9
tutorials/multiphase/interFoam/ras/waterChannel/system/controlDict
View File

@ -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

11
tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes
View File

@ -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;
}

43
tutorials/multiphase/interFoam/ras/waterChannel/system/fvSolution
View File

@ -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;
}
}

1
tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
View File

@ -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;
}