Pythonで２次元弾性骨組幾何学的非線形解析 #Python

プログラム概要

このプログラムは，要素微小変位の仮定において２次元弾性平面骨組の幾何学的非線形解析を行うものです．
要素として，1要素2節点，1節点3自由度（水平変位・鉛直変位・回転）を有する梁要素が用いられています．
荷重として，節点外力増分のみを扱えます．
釣り合い条件を求めるための解析手法として，弧長増分法を用いています．これにより変位は増大するが荷重が減少するような現象も追跡できます．
座標システムにおいては，右方向がx軸，上方向がy軸，反時計回りが正の回転方向です．
連立一次方程式は，numpy.linalg.solve(A, b) で解いています．
入力データファイル，出力データファイルは，空白区切りで記載される必要があり，ファイル名はコマンドライン引数として入力します．

扱える境界条件
項目	説明
節点外力増分	載荷節点数，載荷節点番号，荷重増分値を指定
節点変位拘束	変位拘束節点数，拘束節点番号を指定（完全拘束：変位=0のみ扱えます）

FEM解析プログラム

Program name	Description
py_fem_gfrmAL.py	平面骨組幾何学的非線形構造解析プログラム

FEM解析実行用スクリプト

python3 py_fem_gfrmAL.py inp.txt out.txt nnmax

inp.txt	: 入力データファイル（空白区切りデータ）
out.txt	: 出力データファイル（空白区切りデータ）
nnmax	: 載荷ステップ数

入力データ書式

npoin  nele  nsec  npfix  nlod  # Basic values for analysis
E  A  I                         # Material properties
    ..... (1 to nsec) .....     # 
node_1  node_2  isec            # Element connectivity, material set number
    ..... (1 to nele) .....     # 
x  y                            # Node coordinate of node
    ..... (1 to npoin) .....    # 
lp  fix_x  fix_y  fix_r         # Restricted node number
    ..... (1 to npfix) .....    # 
lp  df_x  df_y  df_r            # Loaded node and loading conditions (load increment)
    ..... (1 to nlod) .....     #     (omit data input if nlod=0)

npoin, nele, nsec	: 節点数，要素数，断面特性数
npfix, nlod	: 拘束節点数，載荷節点数
E, A, I	: 弾性係数，断面積，断面2次モーメント
node_1, node_2, isec	: 節点1，節点2，断面特性番号
x, y	: x座標，y座標
lp, fix_x, fix_y, fix_r	: 節点番号，x・y・回転方向の拘束の有無（0: 自由，1: 完全拘束）
lp, df_x, df_y, df_r	: 節点番号，x・y・回転方向荷重

出力データ書式

npoin  nele  nsec npfix  nlod  nnmax
    (Each value of above)
sec  E  A  I
    sec   : Material number
    E     : Elastic modulus
    A     : Section area
    I     : Moment of inertia
    ..... (1 to nsec) .....
node  x  y  fx  fy  fr  kox  koy  kor
    node   : Node number
    x      : x-coordinate
    y      : y-coordinate
    fx     : Load in x-direction
    fy     : Load in y-direction
    fr     : Moment load
    kox    : Index of restriction in x-direction (0: free, 1: fixed)
    koy    : Index of restriction in y-direction (0: free, 1: fixed)
    kor    : Index of restriction in rotation    (0: free, 1: fixed)
    ..... (1 to npoin) .....
elem  i  j  sec
    elem : Element number
    i    : Node number of start point
    j    : Node number of end point
    sec  : Material number
    ..... (1 to nele) .....
* nnn=0 iii=0 lam=0.0
node  fp-x  fp-y  fp-r  dis-x  dis-y  dis-r  dr-x  dr-y dr-r
    node  : Node number
    fp-x  : Total load in x-direction
    fp-y  : Total load in y-direction
    fp-r  : Total load in rotation
    dis-x : Displacement in x-direction
    dis-y : Displacement in y-direction
    dis-r : Displacement in rotation
    dr-x  : Un-balanced force in x-direction
    dr-y  : Un-balanced force in y-direction
    dr-r  : Un-balanced force in rotation
    ..... (1 to npoin) .....
elem  N_i  S_i  M_i  N_j  S_j  M_j
    elem : Element number
    N_i  : Axial force of node-i
    S_i  : Shear force of node-i
    M_i  : Moment of node-i
    N_j  : Axial force of node-j
    S_j  : Shear force of node-j
    M_j  : Moment of node-j
    ..... (1 to nele) .....
* nnn=1 iii=xx lam=xxx
node  fp-x  fp-y  fp-r  dis-x  dis-y  dis-r  dr-x  dr-y dr-r
    node  : Node number
    fp-x  : Total load in x-direction
    fp-y  : Total load in y-direction
    fp-r  : Total load in rotation
    dis-x : Displacement in x-direction
    dis-y : Displacement in y-direction
    dis-r : Displacement in rotation
    dr-x  : Un-balanced force in x-direction
    dr-y  : Un-balanced force in y-direction
    dr-r  : Un-balanced force in rotation
    ..... (1 to npoin) .....
elem  N_i  S_i  M_i  N_j  S_j  M_j
    elem : Element number
    N_i  : Axial force of node-i
    S_i  : Shear force of node-i
    M_i  : Moment of node-i
    N_j  : Axial force of node-j
    S_j  : Shear force of node-j
    M_j  : Moment of node-j
    ..... (1 to nele) .....

.... .... ....

* nnn=nnmax-1 iii=xx lam=xxx
node  fp-x  fp-y  fp-r  dis-x  dis-y  dis-r  dr-x  dr-y dr-r
    ..... (1 to npoin) .....
elem  N_i  S_i  M_i  N_j  S_j  M_j
    ..... (1 to nele) .....
n=(total degrees of freedom)  time=(calculation time)

fp-x, fp-y, fp-r	: x方向荷重，y方向荷重，回転方向荷重
dis-x, dis-y, dis-r	: x方向変位，y方向変位，回転変位
dr-x, dr-y, dr-r	: x方向不平衡力，y方向不平衡力，回転方向不平衡力
N, S, M	: 軸力，せん断力，モーメント
n	: 全自由度（連立方程式の元）
time	: 計算時間

出力事例

Program name	Description
inp_gfrm_canti.txt	片持梁 FEM解析入力データ
inp_gfrm_cable.txt	ケーブル付き片持梁 FEM解析入力データ
inp_gfrm_arch.txt	アーチ FEM解析入力データ
inp_gfrm_lee.txt	フレーム FEM解析入力データ
py_fig_gfrm.py	荷重ー変位曲線作図プログラム（matplotlib）
py_fig_gfrm_mode.py	変位モード図作成プログラム（matplotlib）

荷重-変位曲線作図プログラム py_fig_gfrm.py

荷重変位曲線作成プログラム py_fig_gfrm.py では，汎用性を出そうとした結果，コマンドラインからの入力がかなりごたつくことになってしまいました．

変位モード図作成プログラム py_fig_gfrm_mode.py

変位モード図作成プログラム py_fig_gfrm_mode.py では，変位モードのみでは汎用性のあるものを比較的平易に実現できますが，境界条件を入れようとするとかなり複雑になってしまいます．よって，作成したいグラフは４枚のみなので，プログラム内で入力ファイルを指定し，ファイル毎に個別の処理をしています．

（荷重を示す矢印）

荷重を示す矢印は，arrowで書いています．arrowでは

ax.arrow(x,y,u,v,head_length=hl,  ....)

のような命令を書き込みます．例えば鉛直の矢印を描く場合，矢印の尻尾から先端までの長さは，( v + head_length ) となることに注意して扱う必要があります．実際の対応は以下の通り．

uu=0.0
vv=(ymax-ymin)*0.1
x1=xx[lnod-1]
y1=yy[lnod-1]+vv
hl=vv*0.4
hw=hl*0.5
ax.arrow(x1,y1,uu,-(vv-hl), lw=2.0,head_width=hw, head_length=hl, fc='#555555', ec='#555555')

（固定端を示す記号）

鉛直片道梁の基部に描いている固定端を示す記号は，fill(... hatch='///') を使って描いています．ここで，fillでは長方形領域を指定していますが，枠線は不要なので，fill の中で，linewidth=0.0 を指定して枠線を描画しないようにしています．実際の対応は以下の通り．

lp=1; x1=x[lp-1]; y1=y[lp-1]
px=[x1-2*scl,x1+2*scl,x1+2*scl,x1-2*scl]
py=[y1,y1,y1-1*scl,y1-1*scl]
ax.fill(px,py,fill=False,linewidth=0.0,hatch='///')
ax.plot([x1-2*scl,x1+2*scl],[y1,y1],color='#000000',linestyle='-',linewidth=1.5)

FEM計算実行・荷重変位曲線作図プログラム実行用スクリプト

# FEM calculation by Arc-Length method
python py_fem_gfrmAL.py inp_gfrm_canti.txt out_gfrm_canti.txt 70
python py_fem_gfrmAL.py inp_gfrm_cable.txt out_gfrm_cable.txt 290
python py_fem_gfrmAL.py inp_gfrm_arch.txt out_gfrm_arch.txt 310
python py_fem_gfrmAL.py inp_gfrm_lee.txt out_gfrm_lee.txt 160

# Drawing of Load-displacement curve
python py_fig_gfrm.py out_gfrm_canti.txt 11 11 -411.07 1000 1000 \$u/L\$ $\v/L\$ \$u/L\$\,$\v/L\$ \$P/P_{cr}\$ LL
python py_fig_gfrm.py out_gfrm_cable.txt 12 11 -1644.3 1000 1000 \$u/L\$ \$v/L\$ \$u/L\$\,\$v/L\$ \$P/P_{cr}\$ LL
python py_fig_gfrm.py out_gfrm_arch.txt 21 21 -666.4 -500 -500 \$u/R\$ \$v/R\$ \$u/R\$\,\$v/R\$ \$P\\cdot\(R^2/EI\)\$ UL
python py_fig_gfrm.py out_gfrm_lee.txt 13 13 -166.6 1000 -1000 \$u/L\$ \$v/L\$ \$u/L\$\,\$v/L\$ \$P\\cdot\(L^2/EI\)\$ LL

# Drawing of displacement mode
python py_fig_gfrm_mode.py

python py_fig_gfrm.py out.txt node-L node-D nd-L nd-u nd-v leg-u leg-v x-Label y-Label loc

out.txt	：FEM出力データファイル（空白区切りデータ）
node-L	：荷重変位曲線を描くための載荷節点番号
node-D	：荷重変位曲線を描くための変位描画節点番号
nd-L	：荷重無次元化のための数値（符号含む）
nd-u	：x方向変位無次元化のための数値（符号含む）
nd-v	：y方向変位無次元化のための数値（符号含む）
leg-u	：x方向変位ラベル（凡例用）
leg-v	：y方向変位ラベル（凡例用）
x-Label	：x軸ラベル（Displacement）
y-Label	：y軸ラベル（Load）
loc	：凡例描画位置

変位モードおよび荷重-変位曲線出力事例

軸圧縮力を受ける片持梁（inp_gfrm_canti.txt）
L=1,000mm, E=200,000MPa, A=100m$^2$, I=833m$^4$ Initial deflection $v_0$=1mm Buckling load $P_{cr}=\pi^2 EI / 4 L^2$=411.07N

ケーブルで先端を引っ張られる片持梁（inp_gfrm_cable.txt）
Column: L=1,000mm, E=200,000MPa, A=100m$^2$, I=833m$^4$ Cable: L=1,000mm, E=200,000MPa, A=28.3$^2$ Initial deflection $v_0$=5mm Buckling load $P_{cr}=\pi^2 EI / L^2$=1644.3NN

鉛直集中荷重を受ける非対称支持されたアーチ（inp_gfrm_arch.txt）
R=500mm, Center angle=215$^\circ$, E=200,000MPa, A=100m$^2$, I=833m$^4$

鉛直集中荷重を受けるフレーム（inp_gfrm_lee.txt）
L=1000mm, E=200,000MPa, A=100m$^2$, I=833m$^4$

以　上