From aadacc5812c7c39be87d98246ac546bb414ba848 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 6 Oct 2010 17:41:40 +0100 Subject: [PATCH 01/15] BUG: corrected laminar clause in kappt wall function --- ...ppatJayatillekeWallFunctionFvPatchScalarField.C | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C index 7aa9fa2fe9..6efef7364b 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C @@ -231,19 +231,17 @@ void kappatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs() scalar P = Psmooth(Prat); scalar yPlusTherm = this->yPlusTherm(P, Prat); - // Evaluate new effective thermal diffusivity - scalar kappaEff = 0.0; - if (yPlus < yPlusTherm) + // Update turbulent thermal conductivity + if (yPlus > yPlusTherm) { - kappaEff = Pr*yPlus; + scalar nu = nuw[faceI]; + scalar kt = nu*(yPlus/(Prt_/kappa_*log(E_*yPlusTherm) + P) - 1/Pr); + kappatw[faceI] = max(0.0, kt); } else { - kappaEff = nuw[faceI]*yPlus/(Prt_/kappa_*log(E_*yPlusTherm) + P); + kappatw[faceI] = 0.0; } - - // Update turbulent thermal diffusivity - kappatw[faceI] = max(0.0, kappaEff - nuw[faceI]/Pr); } fixedValueFvPatchField::updateCoeffs(); From 590e2c3a5253a9b5b69dca192a408cc78e388f73 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 6 Oct 2010 17:47:53 +0100 Subject: [PATCH 02/15] ENH: Updated heat transfer solvers and tutorials to use p_rgh --- .../buoyantBoussinesqPimpleFoam/TEqn.H | 2 +- .../buoyantBoussinesqPimpleFoam/UEqn.H | 9 +- .../buoyantBoussinesqPimpleFoam.C | 2 +- .../createFields.H | 56 ++- .../buoyantBoussinesqPimpleFoam/pEqn.H | 34 +- .../createFields.H | 18 +- .../buoyantBoussinesqSimpleFoam/pEqn.H | 15 +- .../heatTransfer/buoyantPimpleFoam/UEqn.H | 9 +- .../buoyantPimpleFoam/buoyantPimpleFoam.C | 2 +- .../buoyantPimpleFoam/createFields.H | 27 +- .../heatTransfer/buoyantPimpleFoam/hEqn.H | 2 +- .../heatTransfer/buoyantPimpleFoam/pEqn.H | 38 +- .../buoyantSimpleFoam/buoyantSimpleFoam.C | 3 - .../buoyantSimpleFoam/createFields.H | 68 ++- .../heatTransfer/buoyantSimpleFoam/pEqn.H | 9 +- .../buoyantSimpleRadiationFoam.C | 2 +- .../chtMultiRegionFoam/chtMultiRegionFoam.C | 2 + .../chtMultiRegionSimpleFoam/fluid/UEqn.H | 6 +- .../fluid/createFluidFields.H | 59 ++- .../chtMultiRegionSimpleFoam/fluid/hEqn.H | 4 +- .../chtMultiRegionSimpleFoam/fluid/pEqn.H | 53 ++- .../fluid/setRegionFluidFields.H | 5 +- .../fluid/solveFluid.H | 2 +- .../chtMultiRegionFoam/fluid/UEqn.H | 9 +- .../fluid/createFluidFields.H | 41 +- .../chtMultiRegionFoam/fluid/hEqn.H | 13 +- .../chtMultiRegionFoam/fluid/pEqn.H | 42 +- .../fluid/readFluidMultiRegionPISOControls.H | 17 - .../fluid/setRegionFluidFields.H | 14 +- .../chtMultiRegionFoam/fluid/solveFluid.H | 10 + .../fluid/storeOldFluidFields.H | 2 +- .../solid/createSolidFields.H | 6 +- .../chtMultiRegionFoam/solid/solveSolid.H | 12 +- .../lagrangian/coalChemistryFoam/UEqn.H | 2 +- .../coalChemistryFoam/coalChemistryFoam.C | 15 +- .../lagrangian/coalChemistryFoam/hsEqn.H | 2 +- .../lagrangian/coalChemistryFoam/pEqn.H | 30 +- .../hotRoom/0/T.org | 406 +----------------- .../buoyantBoussinesqPimpleFoam/hotRoom/0/p | 15 +- .../hotRoom/0/p_rgh | 45 ++ .../hotRoom/Allclean | 2 +- .../hotRoom/Allrun | 1 + .../hotRoom/system/fvSchemes | 16 +- .../hotRoom/system/fvSolution | 8 +- .../buoyantBoussinesqSimpleFoam/hotRoom/0/p | 15 +- .../hotRoom/Allclean | 2 +- .../hotRoom/Allrun | 1 + .../hotRoom/system/controlDict | 20 + .../hotRoom/system/fvSolution | 6 +- .../iglooWithFridges/0/p | 20 +- .../constant/polyMesh/boundary | 18 +- .../iglooWithFridges/system/controlDict | 20 + .../iglooWithFridges/system/fvSolution | 7 +- .../buoyantPimpleFoam/hotRoom/0/p | 12 +- .../buoyantPimpleFoam/hotRoom/0/p_rgh | 42 ++ .../hotRoom/system/fvSchemes | 4 +- .../hotRoom/system/fvSolution | 12 +- .../buoyantSimpleFoam/buoyantCavity/0/p | 16 +- .../buoyantCavity/system/controlDict | 20 + .../buoyantCavity/system/fvSolution | 11 +- .../buoyantSimpleFoam/hotRoom/0/T | 406 +----------------- .../buoyantSimpleFoam/hotRoom/0/T.org | 406 +----------------- .../buoyantSimpleFoam/hotRoom/0/p | 12 +- .../buoyantSimpleFoam/hotRoom/0/p_rgh | 42 ++ .../hotRoom/system/controlDict | 20 + .../hotRoom/system/fvSchemes | 14 +- .../hotRoom/system/fvSolution | 16 +- .../hotRadiationRoom/0/G | 8 +- .../hotRadiationRoom/0/p | 16 +- .../hotRadiationRoom/0/p_rgh | 48 +++ .../hotRadiationRoom/system/controlDict | 20 + .../hotRadiationRoom/system/fvSchemes | 4 +- .../hotRadiationRoom/system/fvSolution | 14 +- .../hotRadiationRoomFvDOM/0/IDefault | 2 +- .../hotRadiationRoomFvDOM/0/p | 8 +- .../hotRadiationRoomFvDOM/0/p_rgh | 48 +++ .../hotRadiationRoomFvDOM/system/controlDict | 20 + .../hotRadiationRoomFvDOM/system/fvSchemes | 4 +- .../hotRadiationRoomFvDOM/system/fvSolution | 9 +- .../multiRegionHeater/0/epsilon | 1 - .../chtMultiRegionFoam/multiRegionHeater/0/k | 1 - .../multiRegionHeater/0/p_rgh | 29 ++ .../system/bottomAir/changeDictionaryDict | 20 +- .../system/bottomAir/fvSchemes | 6 +- .../system/bottomAir/fvSolution | 19 +- .../multiRegionHeater/system/controlDict | 9 +- .../system/heater/decomposeParDict | 2 +- .../system/heater/fvSolution | 6 + .../system/leftSolid/fvSolution | 6 + .../system/rightSolid/fvSolution | 6 + .../system/topAir/changeDictionaryDict | 39 +- .../multiRegionHeater/system/topAir/fvSchemes | 6 +- .../system/topAir/fvSolution | 43 +- .../snappyMultiRegionHeater/0/alphat | 1 - .../snappyMultiRegionHeater/0/epsilon | 1 - .../snappyMultiRegionHeater/0/k | 1 - .../snappyMultiRegionHeater/0/p_rgh | 29 ++ .../system/bottomAir/changeDictionaryDict | 20 +- .../system/bottomAir/decomposeParDict | 2 +- .../system/bottomAir/fvSchemes | 6 +- .../system/bottomAir/fvSolution | 42 +- .../system/heater/decomposeParDict | 2 +- .../system/heater/fvSolution | 6 + .../system/leftSolid/decomposeParDict | 2 +- .../system/leftSolid/fvSolution | 7 + .../system/rightSolid/decomposeParDict | 2 +- .../system/rightSolid/fvSolution | 7 + .../system/topAir/changeDictionaryDict | 39 +- .../system/topAir/decomposeParDict | 2 +- .../system/topAir/fvSchemes | 6 +- .../system/topAir/fvSolution | 43 +- .../multiRegionHeater/0/U | 2 +- .../multiRegionHeater/0/epsilon | 1 - .../multiRegionHeater/0/k | 1 - .../multiRegionHeater/0/p_rgh | 29 ++ .../system/bottomAir/changeDictionaryDict | 20 +- .../system/bottomAir/decomposeParDict | 2 +- .../system/bottomAir/fvSchemes | 6 +- .../system/bottomAir/fvSolution | 13 +- .../multiRegionHeater/system/controlDict | 7 +- .../multiRegionHeater/system/decomposeParDict | 2 +- .../system/heater/decomposeParDict | 2 +- .../system/leftSolid/decomposeParDict | 2 +- .../system/rightSolid/decomposeParDict | 2 +- .../system/topAir/changeDictionaryDict | 34 +- .../system/topAir/decomposeParDict | 2 +- .../multiRegionHeater/system/topAir/fvSchemes | 6 +- .../system/topAir/fvSolution | 15 +- 128 files changed, 1255 insertions(+), 1780 deletions(-) delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H create mode 100644 tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p_rgh create mode 100644 tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p_rgh create mode 100644 tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/p_rgh create mode 100644 tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p_rgh create mode 100644 tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/p_rgh create mode 100644 tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/p_rgh create mode 100644 tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/p_rgh create mode 100644 tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/p_rgh diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H index dbfc61739f..9a835792a4 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H @@ -12,7 +12,7 @@ ); TEqn.relax(); - TEqn.solve(); + TEqn.solve(mesh.solver(T.select(finalIter))); rhok = 1.0 - beta*(T - TRef); } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H index 35387f4ff4..20a05e5cd4 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H @@ -18,9 +18,10 @@ fvc::reconstruct ( ( - fvc::interpolate(rhok)*(g & mesh.Sf()) - - fvc::snGrad(p)*mesh.magSf() - ) - ) + - ghf*fvc::snGrad(rhok) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ), + mesh.solver(U.select(finalIter)) ); } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C index ebf68d409c..54519906a4 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C @@ -87,7 +87,7 @@ int main(int argc, char *argv[]) if (nOuterCorr != 1) { - p.storePrevIter(); + p_rgh.storePrevIter(); } #include "UEqn.H" diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H index 23d20cfa96..342af079d8 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H @@ -14,12 +14,12 @@ mesh ); - Info<< "Reading field p\n" << endl; - volScalarField p + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh ( IOobject ( - "p", + "p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ, @@ -52,6 +52,18 @@ incompressible::RASModel::New(U, phi, laminarTransport) ); + // Kinematic density for buoyancy force + volScalarField rhok + ( + IOobject + ( + "rhok", + runTime.timeName(), + mesh + ), + 1.0 - beta*(T - TRef) + ); + // kinematic turbulent thermal thermal conductivity m2/s Info<< "Reading field kappat\n" << endl; volScalarField kappat @@ -67,25 +79,41 @@ mesh ); + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("ghf", g & mesh.Cf()); + + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + p_rgh + rhok*gh + ); + label pRefCell = 0; scalar pRefValue = 0.0; setRefCell ( p, + p_rgh, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue ); - - // Kinematic density for buoyancy force - volScalarField rhok - ( - IOobject + if (p_rgh.needReference()) + { + p += dimensionedScalar ( - "rhok", - runTime.timeName(), - mesh - ), - 1.0 - beta*(T - TRef) - ); + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + } + diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H index 21be033f9b..60828e18dc 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H @@ -7,22 +7,23 @@ phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, U, phi); - surfaceScalarField buoyancyPhi = - rUAf*fvc::interpolate(rhok)*(g & mesh.Sf()); - phi += buoyancyPhi; + surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf(); + phi -= buoyancyPhi; for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pEqn + fvScalarMatrix p_rghEqn ( - fvm::laplacian(rUAf, p) == fvc::div(phi) + fvm::laplacian(rUAf, p_rgh) == fvc::div(phi) ); - pEqn.solve + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); + + p_rghEqn.solve ( mesh.solver ( - p.select + p_rgh.select ( ( finalIter @@ -36,17 +37,30 @@ if (nonOrth == nNonOrthCorr) { // Calculate the conservative fluxes - phi -= pEqn.flux(); + phi -= p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector - p.relax(); + p_rgh.relax(); // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rUAf); + U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf); U.correctBoundaryConditions(); } } #include "continuityErrs.H" + + p = p_rgh + rhok*gh; + + if (p_rgh.needReference()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - rhok*gh; + } } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H index 76cc4da8ba..477a322833 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H @@ -96,29 +96,23 @@ p_rgh + rhok*gh ); - label p_rghRefCell = 0; - scalar p_rghRefValue = 0.0; + label pRefCell = 0; + scalar pRefValue = 0.0; setRefCell ( + p, p_rgh, mesh.solutionDict().subDict("SIMPLE"), - p_rghRefCell, - p_rghRefValue + pRefCell, + pRefValue ); - scalar pRefValue = 0.0; - if (p_rgh.needReference()) { - pRefValue = readScalar - ( - mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue") - ); - p += dimensionedScalar ( "p", p.dimensions(), - pRefValue - getRefCellValue(p, p_rghRefCell) + pRefValue - getRefCellValue(p, pRefCell) ); } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H index 5664bb2154..50149ed360 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H @@ -18,17 +18,9 @@ fvm::laplacian(rUAf, p_rgh) == fvc::div(phi) ); - p_rghEqn.setReference(p_rghRefCell, p_rghRefValue); + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); - // retain the residual from the first iteration - if (nonOrth == 0) - { - p_rghEqn.solve(); - } - else - { - p_rghEqn.solve(); - } + p_rghEqn.solve(); if (nonOrth == nNonOrthCorr) { @@ -55,7 +47,8 @@ ( "p", p.dimensions(), - pRefValue - getRefCellValue(p, p_rghRefCell) + pRefValue - getRefCellValue(p, pRefCell) ); + p_rgh = p - rhok*gh; } } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H index d4878d063d..8c6a3f7671 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H @@ -17,8 +17,11 @@ == fvc::reconstruct ( - fvc::interpolate(rho)*(g & mesh.Sf()) - - fvc::snGrad(p)*mesh.magSf() - ) + ( + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ), + mesh.solver(U.select(finalIter)) ); } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C index 0075c1e3ed..cb4c4d34f6 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C @@ -80,7 +80,7 @@ int main(int argc, char *argv[]) if (nOuterCorr != 1) { - p.storePrevIter(); + p_rgh.storePrevIter(); } #include "UEqn.H" diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H index b8ac5595e4..e39d0bb803 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H @@ -53,15 +53,30 @@ ) ); + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("ghf", g & mesh.Cf()); + + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh + ( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + // Force p_rgh to be consistent with p + p_rgh = p - rho*gh; + Info<< "Creating field DpDt\n" << endl; volScalarField DpDt ( "DpDt", fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p) ); - - thermo.correct(); - - dimensionedScalar initialMass = fvc::domainIntegrate(rho); - - dimensionedScalar totalVolume = sum(mesh.V()); diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H index 3125cc3ffa..94537508b3 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H @@ -9,7 +9,7 @@ ); hEqn.relax(); - hEqn.solve(); + hEqn.solve(mesh.solver(h.select(finalIter))); thermo.correct(); } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H index a1c3dbfed5..e625f122a3 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H @@ -1,11 +1,9 @@ { - bool closedVolume = p.needReference(); - rho = thermo.rho(); // Thermodynamic density needs to be updated by psi*d(p) after the // pressure solution - done in 2 parts. Part 1: - thermo.rho() -= psi*p; + thermo.rho() -= psi*p_rgh; volScalarField rUA = 1.0/UEqn.A(); surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); @@ -18,24 +16,23 @@ + fvc::ddtPhiCorr(rUA, rho, U, phi) ); - surfaceScalarField buoyancyPhi = - rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf()); + surfaceScalarField buoyancyPhi = -rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf(); phi += buoyancyPhi; for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pEqn + fvScalarMatrix p_rghEqn ( - fvc::ddt(rho) + psi*correction(fvm::ddt(p)) + fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) + fvc::div(phi) - - fvm::laplacian(rhorUAf, p) + - fvm::laplacian(rhorUAf, p_rgh) ); - pEqn.solve + p_rghEqn.solve ( mesh.solver ( - p.select + p_rgh.select ( ( finalIter @@ -49,34 +46,25 @@ if (nonOrth == nNonOrthCorr) { // Calculate the conservative fluxes - phi += pEqn.flux(); + phi += p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector - p.relax(); + p_rgh.relax(); // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U += rUA*fvc::reconstruct((buoyancyPhi + pEqn.flux())/rhorUAf); + U += rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf); U.correctBoundaryConditions(); } } + p = p_rgh + rho*gh; + // Second part of thermodynamic density update - thermo.rho() += psi*p; + thermo.rho() += psi*p_rgh; DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); #include "rhoEqn.H" #include "compressibleContinuityErrs.H" - - // 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); - thermo.rho() = psi*p; - rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume; - } } diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C index 06fa5d12a0..64cc110cec 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C @@ -62,10 +62,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "hEqn.H" - for (int i=0; i<3; i++) - { #include "pEqn.H" - } } turbulence->correct(); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H index cf23150332..d6fa9acee9 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H @@ -23,20 +23,6 @@ volScalarField& h = thermo.h(); const volScalarField& psi = thermo.psi(); - Info<< "Reading field p_rgh\n" << endl; - volScalarField p_rgh - ( - IOobject - ( - "p_rgh", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -53,7 +39,6 @@ #include "compressibleCreatePhi.H" - Info<< "Creating turbulence model\n" << endl; autoPtr turbulence ( @@ -66,40 +51,39 @@ ) ); + Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); surfaceScalarField ghf("ghf", g & mesh.Cf()); - p = p_rgh + rho*gh; - thermo.correct(); - rho = thermo.rho(); - p_rgh = p - rho*gh; - - label p_rghRefCell = 0; - scalar p_rghRefValue = 0.0; - setRefCell + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh ( - p_rgh, - mesh.solutionDict().subDict("SIMPLE"), - p_rghRefCell, - p_rghRefValue + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh ); + // Force p_rgh to be consistent with p + p_rgh = p - rho*gh; + + + label pRefCell = 0; scalar pRefValue = 0.0; - - if (p_rgh.needReference()) - { - pRefValue = readScalar - ( - mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue") - ); - - p += dimensionedScalar - ( - "p", - p.dimensions(), - pRefValue - getRefCellValue(p, p_rghRefCell) - ); - } + setRefCell + ( + p, + p_rgh, + mesh.solutionDict().subDict("SIMPLE"), + pRefCell, + pRefValue + ); dimensionedScalar initialMass = fvc::domainIntegrate(rho); + dimensionedScalar totalVolume = sum(mesh.V()); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H index 3088b42c69..f1a62d4770 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H @@ -1,11 +1,12 @@ { rho = thermo.rho(); + rho.relax(); volScalarField rUA = 1.0/UEqn().A(); surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); U = rUA*UEqn().H(); - //UEqn.clear(); + UEqn.clear(); phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); bool closedVolume = adjustPhi(phi, U, p_rgh); @@ -20,7 +21,7 @@ fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi) ); - p_rghEqn.setReference(p_rghRefCell, p_rghRefValue); + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.solve(); if (nonOrth == nNonOrthCorr) @@ -42,13 +43,13 @@ p = p_rgh + rho*gh; - // For closed-volume cases adjust the pressure and density levels + // For closed-volume cases adjust the pressure level // to obey overall mass continuity if (closedVolume) { p += (initialMass - fvc::domainIntegrate(psi*p)) /fvc::domainIntegrate(psi); - p_rgh == p - rho*gh; + p_rgh = p - rho*gh; } rho = thermo.rho(); diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C index 0dbe80c67c..6c35536333 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C @@ -58,7 +58,7 @@ int main(int argc, char *argv[]) #include "readSIMPLEControls.H" - p.storePrevIter(); + p_rgh.storePrevIter(); rho.storePrevIter(); // Pressure-velocity SIMPLE corrector diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C index a81b0faaf3..36b3f1c3b0 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C @@ -93,6 +93,8 @@ int main(int argc, char *argv[]) // --- PIMPLE loop for (int oCorr=0; oCorr phiFluid(fluidRegions.size()); PtrList gFluid(fluidRegions.size()); PtrList turbulence(fluidRegions.size()); - PtrList DpDtf(fluidRegions.size()); + PtrList p_rghFluid(fluidRegions.size()); + PtrList ghFluid(fluidRegions.size()); + PtrList ghfFluid(fluidRegions.size()); List initialMassFluid(fluidRegions.size()); List