減衰定数ｈの１質点弾性系について考えます。
総入力エネルギーの速度換算値Veと
固有周期Tの関係を
エネルギースペクトルと定義しています。
それをfortranにより求めていきます。
real h,w,dt,beta,pi,T,Sacc(0:3000),Ve
real acc(0:5000),vel(0:5000),def(0:5000),E
integer i,j
open(1,file='KOBNS.csv',status='old')
h =0.1
dt =0.02
pi =3.14159265358979
do 110 i=0,1499
read(1,*) Sacc(i)
110 continue
acc(0)=-Sacc(0)
vel(0)=0. ![m/s]
def(0)=0. ![m]
beta=1./6.
do 200 j=1,100
T=j*0.05
w=2.*pi/T
E=0.
do 10 i=0,1499
acc(i+1)=-((Sacc(i)+(w**2.)*def(i)+(2.*h*w+(w**2.)*dt)*vel(i))
& +(h*w*dt+(0.5-beta)*(w**2.)*(dt**2.))*acc(i))
& /(1.+h*w*dt+beta*(w**2.)*(dt**2.))
el(i+1)=vel(i)+(acc(i+1)+acc(i))*dt/2.
def(i+1)=def(i)+vel(i)*dt+(0.5-beta)*acc(i)*(dt**2.)
& +beta*acc(i+1)*(dt**2.)
E=E-Sacc(i)*vel(i+1)*dt
10 continue
Ve=sqrt(E*2.)
open(2,file='report3.csv',status='unknown')
write(2,100) j*0.05,Ve
100 format(2f15.6)
200 continue
stop
end
KOBNS.csvという0.02秒刻みの地震波のデータを読み込んでいます。
総入力エネルギーの速度換算値を求める際、
質量ｍは約分されて無くなるので、
総入力エネルギーEには、上のプログラム上、ｍを入れていません。
また、本来、連続数として積分により求めますが、
入力地震派のデータが離散的なので、
その刻みに合わせて加算として求めています。
固有周期を0.05刻みでデータを求めてます。
横軸に固有周期T、縦軸に速度換算値Veとすれば、
エネルギースペクトルの図が出来上がります。