Please help me to find mistake in LBM coding (D2Q9) for flow through channel. Inlet boundary condition is uniform velocity inlet and outlet boundary condition is open boundary conditio. I tried to write code but not getting proper result.
!Program for flow through channel
Program Channel_flow
parameter (n=500, m=40)
real:: f(0:8,0:n,0:m), feq(0:8,0:n,0:m)
real:: rho(0:n,0:m), u(0:n,0:m), v(0:n,0:m)
real:: w(0:8)
real:: cx(0:8), cy(0:8),rowout
integer:: i,j, mstep,kk
real:: sumvelo, dx, dy, dt,alpha, Re, omega, u0, rhow,ssum,rhoo
open (2,file='uvfield')
open (3,file='uvely')
open (4,file='vvelx')
open (8,file='timeu')
open (7,file='fkij')
uo=0.2
sumvelo=0.0
rhoo=1.0
dx=1.0
dy=dx
dt=1.0
alpha=0.02
Re= uo*m/alpha
print*, "Re=", Re
omega=1.0/(3.0*alpha+0.5)
mstep=300
!rowout=1.0
! Weighting function
w(0)=4.0/9.0
do i=1,4
w(i)=1./9.
end do
do i=5,8
w(i)=1./36.
end do
cx(0)=0.0 !Allocate cx
cx(1)=1.0
cx(2)=0.0
cx(3)=-1.0
cx(4)=0.0
cx(5)=1.0
cx(6)=-1.0
cx(7)=-1.0
cx(8)=1.0
cy(0)=0.0 !Allocate cy
cy(1)=0.0
cy(2)=1.0
cy(3)=0.0
cy(4)=-1.0
cy(5)=1.0
cy(6)=1.0
cy(7)=-1.0
cy(8)=-1.0
do j=0,20 !Allocate Initial condition
do i=0,40
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do j=21,m-1 !Allocate Initial condition
do i=0,40
u(i,j)=uo
v(i,j)=0.0
end do
end do
do j=1,m-1 !Allocate Initial condition
do i=41,n
u(i,j)=uo
v(i,j)=0.0
end do
end do
do i=0,n
u(i,0)=0.0
u(i,m)=0.0
end do
do j=0,m ! y- component of velocity is zero
do i=0,n
rho(i,j)=rhoo
v(i,j)=0.0
end do
end do
do kk=1,mstep !Start of main loop
call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
call streaming (f,n,m)
call sfbound(f,n,m,uo)
call rhouv (f,rho,u,v,cx,cy,n,m)
print*,u(0,m/2),v(0,m/2),rho(0,m/2),u(n,m/2),v(n,m/2),rho(n,m/2)
write(8,*) kk, u(n/2,m/2), v(n/2,m/2)
end do ! End of main loop
call result (u,v,rho,uo,n,m)
stop
end
!End of the main program
!_________________________________________________________
subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
real:: f(0:8,0:n,0:m), feq(0:8,0:n,0:m)
real:: rho(0:n,0:m), u(0:n,0:m), v(0:n,0:m)
real:: w(0:8), omega
real:: cx(0:8), cy(0:8), t1, t2
integer:: i,j,k
do i=0,n
do j=0,m
t1= u(i,j)*u(i,j)+v(i,j)*v(i,j)
do k=0,8
t2=u(i,j)*cx(k)+v(i,j)*cy(k)
feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.5*t2*t2-1.50*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
end do
end do
end do
do j=0,m
do i=0,n
sum=0.0
do k=0,8
sum=sum+f(k,i,j)
end do
end do
end do
return
end !end of subroutine collision
!_______________________________________________________
subroutine streaming (f,n,m)
real:: f(0:8,0:n,0:m)
integer:: i,j
!streaming
do j=0,m !Bottom to top
do i=n,1,-1 !Right hand to left
f(1,i,j)=f(1,i-1,j)
end do
do i=0,n-1 !Left hand to Right
f(3,i,j)=f(3,i+1,j)
end do
end do
do j=m,1,-1 ! Top to bottom
do i=0,n
f(2,i,j)=f(2,i,j-1)
end do
do i=n,1,-1
f(5,i,j)=f(5,i-1,j-1)
end do
do i=0,n-1
f(6,i,j)=f(6,i+1,j-1)
end do
end do
do j=0,m-1 ! Bottom to top
do i=0,n
f(4,i,j)=f(4,i,j+1)
end do
do i=0,n-1
f(7,i,j)=f(7,i+1,j+1)
end do
do i=n,1,-1
f(8,i,j)=f(8,i-1,j+1)
end do
end do
return
end
!_________________________________________________________
subroutine sfbound (f,n,m,uo)
real f(0:8,0:n,0:m)
integer:: i,j
real:: uo, rhow,rowout
do j=21,m
!velocity inlet on west boundary
rhow=(f(0,0,j)+f(2,0,j)+f(4,0,j)+2.*(f(3,0,j)+f(6,0,j)+f(7,0,j)))/
&(1.-uo)
!print*, rhow
f(1,0,j)=f(3,0,j)+2.*rhow*uo/3.0
f(5,0,j)=f(7,0,j)+rhow*uo/6.0
f(8,0,j)=f(6,0,j)+rhow*uo/6.0
end do
do i=0,n
!Bounce back on south coundary
f(2,i,0)=f(4,i,0)
f(6,i,0)=f(8,i,0)
f(5,i,0)=f(7,i,0)
end do
do i=0,n
!Bounce back on north boundary
f(4,i,m)=f(2,i,m)
f(8,i,m)=f(6,i,m)
f(7,i,m)=f(5,i,m)
end do
!Open boundary condition at the outlet
do j=1,m
f(1,n,j)= 2.*f(1,n-1,j)-f(1,n-2,j)
f(5,n,j)= 2.*f(5,n-1,j)-f(5,n-2,j)
f(8,n,j)= 2.*f(8,n-1,j)-f(8,n-2,j)
end do
!Obstacle at the inet x=40, y=20
do i=0,40
f(2,i,20)=f(4,i,20)
f(6,i,20)=f(8,i,20)
f(5,i,20)=f(7,i,20)
end do
do j=0,20
f(1,40,j)=f(3,40,j)
f(5,40,j)=f(7,40,j)
f(8,40,j)=f(6,40,j)
end do
return
end
!_______________________________________________________
subroutine rhouv(f,rho,u,v,cx,cy,n,m)
real:: f(0:8,0:n,0:m)
real:: rho(0:n,0:m), u(0:n,0:m), v(0:n,0:m)
real:: cx(0:8), cy(0:8),ssum, usum, vsum
integer:: i,j,k
do j=0,m
do i=0,n
ssum= 0.0
do k=0,8
ssum=ssum+f(k,i,j)
end do
rho(i,j)=ssum
end do
end do
do i=1,n
rho(i,m)=f(0,i,m)+f(1,i,m)+f(3,i,m)+2.*(f(2,i,m)+f(6,i,m)+f(5,i,m
&))
end do
do i=1,n
do j=1,m
usum=0.0
vsum=0.0
do k=0,8
usum=usum+f(k,i,j)*cx(k)
vsum=vsum+f(k,i,j)*cy(k)
end do
u(i,j)=usum/rho(i,j)
v(i,j)=vsum/rho(i,j)
end do
end do
do j=0,m
v(n,j)=0.0
end do
do j=0,20
do i=0,40
u(i,j)=0.0
v(i,j)=0.0
end do
end do
return
end
!__________________________________________________________
subroutine result(u,v,rho,uo,n,m)
real:: u(0:n,0:m),v(0:n,0:m)
real:: rho(0:n,0:m),strf(0:n,0:m),rhoav
integer:: i,j
open(9,file='streamf')
!stream function calculation
strf(0,0)=0.0
do i=0,n
rhoav=0.5*(rho(i-1,0)+rho(i,0))
if(i.ne.0) strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
do j=1,m
rhom=0.5*(rho(i,j)+rho(i,j-1))
strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
end do
end do
!__________________________________________________________
write(2,*)"varaibles= X, Y, U, V, S"
write(2,*)"Zone", "I=", n+1, "J=", m+1,"F= Block"
do j=0,m
write(2,*) (v(i,j), i=0,n)
end do
do j=0,m
write(9,*) (strf(i,j),i=0,n)
end do
do j=0,m
write(3,*) j/float(m),u(n/2,j)/uo,u(n/4,j)/uo,u(3*n/4,j)/uo
end do
do i=0,n
write(4,*) i/float(n),v(i,m/2)/uo
end do
return
end