OpenSeesPyを用いた梁要素の線形座屈解析 -その2線形座屈モードの描画-

2024年4月14日 19:00

はじめに

やぶうちです。
気が向いたので、前回の続きを書きます。
前回はOpenSeesPyを用いて梁要素の線形座屈解析プログラムを構築しました。
詳細は下記リンクをご参照下さい。

OpenSeesPyを用いた梁要素の線形座屈解析 -その1線形座屈荷重の算出-｜やぶうちゆうま (note.com)

プログラム概要

今回は前回構築したプログラムをもとに、Matplotlibを用いて線形座屈モードをプロットしましょう。
解析対象物は前回と同様、下記に示す片側に圧縮力を受ける柱とします。
プロットにあたり、前回の記事で"xi"と定義した線形座屈モードに関するベクトルが必要となります。
したがって、前回構築したプログラムの続きから今回のプログラムを構築することとします。

さっそくプログラムを構築していきましょう。
まず、Matplotlibをimportしましょう。
"xi"という関数をMatplotlibでプロットするための形式へ変換するため、今回もNumpyとOpenSeesPyはimportしておきます。
頭文字が"plt.rcParams"となっている部分では、フォントとフォントサイズの指定を行っています。
今回は図に文字を記述しませんが、僕はMatplotlibをimportするとき原則この3行を宣言するようにしています。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import openseespy.opensees as ops
import openseespy.postprocessing.ops_vis as opsv

plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['mathtext.fontset']='cm'
plt.rcParams['font.size'] = 18

プロットした図はPDFまたはPNG、JPGなど好みの形式で出力できますが、ここでは論文や学会発表の際のパワーポイントに図を掲載することを意識してPDFおよびPNG形式で出力することとします。
図を出力するフォルダー名を"figure"の頭文字をとってここでは"fig"と名付けています。
今回は下記のようにPythonで当該フォルダーを作成していますが、手動で作成していただいても結構です。

import os
try:os.stat('fig/')
except:os.mkdir('fig/'

プログラムの構築と解説

プログラム全体は以下のようになっています。

sfac=0.01
for mode in [1,2,3,4,5]:
	ax=Axes3D(plt.figure(dpi=200,figsize=(6,6)))
	ele_tags = ops.getEleTags()
	for ele_tag in ele_tags:
		nd1,nd2 = ops.eleNodes(ele_tag)
		DOFs_nd1 = ops.nodeDOFs(nd1)
		DOFs_nd2 = ops.nodeDOFs(nd2)
		
		ex = np.array([ops.nodeCoord(nd1)[0], ops.nodeCoord(nd2)[0]])
		ey = np.array([ops.nodeCoord(nd1)[1], ops.nodeCoord(nd2)[1]])
		ez = np.array([ops.nodeCoord(nd1)[2], ops.nodeCoord(nd2)[2]])
		
		# if mode No.
		if DOFs_nd1[0] != -1:ed1 = -xi[DOFs_nd1[0],mode]
		else:ed1=0
		if DOFs_nd1[1] != -1:ed2 = -xi[DOFs_nd1[1],mode]
		else:ed2=0
		if DOFs_nd1[2] != -1:ed3 = -xi[DOFs_nd1[2],mode]
		else:ed3=0
		
		if DOFs_nd2[0] != -1:ed4 = -xi[DOFs_nd2[0],mode]
		else:ed4=0
		if DOFs_nd2[1] != -1:ed5 = -xi[DOFs_nd2[1],mode]
		else:ed5=0
		if DOFs_nd2[2] != -1:ed6 = -xi[DOFs_nd2[2],mode]
		else:ed6=0
		
		ed = np.array([ed1,ed2,ed3,ed4,ed5,ed6])
		ax.plot(ex,ey,ez,color='gray',linestyle='--')
		xdi,ydi = opsv.beam_disp_ends(ex,ey,ed,sfac)
		ydi,zdi = opsv.beam_disp_ends(ey,ez,ed,sfac)
		ax.plot(xdi,ydi,zdi,color='black',linestyle='-')

	ax.view_init(elev=20, azim=90)
	plt.axis('off')
	plt.savefig('fig/mode_'+str(mode)+'.png',transparent=True)
	plt.savefig('fig/mode_'+str(mode)+'.pdf',transparent=True)

順を追って上記のプログラムの説明していきますね。
"sfac"で座屈モードの描画倍率を指定しています。
見やすい倍率に適宜調節して下さい。
なお、座屈モードの単位は[mm]や[m]でなく、無次元量です。
For文を使用しており、"mode"という関数を1から1刻みで5まで繰返し処理させています。
modeは座屈モード次数を示しており、今回は1次から5次までプロットすることとします。
"ax=Axes3D(plt.figure(dpi=200,figsize=(6,6)))"
の部分でプロットする図のベースを作っています。
dpiは解像度でfigsizeは図の縦と横の比率であり、見やすい図となるよう適宜調節します。
ここでは黒色の実線で変形後の構造物、灰色の点線で変形前の構造物をそれぞれプロットしています。
plt.savefig~の部分から図の保存を行っていますが、括弧内の
"transparent=True"
は原則宣言することを推奨します。
これは図の背景を透明にする宣言であり、TeXで論文を書く際に図の位置を調整する際は当該宣言を行った図を使用していないと後々困ることになります。

以上で説明を終わります。
また気が向いたら記事を書く所存です。
先日リリースされたRhino8のGHPythonから、OpenSeesPyやNumpyが使えるようになったらしいので、パラメトリックに構造物の構造解析を行う記事を書こうかなと考えています。

読んでいただきありがとうございました。

この記事が気に入ったらサポートをしてみませんか？