diff --git a/applications/solvers/modules/VoFSolver/VoFSolver.C b/applications/solvers/modules/VoFSolver/VoFSolver.C index 9769c90aa7..9b2052928c 100644 --- a/applications/solvers/modules/VoFSolver/VoFSolver.C +++ b/applications/solvers/modules/VoFSolver/VoFSolver.C @@ -25,8 +25,9 @@ License #include "VoFSolver.H" #include "localEulerDdtScheme.H" -#include "CorrectPhi.H" -#include "geometricZeroField.H" +#include "linear.H" +#include "fvcDiv.H" +#include "fvcMeshPhi.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // diff --git a/applications/solvers/modules/VoFSolver/VoFSolver.H b/applications/solvers/modules/VoFSolver/VoFSolver.H index 35243f885b..80b8317301 100644 --- a/applications/solvers/modules/VoFSolver/VoFSolver.H +++ b/applications/solvers/modules/VoFSolver/VoFSolver.H @@ -169,9 +169,16 @@ protected: //- Return the pressure reference virtual const Foam::pressureReference& pressureReference() const = 0; + //- Is the flow incompressible? + virtual bool incompressible() const = 0; + //- Is the flow divergent? // i.e. compressible or include phase-fraction sources - virtual bool divergent() = 0; + virtual bool divergent() const = 0; + + //- Return the mixture compressibility/density + // Used by CorrectPhi for compressible mixtures following mesh change + virtual tmp psiByRho() const = 0; //- Return the momentum equation stress term virtual tmp divDevTau(volVectorField& U) = 0; diff --git a/applications/solvers/modules/VoFSolver/moveMesh.C b/applications/solvers/modules/VoFSolver/moveMesh.C index a61004238c..9be50fadbe 100644 --- a/applications/solvers/modules/VoFSolver/moveMesh.C +++ b/applications/solvers/modules/VoFSolver/moveMesh.C @@ -65,29 +65,44 @@ void Foam::solvers::VoFSolver::moveMesh() correctUphiBCs(U_, phi, true); - if (divU.valid()) + if (incompressible()) { - CorrectPhi - ( - phi, - U, - p_rgh, - surfaceScalarField("rAUf", fvc::interpolate(rAU())), - divU(), - pressureReference(), - pimple - ); + if (divU.valid()) + { + CorrectPhi + ( + phi, + U, + p_rgh, + surfaceScalarField("rAUf", fvc::interpolate(rAU())), + divU(), + pressureReference(), + pimple + ); + } + else + { + CorrectPhi + ( + phi, + U, + p_rgh, + surfaceScalarField("rAUf", fvc::interpolate(rAU())), + geometricZeroField(), + pressureReference(), + pimple + ); + } } else { CorrectPhi ( phi, - U, p_rgh, + psiByRho(), surfaceScalarField("rAUf", fvc::interpolate(rAU())), - geometricZeroField(), - pressureReference(), + divU(), pimple ); } diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.C b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.C index b97def57d6..f9bdab5580 100644 --- a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.C +++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.C @@ -24,9 +24,9 @@ License \*---------------------------------------------------------------------------*/ #include "compressibleMultiphaseVoF.H" -#include "CorrectPhi.H" #include "geometricZeroField.H" #include "fvcDdt.H" +#include "fvcDiv.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.H b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.H index f2fa2c6d85..3e1a9a759d 100644 --- a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.H +++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.H @@ -128,10 +128,24 @@ protected: return pressureReference_; } - //- Compressible flow is divergent - virtual bool divergent() + //- The flow is incompressible if all phases are incompressible + virtual bool incompressible() const { - return true; + return mixture.incompressible(); + } + + //- The flow is divergent if it is not incompressible + // Mass sources are not currently supported + virtual bool divergent() const + { + return !incompressible(); + } + + //- Return the mixture compressibility/density + // Used by CorrectPhi for compressible mixtures following mesh change + virtual tmp psiByRho() const + { + return mixture.psiByRho(); } //- Return the momentum equation stress term diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.C b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.C index be4046e187..524d7f757b 100644 --- a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.C +++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.C @@ -67,6 +67,20 @@ Foam::compressibleMultiphaseVoFMixture::compressibleMultiphaseVoFMixture // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // +bool Foam::compressibleMultiphaseVoFMixture::incompressible() const +{ + bool incompressible = true; + + forAll(phases_, phasei) + { + incompressible = + incompressible && phases_[phasei].thermo().incompressible(); + } + + return incompressible; +} + + Foam::tmp Foam::compressibleMultiphaseVoFMixture::nu() const { @@ -103,6 +117,25 @@ Foam::tmp Foam::compressibleMultiphaseVoFMixture::nu } +Foam::tmp +Foam::compressibleMultiphaseVoFMixture::psiByRho() const +{ + tmp tpsiByRho + ( + phases_[0]*phases_[0].thermo().psi()/phases_[0].thermo().rho() + ); + + for (label phasei=1; phasei Foam::compressibleMultiphaseVoFMixture::alphaEff ( const volScalarField& nut diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.H b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.H index 5657d023fb..ce98d54da7 100644 --- a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.H +++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.H @@ -104,6 +104,9 @@ public: return rho_; } + //- Return true if all phases are incompressible + bool incompressible() const; + //- Return the kinematic laminar viscosity virtual tmp nu() const; @@ -113,6 +116,9 @@ public: //- Return the face-interpolated dynamic laminar viscosity tmp nuf() const; + //- Return the mixture compressibility/density + virtual tmp psiByRho() const; + //- Correct the thermodynamics of each phase virtual void correctThermo(); diff --git a/applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.C b/applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.C index 5992fb764e..2b0198bddc 100644 --- a/applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.C +++ b/applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.C @@ -203,6 +203,17 @@ Foam::tmp Foam::compressibleTwoPhaseVoFMixture::nu } +Foam::tmp +Foam::compressibleTwoPhaseVoFMixture::psiByRho() const +{ + return + ( + alpha1()*thermo1_->psi()/thermo1_->rho() + + alpha2()*thermo2_->psi()/thermo2_->rho() + ); +} + + bool Foam::compressibleTwoPhaseVoFMixture::read() { if (twoPhaseVoFMixture::read()) diff --git a/applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.H b/applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.H index f40c0d5a94..84078cd0a0 100644 --- a/applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.H +++ b/applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.H @@ -161,12 +161,23 @@ public: return thermo2_->rho(); } + //- The fluid is incompressible if both phases are incompressible + bool incompressible() const + { + return + thermo1_->incompressible() + && thermo2_->incompressible(); + } + //- Return mixture density [kg/m^3] virtual const volScalarField& rho() const { return rho_; } + //- Return the mixture compressibility/density + virtual tmp psiByRho() const; + //- Correct the thermodynamics of each phase virtual void correctThermo(); diff --git a/applications/solvers/modules/compressibleVoF/compressibleVoF.H b/applications/solvers/modules/compressibleVoF/compressibleVoF.H index 817d1e157b..81237d84b4 100644 --- a/applications/solvers/modules/compressibleVoF/compressibleVoF.H +++ b/applications/solvers/modules/compressibleVoF/compressibleVoF.H @@ -148,10 +148,25 @@ protected: return pressureReference_; } - //- Compressible flow is divergent - virtual bool divergent() + //- The fluid is incompressible if both phases are incompressible + virtual bool incompressible() const { - return true; + return mixture.incompressible(); + } + + //- Compressible flow is divergent + virtual bool divergent() const + { + return + !incompressible() + || fvModels().addsSupToField(alpha1.name()); + } + + //- Return the mixture compressibility/density + // Used by CorrectPhi for compressible mixtures following mesh change + virtual tmp psiByRho() const + { + return mixture.psiByRho(); } //- Calculate the alpha equation sources diff --git a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoF.H b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoF.H index 502a902ee1..86ff65eabc 100644 --- a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoF.H +++ b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoF.H @@ -111,13 +111,26 @@ protected: return pressureReference_; } - //- Is the flow divergent? - // i.e. includes phase-fraction sources - virtual bool divergent() + //- The flow is incompressible + virtual bool incompressible() const + { + return true; + } + + //- The flow is not incompressible and hence not divergent + // Mass sources are not currently supported + virtual bool divergent() const { return false; } + //- Return the mixture compressibility/density + // Not required for incompressible fluids + virtual tmp psiByRho() const + { + return tmp(nullptr); + } + //- Return the momentum equation stress term virtual tmp divDevTau(volVectorField& U) { diff --git a/applications/solvers/modules/incompressibleVoF/incompressibleVoF.H b/applications/solvers/modules/incompressibleVoF/incompressibleVoF.H index 054cc4b0d5..27ac4cdd1b 100644 --- a/applications/solvers/modules/incompressibleVoF/incompressibleVoF.H +++ b/applications/solvers/modules/incompressibleVoF/incompressibleVoF.H @@ -112,13 +112,26 @@ protected: return pressureReference_; } + //- The flow is incompressible + virtual bool incompressible() const + { + return true; + } + //- Is the flow divergent? // i.e. includes phase-fraction sources - virtual bool divergent() + virtual bool divergent() const { return fvModels().addsSupToField(alpha1.name()); } + //- Return the mixture compressibility/density + // Not required for incompressible fluids + virtual tmp psiByRho() const + { + return tmp(nullptr); + } + //- Calculate the alpha equation sources virtual void alphaSuSp ( diff --git a/applications/solvers/modules/isothermalFluid/moveMesh.C b/applications/solvers/modules/isothermalFluid/moveMesh.C index b5ebf5c291..1907420a73 100644 --- a/applications/solvers/modules/isothermalFluid/moveMesh.C +++ b/applications/solvers/modules/isothermalFluid/moveMesh.C @@ -61,7 +61,6 @@ void Foam::solvers::isothermalFluid::moveMesh() ( phi, buoyancy.valid() ? p_rgh : p, - rho, thermo.psi(), dimensionedScalar("rAUf", dimTime, 1), divrhoU(), diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/moveMesh.C b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/moveMesh.C index a843cc4de9..589c276d7e 100644 --- a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/moveMesh.C +++ b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/moveMesh.C @@ -40,6 +40,7 @@ void Foam::solvers::multiphaseEuler::moveMesh() if ( (correctPhi || mesh.topoChanged()) + // && divergent() && !divU.valid() ) { diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C index de3c973dee..a3f09fbd66 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C @@ -774,17 +774,65 @@ void Foam::phaseSystem::correctPhi phi_ += alphafs[phasei]*(mesh_.Sf() & phase.UfRef()); } - CorrectPhi - ( - phi_, - movingPhases()[0].U(), - p_rgh, - // surfaceScalarField("rAUf", fvc::interpolate(rAU())), - dimensionedScalar(dimTime/dimDensity, 1), - divU(), - pressureReference, - pimple - ); + if (incompressible()) + { + if (divU.valid()) + { + CorrectPhi + ( + phi_, + movingPhases()[0].U(), + p_rgh, + dimensionedScalar(dimTime/dimDensity, 1), + divU(), + pressureReference, + pimple + ); + } + else + { + CorrectPhi + ( + phi_, + movingPhases()[0].U(), + p_rgh, + dimensionedScalar(dimTime/dimDensity, 1), + geometricZeroField(), + pressureReference, + pimple + ); + } + } + else + { + volScalarField psi + ( + volScalarField::New + ( + "psi", + mesh_, + dimensionedScalar(dimless/dimPressure, 0) + ) + ); + + forAll(phases(), phasei) + { + phaseModel& phase = phases()[phasei]; + const volScalarField& alpha = phase; + + psi += alpha*phase.thermo().psi()/phase.thermo().rho(); + } + + CorrectPhi + ( + phi_, + p_rgh, + psi, + dimensionedScalar(dimTime/dimDensity, 1), + divU(), + pimple + ); + } // Make the flux relative to the mesh motion fvc::makeRelative(phi_, movingPhases()[0].U()); diff --git a/applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFSolver.C b/applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFSolver.C index c7cd0a71f5..6abce0c831 100644 --- a/applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFSolver.C +++ b/applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFSolver.C @@ -25,8 +25,7 @@ License #include "multiphaseVoFSolver.H" #include "localEulerDdtScheme.H" -#include "CorrectPhi.H" -#include "geometricZeroField.H" +#include "fvcAverage.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // diff --git a/applications/solvers/modules/twoPhaseVoFSolver/twoPhaseVoFSolver.C b/applications/solvers/modules/twoPhaseVoFSolver/twoPhaseVoFSolver.C index ae34b5e64f..69f0bd7461 100644 --- a/applications/solvers/modules/twoPhaseVoFSolver/twoPhaseVoFSolver.C +++ b/applications/solvers/modules/twoPhaseVoFSolver/twoPhaseVoFSolver.C @@ -25,8 +25,7 @@ License #include "twoPhaseVoFSolver.H" #include "localEulerDdtScheme.H" -#include "CorrectPhi.H" -#include "geometricZeroField.H" +#include "fvcAverage.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // diff --git a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C index c34e39bf1f..21250c9fae 100644 --- a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C +++ b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -118,7 +118,6 @@ void Foam::CorrectPhi ( surfaceScalarField& phi, const volScalarField& p, - const volScalarField& rho, const volScalarField& psi, const RAUfType& rAUf, const DivRhoUType& divRhoU, diff --git a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H index 715489da43..21781f7845 100644 --- a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H +++ b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2021 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -87,7 +87,6 @@ namespace Foam ( surfaceScalarField& phi, const volScalarField& p, - const volScalarField& rho, const volScalarField& psi, const RAUfType& rAUf, const DivRhoUType& divRhoU, diff --git a/tutorials/modules/multiphaseEuler/mixerVessel2D/system/fvSolution b/tutorials/modules/multiphaseEuler/mixerVessel2D/system/fvSolution index 0291d786e7..06a89edd6c 100644 --- a/tutorials/modules/multiphaseEuler/mixerVessel2D/system/fvSolution +++ b/tutorials/modules/multiphaseEuler/mixerVessel2D/system/fvSolution @@ -25,9 +25,10 @@ solvers p_rgh { solver GAMG; + smoother DIC; + tolerance 1e-7; relTol 0.01; - smoother GaussSeidel; } p_rghFinal @@ -36,20 +37,30 @@ solvers preconditioner { preconditioner GAMG; + smoother GaussSeidel; + tolerance 1e-7; relTol 0; nVcycles 2; - smoother GaussSeidel; } tolerance 1e-7; relTol 0; - maxIter 30; + maxIter 50; } "pcorr.*" { $p_rghFinal; - tolerance 1e-5; + tolerance 0.001; + relTol 0; + } + + MeshPhi + { + solver smoothSolver; + smoother symGaussSeidel; + + tolerance 1e-2; relTol 0; } @@ -57,6 +68,7 @@ solvers { solver smoothSolver; smoother symGaussSeidel; + tolerance 1e-5; relTol 0; minIter 1; @@ -66,6 +78,7 @@ solvers { solver smoothSolver; smoother symGaussSeidel; + tolerance 1e-8; relTol 0; minIter 1; @@ -79,7 +92,7 @@ PIMPLE nCorrectors 3; nNonOrthogonalCorrectors 0; - correctPhi no; + correctPhi yes; correctMeshPhi no; }