ラベル OpenSeesPy の投稿を表示しています。 すべての投稿を表示
ラベル OpenSeesPy の投稿を表示しています。 すべての投稿を表示

2024/11/10

OpenSeesPyによる多質点モデルの時刻歴応答解析④

 では、これまでに設定したモデル、荷重、解析条件に従って、時刻歴応答解析を行い結果を出力します。時刻歴結果の出力にあたっては、各時刻の応答加速度、速度、変位をそれぞれ取得するために以下のコードを記述し、結果出力用のディクショナリを作成します。
  
outputs = {
"time": [],
"rel_accel": [],
"rel_vel": [],
"rel_disp": [],
}

 次に、全解析ステップの解析を解析時間の終了ステップまで指定した時間刻みで実施し、各ステップの解析結果(解析モデルの節点の応答加速度、応答速度、応答変位)をそれぞれ設定したディクショナリに記録します。結果の取得においてはいずれも、OpenSeesPyのOutputcommandを使用しています。

analysis_time = (len(values)) * dt#解析時間の算定
#解析時刻が解析終了時刻になるまで以下をループ
while op.getTime()< analysis_time:
    #現時刻の取得
    curr_time = op.getTime()
    #OpenSeesの構造解析を1ステップ分実行
    op.analyze(1, dt)    #各解析ステップの時刻、節点加速度、節点速度、節点変位の取得(節点4)
    outputs["time"].append(curr_time)
    outputs["rel_disp"].append(op.nodeDisp(4,1)/ratio)
    outputs["rel_vel"].append(op.nodeVel(4,1)/ratio)
    outputs["rel_accel"].append(op.nodeAccel(4, 1) / ratio)

 得られた各節点の応答結果をグラフに書いてみると以下の通りとなります。


 特に節点4の応答時刻歴については、時刻歴の振動外力を入力していることから、入力加速度波形の継続時間(50s)の間、各時刻で変動する波形が得られていることがわかります。また、固定点である基部節点1の応答加速度を見ますと、ゼロになっていることがわかります。
 一方、地震応答解析においては、固定点であったとしても解析上は全節点に地動加速度が入力されているはずです。上記より、Openseesの地震応答解析においては、節点の応答加速度は相対座標系のもとで出力されていることが確認できるかと思います。

 今回は構造減衰や部材の非線形性などは考慮しない極めて単純なモデルで応答解析を実施しました。

 なお、設定した解析モデルの条件を改めて確認する場合、例えば以下のコードを書くことで、OpenSeesPyからのエコーデータを得ることができます。以下のケースではJSON形式でテキストファイルとして出力しています。
 
op.printModel('-JSON', '-file', "model.json")

本例でのエコーデータは以下の通りとなります。エコーデータの単位系はOpenSeesPyの単位系(ヤードポンド法)となります。

{
"StructuralAnalysisModel": {
	"BIM": "unknown",
	"description": "",
	"engineer": "",
	"units": {
		"force": "",
		"length": "",
		"time": "",
		"temperature": ""
	},
	"properties": {
		"uniaxialMaterials": [

		],
		"ndMaterials": [

		],
		"sections": [

		],
		"crdTransformations": [
			{"name": "1", "type": "LinearCrdTransf2d"}
		]
	},
	"geometry": {
		"nodes": [
			{"name": 1, "ndf": 3, "crd": [0, 0]},
			{"name": 2, "ndf": 3, "crd": [0, 157.48], "mass": [0.116455, 0.116455, 0]},
			{"name": 3, "ndf": 3, "crd": [0, 314.961], "mass": [0.116455, 0.116455, 0]},
			{"name": 4, "ndf": 3, "crd": [0, 472.441], "mass": [0.116455, 0.116455, 0]}
		],
		"elements": [
			{"name": 1, "type": "ElasticTimoshenkoBeam2d", "nodes": [1, 2], "E": 29732.7, "G": 11435.7, "A": 474.774, "Avy": 237.387, "Iz": 350013, "massperlength": 0, "crdTransformation": "1"},
			{"name": 2, "type": "ElasticTimoshenkoBeam2d", "nodes": [2, 3], "E": 29732.7, "G": 11435.7, "A": 474.774, "Avy": 237.387, "Iz": 350013, "massperlength": 0, "crdTransformation": "1"},
			{"name": 3, "type": "ElasticTimoshenkoBeam2d", "nodes": [3, 4], "E": 29732.7, "G": 11435.7, "A": 474.774, "Avy": 237.387, "Iz": 350013, "massperlength": 0, "crdTransformation": "1"}
		]
	}
}
}


2024/11/04

OpenSeesPyによる多質点モデルの時刻歴応答解析③

 地震応答解析においては地震荷重や、直接積分手法などの解析条件を定義する必要があり、 ここではそれらをOpenseespyで定義する方法を説明します。

1.時間刻みの設定:

dt = 0.01
ここでは、解析の時間刻み(Δt)を0.01秒に設定しています。この値はシミュレーションの時間解像度を示し、入力する地震波のデータの刻み幅を入力します。

2.外乱となる時刻歴データの読み込み:

motion = pd.read_csv('white_noise.csv')
motion = list(motion['time'])
ratio = 0.393700787 #cm/s2→in/s2
values = [round(float(i)*ratio,4) for i in motion]
op.timeSeries('Path', 1, '-dt', dt, '-values', *values)

地震応答解析においては、外乱データは時刻歴のデータとなります。そこで、その時刻歴のデータを外部csvファイルなどに用意し、そのデータの読み込みをここで定義しています。地震波の加速度入力データは通常gal(cm/s2)で表されることが多いですが、ここでは、単位系はOpenseespyで用いられている単位系(in/s2)に適宜換算する必要があります。ここでは例題としてホワイトノイズ加振を実施するため、適当な乱数に基づいて下図のようなホワイトノイズの時刻歴波形をcsvデータとして用意し入力します。入力データをリスト形式のデータに変換し、単位換算用の係数を乗じた上で、timeSeriesコマンドにて、Path形式でOpenSeesPyに指定した時間刻みの時刻歴荷重として読み込ませています。


3.荷重パターンの設定:
op.pattern('UniformExcitation', 1, 1, '-accel', 1)

一様励振(Uniform Excitation)パターンを荷重パターンとして設定し、入力した構造モデル全体を対象に加振します。ここでは全体座標系X方向に、2で定義した時刻歴荷重データを加速度として入力し加振する設定となっています。

4.時刻歴応答解析の準備:

op.wipeAnalysis()
op.algorithm('RaphsonNewton')
op.system('BandGen')
op.constraints('Plain')
op.integrator('Newmark', 0.5, 0.25,'A')
op.analysis('Transient')

これらの行では、時刻歴応答解析のための設定を行っています。

  • wipeAnalysis() では以前の解析設定をクリアしています。(ここでは以前の解析設定はありませんが、念のため定義します。)
  • algorithm('RaphsonNewton') では、非線形問題を解く際のアルゴリズムを選択定義します。ここではニュートンラプソン法を用いた解法を選択しています。
  • system('BandGen') は、連立方程式を解くための線形ソルバーを定義します。ここでは、とりあえずBandGeneralを定義しています。
  • constraints('Plain') で、制約条件を定義します。
  • integrator('Newmark', 0.5, 0.25,'A') では解析における積分手法を定義しています。ここでは、Newmark、γ=0.5、β=0.25をそれぞれ指定することで、平均加速度法を用いることを定義しています。この手法は通常の動的解析で広く用いられます。
  • analysis('Transient') ではどのような解析を行うかを定義します。Transient Analysisを指定し、時刻歴荷重の入力を想定した時間変化を考慮した解析を行う準備をしています。時刻歴解析を行うことから、結果も時刻歴で出力されることになり、解析結果を記録するコードを書く必要があり、後ほど説明します。
全体として、このコードは地震応答解析の初期設定を行い、地震などの外力に対する構造物の応答を計算するための準備をしています。ここまでで構造モデルと時刻歴解析の荷重条件、解析条件が定義できましたので次は解析を実行し結果を見ていきます。

2022/09/25

OpenSeesPyによる多質点モデルの時刻歴応答解析②

時刻歴応答解析の実施に当たり、まず解析モデルを設定します。解析モデルは図に示すような高さ方向に集中質量を有する質点が並ぶ2次元の多質点モデルとします。各質点の間は、梁要素があるものとし、各要素の断面特性をそれぞれ設定することで、層間の構造特性を定義したモデルとなっています。

対象解析モデル

まず、モデル全体の諸元として、対象は2次元平面にあるモデルで、各節点において、3つの方向の自由度(鉛直、水平、回転方向)を有するモデルであることから、以下のように定義します。

# モデルの初期化
wipe()
# モデルの定義
model('basic', '-ndm', 2, '-ndf', 3)

次に各節点の座標を定義します。2次元のモデルであることから、節点座標は、2つの座標情報にて定義されます。例題のようなモデルを想定する場合、以下のように記述します。
Openseespyでの処理上、単位系はヤードポンド法の長さの単位であるインチに変換し、各節点の座標情報を記述します。
# 節点座標の定義
m_inch=39.3700787#単位系m→inchへの変換
op.node(1,0,0*m_inch)#節点1の座標設定(節点番号,x座標,y座標)
op.node(2,0,4*m_inch)#節点2の座標設定
op.node(3,0,8*m_inch)#節点3の座標設定
op.node(4,0,12*m_inch)#節点4の座標設定
解析モデルにおいては、境界条件を設定する必要があります。今回の場合、解析モデルの基礎部分に対して境界条件を設定するため、基部の節点1において、境界条件(fix command)を設定します。コマンドの引数において、0:非拘束、1:拘束となります。今回の場合、3自由度ともに固定とするので、以下のように設定します。
# 境界条件の定義
op.fix(1,1,1,1)
#節点1の境界条件(対象節点番号,拘束条件(鉛直、水平、回転自由度固定(1,1,1)))
各節点の質量を定義します。インチ系の加速度で規準化された値で入力する必要があるので単位をインチで換算した重力加速度386in/s2で値を割る必要があることに注意が必要です。日本人の感覚だとこの単位系の扱いが慣れず、とても分かりづらく感じます。とりあえず各節点に対し回転慣性質量は考えずに、鉛直方向,水平方向それぞれに質量を以下のように入力します。
# 節点質量の定義
kN_pond=0.22480894387096#重量単位系kN→kipsへの変換
g = 386.0885827#インチ系における重力加速度(in/s2)
op.mass(2,200*kN_pond/g,200*kN_pond/g,0)
#節点2の節点質量設定(対象節点番号,鉛直方向の質量,水平方向の質量,回転方向の質量)
op.mass(3,200*kN_pond/g,200*kN_pond/g,0)#節点3の節点質量設定
op.mass(4,200*kN_pond/g,200*kN_pond/g,0)#節点4の節点質量設定
最後に各質点を結ぶ梁要素の諸元を設定します。
梁要素は、element commandの中で装備されている種々の特性を有する要素を設定できます。ここでは、各要素は弾性梁とし、曲げ、せん断変形を共に考慮できる「Elastic Timoshenko Beam Column Element」を用いることにします。梁要素の剛性としては、とりあえず鋼材を想定してヤング係数、せん断弾性係数を設定します。断面形状はとりあえず直径2m、板厚50mmの円筒(A=0.145686451m2、I=0.306305284m4)を想定してみます。各要素に対して、それぞれ断面積、断面二次モーメント、せん断断面積を定義します。
# 梁要素の定義
E_mod=205000000*kN_pond/((m_ft)**2);G_mod=78846153.8461538*kN_pond/((m_ft)**2)#鋼材を想定した材料定数の設定
op.geomTransf("Linear",1)#要素座標系から全体座標系への変換手法の定義(線形変換)
op.element('ElasticTimoshenkoBeam', 1, 1, 2, E_mod, G_mod, 0.306305284*(m_ft)**2, 0.145686451*(m_ft)**4,0.5*0.306305284*(m_ft)**2,1)#梁要素1の諸元設定
op.element('ElasticTimoshenkoBeam', 2, 2, 3, E_mod, G_mod, 0.306305284*(m_ft)**2, 0.145686451*(m_ft)**4,0.5*0.306305284*(m_ft)**2,1)#梁要素2の諸元設定
op.element('ElasticTimoshenkoBeam', 3, 3, 4, E_mod, G_mod, 0.306305284*(m_ft)**2, 0.145686451*(m_ft)**4,0.5*0.306305284*(m_ft)**2,1)#梁要素1の諸元設定
以上より、対象とする多質点モデルの入力が終わりました。次にこの解析モデルに対して、地震波を入力した時刻歴応答解析を実施してみたいと思います。

2022/09/24

OpenSeesPyによる多質点モデルの時刻歴応答解析①

建物構造物などは多質点モデルなどに簡略化したうえで、地震波を入力した時刻歴応答解析を実施することがよくあります。Openseespyでは構築した種々の構造モデルの時刻歴応答解析が可能ですが、今回から数回に分けて、Openseespyにおける多質点モデルの時刻歴応答解析の実施例を紹介します。なお、対象モデルは多質点の曲げせん断モデルとします。地震応答解析の手順は以下の通りとなります。

①多質点モデルの構築

 これまで紹介してきた応力解析事例と同様に、解析対象となる構造物のモデルを定義します。今回は説明のため、単純な2次元の多質点モデルを想定し、多質点曲げせん断モデルの節点座標、節点重量、梁要素の諸元、境界条件をOpenseespyのコマンドを用いて定義していきます。

②入力荷重の定義

 入力荷重を定義します。時刻歴応答解析では、地震、風荷重などの各時刻で変化する荷重を入力する荷重データとして想定することになります。それらをOpenseespyのコマンド(timeSeries command)を用いて定義します。

③地震応答解析の条件設定

 解析条件として、どのようなソルバーを用いるか、また、構造減衰の設定条件などを定義します。Openseespyには、様々な解析条件の仕様が装備されているようです。

④応答解析結果の出力

 各時刻の解析結果のデータをpythonのリストなどに記録することで、出力します。

以降、何回かに分けて、解析コードとともに、記述方法を説明していく予定です。

2021/06/20

Rhino+GrasshopperでのCPythonの実行環境(GH_CPython)

RhinocerosにおけるビジュアルコーディングツールであるGrasshopperには、python言語でのコーディングが可能なGh pythonがあらかじめ付属しています。しかし、このGh pythonは通常のPython(CPython)環境ではなく、Rhinocerosが動いているNET Framework上で動作するironPythonであるため、Rhinoceros上での描画アドインは動かせるものの、一般的にPython言語で用いられている多くの拡張モジュールをGh python上で通常実行することができません。Python環境において複雑な数値計算や最適化、応力解析などを実行できる無数の拡張モジュール(Numpy、Deep、Openseespyなど)が、一般的には利用できないのです。

一方、こうした数値計算は建築設計や構造解析において多く活用されており、これらをRhino+Grasshopperのプログラム上でもなんとか利用できる環境が構築できると便利です。

上記を可能とするために、Grasshopper上で一般的なpython言語であるCPythonを実行できるライブラリGH_CPythonの活用が考えられます。本ページではそのインストール方法を紹介してみます。なお、GH_CPythonを用いたRhino+Grasshopper上でのCPythonの実行においては、あらかじめPythonライブラリのインストールが必要となります。

なお、Pythonのインストールについては、本ブログの別ページ Kの勉強部屋: OpenSeesPyの導入① (kozoarcheng.blogspot.com) をご覧ください。


<GH_CPythonの実行環境の構築>

1.Githubまたはfood4Rhinoにて、GH_CPythonをダウンロード

GH_CPython | Food4Rhino

2.ダウンロードしたファイルのうち、GH_CPython.gha と FastColoredTextBox.dllを自身のRhino+Grasshopperインストールファイル中のLibrariesフォルダ(もしくはComponentsフォルダ)に移動。(ダウンロードファイルはプロパティにてセキュリティを解除する。)

3.Rhinocerosを管理者アカウントで実行

4.GH_CPythonをGrasshopper上で開き、「Python」→「Choose Interpreter」をクリックし、以下のように自身のPCにインストール済みのpython.exeのディレクトリを指定



基本的には、上記の操作でGrasshopper上にて、CPythonの環境構築ができると思います。

CPythonにおいては、本サイトで使用法を紹介している応力解析アドインOpenseespyなども実行できることから、上記の環境構築により、Rhino+Grasshopperで構築したプログラムフローにおける3Dcadやビジュアルコーディングと構造物の応力解析、時刻歴応答解析などの数値解析を連動させた処理が可能となり、使い方によっては解析上非常に有用なツールとなりえるかと思います。

上記で説明した内容については、作者(MahmoudAbdelRahman氏)のGithubである下記にも記載がありますので、適宜ご参照頂ければと思います。

Home · MahmoudAbdelRahman/GH_CPython Wiki · GitHub

2020/05/06

トラス要素から梁要素への拡張(Openseesによる3次元応力解析)

前回構築したOpenseesによる3次元フレーム要素の弾性応力解析プログラムについて、補足的に説明します。プログラムのソースはあまり3次元トラス要素のプログラム構成と変わりません。そこで、前回までに構築した3次元トラス要素プログラムの一部を修正したポイントを以下にかいていきます。扱う部材がトラス要素→梁要素となったことで以下のように部材の曲げモーメント、せん断力を考慮する必要がでてきます。修正したポイントはその点に関連するところがほとんどです。



・節点固定度、荷重項を3自由度→6自由度に

トラス要素では部材端がピンのみでしたが、梁要素では部材端が剛接の場合も考えられます。そこで、節点については回転方向の固定、作用力も考慮できるようにする必要がありました。部材の回転を考えない場合、部材固定度、節点荷重も3自由度にて表現されますが、回転の概念が入ってくることで、6自由度となります。そこで、プログラムのインプットにおいて6自由度を考慮できるようにリストの構成等を修正しました。

・要素をトラス要素→梁要素に

トラス要素は軸力のみを伝達することを考えるため、断面積のみを設定する仕様ですが、フレーム要素にした場合曲げも伝達することになるため、部材各方向の断面二次モーメント、ねじれモーメントを設定する必要があります。Openseesでは弾性梁要素として、elasticBeamColumnコマンドにて以下のように部材特性を設定します。

element('elasticBeamColumn', eleTag, *eleNodes, Area, E_mod, G_mod, Jxx, Iy, Iz, transfTag)

eleTag:部材番号
eleNodes:材端の座標
Area:断面積
E_mod:弾性係数
G_mod:せん断弾性係数
Jxx:ねじり断面二次モーメント
Iy:部材座標系y軸方向に対する断面二次モーメント
Iz:部材座標系z軸方向に対する断面二次モーメント
transfTag:部材要素の座標変換に関する設定(後述)

・部材要素の座標変換に関する設定の追加

全体座標系→部材座標系の変換に以下のコマンドを新たに入力する必要がありました。(通常トラス要素の応力解析でもこのような座標変換は行われていることから、あくまでもOpenseesならではの仕様なのかもしれません。)

transfTag = 1
vecxz =0,0,1
geomTransf('Linear', transfTag, vecxz)
’Linear’:線形での座標変換 transfTag:座標変換を定義するタグ(それぞれの部材要素に設定) vecxz:ユーザーが指定する基本ベクトルであり、部材座標系のx軸と並行であってはならない。 (ここでは解析の対象モデルに鉛直部材がないことを想定して、全体座標系のZ軸方向の単位ベクトルを設定)

上記の設定を各部材要素に対してタグを設定することにより当てはめていきます。

・解析結果の出力

トラス要素の応力解析では部材に軸力のみが生じましたが、梁要素では曲げ及びせん断力が各部材の材端に生じるため、各要素の作用力を出力するbasicforceコマンドを適用すると、1部材につき、6個の数値がリストとして出力されます。そこで、それらのデータを整理し、軸力、曲げモーメント、せん断力のそれぞれについて、各部材のアウトプットをグラデーションによる色分布で図化出力できるように修正しました。前回投稿で記載したものです。

plotlyによる図化出力においては、各部材の線グラフのデータにhovertextとして作用力の数値データを埋め込み、マウスをhoverさせることにより、各部材の作用力が直感的に確認できるようにしました。
ただ、このplotly、表示機能を盛り込みすぎると、動作速度が大変遅くなってしまうのが難点ではあります。(本来はグラフ描画ツールなので。。)みたい結果をその都度絞りながら、結果を表示させていくことが現実的な使い方かもしれません。


以上を考慮することにより、Openseesを使いながら、3次元の任意形状の構造モデルを対象とした梁要素の応力解析プログラムが構成できました。トラス要素の構造解析プログラムとあまり構成が変わらないので、ソースコードの記載は割愛します。上記のポイントを押さえながら比較的簡単に拡張できると思いますので、試してみてください。

2020/05/05

Grasshopperによる3次元ジオメトリのOpenseesによる応力解析

今回は、Openseesにより構成した3次元の応力解析プログラムの実務的な利用について考えてみます。特に複雑な形状モデルを対象とする場合、3次元の形状データから、目視や手入力による作業で応力解析プログラムのインプットを作成することを考えると、膨大で複雑な数値データの入力が必要となることから、大変面倒な作業が必要となってしまいます。そこで、近年建築分野で普及が進んでいる3DcadであるRhinoceros+Grasshopperにより作成した形状データをもとにプログラムによるOpenseesの解析入力ファイルの作成にトライしてみます。また、Rhinoceros+Grasshopperによる入力データの作成および解析結果のKaramba3Dによる解析結果の簡単な検証を試みます。

・Rhinoceros+Grasshopperとは
RhinocerosはRobert McNeel & Associates社が開発した3次元CADであり、NURBS局面などの複雑なジオメトリを柔軟に扱えることから、近年、建築分野でも多く活用されています。またGrasshopperはRhinocerosをスクリプトベースで自由に操ることができるツールで、視覚的なアイコンをつなぐことで容易にプログラムが構成できるため、コンピュテーショナルデザインなどに活用されています。さらに近年は世界中のエンジニア・プログラマーにより開発された最適化や応力解析などの多くの外部モデュールが公開されており、その多くがフリーで利用可能となっています。そうしたツールは以下のサイトfood4Rhinoで利用可能となっています。

https://www.food4rhino.com/

・Karamba3Dとは
Grasshopperには、多くの外部機能がPythonのモジュールのように提供されており、Karamba3DはRhinoceros上のジオメトリを用いた応力解析を容易におこなうことができる応力解析モジュールです。近年は有料化されていますが、(unlimitedlicenceで1150ユーロ)多くの構造エンジニアにも使用され始めています。また、最適化ツールなどの機能拡張が日々行われています。

https://www.karamba3d.com/

まず、応力解析の対象とするトラス屋根モデルを以下のように懸垂円筒状のトラスモデルとして作成してみました。作成に当たってはトラス要素のジオメトリを簡便に作成できる外部モデュール(Lunchbox)を使用しています。Lunchboxについても先述したfood4Rhinoにてフリーで弾雨ロードすることが可能です。


Grasshopper上では、上記のようなアイコンをつなげるだけの簡単なコーディングで以下のような複雑なトラスの構造モデルが作成できます。立体的な幾何学を扱うにはとっても便利なツールです。



このような複雑な形状の構造解析モデルの座標や要素の情報を一つ一つinputファイルに記述するのは例えExcelファイルによる記述であってもとても大変です。

しかし、幸いなことに、Grasshopper上ではPythonで書いたプログラムをコンポーネントとして扱うことができるため、作成したジオメトリデータ(line、座標データなど)とPythonを連携させて使用することできます。そこで、Grasshopper上で作成したジオメトリデータをOpenseesの構造解析用のインプットデータとしての節点、要素の数値データに置換するプログラムをPythonで作成し、そのプログラムをつかってOpenseesによる3次元解析プログラムのインプットデータを作成しました。プログラムといっても、Rhinoceros上のジオメトリの節点座標と、各梁要素の材端座標を拾って、input用のデータ書式に再構成するだけのシンプルなものとなります。

例題として、作成した解析モデルの頂部の中心節点に鉛直荷重が加わった場合の応力解析をOpenseesにより行います。
前回可視化に利用したPlotlyにより、Openseesによる応力解析結果を視覚化するとこんな感じになります。


軸力図

軸力図については圧縮材が赤色、引張材が青色となるように設定しました。


最大曲げモーメントの分布図


そして変位図


Karambaによる変位図

上記のOpenseesによる解析結果とKarambaによる解析結果とを比較すると、変形性状は概ね対応していることがわかります。上記では、Grasshopper→Openseesによる応力解析により複雑な形状の構造物の応力解析ができました。Pythonによるコーディングにおいて幾何学形状を扱うことは労力がかかります。そうしたことから、このように幾何学の扱いに長けており、デザイナーによるコーディングも可能なRhinoceros+Grasshopperをデザインなどで使用し、それらの分析にOpenseesを活用するといったスキームは、その間をつなぐ環境がしっかりと構築されれば、建築設計実務でも十分活用できるのではと考えます。

2020/04/29

Plotlyによる応力解析結果の出力

前回まではmatplotlibを用いた2次元画像による応力解析結果の可視化を試してみました。
しかし、3次元の構造解析モデルの応力解析結果を確認する場合、汎用解析ツールのUIのように
同じく3次元空間でグリグリとオブジェクトを回しながら結果確認ができると実務的に便利です。
そこで今回は、Pythonにも対応しているインタラクティブなデータ可視化ライブラリPlotlyを用いてOpenseesPyによる3次元トラスの応力解析結果の出力を行います。

Plotlyはオープンソースのデータ可視化ライブラリです。
Pythonだけでなく、JAVAscriptやMATLAB等にも対応しています。
同じデータ可視化ライブラリであるmatplotlibは主にグラフデータの画像出力に対応していましたが、
Plotlyでは出力したグラフデータをブラウザ上でマウスなどを動かしながらインタラクティブな状態で確認できます。
Pythonでの使用にはOpenseesPyの場合と同じく、ライブラリのインストールおよびプログラムへのインポートが必要となります。
インストール法や使用方法の詳細につきましては、以下サイトをご覧ください。

https://plotly.com/python/

今回もこれまでに構築した3次元トラス解析プログラムをベースに検討を進めます。
Plotlyを使ってデータをグラフ化するプロセスを以下に簡単に示します。

・ライブラリのインポート

import plotly.offline as po
import plotly.graph_objs as go
上記によりPlotlyをプログラムにインポートします。

・初期化

po.init_notebook_mode()
上記によりPlotlyを初期化します。

・グラフの定義

graph=go.Scatter3d(x=node_x, y=node_y, z=node_z, name='beam%s'%(int(ele_no[i])), mode='lines+text', line=dict(width=int(area[i]),color=color_tag),
hovertext='AxialForce %s'%(round(temp[0],2)), showlegend=True)
上記スクリプトの形式で描画するグラフを定義します。
なお、上記スクリプトはモデル図兼応力図を描くためのスクリプトになります。

モデルの各梁要素を線グラフとして描画します。線の太さは梁要素の断面積のデータから、線の色は応力解析により得られた梁要素の軸力から決定するように定義しています。なお、軸力についてはグラフ上でカーソルを合わせると確認できる「hovertext」として定義しています。

・レイアウトの定義とグラフの描画

data.append(graph)
layout = go.Layout(xaxis=dict(showgrid=False))
fig2 = dict(data = data, layout=layout)
po.plot(fig2, filename="model")
定義したグラフデータを保存したdataリストとグラフのレイアウト等を定義したlayoutリストをfig2に呼び出して、po.plotにて描画しています。


上記のPlotlyによるデータ図化機能を応用して、結果としてグリグリ動かせるモデル軸力図・変形図が出力できました。
ブラウザ上の画面で動かせる形式になっており、
マウスを梁上に移動させると軸力も表示させることができます。


変位図を出力したり、画面上で選択することで、モデル図応力図と変位図の2つのグラフを同時に出現させたりすることもできます。
マウスで図を回転させたり、拡大したりすることが非常に直感的に行えます。



以上、汎用の構造解析ツールのような操作感覚で、解析結果を確認することができました。
今後さらに複雑な立体トラスモデルの応力解析結果を確認する場合には、非常に役立ちそうなツールだと思います。
最後に今回使用したソースコードを以下に示します。
#plotlyによる描画
po.init_notebook_mode()
trace1 = go.Scatter3d(x=cor_x, y=cor_y, z=cor_z, mode='markers+text', marker_color='blue', marker_size=2, name="Node", text=no, showlegend=True)
data = [trace1]

#モデル図・軸力図の描画
for i in range(len(ele_no)):
    temp2=[]
    temp3=[]
    node_x=[]
    node_y=[]
    node_z=[]
    temp2=nodeCoord(int(n1[i]))
    temp3=nodeCoord(int(n2[i]))
    node_x.append(temp2[0])
    node_x.append(temp3[0])
    node_y.append(temp2[1])
    node_y.append(temp3[1])
    node_z.append(temp2[2])
    node_z.append(temp3[2])
    
    temp=basicForce(int(ele_no[i]))#各部材の応力値をRGB値に変換
    if temp[0] >=0:
        if max_tens == 0:
            color_tag='rgb(0,0,'+str(255*temp[0])+')'
        else:
            color_tag='rgb(0,0,'+str(255*temp[0]/max_tens)+')'     
    else:
        if max_comp == 0:
            color_tag='rgb('+str(-255*temp[0])+',0,0)'
        else:  
            color_tag='rgb('+str(-255*temp[0]/max_comp)+',0,0)'
            
    graph=go.Scatter3d(x=node_x, y=node_y, z=node_z, 
                               name='beam%s'%(int(ele_no[i])), mode='lines+text', line=dict(width=int(area[i]),color=color_tag),
                               hovertext='AxialForce %s'%(round(temp[0],2)), showlegend=True)

    data.append(graph)

#変位図の描画
node_x=[]
node_y=[]
node_z=[]
temp2=[]
temp3=[]
for i in range(len(ele_no)):
    temp2=[x + y * mag_factor for (x, y) in zip(nodeCoord(int(n1[i])),nodeDisp(int(n1[i])))]
    temp3=[x + y * mag_factor for (x, y) in zip(nodeCoord(int(n2[i])),nodeDisp(int(n2[i])))]
    node_x.append(temp2[0])
    node_x.append(temp3[0])
    node_y.append(temp2[1])
    node_y.append(temp3[1])
    node_z.append(temp2[2])
    node_z.append(temp3[2])
graph=go.Scatter3d(x=node_x, y=node_y, z=node_z, 
                               name="Disp", mode='lines', line_color='blue',line=dict(width=2, dash='dot'),
                               showlegend=True)
data.append(graph)

layout = go.Layout(xaxis=dict(showgrid=False))
fig2 = dict(data = data, layout=layout)
po.plot(fig2, filename="model")

2020/04/12

Matplotlibによる応力・変形図の図化出力(3D対応)

前回紹介したMatplotlibによる図化出力ですが、以下の手順で3次元トラスの応力解析結果の図化出力にも拡張できました。蛇足ですが説明します。

図化出力のプログラムの頭に以下のコードを記載します。
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

また、グラフのスクリプトを記述する際には、以下を記載します。
ax1 = fig.add_subplot(2,1,1, projection='3d')

以上を考慮するとMatplotlibによるグラフにおいてz軸の概念が導入され、2次元のグラフ描画と同様のルールでグラフの表示カスタマイズ(軸最大値の設定、表示する角度の設定など)も可能となります。Matplotlibによる3次元グラフ描画の詳細については以下のサイトに記載があります。


上記を踏まえて組んだプログラムにより、以下のように3次元のグラフが描画できました。意外といろいろなことができるんですね。Matplotlib。


参考として、今回の図化出力に用いた3次元トラスの応力解析結果の描画部分のソースコードを以下に示します。

#変位図、応力図を描く    
fig = plt.figure(figsize=(6, 9))
max_x=5
max_y=5#グラフの表示範囲を設定
max_z=5

ax1 = fig.add_subplot(2,1,1, projection='3d')
ax2 = fig.add_subplot(2,1,2, projection='3d')
for e in getEleTags():
    color_tag=0
    nodes = eleNodes(e)
    Nnodes = len(nodes)

    xyz = np.zeros((Nnodes,3), dtype=np.double)
    uu   = np.zeros((Nnodes,3), dtype=np.double)
    for i in range(Nnodes):
        xyz[i,:] = nodeCoord(nodes[i])#return coordinate
        uu[i,:] = nodeDisp(nodes[i])#return current displacement

    x = np.zeros(12,dtype=np.double)
    y = np.zeros(12,dtype=np.double)
    z = np.zeros(12,dtype=np.double)    
    u = np.zeros(12,dtype=np.double)
    v = np.zeros(12,dtype=np.double)
    t = np.zeros(12,dtype=np.double) 

    conec = [0, 1]

    x = xyz[conec,0]
    y = xyz[conec,1]
    z = xyz[conec,2]
    x[2::3] = np.nan
    y[2::3] = np.nan
    z[2::3] = np.nan

    mag_factor = 20.#変形図の拡大率
    u = xyz[conec,0] + mag_factor*uu[conec,0]
    v = xyz[conec,1] + mag_factor*uu[conec,1]
    t = xyz[conec,2] + mag_factor*uu[conec,2]
    u[2::3] = np.nan
    v[2::3] = np.nan
    t[2::3] = np.nan
    
    temp=basicForce(e)#各部材の応力をプロット
    if temp[0] >=0:
        if max_tens == 0:
            color_tag=temp[0]
        else:
            color_tag=temp[0]/max_tens*1       
        ax1.text(np.mean(x),np.mean(y),np.mean(z),round(temp[0],2))
        ax1.plot(x,y,z,color=[0,float(1-color_tag),1],lw=area[e-1])
    else:
        if max_comp == 0:
            color_tag=temp[0]
        else:        
            color_tag=-temp[0]/max_comp*1
        ax1.text(np.mean(x),np.mean(y),np.mean(z),round(temp[0],2))
        ax1.plot(x,y,z,color=[1,0,float(1-color_tag)],lw=area[e-1])

    ax2.plot(x,y,z,marker="o",color="0.5")       
    ax2.plot(u,v,t,color="0.",linestyle="dashed")#変形図をプロット

#各グラフのフォーマットの設定
ax1.set_title("Stress Diagram")
ax1.view_init(elev=40, azim=180)
ax1.set_xlim(-1,max_x)
ax1.set_ylim(-1,max_y)
ax1.set_zlim(-1,max_z)
ax1.grid(False)
ax1.set_xticks([]) 
ax1.set_yticks([]) 
ax1.set_zticks([]) 
ax2.set_title("Displacement Diagram")
ax2.view_init(elev=40, azim=180)
ax2.set_xlim(-1,max_x)
ax2.set_ylim(-1,max_y)
ax2.set_zlim(-1,max_z)
ax2.grid(False)
ax2.set_xticks([]) 
ax2.set_yticks([]) 
ax2.set_zticks([]) 
plt.show()

2020/04/11

Matplotlibによる応力・変形図の図化出力

今回はOpenSeesPyによる解析結果の図化出力して遊んでみます。
図化出力の対象は前回までに構築した2次元のトラス解析プログラムの応力解析による結果とします。OpenSeesPyにより出力される節点変位、部材軸力などの数値データを図上に整理します。

図化出力にはpythonのグラフ描画ライブラリMatplotlibを用いることにします。Matplotlibはグラフのみならず様々な図の出力に応用できるようなのです。
私はまだそんなに玄人ではありませんので、少しずつ勉強しながら使っていくことにします。

今回はとりあえず応力解析結果として、部材の軸力図と変形図を出力することにしました。

軸力図は以下のような定義とします。
軸力図における部材の軸力分布はグラフの線色で表現します。
圧縮軸力は赤線、引張軸力は青線で表現、また各部材の軸力の度合いは色のグラデーションで表現します。
さらに、設定した構造部材の断面積を線の太さで表現することにしました。

これらを踏まえて、圧縮軸力となる部材と引張軸力となる部材でそれぞれ場合分けし、軸力図の出力部分は以下のようなコードとしました。

if temp[0] >=0: if max_tens == 0: color_tag=temp[0] else: color_tag=temp[0]/max_tens*1 ax1.text(np.mean(x),np.mean(y),round(temp[0],2)) ax1.plot(x,y,color=[0,float(1-color_tag),1],lw=area[e-1]) else: if max_comp == 0: color_tag=temp[0] else: color_tag=-temp[0]/max_comp*1 ax1.text(np.mean(x),np.mean(y),round(temp[0],2)) ax1.plot(x,y,color=[1,0,float(1-color_tag)],lw=area[e-1])

変形図については、各節点座標とOpenSeesPyのアウトプットコマンドnodeDispで得られた各節点変位の和をグラフ上にプロットします。変位図のプロットの際には、設定した拡大率(mag_factor)により実際の変位を拡大表示します。

以上、組んだプログラムを用いて以下のような応力図が出力できました。


以下のような2次元の任意形状のトラスの応力解析結果の出力に対応します。


このように、色々な形状を設定しながら遊んでみるとトラス構造物の応力の流れが直感的に把握できます。しかし、煩雑なモデルになればなるほどExcelでちまちまと節点や部材の情報を入力するのが面倒になってきます。

そこで今後、CADなどのプログラムと連携させ、より複雑な形状を簡易にインプットできるようなシステムの構築にも挑戦できたらと思います。

今回紹介した図化出力プログラムのソースコードは以下の通りです。
正直、Matplotlibはまだまだ全然使いこなせている感じがしません。
それでも。。
本来はグラフを作成するツールで、枠線、軸線を消したり、2段にまとめたレイアウトとすることによりシンプルな図化出力ができました。
今後も引き続き勉強していきます。

#変位図、応力図を描く    
fig = plt.figure()
max_x=9
max_y=5#グラフの表示範囲を設定

#グラフレイアウトの設定 
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)

#各要素ごとに出力
for e in getEleTags():
    color_tag=0
    nodes = eleNodes(e)
    Nnodes = len(nodes)

    xyz = np.zeros((Nnodes,2), dtype=np.double)
    uu   = np.zeros((Nnodes,2), dtype=np.double)
    for i in range(Nnodes):
        xyz[i,:] = nodeCoord(nodes[i])#return coordinate
        uu[i,:] = nodeDisp(nodes[i])#return current displacement

    x = np.zeros(12,dtype=np.double)
    y = np.zeros(12,dtype=np.double)
    u = np.zeros(12,dtype=np.double)
    v = np.zeros(12,dtype=np.double)

    conec = [0, 1]

    x = xyz[conec,0]
    y = xyz[conec,1]
    x[2::3] = np.nan
    y[2::3] = np.nan

    mag_factor = 10.#変形図の拡大率
    u = xyz[conec,0] + mag_factor*uu[conec,0]
    v = xyz[conec,1] + mag_factor*uu[conec,1]
    u[2::3] = np.nan
    v[2::3] = np.nan
    
    temp=basicForce(e)#各部材の応力をプロット
    if temp[0] >=0:
        if max_tens == 0:
            color_tag=temp[0]
        else:
            color_tag=temp[0]/max_tens*1       
        ax1.text(np.mean(x),np.mean(y),round(temp[0],2))
        ax1.plot(x,y,color=[0,float(1-color_tag),1],lw=area[e-1])
    else:
        if max_comp == 0:
            color_tag=temp[0]
        else:        
            color_tag=-temp[0]/max_comp*1
        ax1.text(np.mean(x),np.mean(y),round(temp[0],2))
        ax1.plot(x,y,color=[1,0,float(1-color_tag)],lw=area[e-1])

    ax2.plot(x,y,marker="o",color="0.5")       
    ax2.plot(u,v,color="0.",linestyle="dashed")#変形図をプロット

#軸ラベルを消す
ax1.tick_params(labelbottom=False,
                labelleft=False,
                labelright=False,
                labeltop=False,
                bottom=False,
                left=False,
                right=False,
                top=False)
ax2.tick_params(labelbottom=False,
                labelleft=False,
                labelright=False,
                labeltop=False,
                bottom=False,
                left=False,
                right=False,
                top=False)

#各グラフのフォーマットの設定
ax1.set_title("Stress Diagram")
ax1.set_xlim(-1,max_x)
ax1.set_ylim(-1,max_y)
ax1.spines["right"].set_color("none")
ax1.spines["left"].set_color("none")
ax1.spines["top"].set_color("none")
ax1.spines["bottom"].set_color("none")
ax2.set_title("Displacement Diagram")
ax2.set_xlim(-1,max_x)
ax2.set_ylim(-1,max_y)
ax2.spines["right"].set_color("none")
ax2.spines["left"].set_color("none")
ax2.spines["top"].set_color("none")
ax2.spines["bottom"].set_color("none")
plt.show()

2020/04/05

2次元トラス解析プログラムの構築(OpenSeesPy)

OpenSeesなどを用いて建築物などの構造解析を行う場合、様々な形状の構造物の応力解析を行う必要があり、そうした構造物を想定する場合どうしても入力する必要のある情報(構造物の節点、部材構成など)の量が増えることになります。そのため、構造物の情報を整理するプログラムをOpenSeesPyモジュールと併用しながら構成すると便利です。そこで、今回は以前の例題で用いたコードを抽象化し、構造物の情報を外部ファイルからインプットしながら整理しOpenSeesPyモジュールに受け渡すソースを組み立て、任意形状の入力に対応した2次元のトラス解析プログラムを構築します。

想定するプログラムのinput、outputは以下の通りです。
input:構造物の節点座標、部材構成、部材断面、節点荷重の位置、場所
output:各節点の座標、部材応力

前回までに用いたトラス構造物を想定した場合、inputファイルはcsvファイルで以下のように構成します。



<ポイント>

・input.csvの読み込み

作成したcsvファイルを扱うため、Pandasモジュールを使用します。Pandasはpythonにて、Excel、csvファイルのデータフレームを扱えるオープンソースのライブラリです。このモジュールを用いることで、csvファイルのデータを読み込むことが可能になります。m使用の際には、別途PCにPandasモジュールをダウンロードし、プログラムの中でインポートする必要があります。
データの読み込み部分は以下のように構成します。

#csvファイルから構造物のデータを読み込み
ip = pd.read_csv('input.csv', header=1)
no=ip['no']
・・・・

上記1行目では、csvファイルのヘッダーを一行読み飛ばし、2行目から各列の名目とデータを取得しています。
2行目はpythonのリストに'no'列(上図A列)のデータを読み込んでいます。
リストに読み込んでしまえば、リストのデータをOpenSeesPyのコマンドに受け渡すことにより、解析ファイルを生成することができます。

・NaN行の非カウント

OpenSeesPyのインプットを作成する際、データフレーム上のNaN行を読み飛ばす必要があります。そこで各インプットデータ作成ループの前に、NaNがない場合に限り、必要ループ数をカウントする構文を追記しました。NaNの判定は、以下のように同じ値を比較することで行えます。

count1=0
for i in range(len(fix_node)):#NaNの場合をカウントしない
    if fix_node[i] == fix_node[i]:
        count1+=1

・解析結果の出力

以下のようにfor文を用いて、input.csvの節点数、部材数に応じて動的に結果出力を行えるようにしました。

# output_results
print("<Node Displacement>")
for i in range(len(no)):
    print("Node%s x:%s y:%s" %(i,nodeDisp(int(no[i]),1),nodeDisp(int(no[i]),2)))
print("<Axial Force of each element>")
for i in range(count2):
    print("N%s=%s" %(i,basicForce(int(ele_no[i]))))

前回対象としたトラス構造物をインプットとした場合、出力される結果は以下の通りになります。しかし、数値だけだと解析結果のイメージが直感的につかめません。そこで、次回はpythonのグラフ描画モジュールmatplotlibを用いて、解析結果を図化出力することにトライしてみようと思います。


今回の検討に用いたソースは以下の通りです。赤字が今回追記した部分になります。


from openseespy.opensees import *

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# ------------------------------
# Start of model generation
# -----------------------------

# remove existing model
wipe()

# set modelbuilder
model('basic', '-ndm', 2, '-ndf', 2)

#input information
ip = pd.read_csv('input.csv', header=1)
no=ip['no']
cor_x=ip['cor_x']
cor_y=ip['cor_y']
fix_node=ip['fix_node']
fix_x=ip['fix_x']
fix_y=ip['fix_y']
ele_no=ip['ele_no']
n1=ip['n1']
n2=ip['n2']
area=ip['area']
load_point=ip['point']
dir_x=ip['dir_x']
dir_y=ip['dir_y']
dat_no=len(ip)#データフレーム数を取得

# create nodes
count=0
for i in range(dat_no):#NaNの場合をカウントしない
    if no[i] == no[i]:
        count+=1
for i in range(count):
    node(int(no[i]), float(cor_x[i]), float(cor_y[i]))

# set boundary condition
count1=0
for i in range(dat_no):#NaNの場合をカウントしない
    if fix_node[i] == fix_node[i]:
        count1+=1
for i in range(count1):
    fix(int(fix_node[i]), int(fix_x[i]), float(fix_y[i]))

# define materials
uniaxialMaterial("Elastic", 1, 3000.0)

# define elements
count2=0
for i in range(dat_no):#NaNの場合をカウントしない
    if ele_no[i] == ele_no[i]:
        count2+=1
for i in range(count2):
    element("Truss",int(ele_no[i]),int(n1[i]),int(n2[i]),float(area[i]),1)

# create TimeSeries
timeSeries("Linear", 1)

# create a plain load pattern
pattern("Plain", 1, 1)

# Create the nodal load - command: load nodeID xForce yForce
count3=0
for i in range(dat_no):#NaNの場合をカウントしない
    if load_point[i] == load_point[i]:
        count3+=1
for i in range(count3):
    load(int(load_point[i]),float(dir_x[i]),float(dir_y[i]))

# ------------------------------
# Start of analysis generation
# ------------------------------

# create SOE
system("BandSPD")

# create DOF number
numberer("RCM")

# create constraint handler
constraints("Plain")

# create integrator
integrator("LoadControl", 1.0)

# create algorithm
algorithm("Linear")

# create analysis object
analysis("Static")

# perform the analysis
analyze(1)

# output_results
max_tens=0
max_comp=0

print("<Node Displacement>")
for i in range(count):
    print("Node%s x:%s y:%s" %(i,nodeDisp(int(no[i]),1),nodeDisp(int(no[i]),2)))
print("<Axial Force of each element>")
for i in range(count2):
    print("N%s=%s" %(i,basicForce(int(ele_no[i]))))
    temp=basicForce(int(ele_no[i]))
    if temp[0] >= 0:
        if temp[0] > max_tens:
            max_tens = temp[0] 
    else:
        if abs(temp[0]) > max_comp:
            max_comp = abs(temp[0])

2020/03/29

OpenSeesPyによる静定トラスの解析(材料非線形性の考慮)

通常、構造物の材料は材料非線形性を有しており、部材荷重がある降伏荷重を超えると、部材が降伏することで部材の荷重増加率が低減し、部材の変形が大きく増加する状態になります(降伏現象)。Openseesではこうした部材の材料非線形性を考慮した応力解析が可能です。今回は前回と同形状の静定トラスを対象として、その部材が材料非線形性を考慮したモデルを題材とし、各部材の降伏を追跡した応力解析による検討を行います。

今回考慮する部材の材料非線形性は下図に示すような1つの降伏点を有するバイリニアの材料非線形性とします。部材の降伏荷重sYや降伏後の部材剛性低下率αを定義することによって、部材の材料非線形性を解析において考慮できます。


今回は以下に示すパラメータを有する架空の材料を想定します。
ヤング係数: E = 3000
降伏荷重: sY=10
剛性低下率: α=0.05

OpenSeesPy上で上記の材料特性を記載すると、以下のようになります。

E=3000
sY = 10
alpha = 0.05
uniaxialMaterial("Steel01", 1, sY, E, alpha)

()内は部材タイプ、部材特性の番号、降伏荷重、ヤング係数、剛性低下率の順に記載します。上記のコマンドを入力することで、タグ1が定義されている部材全てにおいて、設定した非線形の材料特性が適用されます。
なお、「uniaxialMaterial」はトラスなどの軸方向のみ剛性を有する部材の材料特性を定義するコマンドです。今回導入するバイリニアの材料特性以外にも様々な部材タイプの材料特性が用意されています。それぞれの材料特性に対し設定すべきパラメータが定義されており、それらを設定することにより様々な材料特性を考慮した解析を行うことが可能です。詳しくは以下のリンク先をご覧ください。また、そのうち整理してみたいと思います。(コンクリートなどの材料特性を想定したモデルも用意されています。世界中の研究者が提案している最新の材料特性モデルなども多く揃えられています。)

https://openseespydoc.readthedocs.io/en/latest/src/uniaxialMaterial.html

基本的には上述したコマンドを加えるだけで解析モデルにおける材料非線形性が定義され、部材の材料非線形性を考慮した応力解析が可能です。さらに、今回はトラス構造物における各部材の降伏に伴う節点変位を追跡するため、設定した節点荷重に至るまでの各荷重係数(λ=0~1)を想定した場合の結果を出力する荷重増分解析を行います。そのため、解析実施部分のスクリプトを以下のように定義します。

integrator("LoadControl", 1.0/Nsteps)

integratorコマンドを入力することで静的荷重増分解析を行うことができます。増分解析を荷重制御にて行う場合,上記のように記述し解析の分割数Nsteps=100を想定します。また各stepにおける頂部節点変位と外力を出力するため、以下のようにfor文を記述することにより、解析の実施及び解析結果の出力を行います。

Nsteps=100
data = np.zeros((Nsteps+1,2))
for j in range(Nsteps):
    analyze(1)
    data[j+1,0] = -nodeDisp(2,2)
    data[j+1,1] = getLoadFactor(1)*150

nodedispコマンドはnodeDisp(対象節点番号、変位方向)を入力することで各節点の変位を取得できます。また、getloadfactorコマンドでは各解析時の荷重係数をそれぞれ取得しています。結果リストをグラフ描画モデュールmatplotlibにより処理することで、以下のような荷重変位関係が取得できます。

グラフより節点の鉛直荷重(Vertical Load)が増加することにより、部材の降伏荷重に達し、その後節点の鉛直変形(Vertical Displacement)が大きく増加していることが確認できます。また、対象トラス構造物における部材1,3の降伏やその後の部材2の降伏を荷重増分解析により適切に追跡できていることが確認できます。


以下に、今回の検討に用いたソースコードを示します。基本は前回のソースコードと同様ですが、赤字が主に変更した部分になります。


from openseespy.opensees import *

# ------------------------------
# モデルの作成
# -----------------------------

# モデルの初期化
wipe()

# モデルの定義
model('basic', '-ndm', 2, '-ndf', 2)

# 節点の定義
node(1, 0.0, 0.0)
node(2, 4.04.0)
node(3, 8.00.0)

# 境界条件の定義
fix(1, 1, 1)
fix(3, 0, 1)

# 部材材料の定義
E=3000
sY = 10
alpha = 0.05
uniaxialMaterial("Steel01", 1, sY, E, alpha)

# 部材の定義
element("Truss",1,1,2,5.0,1)
element("Truss",2,1,3,5.0,1)
element("Truss",3,2,3,5.0,1)

# 時間に依存する荷重係数の定義
timeSeries("Linear", 1)

# 荷重パターンの定義
pattern("Plain", 1, 1)

# 節点荷重の定義
load(2, 0, -150)

# ------------------------------
# 解析条件の定義、解析実行
# ------------------------------

# create SOE
system("BandSPD")

# create DOF number
numberer("RCM")

# create constraint handler
constraints("Plain")

# create integrator
Nsteps = 100
integrator("LoadControl", 1.0/Nsteps)

# create algorithm
algorithm("Newton")

# create analysis object
analysis("Static")

# perform the analysis
data = np.zeros((Nsteps+1,2))
for j in range(Nsteps):
    analyze(1)
    data[j+1,0] = -nodeDisp(2,2)
    data[j+1,1] = getLoadFactor(1)*150

# ------------------------------
# 解析結果のグラフ出力(Matplotlibによる)
# ------------------------------
plt.plot(data[:,0], data[:,1])
plt.xlabel('Vertical Displacement')
plt.ylabel('Vertical Load')
plt.show()

2020/03/28

OpenSeesPyによる静定トラスの解析

今回は図に示すように、静定トラスの頂部に静的な荷重が加わった場合の応力解析をOpenSeesPyを用いて行っていきます。


・解析モデルの設定
modelコマンドにて全体の解析モデルを定義します。
ndmは解析モデルの次元数、ndfは自由度数です。
今回は2次元の静定トラスを想定するので、いずれも2に設定します。

model('basic', '-ndm', 2, '-ndf', 2)

・トラスの節点座標の設定
nodeコマンドにて解析モデル各節点の座標を定義します。
2次元モデルの場合、node(節点番号、x座標、y座標)のように定義します。

・固定されている節点の境界条件の設定
fixコマンドにて解析モデルの各節点の境界条件を定義します。
2次元モデルの場合、fix(節点番号、x方向の境界条件、y方向の境界条件(いずれも0:フリー,1:固定)と定義します。

・トラス部材の材料の条件を設定
今回は弾性のトラス部材を想定するため、UniaxialMaterialコマンドにて部材のヤング係数のみを定義します。

・トラス部材を設定
elementコマンドにて各部材の特性を定義します。
トラス部材の場合、以下のように定義できます。
element(elementtype、部材番号、部材端部の節点番号、断面積、定義した材料番号)

・外力を設定
節点に生じる外力をloadコマンドにて以下のように定義します。
load(設定する節点番号、x方向の外力、y方向の外力)

上記が対象とする問題を想定した場合の主な解析モデル、荷重条件の設定項目になり、解析に関する種々の条件を定義したうえでanalyzeコマンドを実行すると、応力解析が実行されます。解析結果は「Output Commands」に用意されている各種コマンドを用いて確認することができます。例えば、確認したい項目が、節点変位であればnodeDispコマンド、部材応力であればbasicForceコマンドの結果をprintすることで確認できます。今回の場合、basicForceコマンドにより出力された値が手計算の値と一致することから、正しく解析が行われていることが確認できました。



以下、今回の検討に用いたソースコードを示します。


from openseespy.opensees import *


# ------------------------------
# モデルの作成
# ----------------------------

# モデルの初期化
wipe()

# モデルの定義
model('basic', '-ndm', 2, '-ndf', 2)


# 節点の定義
node(1, 0.0, 0.0)
node(2, 4.04.0)
node(3, 8.00.0)

# 境界条件の定義
fix(1, 1, 1)
fix(3, 0, 1)

# 部材材料の定義
uniaxialMaterial("Elastic", 1, 3000.0)

# 部材の定義
element("Truss",1,1,2,5.0,1)
element("Truss",2,1,3,5.0,1)
element("Truss",3,2,3,5.0,1)

# 時間に依存する荷重係数の定義
timeSeries("Linear", 1)

# 荷重パターンの定義
pattern("Plain", 1, 1)

# 節点荷重の定義
load(2, 0, -50)

# ------------------------------
# 解析条件の定義、解析実行
# ------------------------------

# create SOE
system("BandSPD")

# create DOF number
numberer("RCM")

# create constraint handler
constraints("Plain")

# create integrator
integrator("LoadControl", 1.0)

# create algorithm
algorithm("Linear")

# create analysis object
analysis("Static")

# perform the analysis
analyze(1)

# ------------------------------
# 解析結果の確認
# ------------------------------

N1=basicForce(1)
N2=basicForce(2)
N3=basicForce(3)

print("N1=",N1,"N2=",N2,"N3=",N3)

OpenSeesPyの導入②

今回は、実際にpythonにおいてOpenSeesPyを扱う上で必要なポイントについて説明したいと思います。

・Pythonプログラムの作成環境
前回インストールしたopensees anacondaにはプログラムエディタとして、JupyterNotebookとSpyderが付属しています。

・JupiterNotebook
Webブラウザ上で動作し実行結果がすぐにみられる対話型実行環境。

・Spyder
科学者、エンジニア向けにデザインされたPythonの統合開発環境。テキストエディタ形式でプログラムの編集および実行が可能。


Pythonはコンパイルを不要とするインタプリタ言語であるため、いずれもプログラム記述後すぐに実行結果がみられます。本ブログの検討では基本的にSpyderをベースとした環境にてプログラムを構築していきます。

・Pythonプログラム上でのOpenSeesPyのインポート
Pythonプログラムには様々な機能を提供する標準モジュールが「ライブラリ」として用意されており、それらを自身の作成プログラムにimportすることにより、シンプルなコーディングで様々な機能を有するプログラムを作成することが可能になります。「OpenSeesPy」もそうしたライブラリの一つであり、プログラムの冒頭で以下のようにimportすることにより、使用が可能になります。


・OpenSeesの単位系について
OpenSeesはアメリカ合衆国のPacific Earthquake Engineering Research Center (PEER) によって開発されたため、基本的にヤードポンド法をベースとしており、使用の際にはSI単位系に換算したパラメータを用いる必要があります。 

2020/03/22

OpenSeesPyの導入①

Openseesはアメリカ合衆国のPacific Earthquake Engineering Research Center (PEER) が開発したオープンソース型の構造解析ツールです。オープンソースであることから、全てのコードは公開されています。またその内容は世界中のエンジニアの手で日々更新されており、更新履歴は全てweb上で公開されています。
静的な応力を対象とした構造解析に加え、地震波を入力した時刻歴応答解析などが可能で、全てフリーで使えることから世界中の研究者が解析ツールとして利用しています。詳細は以下をご覧ください。

https://opensees.berkeley.edu/

近年、pythonなどが気軽に使えるプログラムとして台頭しており、openseesもpythonプログラム上で導入可能なpythonモジュール「OpenSeesPy」として公開されています。
各コマンドの機能を解説したマニュアルについても、英語版は下記リンク先にて公開されています。

今回はまずその導入方法を紹介したいと思います。

①pythonを導入
pythonの主要なモデュールを含んだanacondaを以下ページより導入します。

https://www.anaconda.com/distribution/#download-section

開いたページの「Download」をクリックします。

②Jupyter Notebookを起動
Jupyter Notebookはpythonコードを書くと、その実行結果がすぐにみられるツールです。同ツールを開き、ツール上でそれぞれの場合において以下のコードを記述します。


・新規インストール
python -m pip install openseespy
python -m pip install --user openseespy
・ソフトウェアの更新
python -m pip install --upgrade openseespy
python -m pip install --user --upgrade openseespy
上記のように入力すると自動的にモデュールのインストールが始まります。

③pythonコード上にモデュールをインポート
以下のコードをPythonコードの冒頭に記載することで、python上でOpenSeesPyモデュールの使用が可能となります。

import openseespy.opensees as ops

上記の操作により、Pythonコーディング上でOpenseesの構造解析、時刻歴応答解析ツールを動かすための環境構築が終了しました。次は、Python上でOpenseesを動かす方法について紹介していきたいと思います。