From 74b302d6f8f247107c2f0a870c964c0e8d2d2eab Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Sun, 18 Dec 2022 17:28:11 +0000 Subject: [PATCH] solvers::compressibleVoF: Implemented new energy conservative temperature correction equation In order to ensure temperature consistency between the phases it is necessary to solve for the mixture temperature rather than the mixture energy or phase energies which makes it very difficult to conserve energy. The new temperature equation is a temperature correction on the combined phase energy equations which will conserve the phase and mixture energies at convergence. The heat-flux (Laplacian) term is maintained in mixture temperature form so heat-transfer boundary conditions, in particular for CHT, remain in terms of the mixture kappaEff. The fvModels are applied to the phase energy equations and the implicit part converted into an implicit term in the temperature correction part of the equation to improve convergence and stability. This development has required some change to the alphaEqn.H and interFoam has been updated for consistency in preparation for conversion into the solvers::incompressibleVoF modular module. All compressibleVoF fvModels and tutorial cases have been updated for the above change. Note that two entries are now required for the convection terms in the temperature equation, one for explicit phase energy terms and another for the implicit phase temperature correction terms, e.g. tutorials/modules/compressibleVoF/ballValve div(alphaRhoPhi,e) Gauss limitedLinear 1; div(alphaRhoPhi,T) Gauss upwind; In the above the upwind scheme is selected for the phase temperature correction terms as they are corrections and will converge to a zero contribution. However there may be cases which converge better if the same scheme is used for both the energy and temperature terms, more testing is required. --- .../modules/compressibleVoF/alphaPredictor.C | 72 ++++++++++--------- .../compressibleInterPhaseTransportModel.C | 55 +++----------- .../compressibleInterPhaseTransportModel.H | 10 ++- .../modules/compressibleVoF/compressibleVoF.C | 16 ++++- .../modules/compressibleVoF/compressibleVoF.H | 15 +++- .../fvModels/VoFClouds/VoFClouds.C | 32 +++++++-- .../fvModels/VoFClouds/VoFClouds.H | 13 +++- .../VoFSolidificationMeltingSource.C | 31 ++++---- .../VoFSolidificationMeltingSource.H | 11 ++- .../fvModels/VoFSurfaceFilm/VoFSurfaceFilm.C | 33 ++++++--- .../fvModels/VoFSurfaceFilm/VoFSurfaceFilm.H | 13 +++- .../compressibleVoF/momentumPredictor.C | 3 +- .../compressibleVoF/pressureCorrector.C | 2 - .../solvers/modules/compressibleVoF/rhofs.H | 2 - .../compressibleVoF/thermophysicalPredictor.C | 54 +++++++++++--- .../multiphase/interFoam/alphaEqnSubCycle.H | 50 +++++++++++++ .../incompressibleInterPhaseTransportModel.C | 10 +-- .../incompressibleInterPhaseTransportModel.H | 6 +- src/twoPhaseModels/VoF/alphaEqn.H | 18 ++--- src/twoPhaseModels/VoF/alphaEqnSubCycle.H | 41 ----------- .../immiscibleCompressibleTwoPhaseMixture.C | 1 - .../immiscibleIncompressibleTwoPhaseMixture.C | 1 - .../system/fluid/fvSchemes | 3 +- .../ballValve/system/fvSchemes | 5 +- .../climbingRod/system/fvSchemes | 8 ++- .../compressibleVoF/cylinder/system/fvSchemes | 8 ++- .../compressibleVoF/damBreak/system/fvSchemes | 5 +- .../depthCharge2D/system/fvSchemes | 5 +- .../depthCharge2D/system/fvSolution | 4 +- .../depthCharge3D/system/fvSchemes | 5 +- .../depthCharge3D/system/fvSolution | 2 +- .../plateFilm/system/fvSchemes | 8 ++- .../sloshingTank2D/system/fvSchemes | 5 +- 33 files changed, 320 insertions(+), 227 deletions(-) delete mode 100644 applications/solvers/modules/compressibleVoF/rhofs.H create mode 100644 applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H delete mode 100644 src/twoPhaseModels/VoF/alphaEqnSubCycle.H diff --git a/applications/solvers/modules/compressibleVoF/alphaPredictor.C b/applications/solvers/modules/compressibleVoF/alphaPredictor.C index 13b6562868..335182aa60 100644 --- a/applications/solvers/modules/compressibleVoF/alphaPredictor.C +++ b/applications/solvers/modules/compressibleVoF/alphaPredictor.C @@ -36,8 +36,6 @@ void Foam::solvers::compressibleVoF::alphaPredictor() { #include "alphaControls.H" - volScalarField& alpha2(mixture.alpha2()); - const volScalarField& rho1 = mixture.thermo1().rho(); const volScalarField& rho2 = mixture.thermo2().rho(); @@ -47,30 +45,6 @@ void Foam::solvers::compressibleVoF::alphaPredictor() { dimensionedScalar totalDeltaT = runTime.deltaT(); - talphaPhi1 = new surfaceScalarField - ( - IOobject - ( - "alphaPhi1", - runTime.name(), - mesh - ), - mesh, - dimensionedScalar(alphaPhi1.dimensions(), 0) - ); - - surfaceScalarField rhoPhiSum - ( - IOobject - ( - "rhoPhiSum", - runTime.name(), - mesh - ), - mesh, - dimensionedScalar(rhoPhi.dimensions(), 0) - ); - tmp trSubDeltaT; if (LTS) @@ -79,31 +53,59 @@ void Foam::solvers::compressibleVoF::alphaPredictor() fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles); } + // Create a temporary alphaPhi1 to accumulate the sub-cycled alphaPhi1 + tmp talphaPhi1 + ( + surfaceScalarField::New + ( + "alphaPhi1", + mesh, + dimensionedScalar(alphaPhi1.dimensions(), 0) + ) + ); + + // Sub-cycle on both alpha1 and alpha2 + List alphaPtrs({&alpha1, &alpha2}); + for ( - subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + subCycle alphaSubCycle + ( + alphaPtrs, + nAlphaSubCycles + ); !(++alphaSubCycle).end(); ) { #include "alphaEqn.H" talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1; - rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; } alphaPhi1 = talphaPhi1(); - rhoPhi = rhoPhiSum; } else { #include "alphaEqn.H" } - contErr = - ( - fvc::ddt(rho)()() + fvc::div(rhoPhi)()() - - (fvModels().source(alpha1, rho1)&rho1)() - - (fvModels().source(alpha2, rho2)&rho2)() - ); + mixture.correct(); + + alphaRhoPhi1 = fvc::interpolate(rho1)*alphaPhi1; + alphaRhoPhi2 = fvc::interpolate(rho2)*(phi - alphaPhi1); + + rhoPhi = alphaRhoPhi1 + alphaRhoPhi2; + + contErr1 = + ( + fvc::ddt(alpha1, rho1)()() + fvc::div(alphaRhoPhi1)()() + - (fvModels().source(alpha1, rho1)&rho1)() + ); + + contErr2 = + ( + fvc::ddt(alpha2, rho2)()() + fvc::div(alphaRhoPhi2)()() + - (fvModels().source(alpha2, rho2)&rho2)() + ); } diff --git a/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C b/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C index c81bdefc3e..8c7eb7ba27 100644 --- a/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C +++ b/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C @@ -34,13 +34,17 @@ Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel const surfaceScalarField& phi, const surfaceScalarField& rhoPhi, const surfaceScalarField& alphaPhi1, + const surfaceScalarField& alphaRhoPhi1, + const surfaceScalarField& alphaRhoPhi2, const compressibleTwoPhaseMixture& mixture ) : twoPhaseTransport_(false), mixture_(mixture), phi_(phi), - alphaPhi1_(alphaPhi1) + alphaPhi1_(alphaPhi1), + alphaRhoPhi1_(alphaRhoPhi1), + alphaRhoPhi2_(alphaRhoPhi2) { { IOdictionary momentumTransport @@ -68,38 +72,14 @@ Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel if (twoPhaseTransport_) { - const volScalarField& alpha1(mixture_.alpha1()); - const volScalarField& alpha2(mixture_.alpha2()); - - const volScalarField& rho1 = mixture_.thermo1().rho(); - const volScalarField& rho2 = mixture_.thermo2().rho(); - - alphaRhoPhi1_ = - ( - new surfaceScalarField - ( - IOobject::groupName("alphaRhoPhi", alpha1.group()), - fvc::interpolate(rho1)*alphaPhi1_ - ) - ); - - alphaRhoPhi2_ = - ( - new surfaceScalarField - ( - IOobject::groupName("alphaRhoPhi", alpha2.group()), - fvc::interpolate(rho2)*(phi_ - alphaPhi1_) - ) - ); - momentumTransport1_ = ( phaseCompressible::momentumTransportModel::New ( - alpha1, - rho1, + mixture_.alpha1(), + mixture_.thermo1().rho(), U, - alphaRhoPhi1_(), + alphaRhoPhi1_, phi, mixture.thermo1() ) @@ -109,10 +89,10 @@ Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel ( phaseCompressible::momentumTransportModel::New ( - alpha2, - rho2, + mixture_.alpha2(), + mixture_.thermo2().rho(), U, - alphaRhoPhi2_(), + alphaRhoPhi2_, phi, mixture.thermo2() ) @@ -154,19 +134,6 @@ Foam::compressibleInterPhaseTransportModel::divDevTau } -void Foam::compressibleInterPhaseTransportModel::correctPhasePhi() -{ - if (twoPhaseTransport_) - { - const volScalarField& rho1 = mixture_.thermo1().rho(); - const volScalarField& rho2 = mixture_.thermo2().rho(); - - alphaRhoPhi1_.ref() = fvc::interpolate(rho1)*alphaPhi1_; - alphaRhoPhi2_.ref() = fvc::interpolate(rho2)*(phi_ - alphaPhi1_); - } -} - - void Foam::compressibleInterPhaseTransportModel::predict() { if (twoPhaseTransport_) diff --git a/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H b/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H index bc3d779386..7dc4188e0b 100644 --- a/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H +++ b/applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H @@ -77,10 +77,10 @@ class compressibleInterPhaseTransportModel const surfaceScalarField& alphaPhi1_; //- Phase-1 mass-flux (constructed for two-phase transport) - tmp alphaRhoPhi1_; + const surfaceScalarField& alphaRhoPhi1_; //- Phase-2 mass-flux (constructed for two-phase transport) - tmp alphaRhoPhi2_; + const surfaceScalarField& alphaRhoPhi2_; //- Mixture transport model (constructed for mixture transport) autoPtr mixtureMomentumTransport_; @@ -106,6 +106,8 @@ public: const surfaceScalarField& phi, const surfaceScalarField& rhoPhi, const surfaceScalarField& alphaPhi1, + const surfaceScalarField& alphaRhoPhi1, + const surfaceScalarField& alphaRhoPhi2, const compressibleTwoPhaseMixture& mixture ); @@ -121,10 +123,6 @@ public: //- Return the effective momentum stress divergence tmp divDevTau(volVectorField& U) const; - //- Correct the phase mass-fluxes - // (required for the two-phase transport option) - void correctPhasePhi(); - //- Predict the phase or mixture transport models void predict(); diff --git a/applications/solvers/modules/compressibleVoF/compressibleVoF.C b/applications/solvers/modules/compressibleVoF/compressibleVoF.C index 951e952989..b75c780762 100644 --- a/applications/solvers/modules/compressibleVoF/compressibleVoF.C +++ b/applications/solvers/modules/compressibleVoF/compressibleVoF.C @@ -96,6 +96,7 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh) mixture(U), alpha1(mixture.alpha1()), + alpha2(mixture.alpha2()), alphaRestart ( @@ -169,6 +170,18 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh) phi*fvc::interpolate(alpha1) ), + alphaRhoPhi1 + ( + IOobject::groupName("alphaRhoPhi", alpha1.group()), + fvc::interpolate(mixture.thermo1().rho())*alphaPhi1 + ), + + alphaRhoPhi2 + ( + IOobject::groupName("alphaRhoPhi", alpha2.group()), + fvc::interpolate(mixture.thermo2().rho())*(phi - alphaPhi1) + ), + K("K", 0.5*magSqr(U)), momentumTransport @@ -178,6 +191,8 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh) phi, rhoPhi, alphaPhi1, + alphaRhoPhi1, + alphaRhoPhi2, mixture ), @@ -330,7 +345,6 @@ void Foam::solvers::compressibleVoF::prePredictor() { fvModels().correct(); alphaPredictor(); - momentumTransport.correctPhasePhi(); if (pimple.predictTransport()) { diff --git a/applications/solvers/modules/compressibleVoF/compressibleVoF.H b/applications/solvers/modules/compressibleVoF/compressibleVoF.H index 4f340c6e30..5eea866d6b 100644 --- a/applications/solvers/modules/compressibleVoF/compressibleVoF.H +++ b/applications/solvers/modules/compressibleVoF/compressibleVoF.H @@ -95,9 +95,12 @@ protected: //- The compressible two-phase mixture immiscibleCompressibleTwoPhaseMixture mixture; - //- Reference to the primary phase-fraction + //- Reference to the phase1-fraction volScalarField& alpha1; + //- Reference to the phase2-fraction + volScalarField& alpha2; + //- Switch indicating if this is a restart bool alphaRestart; @@ -137,6 +140,12 @@ protected: // Phase-1 volumetric flux surfaceScalarField alphaPhi1; + // Phase-1 mass-flux + surfaceScalarField alphaRhoPhi1; + + // Phase-2 mass-flux + surfaceScalarField alphaRhoPhi2; + //- Kinetic energy field // Used in the energy equation volScalarField K; @@ -164,7 +173,9 @@ protected: tmp rAU; - tmp contErr; + tmp contErr1; + + tmp contErr2; //- MULES Correction tmp talphaPhi1Corr0; diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C b/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C index 5f91f13f28..0f4a8991be 100644 --- a/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C +++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.C @@ -89,7 +89,7 @@ Foam::fv::VoFClouds::VoFClouds Foam::wordList Foam::fv::VoFClouds::addSupFields() const { - return wordList({thermo_.rho()().name(), "U", "T"}); + return wordList({thermo_.rho()().name(), thermo_.he().name(), "U"}); } @@ -108,7 +108,7 @@ void Foam::fv::VoFClouds::correct() void Foam::fv::VoFClouds::addSup ( - const volScalarField& rho, + const volScalarField& alpha, fvMatrix& eqn, const word& fieldName ) const @@ -122,13 +122,31 @@ void Foam::fv::VoFClouds::addSup { eqn += clouds_.Srho(); } - else if (fieldName == "T") + else { - const volScalarField::Internal Cv(thermo_.Cv()); + FatalErrorInFunction + << "Support for field " << fieldName << " is not implemented" + << exit(FatalError); + } +} - eqn += - clouds_.Sh(eqn.psi())()/Cv - + clouds_.Srho()*(eqn.psi() - thermo_.he()/Cv); + +void Foam::fv::VoFClouds::addSup +( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const word& fieldName +) const +{ + if (debug) + { + Info<< type() << ": applying source to " << eqn.psi().name() << endl; + } + + if (fieldName == thermo_.he().name()) + { + eqn += clouds_.Sh(eqn.psi()); } else { diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H b/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H index ed5c571da1..9a9ae34924 100644 --- a/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H +++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFClouds/VoFClouds.H @@ -126,15 +126,24 @@ public: // Add explicit and implicit contributions to compressible equation - //- Add explicit contribution to compressible enthalpy equation + //- Add explicit contribution to phase continuity virtual void addSup ( + const volScalarField& alpha, + fvMatrix& eqn, + const word& fieldName + ) const; + + //- Add explicit contribution to phase energy equation + virtual void addSup + ( + const volScalarField& alpha, const volScalarField& rho, fvMatrix& eqn, const word& fieldName ) const; - //- Add implicit contribution to compressible momentum equation + //- Add implicit contribution to mixture momentum equation virtual void addSup ( const volScalarField& rho, diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C b/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C index 5d0bba9327..834f765fbd 100644 --- a/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C +++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C @@ -92,6 +92,15 @@ Foam::fv::VoFSolidificationMeltingSource::VoFSolidificationMeltingSource relax_(NaN), Cu_(NaN), q_(NaN), + + thermo_ + ( + mesh().lookupObject + ( + "phaseProperties" + ) + ), + alphaSolid_ ( IOobject @@ -115,12 +124,13 @@ Foam::fv::VoFSolidificationMeltingSource::VoFSolidificationMeltingSource Foam::wordList Foam::fv::VoFSolidificationMeltingSource::addSupFields() const { - return wordList({"U", "T"}); + return wordList({"U", thermo_.thermo1().he().name()}); } void Foam::fv::VoFSolidificationMeltingSource::addSup ( + const volScalarField& alpha, const volScalarField& rho, fvMatrix& eqn, const word& fieldName @@ -131,24 +141,7 @@ void Foam::fv::VoFSolidificationMeltingSource::addSup Info<< type() << ": applying source to " << eqn.psi().name() << endl; } - const compressibleTwoPhaseMixture& thermo - ( - mesh().lookupObject - ( - "phaseProperties" - ) - ); - - const volScalarField CpVoF(thermo.thermo1().Cp()); - - if (eqn.psi().dimensions() == dimTemperature) - { - eqn += L_/CpVoF*(fvc::ddt(rho, alphaSolid_)); - } - else - { - eqn += L_*(fvc::ddt(rho, alphaSolid_)); - } + eqn += L_*(fvc::ddt(rho, alphaSolid_)); } diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H b/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H index 3ad77e8564..deaa07b63f 100644 --- a/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H +++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H @@ -88,6 +88,9 @@ SourceFiles namespace Foam { + +class compressibleTwoPhaseMixture; + namespace fv { @@ -119,6 +122,9 @@ class VoFSolidificationMeltingSource //- Coefficient used in porosity calc - default = 0.001 scalar q_; + //- VoF thermo + const compressibleTwoPhaseMixture& thermo_; + //- Solid phase fraction mutable volScalarField alphaSolid_; @@ -171,15 +177,16 @@ public: // Add explicit and implicit contributions to compressible equation - //- Add explicit contribution to compressible enthalpy equation + //- Add explicit contribution to phase energy equation virtual void addSup ( + const volScalarField& alpha, const volScalarField& rho, fvMatrix& eqn, const word& fieldName ) const; - //- Add implicit contribution to compressible momentum equation + //- Add implicit contribution to mixture momentum equation virtual void addSup ( const volScalarField& rho, diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.C b/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.C index 4c66ed7596..ccadc1f3e9 100644 --- a/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.C +++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.C @@ -78,7 +78,7 @@ Foam::wordList Foam::fv::VoFSurfaceFilm::addSupFields() const { surfaceFilm_.rhoPrimary().name(), surfaceFilm_.UPrimary().name(), - surfaceFilm_.TPrimary().name() + surfaceFilm_.primaryThermo().he().name() } ); } @@ -105,7 +105,7 @@ void Foam::fv::VoFSurfaceFilm::correct() void Foam::fv::VoFSurfaceFilm::addSup ( - const volScalarField& rho, + const volScalarField& alpha, fvMatrix& eqn, const word& fieldName ) const @@ -119,14 +119,31 @@ void Foam::fv::VoFSurfaceFilm::addSup { eqn += surfaceFilm_.Srho(); } - else if (fieldName == surfaceFilm_.TPrimary().name()) + else { - const volScalarField::Internal Cv(surfaceFilm_.primaryThermo().Cv()); + FatalErrorInFunction + << "Support for field " << fieldName << " is not implemented" + << exit(FatalError); + } +} - eqn += - surfaceFilm_.Sh()()/Cv - + surfaceFilm_.Srho() - *(eqn.psi() - surfaceFilm_.primaryThermo().he()/Cv); + +void Foam::fv::VoFSurfaceFilm::addSup +( + const volScalarField& alpha, + const volScalarField& rho, + fvMatrix& eqn, + const word& fieldName +) const +{ + if (debug) + { + Info<< type() << ": applying source to " << eqn.psi().name() << endl; + } + + if (fieldName == surfaceFilm_.primaryThermo().he().name()) + { + eqn += surfaceFilm_.Sh(); } else { diff --git a/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.H b/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.H index 8a6d0c9353..33eadabeb5 100644 --- a/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.H +++ b/applications/solvers/modules/compressibleVoF/fvModels/VoFSurfaceFilm/VoFSurfaceFilm.H @@ -116,15 +116,24 @@ public: // Add explicit and implicit contributions to compressible equation - //- Add explicit contribution to compressible enthalpy equation + //- Add explicit contribution to phase continuity virtual void addSup ( + const volScalarField& alpha, + fvMatrix& eqn, + const word& fieldName + ) const; + + //- Add explicit contribution to phase energy equation + virtual void addSup + ( + const volScalarField& alpha, const volScalarField& rho, fvMatrix& eqn, const word& fieldName ) const; - //- Add implicit contribution to compressible momentum equation + //- Add implicit contribution to mixture momentum equation virtual void addSup ( const volScalarField& rho, diff --git a/applications/solvers/modules/compressibleVoF/momentumPredictor.C b/applications/solvers/modules/compressibleVoF/momentumPredictor.C index 273895bb81..e221c34d3b 100644 --- a/applications/solvers/modules/compressibleVoF/momentumPredictor.C +++ b/applications/solvers/modules/compressibleVoF/momentumPredictor.C @@ -31,7 +31,8 @@ void Foam::solvers::compressibleVoF::momentumPredictor() { tUEqn = ( - fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - fvm::Sp(contErr(), U) + fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + - fvm::Sp(contErr1() + contErr2(), U) + MRF.DDt(rho, U) + momentumTransport.divDevTau(U) == diff --git a/applications/solvers/modules/compressibleVoF/pressureCorrector.C b/applications/solvers/modules/compressibleVoF/pressureCorrector.C index 14b3169a84..5295844ab1 100644 --- a/applications/solvers/modules/compressibleVoF/pressureCorrector.C +++ b/applications/solvers/modules/compressibleVoF/pressureCorrector.C @@ -34,8 +34,6 @@ void Foam::solvers::compressibleVoF::pressureCorrector() { volScalarField& p = mixture.p(); - const volScalarField& alpha2(mixture.alpha2()); - const volScalarField& rho1 = mixture.rho1(); const volScalarField& rho2 = mixture.rho2(); diff --git a/applications/solvers/modules/compressibleVoF/rhofs.H b/applications/solvers/modules/compressibleVoF/rhofs.H deleted file mode 100644 index 67dc10f378..0000000000 --- a/applications/solvers/modules/compressibleVoF/rhofs.H +++ /dev/null @@ -1,2 +0,0 @@ -surfaceScalarField rho1f(fvc::interpolate(rho1)); -surfaceScalarField rho2f(fvc::interpolate(rho2)); diff --git a/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C b/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C index 108d7bfe74..169f1598a9 100644 --- a/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C +++ b/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C @@ -29,29 +29,61 @@ License void Foam::solvers::compressibleVoF::thermophysicalPredictor() { + const volScalarField& rho1(mixture.rho1()); + const volScalarField& rho2(mixture.rho2()); + + const volScalarField& e1(mixture.thermo1().he()); + const volScalarField& e2(mixture.thermo2().he()); + + const fvScalarMatrix e1Source(fvModels().source(alpha1, rho1, e1)); + const fvScalarMatrix e2Source(fvModels().source(alpha2, rho2, e2)); + volScalarField& T = mixture.T(); const volScalarField& p = mixture.p(); - const volScalarField& alpha2(mixture.alpha2()); fvScalarMatrix TEqn ( - fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr(), T) - - fvm::laplacian(thermophysicalTransport.alphaEff(), T) + correction + ( + mixture.thermo1().Cv()() + *( + fvm::ddt(alpha1, rho1, T) + fvm::div(alphaRhoPhi1, T) + - ( + e1Source.hasDiag() + ? fvm::Sp(contErr1(), T) + fvm::Sp(e1Source.A(), T) + : fvm::Sp(contErr1(), T) + ) + ) + + mixture.thermo2().Cv()() + *( + fvm::ddt(alpha2, rho2, T) + fvm::div(alphaRhoPhi2, T) + - ( + e2Source.hasDiag() + ? fvm::Sp(contErr2(), T) + fvm::Sp(e2Source.A(), T) + : fvm::Sp(contErr2(), T) + ) + ) + ) + + + fvc::ddt(alpha1, rho1, e1) + fvc::div(alphaRhoPhi1, e1) + - contErr1()*e1 + + fvc::ddt(alpha2, rho2, e2) + fvc::div(alphaRhoPhi2, e2) + - contErr2()*e2 + + - fvm::laplacian(thermophysicalTransport.kappaEff(), T) + + ( - mixture.totalInternalEnergy() + mixture.totalInternalEnergy() ? - fvc::div(fvc::absolute(phi, U), p)()() // - contErr()/rho*p + fvc::div(fvc::absolute(phi, U), p)()() + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - - (U()&(fvModels().source(rho, U)&U)()) - contErr()*K + - (U()&(fvModels().source(rho, U)&U)()) - (contErr1() + contErr2())*K : p*fvc::div(fvc::absolute(phi, U))()() ) - *( - alpha1()/mixture.thermo1().Cv()() - + alpha2()/mixture.thermo2().Cv()() - ) == - fvModels().source(rho, T) + (e1Source&e1) + + (e2Source&e2) ); TEqn.relax(); diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H new file mode 100644 index 0000000000..a5477a0058 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H @@ -0,0 +1,50 @@ +if (nAlphaSubCycles > 1) +{ + dimensionedScalar totalDeltaT = runTime.deltaT(); + tmp trSubDeltaT; + + if (LTS) + { + trSubDeltaT = + fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles); + } + + // Create a temporary alphaPhi1 to accumulate the sub-cycled alphaPhi1 + tmp talphaPhi1 + ( + surfaceScalarField::New + ( + "alphaPhi1", + mesh, + dimensionedScalar(alphaPhi1.dimensions(), 0) + ) + ); + + List alphaPtrs({&alpha1, &alpha2}); + + for + ( + subCycle alphaSubCycle + ( + alphaPtrs, + nAlphaSubCycles + ); + !(++alphaSubCycle).end(); + ) + { + #include "alphaEqn.H" + talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1; + } + + alphaPhi1 = talphaPhi1(); +} +else +{ + #include "alphaEqn.H" +} + +mixture.correct(); +rho == mixture.rho(); + +// Calculate the mass-flux from the accumulated alphaPhi1 +rhoPhi = (alphaPhi1*(rho1 - rho2) + phi*rho2); diff --git a/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C b/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C index 96a2a22ae4..72ade6824b 100644 --- a/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C +++ b/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C @@ -32,14 +32,14 @@ incompressibleInterPhaseTransportModel ( const volVectorField& U, const surfaceScalarField& phi, - const surfaceScalarField& alphaPhi10, + const surfaceScalarField& alphaPhi1, const incompressibleTwoPhaseMixture& mixture ) : twoPhaseTransport_(false), mixture_(mixture), phi_(phi), - alphaPhi10_(alphaPhi10) + alphaPhi1_(alphaPhi1) { { IOdictionary momentumTransport @@ -75,7 +75,7 @@ incompressibleInterPhaseTransportModel new surfaceScalarField ( IOobject::groupName("alphaPhi", alpha2.group()), - (phi_ - alphaPhi10_) + (phi_ - alphaPhi1_) ) ); @@ -85,7 +85,7 @@ incompressibleInterPhaseTransportModel ( alpha1, U, - alphaPhi10_, + alphaPhi1_, phi, mixture.nuModel1() ) @@ -143,7 +143,7 @@ void Foam::incompressibleInterPhaseTransportModel::correctPhasePhi() { if (twoPhaseTransport_) { - alphaPhi2_.ref() = (phi_ - alphaPhi10_); + alphaPhi2_.ref() = (phi_ - alphaPhi1_); } } diff --git a/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.H b/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.H index e8a7cb09ef..69b9759e4b 100644 --- a/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.H +++ b/applications/solvers/multiphase/interFoam/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.H @@ -70,9 +70,9 @@ class incompressibleInterPhaseTransportModel const surfaceScalarField& phi_; //- Phase volumetric flux - const surfaceScalarField& alphaPhi10_; + const surfaceScalarField& alphaPhi1_; - //- Phase-2 mass-flux (constructed from alphaPhi10_ and phi_) + //- Phase-2 mass-flux (constructed from alphaPhi1_ and phi_) tmp alphaPhi2_; //- Mixture transport model (constructed for mixture transport) @@ -94,7 +94,7 @@ public: ( const volVectorField& U, const surfaceScalarField& phi, - const surfaceScalarField& alphaPhi10, + const surfaceScalarField& alphaPhi1, const incompressibleTwoPhaseMixture& mixture ); diff --git a/src/twoPhaseModels/VoF/alphaEqn.H b/src/twoPhaseModels/VoF/alphaEqn.H index c4ac7d876c..2dbf2ff896 100644 --- a/src/twoPhaseModels/VoF/alphaEqn.H +++ b/src/twoPhaseModels/VoF/alphaEqn.H @@ -215,7 +215,8 @@ alpha2 = 1.0 - alpha1; - mixture.correct(); + // Correct only the mixture interface for the interface compression flux + mixture.interfaceProperties::correct(); } if (alphaApplyPrevCorr && MULESCorr) @@ -228,19 +229,13 @@ talphaPhi1Corr0.clear(); } - #include "rhofs.H" - if ( word(mesh.schemes().ddt("ddt(rho,U)")) - == fv::EulerDdtScheme::typeName - || word(mesh.schemes().ddt("ddt(rho,U)")) - == fv::localEulerDdtScheme::typeName + != fv::EulerDdtScheme::typeName + && word(mesh.schemes().ddt("ddt(rho,U)")) + != fv::localEulerDdtScheme::typeName ) - { - rhoPhi = alphaPhi1*(rho1f - rho2f) + phiCN*rho2f; - } - else { if (ocCoeff > 0) { @@ -248,9 +243,6 @@ alphaPhi1 = (alphaPhi1 - (1.0 - cnCoeff)*alphaPhi1.oldTime())/cnCoeff; } - - // Calculate the end-of-time-step mass flux - rhoPhi = alphaPhi1*(rho1f - rho2f) + phi*rho2f; } Info<< "Phase-1 volume fraction = " diff --git a/src/twoPhaseModels/VoF/alphaEqnSubCycle.H b/src/twoPhaseModels/VoF/alphaEqnSubCycle.H deleted file mode 100644 index bdd4070f87..0000000000 --- a/src/twoPhaseModels/VoF/alphaEqnSubCycle.H +++ /dev/null @@ -1,41 +0,0 @@ -if (nAlphaSubCycles > 1) -{ - dimensionedScalar totalDeltaT = runTime.deltaT(); - surfaceScalarField rhoPhiSum - ( - IOobject - ( - "rhoPhiSum", - runTime.name(), - mesh - ), - mesh, - dimensionedScalar(rhoPhi.dimensions(), 0) - ); - - tmp trSubDeltaT; - - if (LTS) - { - trSubDeltaT = - fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles); - } - - for - ( - subCycle alphaSubCycle(alpha1, nAlphaSubCycles); - !(++alphaSubCycle).end(); - ) - { - #include "alphaEqn.H" - rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; - } - - rhoPhi = rhoPhiSum; -} -else -{ - #include "alphaEqn.H" -} - -rho == mixture.rho(); diff --git a/src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/immiscibleCompressibleTwoPhaseMixture.C b/src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/immiscibleCompressibleTwoPhaseMixture.C index d54d42cbe6..f549ab60bd 100644 --- a/src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/immiscibleCompressibleTwoPhaseMixture.C +++ b/src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/immiscibleCompressibleTwoPhaseMixture.C @@ -40,7 +40,6 @@ immiscibleCompressibleTwoPhaseMixture(const volVectorField& U) void Foam::immiscibleCompressibleTwoPhaseMixture::correct() { compressibleTwoPhaseMixture::correct(); - interfaceProperties::correct(); } diff --git a/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C b/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C index d64172c18a..ce74646565 100644 --- a/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C +++ b/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C @@ -40,7 +40,6 @@ immiscibleIncompressibleTwoPhaseMixture(const volVectorField& U) void Foam::immiscibleIncompressibleTwoPhaseMixture::correct() { incompressibleTwoPhaseMixture::correct(); - interfaceProperties::correct(); } diff --git a/tutorials/modules/CHT/VoFcoolingCylinder2D/system/fluid/fvSchemes b/tutorials/modules/CHT/VoFcoolingCylinder2D/system/fluid/fvSchemes index 8ecc19544b..83cef88b09 100644 --- a/tutorials/modules/CHT/VoFcoolingCylinder2D/system/fluid/fvSchemes +++ b/tutorials/modules/CHT/VoFcoolingCylinder2D/system/fluid/fvSchemes @@ -36,7 +36,8 @@ divSchemes div(rhoPhi,U) Gauss linearUpwind limited; turbulence Gauss limitedLinear 1; - div(rhoPhi,T) $turbulence; + div(alphaRhoPhi,e) $turbulence; + div(alphaRhoPhi,T) Gauss upwind; div(rhoPhi,K) $turbulence; div(phi,p) $turbulence; div(phi,k) $turbulence; diff --git a/tutorials/modules/compressibleVoF/ballValve/system/fvSchemes b/tutorials/modules/compressibleVoF/ballValve/system/fvSchemes index 610de57c32..1bad0a8d9d 100644 --- a/tutorials/modules/compressibleVoF/ballValve/system/fvSchemes +++ b/tutorials/modules/compressibleVoF/ballValve/system/fvSchemes @@ -33,10 +33,11 @@ divSchemes div(phi,alpha) Gauss interfaceCompression vanLeer 1; div(rhoPhi,U) Gauss limitedLinear 1; div(rhoPhi,K) Gauss limitedLinear 1; - div(rhoPhi,T) Gauss limitedLinear 1; + div(phi,p) Gauss upwind; + div(alphaRhoPhi,e) Gauss limitedLinear 1; + div(alphaRhoPhi,T) Gauss upwind; div(rhoPhi,k) Gauss upwind; div(rhoPhi,epsilon) Gauss upwind; - div(phi,p) Gauss upwind; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; } diff --git a/tutorials/modules/compressibleVoF/climbingRod/system/fvSchemes b/tutorials/modules/compressibleVoF/climbingRod/system/fvSchemes index 7ff41630b5..437b790d32 100644 --- a/tutorials/modules/compressibleVoF/climbingRod/system/fvSchemes +++ b/tutorials/modules/compressibleVoF/climbingRod/system/fvSchemes @@ -26,11 +26,13 @@ gradSchemes divSchemes { - div(rhoPhi,U) Gauss linearUpwindV grad(U); - div(rhoPhi,T) Gauss linearUpwind grad(T); - div(phi,alpha) Gauss interfaceCompression vanLeer 1; + div(rhoPhi,U) Gauss linearUpwindV grad(U); + + div(alphaRhoPhi,e) Gauss linearUpwind grad(e); + div(alphaRhoPhi,T) Gauss upwind; + div(phi,p) Gauss upwind; div(rhoPhi,K) Gauss upwind; diff --git a/tutorials/modules/compressibleVoF/cylinder/system/fvSchemes b/tutorials/modules/compressibleVoF/cylinder/system/fvSchemes index 0d7d1c0432..76856d7176 100644 --- a/tutorials/modules/compressibleVoF/cylinder/system/fvSchemes +++ b/tutorials/modules/compressibleVoF/cylinder/system/fvSchemes @@ -26,11 +26,13 @@ gradSchemes divSchemes { - div(rhoPhi,U) Gauss linearUpwindV grad(U); - div(rhoPhi,T) Gauss linearUpwind grad(T); - div(phi,alpha) Gauss interfaceCompression vanLeer 1; + div(rhoPhi,U) Gauss linearUpwindV grad(U); + + div(alphaRhoPhi,e) Gauss linearUpwind grad(e); + div(alphaRhoPhi,T) Gauss upwind; + div(phi,p) Gauss upwind; div(rhoPhi,K) Gauss upwind; diff --git a/tutorials/modules/compressibleVoF/damBreak/system/fvSchemes b/tutorials/modules/compressibleVoF/damBreak/system/fvSchemes index 9dbe1506f2..3b8ab02076 100644 --- a/tutorials/modules/compressibleVoF/damBreak/system/fvSchemes +++ b/tutorials/modules/compressibleVoF/damBreak/system/fvSchemes @@ -29,7 +29,10 @@ divSchemes div(phi,alpha) Gauss interfaceCompression vanLeer 1; div(rhoPhi,U) Gauss linearUpwind grad(U);Gauss upwind; - div(rhoPhi,T) Gauss upwind; + + div(alphaRhoPhi,e) Gauss upwind; + div(alphaRhoPhi,T) Gauss upwind; + div(rhoPhi,K) Gauss upwind; div(phi,p) Gauss upwind; diff --git a/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSchemes b/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSchemes index 27525c071d..25cf33d862 100644 --- a/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSchemes +++ b/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSchemes @@ -29,7 +29,10 @@ divSchemes div(phi,alpha) Gauss interfaceCompression vanLeer 1; div(rhoPhi,U) Gauss upwind; - div(rhoPhi,T) Gauss upwind; + + div(alphaRhoPhi,e) Gauss upwind; + div(alphaRhoPhi,T) Gauss upwind; + div(rhoPhi,K) Gauss upwind; div(phi,p) Gauss upwind; div(phi,k) Gauss upwind; diff --git a/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSolution b/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSolution index fa770e4225..99015d85b1 100644 --- a/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSolution +++ b/tutorials/modules/compressibleVoF/depthCharge2D/system/fvSolution @@ -97,8 +97,8 @@ PIMPLE { momentumPredictor no; transonic no; - nOuterCorrectors 2; - nCorrectors 2; + nOuterCorrectors 3; + nCorrectors 1; nNonOrthogonalCorrectors 0; } diff --git a/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSchemes b/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSchemes index 27525c071d..25cf33d862 100644 --- a/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSchemes +++ b/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSchemes @@ -29,7 +29,10 @@ divSchemes div(phi,alpha) Gauss interfaceCompression vanLeer 1; div(rhoPhi,U) Gauss upwind; - div(rhoPhi,T) Gauss upwind; + + div(alphaRhoPhi,e) Gauss upwind; + div(alphaRhoPhi,T) Gauss upwind; + div(rhoPhi,K) Gauss upwind; div(phi,p) Gauss upwind; div(phi,k) Gauss upwind; diff --git a/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSolution b/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSolution index 1428d681d5..82b21c5e3c 100644 --- a/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSolution +++ b/tutorials/modules/compressibleVoF/depthCharge3D/system/fvSolution @@ -98,7 +98,7 @@ PIMPLE momentumPredictor no; transonic no; nOuterCorrectors 3; - nCorrectors 2; + nCorrectors 1; nNonOrthogonalCorrectors 0; } diff --git a/tutorials/modules/compressibleVoF/plateFilm/system/fvSchemes b/tutorials/modules/compressibleVoF/plateFilm/system/fvSchemes index 0d7d1c0432..76856d7176 100644 --- a/tutorials/modules/compressibleVoF/plateFilm/system/fvSchemes +++ b/tutorials/modules/compressibleVoF/plateFilm/system/fvSchemes @@ -26,11 +26,13 @@ gradSchemes divSchemes { - div(rhoPhi,U) Gauss linearUpwindV grad(U); - div(rhoPhi,T) Gauss linearUpwind grad(T); - div(phi,alpha) Gauss interfaceCompression vanLeer 1; + div(rhoPhi,U) Gauss linearUpwindV grad(U); + + div(alphaRhoPhi,e) Gauss linearUpwind grad(e); + div(alphaRhoPhi,T) Gauss upwind; + div(phi,p) Gauss upwind; div(rhoPhi,K) Gauss upwind; diff --git a/tutorials/modules/compressibleVoF/sloshingTank2D/system/fvSchemes b/tutorials/modules/compressibleVoF/sloshingTank2D/system/fvSchemes index 2271bc4c13..50faa84e82 100644 --- a/tutorials/modules/compressibleVoF/sloshingTank2D/system/fvSchemes +++ b/tutorials/modules/compressibleVoF/sloshingTank2D/system/fvSchemes @@ -29,7 +29,10 @@ divSchemes div(phi,alpha) Gauss interfaceCompression vanLeer 1; div(rhoPhi,U) Gauss vanLeerV; - div(rhoPhi,T) Gauss vanLeer; + + div(alphaRhoPhi,e) Gauss vanLeer; + div(alphaRhoPhi,T) Gauss upwind; + div(rhoPhi,K) Gauss linear; div(phi,p) Gauss linear;