About Thermal LBM

Hello Everybody,
Really I am grateful to this forum. I am new user of LBM.I would like to study about thermal LBM. I have found physical and Lb unit relation between Length, mass, time etc. My question, is there any relation between real and LB temperature? I would be grateful to him also, if any body help me to send any sample Fortran code of thermal LBM. With all the best.
Yours Truly

A Fortran code, no. But Andrea has written a Matlab code which simulates a thermal LBM. It’s published, together with other LB samples, on the following page:


Thank you very much.
Would you please tell me is there any relation between real and LB temperature units?


umm, at least for the case of Rayleigh-Benard Convection, where you can define a Rayleigh number (Ra), you know that your Ra is proportional to the gradient of the temperature between upper and lower wall. Now, if you write the Ra number as ratio of quantities in real units on one side and the Ra given by your quantities in LB units on the other side of your equation… well, I think you should be able to find a relation between the gradient in lattice and real units.


and then (maybe)

deltaT_lu= (deltaTgalphaL^3)(nu_luk_lu)/(nuk)(gr_lualpha_lu*(ly-1)^3)

where gr_lu will depend on delta_t and delta_x of your simulation.
alpha_lu=chosen constant.
(ly-1)= number of lattice spaces in the vertical direction (ly-2 if you use bounce_back)
nu_lu= c2s*(tauNS-0.5) % kinematic viscosity in lattice units
k_lu = c2s*(tauThermal-0.5) % thermal diffusivity

I hope it helps… and that what I wrote make sense :wink:


Thanks for your reply.But may be you discuss about cavity flow.
But I would like to know about channel flow, if wall temperature is 300 K, then how can I convert it to lb units. Supoose, if my Gird size is 50x30, and length is L=0.0005 m. then delx=L/50=0.00001m . Now we easily can find U_real=(delx/delt)U_lb. where delt=delxCs/a.=time step.
If I want to put constant wall temp 300k, then how I can put it? I put directly 300K or I need to convert it to LB units, then how can I do it? As a new user of this method, I have some confused. Again thank you very much.

dear andrea
i am a new user of this forum. would you mention the source or paper that you use to write the program. it would be pleasure that you announce the description the case study(thermal boundary,…)taht is applied to this problem.

best regards.

Hi moslemi,

well… you could look at:

Int. J. Numer. Meth. Fluids 2002; 39:325–342 (DOI: 10.1002/d.337)

in the short Matlab code, instead of the D2Q4 lattice used by Guo for the advection diffusion equation, I use a D2Q5 (we add a rest particle to the D2Q4 lattice) but I would say that it does not change a lot for the “not completely correct advection diffusion scheme” implemented in the matlab code I wrote. To add a rest particle it would be necessary in the case you would like to use the correct advection diffusion scheme proposed by Jonas Latt in his Ph.D. thesis: http://www.tufts.edu/~jlatt01/preprints/jlatt-thesis.pdf

About the boundary condition. Reference? I do not know.
What I know is the macroscopic temperature (T_macro) on the boundary. I also know 4 of the 5 particle distribution functions (tin) I use for the advection diffusion lattice Boltzmann scheme.
I know that for a D2Q5 lattice T_macro= sum_i (tin) where i=1,2,3,4,5 and tin is the particle distr

Let’s write T_macro in a explicit form:
T_macro= tin_1+tin_2+tin_3+tin_4+tin_5
On the equation above now you have only one unknown. The unknown depend on which boundary you want to apply the constant temperature condition.

Look the matlab code for a sample simple case. for instance look at

tIn(5,:,ly) = Tcold-tIn(1,:,ly)-tIn(2,:,ly)-tIn(3,:,ly)-tIn(4,:,ly);
tIn(3,:,1) = Thot-tIn(1,:,1) -tIn(2,:,1) -tIn(4,:,1) -tIn(5,:,1);

I hope I was clear … let me know.

Hi mmatadu,

maybe this paper could help you

A Novel Thermal Model for the Lattice
Boltzmann Method in Incompressible Limit
Xiaoyi He, Shiyi Chen, and Gary D. Doolen

… there’s a small example that could be close to the oen you are trying to resolve.
The example is: Couette flow with a temperature gradient.

is this the problem you are interested about? If it is not… explain me a little bit better your geometry please. Do you have a buoyancy force in your problem? the hotter plate is the upper or lower one? … but anyway … take a look at the paper I cited. it could help.


Dear Andrea,

I’m playing with your thermal.m code. I’m trying to modify it in order to simulate a buoyancy-driven cavity flow. I have some problems with boundary conditions, in particular I don’t know how to impose a gradient temperature in the case of D2Q5 lattice. Could you help me?

Thank you.

Ciao Shaq,

Some people use a finite difference scheme for the adiabatic boundary condition.
See for example:

C. Shu, Y. Peng, Y.T. Chew, Simulation of natural convection in a square cavity by Taylor series expansion- and least square-based lattice Boltzmann method, Int. J. Mod. Phys. C 13 (10) (2002) 1399-1414.

but you should be aware of the special case of corners…


Y.Peng,C.Shu and Y.TChew
Simplified thermal lattice Boltzmann model for incompressible thermal flows. Physical Review E68 2003.

but I’m wondering if a normal bounce-back for the temperature distribution function (as we use to do for the fluid) could be enough…


Dear Andrea,
thank you for your perfect reply.
best regards

a comment maybe. In general when you want to impose a gradient of a macroscopic quantity (temperature or velocity or pressure, let’s say temperature T from now on) you can transform this gradient condition into a fixed temperature wall condition using the nearest neighbors. This is more convenient since in general it is not possible to impose a gradient on a lb code easily. Let’s have small example.

Say you want d_x T(x_0,y)=0, where x_0 is the position of the wall (say it is a left wall, i.e. the outward pointing vector is (-1,0)). You can do a Taylor expansion (up to second order) of the nearest neighbors and you have (i’ll omit the y coordinate)

T(x_0+1)=T(x_0)+d_xT(x_0)+1/2d_x^2 T(x_0)
d_xT(x_0)+2*d_x^2 T(x_0)

You already know T(x_0+1) and T(x_0+2). You want to impose d_x T(x_0)=0. Then with these conditions and the two expansions, you are left with


Therefore you can just apply a normal fixed temperature boundary condition with the temperature equal to the formula above and the scheme should remain of second order accuracy.

I hope this explanation is more or less clear…


I am new to LBM. I would like to solve convection heat transferr in enclsoure (cavity, cylinder) and channels. May anybody send sample code for enclosures or channels
with regards


I’m not completely sure about what I’m writing… but if you use the thermal 2D code for Rayleigh Benard convection written in matalb … and you add an external constant force at each element of your grid in the direction perpendicular to the buoyancy force, well … you should get a kind of Poiseuille flow in a heated up system. Using R-B example you would have no-slip boundary condition at the top-bottom walls and periodic boundary in the x-direction.

Well … you could even switch your heater on (to change boundary condition for the temperature from 0 to 1 for example… look I don’t know, but maybe a trigger for the temperature would be still needed ) once you reach a steady state for your Poiseuille solution… Do you see what I mean? it could be cool… maybe it is far from what you have in mind … but it could be helpful to get confident with Lattice Boltzmann and thermal models.

if you already play with OpenLB it should be pretty easy to make the small change that I proposed using the Rayleigh Benard 2D example. If you don’t use OpenLB, it could be a nice problem to start with. (If you wanna use OpenLB look at forced poiseuille flow and Rayleigh Benard example)


don’t hesitate to ask more if I wasn’t clear (it happens too often)

tovarish Wrote:


umm, at least for the case of Rayleigh-Benard


and then (maybe)


where gr_lu will depend on delta_t and delta_x of
alpha_lu=chosen constant.
(ly-1)= number of lattice spaces in the vertical
nu_lu= c2s*(tauNS-0.5) % kinematic viscosity in
k_lu = c2s*(tauThermal-0.5) % thermal

I hope it helps… and that what I wrote make


I’m a New user of LBM. my case is lid driven cavity flow (not affected by natural convection) with thermal boundary conditions.
Now that I can’t use Ra as my connection between LB and real parameters how can I input my temperatures into LB code?
I would be very thankful for any help .