# Wrong implementation of the Boussinesq assumption in heat transfer solvers

 Register Blogs Members List Search Today's Posts Mark Forums Read

 September 11, 2019, 06:23 Wrong implementation of the Boussinesq assumption in heat transfer solvers #1 Member   Diego Angeli Join Date: Mar 2009 Posts: 31 Rep Power: 17 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 kwardle, arjun, visakhmg and 14 others like this.

 September 13, 2019, 16:22 Wrong implementation of the Boussinesq assumption in heat transfer solvers #2 Senior Member   abdikerim kurbanaliev Join Date: Jun 2010 Location: Kyrgyzstan, Osh Posts: 120 Rep Power: 16 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 damu1414 and Teresa Sun like this.

 February 21, 2020, 10:09 #3 Member   Andrea Di Ronco Join Date: Nov 2016 Location: Milano, Italy Posts: 57 Rep Power: 9 Hello Diego! I'm starting to work on solvers based on the legacy buoyantBoussinesqPimpleFoam and I thought that updating them to the newest version of OF at the beginning of my work could be a good idea. After several days I realized too that the new OF solvers were based on different physics. Also, the new buoyant[SP]impleFoam tutorials present some strange convergence behaviour. I'm not new to OF, but I'm not a great expert of compressible solvers. I'm definitely supposed to work on large, possibly closed domains and heavy liquids. As an expert, do you have some specific advice or would you just recommend to switch back to older releases? Thank you very much. Andrea Last edited by Diro7; February 21, 2020 at 11:15.

 February 22, 2020, 08:54 #4 Senior Member     Herpes Free Engineer Join Date: Sep 2019 Location: The Home Under The Ground with the Lost Boys Posts: 932 Rep Power: 12 To my observation, buoyantBoussinesq*Foams are still present in v1912. Please correct me if I am wrong. __________________ The OpenFOAM community is the biggest contributor to OpenFOAM: User guide/Wiki-1/Wiki-2/Code guide/Code Wiki/Journal Nilsson/Guerrero/Holzinger/Holzmann/Nagy/Santos/Nozaki/Jasak/Primer Governance Bugs/Features: OpenFOAM (ESI-OpenCFD-Trademark) Bugs/Features: FOAM-Extend (Wikki-FSB) Bugs: OpenFOAM.org How to create a MWE New: Forkable OpenFOAM mirror

February 22, 2020, 12:11
#5
Member

Diego Angeli
Join Date: Mar 2009
Posts: 31
Rep Power: 17
Quote:
 Originally Posted by Diro7 Hello Diego! I'm starting to work on solvers based on the legacy buoyantBoussinesqPimpleFoam and I thought that updating them to the newest version of OF at the beginning of my work could be a good idea. After several days I realized too that the new OF solvers were based on different physics. Also, the new buoyant[SP]impleFoam tutorials present some strange convergence behaviour. I'm not new to OF, but I'm not a great expert of compressible solvers. I'm definitely supposed to work on large, possibly closed domains and heavy liquids. As an expert, do you have some specific advice or would you just recommend to switch back to older releases? Thank you very much. Andrea
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

February 22, 2020, 12:14
#6
Member

Diego Angeli
Join Date: Mar 2009
Posts: 31
Rep Power: 17
Quote:
 Originally Posted by HPE To my observation, buoyantBoussinesq*Foams are still present in v1912. Please correct me if I am wrong.
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.

 February 22, 2020, 15:46 #7 Senior Member     Herpes Free Engineer Join Date: Sep 2019 Location: The Home Under The Ground with the Lost Boys Posts: 932 Rep Power: 12 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. __________________ The OpenFOAM community is the biggest contributor to OpenFOAM: User guide/Wiki-1/Wiki-2/Code guide/Code Wiki/Journal Nilsson/Guerrero/Holzinger/Holzmann/Nagy/Santos/Nozaki/Jasak/Primer Governance Bugs/Features: OpenFOAM (ESI-OpenCFD-Trademark) Bugs/Features: FOAM-Extend (Wikki-FSB) Bugs: OpenFOAM.org How to create a MWE New: Forkable OpenFOAM mirror

 June 22, 2020, 10:33 #8 New Member   damu Join Date: Feb 2017 Posts: 24 Rep Power: 9 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

June 22, 2020, 12:31
#9
Member

Andrea Di Ronco
Join Date: Nov 2016
Location: Milano, Italy
Posts: 57
Rep Power: 9
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...

Quote:
 Originally Posted by damu1414 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.
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.

Quote:
 Originally Posted by damu1414 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?
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.

Quote:
 Originally Posted by damu1414 Are there any other issues you have noticed in v7 of OpenFOAM specifically in buoyancy related solvers?
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

June 23, 2020, 13:04
#10
New Member

damu
Join Date: Feb 2017
Posts: 24
Rep Power: 9
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
Attached Files
 caseFiles.zip (119.4 KB, 22 views) cavity_mesh.zip (120.6 KB, 10 views)

 June 25, 2020, 09:14 #11 Member   Andrea Di Ronco Join Date: Nov 2016 Location: Milano, Italy Posts: 57 Rep Power: 9 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

 July 9, 2020, 01:48 #12 New Member   damu Join Date: Feb 2017 Posts: 24 Rep Power: 9 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.

 August 25, 2021, 20:37 #13 Senior Member   Kent Wardle Join Date: Mar 2009 Location: Illinois, USA Posts: 219 Rep Power: 21 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.

 August 27, 2021, 09:31 #14 Senior Member   Join Date: Apr 2020 Location: UK Posts: 674 Rep Power: 14 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.

September 7, 2021, 09:00
#15
Member

Andrea Di Ronco
Join Date: Nov 2016
Location: Milano, Italy
Posts: 57
Rep Power: 9
Quote:
 Originally Posted by kwardle 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.
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

 February 7, 2022, 19:14 #16 Senior Member   abdikerim kurbanaliev Join Date: Jun 2010 Location: Kyrgyzstan, Osh Posts: 120 Rep Power: 16 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

 April 22, 2024, 17:53 codedFvModels Boussinesq buoyancy Force #17 New Member   Felipe Noh Join Date: Aug 2014 Posts: 8 Rep Power: 11 Dear FOAMers, If someone could please advise me on the following: I am implementing buoyancy force with the Boussinesq approximation in the momentum equations as a source term in the "fvModels" file. buoyancyBoussinesqForce { type coded; select all; field U; codeInclude #{ #}; codeAddSup #{ Pout<< "**codeAddSup**" << endl; #}; codeAddRhoSup #{ Pout<< "**codeAddRhoSup**" << endl; const scalar T0 = 293.5; const scalar rho0 = 1.20329; const scalar beta = 0.003407155; const volScalarField& T = mesh().lookupObject("T"); //T field //double buoyancy = -0.001; double buoyancy= -9.81*rho0*(1.0-beta*(T-T0)); //Pout<< "T = "< >’ to ‘double’ in initialization make: *** [/opt/openfoam11/wmake/rules/General/transform:26: Make/linux64GccDPInt32Opt/codedFvModelTemplate.o] Error 1 The error is in: double buoyancy= -9.81*rho0*(1.0-beta*(T-T0)); What is the correct variable type for buoyancy? I'm solving natural convection in a square cavity with OpenFOAM 11. For verification purposes, I want to solve it using the Boussinesq approximation. Thansk a lot Felipe Noh

April 23, 2024, 01:38
#18
Member

Andrea Di Ronco
Join Date: Nov 2016
Location: Milano, Italy
Posts: 57
Rep Power: 9
Quote:
 Originally Posted by fnohpat Dear FOAMers, If someone could please advise me on the following: I am implementing buoyancy force with the Boussinesq approximation in the momentum equations as a source term in the "fvModels" file. buoyancyBoussinesqForce { type coded; select all; field U; codeInclude #{ #}; codeAddSup #{ Pout<< "**codeAddSup**" << endl; #}; codeAddRhoSup #{ Pout<< "**codeAddRhoSup**" << endl; const scalar T0 = 293.5; const scalar rho0 = 1.20329; const scalar beta = 0.003407155; const volScalarField& T = mesh().lookupObject("T"); //T field //double buoyancy = -0.001; double buoyancy= -9.81*rho0*(1.0-beta*(T-T0)); //Pout<< "T = "< >’ to ‘double’ in initialization make: *** [/opt/openfoam11/wmake/rules/General/transform:26: Make/linux64GccDPInt32Opt/codedFvModelTemplate.o] Error 1 The error is in: double buoyancy= -9.81*rho0*(1.0-beta*(T-T0)); What is the correct variable type for buoyancy? I'm solving natural convection in a square cavity with OpenFOAM 11. For verification purposes, I want to solve it using the Boussinesq approximation. Thansk a lot Felipe Noh
After a very quick glance:

- Since T is a volScalarField, your best chance is to cast buoyancy as a volScalarField as well. This is exactly what the error message is suggesting since, from what I remember, volScalarField is shorthand for <Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>>. For sure, casting it as a mere double cannot work.

- Also T - T0 might be problematic since T has physical dimension and T0 is a scalar (i.e. a mere double again) with no dimensions. But I'm not sure about this, I'm not writing openfoam since quite a long time now

April 23, 2024, 13:58
#19
New Member

Felipe Noh
Join Date: Aug 2014
Posts: 8
Rep Power: 11
Quote:
 Originally Posted by Diro7 After a very quick glance: - Since T is a volScalarField, your best chance is to cast buoyancy as a volScalarField as well. This is exactly what the error message is suggesting since, from what I remember, volScalarField is shorthand for >. For sure, casting it as a mere double cannot work. - Also T - T0 might be problematic since T has physical dimension and T0 is a scalar (i.e. a mere double again) with no dimensions. But I'm not sure about this, I'm not writing openfoam since quite a long time now

Thanks Andrea:

I have made some changes but I have more errors.

buoyancyBoussinesqForce
{
type coded;
select all;

field U;

codeInclude
#{
#include "createTime.H"
//#include "createMesh.H"
#};

#{
#};

#{
const scalar T0 = 293.5;
const scalar rho0 = 1.20329;
const scalar beta = 0.003407155;

const volScalarField& T = mesh().lookupObject<volScalarField>("T"); //T field

volScalarField buoyancy
(
IOobject
(
"buoyancy",
runTime.timeName(),
mesh,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar ("buoyancy",dimensionSet(1,-2,-2,0,0,0,0),0)
);

vectorField& USource = eqn.source();

//double buoyancy = -0.001;
buoyancy= -9.81*rho0*(1.0-beta*(T-T0));
//Pout<< "T = "<<T<< endl;

USource += vector(0, buoyancy, 0); //,
#};

#{
#};
}

Could I send you the complete case for you to review?

Thanks

 April 24, 2024, 05:06 #20 Member   Andrea Di Ronco Join Date: Nov 2016 Location: Milano, Italy Posts: 57 Rep Power: 9 No, of course you can't. I don't have time to review your entire case and, as I said, I'm not working with openfoam anymore. If you need someone to do your job for you, I'm sure you'll find plenty of people here who do that for a living. Otherwise, you'll find plenty of free information in the forum or elsewhere to help you tackle your problem. Posting each iteration of your case every few hours or asking direct work from other people is not the proper way to ask for help here. Good luck, Andrea Tobermory likes this.