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_sc_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]+=tauforce[i][j][1dt/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);
}
…