File Exchange

image thumbnail

Lattice Boltzmann immiscible two-phase model (LBM)

version 1.7 (15.7 KB) by

Starting point for newcomers wanting to study RK-type lattice Boltzmann immiscible two-phase model.



View License

This code may be use as a starting point for newcomers wanting to study RK-type lattice Boltzmann immiscible two-phase model. It may give a better understanding of the RK-type lattice Boltzmann two-phase model.
This m-code reproduces some of the results of Table 3 in section 3.2 "Steady Bubble" of the following reference:
R1) Leclaire, S., Reggio, M. and Trépanier, J.-Y. (2012). Numerical evaluation of two recoloring operators for an immiscible two-phase flow lattice Boltzmann model, Applied Mathematical Modelling 36(5): 2237-2252.

The results are stored in the file "results.txt". The code "main.m" will take several hours to complete all simulations.

Moreover, instead of using the second order non-isotropic color gradient, it is possible to use isotropic color gradient of the following reference:
R2) Leclaire, S., Reggio, M. and Trépanier, J.-Y. (2011). Isotropic color gradient for simulating very high-density ratios with a two-phase flow lattice Boltzmann model, Computers and Fluids 48(1): 98-112.

In average, results are better with the isotropic discretisation.

If these m-files were useful in your work, please cite:
1) The above two references R1 and R2 accordingly.


If you have downloaded these m-files, the following references could also be of interest:

Generalization of this model to multiple immiscible components:
R3) Leclaire, S., Reggio, M. and Trépanier, J.-Y. (2013) Progress and investigation on lattice Boltzmann modeling of multiple immiscible fluids or components with variable density and viscosity ratios. Journal of computational physics 246: 318-342.

Additional terms for modeling variable density ratios:
R4) Leclaire, S., Pellerin, N., Reggio, M. and Trépanier, J.-Y. (2013) Enhanced equilibrium distribution functions for simulating immiscible multiphase flows with variable density ratios in a class of lattice Boltzmann models. International Journal of Multiphase Flow 57: 159-168.

Introduction of MRT and model validation for simulating unsteady multiphase flows:
R5) Leclaire, S., Pellerin, N., Reggio, M. and Trépanier, J.-Y. (2014) Unsteady immiscible multiphase flow validation of a multiple-relaxation-time lattice Boltzmann method. Journal of Physics A: Mathematical and Theoretical 47(10): 105501.

Introduction of the cascaded collision (with unit density ratio):
R6) Leclaire, S., Pellerin, N., Reggio, M. and Trépanier, J.-Y. (2014) Multiphase flow modeling of spinodal decomposition based on the cascaded lattice Boltzmann method. Physica A: Statistical Mechanics and its Applications 406: 307-319.

Others related publications may be found there:

Comments and Ratings (9)


Seb (view profile)

@Qingqing Gu
Thank you for your interest and I'm happy you appreciate the code. Please take a look at the comment I gave below the 6 May 2015.

Qingqing Gu

Dear Leclaire,
Many thanks for your code! I am wondering why the gradient weight for standard anisotropic color gradient in the line 55, collisionMultiPhase.m is [0 1 1 1 1 1 1 1 1] instead of [0 1/6 1/6 1/6 1/6 1/6 1/6 1/6 1/6] as mentioned in your R2(Eq19 and Eq20)?
Many thanks in advance.


Seb (view profile)

@Kaige Shi

I do not know references for "pressure" or "velocity" boundary conditions when two or more components are in contact with the boundary.

I know that single-phase flow "pressure" or "velocity" boundary conditions are working with this two-phase flow model as long as only one of the component is in contact with the boundary. This is working because the phases in this model are strongly immiscible. However, this is only a particular situation and is not general.

Maybe a first step would be to check the work of Lou et al. [1] who have study outflow boundary condition for two-phase flows in the lattice Boltzmann method. This reference might be a starting point.

If eventally you are looking for wetting boundary conditions, one starting point would be to check our last work [2].

[1] Lou, Q., Guo, Z., & Shi, B. (2013). Evaluation of outflow boundary conditions for two-phase lattice Boltzmann equation. Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics, 87(6), 063301. doi:10.1103/PhysRevE.87.063301

[2] Leclaire, S., Abahri, K., Belarbi, R., & Bennacer, R. (2016). Modeling of static contact angles with curved boundaries using a multiphase lattice Boltzmann method with variable density and viscosity ratios. International Journal for Numerical Methods in Fluids, n/a–n/a. doi:10.1002/fld.4226

Kaige Shi

Dear Leclaire,

Thanks for your code. I am wondering how to apply pressure or velocity boundary condition to LBM when there are two components. Do you have reference about it?
Thanks in advance.


Rex Li

Rex Li (view profile)

Hi Sébastien,

Many thanks for your reply. Wah didn't realise that relation comes from your newer publication - now I really got a collection of your publications lol..
Thanks for your help



Seb (view profile)

Hi Rex,

Thank you for your interest in the multiphase model I have submitted. You made a good remark that at I also did not understand very well in the beginning. The point is that the original (non-isotropic) color gradient is, in fact, not a gradient approximation. It estimates the gradient orientation OK, but not the gradient norm. So that is why there is multiplication factor between the two color “gradients”. Indeed, if you take a look at the Eq. (14) in the following publication:

Leclaire, S., El-Hachem, M. and Trépanier, J.-Y., Reggio, M. (2014) High Order Spatial Generalization of 2D and 3D Isotropic Discrete Gradient Operators with Fast Evaluation on GPUs. Journal of Scientific Computing 59

You will found out that A=6 with the original gradient discretization and not A=1 as it should be to correctly approximate the gradient norm. From the LBM perspective, this is not really a problem since the surface tension is proportional to the color gradient. So, in the end, a surface tension multiplication coefficient factor can be used to correctly level the strength of the surface tension.

I hope that help,


Rex Li

Rex Li (view profile)

Rex Li

Rex Li (view profile)

Rex Li

Rex Li (view profile)

Hi Leclaire,
Million thanks for your code ! It is very very useful for me to learn the color-gradient model (I came as a fan of Shan-Chen model but gradually I realised it got some limitations killing all my hope).

I am just wondering that in the line 99-101,main.m, you set the options.aR of isotropic color gradient to be 6 times bigger than that of anisotropic color gradient. I was trying to search in your R2 (Computers and Fluids 48(1): 98-112) however and could not find such a relation. Yeah, I sort of realised that by using isotropic color gradient, the absolute magnitude of the color gradient is smaller compared to using anisotropic color gradient, so dose this relation have some effects to put isotropic & anisotropic color gradients on the same footing ?

Many thanks in advance




Update of the description.


Update of the description.


Description correction.


Description update.


I just change the text in the description

MATLAB Release
MATLAB 7.11 (R2010b)

Inspired by: plotclr

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Win prizes and improve your MATLAB skills

Play today