From 3b038c37029b1e51ef18f7b4b67dfb2c1e92890e Mon Sep 17 00:00:00 2001 From: henry Date: Thu, 28 Jan 2010 16:06:08 +0000 Subject: [PATCH] settlingFoam: Changed to solving for pmh (static pressure minus hydrostatic pressure). While p and pmh (pd in OpenFOAM-1.5.?) are equivalent on orthogonal meshes they are not on non-orthogonal meshes and the difference is very important for buoyancy-dominated flows such as settling. settlingFoam is now written in terms of pmh (static pressure minus hydrostatic pressure) which used to be called pd but that confused too may people. --- .../solvers/multiphase/settlingFoam/UEqn.H | 3 +- .../multiphase/settlingFoam/createFields.H | 56 ++++++++++++++++--- .../solvers/multiphase/settlingFoam/pmhEqn.H | 52 +++++++++++++++++ .../multiphase/settlingFoam/settlingFoam.C | 30 +++++----- .../settlingFoam/ras/dahl/0/{p => pmh} | 2 +- .../ras/dahl/constant/polyMesh/boundary | 2 +- .../settlingFoam/ras/dahl/system/fvSchemes | 4 +- .../settlingFoam/ras/dahl/system/fvSolution | 2 +- .../multiphase/settlingFoam/ras/tank3D/0/U | 3 +- .../settlingFoam/ras/tank3D/0/{p => pmh} | 17 ++++-- .../settlingFoam/ras/tank3D/Allclean | 2 +- .../ras/tank3D/system/controlDict | 2 +- .../settlingFoam/ras/tank3D/system/fvSchemes | 4 +- .../settlingFoam/ras/tank3D/system/fvSolution | 12 +++- 14 files changed, 149 insertions(+), 42 deletions(-) create mode 100644 applications/solvers/multiphase/settlingFoam/pmhEqn.H rename tutorials/multiphase/settlingFoam/ras/dahl/0/{p => pmh} (98%) rename tutorials/multiphase/settlingFoam/ras/tank3D/0/{p => pmh} (86%) diff --git a/applications/solvers/multiphase/settlingFoam/UEqn.H b/applications/solvers/multiphase/settlingFoam/UEqn.H index 04e9194363..390eeb0502 100644 --- a/applications/solvers/multiphase/settlingFoam/UEqn.H +++ b/applications/solvers/multiphase/settlingFoam/UEqn.H @@ -22,8 +22,7 @@ == fvc::reconstruct ( - fvc::interpolate(rho)*(g & mesh.Sf()) - - fvc::snGrad(p)*mesh.magSf() + (- ghf*fvc::snGrad(rho) - fvc::snGrad(pmh))*mesh.magSf() ) ); } diff --git a/applications/solvers/multiphase/settlingFoam/createFields.H b/applications/solvers/multiphase/settlingFoam/createFields.H index a085aebb37..6833c90ec9 100644 --- a/applications/solvers/multiphase/settlingFoam/createFields.H +++ b/applications/solvers/multiphase/settlingFoam/createFields.H @@ -1,9 +1,9 @@ - Info<< "Reading field p\n" << endl; - volScalarField p + Info<< "Reading field pmh\n" << endl; + volScalarField pmh ( IOobject ( - "p", + "pmh", runTime.timeName(), mesh, IOobject::MUST_READ, @@ -119,11 +119,6 @@ #include "compressibleCreatePhi.H" - label pRefCell = 0; - scalar pRefValue = 0.0; - setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue); - - Info<< "Calculating field mul\n" << endl; volScalarField mul ( @@ -346,3 +341,48 @@ ), mut + mul ); + + + Info<< "Calculating field (g.h)f\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf = surfaceScalarField("gh", g & mesh.Cf()); + + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + pmh + rho*gh + ); + + label pmhRefCell = 0; + scalar pmhRefValue = 0.0; + setRefCell + ( + pmh, + mesh.solutionDict().subDict("PISO"), + pmhRefCell, + pmhRefValue + ); + + scalar pRefValue = 0.0; + + if (pmh.needReference()) + { + pRefValue = readScalar + ( + mesh.solutionDict().subDict("PISO").lookup("pRefValue") + ); + + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pmhRefCell) + ); + } diff --git a/applications/solvers/multiphase/settlingFoam/pmhEqn.H b/applications/solvers/multiphase/settlingFoam/pmhEqn.H new file mode 100644 index 0000000000..5a9b01a6f2 --- /dev/null +++ b/applications/solvers/multiphase/settlingFoam/pmhEqn.H @@ -0,0 +1,52 @@ +volScalarField rUA = 1.0/UEqn.A(); + +surfaceScalarField rUAf +( + "(rho*(1|A(U)))", + fvc::interpolate(rho)*fvc::interpolate(rUA) +); + +U = rUA*UEqn.H(); +phi = + fvc::interpolate(rho) + *( + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) + ); + +surfaceScalarField phiU("phiU", phi); +phi -= ghf*fvc::snGrad(rho)*rUAf*mesh.magSf(); + +for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) +{ + fvScalarMatrix pmhEqn + ( + fvm::laplacian(rUAf, pmh) == fvc::ddt(rho) + fvc::div(phi) + ); + + pmhEqn.setReference(pmhRefCell, pmhRefValue); + pmhEqn.solve(); + + if (nonOrth == nNonOrthCorr) + { + phi -= pmhEqn.flux(); + } +} + +p == pmh + rho*gh; + +if (pmh.needReference()) +{ + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pmhRefCell) + ); +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +U += rUA*fvc::reconstruct((phi - phiU)/rUAf); +U.correctBoundaryConditions(); diff --git a/applications/solvers/multiphase/settlingFoam/settlingFoam.C b/applications/solvers/multiphase/settlingFoam/settlingFoam.C index 942ed729be..006870828a 100644 --- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C +++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C @@ -64,22 +64,24 @@ int main(int argc, char *argv[]) #include "rhoEqn.H" - #include "calcVdj.H" - - #include "UEqn.H" - - #include "alphaEqn.H" - - #include "correctViscosity.H" - - - // --- PISO loop - for (int corr=0; corr/dev/null +rm log.* diff --git a/tutorials/multiphase/settlingFoam/ras/tank3D/system/controlDict b/tutorials/multiphase/settlingFoam/ras/tank3D/system/controlDict index cd4042ef63..b6d0e5c258 100644 --- a/tutorials/multiphase/settlingFoam/ras/tank3D/system/controlDict +++ b/tutorials/multiphase/settlingFoam/ras/tank3D/system/controlDict @@ -17,7 +17,7 @@ FoamFile application settlingFoam; -startFrom startTime; +startFrom latestTime; startTime 0; diff --git a/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSchemes b/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSchemes index 3802ef8dca..e076870135 100644 --- a/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSchemes +++ b/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSchemes @@ -39,7 +39,7 @@ laplacianSchemes { default none; laplacian(muEff,U) Gauss linear corrected; - laplacian((rho*(1|A(U))),p) Gauss linear corrected; + laplacian((rho*(1|A(U))),pmh) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; laplacian(mut,Alpha) Gauss linear corrected; @@ -58,7 +58,7 @@ snGradSchemes fluxRequired { default no; - p ; + pmh ; } diff --git a/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSolution b/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSolution index 148cdcddc5..ac27adffb4 100644 --- a/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSolution +++ b/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSolution @@ -17,7 +17,7 @@ FoamFile solvers { - p + pmh { solver PCG; preconditioner DIC; @@ -68,9 +68,17 @@ solvers PISO { - nCorrectors 2; + momentumPredictor yes; + nOuterCorrectors 1; + nCorrectors 2; nNonOrthogonalCorrectors 0; } +relaxationFactors +{ + U 1; + k 1; + epsilon 1; +} // ************************************************************************* //