two phase flow by shan chen

Dear guys

  1. I now boundary condition (pressure boundary condition & velocity boundary condition) in Single phase for the lattice Boltzmann implemented that attached below. I implementing this method without any change for each component in MFMC code but I can’t running that.
  2. for calculating the interaction force, the effective mass (e.g. psi(x)) for the neighbor nodes (e.g. psi(x+e) & psi(x-e)) should be known. In the case of inlet velocity, the function psi(x-e) a is out of the computational domain. How can I compute the interaction force for inlet boundary nodes?

for the left boundary condition:
rho_in=2.05
ux_in=1.0d0-(ff(2,1,y)+ff(4,1,y)+ff(0,1,y)+2.0d0*(ff(6,1,y)+ff(3,1,y)+ff(7,1,y)))/rho_in
rho_in=(ff (2,1,y)+ff (4,1,y)+ff1(0,1,y)+2.0d0*(ff (6,1,y)+ff (3,1,y)+ff (7,1,y)))/(1-ux_in)
ff(1,1,y)=ff1(3,1,y)+2.0d0/3.0d0rho_in1ux_in1
ff(5,1,y)=ff(7,1,y)-(ff(2,1,y)-ff(4,1,y))/2.0d0+1.0d0/6.0d0rho_in1ux_in
ff(8,1,y)=ff(6,1,y)+(ff(2,1,y)-ff(4,1,y))/2.0d0+1.0d0/6.0d0rho_in1ux_in
For the right boundary condition:

rho_out=2.0
ux_out=-1.0d0+(ff(2,lx,y)+ff(4,lx,y)+ff(0,lx,y)+2.0d0*(ff(1,lx,y)+ff(5,lx,y)+ff(8,lx,y)))/rho_out
ff(3,lx,y)=ff(1,lx,y)-2.0d0/3.0d0rho_outux_out
ff(6,lx,y)=ff1(8,lx,y)-(ff(2,lx,y)-ff(4,lx,y))/2.0d0-1.0d0/6.0d0rho_outux_out
ff(7,lx,y)=ff(5,lx,y)+(ff(2,lx,y)-ff(4,lx,y))/2.0d0-1.0d0/6.0d0rho_outux_out