From 27ab3edc5ebf0a03f317e543b88a9646277c8075 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 5 Mar 2020 09:22:05 +0000 Subject: [PATCH] buoyantPimpleFoam, chtMultiRegionFoam, rhoReactingBuoyantFoam: Share pEqn.H --- .../rhoReactingBuoyantFoam/Make/options | 1 + .../rhoReactingBuoyantFoam/UEqn.H | 3 +- .../rhoReactingBuoyantFoam/createFields.H | 19 ++- .../rhoReactingBuoyantFoam/pEqn.H | 74 --------- .../rhoReactingBuoyantFoam.C | 5 +- .../heatTransfer/buoyantPimpleFoam/UEqn.H | 3 +- .../buoyantPimpleFoam/createFields.H | 2 + .../heatTransfer/buoyantPimpleFoam/pEqn.H | 154 +++++++++++++----- .../chtMultiRegionFoam/fluid/pEqn.H | 145 ----------------- .../fluid/setRegionFluidFields.H | 4 + .../chtMultiRegionFoam/fluid/solveFluid.H | 2 +- 11 files changed, 145 insertions(+), 267 deletions(-) delete mode 100644 applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options index b403a0268e..57f5c58466 100644 --- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options +++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options @@ -2,6 +2,7 @@ EXE_INC = \ -I. \ -I$(FOAM_SOLVERS)/combustion/reactingFoam \ -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/sampling/lnInclude \ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/UEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/UEqn.H index dbda898873..e6bb14c76b 100644 --- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/UEqn.H +++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/UEqn.H @@ -1,6 +1,6 @@ MRF.correctBoundaryVelocity(U); - fvVectorMatrix UEqn + tmp tUEqn ( fvm::ddt(rho, U) + fvm::div(phi, U) + MRF.DDt(rho, U) @@ -8,6 +8,7 @@ == fvOptions(rho, U) ); + fvVectorMatrix& UEqn = tUEqn.ref(); UEqn.relax(); diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H index 0a4284e82c..7274bf121d 100644 --- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H +++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H @@ -89,6 +89,17 @@ volScalarField p_rgh // Force p_rgh to be consistent with p p_rgh = p - rho*gh; +pressureControl pressureControl +( + p, + p_rgh, + rho, + pimple.dict(), + thermo.incompressible() +); + +mesh.setFluxRequired(p_rgh.name()); + Info<< "Creating field dpdt\n" << endl; volScalarField dpdt ( @@ -105,9 +116,10 @@ volScalarField dpdt Info<< "Creating field kinetic energy K\n" << endl; volScalarField K("K", 0.5*magSqr(U)); +dimensionedScalar initialMass = fvc::domainIntegrate(rho); + multivariateSurfaceInterpolationScheme::fieldTable fields; - forAll(Y, i) { fields.add(Y[i]); @@ -116,3 +128,8 @@ fields.add(thermo.he()); #include "createMRF.H" #include "createFvOptions.H" + + +// This solver does not support moving mesh but it uses the pressure equation +// of one which does, so we need a dummy face-momentum field +autoPtr rhoUf(nullptr); diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H deleted file mode 100644 index 1ef1e4d8a2..0000000000 --- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H +++ /dev/null @@ -1,74 +0,0 @@ -rho = thermo.rho(); - -// Thermodynamic density needs to be updated by psi*d(p) after the -// pressure solution -const volScalarField psip0(psi*p); - -volScalarField rAU(1.0/UEqn.A()); -surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); -volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); - -surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); - -surfaceScalarField phiHbyA -( - "phiHbyA", - ( - fvc::interpolate(rho)*fvc::flux(HbyA) - + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi)) - ) - + phig -); - -MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - -// Update the pressure BCs to ensure flux consistency -constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); - -fvScalarMatrix p_rghDDtEqn -( - fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) - + fvc::div(phiHbyA) - == - fvOptions(psi, p_rgh, rho.name()) -); - -while (pimple.correctNonOrthogonal()) -{ - fvScalarMatrix p_rghEqn - ( - p_rghDDtEqn - - fvm::laplacian(rhorAUf, p_rgh) - ); - - p_rghEqn.solve(); - - if (pimple.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); - K = 0.5*magSqr(U); - } -} - -p = p_rgh + rho*gh; - -// Thermodynamic density update -thermo.correctRho(psi*p - psip0); - -if (thermo.dpdt()) -{ - dpdt = fvc::ddt(p); -} - -#include "rhoEqn.H" -#include "compressibleContinuityErrs.H" diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C index ec3dde4f72..ec9b3cc559 100644 --- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C +++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2012-2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -36,6 +36,7 @@ Description #include "turbulentFluidThermoModel.H" #include "multivariateScheme.H" #include "pimpleControl.H" +#include "pressureControl.H" #include "fvOptions.H" #include "localEulerDdtScheme.H" #include "fvcSmooth.H" @@ -97,7 +98,7 @@ int main(int argc, char *argv[]) // --- Pressure corrector loop while (pimple.correct()) { - #include "pEqn.H" + #include "../../../heatTransfer/buoyantPimpleFoam/pEqn.H" } if (pimple.turbCorr()) diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H index 6ff0c80867..832094b597 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H @@ -2,7 +2,7 @@ MRF.correctBoundaryVelocity(U); - fvVectorMatrix UEqn + tmp tUEqn ( fvm::ddt(rho, U) + fvm::div(phi, U) + MRF.DDt(rho, U) @@ -10,6 +10,7 @@ == fvOptions(rho, U) ); + fvVectorMatrix& UEqn = tUEqn.ref(); UEqn.relax(); diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H index 70453aa3dd..58d825f881 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H @@ -102,6 +102,8 @@ volScalarField dpdt Info<< "Creating field kinetic energy K\n" << endl; volScalarField K("K", 0.5*magSqr(U)); +dimensionedScalar initialMass = fvc::domainIntegrate(rho); + #include "createMRF.H" #include "createRadiationModel.H" #include "createFvOptions.H" diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H index 93057220e4..3499060d41 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H @@ -1,4 +1,4 @@ -if (!pimple.simpleRho()) +if (!mesh.steady() && !pimple.simpleRho()) { rho = thermo.rho(); } @@ -11,17 +11,22 @@ volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh)); +if (pimple.nCorrPiso() <= 1) +{ + tUEqn.clear(); +} + surfaceScalarField phiHbyA ( "phiHbyA", - ( - fvc::interpolate(rho)*fvc::flux(HbyA) - + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf)) - ) + fvc::interpolate(rho)*fvc::flux(HbyA) + + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf)) ); MRF.makeRelative(fvc::interpolate(rho), phiHbyA); +bool adjustMass = false; + surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); phiHbyA += phig; @@ -31,66 +36,128 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); fvc::makeRelative(phiHbyA, rho, U); -fvScalarMatrix p_rghDDtEqn -( - fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) - + fvc::div(phiHbyA) - == - fvOptions(psi, p_rgh, rho.name()) -); +fvScalarMatrix p_rghEqn(p_rgh, dimMass/dimTime); -while (pimple.correctNonOrthogonal()) +if (pimple.transonic()) { - fvScalarMatrix p_rghEqn(p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh)); - - p_rghEqn.setReference + surfaceScalarField phid ( - pressureControl.refCell(), - pressureControl.refValue() + "phid", + (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA ); - p_rghEqn.solve(); + phiHbyA -= fvc::interpolate(psi*p_rgh)*phiHbyA/fvc::interpolate(rho); - if (pimple.finalNonOrthogonalIter()) + fvScalarMatrix p_rghDDtEqn + ( + fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) + + fvc::div(phiHbyA) + fvm::div(phid, p_rgh) + == + fvOptions(psi, p_rgh, rho.name()) + ); + + while (pimple.correctNonOrthogonal()) { - // Calculate the conservative fluxes - phi = phiHbyA + p_rghEqn.flux(); + p_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh); - // Explicitly relax pressure for momentum corrector - p_rgh.relax(); + // Relax the pressure equation to ensure diagonal-dominance + p_rghEqn.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); - K = 0.5*magSqr(U); + p_rghEqn.setReference + ( + pressureControl.refCell(), + pressureControl.refValue() + ); + + p_rghEqn.solve(); + } +} +else +{ + if (mesh.steady()) + { + adjustMass = adjustPhi(phiHbyA, U, p_rgh); + } + + fvScalarMatrix p_rghDDtEqn + ( + fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) + + fvc::div(phiHbyA) + == + fvOptions(psi, p_rgh, rho.name()) + ); + + while (pimple.correctNonOrthogonal()) + { + p_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh); + + p_rghEqn.setReference + ( + pressureControl.refCell(), + pressureControl.refValue() + ); + + p_rghEqn.solve(); } } +phi = phiHbyA + p_rghEqn.flux(); + p = p_rgh + rho*gh; -bool limitedp = pressureControl.limit(p); - -if (limitedp) +if (mesh.steady()) { + #include "incompressible/continuityErrs.H" +} +else +{ + const bool limitedp = pressureControl.limit(p); + + // Thermodynamic density update + thermo.correctRho(psi*p - psip0); + + if (limitedp) + { + rho = thermo.rho(); + } + + #include "rhoEqn.H" + #include "compressibleContinuityErrs.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); +K = 0.5*magSqr(U); + +if (mesh.steady()) +{ + pressureControl.limit(p); +} + +// For steady closed-volume compressible cases adjust the pressure level +// to obey overall mass continuity +if (adjustMass && !thermo.incompressible()) +{ + p += (initialMass - fvc::domainIntegrate(thermo.rho())) + /fvc::domainIntegrate(psi); p_rgh = p - rho*gh; + p.correctBoundaryConditions(); } -// Thermodynamic density update -thermo.correctRho(psi*p - psip0); - -if (limitedp) +if (mesh.steady() || pimple.simpleRho() || adjustMass) { rho = thermo.rho(); } -#include "rhoEqn.H" -#include "compressibleContinuityErrs.H" - -if (pimple.simpleRho()) +if ((mesh.steady() || pimple.simpleRho()) && !pimple.transonic()) { - rho = thermo.rho(); + rho.relax(); } // Correct rhoUf if the mesh is moving @@ -105,3 +172,6 @@ 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/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H deleted file mode 100644 index 9e89597ab4..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H +++ /dev/null @@ -1,145 +0,0 @@ -if (!mesh.steady() && !pimple.simpleRho()) -{ - rho = thermo.rho(); -} - -volScalarField rAU("rAU", 1.0/UEqn.A()); -surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); -volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh)); -if (pimple.nCorrPiso() <= 1) -{ - tUEqn.clear(); -} - -surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); - -surfaceScalarField phiHbyA -( - "phiHbyA", - fvc::interpolate(rho)*fvc::flux(HbyA) - + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi)) -); - -MRF.makeRelative(fvc::interpolate(rho), phiHbyA); - -const bool closedVolume = mesh.steady() && adjustPhi(phiHbyA, U, p_rgh); -const bool adjustMass = closedVolume && !thermo.incompressible(); - -phiHbyA += phig; - -// Update the pressure BCs to ensure flux consistency -constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); - -{ - fvScalarMatrix p_rghEqnComp - ( - fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) - == - fvOptions(psi, p_rgh, rho.name()) - ); - - if (pimple.transonic()) - { - surfaceScalarField phid - ( - "phid", - (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA - ); - - phiHbyA -= fvc::interpolate(psi*p_rgh)*phiHbyA/fvc::interpolate(rho); - - p_rghEqnComp += fvm::div(phid, p_rgh); - } - - // Thermodynamic density needs to be updated by psi*d(p) after the - // pressure solution - tmp psip0(mesh.steady() ? tmp() : psi*p); - - while (pimple.correctNonOrthogonal()) - { - fvScalarMatrix p_rghEqnIncomp - ( - fvc::div(phiHbyA) - - fvm::laplacian(rhorAUf, p_rgh) - ); - - fvScalarMatrix p_rghEqn(p_rghEqnComp + p_rghEqnIncomp); - - p_rghEqn.setReference - ( - pressureControl.refCell(), - pressureControl.refValue() - ); - - p_rghEqn.solve(); - - if (pimple.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_rghEqnIncomp.flux())/rhorAUf); - U.correctBoundaryConditions(); - fvOptions.correct(U); - K = 0.5*magSqr(U); - } - } - - p = p_rgh + rho*gh; - - pressureControl.limit(p); - - // Thermodynamic density update - if (!mesh.steady()) - { - thermo.correctRho(psi*p - psip0); - } -} - -// Update pressure time derivative if needed -if (thermo.dpdt()) -{ - dpdt = fvc::ddt(p); -} - -// Solve continuity -if (!mesh.steady()) -{ - #include "rhoEqn.H" - #include "compressibleContinuityErrs.H" -} -else -{ - #include "incompressible/continuityErrs.H" -} - - -// For closed-volume compressible cases adjust the pressure level -// to obey overall mass continuity -if (adjustMass) -{ - p += (initialMass - fvc::domainIntegrate(thermo.rho())) - /fvc::domainIntegrate(psi); - p_rgh = p - rho*gh; - p.correctBoundaryConditions(); -} - -// Density updates -if (adjustMass || mesh.steady() || pimple.simpleRho()) -{ - rho = thermo.rho(); -} - -if (mesh.steady() && !pimple.transonic()) -{ - rho.relax(); -} - -Info<< "Min/max rho:" << min(rho).value() << ' ' - << max(rho).value() << endl; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H index 0e08ee5249..f0fdc7c2f9 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H @@ -60,3 +60,7 @@ pressureControl& pressureControl = pressureControlFluid[i]; scalar cumulativeContErr = cumulativeContErrs[i]; + + // This solver does not support moving mesh but it uses the pressure + // equation of one which does, so we need a dummy face-momentum field + autoPtr rhoUf(nullptr); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H index 779b857294..137c18d230 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H @@ -16,7 +16,7 @@ else // --- PISO loop while (pimple.correct()) { - #include "pEqn.H" + #include "../../buoyantPimpleFoam/pEqn.H" } if (pimples.pimpleTurbCorr(i))