[Artur]

Good morning Tobi,

Yes, I use buoyantSimpleFoam and without initializing with incompatible solvers so just from uniform fields.

All the best

Artur

[rezika]

hi tobi


can you please post your buoyantSimpleFoam case, i have exactly the same error,



did you find a solution for the buoyantPimpleFOam?



thanks

[Andrea_85]

Hi All,

after have been playing with buoyantSimple/PimpleFoam with no success for a while, i would like to share my case here. I think it is pretty much related to what has been discussed before.

I am simulating wind flow around buildings. The geometry is an isolated street canyon as shown in fig1. The domain has dimension 230x60x250 m^3 and the 2 buildings 10x10x50 m^3. The mesh is gradually refined toward the building surfaces, such that the first cell height on every surface is 10 cm. I am using k-epsilon model for turbulence. At inlet, sides and top i use atmospheric boundary layer conditions (for U, k and epsilon), standard wall function for all street canyon surfaces and a rough wall function for the ground. Temperature is fixed at the inlet, sides and top to 297K, while on the leeward, windward and street surfaces temperature is set to 305K (all the other surfaces have zeroGradient). I use perfectGas (or incompressiblePerfectGas) as equation of state (OF6).

The case run pretty well using buoyantSimpleFoam for convective dominated flow (wind speed at the building height, i.e. 10m, U10 > 2 m/s). Starting from uniform fields, it typically converges in less than 1000 iterations. Problems come for buoyancy-driven flow, when I decrease the wind speed at the inlet. Results below refer to a case with U10 = 0.5 m/s.

- buoyantSimpleFoam
The simulation does not converge after 20000 iters (especially pressure).

Time = 20000

DILUPBiCGStab:  Solving for Ux, Initial residual = 7.79058e-05, Final residual = 7.50131e-06, No Iterations 1
DILUPBiCGStab:  Solving for Uy, Initial residual = 0.000732887, Final residual = 6.29222e-05, No Iterations 1
DILUPBiCGStab:  Solving for Uz, Initial residual = 0.00149779, Final residual = 0.000141864, No Iterations 1
DILUPBiCGStab:  Solving for h, Initial residual = 0.000721349, Final residual = 5.48416e-05, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 0.0160753, Final residual = 8.88376e-05, No Iterations 3
time step continuity errors : sum local = 7.16891e-08, global = 8.96698e-10, cumulative = -9.06306e-06
rho max/min : 1.16767 1.13963
DILUPBiCGStab:  Solving for epsilon, Initial residual = 3.35147e-05, Final residual = 1.10589e-06, No Iterations 1
DILUPBiCGStab:  Solving for k, Initial residual = 7.987e-05, Final residual = 2.87682e-06, No Iterations 1
ExecutionTime = 4518.52 s  ClockTime = 4532 s

I tried almost everything: refining the mesh, different turbulence model, longer domain, lower under-relaxation factors, upwind for velocity. None of them works. I suspect that the buoyancy-driven case is intrinsically unsteady and it is not easy to get a converged solution for steady-state. Or at least i did not manage to find one.


- buoyantPimpleFoam
Since the case is most probably unsteady, I tried to run it transient. My idea was to compare the time-averaged transient solution with the time-average of the not-converged steady solution. Perhaps they are not that different...

I used one of the (not-converged) steady-state solutions to initialize all fields (fig2 shows an example for velocity). However, i am not able to simulate more than a couple of seconds. Then i get very high velocity close to the leeward wall, Co increases and the simulation crashes. The velocity field at the time of the crash is shown in fig3.

Do you have any suggestions to have either a converged steady-state or a stable transient simulation?


ps The only differences between the buoyantSimpleFoam case (that is attached) and the buoyantPimpleFoam case are controlDict and fvSolution, that are reported below. I have also tried to use PIMPLE as suggested in this thread but with no success.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.2                               |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs ( "libOpenFOAM.so" "libgroovyBC.so")

application     buoyantPimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         100;

deltaT          0.001;

writeControl    adjustableRunTime;

writeInterval   2;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression on;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

adjustTimeStep  yes;

maxCo 0.5;

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.2                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "rho.*"
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       0;
        relTol          0;
    }
    p_rgh
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;

        smoother         DICGaussSeidel;

        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator     faceAreaPair;
        mergeLevels      1;
    }
    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    "(U|k|epsilon|h|e)"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.1;
    }

    "(U|h|e|k|epsilon|R)Final"
    {
        $U;
        relTol          0;
    }

}



PIMPLE
{
    nNonOrthogonalCorrectors 1;
    rhoMax          rhoMax [ 1 -3 0 0 0 ] 1.5;
    rhoMin          rhoMin [ 1 -3 0 0 0 ] 0.9;

    momentumPredictor yes;
    nOuterCorrectors 1;
    nCorrectors 2;

    residualControl
    {
        U      1e-4;
        p_rgh  1e-4;
        h      1e-4;
        "(k|epsilon)" 1e-4;
    }
}

/*
relaxationFactors
{
    fields
    {
      p_rgh 0.7;
      rho   0.7;
      "(p_rgh|rho)Final" 1;
    }
    equations
    {
        "(T|k|epsilon)" 0.7;
        U 0.7;
        "(h|e)" 0.7;
        "(h|e|U|T|k|epsilon)Final" 1;
    }
}
*/

Thanks,
Andrea

[Jost K]

Hi All,

since this thread has been dealing with pressure problems in buoyantPimpleFoam a lot, I would like to post this here:

I am working on a buoyantPimpleFoam case with a large water column (150m in vertical direction). Since I was having some trouble with boundary conditions I chose to set up a test case with no convection:


My domain is a simple cuboid (20x20x150m).
My fluid is water.

I initialize the temperature to a given profile.
All the walls are slip walls with fixedFluxPressure for p_rgh and calculated for p (much like in the tutorials) except for the top wall which is a pressureInletOutlet patch with totalPressure = 0 for p_rgh and calculated for p (to simulate a freeSurface / openTop).
I use polynomials for the thermophysical properties and simulate sensibleInternalEnergy.



Now, when I run the case for a couple of seconds I get some flow through the upper boundary while the hydrostatic pressure gradient adjusts (which messes up my initial temperature profile but that is only a side note here.) and then the fluid settles (U=0) with the expected p gradient (from 0 at the top to roughly 1.5e6 at the bottom) and a p_rgh gradient from 0 at the top (z=0) to -1300 at the bottom (z=-150).
This p_rgh gradient confuses me. Am I wrong to expect a uniform p_rgh field when the fluid is at rest?

I'm assuming there is something seriously wrong about the pressure in this solver.

Can anyone explain this?


Thanks,
Jost

[Artur]

Hi Jost,

That's an interesting problem you're working on. May I ask what is the final application?

As for the issue you're having, I'd also expect p_rgh to be zero throughout the water column for a fluid at rest. That being said, are you sure it is not moving at all? The value of 1300 Pa is indeed very small over the distance of 150 m, so it could be due to a residual velocity field that's barely perceptible. One thing to try would be to add a line monitor along the centreline of the domain and confirm the velocity distribution, as well as run the simulation for a longer time (e.g. twice what you have done so far) and see if the pressure gradient changes at least a little bit.

Cheers,

A

[Jost K]

Hi Artur,
Thanks for your quick reply!

My final application will be oceanic nutrient transport and pumping cold depth water to the ocean surface region (There is a discussion on some of the changes I made to buoyantPimpleFoam here( Buoyant environmental flow including salinity in OpenFOAM ). [I made sure the this problem here has nothing to do with those modifications]).

I found the source of my mysterious p_rgh gradient.
It has nothing to do with any velocities and it persists even after simulating for as long as 20000s.
In buoyantPimpleFoam (and the other heat transfer solvers) p is calculated from p_rgh as

(1)

 p=p_rgh+rho*gh

However, rho*gh only represents the hydrostatic pressure if rho is uniform (i.e. Eqn. (1) assumes rho in all cells above an arbitrary cell to be the same as in this cell).
In my case there is a vertical rho gradient which corresponds to the the temperature gradient I initialize (using funkySetFiels).

In this nonuniform case the hydrostatic pressure has to be calculated similar to this:

(2)


(with h being the vertical coordinate direction).
I first thought the assumption of uniform density was fine, because we only use the pressure gradient in our conservation equations, not the absolute value.
But for some reason the difference between the two above approximations (1) and (2) seems to be in p_rgh.
This can easily be shown though hand (excel) calculations in my case.

I cannot quite get my head around why I get this hydrostatic component in p_rgh but maybe someone can give me a hint here.
Anyway, I do not think this problem introduces large errors in the solution (because we only use the pressure gradient) but it certainly makes it impossible to formulate totalPressure or fixedValue boundary conditions for prgh.
I believe this also connects to Tobis problem, when he first started this thread and certainly to this problem ( Correct boundary conditions for p_rgh (special for vertical patches) ).

Last but not least:
has anyone got an idea for a smart way to implement the correct hydrostatic pressure calculation in BuoyantPimpeFoam?
Any function I can use for the integration?

Cheers,
Jost

[Artur]

Hi,

Thanks, that sounds like a very cool problem indeed. I'll gladly have a read through your post about the overall idea for the project.

I'm no expert in this type of problems, but your hypothesis about the origin of the problem and a possible solution you suggested both sound very plausible to me.

From the top of my head, I cannot think of a generic, straightforward way to implement an integrator for the weight of the fluid that would work robustly on complex grids and for a wide variety of problems.

That being said, the problem you are trying to solve should be much more straightforward to deal with since you're working with a water column. Hence, you could invoke a 2D assumption and integrate the weight of the fluid along each horizontal plane of the grid first, and then integrate those values moving from top to bottom of the domain. This should be reasonably accurate, provided density does not change too much in the horizontal plane.

Alternatively, if you are using structured grids without any nested refinement and your density does vary significantly in the horizontal plane, you could do the integration for each vertical "strip" of cells. This would be a bit more tricky to do in parallel, but still doable.

A third approach I can think of is the most crude but also the simplest. You could run a benchmark simulation first, then sample the solution using a uniform grid of points, do the integration externally to OpenFOAM (e.g. Python) and come up with a hydrostatic pressure gradient that way. This could then be applied on the vertical walls using either a nonUniformFixedValue BC or something like swak4foam.

Do any of those sound acceptable to you?

All the best,

A

[Tobi]

Hi all,

as I am investigating the first time deeply into the pressure calculation in OpenFOAM, I found some mistakes I always had in my mind (and therefore, I always made some mistakes in the pressure boundary conditions).

First of all. The part is not the hydrostatic pressure as I was always thinking it is. The main problem I had was the following:

• I thought that the p_rgh field is the static pressure minus the hydrostatic part
• However, actally it is but not in the sense I always thought. The hydrostatic part is always the hight difference from the reference height to the cell center (commonly to the height h = 0 m). Based on the reference height, the field can look complete different and can achieve pressure values that do not match with analytical calculations such as: at position (x) the p_rgh value has to be p - rho g dh. It only does make sense, if the reference hight in the calculation is set properly. However, for the p field itself, it does n ot matter how we set the reference height.
• I always initiat the fields using a transient solver. After a lot of tries, I realized that for such problems, the steady-state solver might go into the wrong direction (even though we have correct boundary conditions).

@Jost. Where is the difference between your equation (1) and (2)? In Foam you have the density based on the temperature and if the temperature changes with respect to the height, it is included.

[Jost K]

Hi All,

first of all I will try to clarify my point:


@Tobi:
I was a little lazy when writing eqn (2).
In my view when using variable density, the hydrostatic pressure at a certain vertical position should be calculated:

This is using the actual value of rho at each vertical position above the position . The way I understand it, calculating the way buoyantPimpleFoam does (along the lines of eqn. (1) from my last post) would result in something like this:

Which means assuming that rho at each vertical position above is equal to rho at the position .
Would you agree that that this is the way buoyantPimpleFoam calculates the hydrostatic pressure?

@Artur:
Thanks again for your quick reply.
I think you are probably right in saying that there is no easy generic way to solve this issue.
I also agree that for my problem there are plenty of reasonable workarounds. I will definitely have a closer look at all of the options you outlined. I am pretty sure they will work fine.


However, I would still like to put a big ! here, since this an issue that may well have a general influence on solutions obtained with buoyantPimpleFoam, depending on the case. Do you agree?

Cheers,
Jost

[Tobi]

Hi Jost,


I agree, FOAM calculates the field p_rgh using:

 p=p_rgh+rho*gh

but the density field is a volume scalar field and thus the quantity can be different in each numerical cell (based on the temperature and the equations you are using for the density). Thus, rho is not constant!

[Artur]

Well, it's not constant, but then again gh is:

https://develop.openfoam.com/Develop...l/include/gh.H

so one is really solving eq 1 and not eq 2 from Jost's post, it seems (i.e. ignoring the integral). I suppose this works fine e.g. for a free-surface flow, where there's a well-defined reference location and the two fluids more-or-less do not mix. But for the water column problem described here, I'm not so sure this is an appropriate way to deal with the issue.

[Jost K]

Hi Tobi,

I am sorry if my point does not come about clearly.
Here is a little example:
Let's consider the uniform grid in the attached figure,
grid.png
where we want to calculate the hydrostatic pressure at point P.

When doing so along the lines of

p=p_rgh + rho*gh

you would end up with


(just try to imagine what the rho and gh field would look like in the above case)

However, I am sure you'd agree that the actual weight of the fluid above the point P which causes the hydrostatic pressure depends on the density in the cells above P so that the calculation should be something similar to:



does this make my problem clearer?

(I apologize for my weird pseudo code, I don't know a better way to put this)

Cheers,
Jost

[Tobi]

Hi Joes,

got it. Hmm... well the p field should have these behavior, shouldn´t it? Interesting point and question that I did not investigate till now. I agree to both of you. For uniform density fields it is clear. For non-uniform-height dependend density field I am not sure right now. A simple test case should tell us how it works.

However, out of the box. If you both are right, this would mean that the static pressure field would be wrong always if there is a temperature gradient in the fluid (lets consider a large one) and this should introduce instability. Hmmm... At the end, this would mean that the physics are not correct.

[Tobi]

Hi all,

I just made a test case and based on that, the hydrostatic part is based on equation (2). Otherwise we should have the same static pressure in both cases, right.

In my case, I just changed the temperature at different locations and thus, the density changed. Do I miss something here?

[Jost K]

Hi Tobi,

the hydrostatic pressure would not be same in both cases if calculation was based on Eqn (1).

Note that the hydrostatic pressure in the areas with lower density (displayed in green and red in the right subfigure of your first figure) would still be calculated based on the value of the density there.

However, the hydrostatic pressure should display some discontinuous behavior (which will probably vanish over time) when you plot it along a vertical line.

The line should look somewhat like this:

\_
   \_
      \

if you know that I mean (density is on the horizontal axis and depth on the vertical one)


Anyway, you are right in the sense that the maximum hydrostatic pressure should be the same in both plots if the calculation was based on (1).
What exactly do you plot as the hydrostatic pressure, is it p - p_rgh or just p?

[Jost K]

Originally Posted by Jost K

What exactly do you plot as the hydrostatic pressure, is it p - p_rgh or just p?

This is really the important question here.
Because, as I mentioned earlier, in my case I get a p_rgh gradient such that p looks like it had been calculated according to Eqn (2) (which it really isn't from my understanding of the code).
(Note that p includes p_rgh.)
It's this unphysical p_rgh gradient that worries me not the fact that p might be wrong.

Did you have a look at p_rhg in your case?

Cheers,
Jost

[Muhammad Omer Mughal]

Dear Andrea

I am not sure if you have found the solution to your problem. However, I am having a similar problem in the flow around buildings. I am also using buoyantSimpleFOAM but my building geometry is a bit more complex than yours. My problem is that the flow field looks weird at a location far from (5H) from the building location. Initially the flow field looks fine but as the time is increased it looks like there might be an upward motion somewhere due to buoyancy. It could be that, at night, there is forced convection at lower surface temperatures. And during daytime, when surface temperatures increase, it becomes buoyant flow. I am attaching my flow fields for time steps 20000,30000,40000,50000,60000,70000. I would appreciate your help

@Tobi it would be great if you can comment as well

https://pasteboard.co/J2ddsiB.png
https://pasteboard.co/J2ddPnr.png
https://pasteboard.co/J2debDG.png
https://pasteboard.co/J2deotd.png
https://pasteboard.co/J2deCtM.png
https://pasteboard.co/J2deQxZ.png

[nipinl]

Originally Posted by Tobi

Hi all,

I just made a test case and based on that, the hydrostatic part is based on equation (2). Otherwise we should have the same static pressure in both cases, right.

In my case, I just changed the temperature at different locations and thus, the density changed. Do I miss something here?

This thread was a nice read. Any conclusion made regarding the calculation of p? Is that really based on equation (2) ? which I hope it should be.

[Jost K]

Hi Nipin,

The pressure is calculated from the pressure equation, which is found by combining the continuity and momentum equations.
For buoyantPimleFoam, a very good description of the pressure Equation is given here (just use google translate in case you don't speak Chinese).
For steady state flow with zero velocity the pressure equation boils down to equality of the explicit source terms found in the momentum equation.
Namely, the buoyancy and pressure gradient terms.
For a general momentum equation this means

which defines p as the hydrostatic pressure, as given in Eqn (2) in my post above.
In the case of a p_rgh solver (like buoyantSimpleFoam and buoyantPimpleFoam) the equality of the source terms is

This clearly shows that p_rgh has a hydrostatic component in flows with nonuniform density ().

This has a couple of consequences for simulations of flows with nonuniform density, most importantly:
1) Do not use zeroGradient as pressure boundary condition on horizontal walls. This can lead to large continuity errors. Use fixedFluxPressure instead.
2) Do not use fixedValue as pressure boundary condition on vertical outlets. This will lead unwanted flow, which is driven by the hydrostatic component of p_rgh. fireFoam uses the hydostatic background pressure ph_rgh, which can be applied as fixed value on pressure outlet boundaries. This seems to be the way to go for vertical outlet boundaries.


Cheers,
Jost

[Tobi]

The ph_rgh boundary condition however, can only be used with fireFoam as this field is not present for, e.g., buoyant*Foam. Hence, my question is still what we should do for vertical outlets? As I already pointed out somewhere in this forum, we probably need a new boundary condition (similar to the ph_rgh guy from firefoam) with other solvers.