目次
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 ノードを使ってアニメーションにするのも面白いと思います。
おわりに
スクリプト単体では使い道がないですが、データ操作やアニメーション出力の箇所はいろいろと応用できそうな気がします。