# Problems about SC multiphase modael

Hi,
Everyone! I’m a beginners for LBM and now I’m coding the Shan and Chen Multiphase model by myself.I used The Lattice Boltzmann Modeling(an introduction to geoscienctists and engineers,by M.C. Sukop,D.T. Thorne, Jr.) as the reference book.The problem I simulated is phase separation with a initial density 200 with (0,1) random density fluctuation.It seems that there isn’t any problem in the program(including the parameters) because I nearly copied the program from the book.But the velocity of the system seems to be increasing very fast at the initial steps and even past the sound velocity.I don’t know why these phenomenon happen.Could anybody help me?
Listed below is the major code I programmed:

int _tmain(int argc, _TCHAR* argv[])
{
int Count;
time=0;
Count=0;
Init();
Time_Set=100000.0dt;
while (time<Time_Set)
{
compute_phase_force();
collision();
flow();
compute_macro_var();
time+=dt;
Count++;
if (Count %50==0) //for monitoring the parameters in the computation
{
printf(“Max_velocity=%.15lf,max_rho=%.6lf,sum_mass=%.5f\n”,max_velocity,max_rho,sum_mass);
writedata_to_file();
}
if (abs((max_velocity-pre_max_velocity)/max_velocity)<1e-5) //check the convergence
{
writedata_to_file();
}
}
return 0;
}
void Init()
{
e[0][0]=0;e[0][1]=0;
e[1][0]=1;e[1][1]=0;
e[2][0]=0;e[2][1]=1;
e[3][0]=-1;e[3][1]=0;
e[4][0]=0;e[4][1]=-1;
e[5][0]=1;e[5][1]=1;
e[6][0]=-1;e[6][1]=1;
e[7][0]=-1;e[7][1]=-1;
e[8][0]=1;e[8][1]=-1;
//initial the 9 opposite direction
op_e[0]=0;
op_e[1]=3;op_e[2]=4;op_e[3]=1;op_e[4]=2;
op_e[5]=7;op_e[6]=8;op_e[7]=5;op_e[8]=6;
//initial the parameters for computation
R=1.0;
dx=1;dt=1;
c=dx/dt;
c_s=c/sqrt(3.0);
T0=c_s
c_s/R;
G=-120.0; //G,psi0 and rho0 are the same as the reference book!!!
psi0=4.0;
rho0=200.0;
tau=1;
vis=(tau-0.5)dxdx/dt/6;

init_macro();
init_micro();
}

void init_macro()
{
int i,j;

srand(1);
for(i=0;i<Nx;i++)
{ for(j=0;j<Ny;j++)
{
u[i][j]=0;
v[i][j]=0;
rho[i][j]=+AverageRandom(0,1); //AverageRandom is a function which produce randon real number between 0 and 1
force[i][j][0]=force[i][j][1]=0;
}
}
}
void init_micro() //???
{
int i,j,k;
for(i=0;i<Nx;i++)
{
for(j=0;j<Ny;j++)
{
compute_feq(u[i][j],v[i][j],rho[i][j],feq);
for(k=0;k<Q;k++)
f[i][j][k]=feq[k];
}
}
}

void compute_feq(double u,double v,double rho,double *feq)
{
double temp_u2, temp_rho, uav, vmu;
int k;

u /= c;
v /= c;
temp_u2 = 1-1.5*(u*u + v*v);
temp_rho = 4.0*rho/9.0;
uav = u+v;
vmu = v-u;

feq[0] = temp_rho * temp_u2;

temp_rho *= 0.25;
feq[1] = temp_rho * (temp_u2 + 3*u + 4.5*u*u);
feq[2] = temp_rho * (temp_u2 + 3*v + 4.5*v*v);
feq[3] = feq[1] - temp_rho * 6*u;
feq[4] = feq[2] - temp_rho * 6*v;

temp_rho *= 0.25;
feq[5] = temp_rho * (temp_u2 + 3*uav + 4.5*uav*uav);
feq[6] = temp_rho * (temp_u2 + 3*vmu + 4.5*vmu*vmu);
feq[7] = feq[5] - temp_rho * 6*uav;
feq[8] = feq[6] - temp_rho * 6*vmu;


}
void compute_macro_var()
{
int i,j,k;
rho_max=0;
max_velocity=0;
sum_mass=0;
for(i=0;i<Nx;i++)
{
for(j=0;j<Ny;j++)
{

   rho[i][j] = 0;
u[i][j] = v[i][j] = 0;
for(k=0; k<Q; k++)
{
u[i][j] += e[k][0] * f[i][j][k];
v[i][j] += e[k][1] * f[i][j][k];
rho[i][j] += f[i][j][k];
}


}
void collision()
{
int i,j,k;
double u_eq,v_eq;
for(i=0; i<Nx; i++)
{
for(j=0; j<Ny; j++)
{
u[i][j]+=tauforce[i][j][0]dt/rho[i][j]; //Add the long range interaction force’s influence to the velocity
v[i][j]+=tau
force[i][j][1
dt/rho[i][j];

        compute_feq(u[i][j],v[i][j], rho[i][j], feq);
for(k=0; k<Q; k++)
ftemp[i][j][k] = f[i][j][k]+(feq[k]-f[i][j][k])/tau;
}


}
}

void flow()
{
int i,j,ip,in,jp,jn,ik,jk,k;
rho_max=0;

//periodic conditions both for x and y directions
for(i=0;i<Nx;i++)
{
in=(i>0)?(i-1):(Nx-1);
ip=(i<Nx-1)?(i+1):(0);
for(j=0;j<Ny;j++)
{
jn=(j>0)?(j-1):(Ny-1);
jp=(j<Ny-1)?(j+1):(0);
f[i][j][0]=ftemp[i][j][0];
f[ip][j][1]=ftemp[i][j][1];
f[i][jp][2]=ftemp[i][j][2];
f[in][j][3]=ftemp[i][j][3];
f[i][jn][4]=ftemp[i][j][4];
f[ip][jp][5]=ftemp[i][j][5];
f[in][jp][6]=ftemp[i][j][6];
f[in][jn][7]=ftemp[i][j][7];
f[ip][jn][8]=ftemp[i][j][8];
}
}


}
void compute_phase_force()
{
int i,j,in,jn,ip,jp;

for(i=0;i<Nx;i++)
for(j=0;j<Ny;j++)
psi[i][j]=4.0*exp(-200/(rho[i][j]));

for(i=0;i<Nx;i++)
{
in=(i>0)?(i-1):(Nx-1);
ip=(i<Nx-1)?(i+1):(0);
for(j=0;j<Ny;j++)
{
jn=(j>0)?(j-1):(Ny-1);
jp=(j<Ny-1)?(j+1):(0);

		force[i][j][0]=0.;
force[i][j][1]=0.;

force[i][j][0]+=4*e[1][0]*psi[ip][j];force[i][j][1]+=4*e[1][1]*psi[ip][j];
force[i][j][0]+=4*e[2][0]*psi[i][jp];force[i][j][1]+=4*e[2][1]*psi[i][jp];
force[i][j][0]+=4*e[3][0]*psi[in][j];force[i][j][1]+=4*e[3][1]*psi[in][j];
force[i][j][0]+=4*e[4][0]*psi[i][jn];force[i][j][1]+=4*e[4][1]*psi[i][jn];
force[i][j][0]+=1*e[5][0]*psi[ip][jp];force[i][j][1]+=1*e[5][1]*psi[ip][jp];
force[i][j][0]+=1*e[6][0]*psi[in][jp];force[i][j][1]+=1*e[6][1]*psi[in][jp];
force[i][j][0]+=1*e[7][0]*psi[in][jn];force[i][j][1]+=1*e[7][1]*psi[in][jn];
force[i][j][0]+=1*e[8][0]*psi[ip][jn];force[i][j][1]+=1*e[8][1]*psi[ip][jn];

force[i][j][0]*=-G*psi[i][j];force[i][j][1]*=-G*psi[i][j];
}


}
}

void writedata_to_file()
{
int i,j;
double pressure,phase_force,s_force,net_force;
FILE *f_result;
f_result=fopen(“data.dat”,“wb”);
fprintf(f_result,“TITLE = “SC model XY Plot”\n”);
fprintf(f_result,“VARIABLES = “x”,“y”,“Velocity_x”,“Velocity_y”,“density”,“Pressure”,“force x”,“force y”,“Phase force”\n”);
fprintf(f_result,“ZONE T=“BIG ZONE”, I=%i, J=%i, F=POINT\n”,Nx,Ny);

for(i=0;i<Nx;i++)
{
for(j=0;j<Ny;j++)
{
pressure=rho[i][j]*R*T0+G*R*T0/2*(4.0*exp(-200.0/rho[i][j]))*(4.0*exp(-200.0/rho[i][j]));
phase_force=sqrt(force[i][j][0]*force[i][j][0]+force[i][j][1]*force[i][j][1]);
fprintf(f_result,"%i\t%i\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",i,j,u[i][j]+tau*force[i][j][0]*dt/rho[i][j]/2,v[i][j]+tau*(force[i][j][1]*dt/rho[i][j]/2,rho[i][j],pressure,force[i][j][0],force[i][j][1],phase_force);
}
}
fclose(f_result);


}

Hi,
firstly I have to say that I don’t any experience with the multiphase (single particle distribution) ShanChen model… actually I play a lot with Multi-Component (multi particle distribution) ShanChen Model. Well… I guess you like a lot Sukop’s book then you should know the difference between the two. (they are both in the book).

## You wrote: --------------- But the velocity of the system seems to be increasing very fast at the initial steps and even past the sound velocity.I don’t know why these phenomenon happen

I would like to ask you what kind of physical problem you are looking at… if the speed bigger than the speed of sound that you get is localized in specific region of your computation domain (for examples closer to the interface between the two phases (is the interface curved by any chance?) ) … Does the qualitative results of the simulation looks good ?
And how do you calculate your velocity field?

if the answer is: YES, the images looks pretty and my huge speed is localized close to the interface between the two phase … maybe your looking at the Spurious velocities (not physical velocity that forms at curved interface due to a not totally clear (yet) lack of isotropy of the interparticle forces … and many other things )

I won’t go trough your code … but what I can tell you is that Sukop had put on the web some of his “LB codes for dummies” (the kind of codes I like given my poor ability in coding) and between them there was also a code for the ShanChen model.
The web page where you can download the codes should be: http://www.fiu.edu/~sukopm/Research/LB.htm (I don’t know why, but today doesn’t work).
In case you could write to Sukop or Thorne. I think that the codes were included in the price book

If you want some more information about the single particle ShanChen model, you could look at a paper of Sukop about the Cavitation problem (this is done with the single component ShanChen model … They give some more values for parameters G and epsilon and so on … )

the paper is:
Phys. Rev. E 71, 046703 (2005)

Maybe one of these days a singleCoponent ShanChen Model will be available in the OpenLB (you have the MultiComponent there ). Then we could cross-check our results…

ciao
Andrea

Dear Robert
hello
This problem is not weird.The Shan and Chen model can have some kind of this problem.
we can not apply high density ratio on this model.high spaturious velocity occurs on boundary…
First of all i wondre which applicatication do you want to model??
I can tell you to what extent you can use this model…
Is this the complete code??? I can rin it at convenient time???
Best regards

Dear Robert,

If you need Shan-Chen with OpenLB. I can send it to you.

Dear Alex
Hello
Is Open LB Project capable to model Shan and Chen theory??
To what extent can you use this model??
Is it capable to model high density ratio??
best regards

Dear moslemi,

I added Shan-Chen to my files. Probably the right word for it is postprocessor. It’s not in the libraries, but in the main.cpp files. You can use this model upto liquid gas density ratio 100. It’s usual Shan-Chen. You can use extended model for 188 - but it’s not OpenLB, but my small code. But there are special equation of states where you can reach upto 1000 and higher. I can recommend you to look into the paper Peng Yuan and Laura Shaefer “Equations of state in a lattice Boltzmann model”. I want to include these models to the libraries of OpenLB and will do it soon.

Alex

Dear Alex
hello
Thank you
Hamid

Dear Hamid,

Alex

Dear Alex
Hello
I hope this message recieves you in good health and spirits.I am trying to develop a code to model Multiphase flow based on He’s model.
Did you see this model?
I cite it below:
A Lattice Boltzmann Scheme for Incompressible
Multiphase Flow and Its Application in
Simulation of Rayleigh–Taylor Instability1
Xiaoyi He, Shiyi Chen, and Raoyang Zhang
Journal of Computational Physics 152, 642–663 (1999)
thank you
Best Regards
Hamid

Hi Hamid,

I tried to implement this model with Shan-Chen equation of state and spend a lot of time trying to make it work. But no success. It’s killing me…

Actually, there are some interesting and non-understandable things for me in this model (my comments):

1. You need your force to be implemented as Fc_i/(cs^2 \rho) f^eq - something like that as I remember. Only because your force term is proportional to the f^eq (equilibrium distribution function) you can obtain the factor ?(U)-?(0) which is proportional to Mach number for pressure gradient. If you try to do forces in other representation view you won’t achieve it or probably I didn’t find the way to implement

2. The most confusing part for me is how they actually introduce the order parameter distribution function and how they state that surface tension and gravity are not changing local density and therefore are not involved in the order parameter equation. It’s confusing for me. If you look through the papers Swift et al. or Landau about Binary Models you will find convection-diffusion but not hydrodynamic equation for the order parameter (!) involving chemical potential and special type of equation of states for order parameters. How do they introduce this order parameter and how it’s connected to real physics - I don’t know, I’d like to hear any comments?

3. I’ve heard that you can’t achieve big density ratios with this model - around 10 or something like that. But it’s what I’ve heard and read from papers.

4. Also confusing part for me - it’s incompressible limit which they are using - you will find somewhere there in the paper how authors take the full derivative d grad\psi /dt that it equals to u \grad \psi - if you do it thoroughly you can see that it’s true only for incompressible equations. But how we can talk about incompressible equations if we have density ratio about 10. Then the binary liquid model should be introduced - but for this you need absolutely different equations.

So I’m confused, I tried to implement it but no success. But I’m still really interested in it. Especially in the equation for order parameter. If you can successfully implement it or any comments from you please let me know.

I will be away till the next Tuesday,