矩形障碍算例（附Fortran计算代码及MATLAB后处理代码）

! 我又来啦，做了个小算例，在100*100个格子的空腔中加了几个矩形障碍

!上代码！西边界为水流进入区

parameter (n=100,m=100)
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='uvfield.txt')
open(3,file='uvely.txt')
open(4,file='vvelx.txt')
open(8,file='timeu.txt')
uo=0.08          !有变动，设Re=400，H=100,L=100, v=0.01，密度1000（水）
sumvelo=0.0
rhoo=1000
dx=1.0
dy=dx
dt=1.0
alpha=0.01
Re=uo*m/alpha
print *, "Re=", Re
omega=1.0/(3.*alpha+0.5)
mstep=4000        !注意修改步数
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
cx(1)=1
cx(2)=0
cx(3)=-1
cx(4)=0
cx(5)=1
cx(6)=-1
cx(7)=-1
cx(8)=1
cy(0)=0
cy(1)=0
cy(2)=1
cy(3)=0
cy(4)=-1
cy(5)=1
cy(6)=1
cy(7)=-1
cy(8)=-1
do j=0,m
do i=0,n
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do j=1,m-1
u(0,j)=uo
v(0,j)=0.0
end do
! main loop
1 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)
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 the 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)
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.50*t2*t2-1.50*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
END DO