[hamed.majeed]

Hi guys,

I am working on lattice boltzmann method LBM and currently following the book LBM by A. A. Mohamad.

I have issues in the FORTRAN code for poiseuille flow for flat plates. The code file is as follows:

!computer code for flow in a channel
parameter (n=1000,m=40)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m), rho(0:n,0:m)
real w(0:8), cx(0:8), cy(0:8)
real u(0:n,0:m),v(0:n,0:m)
integer i
open(2,file='result')
!
uo=0.2
rhoo=1000.00
dx=1.0
dy=dx
dt=1.0
alpha=0.02 !nu lattice
Re=uo*m/alpha !Re=400
omega=1./((alpha*3.*dt/(dx*dx))+0.5)
mstep=20
w(0)=4./9.
do i=1,4
w(i)=1./9.
end do
do i=5,8
w(i)=1./36.
end do
cx(0)=0
cy(0)=0
cx(1)=1
cy(1)=0
cx(2)=0
cy(2)=1
cx(3)=-1
cy(3)=0
cx(4)=0
cy(4)=-1
cx(5)=1
cy(5)=1
cx(6)=-1
cy(6)=1
cx(7)=-1
cy(7)=-1
cx(8)=1
cy(8)=-1
do j=0,m
do i=0,n
rho(i,j)=rhoo !Inital value density is set to rhoo=1000.00
u(i,j)=0.0 !Initial value fo u and v velcities is set to 0
v(i,j)=0.0
do k=0,8
f(k,i,j)=0.0
end do
end do
end do
do j=1,m-1
u(0,j)=uo !Velocity at inlet is set to uo at the left most inlet (U lattice)
v(0,j)=0.0
end do

!Main Loop
do kk=1,mstep
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)
end do
!End of Main Loop
do j=0,m
print *,kk,j,u(300,j)
end do
stop
end
!End of Main Program

subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8),cx(0:8),cy(0:8)
real u(0:n,0:m),v(0:n,0:m)
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.5*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
end do
end do
end do
return
end

subroutine streaming(f,n,m)
real f(0:8,0:n,0:m)
!Streaming
do j=0,m
do i=n,1,-1 !Right to left
f(1,i,j)=f(1,i-1,j)
end do
do i=0,n-1 !Left 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) !Boundary conditions
real f(0:8,0:n,0:m)
real uo
!Boundary condition with known velocity at west wall
do j=0,m
rhow=(f(0,0,j)+f(2,0,j)+f(4,0,j)+2.0*(f(3,0,j)+f(6 ,0,j)+f(7,0,j)))/(1-uo)
f(1,0,j)=f(3,0,j)+2.0*rhow*uo/3.0
f(5,0,j)=f(7,0,j)-0.5*(f(2,0,j)-f(4,0,j))+rhow*uo/6.0
f(8,0,j)=f(6,0,j)+0.5*(f(2,0,j)-f(4,0,j))+rhow*uo/6.0
end do
!Bounce back on south boundary
do i=0,n
f(2,i,0)=f(4,i,0)
f(5,i,0)=f(7,i,0)
f(6,i,0)=f(8,i,0)
end do
!Bounce back on north boundary
do i=0,n
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 outlet
do j=1,m
f(1,n,j)=2.0*f(1,n-1,j)-f(1,n-2,j)
f(5,n,j)=2.0*f(5,n-1,j)-f(5,n-2,j)
f(8,n,j)=2.0*f(8,n-1,j)-f(8,n-2,j)
end do
return
end

subroutine rhouv(f,rho,u,v,cx,cy,n,m) !To find net rho and velocities
real f(0:8,0:n,0:m),rho(0:n,0:m),u(0:n,0:m),v(0:n,0:m), cx(0:8),cy(0:8)
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 j=0,m
do i=0,n
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=1,m
v(n,j)=0.0
end do
return
end
!==========end of program===============

[flotus1]

And what exactly is the problem with your code?

[hamed.majeed]

The "u" velocity plot along the "y" direction shows the velocity to be zero...however, in actual it shouldn't!

[flotus1]

Having a quick look at the code I spotted 2 obvious problems. I assume there are more.

1) You plot the velocity profile 300 lattice spacings away from the inlet after 20 time steps. Whatever you do at the inlet has not propagated that far after 20 time steps.

2) You initialize the velocity distributions functions with zero. That is quite unusual to say the least. Better use equilibrium values based on the macroscopic quantities.

[hamed.majeed]

Dear Alex,

Yes, when I changed the total time to 20000, the values turn out to be okay. And you were right about the whole domain setting to zero. Thank you.

Regards

[hamed.majeed]

The FORTRAN file is available at:
http://1drv.ms/1SA7J5o

[omertar]

Hello , how do you choose the density and the nu?


Sent from my iPhone using CFD Online Forum mobile app

[souai]

I'm working on the simulation of conduction and or convection radiation heat transfer in a 2D parallel channel by lattice Boltzmann method. when I tried to validate our code with a local nusselt number I have a disruption in some lasts values of nusselt. so I appreciate that you help us to explain or correct how to compute the nusselt number with LBM.

[FARKACH YOUNES]

You can find it in this book https:
//drive.google.com/file/d/15ZipkvhgK7MOP-RvoHNWQJ1w1ehTKgAS/view?usp=sharing