[Fortran] 纯文本查看 复制代码
program main !顶盖驱动boltzmann 程序
implicit none
real*8 dx,dy,Lx,Ly,dt,c,cs
integer Nx,Ny
real*8 Re,niu,U0,rho0
real*8,allocatable:: rho(:,:)
real*8,allocatable:: w(:),e(:,:),un(:,:,:),un_1(:,:,:)
real*8 tau
integer Q
real*8,allocatable::Fn(:,:,:),Fn_1(:,:,:)
integer i,j,k,n,n_max
real*8 feq,feq_n,du,u,err_e,err
!初始化111111111111111111111111111111111111111111111111111111111111111111111111111111111111
n_max=20000 !迭代次数
err_e=10**(-6) !终止误差
dx=1.0 !格子的流向长度
dy=1.0 !格子的垂向长度
Nx=256 !流向格子数量
Ny=256
Lx=dx*Nx !流向长度
Ly=dy*Ny
dt=1.0 !时间步长
c=dx/dt !格子速度
cs=c/sqrt(3.0) !格子声速
Re=1000.0 !雷诺数
U0=0.1 !顶盖速度
rho0=1.0 !初始密度
niu=U0*Lx/Re !运动粘性系数
tau=3.0*niu+0.5 !松弛时间
Q=9 !速度离散数量
allocate(rho(Nx,Ny)) !各节点的密度
allocate(w(Q)) !权系数
allocate(e(Q,2)) !离散速度
allocate(un(Nx,Ny,2)) !各节点迭代后的速度
allocate(un_1(Nx,Ny,2)) !各节点迭代前的速度
allocate(Fn(Nx,Ny,Q)) !迭代后的分布函数
allocate(Fn_1(Nx,Ny,Q)) !迭代前的分布函数
rho(:,:)=rho0
w(:)=(/4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36/)
!open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
!write(10,*) w !正常情况下应当屏蔽*******************************************************************
e(:,1)=(/ 0,1,0,-1,0,1,-1,-1,1 /)
e(:,2)=(/ 0,0,1,0,-1,1,1,-1,-1/)
!open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
!write(10,*) (e(i,:),i=1,Q,1) !正常情况下应当屏蔽***************************************************************
un(:,:,:)=0 !各节点速度
un(:,Ny,1)=U0
!open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
!write(10,*) ((un(i,j,1),i=1,Nx,1),j=1,Ny,1) !正常情况下应当屏蔽****************************************
rho(:,:)=0 !各节点密度
do i=1,Nx,1
do j=1,Ny,1
do k=1,Q,1
call sub_equi(rho(i,j),w(k),e(k,:),un(i,j,:),cs,feq)
Fn_1(i,j,k)=feq !分布函数初始化为平衡态
enddo
enddo
enddo
end
下面是我的子程序
subroutine sub_equi(rho,w,e,u,cs,feq)
!该程序用于求解Maxwell 分布函数
!rho 密度
!w 权系数
!e 离散速度矢量
!u 格子宏观速度矢量
!cs 格子声速
!feq Maxwell 分布函数
implicit none
real*8 rho,w,cs,feq,eu,uu
real*8 e(2),u(2)
integer i
open(unit=10,file='inspection.txt') !正常情况下应当屏蔽*******************************************************************
write(10,*) u !正常情况下应当屏蔽***************************************************************
eu=e(1)*u(1)+e(2)*u(2) !向量e、u的点乘结果
uu=u(1)**2+u(2)**2 !向量u、u的向量点乘
feq=rho*w*(1+eu/cs**2+(eu/cs**2)**2/2-uu/(2*cs**2))
end