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')
read(8,*)
read(8,*) g, tau, G
read(8,*)
read(8,*) nts
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