From b36d8fca44936739d90a8f05371b35c7341b19f9 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Thu, 17 Nov 2022 18:08:59 +0000 Subject: [PATCH] solvers::isothermalFluid: Improved the consistency, accuracy and stability of the pressure-work term for moving mesh cases, in both the internal energy and enthalpy equations --- .../solvers/combustion/XiFoam/EaEqn.H | 2 +- .../solvers/combustion/XiFoam/EauEqn.H | 2 +- .../solvers/combustion/XiFoam/PDRFoam/EaEqn.H | 2 +- .../combustion/XiFoam/PDRFoam/EauEqn.H | 2 +- .../modules/fluid/thermophysicalPredictor.C | 2 +- .../isothermalFluid/correctBuoyantPressure.C | 2 +- .../modules/isothermalFluid/correctPressure.C | 2 +- .../modules/isothermalFluid/isothermalFluid.C | 31 +++++++++++++++++++ .../modules/isothermalFluid/isothermalFluid.H | 11 +++++++ .../thermophysicalPredictor.C | 8 ++--- .../lockExchange/system/fvSchemes | 3 +- 11 files changed, 55 insertions(+), 12 deletions(-) diff --git a/applications/solvers/combustion/XiFoam/EaEqn.H b/applications/solvers/combustion/XiFoam/EaEqn.H index 6827211f60..3c7917ddbd 100644 --- a/applications/solvers/combustion/XiFoam/EaEqn.H +++ b/applications/solvers/combustion/XiFoam/EaEqn.H @@ -7,7 +7,7 @@ + fvc::ddt(rho, K) + fvc::div(phi, K) + ( hea.name() == "ea" - ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho) + ? fvc::div(fvc::absolute(phi, rho, U), p/rho) : -dpdt ) + thermophysicalTransport.divq(hea) diff --git a/applications/solvers/combustion/XiFoam/EauEqn.H b/applications/solvers/combustion/XiFoam/EauEqn.H index b67022e7a7..23066e356c 100644 --- a/applications/solvers/combustion/XiFoam/EauEqn.H +++ b/applications/solvers/combustion/XiFoam/EauEqn.H @@ -8,7 +8,7 @@ if (ign.ignited()) + (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou() + ( heau.name() == "eau" - ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho) + ? fvc::div(fvc::absolute(phi, rho, U), p/rho) *rho/thermo.rhou() : -dpdt*rho/thermo.rhou() ) diff --git a/applications/solvers/combustion/XiFoam/PDRFoam/EaEqn.H b/applications/solvers/combustion/XiFoam/PDRFoam/EaEqn.H index 3aae440578..92553db474 100644 --- a/applications/solvers/combustion/XiFoam/PDRFoam/EaEqn.H +++ b/applications/solvers/combustion/XiFoam/PDRFoam/EaEqn.H @@ -7,7 +7,7 @@ + betav*fvc::ddt(rho, K) + fvc::div(phi, K) + ( hea.name() == "ea" - ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho) + ? fvc::div(fvc::absolute(phi, rho, U), p/rho) : -betav*dpdt ) - fvm::laplacian(Db, hea) diff --git a/applications/solvers/combustion/XiFoam/PDRFoam/EauEqn.H b/applications/solvers/combustion/XiFoam/PDRFoam/EauEqn.H index a8f988b9f7..1ccef6b657 100644 --- a/applications/solvers/combustion/XiFoam/PDRFoam/EauEqn.H +++ b/applications/solvers/combustion/XiFoam/PDRFoam/EauEqn.H @@ -8,7 +8,7 @@ if (ign.ignited()) + (betav*fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou() + ( heau.name() == "eau" - ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho) + ? fvc::div(fvc::absolute(phi, rho, U), p/rho) *rho/thermo.rhou() : -betav*dpdt*rho/thermo.rhou() ) diff --git a/applications/solvers/modules/fluid/thermophysicalPredictor.C b/applications/solvers/modules/fluid/thermophysicalPredictor.C index c558660a0b..489aafa0c0 100644 --- a/applications/solvers/modules/fluid/thermophysicalPredictor.C +++ b/applications/solvers/modules/fluid/thermophysicalPredictor.C @@ -37,7 +37,7 @@ void Foam::solvers::fluid::thermophysicalPredictor() + fvc::ddt(rho, K) + fvc::div(phi, K) + ( he.name() == "e" - ? fvc::div(fvc::absolute(phi, rho, U), p/rho) + ? addMotionWork(fvc::div(phi, p/rho)) : -dpdt ) + thermophysicalTransport->divq(he) diff --git a/applications/solvers/modules/isothermalFluid/correctBuoyantPressure.C b/applications/solvers/modules/isothermalFluid/correctBuoyantPressure.C index bdb8dfa504..7b0f57240c 100644 --- a/applications/solvers/modules/isothermalFluid/correctBuoyantPressure.C +++ b/applications/solvers/modules/isothermalFluid/correctBuoyantPressure.C @@ -265,7 +265,7 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() if (mesh.moving()) { - dpdt -= fvc::div(fvc::meshPhi(rho, U), p); + dpdt -= motionWork(); } } } diff --git a/applications/solvers/modules/isothermalFluid/correctPressure.C b/applications/solvers/modules/isothermalFluid/correctPressure.C index 792f6b009d..5a2dab0319 100644 --- a/applications/solvers/modules/isothermalFluid/correctPressure.C +++ b/applications/solvers/modules/isothermalFluid/correctPressure.C @@ -241,7 +241,7 @@ void Foam::solvers::isothermalFluid::correctPressure() if (mesh.moving()) { - dpdt -= fvc::div(fvc::meshPhi(rho, U), p); + dpdt -= motionWork(); } } } diff --git a/applications/solvers/modules/isothermalFluid/isothermalFluid.C b/applications/solvers/modules/isothermalFluid/isothermalFluid.C index 93dfa6518f..2f132eb4fc 100644 --- a/applications/solvers/modules/isothermalFluid/isothermalFluid.C +++ b/applications/solvers/modules/isothermalFluid/isothermalFluid.C @@ -54,6 +54,37 @@ void Foam::solvers::isothermalFluid::continuityErrors() } +// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // + +Foam::tmp +Foam::solvers::isothermalFluid::motionWork() const +{ + return fvc::div + ( + fvc::interpolate(rho)*fvc::meshPhi(rho, U), + p/rho, + "div(phi,(p|rho))" + ); +} + + +Foam::tmp +Foam::solvers::isothermalFluid::addMotionWork +( + const tmp& work +) const +{ + if (mesh.moving()) + { + return work + motionWork(); + } + else + { + return work; + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::solvers::isothermalFluid::isothermalFluid diff --git a/applications/solvers/modules/isothermalFluid/isothermalFluid.H b/applications/solvers/modules/isothermalFluid/isothermalFluid.H index eff45cb6c1..6e211c6aba 100644 --- a/applications/solvers/modules/isothermalFluid/isothermalFluid.H +++ b/applications/solvers/modules/isothermalFluid/isothermalFluid.H @@ -193,6 +193,17 @@ private: void correctBuoyantPressure(); +protected: + + //- Return the mesh-motion work + // Used to correct dpdt of the enthalpy equation + // or added to the div(Up) term of the internal energy equation + tmp motionWork() const; + + //- Adds the mesh-motion work to the work term provided + tmp addMotionWork(const tmp&) const; + + public: //- Runtime type information diff --git a/applications/solvers/modules/multicomponentFluid/thermophysicalPredictor.C b/applications/solvers/modules/multicomponentFluid/thermophysicalPredictor.C index 778ef2173c..e6c613d6b3 100644 --- a/applications/solvers/modules/multicomponentFluid/thermophysicalPredictor.C +++ b/applications/solvers/modules/multicomponentFluid/thermophysicalPredictor.C @@ -51,9 +51,9 @@ void Foam::solvers::multicomponentFluid::thermophysicalPredictor() fvScalarMatrix YiEqn ( fvm::ddt(rho, Yi) - + mvConvection->fvmDiv(phi, Yi) - + thermophysicalTransport->divj(Yi) - == + + mvConvection->fvmDiv(phi, Yi) + + thermophysicalTransport->divj(Yi) + == reaction->R(Yi) + fvModels().source(rho, Yi) ); @@ -79,7 +79,7 @@ void Foam::solvers::multicomponentFluid::thermophysicalPredictor() + fvc::ddt(rho, K) + fvc::div(phi, K) + ( he.name() == "e" - ? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho) + ? addMotionWork(mvConvection->fvcDiv(phi, p/rho)) : -dpdt ) + thermophysicalTransport->divq(he) diff --git a/tutorials/modules/multicomponentFluid/lockExchange/system/fvSchemes b/tutorials/modules/multicomponentFluid/lockExchange/system/fvSchemes index d023d0f44d..78201e10a4 100644 --- a/tutorials/modules/multicomponentFluid/lockExchange/system/fvSchemes +++ b/tutorials/modules/multicomponentFluid/lockExchange/system/fvSchemes @@ -35,7 +35,8 @@ divSchemes sludge vanLeer; ".*" upwind; }; - div(phi,K) Gauss limitedLinear 1; + div(phi,K) Gauss upwind; + div(phi,(p|rho)) Gauss upwind; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; }