dsmc_to_couetee.f90

DSMC_to_couetee.rar是一个直接模拟的流体流动的源程序,用于计算couetee问题.

program  couette
! couette (DSMC,*IP)

	implicit doubleprecision (a-h,o-z)

	parameter (no_cell=50) !no_cell,网格数

	parameter (no_molecule_each_cell=30) !no_molecule_each_cell,每个网格分子数
	
	parameter (no_molecule=1500) !!no_molecule,分子总数

	parameter (jnis=1000,nloop=11000) !jnis_开始取样循环数，nloop,总循环数

	common/cons/pi,boltz,eta51,y_length,dtm,area,drag_ns, &
			cell_height,rmass,tkom,amda,rkn !see subc1
	common/init/t_ini,pressure,fnd,t_wall,u_wall,vm_ini,vm_wall,vrm_ini,vrmax_eta51 !see subc1

	dimension p(4,no_molecule),vmean(3,no_molecule),lcr(no_molecule)
	!p(i,n),i=1,2,3为分子速度的x,y,z方向分量，p(4,n)为分子的y方向位置
	!vmean(i,n)，i=1,2,3 is分子信息速度的x,y,z方向分量
	!lcr(n)中放置新位置编码为n的分子的原始号码

	dimension yc(no_cell),ic(2,no_cell),coll_remain(no_cell)

	!yc(n) is 第i格网格中心的坐标
	!ic(1,I)是第I网格的分子数
	!IC(2,I)是第I网格中第1分子的序号-1
	!COLL_REMAIN(I)是第I网格中DTM时间剩余的不足1的碰撞数

	dimension wall(7),up_wall(7),field(7,no_cell)

	!wall(I) is在下平板累计的:(1)入射分子数;(2)(5)入射和反射切向动量,
	!(3)(6)入射和反射法向动量;(4)(7)入射和反射动能
	!up_wall(I)是在上平板累计的以上各量
	!field(i,no_cell):流场量求和,(1,n)是第N网格中累计分子数累计
								!(2,n),(3,n)是第N网格中速度(*信息速度)1,2分量累计
								!(4,n),(5,n)(6,n)是第n网格中的速度分量平方累计
								
	dimension op(8)
	!op(8) 输出数据单元
	!分子类型Ar
	
	rkn=10.0d0 !knudsen倒数

	open(9,file='qcoul.dat')

	

	call subc1(no_cell,no_molecule,yc,coll_remain) !置常数与初值

	iloop=0 !循环数
	
	do ijk=1,7
		up_wall(ijk)=0.0d0
		wall(ijk)=0.0d0
			do icell=1,no_cell
				field(ijk,icell)=0.0d0
			enddo
	enddo

	

	call subc2(no_molecule,no_cell,no_molecule_each_cell,vm_ini, &
				cell_height,yc,p,vmean,y_length,u_wall) !置分初始速度

100	continue

	iloop=iloop+1

	if(mod(iloop,100).eq.0) write(*,*) 'iloop=',iloop
	
	
	
		call subc3(iloop,jnis,no_molecule,p,vmean,up_wall,wall) !计算分子运动和表面的反射
		call subc4(no_molecule,no_cell,cell_height,p,vmean,ic,lcr) !分子编号
		call subc5(no_cell,no_molecule,coll_remain,ic,lcr,p,vmean) !计算碰撞

		if(iloop.gt.jnis)then
		call subc6(no_cell,no_molecule,ic,lcr,p,vmean,field) !当前循环数大于开始循环数时,进行壁面和流场取样
		endif


	if (iloop.le.nloop) goto 100

	call subc7(nloop,jnis,no_cell,yc,no_molecule_each_cell,wall,up_wall,field,op)

	close(3)
	close(26)

24	format(1x,i4,f10.8,f15.5,7f10.5)
	stop
	end

!------------------------------------------------------------------------------------!
	subroutine subc1(no_cell,no_molecule,yc,coll_remain) !置常数与初值

	implicit doubleprecision(a-h,o-z)
	common /cons/pi,boltz,eta51,y_length,dtm,area,drag_ns, &
			cell_height,rmass,tkom,amda,rkn 
	common/init/t_ini,pressure,fnd,t_wall,u_wall,vm_ini,vm_wall,vrm_ini,vrmax_eta51 
	
	dimension yc(no_cell),coll_remain(no_cell)

	pi=3.1415926d0
	boltz=1.3805d-23 !boltzmann,(J/K)
	viscosity=2.117d-5 !NS/m^2-附表1,2
	tref=273.0
	rmass=6.63d-26 !(kg)原子质量
	diaref=3.659d-10 !(m) 分子半径,如用VHS,diaref=4.283d-10
	power=999999.9d0 !=eta(2-78),如用VHS,power=7.451d0
	eta51=(power-5.0d0)/(power-1.0d0) !????
	ksai=2.0d0/(power-1.0d0) !(2-99)????
	amr=rmass/2.0d0 !折合质量
	cxsref=pi*diaref**2 
	!P.292

	vhs_coe=1.0
	tkom=cxsref*(2.0*boltz*tref/amr)**ksai/vhs_coe


	t_ini=273.0d0
	pressure=101325d0
	fnd=pressure/(boltz*t_ini)

	t_wall=273.0d0
	u_wall=100.0d0
!*	u_wall=1.0d0

	vm_ini=dsqrt(2.0e0*boltz*t_ini/rmass)
	vm_wall=dsqrt(2.0e0*boltz*t_wall/rmass)
	vrm_ini=2.0*dsqrt(2.0/pi)*vm_ini

	vrmax_eta51=2.0d0*(vrm_ini)**eta51


	amda=viscosity*16.0d0/(5.0d0*dsqrt(pi)*rmass*fnd*vm_ini) !lambda,(2-201)
	write(*,*) 'amda=',amda
	y_length=rkn*amda
	dtm=0.23*amda/vm_ini
	!以上为流场范围内和时间步长

	area=no_molecule/(fnd*y_length)

	cell_height=y_length/dble(no_cell) !网格高度

	do n=1,no_cell
		yc(n)=(n-0.5d0)*cell_height
		coll_remain(n)=rf(0)
	enddo

	write(9,'(2x,4e16.6)') rmass,eta51,tkom,dtm,amda,y_length,area,cell_height
	write(9,'(2x,4e16.6)') t_ini,t_wall,pressure,fnd,u_wall,vm_ini,vm_wall

	return
	end

!---------------------------------------------------------------------------------------!
	subroutine subc2(no_molecule,no_cell,no_molecule_each_cell, &
					vm_ini,cell_height,yc,p,vmean,y_length,u_wall)
	! 
	! 置分子初始速度与位置
	!

	implicit doubleprecision (a-h,o-z)
	dimension p(4,no_molecule),wmean(3,no_molecule),yc(no_cell)

	pi=3.1415926d0

	do n=1,no_cell
		do m=1,no_molecule_each_cell
			iall=(n-1)*no_molecule_each_cell+m
				do ijk=1,3
					abc1=dsqrt(-dlog(rf(0)))
					abc2=2.0d0*pi*rf(0)
					p(ijk,iall)=abc1*dsin(abc2)*vm_ini
!**					vmean(ijk,iall)=0.0d0
				end do
!	这一循环产生模拟分子的初始分子速度，即其初始宏观速度分量加上平衡气体中热运动速度分量，
!   参见（3-18）和（3-19）及以下讨论。初始信息速度为0
			if(rf(0).gt.0.5d0)then
				p(4,iall)=yc(n)+0.5*rf(0)*cell_height
			else
				p(4,iall)=yc(n)-0.5*rf(0)*cell_height
!	给出模拟分子的在网格中的初始随机位置，这里考虑了减少方差原则，参见附录（III-5）
!	与（III-5）’，
			end if 
		end do
	end do
	return
	end

!----------------------------------------------------------------------------------!

	subroutine subc3(iloop,jnis,no_molecule,p,vmean,up_wall,wall)
!	计算分子运动和在上下表面的反射和取样
	implicit doubleprecision(a-h,o-z)
	common/cons/pi,boltz,eta51,y_length,dtm,&
	  area,drag_ns,cell_height, rmass,tkom,amda,rkn
	 common/init/t_ini,pressure,fnd,t_wall,u_wall,vm_ini, &
	   vm_wall,vrm_ini,vrmax_eta51
	dimension p(4,no_molecule),vmean(3,no_molecule)
	dimension wall(7),up_wall(7)
	do 180 m=1,no_molecule
		y=p(4,m)+p(2,m)*dtm
!	Y,模拟分子以P（2，M）运动经DTM时间后到达的位置
		if(y.ge.y_length)then
!如果Y大于Y_LENGHT，分子在上壁面反射
			dtr=(y-y_length)/p(2,m)
	!反射余的分子运动时间
			if(iloop.ge.jnis) then 
				p(1,m)=p(1,m)+u_wall/2.0d0
!**	vmean(1,m)=vneam(1,m)+u_wall/2.0d0
	!这里将参考系转换为与上平板（以-U_wal/2.0d0运动）相联系
				up_wall(1)=up_wall(1)+1.0d0
				up_wall(2)=up_wall(2)+p(1,m)
!*	up_wall(2)=up_wall(2)+vmean(1,m)
!在IP方法中UP_WALL（2）中累计的是入射分子的IP速度，见(7-75)
				up_wall(3)=up_wall(3)-p(2,m)
				up_wall(4)=up_wall(4)+0.5*(p(1,m)**2+p(2,m)**2+p(3,m)**2)
			end if 
! 当当前循环数ILOOP大于开始取样循环数JNIS时，进行了上平板取样（UP_WALL）
! 累计了入射分子贡献；（1）计算（2）切向动量（3）法向动量（4）动能。
			abc1=dsqrt(-dlog(rf(0)))*vm_wall
			abc2=2.0*pi*rf(0)
			p(1,m)=abc1*dsin(abc2)
			p(3,m)=abc1*dcos(abc2)
			p(2,m)=-dsqrt(-dlog(rf(0)))*vm_wall
!** vmean(1,m)=0.0d0
!** vmean(2,m)=0.0d0
!** vmean(3,m)=0.0d0

!   计算了在平板温反射后的分子的速度，参见（3-18）（3-19）（3-15）
!	（3-20）。反射后的IP速 度在与平板相联系的坐标系中为零。
			if(iloop.ge.jnis) then
				up_wall(5)=up_wall(5)+p(1,m)
				up_wall(5)=up_wall(5)+vmean(1,m)
!在IP方法中UP_wall(5)中累计的是反射分子的IP速度，见（7-75）
				up_wall(6)=up_wall(6)+p(2,m)
				up_wall(7)=up_wall(7)+0.5*(p(1,m)**2+p(2,m)**2+p(3,m)**2)
			end if 

!以上累计子反射分子贡献；（5）切向动量（为0），（6）法向动量，（7）动能

			p(1,m)=p(1,m)-u_wall/2.0d0 
!** vmean(1,m)=-u_wall/2.0d0
!** vmean(2,m)=0.0d0

			p(4,m)=y_length+p(2,m)*dtr

!以上从平板反射的切向分子速度加上了平板速度，切向和法向信息速度赋值。分子运动到新位置。

			goto 180
		end if 

		if(y.le.0.0d0) then
!如果Y小于0，分子在下壁面反射。分子的反射和取样与上板类似

			dtr=y/p(2,m)
			if(dtr.lt.1.0d-3*dtm) dtr=1.0d-3*dtm
				if(iloop.ge.jnis) then
					p(1,m)=p(1,m)-u_wall/2.0d0
!**  vmean(1,m)=vmean(1,m)-u_wall/2.0d0
					wall(1)=wall(1)+1.0d0
					wall(2)=wall(2)+p(1,m)
!*	wall(2)=wall(2)+vmean(1,m)
					wall(3)=wall(3)-p(2,m)
					wall(4)=wall(4)+0.5*(p(1,m)**2+p(2,m)**2+p(3,m)**2)
				end if 
				abc1=dsqrt(-dlog(rf(0)))*vm_wall
				abc2=2.0d0*pi*rf(0)
				p(1,m)=abc1*dsin(abc2)
				p(3,m)=abc1*dcos(abc2)
				p(2,m)=dsqrt(-dlog(rf(0)))*vm_wall
!** vmean(1,m)=0.0d0
!** vmean(2,m)=0.0d0
!** vmean(3,m)=0.0d0
				if (iloop.ge.jnis) then
					wall(5)=wall(5)+p(1,m)
!*	wall(5)=wall(5)+vmean(1,m)
					wall(6)=wall(6)+p(2,m) !missing
					wall(7)=wall(7)+0.5*(p(1,m)**2+p(2,m)**2+p(3,m)**2)
				end if 
				p(1,m)=p(1,m)+u_wall/2.0d0
!** vmean(1,m)=u_wall/2.0d0
!** vmean(2,m)=0.0d0
				p(4,m)=p(2,m)*dtr
				goto 180
			end if 
	p(4,m)=y
180 continue
    return
	end

!-------------------------------------------------------------------------------------!

	subroutine subc4(no_molecule,no_cell,cell_height,p,vmean,ic,lcr)
!原始号码为M的分子根据其运动后所在网格重新编号为K，M存于LCR（K）中
	implicit doubleprecision (a-h,o-z)
	dimension p(4,no_molecule),vmean(3,no_molecule),lcr(no_molecule),ic(2,no_cell)

150 do icell=1,no_cell
		ic(1,icell)=0
	end do

	do m=1,no_molecule
		ncell=p(4,m)/cell_height+0.9999999d0
		if(ncell.eq.0) ncell=1
		if(ncell.ge.no_cell) ncell=no_cell
 !ncell 是M分子所在的新的网格编号
		ic(1,ncell)=ic(1,ncell)+1
	end do 

	!这里数好了所有NO-CELL个网格中的分子数IC（1，I），为得到IC（2，I）做准备
	ic2=0

	do icell=1,no_cell
		ic(2,icell)=ic2