応答スペクトル

今回は応答スペクトルを計算するプログラムのご紹介です。言語はPython3とします。

トリパタイトに関する記事はこちら

簡単なサンプル

設定データ

応答スペクトル計算に関する設定データを用意します。対象とする固有周期の範囲減衰定数波形ファイルパスなどを設定します。
{
  "condition":
    {
      "period_begin" : 0.01,
      "period_end"   : 10.0,
      "period_point" : 1000,
      "damp_factor"  : 0.05,
      "dt"   : 0.02
    },
  "wave":
    {
      "path"  : "./samplewave.csv",
      "dt"    : 0.02,
      "factor": 0.01
    }
}
### 解説(各項目の説明) ###
{
  "condition":
    {
      "period_begin" : 0.01,    # 最小固有周期
      "period_end"   : 10.0,    # 最大固有周期
      "period_point" : 1000,  # 計算ポイント数
      "damp_factor"  : 0.05,    # 減衰定数
      "dt"   : 0.02             # 解析時間刻み
    },
  "wave":
    {
      "path"  : "./samplewave.csv", # 波形ファイルパス
      "dt"    : 0.02,                 # 波形刻み
      "factor": 0.01               # 波形倍率
    }
}

描画結果

ある地震波について応答スペクトルを計算して描画した結果を示します。

 

サンプルコード

上記のサンプルコードを以下に示します。設定ファイルに記述した地震波、減衰定数に対して固有周期毎にNigam・Jennings法(*1)を用いて応答計算を行い最大応答値をプロットしています。並列計算が可能な内容のため、マルチプロセスで高速化を図っています。描画はplotlyを利用しています。   *1:柴田明徳(2003)『最新 耐震構造解析 第2版』森北出版株式会社
# coding:utf-8
import sys, os, json, math
import numpy as np
import pandas as pd
from docopt import docopt
from scipy import interpolate
from concurrent.futures import ProcessPoolExecutor, Future, wait
import plotly.offline as of
import plotly.graph_objs as go

# 解析条件クラス
class AnalysisCondition:
    def __init__(self, config):
        self.period_begin = float(config['condition']['period_begin'])
        self.period_end = float(config['condition']['period_end'])
        self.period_point = int(config['condition']['period_point'])
        self.damp_factor = float(config['condition']['damp_factor'])
        self.dt = float(config['condition']['dt'])

# 波形クラス
class Wave:
    def __init__(self, config):
        self.path = config['wave']['path']
        self.dt = float(config['wave']['dt'])
        self.factor = float(config['wave']['factor'])

# 初期化
def init(config_file_path):
    with open(config_file_path, mode='r') as f:
        config = json.load(f)
    condition = AnalysisCondition(config)
    wave = Wave(config)
    return condition, wave

# 波形データ(加速度時刻歴)取得
def get_wave_data(condition, wave):
    with open(wave.path, 'r') as f:
        df = pd.read_csv(f)
    time = df['t'].values.tolist()
    acc = list(map(lambda x: x * wave.factor, df['acc'].values.tolist()))
    step = int((len(time)-1) * wave.dt / condition.dt)
    func_interpolate = interpolate.interp1d(time, acc) # 線形補間関数
    new_time = np.arange(0.0, step*condition.dt, condition.dt)
    new_acc = func_interpolate(new_time)
    return new_acc

# Nigam・Jennings法による応答計算
def nigam_jennings(period, acc0, dt, h):
    w = 2.0 * math.pi / period
    w_dt = w * dt
    dw = math.sqrt(1.0-h**2.0) * w
    dw_dt = dw * dt
    c = math.cos(dw_dt)
    s = math.sin(dw_dt)
    e = math.e ** (-1.0*h*w_dt)
    c0 = h / math.sqrt(1.0-h**2.0)
    c1 = c - c0 * s
    c2 = (2.0 * h**2.0 - 1.0) / (w**2.0 * dt)
    c3 = s / dw
    c4 = (2.0 * h) / (w**3.0 * dt)
    c5 = dw * s + h * w * c
    c6 = 1.0 / (w**2.0 * dt)
    c7 = 1.0 / (w**2.0)
    a11 = e * (c0 * s + c)
    a12 = e * s / dw
    a21 = -1.0 * w * e * s / math.sqrt(1.0-h**2.0)
    a22 = e * c1
    b11 = e * ((c2+h/w)*c3 + (c4+c7)*c) - c4
    b12 = -1.0 * e * (c2*c3 + c4*c) - c7 + c4
    b21 = e * ((c2+h/w)*c1 - (c4+c7)*c5) + c6
    b22 = -1.0 * e * (c2*c1-c4*c5) - c6
    dis = 0.0
    vel = 0.0
    acc = -1.0*(2.0*h*w*vel+w**2.0*dis)
    dis_max = 0.0
    vel_max = 0.0
    acc_max = 0.0
    for n in range(len(acc0)-1):
        d = a11*dis + a12*vel + b11*acc0[n] + b12*acc0[n+1]
        v = a21*dis + a22*vel + b21*acc0[n] + b22*acc0[n+1]
        a = -1.0 * (2.0*h*w*v + w**2.0*d)
        dis = d
        vel = v
        acc = a
        if abs(dis_max) < abs(d) : dis_max = abs(d)
        if abs(vel_max) < abs(v) : vel_max = abs(v)
        if abs(acc_max) < abs(a) : acc_max = abs(a)
    return period, dis_max, vel_max, acc_max

# 応答スペクトル計算
def spectrum(condition, wave):
    # 初期化
    df = pd.DataFrame([], columns=['period', 'acc', 'vel', 'dis'])
    periods = np.logspace(np.log10(condition.period_begin), np.log10(condition.period_end), num=condition.period_point)
    h = condition.damp_factor
    dt = condition.dt
    acc0 = get_wave_data(condition, wave)
    # 並列化
    with ProcessPoolExecutor(max_workers=os.cpu_count()) as executor:
        print('並列計算:CPU', os.cpu_count())
        futures = [executor.submit(nigam_jennings, period, acc0, dt, h) for period in periods]
        df['period'] = [future.result()[0] for future in futures]
        df['dis'] = [future.result()[1] for future in futures]
        df['vel'] = [future.result()[2] for future in futures]
        df['acc'] = [future.result()[3] for future in futures]
    return df

# 結果出力
def output(out_file_path, df):
    with open(out_file_path, 'w', newline='') as f:
        df.to_csv(f)

# グラフ描画
def draw(df):
    dis = go.Scatter(x=df['period'], y=df['dis'], name='Dis', showlegend=True, line=dict(color='green'))
    vel = go.Scatter(x=df['period'], y=df['vel'], name='Vel', showlegend=True, line=dict(color='blue'))
    acc = go.Scatter(x=df['period'], y=df['acc'], name='Acc', showlegend=True, line=dict(color='red'))
    fig_dis = dict(data=[dis], layout=dict(xaxis=dict(type='log', title=dict(text='固有周期(s)')), yaxis=dict(title=dict(text='応答変位(m)'))))
    fig_vel = dict(data=[vel], layout=dict(xaxis=dict(type='log', title=dict(text='固有周期(s)')), yaxis=dict(title=dict(text='応答速度(m/s)'))))
    fig_acc = dict(data=[acc], layout=dict(xaxis=dict(type='log', title=dict(text='固有周期(s)')), yaxis=dict(title=dict(text='応答加速度(m/s2)'))))
    of.plot(fig_dis)
    of.plot(fig_vel)
    of.plot(fig_acc)

if __name__ == "__main__":
    __doc__ = """
    Usage:
        main.py <config_file_path> <out_file_dir>
        main.py -h | --help
    Options:
        -h --help  show this help message and exit
    """
    args = docopt(__doc__)
    config_file_path = args['<config_file_path>']
    out_file_dir = args['<out_file_dir>']
    condition, wave = init(config_file_path)
    df_spectrum = spectrum(condition, wave)
    output(out_file_dir + 'out_spectrum.csv', df_spectrum)
    draw(df_spectrum)

まとめ

今回は応答スペクトル計算のご紹介でした。    

※本プログラムによるいかなるトラブル、損失、損害について弊社は一切責任を負いません。本記事に掲載している図やプログラムコードはあくまでも一例です。

返信を残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です