From 7b825d0817e14e1dab4c75a4b33d200b7cd77e41 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Fri, 12 May 2017 11:04:38 +0100 Subject: [PATCH] rhoPimpleFoam: Updated transonic option to be consistent with sonicFoam Improves stability on start-up and allows running at slightly larger time-steps. --- .../solvers/compressible/rhoPimpleFoam/pEqn.H | 15 +- .../compressible/rhoPimpleFoam/pcEqn.H | 14 +- .../solvers/compressible/rhoSimpleFoam/pEqn.H | 197 +++++++++--------- .../compressible/rhoSimpleFoam/pcEqn.H | 14 +- .../general/pressureControl/pressureControl.C | 22 +- .../general/pressureControl/pressureControl.H | 4 +- 6 files changed, 145 insertions(+), 121 deletions(-) diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H index 5aab536239..9eecd01d99 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H @@ -1,3 +1,5 @@ +rho = thermo.rho(); + volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); @@ -10,7 +12,7 @@ if (pimple.nCorrPISO() <= 1) surfaceScalarField phiHbyA ( "phiHbyA", - fvc::flux(rho*HbyA) + fvc::interpolate(rho)*fvc::flux(HbyA) + rhorAUf*fvc::ddtCorr(rho, U, phi) ); @@ -26,7 +28,8 @@ if (pimple.transonic()) "phid", (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA ); - phiHbyA -= fvc::interpolate(p)*phid; + + phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho); while (pimple.correctNonOrthogonal()) { @@ -84,9 +87,11 @@ U.correctBoundaryConditions(); fvOptions.correct(U); K = 0.5*magSqr(U); -pressureControl.limit(p); -p.correctBoundaryConditions(); -rho = thermo.rho(); +if (pressureControl.limit(p)) +{ + p.correctBoundaryConditions(); + rho = thermo.rho(); +} if (!pimple.transonic()) { diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H index afbc2851e4..34657a1892 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H @@ -1,3 +1,5 @@ +rho = thermo.rho(); + volScalarField rAU(1.0/UEqn.A()); volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1())); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); @@ -11,7 +13,7 @@ surfaceScalarField phiHbyA ( "phiHbyA", ( - fvc::flux(rho*HbyA) + fvc::interpolate(rho)*fvc::flux(HbyA) + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi) ) ); @@ -33,7 +35,7 @@ if (pimple.transonic()) phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf() - - fvc::interpolate(p)*phid; + - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho); HbyA -= (rAU - rAtU)*fvc::grad(p); @@ -96,9 +98,11 @@ U.correctBoundaryConditions(); fvOptions.correct(U); K = 0.5*magSqr(U); -pressureControl.limit(p); -p.correctBoundaryConditions(); -rho = thermo.rho(); +if (pressureControl.limit(p)) +{ + p.correctBoundaryConditions(); + rho = thermo.rho(); +} if (!pimple.transonic()) { diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index c7edbb9e87..0f290ed273 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -1,109 +1,110 @@ +volScalarField rAU(1.0/UEqn.A()); +surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); +volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); +tUEqn.clear(); + +bool closedVolume = false; + +surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA)); +MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + +// Update the pressure BCs to ensure flux consistency +constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); + +if (simple.transonic()) { - volScalarField rAU(1.0/UEqn.A()); - surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); - volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); - tUEqn.clear(); + surfaceScalarField phid + ( + "phid", + (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA + ); - bool closedVolume = false; + phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho); - surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA)); - MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - - // Update the pressure BCs to ensure flux consistency - constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); - - if (simple.transonic()) + while (simple.correctNonOrthogonal()) { - surfaceScalarField phid + fvScalarMatrix pEqn ( - "phid", - (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA + fvc::div(phiHbyA) + + fvm::div(phid, p) + - fvm::laplacian(rhorAUf, p) + == + fvOptions(psi, p, rho.name()) ); - phiHbyA -= fvc::interpolate(p)*phid; - while (simple.correctNonOrthogonal()) + // Relax the pressure equation to ensure diagonal-dominance + pEqn.relax(); + + pEqn.setReference + ( + pressureControl.refCell(), + pressureControl.refValue() + ); + + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) { - fvScalarMatrix pEqn - ( - fvc::div(phiHbyA) - + fvm::div(phid, p) - - fvm::laplacian(rhorAUf, p) - == - fvOptions(psi, p, rho.name()) - ); - - // Relax the pressure equation to ensure diagonal-dominance - pEqn.relax(); - - pEqn.setReference - ( - pressureControl.refCell(), - pressureControl.refValue() - ); - - pEqn.solve(); - - if (simple.finalNonOrthogonalIter()) - { - phi = phiHbyA + pEqn.flux(); - } + phi = phiHbyA + pEqn.flux(); } } - else - { - closedVolume = adjustPhi(phiHbyA, U, p); - - while (simple.correctNonOrthogonal()) - { - fvScalarMatrix pEqn - ( - fvc::div(phiHbyA) - - fvm::laplacian(rhorAUf, p) - == - fvOptions(psi, p, rho.name()) - ); - - pEqn.setReference - ( - pressureControl.refCell(), - pressureControl.refValue() - ); - - pEqn.solve(); - - if (simple.finalNonOrthogonalIter()) - { - phi = phiHbyA + pEqn.flux(); - } - } - } - - - #include "incompressible/continuityErrs.H" - - // Explicitly relax pressure for momentum corrector - p.relax(); - - U = HbyA - rAU*fvc::grad(p); - U.correctBoundaryConditions(); - fvOptions.correct(U); - - pressureControl.limit(p); - - // For closed-volume cases adjust the pressure and density levels - // to obey overall mass continuity - if (closedVolume) - { - p += (initialMass - fvc::domainIntegrate(psi*p)) - /fvc::domainIntegrate(psi); - } - - p.correctBoundaryConditions(); - - rho = thermo.rho(); - - if (!simple.transonic()) - { - rho.relax(); - } +} +else +{ + closedVolume = adjustPhi(phiHbyA, U, p); + + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvc::div(phiHbyA) + - fvm::laplacian(rhorAUf, p) + == + fvOptions(psi, p, rho.name()) + ); + + pEqn.setReference + ( + pressureControl.refCell(), + pressureControl.refValue() + ); + + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } +} + +#include "incompressible/continuityErrs.H" + +// Explicitly relax pressure for momentum corrector +p.relax(); + +U = HbyA - rAU*fvc::grad(p); +U.correctBoundaryConditions(); +fvOptions.correct(U); + +bool pLimited = pressureControl.limit(p); + +// For closed-volume cases adjust the pressure and density levels +// to obey overall mass continuity +if (closedVolume) +{ + p += (initialMass - fvc::domainIntegrate(psi*p)) + /fvc::domainIntegrate(psi); +} + +if (pLimited || closedVolume) +{ + p.correctBoundaryConditions(); +} + +rho = thermo.rho(); + +if (!simple.transonic()) +{ + rho.relax(); } diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H index c7dc0f864d..51560aab5f 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H @@ -1,3 +1,5 @@ +rho = thermo.rho(); + volScalarField rAU(1.0/UEqn.A()); volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1())); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); @@ -5,7 +7,7 @@ tUEqn.clear(); bool closedVolume = false; -surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA)); +surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA)); MRF.makeRelative(fvc::interpolate(rho), phiHbyA); volScalarField rhorAtU("rhorAtU", rho*rAtU); @@ -23,7 +25,7 @@ if (simple.transonic()) phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf() - - fvc::interpolate(p)*phid; + - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho); HbyA -= (rAU - rAtU)*fvc::grad(p); @@ -98,7 +100,7 @@ U = HbyA - rAtU*fvc::grad(p); U.correctBoundaryConditions(); fvOptions.correct(U); -pressureControl.limit(p); +bool pLimited = pressureControl.limit(p); // For closed-volume cases adjust the pressure and density levels // to obey overall mass continuity @@ -108,9 +110,11 @@ if (closedVolume) /fvc::domainIntegrate(psi); } -p.correctBoundaryConditions(); +if (pLimited || closedVolume) +{ + p.correctBoundaryConditions(); +} -// Recalculate density from the relaxed pressure rho = thermo.rho(); if (!simple.transonic()) diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C index aca22e7699..abd71e60ac 100644 --- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C +++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C @@ -207,14 +207,24 @@ Foam::pressureControl::pressureControl // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -void Foam::pressureControl::limit(volScalarField& p) const +bool Foam::pressureControl::limit(volScalarField& p) const { - Info<< "pressureControl: p max/min " - << max(p).value() << " " - << min(p).value() << endl; + scalar pMax = max(p).value(); + scalar pMin = min(p).value(); - p = max(p, pMin_); - p = min(p, pMax_); + if (pMax > pMax_.value() || pMin < pMin_.value()) + { + Info<< "pressureControl: p max/min " << pMax << " " << pMin << endl; + + p = max(p, pMin_); + p = min(p, pMax_); + + return true; + } + else + { + return false; + } } diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H index 61f510ac60..0bb49c99a5 100644 --- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H +++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H @@ -89,8 +89,8 @@ public: //- Return the pressure reference level inline scalar refValue() const; - //- Limit the pressure - void limit(volScalarField& p) const; + //- Limit the pressure if necessary and return true if so + bool limit(volScalarField& p) const; };