目次

1. はじめに
 2. ライブラリのインストール
 3. 地震応答解析
 4. グラフの出力
 5. おわりに

はじめに

構造系ではGrasshopperが人気ですが、Dynamoでも簡単な応答解析ができないか試してみます。

ライブラリのインストール

Pythonで使う外部ライブラリをインストールします。

1. はじめに Dynamo を起動し、下記のコードで CPython のバージョンを確認します。

PythonScript

import sys
OUT = sys.version

2. Anaconda をダウンロードし、Anaconda Navigator から「CMD.exe Prompt」を実行します。

3. 下記のコードを入力して Dynamo 用の環境をつくります。
3.8.3 は最初に確認した CPython のバージョンの番号です。

Prompt

conda create --name Dynamo383 python=3.8.3

4. Proceed ([y]/n)?と聞かれるので Y を入力します。

5. いま作成した環境に切り替えます。

Prompt

conda activate Dynamo383

6. 必要なライブラリをまとめてインストールします。

Prompt

conda install numpy pandas scipy Pillow matplotlib

地震応答解析

せん断変形のみを考慮した多質点系モデルの振動解析を実装します。
波形データは東京都建設局で公開されているCSVデータを利用します。

条件値と波形データの入力

ライブラリをインポートします。

PythonScript

#ライブラリをロード
sys.path.append(r'C:\Users\shota\.conda\envs\Dynamo383\Lib\site-packages')
import sys, math, io
import numpy as np
import pandas as pd
from scipy import interpolate
import matplotlib.pyplot as plt
from PIL import Image as Img
import System.Drawing
from System.Drawing import *
from System.Drawing.Imaging import *
from System.IO import MemoryStream

コード中に IN[0] とかあると分かりにくいので入力値は Inputクラス にまとめておきます。

PythonScript

# インプット
class Input:
    def __init__(self):
        self.g = IN[0] # 重力加速度
        self.dt = IN[1] # 積分時間刻み
        self.beta = IN[2] # Newmarkβ法のβ
        self.damp_factor = IN[3] # 減衰定数
        self.natural_frequency = IN[4] # 振動数
        self.wave_path = IN[7] # CSVファイルのパス
        self.arr_m = IN[5] # 重量
        self.arr_k = IN[6] # 剛性

CSV形式の波形データを読み込みます。

PythonScript

# 波形データ
def WaveData(input):
    with open(input.wave_path, 'r') as f:
        df = pd.read_csv(f)
    return list(df['acc'].values.tolist())

マトリクス生成

質量マトリクス、減衰マトリクス、剛性マトリクスをつくります。

PythonScript

# 質量マトリクス
def MassMatrix(arr_m, input):    
    return np.diag(list(map(lambda x: x / input.g, arr_m)))
    
# 剛性マトリクス
def StiffMatrix(arr_k):
    size = len(arr_k)
    full = np.zeros((size, size))
    for i in range(size):
        k = arr_k[i]
        part = [[k , -k], [-k,   k]]
        for i_row in range(len(part)):
            for i_col in range(len(part)):
                target_row = i-1 + i_row
                target_col = i-1 + i_col
                if 0 <= target_row and 0 <= target_col:
                    full[target_row][target_col] += part[i_row][i_col]
    return full
    
# 減衰マトリクス
def DampMatrix(k, input):
    omega = 2 * math.pi * input.natural_frequency
    return 2 * input.damp_factor * k / omega

質点系モデル作成

前項で生成したマトリクスを用いて質点系モデルをつくります。

PythonScript

# 質点系モデル
class Model:
    def __init__(self,input):
        arr_m = input.arr_m
        arr_k = input.arr_k
        self.size = len(arr_m)
        arr_k.reverse()
        arr_m.reverse()
        self.m = MassMatrix(arr_m, input)
        self.k = StiffMatrix(arr_k)
        self.c = DampMatrix(self.k, input)

解析

Newmarkβ法を用いて数値解析を行い、入力波形のステップごとに加速度・速度・変位を求めます。

PythonScript

def dynamic_analysis(model, input):
# 初期化
    m = model.m
    k = model.k
    c = model.c
    acc0 = WaveData(input)    
    unit_vector = np.ones(model.size) 
    pre_acc0 = 0.0
    time = 0.0
    dis = np.zeros(model.size)
    vel = np.zeros(model.size)
    acc = np.zeros(model.size)
    ddis = np.zeros(model.size)
    dvel = np.zeros(model.size)
    dacc = np.zeros(model.size)
    dis_his = {}
    vel_his = {}
    acc_his = {}
    for i in range(0, model.size):
        dis_his[i] = []
        vel_his[i] = []
        acc_his[i] = []
    time_his = []
    # Newmarkβ法で解析
    for i in range(0, len(acc0)):
        kbar = k + (1.0/(2.0*input.beta*input.dt)) * c + (1.0/(input.beta*input.dt**2.0)) * m
        dp1 = -1.0 * m.dot(unit_vector) * (acc0[i] - pre_acc0)
        dp2 = m.dot((1.0/(input.beta*input.dt))*vel + (1.0/(2.0*input.beta))*acc)
        dp3 = c.dot((1.0/(2.0*input.beta))*vel + (1.0/(4.0*input.beta)-1.0)*acc*input.dt)
        dp = dp1 + dp2 + dp3
        ddis = np.linalg.inv(kbar).dot(dp)
        dvel = (1.0/(2.0*input.beta*input.dt))*ddis - (1.0/(2.0*input.beta))*vel - ((1.0/(4.0*input.beta)-1.0))*acc*input.dt
        dacc = (1.0/(input.beta*input.dt**2.0))*ddis - (1.0/(input.beta*input.dt))*vel - (1.0/(2.0*input.beta))*acc
        dis += ddis
        vel += dvel
        acc += dacc
        acc_abs = acc + [acc0[i] for n in range(1,model.size)]
        [dis_his[i].append(x) for i, x in enumerate(dis)]
        [vel_his[i].append(x) for i, x in enumerate(vel)]
        [acc_his[i].append(x) for i, x in enumerate(acc_abs)]
        time_his.append(time)
        time += input.dt
        pre_acc0 = acc0[i]
    return time_his, dis_his, vel_his, acc_his

結果出力

前項で作成した__dynamic_analysis__を実行して結果を求めます。

PythonScript

# メイン
input = Input()
model = Model(input)
time_his, dis_his, vel_his, acc_his = dynamic_analysis(model, input)

グラフの出力

matplotlib でグラフを作成します。
グラフを Watch Image ノードで表示するために、System.Drawing.Bitmap に変換します。

PythonScript

# 画像変換
def convertToBitmap(fig):
    rgba_buf = fig.canvas.buffer_rgba()
    (w,h) = fig.canvas.get_width_height()
    npImgArray = np.frombuffer(rgba_buf, dtype=np.uint8).reshape((h,w,4))   
    bitmap_ = None
    # alphaチャンネルを削除
    if npImgArray.ndim == 3 and npImgArray.shape[-1] == 4:
        npImgArray = npImgArray[:, :, :-1]
    img = Img.fromarray(npImgArray, "RGB")
    # 画像をバイナリデータに変換し、メモリ上に保存
    byteIO = io.BytesIO()
    img.save(byteIO, format='BMP')
    byteArr = byteIO.getvalue()
    # convert to Net ByteArray
    netBytes = System.Array[System.Byte](byteArr)
    with MemoryStream(netBytes) as ms:
        bitmap_ = Bitmap(ms)
    return bitmap_
    
# グラフ
fig = plt.figure(figsize=(10,10))

ax1 = fig.add_subplot(3,1,1)
ax1.set_ylabel("acc") # 縦軸のタイトル
ax1.plot(time_his, acc_his[1], label="1F", color="darkturquoise") # 1階のグラフ
ax1.plot(time_his, acc_his[0], label="2F", color='limegreen') # 2階のグラフ
ax1.legend(loc='lower right') # 凡例の位置
ax1.grid(True) # グリッド

ax2 = fig.add_subplot(3,1,2)
ax2.set_ylabel("vel") # 縦軸のタイトル
ax2.plot(time_his, vel_his[1], label="1F", color="darkturquoise") # 1階のグラフ
ax2.plot(time_his, vel_his[0], label="2F", color='limegreen') # 2階のグラフ
ax2.legend(loc='lower right') # 凡例の位置
ax2.grid(True) # グリッド

ax3 = fig.add_subplot(3,1,3)
ax3.set_xlabel("time") # 横軸のタイトル
ax3.set_ylabel("dis") # 縦軸のタイトル
ax3.plot(time_his, dis_his[1], label="1F", color="darkturquoise") # 1階のグラフ
ax3.plot(time_his, dis_his[0], label="2F", color='limegreen') # 2階のグラフ
ax3.legend(loc='lower right') # 凡例の位置
ax3.grid(True) # グリッド

fig.canvas.draw()
OUT = convertToBitmap(fig)

DateTime.Now ノードを使ってアニメーションにするのも面白いと思います。

おわりに

スクリプト単体では使い道がないですが、データ操作やアニメーション出力の箇所はいろいろと応用できそうな気がします。