はじめに
1月1日の能登半島地震では多くの被害が出たことと思います。
一刻も早く日常に戻れることをお祈りしています。
防災を生業とするエンジニアとして、今回の地震がどんな地震だったか調査したので結果をお伝えします。今回使用したデータは防災科学技術研究所と神戸海洋気象台のデータを使わせていただきました。関係各位にお礼申し上げます。
なお、今回は速報のため、間違いがありましたらご指摘ください。
加速度波形
図は今回の地震で最も加速度が大きかったK-NET富来のEW（東西方向）成分の加速度波形です。ピークが複数あり、ニュースで言われているような群発地震であることがわかります。最大加速度が2700gal程度あり、東北地方太平洋沖地震と同程度の地震と思われます。
加速度応答スペクトル
加速度応答スペクトルは各周期ごとの最大加速度応答をプロットしたものです。建物には固有周期があり、固有周期と地震動の最大加速度応答が大きい周期が一致すると、被害が大きくなります。
上図の青線が先程のK-NET富来の結果です。比較しているのは、兵庫県南部地震の神戸海洋気象台観測記録（赤）、東北地方太平洋沖地震のK-NET築館の観測記録（緑）、熊本地震のK-NET益城の観測記録（オレンジ）です。富来と築館は短周期側が非常に大きくなっている海溝型地震の特徴を示しています。一方、神戸海洋気象台とK-NET益城は周期0.2〜0.5程度が大きくなる内陸型地震の特徴を示しています。建物の固有周期は一般家屋で0.5程度で高層ビルになるほど固有周期は大きくなります。よって、建物への影響は兵庫県南部地震や熊本地震の方が大きいといえます。
非線形地震応答解析
地震応答解析はその名のとおり地震動に対して運動方程式を解き、構造物がどのような応答をするか調べるものです。ここでは、下図のようなモデルで降伏震度を設定する非線形解析を行います。降伏震度を超えると残留変位が発生します。今回は木造建築相当を想定し、固有周期0.5sec、建物の減衰率5%、降伏震度0.4復元力モデルにはバイリニアモデルを用いました。
上図の結果から、変位は小さいですが、降伏震度を超え残留変位が生じていることがわかります。最初のピークで降伏震度を超えており、それ以降のピークでは残留変位は大きくなっていないようです。以上より、家屋などを想定すると、地震動による被害は兵庫県南部地震や熊本地震の方が被害は大きいと思います。しかし、今回は火災と津波という兵庫県南部地震と東北地方太平洋沖地震の両方の要素があり、被害が大きくなったようです。
おわりに
もう少しすると詳細調査の結果が出てくると思います。その成果は今後に生かされるとして、今はとにかくこれ以上被害が出ないことを願っています。
最後に今回使用したPythonプログラムを載せておきます。
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import time
import math
from tqdm import tqdm
def nmk(beta, am, ac, ak, dt, ff, acc1, vel1, dis1):
acc2 = (ff-ac*(vel1+0.5*dt*acc1)-ak*(dis1+dt*vel1 +
(0.5-beta)*dt*dt*acc1))/(am+ac*0.5*dt+ak*beta*dt*dt)
vel2 = vel1+0.5*dt*(acc1+acc2)
dis2 = dis1+dt*vel1+(0.5-beta)*dt*dt*acc1+beta*dt*dt*acc2
return acc2, vel2, dis2
def inputdt():
file = pd.read_csv("KOBEJMA.csv")
t = file['TIME'].values
acc = file['NS'].values
return t, acc
def main_sdof():
fout = open('output.csv', 'w')
print('Period', ',', 'ResAcc.', file=fout)
start = time.time()
# Definition of time and acceleration
t, acc = inputdt()
ainp = np.array(acc)
dt = 0.02
# Parameter setting
h = 0.05
am = 1.0
f = -am*ainp
# Numerical integration
md = 2
ddt = dt/float(md)
beta = 1.0/6.0
rac = np.zeros(len(t))
vel = np.zeros(len(t))
dis = np.zeros(len(t))
T = np.zeros(100)
AMAX = np.zeros(100)
acc1 = 0.0
vel1 = 0.0
dis1 = 0.0
pa = math.log10(0.1)
pb = math.log10(5.0)
pc = (pb-pa)/99.0
pd = pa
for k in tqdm(range(0, 100)):
time.sleep(0.01)
pe = 10.0**pd
T1 = pe
pd = pd+pc
ak = 4.0*np.pi**2.0*am/T1**2.0
# print(ak)
ac = 2.0*h*np.sqrt(ak*am)
T[k] = T1
# print(ac)
for i in range(0, len(t)):
f1 = 0.0
if 1 <= i:
f1 = f[i-1]
f2 = f[i]
fm = f1+(f2-f1)
acc2, vel2, dis2 = nmk(beta, am, ac, ak, ddt, fm, acc1, vel1, dis1)
acc1 = acc2
vel1 = vel2
dis1 = dis2
# fout=open(fnameW,'a')
# print(acc1,file=fout)
# fout.close()
rac[i] = acc1+ainp[i]
vel[i] = vel1
dis[i] = dis1
AMAX[k] = np.max(rac)
print(T[k], ',', AMAX[k], file=fout)
print('{0:10.3f}(sec)'.format(time.time()-start))
if __name__ == '__main__':
main_sdof()
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import time
import math
def resilience(fy, delty, ak, gm, dis1, ndis2, qy1):
ql = qy1+ak*(ndis2-dis1)
qu = fy+gm*ak*(ndis2-delty)
qd = -fy+gm*ak*(ndis2+delty)
if ql > qu:
qy2 = qu
elif ql < qd:
qy2 = qd
else:
qy2 = ql
return qy2
def osmethod(fy, delty, gm, am, ac, ak, dt, fm, acc1, vel1, dis1, ndis1, qy1):
dndis = dt*vel1+dt*dt*acc1/4.0
ndis2 = dis1+dndis
qy2 = resilience(fy, delty, ak, gm, dis1, ndis2, qy1)
dnf = qy2-qy1-ak*(ndis2-ndis1)
ddis = ((am*(acc1+4.0*vel1/dt)+ac*vel1)-qy1-ak *
(dis1-ndis1)-dnf+fm)/(4.0*am/dt/dt+2.0*ac/dt+ak)
dvel = -2.0*vel1+2.0*ddis/dt
dacc = -2.0*acc1-4.0*vel1/dt+4.0*ddis/dt/dt
dis2 = dis1+ddis
vel2 = vel1+dvel
acc2 = acc1+dacc
return ndis2, qy2, dis2, vel2, acc2
def input():
file = pd.read_csv("Noto_EW.csv")
t = file['Time'].values
acc = file['Acc.'].values
t_c = np.zeros(len(t) * 10)
acc_c = np.zeros(len(acc) * 10)
for i in range(0, len(t)):
for j in range(0, 10):
k = i * 10 + j
t_c[k] = 0.001 * k + 0.001
acc_c[k] = acc[i] / 10.0 * j
return t_c, acc_c
def main_sdof():
start = time.time()
# Definition of time and acceleration
t, acc = input()
ainp = np.array(acc)
dt = 0.001
# Parameter setting
g = 9.80665 # Gravitational acceleration
h = 0.05 # Damping factor
am = 100.0 # Mass
T = 0.5 # Natural period
gm = 0.02 # Stiffiness deterioration
ak = 4.0*np.pi**2.0*am/T**2.0 # Initial stiffiness
ac = 2.0*h*np.sqrt(ak*am) # Damping coefficient
khy = 0.4 # Yield seismic intensity
fy = am*g*khy # Yield stress
delty = fy/ak # Yield displacement
# Numerical integration
rac = np.zeros(len(t))
vel = np.zeros(len(t))
dis = np.zeros(len(t))
ndis = np.zeros(len(t))
qy = np.zeros(len(t))
acc1 = 0.0
vel1 = 0.0
dis1 = 0.0
ndis1 = 0.0
qy1 = 0.0
for i in range(0, len(t)):
fm = -am*ainp[i]*0.01
ndis2, qy2, dis2, vel2, acc2 = osmethod(
fy, delty, gm, am, ac, ak, dt, fm, acc1, vel1, dis1, ndis1, qy1)
acc1 = acc2
vel1 = vel2
dis1 = dis2
ndis1 = ndis2
qy1 = qy2
rac[i] = acc1+ainp[i]
vel[i] = vel1
dis[i] = dis1
ndis[i] = ndis1
qy[i] = qy1
print('{0:10.3f}(sec)'.format(time.time()-start))
# Show drawing
fig = plt.figure()
plt.subplot(2, 1, 1)
plt.plot(t, ainp, linewidth=0.5, color='black')
plt.ylabel('Acceleration (gal)')
plt.subplot(2, 1, 2)
plt.plot(t, dis, linewidth=0.5, color='red')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.show()
fig = plt.figure()
plt.subplot(1, 1, 1)
plt.plot(dis, qy, 'b-')
plt.xlabel('Displacement (m)')
plt.ylabel('Force (kN)')
plt.show()
if __name__ == '__main__':
main_sdof()