
[Sponsors] 
Wrong implementation of the Boussinesq assumption in heat transfer solvers 

LinkBack  Thread Tools  Search this Thread  Display Modes 
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 foamextend version), and as a researcher and DNS code developer dealing with buoyancyinduced 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 lowMach 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 illconditioned 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 RayleighBé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 

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 

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/Wiki1/Wiki2/Code guide/Code Wiki/Journal Nilsson/Guerrero/Holzinger/Holzmann/Nagy/Santos/Nozaki/Jasak/Primer Governance Bugs/Features: OpenFOAM (ESIOpenCFDTrademark) Bugs/Features: FOAMExtend (WikkiFSB) 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:
which cases are you targeting? In the Boussinesq domain I would stick with the legacy solvers if you have to do fluidonly simulations, while I'd use the foamextend 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:
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 .orgMantis, otherwise. And/or in .comGitLab even though the solvers exist there. Thank you for your contributions.
__________________
The OpenFOAM community is the biggest contributor to OpenFOAM: User guide/Wiki1/Wiki2/Code guide/Code Wiki/Journal Nilsson/Guerrero/Holzinger/Holzmann/Nagy/Santos/Nozaki/Jasak/Primer Governance Bugs/Features: OpenFOAM (ESIOpenCFDTrademark) Bugs/Features: FOAMExtend (WikkiFSB) 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:
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:
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:
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 

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 inhouse 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:
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<volScalarField>("T"); //T field //double buoyancy = 0.001; double buoyancy= 9.81*rho0*(1.0beta*(TT0)); //Pout<< "T = "<<T<< endl; vectorField& USource = eqn.source(); USource += vector(0, buoyancy, 0); //, #}; codeAddAlphaRhoSup #{ Pout<< "**codeAddAlphaRhoSup**" << endl; #}; } But, I have the following error: error: cannot convert ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ 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.0beta*(TT0)); 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:
 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:
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" #}; 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<volScalarField>("T"); //T field volScalarField buoyancy ( IOobject ( "buoyancy", runTime.timeName(), mesh, IOobject::NO_READ, 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.0beta*(TT0)); //Pout<< "T = "<<T<< endl; USource += vector(0, buoyancy, 0); //, #}; codeAddAlphaRhoSup #{ Pout<< "**codeAddAlphaRhoSup**" << endl; #}; } 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 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Radiation heat transfer boundary condition  natrask  OpenFOAM Programming & Development  0  February 8, 2015 09:05 
Different solvers for heat transfer problems  aylalisa  OpenFOAM Running, Solving & CFD  7  March 15, 2013 04:11 
Question about heat transfer simulation  Anna Tian  Main CFD Forum  0  January 25, 2013 18:53 
Water subcooled boiling  Attesz  CFX  7  January 5, 2013 03:32 
How can I increase Heat Transfer at Domain Interf?  B.Simon  CFX  3  October 28, 2008 18:53 