From 730481aed3e00a66a8a85ecee4643824e42db5fd Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 25 Oct 2010 17:49:56 +0100
Subject: [PATCH] LTSInterFoam: Initial version of interFoam supporting local
time-stepping for acceleration to steady-state
---
.../solvers/multiphase/interFoam/Allwclean | 1 +
.../solvers/multiphase/interFoam/Allwmake | 1 +
.../interFoam/LTSInterFoam/LTSInterFoam.C | 101 +++
.../multiphase/interFoam/LTSInterFoam/MULES.C | 83 +++
.../multiphase/interFoam/LTSInterFoam/MULES.H | 136 ++++
.../interFoam/LTSInterFoam/MULESTemplates.C | 602 ++++++++++++++++++
.../interFoam/LTSInterFoam/Make/files | 4 +
.../interFoam/LTSInterFoam/Make/options | 15 +
.../interFoam/LTSInterFoam/alphaEqn.H | 36 ++
.../interFoam/LTSInterFoam/alphaEqnSubCycle.H | 35 +
.../LTSInterFoam/setInitialrDeltaT.H | 30 +
.../interFoam/LTSInterFoam/setrDeltaT.H | 109 ++++
12 files changed, 1153 insertions(+)
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
diff --git a/applications/solvers/multiphase/interFoam/Allwclean b/applications/solvers/multiphase/interFoam/Allwclean
index 350d4b268b0..2c54bfa3218 100755
--- a/applications/solvers/multiphase/interFoam/Allwclean
+++ b/applications/solvers/multiphase/interFoam/Allwclean
@@ -6,5 +6,6 @@ wclean
wclean interDyMFoam
wclean MRFInterFoam
wclean porousInterFoam
+wclean LTSInterFoam
# ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/interFoam/Allwmake b/applications/solvers/multiphase/interFoam/Allwmake
index b25e4440e60..8044426582e 100755
--- a/applications/solvers/multiphase/interFoam/Allwmake
+++ b/applications/solvers/multiphase/interFoam/Allwmake
@@ -6,5 +6,6 @@ wmake
wmake interDyMFoam
wmake MRFInterFoam
wmake porousInterFoam
+wmake LTSInterFoam
# ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
new file mode 100644
index 00000000000..9fcb19a7cc2
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@@ -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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
new file mode 100644
index 00000000000..e281ae575cb
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
@@ -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
+ );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
new file mode 100644
index 00000000000..729b041aa62
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
@@ -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
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
new file mode 100644
index 00000000000..a6d3a3f6b24
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
@@ -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>());
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
new file mode 100644
index 00000000000..205812e978d
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
@@ -0,0 +1,4 @@
+LTSInterFoam.C
+MULES.C
+
+EXE = $(FOAM_APPBIN)/LTSInterFoam
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options
new file mode 100644
index 00000000000..24349f694e0
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options
@@ -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
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
new file mode 100644
index 00000000000..a4cb6e7688c
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
@@ -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;
+}
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
new file mode 100644
index 00000000000..9aae37a8bef
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
@@ -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;
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
new file mode 100644
index 00000000000..1f2f082fb36
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
@@ -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)
+);
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
new file mode 100644
index 00000000000..f5b8450ddd9
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
@@ -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;
+}
--
GitLab