smdsam0.for

给予DSMC的二维CFD模拟FORTRAN代码

	SUBROUTINE SMDSAM(K,DDT,X,Y,Z,VX,VY,VZ)
	INCLUDE 'dsmc0.h'

*--SMSMCSS(N,M,L) sampled info. on the molecules striking the boundaries
*----M is the code number of the element; L is the species
*----N=1 the number sum
*----N=2 the sum of the normal momentum of the incident molecules
*----N=3 the sum of the normal momentum for the reflected molecules
*----N=4 the sum of the incident tangential momentum
*----N=5 the sum of the incident translational energy
*----N=6 the sum of the reflected translational energy
*----N=7 the sum of the incident rotational energy
*----N=8 the sum of the reflected rotational energy
*----N=9 the sum of the reflected tangential momentum
*
        VMP=SQRT(2.*BOLTZ*FTMP/SP(5,1))
	L =1
	SMT0=DDT
	x0=x
	Y0=Y
	Z0=Z
	vx0=vx
	vY0=vY
	vZ0=vZ

	VY=VY-VFY(K)
	CALL SURF_MD(K,DDT,X,Y,Z,VX,VY,VZ,DP)
CCC	if (abs(vx).gt.10.0d0*vmp)then
CCC	 write (*,*)'too large vx,at 3rd'
CCC	 write (*,*)k,'surface',nsurfcoll(k)
CCC	 write (*,*)x,y,z
CCC	 write (*,*)vx,vy,vz
CCC	 write (*,*)ddt
CCC	 pause
CCC	endif
	if (abs(vx).gt.3.0d0*vmp)then
	       NVELMOD(K)=	NVELMOD(K)+1
*  速度大得变态时，放弃MD计算，转而漫反射确定分子速度
             VX = SQRT(-LOG(RF(0)))*VMP
	       IF (K.EQ.2)VX=-VX
             CALL RVELC(VY,VZ,VMP)
	endif
	VY=VY+VFY(K)
	SMDDIS(K,1)=SMDDIS(K,1)+DP
	SMDDIS(K,2)=SMDDIS(K,2)+(Y-Y0)
	SMDDIS(K,3)=SMDDIS(K,3)+(Z-Z0)
	SMDTIME(K)=SMDTIME(K)+SMT0-DDT

	IF (DDT .GT. 0.0D0)THEN
*   假设分子不会连续两次与壁面碰撞
	X=X+VX*DDT
	Y=Y+VY*DDT
	Z=Z+VZ*DDT
       IF (K.EQ.1) THEN
        SMCSS(3,K,L)=SMCSS(3,K,L)+VX
       END IF
       IF (K.EQ.2) THEN
        SMCSS(3,K,L)=SMCSS(3,K,L)-VX
       END IF
       SMCSS(8,K,L)=SMCSS(8,K,L)+VY
       SMCSS(9,K,L)=SMCSS(9,K,L)+VZ
       SMCSS(6,K,L)=SMCSS(6,K,L)+0.5
     &           *(VX*VX+VY*VY+VZ*VZ)
	ENDIF


	RETURN
	END