dsmc_airflowinmicro.f90

直接蒙特卡洛方法计算微管气体流动

program main

implicit real(kind=8)(a-h,o-z)
parameter (no_cellx=10,no_celly=5,no_cell=50)
parameter (no_molecule_each_cell=30)
parameter (no_molecule=1500)
parameter (nloop=10)

common /cons/ pi,boltz,eta51,x_length,y_length,dtm,cellx_height,celly_height,vrmax_eta51,vole,t_ini,v_max
common /init/ pressurein,pressureout,tkom,rmass
dimension p(5,no_molecule),lcr(no_molecule)
dimension yc(no_celly),xc(no_cellx),coll_remain(no_cell)
dimension field(7,no_cell)
dimension kgk(no_molecule) 
dimension js(2,no_cell)

call subc1(no_cellx,no_celly,no_cell,coll_remain,xc,yc,no_molecule)
call subc2(no_cellx,no_celly,no_molecule_each_cell,no_molecule,xc,yc,p)
iloop=1
do while(iloop<nloop)
  write(*,*) iloop
  call subc3(p,no_molecule,kgk,iloop)
  call subc4(no_molecule,no_cellx,no_celly,no_cell,kgk,lcr,p,xc,yc,js,iloop)
  call subc5(no_cell,no_molecule,js,p,coll_remain,lcr,iloop)
  call subc6(no_cell,no_molecule,js,p,field,no_cellx,no_celly,iloop)
  if(iloop>50)then
    if(mod(iloop,20)==0)then	
      call subc7(p,field,no_molecule,no_cell,js,no_celly,no_cellx,iloop)
	endif
  endif
  iloop=iloop+1
enddo
stop
end

subroutine subc1(no_cellx,no_celly,no_cell,coll_remain,xc,yc,no_molecule)
implicit real(kind=8)(a-h,o-z)
common /cons/ pi,boltz,eta51,x_length,y_length,dtm,cellx_height,celly_height,vrmax_eta51,vole,t_ini,v_max
common /init/ pressurein,pressureout,tkom,rmass
dimension yc(no_celly),xc(no_cellx),coll_remain(no_cell)
pi=3.1415926d0
boltz=1.3805d-23
viscosity=2.117d-5
tref=300.0d0
rmass=6.63d-26
diaref=3.659d-10
power=999999.9d0
eta51=(power-5.0d0)/(power-1.0d0)
aksai=2.0d0/(power-1.0d0)
amr=rmass/2.0d0
cxsref=pi*diaref**2
vhs_coe=1.0d0
tkom=cxsref*(2.0d0*boltz*tref/amr)**aksai/vhs_coe
t_ini=300.0d0
pressurein=101325d0
pressureout=0.3*101325d0
amda=boltz*t_ini/dsqrt(2.0d0)/pi/diaref**2/(pressurein*0.5)
v_max=dsqrt(2.0d0*boltz*t_ini/rmass)
vrmax_eta51=2.0d0*(v_max)**eta51

y_length=5*5.0d-8
x_length=5.0d-7

dtm=0.23d0*amda/v_max
cellx_height=x_length/dble(no_cellx)
celly_height=y_length/dble(no_celly)
open(15,file='changshu.txt')
write(15,54)cellx_height,celly_height,pressurein,dtm,v_max,amda
54 format(7e25.16)
close(15)
vole=pi*cellx_height*celly_height**2 
do m=1,no_cellx
   xc(m)=(m-0.5d0)*cellx_height
enddo
do n=1,no_celly
   yc(n)=(n-0.5d0)*celly_height
enddo
do mm=1,no_cell
	 coll_remain(mm)=rf(0)   
enddo
return
end

subroutine subc2(no_cellx,no_celly,no_molecule_each_cell,no_molecule,xc,yc,p)
implicit real(kind=8)(a-h,o-z)
common /cons/ pi,boltz,eta51,x_length,y_length,dtm,cellx_height,celly_height,vrmax_eta51,vole,t_ini,v_max
common /init/ pressurein,pressureout,tkom,rmass
dimension xc(no_cellx),yc(no_celly)
dimension p(5,no_molecule)
do i=1,no_cellx
  do j=1,no_celly
    do n=1,no_molecule_each_cell
	  iall=((i-1)*no_celly+j-1)*no_molecule_each_cell+n
	  do ijk=1,3
	    abc1=dsqrt(-dlog(rf(0)))
		abc2=2.0d0*pi*rf(0)
		if(ijk==1)then
		   p(ijk,iall)=abc1*dsin(abc2)*v_max+1.0d0
		else
		   p(ijk,iall)=abc1*dsin(abc2)*v_max
		endif
	  enddo
	  rfx=rf(0)
	  rfy=rf(0)
	  if(rfx>0.5d0)then
	    if(rfy>=0.5d0)then
		  p(4,iall)=xc(i)+0.5d0*rfx*cellx_height
		  p(5,iall)=yc(j)+0.5d0*rfy*celly_height
		else
		  p(4,iall)=xc(i)+0.5d0*rfx*cellx_height
		  p(5,iall)=yc(j)-0.5d0*rfy*celly_height
		endif
	  elseif(rfx<=0.5d0)then
	    if(rfy>0.5d0)then
		  p(4,iall)=xc(i)-0.5d0*rfx*cellx_height
		  p(5,iall)=yc(j)+0.5d0*rfy*celly_height
		else
		  p(4,iall)=xc(i)-0.5d0*rfx*cellx_height
		  p(5,iall)=yc(j)-0.5d0*rfy*celly_height
		endif
	  endif
	enddo
  enddo
enddo
return
end

subroutine subc3(p,no_molecule,kgk,iloop)
implicit real(kind=8)(a-h,o-z)
common /cons/ pi,boltz,eta51,x_length,y_length,dtm,cellx_height,celly_height,vrmax_eta51,vole,t_ini,v_max
common /init/ pressurein,pressureout,tkom,rmass
dimension kgk(no_molecule)
dimension p(5,no_molecule)
kgk(1:no_molecule)=0
do 180 m=1,no_molecule
       x=p(4,m)+p(1,m)*dtm
	   y=p(5,m)+p(2,m)*dtm
	   if(x<=0.or.x>=x_length)then
	     p(1:5,m)=0.0d0
		 kgk(m)=5
		 goto 180
	   endif
	   if(y>=y_length)then
	      dtr=(y-y_length)/p(2,m)
		  p(4,m)=p(4,m)+p(1,m)*(dtm-dtr)
		  abc1=dsqrt(-dlog(rf(0)))
		  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)))*v_max
		  p(4,m)=p(4,m)+p(1,m)*dtr
		  p(5,m)=y_length+p(2,m)*dtr
		  goto 180
		endif
		if(y<=0.0d0)then
		 if(kgk(m)-3.0d0<0)then
           p(1,m)=p(1,m)
		   p(2,m)=-p(2,m)
		   p(3,m)=p(3,m)
		   p(4,m)=x
		   p(5,m)=-y
           goto 180
		  endif
		endif 
		 p(4,m)=x
		 p(5,m)=y
180 continue
return
end	
	   
subroutine subc4(no_molecule,no_cellx,no_celly,no_cell,kgk,lcr,p,xc,yc,js,iloop)
implicit real(kind=8)(a-h,o-z)
common /cons/ pi,boltz,eta51,x_length,y_length,dtm,cellx_height,celly_height,vrmax_eta51,vole,t_ini,v_max
common /init/ pressurein,pressureout,tkom,rmass
dimension p(5,no_molecule),xc(no_cellx),yc(no_celly)
dimension lcr(no_molecule),ic(no_cellx,no_celly),kgk(no_molecule),js(2,no_cell) !,lcrr(no_molecule)
dimension pp(5,no_molecule),ppp(5,no_molecule)
dimension jss(2,no_cell) 
do i=1,no_celly
   do j=1,no_cellx
       ic(j,i)=0
   enddo
enddo
ncounter=0
do n=1,no_molecule
   if(kgk(n)-1.0d0>0)then
      ncounter=ncounter+1
   endif
enddo
do 110 m=1,no_molecule
   if(kgk(m)-1>0) goto 110
   ncellx=int(p(4,m)/cellx_height+0.99999999d0)
   ncelly=int(p(5,m)/celly_height+0.99999999d0)
   if(ncellx==0) ncellx=1
   if(ncelly==0) ncelly=1
   if(ncellx>=no_cellx) ncellx=no_cellx
   if(ncelly>=no_celly) ncelly=no_celly
   ic(ncellx,ncelly)=ic(ncellx,ncelly)+1
110 continue 
do i=1,no_cellx
   do j=1,no_celly
	  k=(i-1)*no_celly+j
	  js(1,k)=ic(i,j)
   enddo
enddo
ic2=0
do icell=1,no_cell
   js(2,icell)=ic2 
   ic2=ic2+js(1,icell)
   js(1,icell)=0
enddo
do i=1,no_celly
   do j=1,no_cellx
       ic(j,i)=0
   enddo
enddo

do 112 m=1,no_molecule
   if(kgk(m)-1>0) goto 112
   ncellx=int(p(4,m)/cellx_height+0.99999999d0)
   ncelly=int(p(5,m)/celly_height+0.99999999d0)
   if(ncellx==0) ncellx=1
   if(ncelly==0) ncelly=1
   if(ncellx>=no_cellx) ncellx=no_cellx
   if(ncelly>=no_celly) ncelly=no_celly
   ic(ncellx,ncelly)=ic(ncellx,ncelly)+1
   k=js(2,(ncellx-1)*no_celly+ncelly)+ic(ncellx,ncelly)   
   pp(1:5,k)=p(1:5,m)
112 continue
do i=1,no_cellx
   do j=1,no_celly
	  k=(i-1)*no_celly+j
	  js(1,k)=ic(i,j)
   enddo
enddo

if(ncounter==0) goto 1199

do k=1,ncounter
   jishu1=int(no_celly*rf(0)+0.99999999d0)
   jibh=js(2,jishu1)+1
   jifzs=no_molecule-ncounter+k-1
   do 99 i=1,jifzs
      if(i<jibh)then
	     ppp(1:5,i)=pp(1:5,i)
	     goto 99
	  endif
	  if(i>=jibh)then
         ppp(1:5,i+1)=pp(1:5,i)
		 goto 99
	  endif
   99 continue
   do 66 jj=1,no_cell
      if(jj<=jishu1) jss(2,jj)=js(2,jj)
      if(jj>jishu1)  jss(2,jj)=js(2,jj)+1
	  if(jj==jishu1)then
	      jss(1,jishu1)=js(1,jishu1)+1
		  goto 66
	  endif
	  jss(1,jj)=js(1,jj)
   66 continue
    do 2112 ijk=1,3
	    abc1=dsqrt(-dlog(rf(0)))
		abc2=2.0d0*pi*rf(0)
		if(ijk==1)then
		   ppp(ijk,jibh)=abc1*dsin(abc2)*v_max+5.0d0
		   goto 2112
		endif
		   ppp(ijk,jibh)=abc1*dsin(abc2)*v_max
   2112 continue
	rfx=rf(0)
	rfy=rf(0)
	if(rfx>0.5d0)then
	    if(rfy>=0.5d0)then
		  ppp(4,jibh)=xc(1)+0.5d0*rfx*cellx_height
		  ppp(5,jibh)=yc(jishu1)+0.5d0*rfy*celly_height
		else
		  ppp(4,jibh)=xc(1)+0.5d0*rfx*cellx_height
		  ppp(5,jibh)=yc(jishu1)-0.5d0*rfy*celly_height
		endif
	elseif(rfx<=0.5d0)then
	    if(rfy>0.5d0)then
		  ppp(4,jibh)=xc(1)-0.5d0*rfx*cellx_height
		  ppp(5,jibh)=yc(jishu1)+0.5d0*rfy*celly_height
		else
		  ppp(4,jibh)=xc(1)-0.5d0*rfx*cellx_height
		  ppp(5,jibh)=yc(jishu1)-0.5d0*rfy*celly_height
		endif
	endif
	do i=1,jifzs+1
	   pp(1:5,i)=ppp(1:5,i)
	enddo
	do i=1,no_cell
	   js(1,i)=jss(1,i)
	   js(2,i)=jss(2,i)
	enddo
enddo
do i=1,no_molecule
   p(1:5,i)=ppp(1:5,i)
enddo

1199 ic(1:no_cellx,1:no_celly)=0
do 212 m=1,no_molecule
   ncellx=int(p(4,m)/cellx_height+0.99999999d0)
   ncelly=int(p(5,m)/celly_height+0.99999999d0)
   if(ncellx==0) ncellx=1
   if(ncelly==0) ncelly=1
   if(ncellx>=no_cellx) ncellx=no_cellx
   if(ncelly>=no_celly) ncelly=no_celly
   ic(ncellx,ncelly)=ic(ncellx,ncelly)+1
   k=js(2,(ncellx-1)*no_celly+ncelly)+ic(ncellx,ncelly)   
   lcr(k)=m 
212 continue
return
end

subroutine subc5(no_cell,no_molecule,js,p,coll_remain,lcr,iloop)
implicit real(kind=8)(a-h,o-z)
common /cons/ pi,boltz,eta51,x_length,y_length,dtm,cellx_height,celly_height,vrmax_eta51,vole,t_ini,v_max
common /init/ pressurein,pressureout,tkom,rmass

dimension p(5,no_molecule),vrc(3),vccm(3)
dimension coll_remain(no_cell)  
dimension lcr(no_molecule),js(2,no_cell)
dimension area(no_cell)

190 do 140 m=1,no_cell 
      area(m)=dble(js(1,m))*boltz*t_ini/(pi*pressurein*y_length**2)
	  if(js(1,m)<2) goto 140
	    vaver=0.0d0
	    icontr=js(1,m)
		if(js(1,m)==2) icontr=1
		if(js(1,m)==3) icontr=2
	  do icoll=1,icontr 
		 k=int(rf(0)*js(1,m)+js(2,m)+0.99999999d0)
		 if(k==js(2,m)) k=k+1
		 l1=lcr(k)
        210 j=int(rf(0)*js(1,m)+js(2,m)+0.99999999d0)
		 if(j==js(2,m)) j=j+1
		 if(j==k) goto 210
		 l2=lcr(j)
		 do k=1,3
			vrc(k)=p(k,l1)-p(k,l2)
		 enddo
		 vr=dsqrt(vrc(1)**2+vrc(2)**2+vrc(3)**2)
	 	 vaver=vaver+vr**eta51
	  enddo
		vaver=vaver/icontr
		cnoic=js(1,m)*(js(1,m)/(area(m)*cellx_height*celly_height))*tkom*vaver*dtm/2.0d0  
		cnoic=cnoic+coll_remain(m)
		ncoll=anint(cnoic)
		coll_remain(m)=cnoic-ncoll
		if(ncoll<1) goto 140
		nacoll=0
	300 k=int(rf(0)*js(1,m)+js(2,m)+0.99999999d0)
		if(k==js(2,m)) k=k+1
		l1=lcr(k)
	310 j=int(rf(0)*js(1,m)+js(2,m)+0.99999999d0)
        if(j==js(2,m)) j=j+1
		if(j==k)goto 310
		l2=lcr(j)
		do k=1,3
			vrc(k)=p(k,l1)-p(k,l2)
		enddo
		vr=dsqrt(vrc(1)**2+vrc(2)**2+vrc(3)**2)
		if(vr**eta51>vrmax_eta51) vrmax_eta51=vr**eta51
		a=vr**eta51/vrmax_eta51
		b=rf(0)
		if(a<b) goto 300
		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
		do k=1,3
			vccm(k)=0.5d0*(p(k,l1)+p(k,l2))
			p(k,l1)=vccm(k)+vrc(k)*0.5d0
            p(k,l2)=vccm(k)-vrc(k)*0.5d0
		enddo
		nacoll=nacoll+1
		if(nacoll<ncoll) goto 300
	140 continue
do i=1,no_molecule
   p(4,i)=p(4,i)
   p(5,i)=p(5,i)
enddo
return
end

subroutine subc6(no_cell,no_molecule,js,p,field,no_cellx,no_celly,iloop)
implicit real(kind=8)(a-h,o-z)
common /cons/ pi,boltz,eta51,x_length,y_length,dtm,cellx_height,celly_height,vrmax_eta51,vole,t_ini,v_max
common /init/ pressurein,pressureout,tkom,rmass

dimension field(7,no_cell),p(5,no_molecule),density(no_cellx,no_celly)
dimension nn(no_celly),js(2,no_cell) 
dimension pressure(no_cellx,no_celly),tt(no_cellx,no_celly) 
dimension velocity(no_cellx,no_celly),dpressure(no_cellx,no_celly),sudu(3,no_cell)
character(len=20) :: ifilename,iform
field=0.0d0
do i=1,no_celly
   vole2=(2*i-1)*vole
   nn(i)=vole2*pressurein/(boltz*t_ini*js(1,i))
   do j=1,no_cellx
      icell=(j-1)*no_celly+i
      do molecule=1,js(1,icell)
         k=js(2,icell)+molecule
	     field(1,icell)=field(1,icell)+1.0d0*dble(nn(i)) 
	     field(2,icell)=field(2,icell)+p(1,k)*dble(nn(i))
	     field(3,icell)=field(3,icell)+p(2,k)*dble(nn(i))
	     field(4,icell)=field(4,icell)+p(3,k)*dble(nn(i))
	     do ijk=1,3
	       field(4+ijk,icell)=field(4+ijk,icell)+dble(nn(i))*p(ijk,k)**2  
	    enddo
	  enddo
	  sudu(1,icell)=field(2,icell)/field(1,icell)
	  sudu(2,icell)=field(3,icell)/field(1,icell)
	  sudu(3,icell)=field(4,icell)/field(1,icell)
	enddo
enddo

do i=1,no_celly
   vole2=(2*i-1)*vole
   nn(i)=vole2*pressurein/(boltz*t_ini*js(1,i))
   do j=1,no_cellx
      incell=(j-1)*no_celly+i   
         tt(j,i)=rmass*(field(5,incell)+field(6,incell)+field(7,incell))/field(1,incell)/boltz/3 
		 velocity(j,i)=3331.45d0*dsqrt(1+tt(j,i)/273.15d0)
		 pressure(j,i)=js(1,incell)*nn(i)/vole2*boltz*tt(j,i)
		 density(j,i)=pressure(j,i)/tt(j,i)/208.1303d0 
	enddo
enddo
do j=1,no_celly
    dpressure(1,j)=pressurein-pressure(1,j)
	do k=1,js(1,j)
	   kk=js(2,j)+k
	   p(1,kk)=p(1,kk)+dpressure(1,j)/(density(1,j)*velocity(1,j))
	enddo
enddo
do j=1,no_celly  
    dpressure(no_cellx,j)=pressure(no_cellx,j)-pressureout
	jj=no_celly*(no_cellx-1)+j
	do k=1,js(1,jj)
	   kk=js(2,jj)+k
	   p(1,kk)=p(1,kk)+dpressure(no_cellx,j)/(density(no_cellx,j)*velocity(no_cellx,j))
	enddo
enddo
do i=1,no_molecule
   do j=2,5
      p(j,i)=p(j,i)
   enddo
enddo

if(iloop>50)then
  if(mod(iloop,20)==0)then
    write(iform,'(i4)') iloop
    write(ifilename,*) "field",trim(iform),".txt"
    open(34,file=ifilename)
    do i=1,no_cell
      write(34,332)i,field(1,i),field(2,i),field(5,i)
      332 format(1x,i6,3e25.15)
    enddo
    close(34)
    write(iform,'(i4)') iloop
    write(ifilename,*) "sudu",trim(iform),".txt"
    open(78,file=ifilename)
    do i=1,no_cell
      write(78,241)i,sudu(1,i),sudu(2,i),sudu(3,i)
	  241 format(1x,i6,3e25.15)
    enddo
    close(78)
    write(iform,'(i4)') iloop
    write(ifilename,*) "pressure",trim(iform),".txt"
    open(13,file=ifilename)
    do i=1,no_cellx
      do j=1,no_celly
        write(13,536)i,j,pressure(i,j)
        536 format(1x,2i6,e25.15)
      enddo
    enddo
    close(13)
    write(iform,'(i4)') iloop
    write(ifilename,*) "ppppp",trim(iform),".txt"
	open(33,file=ifilename)
    do jj=1,no_molecule
      write(33,67)jj,p(1,jj),p(2,jj),p(3,jj),p(4,jj),p(5,jj)
      67 format(1x,i4,3f20.12,2e20.12)
    enddo
    close(33)
  endif
endif

return
end

subroutine subc7(p,field,no_molecule,no_cell,js,no_celly,no_cellx,iloop)
implicit real(kind=8)(a-h,o-z)
common /cons/ pi,boltz,eta51,x_length,y_length,dtm,cellx_height,celly_height,vrmax_eta51,vole,t_ini,v_max
common /init/ pressurein,pressureout,tkom,rmass

dimension p(5,no_molecule),field(7,no_cell),js(2,no_cell),ann(no_celly)
character(len=20) :: ifilename,iform
zhiliangliu=0.0d0
do i=1,no_celly
   vole2=(2*i-1)*vole
   ann(i)=vole2*pressurein/(boltz*t_ini*js(1,i))
   jj=(no_cellx-1)*no_celly+i
   zhiliangliu=zhiliangliu+field(1,jj)*ann(i)*rmass/celly_height
enddo
write(iform,'(i4)') iloop
write(ifilename,*) "zhiliangliu",trim(iform),".txt"
open(23,file=ifilename)
write(23,*) zhiliangliu
close(23)
return
end



		  
function rf(idum)
implicit real(kind=8)(a-h,o-z)
save ma,inext,inextp
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
dimension ma(55)
data iff/0/
if(idum<0.or.iff==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<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)<mz) ma(i)=ma(i)+mbig
    60	continue
100 continue
inext=0
inextp=31
endif
200  inext=inext+1
	 if(inext==56) inext=1
	inextp=inextp+1
	if(inextp==56) inextp=1
	mj=ma(inext)-ma(inextp)
	if(mj<mz) mj=mj+mbig
	ma(inext)=mj
	rf=mj*fac
	if(rf>1.e-8.and.rf<0.999999999) return
goto 200
end