From c22a5a7aa6b33ab9903a7236cb1f4a396fdd4754 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Thu, 22 Sep 2022 15:05:33 +0100 Subject: [PATCH] isothermalFluid::correctBuoyantPressure: Relax the net force momentum equation source rather than p_rgh which introduces an imbalance between the pressure and buoyancy forces. The relaxation factor for p_rgh specified in fvSolution is used to relax the net force as the intent is to relax the pressure and this provides convenient usage and backwards-compatibility. Optional relaxation for the thermodynamic pressure p is also available for case this provides convergence benefit for steady cases by relaxing the pressure work term is the energy equation. Note these changes only relate to the operation of the isothermalFluid solver module for buoyant cases. --- .../isothermalFluid/correctBuoyantPressure.C | 36 ++++++---- .../fluid/isothermalFluid/isothermalFluid.C | 15 ++++ .../fluid/isothermalFluid/isothermalFluid.H | 4 ++ .../fluid/isothermalFluid/momentumPredictor.C | 8 +-- .../GeometricField/GeometricField.C | 72 +++++++++++++++---- .../GeometricField/GeometricField.H | 26 +++++-- 6 files changed, 123 insertions(+), 38 deletions(-) diff --git a/applications/solvers/modules/fluid/isothermalFluid/correctBuoyantPressure.C b/applications/solvers/modules/fluid/isothermalFluid/correctBuoyantPressure.C index 291831329e..cf0bd861ed 100644 --- a/applications/solvers/modules/fluid/isothermalFluid/correctBuoyantPressure.C +++ b/applications/solvers/modules/fluid/isothermalFluid/correctBuoyantPressure.C @@ -102,8 +102,9 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() *fvc::relative(phiHbyA, rho, U) ); - const fvScalarMatrix divPhidp(fvm::div(phid, p)); - phiHbyA -= divPhidp.flux(); + // Subtract the compressible part + // The resulting flux will be zero for a perfect gas + phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho); if (pimple.consistent()) { @@ -120,7 +121,7 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() fvScalarMatrix p_rghDDtEqn ( fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) - + fvc::div(phiHbyA) + divPhidp + + fvc::div(phiHbyA) + fvm::div(phid, p) == fvModels().source(psi, p_rgh, rho.name()) ); @@ -183,6 +184,21 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() phi = phiHbyA + p_rghEqn.flux(); + // Calculate and relax the net pressure-buoyancy force + netForce.ref().relax + ( + fvc::reconstruct((ghGradRhof + p_rghEqn.flux()/rhorAAtUf)), + p_rgh.relaxationFactor() + ); + + // Correct the momentum source with the pressure gradient flux + // calculated from the relaxed pressure + U = HbyA + rAAtU*netForce(); + U.correctBoundaryConditions(); + fvConstraints().constrain(U); + + K = 0.5*magSqr(U); + if (!mesh.schemes().steady()) { p = p_rgh + rho*gh + pRef; @@ -202,19 +218,8 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() continuityErrors(); - // Explicitly relax pressure for momentum corrector - p_rgh.relax(); - p = p_rgh + rho*gh + pRef; - // Correct the momentum source with the pressure gradient flux - // calculated from the relaxed pressure - U = HbyA + rAAtU*fvc::reconstruct((ghGradRhof + p_rghEqn.flux()/rhorAAtUf)); - U.correctBoundaryConditions(); - fvConstraints().constrain(U); - - K = 0.5*magSqr(U); - if (mesh.schemes().steady()) { if (fvConstraints().constrain(p)) @@ -234,6 +239,9 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() p_rgh.correctBoundaryConditions(); } + // Optionally relax pressure for the thermophysics + p.relax(); + if (mesh.schemes().steady() || pimple.simpleRho() || adjustMass) { rho = thermo.rho(); diff --git a/applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.C b/applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.C index 3ccbff8fd5..6da3efdead 100644 --- a/applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.C +++ b/applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.C @@ -170,6 +170,21 @@ Foam::solvers::isothermalFluid::isothermalFluid thermo, pimple.dict() ); + + netForce = new volVectorField + ( + IOobject + ( + "netForce", + runTime.timeName(), + mesh + ), + fvc::reconstruct + ( + (-buoyancy->ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh)) + *mesh.magSf() + ) + ); } if (mesh.dynamic()) diff --git a/applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.H b/applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.H index a20550be75..f11aeaa1ee 100644 --- a/applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.H +++ b/applications/solvers/modules/fluid/isothermalFluid/isothermalFluid.H @@ -147,6 +147,10 @@ protected: // Cached temporary fields + //- Momentum equation net force source term + // Used for buoyant simulations only + tmp netForce; + //- Pointer to the surface momentum field // used to recreate the flux after mesh-change autoPtr rhoUf; diff --git a/applications/solvers/modules/fluid/isothermalFluid/momentumPredictor.C b/applications/solvers/modules/fluid/isothermalFluid/momentumPredictor.C index e23eca03d0..75452a97b0 100644 --- a/applications/solvers/modules/fluid/isothermalFluid/momentumPredictor.C +++ b/applications/solvers/modules/fluid/isothermalFluid/momentumPredictor.C @@ -51,13 +51,7 @@ void Foam::solvers::isothermalFluid::momentumPredictor() ( UEqn == - fvc::reconstruct - ( - ( - - buoyancy->ghf*fvc::snGrad(rho) - - fvc::snGrad(p_rgh) - )*mesh.magSf() - ) + netForce() ); } else diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C index 4281adf868..036475264e 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C @@ -1274,18 +1274,23 @@ bool Foam::GeometricField::needReference() const template class PatchField, class GeoMesh> void Foam::GeometricField::relax(const scalar alpha) { - if (debug) + if (alpha < 1) { - InfoInFunction - << "Relaxing" << endl << this->info() << " by " << alpha << endl; - } + if (debug) + { + InfoInFunction + << "Relaxing" << endl << this->info() + << " by " << alpha << endl; + } - operator==(prevIter() + alpha*(*this - prevIter())); + operator==(prevIter() + alpha*(*this - prevIter())); + } } template class PatchField, class GeoMesh> -void Foam::GeometricField::relax() +Foam::scalar +Foam::GeometricField::relaxationFactor() const { if ( @@ -1297,18 +1302,61 @@ void Foam::GeometricField::relax() && this->mesh().solution().relaxField(this->name() + "Final") ) { - relax + return this->mesh().solution().fieldRelaxationFactor ( - this->mesh().solution().fieldRelaxationFactor - ( - this->name() + "Final" - ) + this->name() + "Final" ); } else if (this->mesh().solution().relaxField(this->name())) { - relax(this->mesh().solution().fieldRelaxationFactor(this->name())); + return this->mesh().solution().fieldRelaxationFactor(this->name()); } + else + { + return 1; + } +} + + +template class PatchField, class GeoMesh> +void Foam::GeometricField::relax() +{ + relax(relaxationFactor()); +} + + +template class PatchField, class GeoMesh> +void Foam::GeometricField::relax +( + const tmp>& tgf, + const scalar alpha +) +{ + if (alpha < 1) + { + if (debug) + { + InfoInFunction + << "Relaxing" << endl << this->info() + << " by " << alpha << endl; + } + + operator==(*this + alpha*(tgf - *this)); + } + else + { + operator==(tgf); + } +} + + +template class PatchField, class GeoMesh> +void Foam::GeometricField::relax +( + const tmp>& tgf +) +{ + relax(tgf, relaxationFactor()); } diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H index 932d26615a..a2c00276a4 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H @@ -462,15 +462,31 @@ public: tmp> T() const; //- Relax field (for steady-state solution). - // alpha = 1 : no relaxation - // alpha < 1 : relaxation - // alpha = 0 : do nothing + // alpha >= 1 : no relaxation + // alpha < 1 : relaxation void relax(const scalar alpha); - //- Relax field (for steady-state solution). - // alpha is read from controlDict + //- Return the field relaxation factor read from fvSolution + // or 1 if not specified + scalar relaxationFactor() const; + + //- Relax current field with respect to the cached previous iteration. + // Relaxation factor is read from fvSolution void relax(); + //- Relax given field with respect to the current field + // and reset the field to the result + void relax + ( + const tmp>&, + const scalar alpha + ); + + //- Relax given field with respect to the current field + // and reset the field to the result + // Relaxation factor is read from fvSolution + void relax(const tmp>&); + //- Select the final iteration parameters if `final' is true // by returning the field name + "Final" // otherwise the standard parameters by returning the field name