diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H index 29cc3274cb..90392ed475 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H @@ -40,7 +40,13 @@ volVectorField U #include "compressibleCreatePhi.H" -pressureControl pressureControl(p, rho, pimple.dict(), false); +pressureControl pressureControl +( + p, + rho, + pimple.dict(), + thermo.incompressible() +); mesh.setFluxRequired(p.name()); diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H index 6822efa11a..1d525d7d76 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H +++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H @@ -38,7 +38,13 @@ volVectorField U #include "compressibleCreatePhi.H" -pressureControl pressureControl(p, rho, simple.dict()); +pressureControl pressureControl +( + p, + rho, + simple.dict(), + thermo.incompressible() +); mesh.setFluxRequired(p.name()); diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index 064c46cb5e..e7571bfb26 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -27,11 +27,16 @@ volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); tUEqn.clear(); -bool closedVolume = false; +surfaceScalarField phiHbyA +( + "phiHbyA", + fvc::interpolate(rho)*fvc::flux(HbyA) +); -surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA)); MRF.makeRelative(fvc::interpolate(rho), phiHbyA); +bool closedVolume = false; + // Update the pressure BCs to ensure flux consistency constrainPressure(p, rho, U, phiHbyA, rhorAAtUf, MRF); @@ -127,7 +132,7 @@ pressureControl.limit(p); // For closed-volume cases adjust the pressure and density levels // to obey overall mass continuity -if (closedVolume) +if (closedVolume && !thermo.incompressible()) { p += (initialMass - fvc::domainIntegrate(psi*p)) /fvc::domainIntegrate(psi); diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H index 3499060d41..d29960b306 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H @@ -7,8 +7,8 @@ if (!mesh.steady() && !pimple.simpleRho()) // pressure solution const volScalarField psip0(psi*p); -volScalarField rAU(1.0/UEqn.A()); -surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); +const volScalarField rAU("rAU", 1.0/UEqn.A()); +const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh)); if (pimple.nCorrPiso() <= 1) @@ -25,7 +25,7 @@ surfaceScalarField phiHbyA MRF.makeRelative(fvc::interpolate(rho), phiHbyA); -bool adjustMass = false; +bool adjustMass = pimple.transonic() ? false : adjustPhi(phiHbyA, U, p_rgh); surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); @@ -172,6 +172,3 @@ if (thermo.dpdt()) dpdt -= fvc::div(fvc::meshPhi(rho, U), p); } } - -Info<< "Min/max rho:" << min(rho).value() << ' ' - << max(rho).value() << endl; diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options b/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options index 5170ca5d57..2d14ae533c 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options @@ -1,5 +1,6 @@ EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ -I$(LIB_SRC)/sampling/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/transportModels/compressible/lnInclude \ diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C index f5f2f099c4..62f5b071c4 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -35,6 +35,7 @@ Description #include "turbulentFluidThermoModel.H" #include "radiationModel.H" #include "simpleControl.H" +#include "pressureControl.H" #include "fvOptions.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H index 72c828d43b..02b0811887 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H @@ -70,21 +70,18 @@ volScalarField p_rgh // Force p_rgh to be consistent with p p_rgh = p - rho*gh; -label pRefCell = 0; -scalar pRefValue = 0.0; -setRefCell +pressureControl pressureControl ( p, p_rgh, + rho, simple.dict(), - pRefCell, - pRefValue + thermo.incompressible() ); mesh.setFluxRequired(p_rgh.name()); dimensionedScalar initialMass = fvc::domainIntegrate(rho); -dimensionedScalar totalVolume = sum(mesh.V()); #include "createMRF.H" #include "createRadiationModel.H" diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H index 02bf5c9614..4a08c51ad8 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H @@ -1,70 +1,110 @@ +rho = thermo.rho(); + +const volScalarField rAU("rAU", 1.0/UEqn.A()); +const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); +volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh)); + +tUEqn.clear(); + +surfaceScalarField phiHbyA +( + "phiHbyA", + fvc::interpolate(rho)*fvc::flux(HbyA) +); + +MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + +bool closedVolume = simple.transonic() ? false : adjustPhi(phiHbyA, U, p_rgh); + +surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + +phiHbyA += phig; + +// Update the pressure BCs to ensure flux consistency +constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); + +fvScalarMatrix p_rghEqn(p_rgh, dimMass/dimTime); + +if (simple.transonic()) { - rho = thermo.rho(); - rho.relax(); - - volScalarField rAU("rAU", 1.0/UEqn.A()); - surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); - volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh)); - tUEqn.clear(); - - surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); - - surfaceScalarField phiHbyA + surfaceScalarField phid ( - "phiHbyA", - fvc::interpolate(rho)*fvc::flux(HbyA) + "phid", + (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA ); - MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - - bool closedVolume = adjustPhi(phiHbyA, U, p_rgh); - - phiHbyA += phig; - - // Update the pressure BCs to ensure flux consistency - constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); + phiHbyA -= fvc::interpolate(psi*p_rgh)*phiHbyA/fvc::interpolate(rho); while (simple.correctNonOrthogonal()) { - fvScalarMatrix p_rghEqn + p_rghEqn = + fvc::div(phiHbyA) + fvm::div(phid, p_rgh) + - fvm::laplacian(rhorAUf, p_rgh) + == + fvOptions(psi, p_rgh, rho.name()); + + // Relax the pressure equation to ensure diagonal-dominance + p_rghEqn.relax(); + + p_rghEqn.setReference ( - fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA) + pressureControl.refCell(), + pressureControl.refValue() ); - p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.solve(); - - if (simple.finalNonOrthogonalIter()) - { - // Calculate the conservative fluxes - phi = phiHbyA - p_rghEqn.flux(); - - // Explicitly relax pressure for momentum corrector - p_rgh.relax(); - - // Correct the momentum source with the pressure gradient flux - // calculated from the relaxed pressure - U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf); - U.correctBoundaryConditions(); - fvOptions.correct(U); - } } - - #include "continuityErrs.H" - - p = p_rgh + rho*gh; - - // For closed-volume cases adjust the pressure level - // to obey overall mass continuity - if (!thermo.incompressible() && closedVolume) - { - p += (initialMass - fvc::domainIntegrate(psi*p)) - /fvc::domainIntegrate(psi); - p_rgh = p - rho*gh; - } - - rho = thermo.rho(); - rho.relax(); - Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() - << endl; +} +else +{ + while (simple.correctNonOrthogonal()) + { + p_rghEqn = + fvc::div(phiHbyA) + - fvm::laplacian(rhorAUf, p_rgh) + == + fvOptions(psi, p_rgh, rho.name()); + + p_rghEqn.setReference + ( + pressureControl.refCell(), + pressureControl.refValue() + ); + + p_rghEqn.solve(); + } +} + +phi = phiHbyA + p_rghEqn.flux(); + +p = p_rgh + rho*gh; + +#include "incompressible/continuityErrs.H" + +// Explicitly relax pressure for momentum corrector +p_rgh.relax(); + +// Correct the momentum source with the pressure gradient flux +// calculated from the relaxed pressure +U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf); +U.correctBoundaryConditions(); +fvOptions.correct(U); + +pressureControl.limit(p); + +// For closed-volume compressible cases adjust the pressure level +// to obey overall mass continuity +if (closedVolume && !thermo.incompressible()) +{ + p += (initialMass - fvc::domainIntegrate(thermo.rho())) + /fvc::domainIntegrate(psi); + p_rgh = p - rho*gh; + p.correctBoundaryConditions(); +} + +rho = thermo.rho(); + +if (!simple.transonic()) +{ + rho.relax(); }