Determining physical properties in multiphase D2Q9 lbm

I am simulating a problem in which water leaves the system and i measure the amount left. I am able to measure the density of the water that has left the system in each iteration in lattice density units. I tried to use the critical density to co-relate the rho_lb to rho_p all the while calculating delt by delt=(nu_lb/nu_p)*delx^2. But now i realize that the density i have is in two dimension, i.e. i have mu/lu^2 while the real density is kg/m^3. So how exactly do i calculate the density in physical units? Also, the equation for delt given above would lead to multiple values of delt when multi-component simulations are run. So is there any other relation that has to be used?

And finally, i have another doubt not so much related to the other two. I have done single phase thermal simulations in which i took thot=1 and tcold=0. But if i want to simulate boiling of water, how do i relate the temperature to the boiling point of water and how do i convert the temperature in lattice units back to physical units?

Please feel free to point out any mistakes in the logic i have presented above and any reply will be appreciated.
Thank you in advance :slight_smile:

Dear Githin,

to convert the density to physical units, the trick is to assume that the domain has dimensions Nx x Ny x 1. This way, you can simply calculate the 3D conversion factor and use it in your 2D system.

In multicomponent systems, you can do the whole unit conversion for one component to get the scaling factors and then apply these to obtain the numerical parameters for the second/third etc. components. Since the conversion factor for time (in your case Delt) will then depend on the relaxation parameter and the physical viscosity of this first component, you will then end up with a different relaxation time for the second component. The timestep (Delt) will of course remain the same.

Concerning non-dimensionalisation, I would recommend the book

Krüger et al.: “The Lattice Boltzmann Method: Principles and Practice”, Springer (2016)

which contains a very detailed description on how to convert units with several examples.



Dear kk,

Thank you for the quick reply. It was really helpful and i think i understand. But i have another doubt. I need to simulate air and water and i have been doing so assuming the air-air and air-water interactions are minimum compared to that of water. So my code calculates inter-particle forces only for water and not for air. But here again i have the problem of density. I am simulating the atmosphere in contact with liquid water. So i am confused about what the density of air and the density of water vapor should be. the density of liquid water is about 528. i did some calculations and got a value for air as about 9.5 lattice density units. This is fine but using critical density and actual humidity, i get the value of atmospheric water vapor density as around 0.6 lattice density units which is creating a huge suction effect as water is at 528. Please suggest a way to fix this.
Also, if i were to include inter-particle forces for air, can i still take any arbitrary value for G_air or does it also depend on G_water (currrently G_water=-120)? also do i calculate the value of G_ads_air in terms of contact angle like given in the book by sukop and thorne?

Thanks in advance

Dear Githin,

So, just so that I understand correctly. You are using a shan-chen multicomponent approach to simulate a water-air system. To adjust the densities of each component, you modify the equation of state of each component. So you have forces acting between different species (Fij) and “inside” of each species (Fii) dependent on the parameters Gij and Gii.

  1. If I understand correctly, the forces between different species Fij only act on your water phase. This would mean the air-phase is decoupled from the water-phase, which is imho not physical. If water moves, it should have an effect on the air-phase (it should displace air).

  2. Suction: I am assuming that the forces between species Fij are small. This means that if you choose a Gii for each phase, you will have (inside this component) a vapor and liquid density, which are rather close to predicitions from the maxwell-construction.
    I do not know which equationof state you are using, but if the densities you choose for this component are far away from the equilibrium densities, then the density values will change. E.g. if your liquid density of water is too small and the vapor density of water is too large, then mass from the vapor phase will be transferred from vapor to liquid, which changes the density and also the number of cells occupied by each phase. This can happen in both ways and would explain this suction effect.

  3. Modelling contact angles: You can of course implement the contact angle forces as described in the sukop/thorne book and try this out. However, I do not think that their analytical expression for predicting the contact angles applies to your rather complex case.

Here are two publications, where the model has been used in a similar way (especially the first one might be interesting for you)

Stiles, Christopher D.; Xue, Yongqiang:
“High density ratio lattice Boltzmann method simulations of multicomponent multiphase transport of H2O in air”
Comput. Fluids 131, 81-90 (2016)

Liu, Malin; Yu, Zhao; Wang, Tiefeng; Wang, Jinfu; Fan, Liang-Shih:
“A modified pseudopotential for a lattice Boltzmann simulation of bubbly flow”
Chem. Eng. Sci. 65, 5615–5623 (2010)



Dear kk,

Thank you for the quick reply.

If i am understanding it right, u r saying that in order to get proper physical resemblance, i need to take into account all the G values, i.e. air-air, water-water, air-water adn water-air. Currently i am only taking water-water and water-solid into consideration and the effect of air flow is purely due to its effect on overall velocity. But if i am to take all those into consideration, how do i decide on the values of Gij for air-water and water-air?

i am using the G_water-water value of -120, rho0 of 200 and psi0 of 4. I am using the value of liquid density the same as i got in a simple equilibrium simulation i got, which is 529 and 89. But now i am confused. The critical density of water at STP is 322 kg/m^3 and in lattice units it is 200. So, liquid equilibrium density should be 1000*(200/322)=621 and vapor should be 0.02*(322/200)=0.0124. Obviously my method is wrong. So how do i actually do this calculation?

Thanks for the recommended publications. I will check them out as soon as possible
Thanks in advance

Dear kk,

Please disregard my last post. The papers you suggested have been of great help and most of my doubts have been cleared. One final doubt i have is how they decide on the parameters in the equation of states. When they incorporate an equation of state into lbm, it usually has atleast two parameters (a and b), if not more, and these parameters are linked to the critical temperatures. I get that since we need to find the critical temperatures we substitute for a and b in order to get critical properties. But where do they get the values of a and b to substitute? It has to be somewhat related to the species present. But there is such contrast between the a and b values in real life vs and b values in the lbm units which they substitute. So basically my doubt is how do they decide on what values of a and b to substitute for different equations of state.

Thanks in advance

Dear Githin,

most of the time, these equations of state are used, people apply these values for the coefficients. These are the values that were used in the first publication, where the usage of different equations of state was introduced (Yuan & Schaefer: “Equations of state in a lattice Boltzmann model”. Phys. Fluids 18, 042101 (2006) ). So, probably everyone uses those because they are the most “analysed” values. Also, in case of the Carnahan-Starling equation of state, the choice of parameters a and b can strongly simplify the equation (in case of a=1, b=4). You could of course choose different parameters, but you would have to look at how the resulting potential function behaves numerically.



Dear kk,

Thank you for your quick reply.
i still have a few doubts regarding these. First, wouldnt the choice of a and b be what determines which fluid it is? I mean, when using the real EOS, we use values of Tr and Pr according to the compound and so get different values of a and b. In this case, he is simulating water so i guess the values he gave are for them so it is what i need but generally how do we simulate different fluids?

Secondly and most importantly, I wrote a code trying to use other equations of state but i just cant make it work. Either after a few iterations the value of psi goes complex and so MATLAB cant plot it or phase separation never occurs and everything just diffuses away. I found an existing post in palabos with the same problem. his code is very similar to mine. All the relevent parts are the same but there is no proper solution mentioned. The link to the code is given below. Please go through the code and help me find out what im doing wrong.,7913,7913#msg-7913

Thanks in advance

Dear Githin,

I was thinking you want to simulate a system that contains multiple components where you want to change the equation of state to get different densities in each component. Looking at your previous post, it looks like you are simulating a system containing a SINGLE component which is present in two phases (liquid and vapor). At least, the code in the link is used for SINGLE component multiphase flows.

To your problem: When you change the equation of state, the value of G becomes redundant. Only its sign is important to keep the radicand in the “modified” version of psi positive. So if psi goes complex, you might have to switch the sign of G, e.g. from 1 to -1.
Assuming your implementation is correct, if no phase separation occurs, then the temperature you are using might be too large. To enforce phase separation you would have to use a temperature that is below the critical one.



Dear kk,

I am doing multi component simulation. I was doing it using the simple SC EOS and so i am trying to convert it into a better EOS. I thought it would be best to try it on a simple SCMP problem first and then move on to MCMP. the code i shared from palabos had the same problem i had, just was in SCMP instead of MCMP.

I tried what you asked. Changing G to -1 helped make it run a bit longer but i still cannot observe any phase separation for any value of T. I will give my SCMP code below. i am not able to pinpoint the error.

I adapted the code from another one which i had written some time ago. So ignore the code written for the walls and the wall interaction as i did not place any walls.

clear all
close all

nx=400; %matrix x dimension
ny=400; %matrix y dimension
density=20; %density


to = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
ex = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
ey = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
oppf = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];


f=ones(9,ny,nx); %density distribution function for air

rho=ones(ny,nx); %density of air
fx=zeros(ny,nx); %intra fluid x force
fy=zeros(ny,nx); %intra fluid y force
fwx=zeros(ny,nx); %intra fluid x force
fwy=zeros(ny,nx); %intra fluid y force
ran=rand([ny nx]);
for i=1:9

for i=1:ny
for j=1:nx
if solid(i,j)==1

for z=1:tend
%calculating macrovariables

fx=to(2)*circshift(psi,[0 -1])+to(6)*circshift(psi,[-1 -1])+to(9)*circshift(psi,[1 -1])-to(4)*circshift(psi,[0 1])-to(8)*circshift(psi,[1 1])-to(7)*circshift(psi,[-1 1]);
fy=to(3)*circshift(psi,[-1 0])+to(6)*circshift(psi,[-1 -1])+to(7)*circshift(psi,[-1 1])-to(5)*circshift(psi,[1 0])-to(8)*circshift(psi,[1 1])-to(9)*circshift(psi,[1 -1]);

% fx=to(2).*2.*circshift(psi,[0 -1])+to(6)./2.*circshift(psi,[-1 -1])+to(9)./2.*circshift(psi,[1 -1])-to(4).*2.*circshift(psi,[0 1])-to(8)./2.*circshift(psi,[1 1])-to(7)./2.*circshift(psi,[-1 1]);
% fy=to(3).*2.*circshift(psi,[-1 0])+to(6)./2.*circshift(psi,[-1 -1])+to(7)./2.*circshift(psi,[-1 1])-to(5).*2.*circshift(psi,[1 0])-to(8)./2.*circshift(psi,[1 1])-to(9)./2.*circshift(psi,[1 -1]);

% not sure if it is this one or the other (commented above)

fwx=to(2)*circshift(solid,[0 -1])+to(6)*circshift(solid,[-1 -1])+to(9)*circshift(solid,[1 -1])-to(4)*circshift(solid,[0 1])-to(8)*circshift(solid,[1 1])-to(7)*circshift(solid,[-1 1]);
fwy=to(3)*circshift(solid,[-1 0])+to(6)*circshift(solid,[-1 -1])+to(7)*circshift(solid,[-1 1])-to(5)*circshift(solid,[1 0])-to(8)*circshift(solid,[1 1])-to(9)*circshift(solid,[1 -1]);

%above wall interaction will be zero as solid=0


%     calculating equilibrium conditions and collision

for i=1:9
    cuNS = 3*(ex(i)*ux+ey(i)*uy);
    feq(i,:,:) = rho .* to(i) .* ...
        ( 1 + cuNS + 1/2*(cuNS.*cuNS) - 3/2*(ux.^2+uy.^2) );
    for j=1:ny
        for k=1:nx
            if solid(j,k)==1
    f(i,:,:) = f(i,:,:) -  (f(i,:,:)-feq(i,:,:))./tau;

%solid collision
for i=1:ny
    for j=1:nx
        if solid(i,j)==1
for i=1:9
    f(i,:,:)=circshift(f(i,:,:),[0,ey(i) ex(i)]);

axis equal off;


I greatly appreciate your help.

Dear Githin,

At a first glance: When calculating psi, you are using the absolute value of the radicand. Remove the abs(…) function and ensure positivity of the radicand only by adjusting the sign of the parameter G. To be consistent, you must use G in the calculation of the force terms fx and fy as well (the value does not matter, but the sign must be taken into account!). You cannot just neglect G.

I also do not know if the critical values are correct. I have no experience in using the PR-EOS, but why are you using an initial density of 20? Also alpha should be a temperature dependent variable.



Dear kk,

Thanks for the quick reply.
the thing about the G was an accident. I was just trying it out and forgot to change it back. I do always use G in both the psi term as well as the fx fy term.
Oh and btw im using RKS EOS.
Regarding alpha, I calculated it based on the value of T i am taking and input it as a constant.
I have cross checked the values of critical variables multiple times and have not noticed any mistakes.

Also, i am using different values of density inorder to try to achieve phase separation. i am not able to achieve it for any temperature and 20 was just one of the values i tried. I tried values ranging between 1 to 500. nothing worked

Any help regarding this problem is appreciated.

Hi Githin,

can u provide the code of the MCMP to me?

kind regards,