Wrong implementation of the Boussinesq assumption in heat transfer solvers
 

Dear FOAMers,

I don't know if this issue has ever been discussed before, but nevertheless I feel that it's important to raise it up, hoping that core OF developers can have a look at the matter.

According to the release notes of Version 7 (https://openfoam.org/release/7/) which I only read today, "The Boussinesq equation of state can now be applied to any buoyancy solver, deprecating the specialized buoyantBoussinesq[SP]impleFoam solvers".

As a convinced (and stubborn) user of the buoyantBoussinesq[SP]impleFoam solvers (and of the conjugateHeatFoam in the foam-extend version), and as a researcher and DNS code developer dealing with buoyancy-induced flows since 15 years, I am quite disappointed by the news, since, in my opinion, this corresponds to the end of the possibility of simulating buoyant flows with a correct implementation of the Boussinesq assumption in OpenFOAM.

To put it shortly:
- the Boussinesq assumption is not an equation of state
- in the Boussinesq assumption, density varies linearly with temperature only in the buoyancy term
- the usefulness of the Boussinesq assumption from a numerical point of view is the possiblility to exploit the advantages of the incompressible framework to simulate what, as a matter of fact, is a low-Mach number compressibility effect, using a gauge pressure to avoid numerical instability

In the OpenFOAM implementation (please do correct me if I am wrong):
- rho = rho0*(1 - beta*(T - T0)) is treated as an equation of state
- this is applied to all terms in the conservation equations
- now all heat transfer solvers are based on a compressible framework, and pressure is absolute. This could lead to ill-conditioned linear systems when, e.g. free convection problems, especially with open boundaries or with big, closed domains mimicking open spaces (think of, say, flow originating from an unconfined heated wire), are studied.

I am aware of the fact that there was a need to rationalize the heat transfer compartment in OpenFOAM, but I suggest finding a way to make the Boussinesq assumption a variant of the rhoConst equationOfState, only affecting the buoyancy term in the momentum equation, rather than an equationOfState as it is now. This would not only be more correct and coherent with the theory, but it also would restore the exact implementation of the buoyantBoussinesq[SP]impleFoam solvers.

One more remark: the "BernardCells" of the tutorial are misspelt. The phenomenon is known worldwide as Rayleigh-Bénard convection, the second surname being that of the French physicist Henri Bénard. I believe that a tutorial of such a fantastic code like OpenFOAM should have correct spelling of scientists' names.

Thanks a lot by now for anyone's attention, and replies/remarks

Diego

Wrong implementation of the Boussinesq assumption in heat transfer solvers
 

Dear All,
I was slightly surprised by the name of folder BernardCells when testing buoyantPimpleFoam. After visualizing the calculation results in Paraview, I realized that this is Bénard – Rayleigh convection in Bénard cells!

Kerim

Hello Diego!

To my observation, buoyantBoussinesq*Foams are still present in v1912. Please correct me if I am wrong.

Dear Andrea,

which cases are you targeting?

In the Boussinesq domain I would stick with the legacy solvers if you have to do fluid-only simulations, while I'd use the foam-extend conjugateHeat(Simple)Foam solvers for conjugate heat transfer.

Best regards

Diego

HPE, thanks for the info!

This could surely be useful to "patch back" the Boussinesq solvers in the .org version. However, my point is more conceptual about what is and what is not the Boussinesq approximation... and this seems to be ignored by the developers of the biggest open source continuum mechanics library.

Hi,

Might the org developers dont know the issue or this thread? I would suggest to submit an issue ticket in .org-Mantis, otherwise.

And/or in .com-GitLab even though the solvers exist there.

Thank you for your contributions.

Dear Andrea,

I have been trying to get a converged solution for natural convection in square cavity case using buoyantSimpleFoam. This should be one of the simplest problems however, as you said, I am unsure about the default convergence criteria set in the solver. I would be grateful to have a discussion here. Apart from the convergence criteria, I would like to know if using buoyantPimpleFoam is a better choice for steady state problems. True that piso and pimple are transient solvers but, how bad is using that for steady state? Are there any other issues you have noticed in v7 of OpenFOAM specifically in buoyancy related solvers?

FYI, I am new to OF and have been using it for a year :)

Dear Diego and Kerim, request you too to join.

Thank you

Dear damu1414,

I'm not sure I'm the best one around here to give general advice but I'll try to give you some hints that could be helpful with the discussion.
Also, I'm not so sure this thread is the most appropriate one for your questions, anyway... ;)

I do not have great experience with density-based solvers, as after my initial troubles I switched back to OF6 to be able to use the Boussinesq-type ones (since my problem at hand hardly needed compressible fluid modelling).
Generally speaking, since you are talking about convergence it should be useful to include your case setup, and somewhat more precise info on results (maybe residuals etc).
If everything is set "correctly", or without further information, the only advice I can give you is to try to play with relaxation factors and see what appens.

I can't see any particular reason why you shouldn't at all use transient solvers for steady state calculations, apart from computational efficiency and some slightly more complex case setup.
As a matter of fact, in my (very limited) experience some problems can be quite tough to be brought to convergence with steady state algorithms.
Transient solvers may be helpful in such cases as stability can be enforced on more "physical" grounds (by means of Courant control etc) instead of a bare iterative procedure and somewhat arbitrary relaxation factors.

As I said before, my adventures with OF7 and buoyant solvers came to an end quite quickly, so I'm afraid I'm of no help with regard to this :)

Best,

Andrea

Dear Andrea

Many thanks for the response. Please see attached the case files. This is for a square cavity with left and right walls kept at hot and cold temperatures respectively and top and bottom being adiabatic. I have used the numericals to match Ra=10^4.
Physical properties of air taken at 1atm and 100degC. I have also included the mesh file(done in gambit).

FYI, I am trying to validate my results(Nusselt number to be specific) with the literature NATURAL CONVECTION OF AIR IN A SQUARE CAVITY : A BENCH MARK NUMERICAL SOLUTION by G. DE VAHL DAVIS

I would request you to have a look and please let me know if any further info is required.

Thank you

Two quick preliminary considerations:

1) Why didn't you use blockMesh for mesh generation? Your geometry is so simple and already available in lots of OpenFOAM tutorials

2) Did you notice that under the buoyantSimpleFoam tutorials folder there is a buoyantCavity tutorial which is very similar to your case? It seems to run fine

Andrea

Dear Andrea,

I was looking for an update before responding to you. And yes, there is.

I switched to OF5(I had to downgrade to Ubuntu 18.04 though) and applied the buoyantBoussinesqPimpleFoam solver for my case. To my surprise, it gave me results matching with the existing literature. As discussed earlier in this thread, it seems to be an issue with OF7. Am not sure.

To your question of why am not using blockMesh. I have been using gambit for quite sometime and some old habits die hard. :). Also, am new to OF and learning in-house meshing tools.

Thanks a ton for you and Diego once again.

I am interested in the original intent of this thread. Any update on reverting this change for Boussinesq solvers? I have certainly observed challenges with instability using the new solver that I never saw previously.

Sounds like the answer is "no" - no change. Possible workarounds are:
1. use rhoConst EOS in the solver (with rho = rho0) , and then add your own momentum source term for the Boussinesq gravity term (rho0*beta*(T - T0))
2. or, build you own Boussinesq solver from [sp]impleFoam, adding in the temperature equation and the Boussinesq gravity term.

Option 2 would be cleaner, since you avoid having to deal with p_rgh, but is obviously a little bit more effort. You can reuse much of the coding from buoyant[SP]impleFoam however, for the T equation. Good luck.

You can rather easily port the boussinesq solvers from OF6 to newer code releases. Some parts of the OF libraries needed by the legacy boussinesq solvers have been moved/renamed, but AFAIK all the functionality is still present somewhere and the exact behaviour can be reproduced.

Andrea

Dear damu1414,
I am sorry for late joining to this discussing.
First of all there is inconsistency between blockMeshDict (3D) and the files (2D) from 0 folder in your caseFiles.

I would like to repeat your OF5 results. Could you share your right caseFiles?

In the meantime, if you're interested, please see papre, where we've made a few corrections to buoyantCavity tutorial mentioned by Duro:
https://aip.scitation.org/doi/abs/10.1063/5.0071571