View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0002132OpenFOAMBugpublic2016-06-27 16:462017-10-20 11:59
Reporterfabian_roesler Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityalways
Status closedResolutionsuspended 
PlatformLinuxOSRed HatOS Version(please specify)
Product Versiondev 
Summary0002132: Wrong inizialisation of pressure in thermal solvers
DescriptionIn all thermal solvers, p_rgh is initialized by using the p initial p field as follows:

    //Force p_rgh to be consistent with p
    p_rgh = p – rho*gh;

Most users will define p_rgh and leave p boundaries to be calculated by the solver. However, this will lead to a inverse (wrong) initialized p_rgh field with increasing pressure along the z axis.
For cases with small z dimension this wrong initialization is accounted for in the first solver iterations leading to artificially high velocities. With growing z dimension, the solver can not compensate the wrong pressure field and will diverge in the second iteration.
For most users pressure p should be forced to be consistent with p_rgh iteratively as rho might depend on p_rgh:

    //Force p to be consistent with p_rgh
    p = p_rgh + rho*gh;

In the solver fireFoam this hydrostaticInitialization is selectable via a new switch in pimpleDict. This method phrghEqn.H should be made a comprehensive method to be used with all thermal solvers:

    1 if (pimple.dict().lookupOrDefault<bool>("hydrostaticInitialization", false))
    2 {
    3 volScalarField& ph_rgh = regIOobject::store
    4 (
    5 new volScalarField
    6 (
    7 IOobject
    8 (
    9 "ph_rgh",
   10 "0",
   11 mesh,
   12 IOobject::MUST_READ,
   13 IOobject::NO_WRITE
   14 ),
   15 mesh
   16 )
   17 );
   18
   19 if (equal(runTime.value(), 0))
   20 {
   21 p = ph_rgh + rho*gh + pRef;
   22 thermo.correct();
   23 rho = thermo.rho();
   24
   25 label nCorr
   26 (
   27 pimple.dict().lookupOrDefault<label>("nHydrostaticCorrectors", 5)
   28 );
   29
   30 for (label i=0; i<nCorr; i++)
   31 {
   32 surfaceScalarField rhof("rhof", fvc::interpolate(rho));
   33
   34 surfaceScalarField phig
   35 (
   36 "phig",
   37 -rhof*ghf*fvc::snGrad(rho)*mesh.magSf()
   38 );
   39
   40 // Update the pressure BCs to ensure flux consistency
   41 constrainPressure(ph_rgh, rho, U, phig, rhof);
   42
   43 fvScalarMatrix ph_rghEqn
   44 (
   45 fvm::laplacian(rhof, ph_rgh) == fvc::div(phig)
   46 );
   47
   48 ph_rghEqn.solve();
   49
   50 p = ph_rgh + rho*gh + pRef;
   51 thermo.correct();
   52 rho = thermo.rho();
   53
   54 Info<< "Hydrostatic pressure variation "
   55 << (max(ph_rgh) - min(ph_rgh)).value() << endl;
   56 }
   57
   58 ph_rgh.write();
   59
   60 p_rgh = ph_rgh;
   61 }
   62 }
Steps To ReproduceBuild a simple 2D blockMesh case with huge z dimensions and run a thermal solver e.g. buoyantSimpleFoam, rhoPimpleFoam, ...
The case will diverge in the first iterations.
TagsNo tags attached.

Activities

fabian_roesler

2016-06-27 16:50

reporter   ~0006467

One more thing to mention. The switch for hydrostatic initialization should better be placed into thermophysicalProperties, as only those solvers are affected by this wrong initialization. However, this is a matter of taste.

henry

2016-06-27 17:02

manager   ~0006468

You refer to product version 3.0.x but the experimental hydrostatic initialization is in OpenFOAM-dev.

I do plan to expand the concept of hydrostatic initialization across all p_rgh based solvers once all the details have been worked out and tested. This will take a bit of time but in the meantime you can test the code I put into fireFoam in any of the other p_rgh solvers and report your experience.

> The switch for hydrostatic initialization should better be placed into thermophysicalProperties

I don't think this is the correct place for this initialization switch as it relates to the solution algorithm not the thermophysical properties.

fabian_roesler

2016-06-27 19:42

reporter   ~0006469

> You refer to product version 3.0.x but the experimental hydrostatic initialization is in OpenFOAM-dev.

I'm sorry, I was mixing up all the version I am using :)

> I don't think this is the correct place for this initialization switch as it relates to the solution algorithm not the thermophysical properties.

You're right on that. As I said, it's a matter of taste.

I will continue testing the code with the other solvers and will report back here as soon as I am done.
I'm glad that you plan to merge the code with the other solvers.

Cheers

Fabian

henry

2017-10-20 11:59

manager   ~0008906

Awaiting feedback.

Issue History

Date Modified Username Field Change
2016-06-27 16:46 fabian_roesler New Issue
2016-06-27 16:50 fabian_roesler Note Added: 0006467
2016-06-27 16:57 henry Priority high => normal
2016-06-27 17:02 henry Note Added: 0006468
2016-06-27 17:03 henry Product Version 3.0.x => dev
2016-06-27 19:42 fabian_roesler Note Added: 0006469
2017-10-20 11:59 henry Assigned To => henry
2017-10-20 11:59 henry Status new => closed
2017-10-20 11:59 henry Resolution open => suspended
2017-10-20 11:59 henry Note Added: 0008906