[SadBoySquad]

Hello!

I am trying to simulate a buoyant flow inside a building.
The buoyant currents are created by heatFluxes or higher than ambient temperatures on some of the walls.

Realistic configurations most of the time require heat fluxes, however fixed temperatures can also be used for testing purposes.
The building has various openings that the flow can enter and exit from depending on the currents.

The premise

1. boundary condition setup

For openings, I am using the following boundary conditions:

• p_rgh is prghTotalHydrostaticPressure
• U is pressureInletOutletVelocity
• k|epsilon|omega|T are inletOutlet

The rest of the boundaries are solid walls.

Some of the walls are adiabatic whereas others may have a fixedValue or externalWallHeatFlux for T.

2. initializing hydrostatic pressure

Since I am using the prghTotalHydrostaticPressure boundary condition, I need to initialize p_rgh properly.
I achieve this by utilizing the hydrostaticPressure function.

In my controlDict file, I include:

    initPrgh
    {   
        type            hydrostaticPressure;
        libs            (initialisationFunctionObjects);
        pRef            none;
    }

3. buoyantPimpleFoam with perfectGas equation of state works well and gives reasonable results

I have successfully simulated the flow with buoyantPimpleFoam and perfectGas equation of state.

Hydrostatic pressure is initialized correctly.
Flow enters from lower openings and exits from upper openings (since buoyant currents drive the flow upwards).
So far so good.

The problems

I want to simulate the exact same problem with boussinesq approximation, in order to compare the results with a different solver.

1. buoyantPimpleFoam with boussinesq equation of state crashes on the first iteration

As far as I can tell, p_rgh is not initialized correctly by the hydrostaticPressure function; p_rgh in 0 time has a uniform 1e05 value.
In the perfectGas solution, p_rgh did not have a uniform value. It was correctly initialized by the hydrostaticPressure function.

The only thing I change is the relevant entry in the thermophysicalProperties dictionary and also adding the equationOfState object with rho0, T0, and beta.
Temperature BCs don't matter; I have even tested with 1K temperature difference between ambient temperature and fixedValue in "hot" walls.

Solver always crashes on first iteration.

Is there anything else that I need to change?

2. buoyantBoussinesqPimpleFoam does not crash but also does not initialize p_rgh correctly, leading to wrong results.

When trying to simulate the same problem with buoyantBoussinesqPimpleFoam, I am not able to utilize the hydrostaticPressure function.

I get the following warnings:

--> FOAM Warning : 

    failed lookup of thermophysicalProperties (objectRegistry region0)
    available objects of type basicThermo:
0()


    From const Type& Foam::objectRegistry::lookupObject(const Foam::word&, bool) const [with Type = Foam::basicThermo]
    in file ./src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 571.
--> loading function object 'initPrgh'

The solver does not crash, but the resulting flow is not correct.
This is expected, as p_rgh at the start of the simulation does not have the correct values for hydrostatic equilibrium.

Is the hydrostaticPressure function available in buoyantBoussinesqPimpleFoam?
What should I do in order to initialize p_rgh correctly?