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