見出し画像
Photo by io28

OpenseesPy静的解析時の荷重・荷重制御設定のメモ

watada

OpenseesPyでの荷重の設定や途中での変更方法についてもやもやし、調べた部分の個人用メモ。OpenSees、OpenSeesPyマニュアルは以下。


前提

・analysisは'Static'。patternは'Plain'。
・analyze(1)実行することを1解析ステップと呼ぶ。$${i}$$は1からはじまりanalyze実行毎に1加算される解析ステップ数とする。
・呼称・記号などは公式ドキュメントなどとは整合していません。説明のためここで独自に定義したものも多く、その明示もしていません。
・誤りなども含むと思います。参考にされる場合は自己責任で。

荷重の重複定義について

loadは設定済loadからの差分（加算式）

load(nodeTag, *loadValues)は設定済loadからの差分$${\Delta \bar{\boldsymbol{F}}_i}$$。コード中に同じnodeTagへの重複定義があった場合、それまでにそのnodeTagに定義していたすべてのloadValuesに加算し、次解析ステップ用の基準荷重ベクトル$${\bar{\boldsymbol {F}}_i}$$を設定する（上書きではない）。

spは上書き

loadの代わりにsp（single-point constraint）による強制変位を使用する場合、コード中に同じnodeTag, dofへのspの重複定義があった場合には、それまでに当該nodeTag, dofに定義していたdofValueを上書きし、次解析の強制変位設定時の（基準）指定移動先座標$${\bar{\boldsymbol{x}}_{target,i}}$$として使用する。
なお、spを使用する場合はconstraintsは'Transformation'などを使用する（'Plain'不可）。

荷重制御方法（integrator），time，timeseriesについて

integrator('LoadControl', incr)のincrは次解析ステップでのtimeの増分（荷重係数の増分ではない）

integrator('LoadControl', incr)のincrは直前のanalyze実行状態からのtime$${t_i}$$の増加分$${\Delta t_i}$$。ここから次ステップ解析時の時間$${t_i=\sum \Delta t_i}$$が定まる。recorder('Node')等'-time'指定時に出力される1列目の値はこの$${t_i}$$。

timeとtimeseriesの設定で、次解析時の荷重係数lambdaが定まる

timeseriesはこの$${t_i}$$と次解析時の荷重係数$${\lambda_i}$$を関連付ける関数$${f}$$。
・'Linear': $${\lambda_i =f(t_i)=cFactor * (t_i -tStart)}$$
・'Constant': $${\lambda_i=f(t_i)=C}$$
　他　https://openseespydoc.readthedocs.io/en/latest/src/timeSeries.html

次解析ステップで実際に構造物に作用させる荷重$${\boldsymbol{F}_i}$$は、$${\lambda_i}$$と基準荷重ベクトル$${\bar{\boldsymbol {F}}_i}$$の積（ $${\boldsymbol{F}_i=\lambda_i \bar{\boldsymbol {F}}_i}$$）。

荷重としてspを指定した場合、前述の基準移動先座標$${\bar{\boldsymbol{x}}_{target}}$$の各成分を$${\lambda_i}$$倍した$${\boldsymbol{x}_{target,i}=\lambda_i \bar{\boldsymbol{x}}_{target,i}}$$を次解析時の指定移動先座標とする強制変位を与える（前ステップ座標との内外分値では無く0との内外分値になるため、注意が必要）。

integrator('DisplacementControl')使用時のtime, timeseries, lambdaについて

integrator('DisplacementControl', nodeTag, dof, incr)のincrは前解析ステップ終了時からの変位増分の制御量$${\Delta U_{dof}}$$で、そのときの荷重ベクトル$${\bar{\boldsymbol {F}}_i}$$は設定済のops.loadの総和。次解析ステップでは指定の$${\Delta U_{dof}}$$を達成する荷重$${\boldsymbol {F}_i=\lambda_i \bar{\boldsymbol {F}}_i}$$を定めるべく荷重係数$${\lambda_i}$$の収束計算を行う。
なお、この場合の各ステップ毎のtime $${t_i}$$はtimeseriesで指定した$${\lambda_i=f(t_i)}$$から逆算される（物理的な意味はないpseudotime）。timeseriesを'Constant'にしている場合など、対応する$${t_i}$$が逆算により定まらない場合エラーになるっぽい。

timeの設定・リセットにはsetTime、確認にはgetTimeを用いる

integratorに'LoadControl'を指定する場合、timeseriesに'Constant'を指定している場合を除き、上記のとおり、time $${t_i}$$の設定が荷重係数$${\lambda_i}$$に影響を及ぼすため、途中で載荷方向や条件を変更する場合、timeseriesのtStartを適切に設定するか、都度timeを変更・リセットするなどが必要。timeの変更にはsetTime(pseudoTime)を用いる（リセットする場合はpseudoTime=0とする）。途中でその時点でのtimeの値を確認するには、getTime()を用いる。

荷重のPatternについて

新規のpattern('Plain', patternTag, tsTag)を宣言した場合、そのスコープは次のpattern('Plain', patternTag, tsTag)が宣言されるまでに追加されるload,spに有効。また、一度あるpatternTagで定義されたloadは、別のpatternTagが宣言された後も、引き続き作用しつづける（消えずに加算になる）。削除するにはremove('loadPattern', patternTag)を指定する。これを指定すると、次解析以降、指定したpatternTagで定義されていた荷重は削除される。
例えば、2次元問題で節点nodeTagに鉛直荷重[0,-100]（patternTag=1, timeseriesは'Constant'を指定）をステップ1-5の5ステップかけて線形載荷させたあと、追加で水平荷重[20,0]（patternTag=2, timeseriesは'Linear'を指定）を10ステップ（ステップ6-15）かけて漸増させる場合には、以下のようにすればよい。

import openseespy.opensees as ops
import numpy as np

#(中略)

# ステップ1-5
tsTagC=1
tsTagL=2
ops.timeSeries('Constant', tsTagC)
ops.timeSeries("Linear", tsTagL)

stepv=5
patternTag1=1
ops.pattern("Plain", patternTag1, tsTagC)
ops.integrator('LoadControl', 1.0/stepv)
for i in range(1, stepv+1):
    ops.load(nodeTag, *[0, -100/stepv, 0])
    ops.analyze(1)

# ステップ6-15
ops.setTime(0.0)
steph=10
patternTag2=2
ops.pattern("Plain", patternTag2, tsTagL)
ops.load(nodeTag, *[20, 0, 0])
ops.integrator('LoadControl', 1.0/steph)
for i in range(1, steph+1):
    ops.analyze(1)

ステップ1-5では'Constant'を指定しているので、ops.integrator('LoadControl', 1.0/stepv)の第2引数は結果には影響しないが、経過するpseudotimeには影響する。上の例では5ステップ終了時点で出力されるpseudotime=1.0となる。
一方、ステップ6-15の間に指定した荷重はすべて'Linear'で指定しているため、このあと（ステップ16以降など）にさらに処理を追加する場合、その時点のtimeの影響を受ける。removeなどを利用して、すべて'Constant'のpatternに置き換えておくほうが扱いやすい。例えば、上の処理のあとに続けて、ステップ16-30でY方向の変位を0まで変位制御で戻すには、以下の処理を追記する。

# ステップ16-30
# ここまでに作用している荷重を一旦削除
ops.setTime(0.0)
ops.remove('loadPattern', patternTag1)
ops.remove('loadPattern', patternTag2)

# 削除した荷重を'Constant'に置き換え再定義
patternTag3=3
ops.pattern('Plain', patternTag3, tsTagC)
ops.load(nodeTag, *[20,-100,0]) # ここまでに作用している荷重を'Constant'に置き換え

# Y方向変位除去用基準荷重を'Linear'パターンで追加定義
patternTag4=4
ops.pattern('Plain', patternTag4, tsTagL)
ops.load(nodeTag,*[0,1.0,0]) 

peak_disp=ops.nodeDisp(nodeTag,2) # ステップ15終了時点のY方向変位を取得

stepv2=15
ops.integrator('DisplacementControl', nodeTag, 2, -peak_disp/stepv2)
for i in range(1,1+stepv2):
    ops.analyze(1)

同様に、ステップ31-38で節点nodeTagの座標をspで[1,0]から単位円周上に左回りに回転させるには、以下の処理を追記する。

# ステップ31-38
patternTag5=5
ops.pattern('Plain', patternTag5, tsTagC) # 'timeSeriesは'Constant'
ops.load(nodeTag,*[20*ops.getLoadFactor(3),-100*ops.getLoadFactor(3),0])
ops.load(nodeTag,*[0,1.0*ops.getLoadFactor(4),0])
ops.remove(patternTag3)
ops.remove(patternTag4)
ops.setTime(0.0)

stepr=8
ops.integrator('LoadControl', 1.0/stepr)
for i in range(stepr):
    ops.sp(nodeTag, 1, np.cos(2*np.pi*i/stepr))
    ops.sp(nodeTag, 2, np.sin(2*np.pi*i/stepr))
    ops.analyze(1)

以上の処理を含むOpenseesPyのサンプルコードをまとめると以下のとおりである。

import openseespy.opensees as ops
import numpy as np

# モデル定義と基本設定
ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3)  # 2次元モデル, 節点自由度は3 (UX, UY, ROTZ)

# 節点の設定
nodesupTag=1
nodeTag=2
ops.node(nodesupTag, 0.0, 0.0)
ops.node(nodeTag, 1.0, 0.0)

# 境界条件
ops.fix(nodesupTag, *[1, 1, 1])  # 節点1は全自由度固定

# 材料の定義
ops.uniaxialMaterial("Elastic", 1, 210000)  # 材料ID, ヤング率（kN/m^2）

# 座標変換の定義
ops.geomTransf("Linear", 1)  # 線形座標変換、タグ1

# 断面の定義
A = 0.02  # m^2
Iz = 0.0001  # m^4
Av = 0.01  # m^2, せん断面積

# 要素の定義
ops.element("elasticBeamColumn", 1, 1, 2, A, 210000, Iz, 1, '-mass', 0.0, '-cMass')

# Load Patternの定義
tsTagC=1
tsTagL=2
ops.timeSeries('Constant', tsTagC)
ops.timeSeries("Linear", tsTagL)

# 解析の設定
ops.system("BandGeneral")
ops.numberer("Plain")
ops.constraints("Transformation")
ops.test('NormDispIncr', 1.0e-4, 30, 0) 
ops.algorithm('Newton') 
ops.analysis('Static')

ops.recorder('Node', '-file', './reac.txt', '-time', '-node', *[nodesupTag], '-dof', *[1,2], 'reaction')
ops.recorder('Node', '-file', './disp.txt', '-time', '-node', *[nodeTag], '-dof', *[1,2], 'disp')

# ステップ1-5
stepv=5
patternTag1=1
ops.pattern("Plain", patternTag1, tsTagC)
ops.integrator('LoadControl', 1.0/stepv)
for i in range(1, stepv+1):
    ops.load(nodeTag, *[0, -100/stepv, 0])
    ops.analyze(1)

# ステップ6-15
ops.setTime(0.0)
steph=10
patternTag2=2
ops.pattern("Plain", patternTag2, tsTagL)
ops.load(nodeTag, *[20, 0, 0])
ops.integrator('LoadControl', 1.0/steph)
for i in range(1, steph+1):
    ops.analyze(1)

# ステップ16-30
# ここまでに作用している荷重を一旦削除
ops.setTime(0.0)
ops.remove('loadPattern', patternTag1)
ops.remove('loadPattern', patternTag2)

# 削除した荷重を'Constant'に置き換え再定義
patternTag3=3
ops.pattern('Plain', patternTag3, tsTagC)
ops.load(nodeTag, *[20,-100,0]) # ここまでに作用している荷重を'Constant'に置き換え

# Y方向変位除去用基準荷重を'Linear'パターンで追加定義
patternTag4=4
ops.pattern('Plain', patternTag4, tsTagL)
ops.load(nodeTag,*[0,1.0,0]) 

peak_disp=ops.nodeDisp(nodeTag,2) # ステップ15終了時点のY方向変位を取得

stepv2=15
ops.integrator('DisplacementControl', nodeTag, 2, -peak_disp/stepv2)
for i in range(1,1+stepv2):
    ops.analyze(1)

# ステップ31-38
patternTag5=5
ops.pattern('Plain', patternTag5, tsTagC) # 'timeSeriesは'Constant'
ops.load(nodeTag,*[20*ops.getLoadFactor(3),-100*ops.getLoadFactor(3),0])
ops.load(nodeTag,*[0,1.0*ops.getLoadFactor(4),0])
ops.remove(patternTag3)
ops.remove(patternTag4)
ops.setTime(0.0)

stepr=8
ops.integrator('LoadControl', 1.0/stepr)
for i in range(stepr):
    ops.sp(nodeTag, 1, np.cos(2*np.pi*i/stepr))
    ops.sp(nodeTag, 2, np.sin(2*np.pi*i/stepr))
    ops.analyze(1)

これを実行して得られるnodeTagの変位と支点の反力を出力したものは以下のようになる。いずれも行番号が解析ステップ番号$${i}$$、1列目の値がtime $${t_i}$$を表す。

なおOpenseesPyのマニュアルによれば、setTime(0.0)を指定するかわりにops.timeSeries('Linear', tsTagL)にtStartオプションを指定する使い方もあると思うが、手元ではうまくいかない。原因は不明。

まとめ：どうするのが良い？

静的解析の場合、以下の組み合わせなどが考えられる。
※参考（呼称もここから）https://www.youtube.com/watch?v=O2f8ihoN7Ts

Load control with force / LcF

荷重：load(nodeTag, *loadValues)
制御：integrator('LoadControl', incr)
※timeseriesの設定に注意

Load control with displacement / LcD

荷重：sp(nodeTag, dof, dofValue)
制御：integrator('LoadControl', incr)
※移動前のdofValueが0でない場合、incr=1.0, timeseriesは'Const'を用いる

Displacement control with force / DcF or DcL

荷重：load(nodeTag, *loadValues)
制御：integrator('DisplacementControl', nodeTag, dof, incr, dUmin=incr, dUmax=incr)
※Cyclic Analysisなど変位の進行方向の正負が入れ替わる点では、incrの符号やdUmin, dUmaxの値は変更する。
※timeseriesは'Linear'としておく（timeに物理的な意味はない）。

備考・その他

・途中で荷重条件を変更する場合、setTimeで時間をリセットするほうが扱いやすい。
・eleloadや他のpattern、Arclengthなど他の制御方法、動的解析については未確認。確認できたら追記予定。
・その他、参考サイト：
https://portwooddigital.com/2021/01/10/load-patterns-and-time-series/


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