[cckohrs]

I have a large domain axi-symmetrical (100 m x 100 m) atmospheric simulation with buoyant flow. I need boundary conditions on the inlet (side) and outlet (top) for which both the direction (normal to the patches, -x and +y respectively) and sign (only inlet -x on inlet, only outlet +y on outlet) are specified.

I have tried combinations of inletOutlet and outletInlet, pressureInletVelocity, pressureNormalInletOutletVelocity and pressureDirectedInletVelocity but I keep getting reverse flow in both cases. Can this be done without programming a custom BC?

Because this is atmospheric with a large domain (and based on other people's work in Fluent) reverse flow is not expected.

[Alczem]

Hey


Large domains and proper open boundary conditions for natural convection seems to be one of the recurrent issues with OpenFOAM. I have been looking for a correct setup for quite some time, and I am not the only one take a look at this thread, maybe it will give you some hints and you will figure it out (I have not yet) :


Common setup for p_rgh in complete open domains



Hopefully we will get a solution soon!

[cckohrs]

Thank you. That thread helped me a while ago to get to my current solution, which has steady and normal-looking behaviour (apart from my mentioned problem). If you are interested I can share my settings, maybe it helps you.

For my current problem: I got half a solution by reducing the domain size (both height and width), which solves the back-flow problem, but I am not fully happy with this work-around and will probably revisit it at some point.

[Alczem]

If you are interested I can share my settings, maybe it helps you.

If you can share your case it will definitely help.

I got half a solution by reducing the domain size (both height and width)

Huh this is weird, this is the opposite of the usual recommendations about setting the boundaries far from the region of interest, but hey, if it works

[cckohrs]

We can discuss what I can share, will DM you.

About the reduced domain size: I'm almost sure there is some issue with the BCs causing this and my "remedy" is just ignoring the issue while causing some unknown inaccuracies. So in the long run I will need a proper solution.

[cckohrs]

Oh never mind. The cause of my issues seem to be that all the fluid in my domain starts dropping as soon as the simulation starts (with thermal and cloud effects removed). Experiments so far show that it is dependent of molar mass (drops less with lower molar mass, but doesn't rise), and is independent of nHydrostaticCorrectors and vapour content (I'm working with humid atmospheric air).

Any suggestions?

[Alczem]

I remember having better stability with incompressiblePerfectGas. The fluid starting to "fall" during the first iterations are probably due to the density and pressure variations right after initialisation. If there are zeroGradient conditions, it might not help much to correct this behavior, not sure about this


Did you try with other equations of state?

[cckohrs]

It appears that this was due to my p_rgh boundary conditions; setting everything solid to fixedFluxPressure solves the falling problem (I was able to see this by writing relatively small timesteps, the pressures near my solid components started dropping as soon as time began and remain there, this was the clue that they were the issue rather than my inlet or outlet or indeed the equation of state). Combined with totalPressure on the outlet (top), I'm getting seemingly good behaviour there (the flow is still in reverse, like my image above, but it flow buoyant flow at least exits the domain instead of simply recirculating (a problem I've had in the past). Nothing I'm doing is, however, giving the expected inlet (side) behaviour.

I saw on one of your other posts that you were struggling to get steady-state solution, taking 5000 iterations to solve. What has worked for me is to solve with Euler until I have nearly steady flow and then changing to localEuler (this has worked better than steadyState).

[cckohrs]

So turns out specifying an hRef is very important, although I don't understand the function of hRef. I did not have one before. Now, when I set it to the elevation of the top of my domain, I get the flow exiting correctly at the top and no falling behaviour. Nearly perfect now with just some small steady movement to figure out (when I neglect the heat transfer effects).

[Alczem]

Hey


It sounds like you are making good progress!

What has worked for me is to solve with Euler until I have nearly steady flow and then changing to localEuler (this has worked better than steadyState).

Yep, I often end up doing that, but it is rather tedious especially in CHT cases with solids taking a long time to reach a steady state (I tried to set their schemes to steadyState while the fluid would remain with Euler but it was really unstable, and decreasing their specific heat so they would heat up faster was also shaky ).

So turns out specifying an hRef is very important

Aah, I could not figure out how this file would affect the simulation. I will have to try your way. How do you think it would behave with only vertical patches? Should we set hRef to the top of the domain anyway?


Could you post some pictures of your flow with a different hRef if you have some time?

[cckohrs]

Originally Posted by Alczem

Hey
Yep, I often end up doing that, but it is rather tedious especially in CHT cases with solids taking a long time to reach a steady state (I tried to set their schemes to steadyState while the fluid would remain with Euler but it was really unstable, and decreasing their specific heat so they would heat up faster was also shaky ).

OK. Haven't done CHT yet. Perhaps something that may help (not sure about solids though): I find that with buouyantReactingFoam I can increase the max courant number (maxCo) in controlDict up to 30 and maintain stability. I'm not sure why this works, but it does. Also, you did not mention localEuler; have you used this? I've found much better behaviour than steadyState.

Originally Posted by Alczem

Aah, I could not figure out how this file would affect the simulation. I will have to try your way. How do you think it would behave with only vertical patches? Should we set hRef to the top of the domain anyway?

Could you post some pictures of your flow with a different hRef if you have some time?

Not sure what its effect is exactly so it's difficult to predict what will happen with vertical patches. I will re-run my sim with hRef = 0, hRef = 50 (middle) and hRef = 100 (top) and post pictures.

[Alczem]

Also, you did not mention localEuler; have you used this?

I did with good results, there were less spurious currents near the boundaries, and the convective flow was established in around 1500-2000 iterations. It seems better suited for natural convection and I think it makes sense since localEuler is a "fast" Euler (if I am not mistaken). I had to use prghTotalHydrostaticPressure on the open boundaries and the hydrostatic initialization with hRef set to the top boundary. Weirdly, at T=0, even with the hydrostatic initialization, the pressure values at the boundaries were constant, and the p_rgh field looked inverted (with values increasing with height).


My tests were run with a maxCo = 5 in the fvSolution.PIMPLE dictionary so I expect you can ramp it up to 10 or maybe 20 depending on your case. I am not sure how to set up relaxation factors with localEuler though: do we need to relax UFinal, pFinal etc or just U, p etc?

[cckohrs]

Originally Posted by Alczem

I did with good results, there were less spurious currents near the boundaries, and the convective flow was established in around 1500-2000 iterations. It seems better suited for natural convection and I think it makes sense since localEuler is a "fast" Euler (if I am not mistaken). I had to use prghTotalHydrostaticPressure on the open boundaries and the hydrostatic initialization with hRef set to the top boundary. Weirdly, at T=0, even with the hydrostatic initialization, the pressure values at the boundaries were constant, and the p_rgh field looked inverted (with values increasing with height).

I had this issue with the inverse pressure field and I can't remember for certain how I managed to solve it. I think I may have had the p and p_rgh fields swapped. So what's working for me now is p is all calculated except for empty or wedge fields; ph_rgh is all fixedFluxPressure; p_rgh is all fixedFluxPressure except for the side (inlet) which is totalPressure. I'm starting to thing that prghTotalHydrostaticPressure and its friends should be used in cases where p and p_rgh are not split (but I havent run cases like that) so I stand to be corrected)

Originally Posted by Alczem

My tests were run with a maxCo = 5 in the fvSolution.PIMPLE dictionary so I expect you can ramp it up to 10 or maybe 20 depending on your case. I am not sure how to set up relaxation factors with localEuler though: do we need to relax UFinal, pFinal etc or just U, p etc?

If I remember correctly localEuler (or local timestepping) calculates a local Courant number and adjusts the timestep by cell to the maximum. Since it's not really steadystate, I don't know if it can be underrelaxed, with the local timestep acting like the relaxation in a sense. I also don't know how the maximum local Courant number is set i OpenFoam.

You make me doubt where to specify the maxCo. It can be specified in fvSolution or in controlDict... So which one is the governing one.

About the pictures: I ran the simulations but cfd-online won't allow me to upload them. Will try again at some point.

[Alczem]

So what's working for me now is p is all calculated except for empty or wedge fields; ph_rgh is all fixedFluxPressure; p_rgh is all fixedFluxPressure except for the side (inlet) which is totalPressure.

I usually have a similar setup.

I'm starting to thing that prghTotalHydrostaticPressure and its friends should be used in cases where p and p_rgh are not split

What do you mean by "p and p_rgh are not split"?

If I remember correctly localEuler (or local timestepping) calculates a local Courant number and adjusts the timestep by cell to the maximum. Since it's not really steadystate, I don't know if it can be underrelaxed, with the local timestep acting like the relaxation in a sense. I also don't know how the maximum local Courant number is set i OpenFoam.

You make me doubt where to specify the maxCo. It can be specified in fvSolution or in controlDict... So which one is the governing one.

That is my understanding as well about localEuler, it speeds up a transient simulation by allowing each cell to use the biggest timestep possible, so you lose the time accurate evolution of your system, but the final state is reached faster. I am pretty sure the maxCo for localEuler is the one defined in fvSolution.PIMPLE. I did several tests where I increased it and the simulation eventually crashed. But it ran fine if I changed the maxCo in controlDict. There are also more settings besides maxCo such as damping and sweep coefficients but they are a bit mysterious


Besides this weird inverted pressure field at the initialization, I am still satisfied with the hydrostatic initialization, localEuler and prghTotalHydrostaticPressure. I was previously reluctant to use localEuler solely because I thought natural convection should be achievable by a steady state simulation, but I don't think there are many other options at the moment!

[cckohrs]

Originally Posted by Alczem

I usually have a similar setup.
What do you mean by "p and p_rgh are not split"?

I thought prghTotalHydrostaticPressure provides a BC that varies with height. But we don't want p_rgh to vary with height since p_rgh=p-rgh? And p is (in my model at least) a value that varies with height (I'm not sure if this is due to Hydrostatic Initialisation in fvSolution, or due to anything I'm doing in ph_rgh, or both). Point is, I'm getting what I believe are desirable results without using prghTotalHydrostaticPressure or similar. I was thinking about a case like a dambreak but with an inlet wall, or perhaps just any undersea case where pressure variation with depth is important but the domain is part of an infinite field (there may be tutorials like that, I haven't looked). The domain wall pressure needs to vary with depth, but you maybe won't split this into p and p_rgh fields. My experience with OpenFoam is quite limited so I don't know.

Originally Posted by Alczem

That is my understanding as well about localEuler, it speeds up a transient simulation by allowing each cell to use the biggest timestep possible, so you lose the time accurate evolution of your system, but the final state is reached faster. I am pretty sure the maxCo for localEuler is the one defined in fvSolution.PIMPLE. I did several tests where I increased it and the simulation eventually crashed. But it ran fine if I changed the maxCo in controlDict. There are also more settings besides maxCo such as damping and sweep coefficients but they are a bit mysterious

Thanks for this. I'll experiment some more.

Originally Posted by Alczem

Besides this weird inverted pressure field at the initialization, I am still satisfied with the hydrostatic initialization, localEuler and prghTotalHydrostaticPressure. I was previously reluctant to use localEuler solely because I thought natural convection should be achievable by a steady state simulation, but I don't think there are many other options at the moment!

I'm still worried about the inverted field but it's probably my lack of understanding. Does p also vary with height in your simulation?

Probably a silly question but I only ask due to personal experience Is your g in the right direction.

[Alczem]

I thought prghTotalHydrostaticPressure provides a BC that varies with height. But we don't want p_rgh to vary with height since p_rgh=p-rgh? And p is (in my model at least) a value that varies with height (I'm not sure if this is due to Hydrostatic Initialisation in fvSolution, or due to anything I'm doing in ph_rgh, or both). Point is, I'm getting what I believe are desirable results without using prghTotalHydrostaticPressure or similar. I was thinking about a case like a dambreak but with an inlet wall, or perhaps just any undersea case where pressure variation with depth is important but the domain is part of an infinite field (there may be tutorials like that, I haven't looked). The domain wall pressure needs to vary with depth, but you maybe won't split this into p and p_rgh fields. My experience with OpenFoam is quite limited so I don't know.

I agree with the fact that p_rgh should be relatively constant with height. But I rarely work with domains tall enough to affect significatively the hydrostatic field. Not sure how prghTotalHydrostaticPressure would handle multiphase flows, but I think that you would need both p and p_rgh for the example you discuss (water and depth).

I'm still worried about the inverted field but it's probably my lack of understanding. Does p also vary with height in your simulation?

Probably a silly question but I only ask due to personal experience Is your g in the right direction.

p does vary with height, only in the wrong direction maybe it is due to hRef, I would have to check my latest tests. And yeah, gravity is properly set (not a silly question at all, I stopper counting how many times I made that mistake ) and the flow is eventually moving upwards as expected.


Using localEuler is definitely more stable and the end result looks great (using prghTotalHydrosaticPressure, fixedFluxPressure for ph_rgh and hRef set at the top boundary) so I will settle for that for now. The steadyState gives more oscillations and spurious currents with slower convergence.