dsmc_to_couetee.f90

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

	!现在IC（2，I）已是第I网格第1个分子的序号-1
		ic2=ic2+ic(1,icell)
		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
		ic(1,ncell)=ic(1,ncell)+1
		k=ic(2,ncell)+ic(1,ncell)
		lcr(k)=m
	end do

	!此循环又数出了IC（1，I），并得到了原始****M的分子根椐其新位置的编号K。
	!将M置于LCR（K）中
	return
	end
!-------------------------------------------------------------------------------------!
	subroutine subc5(no_cell,no_molecule,coll_remain, ic, lcr,p,vmean)
	!计算碰撞
	implicit doubleprecision(a-h,o-z)
	common/cons/pi,boltz,eta51,y_lenght,dtm, &
			area,drag_ns,cell_height,rmass,tkom,amda,rkn
	common/init/ti_ini,pressure,fnd,t_wall,u_wall,vm_ini,vm_wall,vrm_ini,vrmax_eta51
	dimension p(4,no_molecule),vmean(3,no_molecule),lcr(no_molecule)
	dimension ic(2,no_cell),coll_remain(no_cell)
	dimension vrc(3),vccm(3)

	!dimension vrc(3),vccm(3),vrcip(3),vccmip(3)

190 do 140 m=1,no_cell
		if(ic(1,m).lt.2) goto 140
	!如果M网格中无分子或仅有一个分子，跳过碰撞
		vaver=0.0
	!VAVER是网格中CR**EAT51累计
		icontr=ic(1,m)
		if (ic(1,m).eq.2) icontr=1
		if (ic(1,m).eq.3) icontr=2
	!IConTR是需要随机选中分子对数，以计算网格中DTM内的碰撞数。
	!见7。2节，RSF
!
		do icoll=1,icontr
			k=int(rf(0)*ic(1,m)+ic(2,m)+0.9999999d0)
			if(k.eq.ic(2,m)) k=k+1
			l1=lcr(k) !在M网格中随机选中一分子l1
210			j=int(rf(0)*ic(1,m)+ic(2,m)+0.9999999d0)
			if(j.eq.ic(2,m)) j=j+1
			if(j.eq.k) goto 210
			l2=lcr(j) !!在M网格中随机选中另一分子l2

			do k=1,3
				vrc(k)=p(k,l1)-p(k,l2)
			end do

			vr=dsqrt(vrc(1)**2+vrc(2)**2+vrc(3)**2)
			vaver=vaver+vr**eta51
		enddo
	!vr**eta51*tkom=S(CR)

		vaver=vaver/icontr
		cnoic=ic(1,m)*(ic(1,m)/(area*cell_height))*tkom*vaver*dtm/2.0d0
!大括号为流场数密度n，tkom*vaver为S（CR）的平均值
!cnoic为DTM内每网格发生的碰撞数Nnsf
		cnoic=cnoic+coll_remain(m)
		ncoll=anint(cnoic)
		coll_remain(m)=cnoic-ncoll
		if(ncoll.lt.1) goto 140
		nacoll=0
	!nacoll是用来统计实际发生的碰撞数

300		k=int(rf(0)*ic(1,m)+ic(2,m)+0.9999999d0)
	 	if(k.eq.ic(2,m)) k=k+1
		l1=lcr(k)

310		j=int(rf(0)*ic(1,m)+ic(2,m)+0.9999999d0)
		if(j.eq.ic(2,m)) j=j+1
		if(j.eq.k) goto 310
		l2=lcr(j)
	 !从网格中随机选两个分子l1,l2

		do k=1,3
			vrc(k)=p(k,l1)-p(k,l2)
		enddo
		vr=dsqrt(vrc(1)**2+vrc(2)**2+vrc(3)**2)
	!vr为碰撞对的相对速度
		if(vr**eta51.gt.vrmax_eta51) vrmax_eta51=vr**eta51
		a=vr**eta51/vrmax_eta51
		b=rf(0)
		if(a.lt.b) goto 300
		!利用取舍法判断分子是否中选，（7-1），VR**ETA51在ETA51合适
		!赋值时，可民包括VHS，VSS模型
		b=1.0-2.0*rf(0)
		a=dsqrt(1.0-b*b)
		vrc(1)=b*vr
		b=2.0*pi*rf(0)
		vrc(2)=a*dcos(b)*vr
		vrc(3)=a*dsin(b)*vr
		!VRC为碰撞后的相对速度的三分量，见（2-111）（2-112），适合于HS和VHS模型
		do k=1,3
			vccm(k)=0.5*(p(k,l1)+p(k,l2))
			p(k,l1)=vccm(k)+vrc(k)*0.5
			p(k,l2)=vccm(k)-vrc(k)*0.5
			!VCCM为分子对的质心速 度分量，这里得到了L1，L2分子碰撞后的速度（2-41）
	!**	ve_macro=0.5*(vmean(k,l1)+vmean(k,l2))
	!** vmean(k,l1)=ve_macro
	!**	vmean(k,l2)=ve_macro
	!动量守恒并均分给两分子，（7-73）
		end do
		nacoll=nacoll+1
		if(nacoll.lt.ncoll) goto 300
		!重复计算碰撞直至碰撞数达到NRSF，（7-10），再进行下一网格计算
140	continue
	return
	end

!***************************************************************************************!

	subroutine subc6(no_cell,no_molecule,ic,lcr,p,vmean,field)
	!流场参数求和,(1)分子数(2)(3)宏观速度(4)(5)(6)分子速度平方累计
	implicit doubleprecision(a-h,o-z)
	dimension p(4,no_molecule),vmean(3,no_molecule),lcr(no_molecule), &
				ic(2,no_cell),field(7,no_cell)

	do icell=1,no_cell
		do molecule=1,ic(1,icell)
			k=ic(2,icell)+molecule
			l=lcr(k)
			field(1,icell)=field(1,icell)+1.0d0
			field(2,icell)=field(2,icell)+p(1,l)
!*			field(2,icell)=field(2,icell)+vmean(1,l)
			field(3,icell)=field(3,icell)+p(2,l)
!*			field(3,icell)=field(3,icell)+vmean(2,l)
			do ijk=1,3
				field(3+ijk,icell)=field(3+ijk,icell)+P(ijk,1)**2
			enddo
		enddo
	enddo
	return
	end
!*******************************************************************************************!

	subroutine subc7(nloop,jnis,no_cell,yc,no_molecule_each_cell,wall,up_wall,field,op)
	!输入表面参量，记入surf_cous.dat,输入流声参量，记u_coue.dat

	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 wall(7),up_wall(7),field(7,no_cell),op(8),yc(no_cell)

	open(26,file='surf_coue.dat')
!以下输出上表面数据
	skn=1.0d0/rkn
	sma=dsqrt(2.0d0/1.667d0)*u_wall/vm_ini
	write(26,10)skn,sma
10	format(1x,'knudsen number=',f6.3,'    Mach number=',f6.3/)
	write(26,18)
18	format(1x,'properties on upper surface'/)
	write(26,19)
19	format(1x,'      sample size number flux    pressure(in,re) &
			shear stress(in,re)         heat flux(in,re)')

	op(1)=wall(1)
	!(1)分子计数,样本大小，（nloop-jnis)*dtm=t时间内打到上来板有分子数
	a=(nloop-jnis)*dtm*(fnd*area*vm_ini)
	!t*n*A*vm
	op(2)=wall(1)/a
	!(2)用n*vm无量纲化的数目通量wall(1)/t*A
	op(3)=wall(3)/(vm_ini*a)
	op(4)=wall(6)/(vm_ini*a)
	

	op(5)=2.0*dsqrt(pi)*wall(2)/(u_wall*a)
	op(6)=2.0*dsqrt(pi)*wall(5)/(u_wall*a)

	op(7)=wall(4)/(vm_ini**2*a)
	op(8)=wall(7)/(vm_ini**2*a)


	write(26,21)op
21	format(1x,e13.5,7e15.5)
	drag=rmass*wall(2)/((nloop-jnis)*dtm*area)
	write(26,*)
	write(26,987)drag
987	format(1x,'     drag=',e10.4)


!c以下输出下表数据
	write(26,*)
	write(26,81)
81  format(1x,'properties on lower suface'/)
    write(26,19)
	op(1)=up_wall(1)
	!(1)分子计数,样本大小，（nloop-jnis)*dtm=t时间内打到上来板有分子数
	a=(nloop-jnis)*dtm*(fnd*area*vm_ini)
	!t*n*A*vm
	op(2)=up_wall(1)/a
	!(2)用n*vm无量纲化的数目通量wall(1)/t*A
	op(3)=up_wall(3)/(vm_ini*a)
	op(4)=up_wall(6)/(vm_ini*a)
	op(5)=2.0*dsqrt(pi)*up_wall(2)/(u_wall*a)
	op(6)=2.0*dsqrt(pi)*up_wall(5)/(u_wall*a)
	op(7)=up_wall(4)/(vm_ini**2*a)
	op(8)=up_wall(7)/(vm_ini**2*a)


	write(26,21) op

	drag=rmass*up_wall(2)/((nloop-jnis)*dtm*area)
	write(26,*)
	write(26,987)drag


	open(3,file='u_coue.dat')
	!以下输入流场数据
	write(3,123)
123 format(1x,'   N    y     sample size   rho     u     v   Tx &
			Ty     Tz       T')
	do 270 n=1,no_cell
	 op(1)=field(1,n) !网格中分子计算，样本大小
	 op(2)=op(1)/((nloop-jnis)*no_molecule_each_cell)!无量纲密度
	 op(3)=field(2,n)/op(1)
	 op(4)=field(3,n)/op(1)
	 !（3）（4）平均切向，法向速度，*平均切向，法向信息速度
	 op(5)=rmass*(field(4,n)/op(1)-op(3)**2)/boltz/t_ini
	 op(6)=rmass*(field(5,n)/op(1)-op(4)**2)/boltz/t_ini
	 op(7)=rmass*field(6,n)/op(1)/boltz/t_ini
	 !(5)(6)(7)xyz方向温度
	 op(8)=(op(5)+op(6)+op(7))/3
	 write(3,111) n,yc(n)/y_length,op
111	format(1x,I4,F9.4,F12.1,2F9.4,F10.4,4F9.4)
270	continue
	return
	end
!------------------------------------------------------------------------------------!
	function rf(idum) !产生随机数
	implicit doubleprecision(a-h,o-z)
	save ma,inext,inextp
	parameter(mbig=1000000000,mseed=161803398,mz=0,fac=1.0e-9)
	dimension ma(55)
	data iff/0/
	if (idum.lt.0.or.iff.eq.0) then
		iff=1
		mj=mseed-iabs(idum)
		mj=mod(mj,mbig)
		ma(55)=mj
		mk=1
	 do 50 I=1,54
		ii=mod(21*I,55)
		ma(ii)=mk
		mk=mj-mk
		if (mk.lt.mz) mk=mk+mbig
		mj=ma(ii)
50	 continue
	 do 100 k=1,4
		do 60 i=1,55
			ma(i)=ma(i)-ma(1+mod(i+30,55))
			if(ma(i).lt.mz) ma(i)=ma(i)+mbig
60		continue
100	 continue
	 inext=0
	 inextp=31
	endif
200	inext=inext+1
	if(inext.eq.56) inext=1
	inextp=inextp+1
	if(inextp.eq.56) inextp=1
	mj=ma(inext)-ma(inextp)
	if(mj.lt.mz) mj=mj+mbig
	ma(inext)=mj
	rf=mj*fac
	if(rf.gt.1.0e-8.and.rf.lt.0.99999999) return
	goto 200
	end