# Lattice Boltzman model in fortran

Hello,
I am new to lattice Boltzmann model . I am doing my thesis in this field and urgently need help. Currently I am working on the Shan Chan model . But I am finding it very difficult to implement the code. Here is what I have done till now . Any help will be highly appreciated

``````  integer          ly, lx
parameter        (ly=100, lx=100)
integer          i,j,k,l,in,ip,jn,jp
integer          nts
double precision g,tau,f1,f2,f3,rt0,rt1,rt2,temp
double precision ueqxij,ueqyij,uxsq,uysq,uxuy5,uxuy6
double precision uxuy7,uxuy8,usq
integer          is_solid_node(ly,lx)
double precision f(ly,lx,9),ftemp(ly,lx,9),feq(ly,lx,9)
double precision rho(ly,lx),u_x(ly,lx),u_y(ly,lx)
double precision x(ly,lx),y(ly,lx), density(ly,lx)
double precision ex(9),ey(9), psi(ly,lx), pressure(ly,lx),  t
``````

``````  open(unit=8,file='lbm.dat')
close(unit=8)
``````

• ``````.. initializing the nodes as fluid nodes ..
``````

``````  do 11 j = 1, ly, 1
do 12 i = 1, lx, 1
is_solid_node(i,j) = 0
``````

12 continue
11 continue

• ``````.. setting solid notes at walls on the top and bottom ..
``````

``````  do 13 i = 1,lx,1
is_solid_node(1,i) = 1
is_solid_node(ly,i) = 1
``````

13 continue

• ``````.. initializing the densities  ..
``````

``````  do 14 j = 1, ly/2, 1
do 15 i = 1, lx, 1
rho(j,i) = 500.0d0
``````

15 continue
14 continue

``````  do 141 j = ly/2, ly, 1
do 151 i = 1, lx, 1
rho(j,i) = 80.0d0
``````

151 continue
141 continue

• ``````.. setting the particle distribution fn ..
``````

``````  do 16 j = 1,ly,1
do 17 i = 1,lx,1
f(j,i,1) = (4.d0/9.0d0)*rho(j,i)
f(j,i,2) = (1.0d0/9.0d0)*rho(j,i)
f(j,i,3) = (1.0d0/9.0d0)*rho(j,i)
f(j,i,4) = (1.0d0/9.0d0)*rho(j,i)
f(j,i,5) = (1.0d0/9.0d0)*rho(j,i)
f(j,i,6) = (1.0d0/36.0d0)*rho(j,i)
f(j,i,7) = (1.0d0/36.0d0)*rho(j,i)
f(j,i,8) = (1.0d0/36.0d0)*rho(j,i)
f(j,i,9) = (1.0d0/36.0d0)*rho(j,i)
``````

17 continue
16 continue

• ``````.. defining velocity vectors ..
``````

``````  ex(1)= 0.0d0
ey(1)= 0.0d0
ex(2)= 1.0d0
ey(2)= 0.0d0
ex(3)= 0.0d0
ey(3)= 1.0d0
ex(4)=-1.0d0
ey(4)= 0.0d0
ex(5)= 0.0d0
ey(5)=-1.0d0
ex(6)= 1.0d0
ey(6)= 1.0d0
ex(7)=-1.0d0
ey(7)= 1.0d0
ex(8)=-1.0d0
ey(8)=-1.0d0
ex(9)= 1.0d0
ey(9)=-1.0d0
``````

• ``````.. starting the time loop ..
``````

``````  do 18 k = 1,nts,1
write(*, *) k
``````
• `````` .. computing macroscopic density, rho, and velocity
do 19 j = 1,ly,1
do 20 i = 1,lx,1
u_x(j,i) = 0.0d0
u_y(j,i) = 0.0d0

rho(j,i) = 0.0d0
if(is_solid_node(j,i).eq.0) then
do 21 l = 1,9,1
rho(j,i) = rho(j,i) + f(j,i,l)
u_x(j,i) = u_x(j,i) + ex(l)*f(j,i,l)
u_y(j,i) = u_y(j,i) + ey(l)*f(j,i,l)
``````

21 continue
u_x(j,i) = u_x(j,i)/rho(j,i)
u_y(j,i) = u_y(j,i)/rho(j,i)
psi(j,i) = 4.0d0*exp(-200.0d0/rho(j,i))

``````     endif
x(j,i) = dble(i)
y(j,i) = dble(j)
``````

20 continue
19 continue

• ``````.. Computing Interaction force ..
``````

``````   f1 = 3.0d0
f2 = 9.0d0/2.0d0
f3 = 3.0d0/2.0d0
do 99 j = 1,ly,1
do 22 i = 1,lx, 1
if(j.gt.1)then
jn = j-1
else
jn = ly
endif
if(j.lt.ly)then
jp = j+1
else
jp = 1
endif

if(i.gt.1)then
in = i-1
else
in = lx
endif
if(i.lt.lx)then
ip = i+1
else
ip = 1
endif

Fx = 0.0d0;
Fy = 0.0d0;

if(is_solid_node(j,i).eq.0) then
Fx = Fx + 1.0d0/9.0d0*ex(2)*psi(j ,ip)
Fy = Fy + 1.0d0/9.0d0*ey(2)*psi(j ,ip)
Fx = Fx + 1.0d0/9.0d0*ex(3)*psi(jp,i )
Fy = Fy + 1.0d0/9.0d0*ey(3)*psi(jp,i )
Fx = Fx + 1.0d0/9.0d0*ex(4)*psi(j ,in)
Fy = Fy + 1.0d0/9.0d0*ey(4)*psi(j ,in)
Fx = Fx + 1.0d0/9.0d0*ex(5)*psi(jn,i )
Fy = Fy + 1.0d0/9.0d0*ey(5)*psi(jn,i )
Fx = Fx + 1.0d0/36.0d0*ex(6)*psi(jp,ip)
Fy = Fy + 1.0d0/36.0d0*ey(6)*psi(jp,ip)
Fx = Fx + 1.0d0/36.0d0*ex(7)*psi(jp,in)
Fy = Fy + 1.0d0/36.0d0*ey(7)*psi(jp,in)
Fx = Fx + 1.0d0/36.0d0*ex(8)*psi(jn,in)
Fy = Fy + 1.0d0/36.0d0*ey(8)*psi(jn,in)
Fx = Fx + 1.0d0/36.0d0*ex(9)*psi(jn,ip)
Fy = Fy + 1.0d0/36.0d0*ey(9)*psi(jn,ip)

Fx = -G * psi(j,i) * Fx
Fy = -G * psi(j,i) * Fy
``````

• ``````.. Computing Equilibrium distribution ..
``````

``````      rt0 = (4.0d0/9.0d0)*rho(j,i)
rt1 = (1.0/9.0d0)*rho(j,i)
rt2 = (1.0d0/36.0d0)*rho(j,i)

ueqxij =  u_x(j,i)+ tau*Fx/rho(j,i)
ueqyij =  u_y(j,i)+ tau*Fy/rho(j,i)
uxsq   =  ueqxij * ueqxij
uysq   =  ueqyij * ueqyij
uxuy5  =  ueqxij + ueqyij
uxuy6  = -ueqxij + ueqyij
uxuy7  = -ueqxij - ueqyij
uxuy8  =  ueqxij - ueqyij
usq    =  uxsq + uysq

feq(j,i,1) = rt0*(1.0d0 - f3*usq)
feq(j,i,2) = rt1*(1.0d0 + f1*ueqxij + f2*uxsq - f3*usq)
feq(j,i,3) = rt1*(1.0d0 + f1*ueqyij + f2*uysq - f3*usq)
feq(j,i,4) = rt1*(1.0d0 - f1*ueqxij + f2*uxsq - f3*usq)
feq(j,i,5) = rt1*(1.0d0 - f1*ueqyij + f2*uysq - f3*usq)
feq(j,i,6) = rt2*(1.0d0 + f1*uxuy5 + f2*uxuy5*uxuy5-f3*usq)
feq(j,i,7) = rt2*(1.0d0 + f1*uxuy6 + f2*uxuy6*uxuy6-f3*usq)
feq(j,i,8) = rt2*(1.0d0 + f1*uxuy7 + f2*uxuy7*uxuy7-f3*usq)
feq(j,i,9) = rt2*(1.0d0 + f1*uxuy8 + f2*uxuy8*uxuy8-f3*usq)
endif
``````

22 continue
99 continue
do 23 j = 1,ly,1
do 24 i = 1,lx,1

• `````` .. bounce back at the solid walls ..
``````

``````    if(is_solid_node(j,i).eq.1)then
temp = f(j,i,2)
f(j,i,2) = f(j,i,4)
f(j,i,4) = temp
temp= f(j,i,3)
f(j,i,3) = f(j,i,5)
f(j,i,5) = temp
temp= f(j,i,6)
f(j,i,6) = f(j,i,8)
f(j,i,8) = temp
temp= f(j,i,7)
f(j,i,7) = f(j,i,9)
f(j,i,9) = temp
else
``````

• `````` .. regular collision at fluid nodes ..
``````

``````     do 25 l = 1,9,1
f(j,i,l) = f(j,i,l)-(f(j,i,l) - feq(j,i,l))/tau
``````

25 continue
endif
24 continue
23 continue

• `````` .. streaming steps ..
``````

``````   do 26 j= 1,ly,1
if(j.gt.1)then
jn = j-1
else
jn = ly
endif
if(j.lt.ly)then
jp = j+1
else
jp = 1
endif
do 27 i = 1,lx,1
if(i.gt.1)then
in = i-1
else
in = lx
endif
if(i.lt.lx)then
ip = i+1
else
ip = 1
endif
ftemp(j,i,1)  = f(j,i,1)
ftemp(j,ip,2) = f(j,i,2)
ftemp(jp,i,3) = f(j,i,3)
ftemp(j,in,4) = f(j,i,4)
ftemp(jn,i,5) = f(j,i,5)
ftemp(jp,ip,6) = f(j,i,6)
ftemp(jp,in,7) = f(j,i,7)
ftemp(jn,in,8) = f(j,i,8)
ftemp(jn,ip,9) = f(j,i,9)
``````

27 continue
26 continue

• ``````.. updating the f array ..
``````

``````   do 28 l = 1,9,1
do 29 j = 1,ly,1
do 30 i = 1,lx,1
f(j,i,l) = ftemp(j,i,l)
``````

30 continue
29 continue
28 continue

• ``````.. end of time loop ..
``````

18 continue

``````  open(unit=11,file='output.dat',status='unknown')
do 31 j = ly/2,ly,1
do 32 i = 1,lx,1
Pressure(j,i) = rho(j,i)/3.0d0 + G*psi(j,i)*psi(j,i)/6.0d0
write(11,*) j, pressure(j,i)
``````

32 continue
write(11,)
write(11,
)
31 continue
close(unit=12)

``````  open(unit=12,file='output1.dat',status='unknown')
do 33 j = ly/2,ly,1
do 34 i = 1,lx,1
Pressure(j,i) = rho(j,i)/3.0d0 + G*psi(j,i)*psi(j,i)/6.0d0
write(12,*) j, pressure(j,i)
``````

34 continue
write(12,)
write(12,
)
33 continue
close(unit=12)

``````  stop
end
``````