From 99c992d65c02cbd58e3927f8777d9eb11cceefcf Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Tue, 28 Feb 2017 11:14:59 +0000 Subject: [PATCH] rhoPimpleFoam: Added support for transonic flow of liquids and real gases Both stardard SIMPLE and the SIMPLEC (using the 'consistent' option in fvSolution) are now supported for both subsonic and transonic flow of all fluid types. rhoPimpleFoam now instantiates the lower-level fluidThermo which instantiates either a psiThermo or rhoThermo according to the 'type' specification in thermophysicalProperties, see also commit 655fc7874808927d14916307a2230a8965bdb860 --- .../compressible/rhoPimpleFoam/createFields.H | 28 +-- .../solvers/compressible/rhoPimpleFoam/pEqn.H | 65 +++--- .../compressible/rhoPimpleFoam/pcEqn.H | 76 +++--- .../rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H | 64 +++-- .../rhoPimpleDyMFoam/rhoPimpleDyMFoam.C | 10 +- .../rhoPimpleFoam/rhoPimpleFoam.C | 5 +- .../compressible/rhoSimpleFoam/createFields.H | 4 +- .../solvers/compressible/rhoSimpleFoam/pEqn.H | 18 +- .../compressible/rhoSimpleFoam/pcEqn.H | 28 +-- .../rhoPorousSimpleFoam/createFields.H | 59 ----- .../general/pressureControl/pressureControl.C | 220 ++++++++++++------ .../general/pressureControl/pressureControl.H | 3 +- .../annularThermalMixer/system/fvSolution | 5 +- .../LES/pitzDaily/system/fvSolution | 5 +- .../constant/thermophysicalProperties | 2 +- .../RAS/angledDuct/system/fvSchemes | 3 +- .../RAS/angledDuct/system/fvSolution | 14 +- .../RAS/angledDuctLTS/system/fvSolution | 4 +- .../RAS/cavity/system/fvSolution | 5 +- .../RAS/mixerVessel2D/system/fvSolution | 15 +- .../helmholtzResonance/system/fvSolution | 5 +- 21 files changed, 295 insertions(+), 343 deletions(-) delete mode 100644 applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H index a7ee3eca4..84f70f2cd 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H @@ -2,11 +2,11 @@ Info<< "Reading thermophysical properties\n" << endl; -autoPtr pThermo +autoPtr pThermo ( - psiThermo::New(mesh) + fluidThermo::New(mesh) ); -psiThermo& thermo = pThermo(); +fluidThermo& thermo = pThermo(); thermo.validate(args.executable(), "h", "e"); volScalarField& p = thermo.p(); @@ -40,27 +40,7 @@ volVectorField U #include "compressibleCreatePhi.H" -dimensionedScalar rhoMax -( - dimensionedScalar::lookupOrDefault - ( - "rhoMax", - pimple.dict(), - dimDensity, - GREAT - ) -); - -dimensionedScalar rhoMin -( - dimensionedScalar::lookupOrDefault - ( - "rhoMin", - pimple.dict(), - dimDensity, - 0 - ) -); +pressureControl pressureControl(p, rho, pimple.dict(), false); Info<< "Creating turbulence model\n" << endl; autoPtr turbulence diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H index ac7107acf..73ecafd89 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H @@ -1,8 +1,3 @@ -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); -rho.relax(); - volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); @@ -12,55 +7,54 @@ if (pimple.nCorrPISO() <= 1) tUEqn.clear(); } +surfaceScalarField phiHbyA +( + "phiHbyA", + ( + fvc::flux(rho*HbyA) + + rhorAUf*fvc::ddtCorr(rho, U, phi) + ) +); + +MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + +// Update the pressure BCs to ensure flux consistency +constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); + if (pimple.transonic()) { surfaceScalarField phid ( "phid", - fvc::interpolate(psi) - *( - fvc::flux(HbyA) - + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho) - ) + (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA ); - - MRF.makeRelative(fvc::interpolate(psi), phid); + phiHbyA -= fvc::interpolate(p)*phid; while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + + 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.solve(mesh.solver(p.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { - phi == pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } else { - surfaceScalarField phiHbyA - ( - "phiHbyA", - ( - fvc::flux(rho*HbyA) - + rhorAUf*fvc::ddtCorr(rho, U, phi) - ) - ); - - MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - - // Update the pressure BCs to ensure flux consistency - constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn @@ -87,19 +81,20 @@ else // Explicitly relax pressure for momentum corrector p.relax(); -// Recalculate density from the relaxed pressure -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); -rho.relax(); -Info<< "rho max/min : " << max(rho).value() - << " " << min(rho).value() << endl; - U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); fvOptions.correct(U); K = 0.5*magSqr(U); +pressureControl.limit(p); +p.correctBoundaryConditions(); +rho = thermo.rho(); + +if (!pimple.transonic()) +{ + rho.relax(); +} + if (thermo.dpdt()) { dpdt = fvc::ddt(p); diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H index 713f443fc..afbc2851e 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H @@ -1,8 +1,3 @@ -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); -rho.relax(); - volScalarField rAU(1.0/UEqn.A()); volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1())); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); @@ -12,72 +7,64 @@ if (pimple.nCorrPISO() <= 1) tUEqn.clear(); } +surfaceScalarField phiHbyA +( + "phiHbyA", + ( + fvc::flux(rho*HbyA) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi) + ) +); + +MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + +volScalarField rhorAtU("rhorAtU", rho*rAtU); + +// Update the pressure BCs to ensure flux consistency +constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF); + if (pimple.transonic()) { surfaceScalarField phid ( "phid", - fvc::interpolate(psi) - *( - fvc::flux(HbyA) - + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi) - /fvc::interpolate(rho) - ) + (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA ); - MRF.makeRelative(fvc::interpolate(psi), phid); - - surfaceScalarField phic - ( - "phic", + phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf() - ); + - fvc::interpolate(p)*phid; HbyA -= (rAU - rAtU)*fvc::grad(p); - volScalarField rhorAtU("rhorAtU", rho*rAtU); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + + fvc::div(phiHbyA) + fvm::div(phid, p) - + fvc::div(phic) - fvm::laplacian(rhorAtU, p) == fvOptions(psi, p, rho.name()) ); + // Relax the pressure equation to ensure diagonal-dominance + pEqn.relax(); + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { - phi == phic + pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } else { - surfaceScalarField phiHbyA - ( - "phiHbyA", - ( - fvc::flux(rho*HbyA) - + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi) - ) - ); - - MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf(); HbyA -= (rAU - rAtU)*fvc::grad(p); - volScalarField rhorAtU("rhorAtU", rho*rAtU); - - // Update the pressure BCs to ensure flux consistency - constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn @@ -109,19 +96,16 @@ U.correctBoundaryConditions(); fvOptions.correct(U); K = 0.5*magSqr(U); -if (thermo.dpdt()) -{ - dpdt = fvc::ddt(p); -} - -// Recalculate density from the relaxed pressure +pressureControl.limit(p); +p.correctBoundaryConditions(); rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); if (!pimple.transonic()) { rho.relax(); } -Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; +if (thermo.dpdt()) +{ + dpdt = fvc::ddt(p); +} diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H index 7bd540df4..0f8556261 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H @@ -1,8 +1,3 @@ -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); -rho.relax(); - volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); @@ -12,55 +7,53 @@ if (pimple.nCorrPISO() <= 1) tUEqn.clear(); } +surfaceScalarField phiHbyA +( + "phiHbyA", + fvc::flux(rho*HbyA) + + rhorAUf*fvc::ddtCorr(rho, U, rhoUf) +); + +fvc::makeRelative(phiHbyA, rho, U); +MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + +// Update the pressure BCs to ensure flux consistency +constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); + if (pimple.transonic()) { surfaceScalarField phid ( "phid", - fvc::interpolate(psi) - *( - fvc::flux(HbyA) - + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho) - ) + (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA ); - - fvc::makeRelative(phid, psi, U); - MRF.makeRelative(fvc::interpolate(psi), phid); + phiHbyA -= fvc::interpolate(p)*phid; while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + + 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.solve(mesh.solver(p.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { - phi == pEqn.flux(); + phi = phiHbyA + pEqn.flux(); } } } else { - surfaceScalarField phiHbyA - ( - "phiHbyA", - fvc::flux(rho*HbyA) - + rhorAUf*fvc::ddtCorr(rho, U, rhoUf) - ); - - fvc::makeRelative(phiHbyA, rho, U); - MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - - // Update the pressure BCs to ensure flux consistency - constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); - while (pimple.correctNonOrthogonal()) { // Pressure corrector @@ -88,19 +81,20 @@ else // Explicitly relax pressure for momentum corrector p.relax(); -// Recalculate density from the relaxed pressure -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); -rho.relax(); -Info<< "rho max/min : " << max(rho).value() - << " " << min(rho).value() << endl; - U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); fvOptions.correct(U); K = 0.5*magSqr(U); +pressureControl.limit(p); +p.correctBoundaryConditions(); +rho = thermo.rho(); + +if (!pimple.transonic()) +{ + rho.relax(); +} + { rhoUf = fvc::interpolate(rho*U); surfaceVectorField n(mesh.Sf()/mesh.magSf()); diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C index c7046bb9b..1d20c3f10 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -22,10 +22,7 @@ License along with OpenFOAM. If not, see . Application - rhoPimpleFoam - -Group - grpCompressibleSolvers grpMovingMeshSolvers + rhoPimpleDyMFoam Description Transient solver for turbulent flow of compressible fluids for HVAC and @@ -38,10 +35,11 @@ Description #include "fvCFD.H" #include "dynamicFvMesh.H" -#include "psiThermo.H" +#include "fluidThermo.H" #include "turbulentFluidThermoModel.H" #include "bound.H" #include "pimpleControl.H" +#include "pressureControl.H" #include "CorrectPhi.H" #include "fvOptions.H" #include "localEulerDdtScheme.H" diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C index 108ad9aa0..f30dae76a 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -34,10 +34,11 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "psiThermo.H" +#include "fluidThermo.H" #include "turbulentFluidThermoModel.H" #include "bound.H" #include "pimpleControl.H" +#include "pressureControl.H" #include "fvOptions.H" #include "localEulerDdtScheme.H" #include "fvcSmooth.H" diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H index 4671347b6..c7914f89d 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H @@ -7,6 +7,8 @@ autoPtr pThermo fluidThermo& thermo = pThermo(); thermo.validate(args.executable(), "h", "e"); +volScalarField& p = thermo.p(); + volScalarField rho ( IOobject @@ -20,8 +22,6 @@ volScalarField rho thermo.rho() ); -volScalarField& p = thermo.p(); - Info<< "Reading field U\n" << endl; volVectorField U ( diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index 5f092ada2..c7edbb9e8 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -6,16 +6,18 @@ bool closedVolume = false; + 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()) { - surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA)); - surfaceScalarField rhof(fvc::interpolate(rho)); - MRF.makeRelative(rhof, phiHbyA); - surfaceScalarField phid ( "phid", - (fvc::interpolate(psi)/rhof)*phiHbyA + (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA ); phiHbyA -= fvc::interpolate(p)*phid; @@ -49,14 +51,8 @@ } else { - surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA)); - MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - closedVolume = adjustPhi(phiHbyA, U, p); - // Update the pressure BCs to ensure flux consistency - constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); - while (simple.correctNonOrthogonal()) { fvScalarMatrix pEqn diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H index af2005feb..c7dc0f864 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H @@ -5,16 +5,20 @@ tUEqn.clear(); bool closedVolume = false; +surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA)); +MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + +volScalarField rhorAtU("rhorAtU", rho*rAtU); + +// Update the pressure BCs to ensure flux consistency +constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF); + if (simple.transonic()) { - surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA)); - surfaceScalarField rhof(fvc::interpolate(rho)); - MRF.makeRelative(rhof, phiHbyA); - surfaceScalarField phid ( "phid", - (fvc::interpolate(psi)/rhof)*phiHbyA + (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA ); phiHbyA += @@ -23,14 +27,12 @@ if (simple.transonic()) HbyA -= (rAU - rAtU)*fvc::grad(p); - volScalarField rhorAtU("rhorAtU", rho*rAtU); - while (simple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( - fvm::div(phid, p) - + fvc::div(phiHbyA) + fvc::div(phiHbyA) + + fvm::div(phid, p) - fvm::laplacian(rhorAtU, p) == fvOptions(psi, p, rho.name()) @@ -55,19 +57,11 @@ if (simple.transonic()) } else { - surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA)); - MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - closedVolume = adjustPhi(phiHbyA, U, p); phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf(); HbyA -= (rAU - rAtU)*fvc::grad(p); - volScalarField rhorAtU("rhorAtU", rho*rAtU); - - // Update the pressure BCs to ensure flux consistency - constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF); - while (simple.correctNonOrthogonal()) { fvScalarMatrix pEqn diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H deleted file mode 100644 index 4671347b6..000000000 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H +++ /dev/null @@ -1,59 +0,0 @@ -Info<< "Reading thermophysical properties\n" << endl; - -autoPtr pThermo -( - fluidThermo::New(mesh) -); -fluidThermo& thermo = pThermo(); -thermo.validate(args.executable(), "h", "e"); - -volScalarField rho -( - IOobject - ( - "rho", - runTime.timeName(), - mesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - thermo.rho() -); - -volScalarField& p = thermo.p(); - -Info<< "Reading field U\n" << endl; -volVectorField U -( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -#include "compressibleCreatePhi.H" - -pressureControl pressureControl(p, rho, simple.dict()); - -mesh.setFluxRequired(p.name()); - -Info<< "Creating turbulence model\n" << endl; -autoPtr turbulence -( - compressible::turbulenceModel::New - ( - rho, - U, - phi, - thermo - ) -); - -dimensionedScalar initialMass = fvc::domainIntegrate(rho); - -#include "createMRF.H" diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C index 101126a5e..aca22e769 100644 --- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C +++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C @@ -32,105 +32,171 @@ Foam::pressureControl::pressureControl ( const volScalarField& p, const volScalarField& rho, - const dictionary& dict + const dictionary& dict, + const bool pRefRequired ) : - refCell_(0), + refCell_(-1), refValue_(0), pMax_("pMax", dimPressure, 0), pMin_("pMin", dimPressure, GREAT) { - if (setRefCell(p, dict, refCell_, refValue_)) + bool pLimits = false; + + // Set the reference cell and value for closed domain simulations + if (pRefRequired && setRefCell(p, dict, refCell_, refValue_)) { + pLimits = true; + pMax_.value() = refValue_; pMin_.value() = refValue_; } - const volScalarField::Boundary& pbf = p.boundaryField(); - const volScalarField::Boundary& rhobf = rho.boundaryField(); - - scalar rhoRefMax = -GREAT; - scalar rhoRefMin = GREAT; - bool rhoLimits = false; - - forAll(pbf, patchi) - { - if (pbf[patchi].fixesValue()) - { - rhoLimits = true; - - pMax_.value() = max(pMax_.value(), max(pbf[patchi])); - pMin_.value() = min(pMin_.value(), min(pbf[patchi])); - - rhoRefMax = max(rhoRefMax, max(rhobf[patchi])); - rhoRefMin = min(rhoRefMin, min(rhobf[patchi])); - } - } - - if (dict.found("pMax")) + if (dict.found("pMax") && dict.found("pMin")) { pMax_.value() = readScalar(dict.lookup("pMax")); - } - else if (dict.found("pMaxFactor")) - { - const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor"))); - pMax_ *= pMaxFactor; - } - else if (dict.found("rhoMax")) - { - // For backward-compatibility infer the pMax from rhoMax - - IOWarningInFunction(dict) - << "'rhoMax' specified rather than 'pMax' or 'pMaxFactor'" << nl - << " This is supported for backward-compatibility but " - "'pMax' or 'pMaxFactor' are more reliable." << endl; - - if (!rhoLimits) - { - FatalIOErrorInFunction(dict) - << "'rhoMax' specified rather than 'pMaxFactor'" << nl - << " but the corresponding reference density cannot" - " be evaluated from the boundary conditions." << nl - << "Please specify 'pMaxFactor' rather than 'rhoMax'" - << exit(FatalError); - } - - dimensionedScalar rhoMax("rhoMax", dimDensity, dict); - - pMax_ *= max(rhoMax.value()/rhoRefMax, 1); - } - - if (dict.found("pMin")) - { pMin_.value() = readScalar(dict.lookup("pMin")); } - else if (dict.found("pMinFactor")) + else { - const scalar pMinFactor(readScalar(dict.lookup("pMinFactor"))); - pMin_ *= pMinFactor; - } - else if (dict.found("rhoMin")) - { - // For backward-compatibility infer the pMin from rhoMin + const volScalarField::Boundary& pbf = p.boundaryField(); + const volScalarField::Boundary& rhobf = rho.boundaryField(); - IOWarningInFunction(dict) - << "'rhoMin' specified rather than 'pMin' or 'pMinFactor'" << nl - << " This is supported for backward-compatibility but" - "'pMin' or 'pMinFactor' are more reliable." << endl; + scalar rhoRefMax = -GREAT; + scalar rhoRefMin = GREAT; + bool rhoLimits = false; - if (!rhoLimits) + forAll(pbf, patchi) { - FatalIOErrorInFunction(dict) - << "'rhoMin' specified rather than 'pMinFactor'" << nl - << " but the corresponding reference density cannot" - " be evaluated from the boundary conditions." << nl - << "Please specify 'pMinFactor' rather than 'rhoMin'" - << exit(FatalError); + if (pbf[patchi].fixesValue()) + { + pLimits = true; + rhoLimits = true; + + pMax_.value() = max(pMax_.value(), max(pbf[patchi])); + pMin_.value() = min(pMin_.value(), min(pbf[patchi])); + + rhoRefMax = max(rhoRefMax, max(rhobf[patchi])); + rhoRefMin = min(rhoRefMin, min(rhobf[patchi])); + } } - dimensionedScalar rhoMin("rhoMin", dimDensity, dict); + reduce(rhoLimits, andOp()); + if (rhoLimits) + { + reduce(pMax_.value(), maxOp()); + reduce(pMin_.value(), minOp()); - pMin_ *= min(rhoMin.value()/rhoRefMin, 1); + reduce(rhoRefMax, maxOp()); + reduce(rhoRefMin, minOp()); + } + + if (dict.found("pMax")) + { + pMax_.value() = readScalar(dict.lookup("pMax")); + } + else if (dict.found("pMaxFactor")) + { + if (!pLimits) + { + FatalIOErrorInFunction(dict) + << "'pMaxFactor' specified rather than 'pMax'" << nl + << " but the corresponding reference pressure cannot" + " be evaluated from the boundary conditions." << nl + << " Please specify 'pMax' rather than 'pMaxFactor'" + << exit(FatalIOError); + } + + const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor"))); + pMax_ *= pMaxFactor; + } + else if (dict.found("rhoMax")) + { + // For backward-compatibility infer the pMax from rhoMax + + IOWarningInFunction(dict) + << "'rhoMax' specified rather than 'pMax' or 'pMaxFactor'" + << nl + << " This is supported for backward-compatibility but " + "'pMax' or 'pMaxFactor' are more reliable." << endl; + + if (!pLimits) + { + FatalIOErrorInFunction(dict) + << "'rhoMax' specified rather than 'pMax'" << nl + << " but the corresponding reference pressure cannot" + " be evaluated from the boundary conditions." << nl + << " Please specify 'pMax' rather than 'rhoMax'" + << exit(FatalIOError); + } + + if (!rhoLimits) + { + FatalIOErrorInFunction(dict) + << "'rhoMax' specified rather than 'pMaxFactor'" << nl + << " but the corresponding reference density cannot" + " be evaluated from the boundary conditions." << nl + << " Please specify 'pMaxFactor' rather than 'rhoMax'" + << exit(FatalIOError); + } + + dimensionedScalar rhoMax("rhoMax", dimDensity, dict); + + pMax_ *= max(rhoMax.value()/rhoRefMax, 1); + } + + if (dict.found("pMin")) + { + pMin_.value() = readScalar(dict.lookup("pMin")); + } + else if (dict.found("pMinFactor")) + { + if (!pLimits) + { + FatalIOErrorInFunction(dict) + << "'pMinFactor' specified rather than 'pMin'" << nl + << " but the corresponding reference pressure cannot" + " be evaluated from the boundary conditions." << nl + << " Please specify 'pMin' rather than 'pMinFactor'" + << exit(FatalIOError); + } + + const scalar pMinFactor(readScalar(dict.lookup("pMinFactor"))); + pMin_ *= pMinFactor; + } + else if (dict.found("rhoMin")) + { + // For backward-compatibility infer the pMin from rhoMin + + IOWarningInFunction(dict) + << "'rhoMin' specified rather than 'pMin' or 'pMinFactor'" << nl + << " This is supported for backward-compatibility but" + "'pMin' or 'pMinFactor' are more reliable." << endl; + + if (!pLimits) + { + FatalIOErrorInFunction(dict) + << "'rhoMin' specified rather than 'pMin'" << nl + << " but the corresponding reference pressure cannot" + " be evaluated from the boundary conditions." << nl + << " Please specify 'pMin' rather than 'rhoMin'" + << exit(FatalIOError); + } + + if (!rhoLimits) + { + FatalIOErrorInFunction(dict) + << "'rhoMin' specified rather than 'pMinFactor'" << nl + << " but the corresponding reference density cannot" + " be evaluated from the boundary conditions." << nl + << " Please specify 'pMinFactor' rather than 'rhoMin'" + << exit(FatalIOError); + } + + dimensionedScalar rhoMin("rhoMin", dimDensity, dict); + + pMin_ *= min(rhoMin.value()/rhoRefMin, 1); + } } Info<< "pressureControl" << nl diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H index 19bd053d4..61f510ac6 100644 --- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H +++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H @@ -76,7 +76,8 @@ public: ( const volScalarField& p, const volScalarField& rho, - const dictionary& dict + const dictionary& dict, + const bool pRefRequired = true ); diff --git a/tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/system/fvSolution b/tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/system/fvSolution index 66786f377..c771e4cfe 100644 --- a/tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/system/fvSolution +++ b/tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/system/fvSolution @@ -63,8 +63,9 @@ PIMPLE nOuterCorrectors 3; nCorrectors 1; nNonOrthogonalCorrectors 0; - rhoMin 0.5; - rhoMax 2.0; + + pMaxFactor 1.2; + pMinFactor 0.8; } relaxationFactors diff --git a/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSolution index 4eddb3872..b2dc2e505 100644 --- a/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSolution +++ b/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSolution @@ -52,8 +52,9 @@ PIMPLE nOuterCorrectors 3; nCorrectors 1; nNonOrthogonalCorrectors 0; - rhoMin 0.5; - rhoMax 2.0; + + pMinFactor 0.5; + pMaxFactor 2.0; } relaxationFactors diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/constant/thermophysicalProperties b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/constant/thermophysicalProperties index 32542b344..7fc31fc96 100644 --- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/constant/thermophysicalProperties +++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/constant/thermophysicalProperties @@ -23,7 +23,7 @@ thermoType thermo hConst; equationOfState perfectGas; specie specie; - energy sensibleEnthalpy; + energy sensibleInternalEnergy; } mixture diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes index 359d6ad3e..c0522b83d 100644 --- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes +++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes @@ -30,8 +30,9 @@ divSchemes default none; div(phi,U) Gauss upwind; div(phid,p) Gauss upwind; + div(phiv,p) Gauss linear; div(phi,K) Gauss linear; - div(phi,h) Gauss upwind; + div(phi,e) Gauss upwind; div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind; div(phi,R) Gauss upwind; diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution index ff2440678..3c7d9238f 100644 --- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution +++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution @@ -28,11 +28,10 @@ solvers pFinal { $p; - tolerance 1e-07; relTol 0; } - "(rho|U|h|k|epsilon|omega)" + "(rho|U|e|k|epsilon|omega)" { solver smoothSolver; smoother symGaussSeidel; @@ -40,10 +39,9 @@ solvers relTol 0.1; } - "(rho|U|h|k|epsilon|omega)Final" + "(rho|U|e|k|epsilon|omega)Final" { $U; - tolerance 1e-06; relTol 0; } } @@ -57,8 +55,8 @@ PIMPLE nNonOrthogonalCorrectors 0; consistent yes; - rhoMin 0.4; - rhoMax 2.0; + pMaxFactor 1.5; + pMinFactor 0.9; residualControl { @@ -76,13 +74,13 @@ relaxationFactors { fields { - "p.*" 0.9; + "p.*" 1; "rho.*" 1; } equations { "U.*" 0.9; - "h.*" 0.7; + "e.*" 0.7; "(k|epsilon|omega).*" 0.8; } } diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuctLTS/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuctLTS/system/fvSolution index 83e599258..8b4068e0d 100644 --- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuctLTS/system/fvSolution +++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuctLTS/system/fvSolution @@ -56,8 +56,8 @@ PIMPLE nCorrectors 1; nNonOrthogonalCorrectors 0; - rhoMin 0.5; - rhoMax 2.0; + pMaxFactor 1.5; + pMinFactor 0.9; maxCo 0.2; rDeltaTSmoothingCoeff 0.1; diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/cavity/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/cavity/system/fvSolution index e7c64b955..17ed812b9 100644 --- a/tutorials/compressible/rhoPimpleFoam/RAS/cavity/system/fvSolution +++ b/tutorials/compressible/rhoPimpleFoam/RAS/cavity/system/fvSolution @@ -59,8 +59,9 @@ PIMPLE nOuterCorrectors 1; nCorrectors 2; nNonOrthogonalCorrectors 0; - rhoMin 0.5; - rhoMax 2.0; + + pMax 1.2e5; + pMin 0.8e5; } diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/fvSolution index d4a55d13c..02ea29588 100644 --- a/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/fvSolution +++ b/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/fvSolution @@ -84,14 +84,13 @@ solvers PIMPLE { - nOuterCorrectors 1; - nCorrectors 2; - nNonOrthogonalCorrectors 0; - momentumPredictor yes; - rhoMin 0.5; - rhoMax 2.0; - pRefCell 0; - pRefValue 1e5; + nOuterCorrectors 1; + nCorrectors 2; + nNonOrthogonalCorrectors 0; + momentumPredictor yes; + + pMax 1.2e5; + pMin 0.8e5; } relaxationFactors diff --git a/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSolution index 4eddb3872..b2dc2e505 100644 --- a/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSolution +++ b/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSolution @@ -52,8 +52,9 @@ PIMPLE nOuterCorrectors 3; nCorrectors 1; nNonOrthogonalCorrectors 0; - rhoMin 0.5; - rhoMax 2.0; + + pMinFactor 0.5; + pMaxFactor 2.0; } relaxationFactors