１質点系完全弾塑性モデルでの復元力特性を求めるfortranプログラム。
ご参考までに。
real h,w,dt,Sacc(0:3000),beta,pi,T
real acc(0:300000),vel(0:300000),def(0:300000)
real Q(0:3000000),m,k,k1,k2,Qy,dy(0:300000)
real Sacc2(0:300000),dt2,ka
integer i,j
open(1,file='KOBNS.csv',status='old')
dt =0.02
pi =3.14159265358979
beta=1./6.
T=1
h=0.
m=1.
w=2.*pi/T
k1=(w**2.)*m
k2=0
Qy=0.1*m*980. !α=0.1として
Sacc(0)=0.
dt2=dt/20. !/分割数
do 110 i=0,1449
read(1,*) Sacc(i)
110 continue
Sacc2(0)=Sacc(0)
do 200 i=0,1499
ka=(Sacc(i+1)-Sacc(i))/dt
do 210 j=0,19 !分割考慮
i2=i*20+j !分割考慮
Sacc2(i2+1)=ka*dt2+Sacc2(i2)
210 continue
200 continue
acc(0)=-Sacc(0)
vel(0)=0.
def(0)=0.
do 10 i=0,1449*20. !*分割数
acc(i+1)=-(Sacc(i+1)-Sacc(i)+((k/m)*dt2)*vel(i)
& +(h*w*dt2+(0.5-beta)*(k/m)*(dt2**2.)-1.)*acc(i))
& /(1.+h*w*dt2+beta*(k/m)*(dt2**2.))
vel(i+1)=vel(i)+(acc(i+1)+acc(i))*dt2*0.5
def(i+1)=def(i)+vel(i)*dt2+(0.5-beta)*acc(i)*(dt2**2.)
& +beta*acc(i+1)*(dt2**2.)
dy(i+1)=def(i+1)-def(i)
Q(i+1)=Q(i)+dy(i)*k
if(abs(Q(i+1)).lt.Qy) then
k=k1
else if(Q(i+1)*dy(i).lt.0.0)then
k=k1
else
k=k2
endif
open(12,file='report.csv',status='unknown')
write(12,100) def(i),Q(i+1)/100.
100 format(2f15.6)
10 continue
stop
end
記号や単位は察してください。（忘れちゃったので
おそらく、最後のif文がキモなのではないかと思います。
大まかに、弾塑性判定の方法は３つあって、
これは、せん弾力が、ある一定の値を超えたら、
その値を維持したまま剛性を０にするやり方です。
他には、せん弾力がある一定の値を超えたら、
その値に戻って、その値で剛性が０で変形させる、というやり方があります。
前者はエネルギー的な観点で見たとき、釣り合うのですが
後者は、値をいじっているので釣り合わない。
結果、エネルギーについて考えたとき、それがあっているのかどうか確認ができない。
ということで、前者の方法で判定しています。
あとは、収束するまで、地震動を細かく分割させるやり方。。。
そもそも、完全弾塑性モデルが、大雑把なモデルなので、
そのモデルで細かくやり過ぎても意味が無いと言う事で、やってません。
ちなみに、できたファイルをエクセル等でグラフ化すると
こんな感じになります。
参考文献