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

LTSInterFoam: Initial version of interFoam supporting local time-stepping for acceleration to steady-state

Browse Source
This commit is contained in:
Henry
2010-10-25 17:49:56 +01:00
parent a67a967d9a
commit 730481aed3
12 changed files with 1153 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/solvers/multiphase/interFoam/Allwclean
View File

@ -6,5 +6,6 @@ wclean
wclean interDyMFoam wclean interDyMFoam
wclean MRFInterFoam wclean MRFInterFoam
wclean porousInterFoam wclean porousInterFoam
wclean LTSInterFoam
# ----------------------------------------------------------------- end-of-file # ----------------------------------------------------------------- end-of-file

1
applications/solvers/multiphase/interFoam/Allwmake
View File

@ -6,5 +6,6 @@ wmake
wmake interDyMFoam wmake interDyMFoam
wmake MRFInterFoam wmake MRFInterFoam
wmake porousInterFoam wmake porousInterFoam
wmake LTSInterFoam
# ----------------------------------------------------------------- end-of-file # ----------------------------------------------------------------- end-of-file

101
applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C Normal file
View File

@ -0,0 +1,101 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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/>.
Application
interFoam
Description
Solver for 2 incompressible, isothermal immiscible fluids using a VOF
(volume of fluid) phase-fraction based interface capturing approach.
The momentum and other fluid properties are of the "mixture" and a single
momentum equation is solved.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
For a two-fluid approach see twoPhaseEulerFoam.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "turbulenceModel.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialrDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readPISOControls.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "setrDeltaT.H"
twoPhaseProperties.correct();
#include "alphaEqnSubCycle.H"
turbulence->correct();
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

83
applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C Normal file
View File

@ -0,0 +1,83 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 "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"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::MULES::explicitLTSSolve
(
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
const scalar psiMax,
const scalar psiMin
)
{
explicitLTSSolve
(
geometricOneField(),
psi,
phi,
phiPsi,
zeroField(), zeroField(),
psiMax, psiMin
);
}
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
);
}
// ************************************************************************* //

136
applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H Normal file
View File

@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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
MULES
Description
MULES: Multidimensional universal limiter with explicit solution.
Solve a convective-only transport equation using an explicit universal
multi-dimensional limiter.
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
\*---------------------------------------------------------------------------*/
#ifndef MULES_H
#define MULES_H
#include "volFields.H"
#include "surfaceFieldsFwd.H"
#include "primitiveFieldsFwd.H"
#include "zeroField.H"
#include "geometricOneField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace MULES
{
template<class RhoType, class SpType, class SuType>
void explicitLTSSolve
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phiBD,
surfaceScalarField& phiPsi,
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
);
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
);
template<class RhoType, class SpType, class SuType>
void limiter
(
scalarField& allLambda,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phiBD,
const 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"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

602
applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C Normal file
View File

@ -0,0 +1,602 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 "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"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitLTSSolve
(
const RhoType& rho,
volScalarField& psi,
const surfaceScalarField& phi,
surfaceScalarField& phiPsi,
const SpType& Sp,
const SuType& Su,
const scalar psiMax,
const scalar psiMin
)
{
Info<< "MULES: Solving for " << 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);
if (mesh.moving())
{
psiIf =
(
mesh.Vsc0()*rho.oldTime()*psi0*rDeltaT/mesh.Vsc()
+ Su.field()
- psiIf
)/(rho*rDeltaT - Sp.field());
}
else
{
psiIf =
(
rho.oldTime()*psi0*rDeltaT
+ Su.field()
- psiIf
)/(rho*rDeltaT - Sp.field());
}
psi.correctBoundaryConditions();
}
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"))
);
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;
}
template<class RhoType, class SpType, class SuType>
void Foam::MULES::limiter
(
scalarField& allLambda,
const RhoType& rho,
const volScalarField& psi,
const surfaceScalarField& phiBD,
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 scalarField& psi0 = psi.oldTime();
const fvMesh& mesh = psi.mesh();
const unallocLabelList& owner = mesh.owner();
const unallocLabelList& neighb = mesh.neighbour();
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();
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 sumPhiBD(psiIf.size(), 0.0);
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]);
sumPhiBD[own] += phiBDIf[facei];
sumPhiBD[nei] -= phiBDIf[facei];
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& phiBDPf = phiBDBf[patchi];
const scalarField& phiCorrPf = phiCorrBf[patchi];
const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
if (psiPf.coupled())
{
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];
sumPhiBD[pfCelli] += phiBDPf[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);
if (mesh.moving())
{
tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
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;
}
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 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]);
}
}
}
syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
}
}
// ************************************************************************* //

4
applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files Normal file
View File

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

15
applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options Normal file
View File

@ -0,0 +1,15 @@
EXE_INC = \
-I.. \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume

36
applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H Normal file
View File

@ -0,0 +1,36 @@
{
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, scalar(1) - alpha1, alpharScheme),
alpha1,
alpharScheme
);
MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
//MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
}
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}

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

@ -0,0 +1,35 @@
label nAlphaCorr
(
readLabel(piso.lookup("nAlphaCorr"))
);
label nAlphaSubCycles
(
readLabel(piso.lookup("nAlphaSubCycles"))
);
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
# include "alphaEqn.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
# include "alphaEqn.H"
}
interface.correct();
rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;

30
applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H Normal file
View File

@ -0,0 +1,30 @@
scalar maxDeltaT
(
piso.lookupOrDefault<scalar>("maxDeltaT", GREAT)
);
volScalarField rDeltaT
(
IOobject
(
"rDeltaT",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
);
volScalarField rSubDeltaT
(
IOobject
(
"rSubDeltaT",
runTime.timeName(),
mesh
),
mesh,
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
);

109
applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H Normal file
View File

@ -0,0 +1,109 @@
{
scalar maxCo
(
piso.lookupOrDefault<scalar>("maxCo", 0.9)
);
scalar maxAlphaCo
(
piso.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
);
scalar rDeltaTSmoothingCoeff
(
piso.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
);
label nAlphaSpreadIter
(
piso.lookupOrDefault<label>("nAlphaSpreadIter", 1)
);
label nAlphaSweepIter
(
piso.lookupOrDefault<label>("nAlphaSweepIter", 5)
);
scalar rDeltaTDampingCoeff
(
piso.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
);
scalar maxDeltaT
(
piso.lookupOrDefault<scalar>("maxDeltaT", GREAT)
);
volScalarField rDeltaT0 = rDeltaT;
// Set the reciprocal time-step using an effective maximum Courant number
rDeltaT = max
(
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
fvc::surfaceSum
(
mag(rhoPhi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf())
)/rho
);
// Limit the time-step further in the region of the interface
{
surfaceScalarField alphaf = fvc::interpolate(alpha1);
surfaceScalarField SfUfbyDelta =
pos(alphaf - 0.01)*pos(0.99 - alphaf)
*mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
rDeltaT = max
(
rDeltaT,
fvc::surfaceSum(mag(SfUfbyDelta/(maxAlphaCo*mesh.magSf())))
);
}
Info<< "Flow time scale min/max = "
<< gMin(1/rDeltaT.internalField())
<< ", " << gMax(1/rDeltaT.internalField()) << endl;
if (rDeltaTSmoothingCoeff < 1.0)
{
fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
}
if (nAlphaSpreadIter > 0)
{
fvc::spread(rDeltaT, alpha1, nAlphaSpreadIter);
}
if (nAlphaSweepIter > 0)
{
fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter);
}
Info<< "Flow time scale min/max = "
<< gMin(1/rDeltaT.internalField())
<< ", " << gMax(1/rDeltaT.internalField()) << endl;
// Limit rate of change of time scale
// - reduce as much as required
// - only increase at a fraction of old time scale
if
(
rDeltaTDampingCoeff < 1.0
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
{
Info<< "Damping rDeltaT" << endl;
rDeltaT = rDeltaT0*max(rDeltaT/rDeltaT0, 1.0 - rDeltaTDampingCoeff);
}
Info<< "Flow time scale min/max = "
<< gMin(1/rDeltaT.internalField())
<< ", " << gMax(1/rDeltaT.internalField()) << endl;
label nAlphaSubCycles
(
readLabel(piso.lookup("nAlphaSubCycles"))
);
rSubDeltaT = rDeltaT*nAlphaSubCycles;
}