Hello guys,
Recently I modified Shan-Chen’s model (1993) following Yuan and Schaefer’s work, where high density ratio is achieved by applied proper eos (PoF, 042101, 2006).
First, I repeat Sukop’s work and the results are correct. The weight factors for psi gradient are w_i=1/36 (i=5,6,7,8), w_i=1/9 (i=1,2,3,4).
Then, the effecitive density, psi, is defined as, alternatively,
psi = sqrt( (p-rho/3.)6/G)
p=aRT[1+brho/4+(brho/4)^2-(brho/4)^3] / (1-brho/4)^3 - arho^2 (C-S eos, a=1, b=4, R=1, T=0.055)
= rho/3.+G/6.0(psi^2)
However the results are different with Yuan’s, where the spurious current is 10^-14 for flat interface, while mine are 0.2. Are there anybody have idea on it, OR do I make any mistakes on it? thank you very much!
Chen
My C codes read,
// lattice density…
for(y=0;y<max_y;y++)
for(x=0;x<max_x;x++)
rho_global[y][x]=f0[y][x]+f1[y][x]+f2[y][x]+f3[y][x]+f4[y][x]+f5[y][x]+f6[y][x]+f7[y][x]+f8[y][x];
// calculate psi…
for(y=0;y<max_y;y++){
for(x=0;x<max_x;x++){
// Psi[y][x]=Psi_0exp(-rho_0/rho_global[y][x]);
dummy=rho_global[y][x];
temp=1./pow((1.-dummy), 3);
temp=RgT(1.+dummy+pow(dummy,2)-pow(dummy,3));
temp-=a*rho_global[y][x];
Psi[y][x] = sqrt((temp-1./3.)*rho_global[y][x]*6./G);
}
}
//peroidic boundary condition
for(y=0;y<max_y;y++){
yp = (y<max_y-1)?(y+1):(0);
yn = (y>0)?(y-1):(max_y-1);
for(x=0;x<max_x;x++){
xp = (x<max_x-1)?(x+1):(0);
xn = (x>0)?(x-1):(max_x-1);
// psi-gradient
Fx=0.0; Fy=0.0;
Fx+= (Psi[y][xp]-Psi[y][xn])4.0; / 1st neighbor*/
Fy+= (Psi[yp][x]-Psi[yn][x])4.0; / 2nd neighbor*/
Fx+= (Psi[yp][xp]-Psi[yp][xn]-Psi[yn][xn]+Psi[yn][xp]);
Fy+= (Psi[yp][xp]+Psi[yp][xn]-Psi[yn][xn]-Psi[yn][xp]);
// interaction forces
Fx = -GPsi[y][x]/36.; Fy = -GPsi[y][x]/36.;
vx=(f1[y][x]-f3[y][x]+f5[y][x]+f8[y][x]-f6[y][x]-f7[y][x])/rho_global[y][x];
vy=(f2[y][x]-f4[y][x]+f5[y][x]+f6[y][x]-f7[y][x]-f8[y][x])/rho_global[y][x];
vx+=Fx*tau/rho_global[y][x]; vy+=Fy*tau/rho_global[y][x];
// … … … … to calculate equilibrium distribution
I have the same questions! first i think the initial densities is important for simulation. but so far i have no idea.how about you now?