粉体工学会

「2次元2成分系球形粒子充填プログラム」の修正ならびに訂正

「粉体シミュレーション入門」に掲載の粒子要素法プログラムの訂正版"i-pem1.f"(このページで1999年4月21日から6月21日まで公開)にアルゴリズム上の間違いおよび 追加の訂正が見つかりました.再訂正個所再修正プログラム"i-pem2.f"を以下に示します.ご利用の皆様には度々の修正申し訳ございません.ご検討いただき、御意見等ございましたらご連絡ください.ashimosa@mail.doshisha.ac.jp

訂正個所(1999/4/21)

  1. 接線方向のばね定数の設定に縦弾性定数とせん断弾性定数との比soを用いていましたが、Mindlin Force(接触部での局部的なすべりはないケースでの近似式)に基づいた設定へ変更
  2. 速度から変位増分への差分近似に用いる修正子の変更の2点を加え,より近似の精度を上げました.
        subroutine init
line 157   commom/dpm/u(ni+3),v(ni+3),f(ni+3)                変位増分値の初期化を追加
line 167   do 5 i=1,n+3
    u(i)=0.d0
    v(I)=0.d0
    f(i)=0.d0
 5 continue

subroutine nposit
line 459   bv0=v0(i)                                     前ステップでの速度の記憶、追加
 bu0=u0(i)
 bf0=f0(i)

line462   v(i)=(1.5d0*v0(i)-0.5d0*bv0)*dt                                  修正子の変更 
 u(i)=(1.5d0*u0(i)-0.5d0*bu0)*dt                          
 f(i)=(1.5d0*f0(i)-0.5d0*bf0)*dt        

line 470   real(n)              dfloat(n)                                   倍精度へ変更

subroutine actf
line 495    .lt.                 .le.                                               変更 

line 510    kss=8.d0*b1*so*e/2.d0/(2.d0-po)           Mindlin Forceによるばね定数に変更
line 512    vss=vnn*dsqrt(kss/knn)
line 517    kss=8.d0*b1*so*e/(2.d0-po)
line 519    vss=vnn*dsqrt(kss/knn)       
   
line 532   kss *un               kss*us                                        変更

line 533   if (i.eq.4) then                                                        削除
          endif

line550    .gt.0               .gt.0.d0                                     実数へ変更

line 559   .lt.                .le.                                                 変更 
      

再訂正個所(1999/6/21)

1. subroutine actf 内 Line453

i-pem1.f
      en(i,jk)=knn*dis**1.5d0
  i-pem2.f
      en(i,jk)=en(i,jk)+knn*un

2. subroutine actf 内 Line470~

i-pem1.f
      if ((dabs(es(i,jk))-frc*en(i,jk)).gt.0.d0) then           
        es(i,jk)=frc*dsign(en(i,jk),es(i,jk))                   
        ds=0.d0
      endif              
        hn=en(i,jk)+dn                                          
        hs=es(i,jk)+ds                    
  i-pem2.f
      hn=en(i,jk)+dn                                          
      hs=es(i,jk)+ds                                                
      if ((dabs(hs)-frc*hn).gt.0.d0) then           
        hs=frc*dsign(hn,hs)                   
      endif 

3.subroutine fposit内

Line108     stop operatorの挿入
  Line119~   array qq(n)  初期化の追加

4.subroutine pcont内

Line303     stop operatorの挿入
  Line310~   if ((j.eq.0).or.(j.eq.i)) goto 80   へ変更

5.subroutine actf内

Line406      common /pos/x0(ni),z0(ni),qq(ni)  ブロックの削除
  Line426      if (j.le.n) then  へ変更

6.subroutine bfout内

Line536      array qq(ni)  の保存追加
  Line537      write(13,*) (u(i),v(i),f(i),i=1,n+3)    n+3へ変更

再訂正個所1,2,は京大(人環研)早川氏、3,4,5,6 は北大(工)Boris Golman氏よりご指摘いただきました.有難うございました..
特に2についてはアルゴリズム上の訂正ですので重要です。御検討ください.

訂正したプログラム

c********************************************************************
c***                                                              ***
c***           two-dimensional particle element method             ***  
c***              sphere model  ......   i-pem2.f                    ***
c***                                              june 18,1999    ***
c********************************************************************
      implicit real*8(a-h,k,m,o-z)                          
      parameter (ni=1000,nj=13,nc=20000)
      common /con/dt,fri,frw,e,ew,po,pow,so,g,de,pi
      common /wep/rr(ni),wei(ni),pmi(ni)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      common /pos/x0(ni),z0(ni),qq(ni)
      common /vel/u0(ni),v0(ni),f0(ni)      
      common /for1/xf(ni),zf(ni),mf(ni)        
      common /for2/en(ni,nj),es(ni,nj),je(ni,nj)     
      common /dpm/u(ni+3),v(ni+3),f(ni+3)         
      data maxstep/10000000/
c
c---------setting up the first position and initial condition
      call fposit(rmax)
c
      call inmat
c
c--------- initial setting------------------------------------------  
      call init
c
      t=0.d0
c
c---------iteration for each step-----------------------------------
      do 120 it=1,maxstep                                     
        t=t+dt 
c                                     
        call ncel
c                             
        do 100 i=1,n                                        
          xf(i)=0.d0                                           
          zf(i)=0.d0                                            
          mf(i)=0.d0
c  
c----------calculation of contact force between particle and wall---     
          call wcont(i)
c
c-----------calculation of contact force between particles---------- 
          call pcont(i,rmax)
c                                          
  100   continue 
c                                           
c----------superposition of incremental displacement ---------------
        call nposit(judge,it)
c
c----------judgement of static state
        if (1-judge) 200,200,199              
c
  199   if (mod(it,10000).eq.0) then
          write (6,*) 'time=',t,x0(1),z0(1),f0(1),f0(3)
        endif
c----------output of graphic data----------------------------------
        
        if (it.eq.1 .or. mod(it,10000).eq.0) then
          call gfout(it,t,rmax)
        endif
c                   
  120 continue                             
c---------output of back up data-----------------------------------
  200 call bfout(t,rmax)
c
      close(10)
      close(11)
      close(30)
c
      stop                                                     
      end                                                      
c  
c******* fposit *******************************************
c                                                             
      subroutine fposit(rmax)
      implicit real*8(a-h,k,m,o-z)
      parameter (ni=1000,nj=13,nc=20000)
      common /con/dt,fri,frw,e,ew,po,pow,so,g,de,pi
      common /wep/rr(ni),wei(ni),pmi(ni)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      common /pos/x0(ni),z0(ni),qq(ni)
      data ii/584287/
      data r1,r2/1.d-2,5.d-3/
      data w,ipz/6.006d-1,6/
c
      rmax=r1
      rmin=r2
      rn=rmax+1.d-5
      ipx=idint(w/2.d0/rn) 
      n=0
      do 20 i=1,ipz
        if (mod(i,2).eq.0)then
           dx=2.d0*rn
           ip=ipx-1
        else
           dx=rn
           ip=ipx
        endif
        do 10 j=1,ip
          call random(ii,ru)
          if (ru.lt.2.d-1) goto 10
          n=n+1
          if (n.gt.ni) then
            write(6,*) 'number of particles is more than ',ni
            stop
          endif
          x0(n)=2.d0*rn*(j-1)+dx
          z0(n)=2.d0*rn*(i-1)+rn 
c
          call random(ii,ru) 
          if (ru.lt.5.d-1) then
            rr(n)=r1
          else
            rr(n)=r2
          endif          
   10   continue
   20 continue 
      do 30 i=1,n
        qq(i)=0.d0
   30 continue
      write (6,*) 'number of particles:',n
      c=rmin*1.35d0
      idx=idint(w/c)+1
      idz=idint(z0(n)/c)+10
      if (idx*idz.gt.nc) then
        write(6,*) 'ncl is over!!',idx*idz
        stop
      endif
      write(6,*) idx,idz
      do 111 i= 1,n
       write(6,*) x0(i),z0(i)
  111 continue 
      return                                                  
      end
c
c******* inmat ******************************************
c                                              
      subroutine inmat
      implicit real*8(a-h,k,m,o-z)
      parameter (ni=1000,nj=13,nc=20000)
      common /con/dt,fri,frw,e,ew,po,pow,so,g,de,pi
      common /wep/rr(ni),wei(ni),pmi(ni)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      data dt/5.d-7/,pi/3.14159d0/,g/9.80665d0/
      data de,e,ew,po,pow/2.48d3,4.9d10,3.9d9,2.3d-1,2.5d-1/
      data fri,frw/2.5d-1,1.7d-1/
c
      so=1.d0/2.d0/(1.d0+po)
      do 10 i=1,n
        wei(i)=4.d0/3.d0*pi*rr(i)*rr(i)*rr(i)*de
        pmi(i)=8.d0/15.d0*de*pi*(rr(i))**5
   10 continue
c
      return
      end
c  
c******* init *******************************************     
c                                                             
      subroutine init
      implicit real*8(a-h,k,m,o-z)
      parameter (ni=1000,nj=13,nc=20000)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)       
      common /for2/en(ni,nj),es(ni,nj),je(ni,nj)
      common /dpm/u(ni+3),v(ni+3),f(ni+3)    
c
      do 3 i=1,n 
        do 2 j=1,13                                           
          en(i,j)=0.d0                                          
          es(i,j)=0.d0                                          
    2   continue                                             
        do 4 ij=1,13                                          
         je(i,ij)=0                                            
    4   continue                                              
    3 continue  
      do 5 i=1,n+3
        u(i)=0.d0
        v(i)=0.d0
        f(i)=0.d0
    5 continue                                    
      return                                                  
      end
c  
c******* ncel *******************************************    
c                                                             
      subroutine ncel
      implicit real*8(a-h,k,m,o-z)
      parameter (ni=1000,nj=13,nc=20000)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      common /pos/x0(ni),z0(ni),qq(ni)
      common /vel/u0(ni),v0(ni),f0(ni)      
      common /for1/xf(ni),zf(ni),mf(ni)        
      common /for2/en(ni,nj),es(ni,nj),je(ni,nj)     
      common /dpm/u(ni+3),v(ni+3),f(ni+3) 
c
      do 10 ib=1,(idx*idz)                                 
        ncl(ib)=0                                          
   10 continue                                              
      do 20 i=1,n                                           
         nncl(i)=0
         ib=idint(z0(i)/c)*idx+idint(x0(i)/c)+1
         ncl(ib)=i
         nncl(i)=ib                                             
   20 continue           
      return                                                  
      end
c  
c******* wcont *******************************************
c                                                             
      subroutine wcont(i)
      implicit real*8(a-h,k,m,o-z)
      parameter (ni=1000,nj=13,nc=20000)
      common /wep/rr(ni),wei(ni),pmi(ni)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      common /pos/x0(ni),z0(ni),qq(ni)
      common /vel/u0(ni),v0(ni),f0(ni)      
      common /for1/xf(ni),zf(ni),mf(ni)        
      common /for2/en(ni,nj),es(ni,nj),je(ni,nj)     
      common /dpm/u(ni+3),v(ni+3),f(ni+3)  
c
      xi=x0(i)
      zi=z0(i)
      ri=rr(i)
c
c  --- left wall ---
c
      jk=11                                                 
      j=n+1                                                 
      if (xi.lt.ri) then                                    
        as=0.d0                                              
        ac=-1.d0
        gap=dabs(xi)
        je(i,jk)=n+1                                               
        call actf(i,j,jk,as,ac,gap)                 
      else                                               
        en(i,jk)=0.d0                                         
        es(i,jk)=0.d0 
        je(i,jk)=0
      endif                                        
c 
c  --- under wall ---
c
      jk=12                                                 
      j=n+2                                                 
      if (zi.lt.ri) then
        as=-1.d0                                               
        ac=0.d0
        gap=dabs(zi)                                               
        je(i,jk)=n+2
        call actf(i,j,jk,as,ac,gap)
      else
        en(i,jk)=0.d0
        es(i,jk)=0.d0
        je(i,jk)=0
      endif
                             
c
c   --- right wall ---
c
      jk=13
      j=n+3
      if (xi+ri.gt.w) then
        as=0.d0
        ac=1.d0
        gap=dabs(xi-w)
        je(i,jk)=n+3
        call actf(i,j,jk,as,ac,gap)
      else
        en(i,jk)=0.d0
        es(i,jk)=0.d0
        je(i,jk)=0
      endif
      return                                                  
      end
c
c  
c******* pcont *******************************************
c                                                             
      subroutine pcont(i,rmax)
      implicit real*8(a-h,k,m,o-z)
      parameter (ni=1000,nj=13,nc=20000)
      common /wep/rr(ni),wei(ni),pmi(ni)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      common /pos/x0(ni),z0(ni),qq(ni)
      common /vel/u0(ni),v0(ni),f0(ni)      
      common /for1/xf(ni),zf(ni),mf(ni)        
      common /for2/en(ni,nj),es(ni,nj),je(ni,nj)     
      common /dpm/u(ni+3),v(ni+3),f(ni+3) 
c
      xi=x0(i)
      zi=z0(i)
      ri=rr(i)
c
      lup=idint((zi+2.d0*rmax)/c)                            
      lun=idint((zi-2.d0*rmax)/c)                            
      llf=idint((xi-2.d0*rmax)/c)                            
      lrg=idint((xi+2.d0*rmax)/c)                           
      if (lup.ge.idz) lup=idz-1                           
      if (lun.lt.0)   lun=0                                 
      if (llf.lt.0)   llf=0 
      if (lrg.ge.idx) lrg=idx-1                                 
c
      if(lup.le.lun) then
         write(6,*) i,'lup=',lup,'lun=',lun,'c=',c,
     &              'rmax=',rmax,xi,zi,idz
      stop
      endif
c
      do 90 lz=lun,lup                                    
        do 80 lx=llf,lrg                                    
         ib=lz*idx+lx+1                                      
         j=ncl(ib)                                           
         if ((j.eq.0).or.(j.eq.i)) goto 80                              
         do 11 jj=1,10                                      
           if (je(i,jj).eq.j) then                             
             jk=jj                                               
             goto 70                                             
           end if                                              
   11   continue                                            
        do 12 jj=1,10                                       
          if (je(i,jj).eq.0) then                             
            jk=jj                                               
            je(i,jj)=j                                          
            goto 70                                             
          end if                                              
   12   continue 
   70   xj=x0(j)                                            
        zj=z0(j)
        rj=rr(j)                                             
        gap=dsqrt((xi-xj)*(xi-xj)+(zi-zj)*(zi-zj))               
        if (gap.lt.(ri+rj)) then
          if (i.gt.j) then  
            ac=(xj-xi)/(gap)                          
            as=(zj-zi)/(gap)                                   
            j0=0
           do 555 jj=1,10
            if (je(j,jj).eq.i) then
                j0=jj
                goto 554
           endif
  555       continue   
  554       call actf(i,j,jk,as,ac,gap)                    
            en(j,j0)=en(i,jk)
            es(j,j0)=es(i,jk)
            j0=0
         endif
        else            
   85     en(i,jk)=0.d0                                       
          es(i,jk)=0.d0                                       
          je(i,jk)=0  
        endif                                        
   80   continue                                            
   90 continue  
      return                                                  
      end
c  
c******* nposit *******************************************
c                              
      subroutine nposit(judge,it)
      implicit real*8(a-h,k,m,o-z)
      parameter (ni=1000,nj=13,nc=20000)
      common /con/dt,fri,frw,e,ew,po,pow,so,g,de,pi
      common /wep/rr(ni),wei(ni),pmi(ni)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      common /pos/x0(ni),z0(ni),qq(ni)
      common /vel/u0(ni),v0(ni),f0(ni)      
      common /for1/xf(ni),zf(ni),mf(ni)        
      common /dpm/u(ni+3),v(ni+3),f(ni+3) 
c
      sum=0.d0                                            
      do 110 i=1,n
        bv0=v0(i)
        bu0=u0(i)
        bf0=f0(i)
        v0(i)=v0(i)+(zf(i)-wei(i)*g)/wei(i)*dt             
        u0(i)=u0(i)+xf(i)/wei(i)*dt             
        f0(i)=f0(i)+mf(i)/pmi(i)*dt                            
        v(i)=(1.5d0*v0(i)-0.5d0*bv0)*dt                           
        u(i)=(1.5d0*u0(i)-0.5d0*bu0)*dt                          
        f(i)=(1.5d0*f0(i)-0.5d0*bf0)*dt                         
        z0(i)=z0(i)+v(i)                                    
        x0(i)=x0(i)+u(i)                                    
        qq(i)=qq(i)+f(i)
        sum=dabs(u(i))+dabs(v(i))+sum
  110 continue                                            
      av=sum/dfloat(n)/2.d0 
      if (mod(it,10000).eq.0) then
         write(6,*) 'average= ',av
      endif                                
      if (av.lt.(dt*dt*g*1.0d-2)) then
        judge=1 
      else
        judge=0
      endif
            
      return                                                  
      end  
c  
c******* actf *******************************************                    
c                                                             
      subroutine actf(i,j,jk,as,ac,gap)              
      implicit real*8(a-h,k,m,o-z)
      parameter (ni=1000,nj=13,nc=20000)
      common /con/dt,fri,frw,e,ew,po,pow,so,g,de,pi
      common /wep/rr(ni),wei(ni),pmi(ni)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      common /vel/u0(ni),v0(ni),f0(ni)      
      common /for1/xf(ni),zf(ni),mf(ni)        
      common /for2/en(ni,nj),es(ni,nj),je(ni,nj)            
      common /dpm/u(ni+3),v(ni+3),f(ni+3) 
c
      ri=rr(i)
      if (j.le.n) then 
        rj=rr(j)
        dis=ri+rj-gap
        wei3=2.d0*wei(i)*wei(j)/(wei(i)+wei(j)) 
      else
        rj=0.d0
        dis=ri-gap
        wei3=wei(i)
      endif 
      enn=en(i,jk)
      if (enn.le.0.d0) then
        enn=1.d-3
      endif 
      if (j.le.n) then                                       
        b1=(3.d0/2.d0/e*ri*rj/(ri+rj)*(1.d0-po*po)
     &  *enn)**(1.d0/3.d0)
        knn=2.d0/3.d0*b1*e/(1.d0-po*po)
c        kss=knn*so
        kss=8.d0*b1*so*e/2.d0/(2.d0-po)
        vnn=dsqrt(2.d0*wei3*knn)
        vss=vnn*dsqrt(kss/knn)
      else
        b1=((3.d0/4.d0*ri*((1.d0-po*po)/e+(1.d0-pow*pow)/ew))
     &  *enn)**(1.d0/3.d0)
        knn=4.d0/3.d0*b1*e*ew/((1.d0-po*po)*ew+(1.d0-pow*pow)*e)
c        kss=knn*so
        kss=8.d0*b1*so*e/(2.d0-po)
        vnn=dsqrt(2.d0*wei3*knn)
        vss=vnn*dsqrt(kss/knn)
      endif
      ddt=1.d-1*dsqrt(wei3/knn)
      if (ddt.lt.1.d-6) write (6,*) 'over!! ddt=',ddt,i,j,jk,knn,wei(i)
      un=(u(i)-u(j))*ac+(v(i)-v(j))*as                        
      us=-(u(i)-u(j))*as+(v(i)-v(j))*ac+(ri*f(i)+rj*f(j))
c----------------------------------------------------------------------
      if(en(i,jk).eq.0.d0) then
        if (un.ne.0.d0) us=us*dis/un
        un=dis
      endif
c----------------------------------------------------------------------
      en(i,jk)=en(i,jk)+knn*un 
      es(i,jk)=es(i,jk)+kss*us
        
      dn=vnn*un/dt                                                             
      ds=vss*us/dt
      if (en(i,jk).lt.0.d0) then                               
        en(i,jk)=0.d0                                           
        es(i,jk)=0.d0                                           
        dn=0.d0                                                 
        ds=0.d0                                                 
        je(i,jk)=0                                                    
      return                                                  
      elseif ((j-n).le.0) then                                  
        frc=fri                                                 
      else                                                
        frc=frw
      endif 
      hn=en(i,jk)+dn                                          
      hs=es(i,jk)+ds                                                
      if ((dabs(hs)-frc*hn).gt.0.d0) then           
        hs=frc*dsign(hn,hs)                   
      endif              
                                       
        xf(i)=-hn*ac+hs*as+xf(i)                                
        zf(i)=-hn*as-hs*ac+zf(i)                                
        mf(i)=mf(i)-ri*hs 
        if (jk.le.10) then                                        
          xf(j)=hn*ac-hs*as+xf(j)                                
          zf(j)=hn*as+hs*ac+zf(j)                                
          mf(j)=mf(j)-rj*hs
        endif 
 
      return                                                  
      end  
c   
c******* gfout *******************************************
c                                                             
      subroutine gfout(it,t,rmax)
      implicit real*8(a-h,k,m,o-z)
      parameter (ni=1000,nj=13,nc=20000)
      common /wep/rr(ni),wei(ni),pmi(ni)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      common /pos/x0(ni),z0(ni),qq(ni)
      common /vel/u0(ni),v0(ni),f0(ni)            
      common /for2/en(ni,nj),es(ni,nj),je(ni,nj) 
c
      if (it.eq.1) then
        open(10,file='graph12.d')
        open(11,file='graph22.d')
        write(10,*) n,t,w,rmax
      else
        write(10,*) (sngl(x0(i)),sngl(z0(i)),sngl(rr(i)),i=1,n)
        write(10,*) (sngl(u0(i)),sngl(v0(i)),sngl(f0(i)),i=1,n)

        write(11,*) ((sngl(es(i,j)),sngl(en(i,j)),j=1,13),i=1,n)
        write(11,*) ((je(i,j),j=1,13),i=1,n)
      endif

      return                                                  
      end

c
c******* bfout *******************************************                    
c                                                             
      subroutine bfout(t,rmax)
      implicit real*8(a-h,k,m,o-z)                          
      parameter (ni=1000,nj=13,nc=20000)
      common /con/dt,fri,frw,e,ew,po,pow,so,g,de,pi
      common /wep/rr(ni),wei(ni),pmi(ni)
      common /cel/n,idx,idz,ipz,w,c,ncl(nc),nncl(ni)
      common /pos/x0(ni),z0(ni),qq(ni)
      common /vel/u0(ni),v0(ni),f0(ni)      
      common /for1/xf(ni),zf(ni),mf(ni)        
      common /for2/en(ni,nj),es(ni,nj),je(ni,nj)     
      common /dpm/u(ni+3),v(ni+3),f(ni+3)       
c      
      open(13,file='back200.d')
      write(13,*) n,idx,idz,ipz
      write(13,*) rmax,t,w,c,dt                                
      write(13,*) de,fri,frw,g,pi 
      write(13,*) e,ew,po,pow,so                      
      write(13,*) (wei(i),pmi(i),i=1,n)   
      write(13,*) (x0(i),z0(i),qq(i),rr(i),i=1,n)
      write(13,*) (u(i),v(i),f(i),i=1,n+3)
      write(13,*) (u0(i),v0(i),f0(i),i=1,n)
      write(13,*) ((es(i,j),en(i,j),j=1,13),i=1,n)
      write(13,*) ((je(i,j),j=1,13),i=1,n)
      close(13)

      return                                                  
      end

c
c
c******* random ******************************************                  
c                                              
      subroutine random(ii,ru)
      implicit real*8(a-h,k,m,o-z)
c
      ii=ii*65539
      if (ii.lt.0) ii=(ii+2147483647)+1
      ru=ii*0.4656613d-9
      return
      end