はじめに

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()