diff --git a/applications/solvers/multiphase/interFoam/Allwclean b/applications/solvers/multiphase/interFoam/Allwclean index 350d4b268b..2c54bfa321 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 b25e4440e6..8044426582 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 0000000000..9fcb19a7cc --- /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 . + +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. + +\*---------------------------------------------------------------------------*/ + +#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 0000000000..729b041aa6 --- /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 . + +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 +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 +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 +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 0000000000..a6d3a3f6b2 --- /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 . + +\*---------------------------------------------------------------------------*/ + +#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 +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(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("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 +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 CDs(mesh); + upwind UDs(mesh, phi); + //fv::uncorrectedSnGrad snGrads(mesh); + + fvScalarMatrix psiConvectionDiffusion + ( + fvm::ddt(rho, psi) + + fv::gaussConvectionScheme(mesh, phi, UDs).fvmDiv(phi, psi) + //- fv::gaussLaplacianScheme(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("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 +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 tVsc = mesh.Vsc(); + const scalarField& V = tVsc(); + + const volScalarField& rDeltaT = + mesh.objectRegistry::lookupObject("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 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 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()); + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files new file mode 100644 index 0000000000..205812e978 --- /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 0000000000..24349f694e --- /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 0000000000..a4cb6e7688 --- /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 1) +{ + dimensionedScalar totalDeltaT = runTime.deltaT(); + surfaceScalarField rhoPhiSum = 0.0*rhoPhi; + + for + ( + subCycle 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 0000000000..1f2f082fb3 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H @@ -0,0 +1,30 @@ +scalar maxDeltaT +( + piso.lookupOrDefault("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 0000000000..f5b8450ddd --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H @@ -0,0 +1,109 @@ +{ + scalar maxCo + ( + piso.lookupOrDefault("maxCo", 0.9) + ); + + scalar maxAlphaCo + ( + piso.lookupOrDefault("maxAlphaCo", 0.2) + ); + + scalar rDeltaTSmoothingCoeff + ( + piso.lookupOrDefault("rDeltaTSmoothingCoeff", 0.1) + ); + + label nAlphaSpreadIter + ( + piso.lookupOrDefault