[Fortran] 纯文本查看 复制代码
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)
real ssum, usum, vsum
integer i
open(2,file='uvfield')
open(3,file='uvely')
open(4,file='vvelx')
open(8,file='timeu')
uo=0.20
sumvelo=0.0
rhoo=5.00 ! 虽然不知道为什么在后边密度的计算没有加上rhoo
dx=1.0
dy=dx
dt=1.0
alpha=0.02
Re=uo*m/alpha
print *, "Re=", Re
omega=1.0/(3.*alpha+0.5)
mstep=40000 ! 时间步长的确定需要考虑
cx(:)=(/0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0/)
cy(:)=(/0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0/)
w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)
! 初始条件的确定 对于初始条件的设定，是流体没有进入空间之前的设定
do j=21,m
do i=0,40
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do j=0,m
do i=41,n
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
! 边界条件设定 这个是流体没有进入前，也就是流体状态的设定
do j=21,m
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
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 j=21,m ! j=20是代表障碍物，21代表了流体
do i=0,40
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
END DO
END DO
do j=0,m
do i=41,n
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
END DO
END DO
return
end
subroutine streaming(f,n,m)
real f(0:8,0:n,0:m)
! 没有障碍物的流情况 这一部分存在着命名的重叠，先写后改 这一部分的节点是有问题的
! 存在障碍物的部分流
do j=22,m
do i=40,1,-1
f(1,i,j)=f(1,i-1,j)
end do
do i=0,39
f(3,i,j)=f(3,i+1,j)
end do
end do
do j=40,22,-1
do i=0,40
f(2,i,j)=f(2,i,j-1)
end do
do i=40,1,-1
f(5,i,j)=f(5,i-1,j-1)
end do
do i=0,39
f(6,i,j)=f(6,i+1,j-1)
end do
end do
do j=21,39
do i=0,40
f(4,i,j)=f(4,i,j+1)
end do
do i=0,39
f(7,i,j)=f(7,i+1,j+1)
end do
do i=40,1,-1
f(8,i,j)=f(8,i-1,j+1)
end do
end do
DO j=0,m
DO i=n,42,-1 !RIGHT TO LEFT
f(1,i,j)=f(1,i-1,j)
END DO
DO i=41,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=41,n
f(2,i,j)=f(2,i,j-1)
END DO
DO i=n,42,-1
f(5,i,j)=f(5,i-1,j-1)
END DO
DO i=41,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=41,n
f(4,i,j)=f(4,i,j+1)
END DO
DO i=41,n-1
f(7,i,j)=f(7,i+1,j+1)
END DO
DO i=n,42,-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)
real uo,rhow
! 有障碍物和没有障碍物的底边条件 这一部分没验证输出数据
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
! 存在障碍物北边界条件
do i=0,40
f(2,i,20)=f(4,i,20)
f(5,i,20)=f(7,i,20)
f(6,i,20)=f(8,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
! 左边没有障碍物的已知速度边界条件
do j=21,m
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)
f(1,0,j)=f(3,0,j)+2.*rhow*uo/3.
f(5,0,j)=f(7,0,j)+rhow*uo/6.
f(8,0,j)=f(6,0,j)+rhow*uo/6.
end do
! 上方边界
do i=0,n
f(4,i,m)=f(2,i,m)
f(7,i,m)=f(5,i,m)
f(8,i,m)=f(6,i,m)
end do
! 没有障碍物的右边界
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
return
end
subroutine rhouv(f,rho,u,v,cx,cy,n,m)
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)
integer :: i, j, k
real :: ssum, usum, vsum
! 有障碍物的一部分
do j=21,m
do i=0,40
ssum=0.0
do k=0,8
ssum=ssum+f(k,i,j)
end do
rho(i,j)=ssum
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)
!print*,usum! 为什么？？？？？？
end do
u(i,j)=usum/rho(i,j)
v(i,j)=vsum/rho(i,j)
!print*,u(i,j),v(i,j)
end do
end do
! 这一部分错了
do j=0,m
do i=41,n
ssum=0.0
usum=0.0
vsum=0.0
do k=0,8
ssum=ssum+f(k,i,j)
end do
rho(i,j)=ssum
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
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),rho(0:n,0:m),strf(0:n,0:m)
integer :: i, j
real :: rhoav, rhom
open(5, file='streamf')
strf(0,0)=0.
do i=1,n
rhoav=0.5*(rho(i-1,0)+rho(i,0))
strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
end do
do i=0,40
do j=21,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
do i=41,n
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,*)"VARIABLES =X, Y, U, V, S"
write(2,*)"ZONE ","I=",n+1,"J=",m+1,",","F=BLOCK"
do j=0,m
do i=0,n
write(2,*)i,j,u(i,j),v(i,j),strf(i,j)
end do
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
end subroutine result