粉体シミュレーション入門｣に掲載の粒子要素法プログラムの訂正版"i-pem1.f"(このページで1999年4月21日から6月21日まで公開)にアルゴリズム上の間違いおよび 追加の訂正が見つかりました．再訂正個所と再修正プログラム"i-pem2.f"を以下に示します．ご利用の皆様には度々の修正申し訳ございません．ご検討いただき、御意見等ございましたらご連絡ください．ashimosa@mail.doshisha.ac.jp
訂正個所(1999/4/21)
- 接線方向のばね定数の設定に縦弾性定数とせん断弾性定数との比ｓｏを用いていましたが、Mindlin Force(接触部での局部的なすべりはないケースでの近似式)に基づいた設定へ変更
- 速度から変位増分への差分近似に用いる修正子の変更の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へ変更
再訂正個所１，２，は京大（人環研）早川氏、３，４，５，６ は北大（工）Boris Golman氏よりご指摘いただきました．有難うございました.．
特に２についてはアルゴリズム上の訂正ですので重要です。御検討ください．
訂正したプログラム
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