大きな修正

(2017.05.18) 座標変換行列修正（部材軸が全体座標系Z軸と平行な場合とそうでない場合に場合分け）
(2018.08.08) 最新版プログラムはこちら：Python：FEMでの疎行列計算利用による高速化（３次元骨組構造解析での事例）
(2018.09.17) 最新版プログラム対応で記事の内容を修正．

概要

前から欲しかった立体骨組構造解析プログラムをPythonで作ってみました．
以前に作成した平面骨組解析プログラムにおいて，主に，剛性行列，座標変換行列を変更し，１節点６自由度の梁要素対応にしたものです．

~~まだテストが十分ではありませんが，~~ 下の記事に刺激されて投稿してみました．
http://qiita.com/sasaco/items/917429a497c6e9a08401

プログラム作成は，Jupyter上で行っていますので，この記事も自分のJupyterの記載に沿って記述しています．
Pythonによる解析プログラムは３分割してあります．
いつも気になるのですが，ファイル書き出しの部分がごちゃごちゃしすぎている．もっとうまく書けないものか．．．

~~計算はいちおう回っています．．．~~

プログラミング環境は
Machine: MacBook Pro (Retina, 13-inch, Mid 2014)
OS: macOS Sierra
Python version: Python 3.6.1

立体骨組構造解析の理論

基本事項

理論展開およびプログラムには，以下の条件が入っています．

線形弾性材料の微小変形を仮定している
物体力として慣性力を考慮する
初期ひずみとして温度変化によるひずみを考慮する

このプログラムでは，モデルが大規模になることを想定し，連立１次方程式を解く際，疎行列処理を導入しています．

要素剛性方程式

要素剛性方程式の一般形は，２次元応力解析におけるものと同一である．

変位を求めるための方程式

\begin{equation}
[\boldsymbol{k}]\{\boldsymbol{u_{nd}}\}=\{\boldsymbol{f}\}+\{\boldsymbol{f_t}\}+\{\boldsymbol{f_b}\}
\end{equation}

\begin{align}
&[\boldsymbol{k}]=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dV  & & \text{（要素剛性行列）}\\
&\{\boldsymbol{f}\}=\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA  & & \text{（節点外力ベクトル）}\\
&\{\boldsymbol{f_t}\}=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV  & & \text{（節点初期歪ベクトル）}\\
&\{\boldsymbol{f_b}\}=\gamma\left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} & & \text{（節点慣性力ベクトル）}
\end{align}

要素内部応力計算式

\begin{equation}
\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\}
\end{equation}

\begin{align}
&\{\boldsymbol{\epsilon}\}  & &\text{：剛性方程式を解いて求めた変位から求めた要素のひずみベクトル}\\
&\{\boldsymbol{\epsilon_0}\} & &\text{：温度ひずみなどの初期ひずみ}\\
&[\boldsymbol{D_e}]          & &\text{：応力ーひずみ関係行列}
\end{align}

立体骨組構造の要素剛性マトリックス

立体骨組構造の要素剛性マトリックスは 12 x 12 の正方行列であり，
これはねじれを考慮すること以外は平面骨組に対するものを拡張することにより導くことができるため，
導出過程は省略し，ここでは結果のみを記載する．

要素剛性方程式（最も単純な表現）

\begin{equation}
\boldsymbol{\{f\}}=\boldsymbol{[k]}\boldsymbol{\{u\}}
\end{equation}

節点変位ベクトル

\begin{equation}
\boldsymbol{\{u\}}=
\begin{Bmatrix}
u_i & v_i & w_i & \theta_{xi} & \theta_{yi} & \theta_{zi} & u_j & v_j & w_j & \theta_{xj} & \theta_{yj} & \theta_{zj}
\end{Bmatrix}^T
\end{equation}

\begin{align}
&u & &\text{：x軸方向変位} & &\theta_x& &\text{：x軸回り回転量（ねじり角）}\\
&v & &\text{：y軸方向変位} & &\theta_y& &\text{：y軸回り回転量}\\
&w & &\text{：z軸方向変位} & &\theta_z& &\text{：z軸回り回転量}
\end{align}

節点力ベクトル

\begin{equation}
\boldsymbol{\{f\}}=
\begin{Bmatrix}
N_i & S_{yi} & S_{zi} & M_{xi} & M_{yi} & M_{zi} & N_j & S_{yj} & S_{zj} & M_{xj} & M_{yj} & M_{zj}
\end{Bmatrix}^T
\end{equation}

\begin{align}
&N  & &\text{：軸力（x軸方向力）} & &M_x& &\text{：x軸回りモーメント（ねじりモーメント）}\\
&S_y& &\text{：y軸方向せん断力}   & &M_y& &\text{：y軸回り曲げモーメント}\\
&S_z& &\text{：z軸方向せん断力}   & &M_z& &\text{：z軸回り曲げモーメント}
\end{align}

要素剛性行列

\begin{equation}
\boldsymbol{[k]}=
\begin{bmatrix}
\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} & 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} \\
0 & 0 & \cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 & 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} \\
-\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} & 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} \\
0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 & 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} \\
\end{bmatrix}
\end{equation}

\begin{align}
&L& &\text{：要素長}        & &J& &\text{：ねじり定数}\\
&E& &\text{：弾性係数}      & &I_y& &\text{：y軸回り断面２次モーメント}\\
&G& &\text{：せん断弾性係数}& &I_z& &\text{：z軸回り断面２次モーメント}\\
&A& &\text{：断面積}        & &
\end{align}

立体骨組要素の節点慣性力ベクトル

物体力として慣性力を考え，節点ベクトルとして「加速度」を与える．
ここでは，最も単純な方法として，要素に作用する慣性力を単純に節点力に配分する．
要素に作用する全慣性力は，全体座標系での値として以下のとおりとなる．

\begin{align}
F_{m(x)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_x \cdot g) \\
F_{m(y)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_y \cdot g) \\
F_{m(z)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_z \cdot g) \\
\end{align}

よって，以下のように慣性力による節点力 ${\boldsymbol{f_{b}}}$ を定めることができる．

\begin{equation*}
\{\boldsymbol{F_{b}}\}=\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot A \cdot L}{2}
\begin{Bmatrix}
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0   \\
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0
\end{Bmatrix}
\end{equation*}

ここに，$\gamma$ は要素の単位体積重量，$A$ は要素断面積，$L$ は要素長さ，$g$ は重力加速度である．

立体骨組要素の節点温度荷重ベクトル

温度変化は部材の軸力のみに影響するとする．
ここで温度上昇を正の温度変化とすれば，温度変化による部材座標系での節点力ベクトルは以下のとおりとなる．

\begin{equation*}
\{\boldsymbol{f_{t}}\}=EA\cdot\alpha\cdot\Delta T
\begin{Bmatrix}
-1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0
\end{Bmatrix}
\end{equation*}

ここに，$EA$ は部材の軸剛性，$\alpha$ は部材の線膨張係数，$\Delta T$ は部材の温度変化量である．

これを座標変換することにより，全体座標系での節点力ベクトルが定まる．

\begin{equation*}
\{\boldsymbol{F_t}\}=[\boldsymbol{T}]^T\{\boldsymbol{f_t}\}
\end{equation*}

ここに $[\boldsymbol{T}]^T$ は座標変換マトリックス $[\boldsymbol{T}]$ の転置マトリックス（この場合は逆行列に等しい）である．

立体骨組要素の節点外力ベクトル

節点外力ベクトルについては，ここで紹介するプログラムでは，直接，等価節点外力を入力データファイルから入力することとしており，特別な計算処理は行っていない．また荷重は，全体座標系における水平方向力，鉛直方向力，モーメントとしている．

座標変換マトリックス

これまでに求めた要素剛性方程式は，要素に固定された「要素座標系」でのものである．
実際の構造物では要素は色々な方向を向いているので，統一的に定義した「全体座標系」との関係を求めておく必要がある．

ここで，全体座標系変位を要素座標系変位に変換する座標変換マトリックスを $[\boldsymbol{T}]$ とおくと，
全体座標系での剛性マトリックス $[\boldsymbol{K}]$ と要素座標系での剛性マトリックス $[\boldsymbol{k}]$ には次の関係がある．

\begin{equation}
[\boldsymbol{K}]=[\boldsymbol{T}]^T[\boldsymbol{k}][\boldsymbol{T}]
\end{equation}

この関係は，要素座標系変位 ${\boldsymbol{u}}$ と全体座標系変位 ${\boldsymbol{U}}$ の関係が，以下の通りとなっていることより理解できる．

\begin{equation}
\{\boldsymbol{u}\}=[\boldsymbol{T}]\{\boldsymbol{U}\} \qquad\qquad \{\boldsymbol{u}\}^T=\{\boldsymbol{U}\}^T [\boldsymbol{T}]^T
\end{equation}

ここに $[\boldsymbol{T}]$ は，全体座標系変位・節点力を要素座標系変位・節点力に変換する座標変換マトリックスである．

座標変換マトリックスの具体的表現

部材軸（ｘ軸）の方向余弦

要素を構成する節点を $i$ および $j$ とし，その座標を $(x_i,y_i,z_i)$，$(x_j,y_j,z_j)$とする．

\begin{gather}
l=\cfrac{x_j-x_i}{\ell} \qquad m=\cfrac{y_j-y_i}{\ell} \qquad n=\cfrac{z_j-z_i}{\ell} \\
\ell=\sqrt{(x_j-x_i)^2+(y_j-y_i)^2+(z_j-z_i)^2}
\end{gather}

座標変換小行列

$\boldsymbol{[t]}$ は，全体座標系を部材座標系に変換する座標変換行列 (12 x 12) の小行列 (3 x 3)．

\begin{equation}
\boldsymbol{\{x\}}=\boldsymbol{[t]}\boldsymbol{\{X\}}
\end{equation}

\begin{align}
&\boldsymbol{\{x\}}=\{x,y,z\}^T  & &\text{：部材座標（local coordinate）}\\
&\boldsymbol{\{X\}}＝\{X,Y,Z\}^T & &\text{：全体座標（global coordinate）}
\end{align}

部材長軸（部材座標系x軸）が全体座標系Z軸と平行でない場合

\begin{equation}
\boldsymbol{[t]}=
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos\theta & \sin\theta \\
0 & -\sin\theta & \cos\theta
\end{bmatrix}
\cdot
\begin{bmatrix}
l & m & n \\
-\cfrac{m}{\sqrt{l^2+m^2}} & \cfrac{l}{\sqrt{l^2+m^2}} & 0 \\
-\cfrac{l \cdot n}{\sqrt{l^2+m^2}} & -\cfrac{m \cdot n}{\sqrt{l^2+m^2}} & \sqrt{l^2+m^2}
\end{bmatrix}
\qquad \text{($\theta$ : chord angle)}
\end{equation}

部材長軸（部材座標系x軸）が全体座標系Z軸と平行な場合
この時，$x_j-x_i=0$，$y_j-y_i=0$ であり，$l=0$，$m=0$ であることから，$\boldsymbol{[t]}$を書き直す必要がある．

\begin{equation}
\boldsymbol{[t]}=
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos\theta & \sin\theta \\
0 & -\sin\theta & \cos\theta
\end{bmatrix}
\cdot
\begin{bmatrix}
0 & 0 & n \\
n & 0 & 0 \\
0 & 1 & 0
\end{bmatrix}
\qquad \text{($\theta$ : chord angle)}
\end{equation}

座標変換行列

\begin{equation}
\boldsymbol{[T]}=
\begin{bmatrix}
\boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0} \\\
\boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0} \\\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0} \\\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0}   &\boldsymbol{[t]}
\end{bmatrix}
\end{equation}

コードアングルが０となる部材座標系（x-y-z）と全体座標系（X-Y-Z）の関係

全体剛性方程式

以上より，全体座標系での剛性方程式は，以下のとおり表現される．

\begin{equation}
[\boldsymbol{K}]\{\boldsymbol{U}\}=\{\boldsymbol{F}\}+\{\boldsymbol{F_b}\}+\{\boldsymbol{F_t}\}
\end{equation}

\begin{align}
&\{\boldsymbol{U}\}                          & & \text{全体座標系での節点変位ベクトル}\\
&\{\boldsymbol{F}\}                          & & \text{全体座標系での節点外力ベクトル}\\
&\{\boldsymbol{F_b}\}=\{\boldsymbol{f_b}\}   & & \text{全体座標系での節点慣性力ベクトル}\\
&\{\boldsymbol{F_t}\}=[\boldsymbol{T}]^T\{\boldsymbol{f_t}\}  & & \text{全体座標系での節点温度荷重ベクトル}\\
&[\boldsymbol{K}]=[\boldsymbol{T}]^T[\boldsymbol{k}][\boldsymbol{T}]   & & \text{全体座標系での剛性マトリックス}\\
&[\boldsymbol{k}]                            & & \text{要素座標系での剛性マトリックス}\\
&\{\boldsymbol{f_b}\}                          & & \text{要素座標系での節点慣性力ベクトル}\\
&\{\boldsymbol{f_t}\}                          & & \text{要素座標系での節点温度荷重ベクトル}\\
&[\boldsymbol{T}]                            & & \text{座標変換行列}
\end{align}

上式を解くことにより，全体座標系での節点変位 ${\boldsymbol{U}}$ を求めることができる．

要素断面力

要素断面力は，上記で得られた全体座標系変位 ${\boldsymbol{U}}$ を，座標変換マトリックス $[\boldsymbol{T}]$ により要素座標系変位に変換し，これに要素剛性マトリックス $[\boldsymbol{k}]$ を乗じて算定する．すなわち，要素断面力 ${\boldsymbol{f_{sec}}}$ は，以下により求める．

\begin{equation}
\{\boldsymbol{f_{sec}}\}=[\boldsymbol{k}]\cdot ([\boldsymbol{T}]\{\boldsymbol{U}\})
\end{equation}

上記で計算した計算した軸力は，温度変化による効果が含まれていないため，以下のように軸力の温度効果補正を行う必要が有ることに注意する．

\begin{align}
N_1'=&N_1+EA\cdot\alpha\cdot\Delta T \\
N_2'=&N_2-EA\cdot\alpha\cdot\Delta T
\end{align}

ここに，$N_1'$，$N_2'$ は温度変化を考慮した節点 1 および 2 での軸力，$N_1$，$N_2$ は連立方程式の解としての変位より計算した軸力，
$EA$ は部材の軸剛性，$\alpha$ は部材の線膨張係数，$\Delta T$ は部材の温度変化量である．

立体骨組構造解析プログラム

プログラムの概要

プログラミング言語はPythonとする
立体骨組弾性解析を行う
要素は，１要素２節点，１節点６自由度の梁要素とする
全体座標システムは，水平面平面上の直行座標系右方向を x 軸正方向，上方向を y 軸正方向とする
荷重として，節点外力，節点強制変位（ゼロ変位含む），節点温度変化，要素への慣性力を扱える
外力は，節点に作用するもののみ扱える．要素中間への荷重は扱えない
連立一次方程式は，その大きさ（１節点自由度ｘ総節点数）を変えずに，節点強制変位導入による処理を行い解いている
連立一次方程式は，CSR 変換後，scipy.sparse.linalg.spsolve() により解く
出力は，節点変位および要素断面力とする
入力データファイル，出力データファイルの文字区切りは空白とする
入力データファイル名，出力データファイル名は，コマンドライン引数として入力する

扱える境界条件の説明

項目	説明
節点外力	載荷節点数，積荷節点番号，荷重値を指定
節点温度変化	全節点の温度変化量を指定．温度上昇を正値とする
要素慣性力	材料特性データの中で x / y / z方向の加速度を重力加速度に対する比で指定
節点強制変位	強制変位を与える節点数，節点番号，変位値を指定

プログラムの構成（プログラム名：py_fem_3dfrmx.py）

   1. モジュール読み込み
   2. 関数：データ入力 (def inpdata_3dfrm)
       (1) 基本数値読み込み
       (2) 配列宣言
       (3) 材料特性データ読み込み
       (4) 要素構成節点番号および材料番号読み込み
       (5) 節点座標および節点温度変化量読み込み
       (6) 境界条件読み込み
       (7) 節点荷重データ読み込み
   3. 関数：入力データ書き出し (def prinp_3dfrm)
   4. 関数：計算結果書き出し (def prout_3dfrm)
   5. 関数：要素剛性マトリックス作成 (def sm_3dfrm)
   6. 関数：座標変換マトリックス作成 (def tm_3dfrm)
   7. 関数：節点慣性力ベクトル作成 (def bfvec_3dfrm)
   8. 関数：節点温度荷重ベクトル作成 (def tfvec_3dfrm)
   9. 関数：要素断面力計算 (def calsecf_3dfrm)
  10. 関数：メインルーチン (def main_3dfrm)
       (1) コマンドラインより入出力ファイル名取得
       (2) データ入力
       (3) 入力データ書き出し
       (4) 剛性マトリックス組み立て
       (5) 境界条件処理（変位既知自由度の処理）
       (6) 剛性方程式を解く（変位計算）
       (7) 境界条件処理（既知変位の組み込み）
       (8) 要素断面力計算
       (9) 計算結果書き出し
       (10) 計算時間等情報表示
  11. メインルーチン実行

実行コマンドと入出力データ

解析実行コマンド

python3 py_fem_3dfrmx.py inp.txt out.txt

項目	説明
`py_fem_3dfrmx.py`	Python による FEM プログラム名
`inp.txt`	入力データファイル名（空白区切りデータ）
`out.txt`	出力データファイル名（空白区切りデータ）

入力データファイルフォーマット

npoin  nele  nsec  npfix  nlod                   # Basic values for analysis
E  po  A  Ix  Iy  Iz  theta  alpha  gamma  gkX  gkY  gkZ
    ..... (1 to nsec) .....                      # Material properties
node_1  node_2  isec                             # Element connectivity, material set number
    ..... (1 to nele) .....                      #
x  y  z  deltaT                                  # Node coordinate, temperature change of node
    ..... (1 to npoin) .....                     #
lp  kox  koy  koz  kmx  kmy  kmz  rdis_x  rdis_y  rdis_z  rrot_x  rrot_y  rrot_z
    ..... (1 to npfix) .....                     # Restricted node and restrict conditions
lp  fp_x  fp_y  fp_z  mp_x  mp_y  mp_z           # Loaded node and loading conditions
    ..... (1 to nlod) .....                      #     (omit data input if nlod=0)

項目	説明
`npoin, nele, nsec`	節点数，要素数，断面特性数
`npfix, nlod`	拘束節点数，載荷節点数
`E, po, A`	要素弾性係数，要素ポアソン比，要素断面積
`Ix, Iy, Iz`	ねじり定数，y軸周り断面２次モーメント，z軸周り断面２次モーメント
`theta, alpha`	コードアングル，要素線膨張係数
`gamma, gkX, gkY, gkZ`	単位体積重量，x / y / z方向加速度（gの比）
`node_1, node_2, isec`	節点１，節点２，断面特性番号
`x, y, z, deltaT`	節点x座標，節点y座標，節点z座標，節点温度変化
`lp, kox, koy, koz`	拘束節点番号，x / y / z 方向変位拘束の有無（拘束: 1，自由: 0）
`kmx, kmy, kmz`	x / y / z 軸周り回転拘束の有無（拘束: 1，自由: 0）
`r_disx, r_disy, rdis_z`	x / y / z方向既知変位（無拘束でも0を入力）
`rrot_x, rrot_y, rrot_z`	x / y / z軸周り既知回転量（無拘束でも0を入力）
`lp, fp_x, fp_y, fp_z`	載荷重節点番号，x方向荷重，y方向荷重，z方向荷重
`mp_x, mp_y, mp_z`	x軸周りモーメント（ねじり），y軸周りモーメント，z軸周りモーメント

出力データファイルフォーマット

npoin  nele  nsec npfix  nlod
    (Each value of above)
sec  E      po     A    J    Iy    Iz    theta
sec  alpha  gamma  gkX  gkY  gkZ
    sec   : Material number
    E     : Elastic modulus
    po    : Poisson's ratio
    A     : Section area
    J     ; Torsional constant
    Iy    : Moment of inertia around y-axis
    Iz    : Moment of inertia around z-axis
    theta : Chord angle
    alpha : Thermal expansion coefficient
    gamma : Unit weight
    gkX   : Acceleration in X-direction (ratio to g)
    gkY   : Acceleration in Y-direction (ratio to g)
    gkZ   : Acceleration in Z-direction (ratio to g)
    ..... (1 to nsec) .....
node   x   y   z   fx   fy   fz   mx   my   mz  deltaT
    node   : Node number
    x      : x-coordinate
    y      : y-coordinate
    fx     : Load in X-direction
    fy     : Load in Y-direction
    fz     : Load in Z-direction
    mx     : Moment load around X-axis
    my     : Moment load around Y-axis
    mz     : Moment load around Z-axis
    deltaT : Temperature change of node
    ..... (1 to npoin) .....
node   kox   koy   koz   kmx   kmy   kmz   rdis_x   rdis_y   rdis_z   rrot_x   rrot_y   rrot_z
    node    : Node number
    kox     : Index of restriction in X-direction (0: free, 1: fixed)
    koy     : Index of restriction in Y-direction (0: free, 1: fixed)
    koz     : Index of restriction in Z-direction (0: free, 1: fixed)
    kmx     : Index of restriction in rotation around X-axis  (0: free, 1: fixed)
    kmy     : Index of restriction in rotation around Y-axis  (0: free, 1: fixed)
    kmz     : Index of restriction in rotation around Z-axis  (0: free, 1: fixed)
    rdis_x  : Given displacement in X-direction
    rdis_y  : Given displacement in Y-direction
    rdis_z  : Given displacement in Z-direction
    rrot_x  : Given displacement in rotation around X-axis
    rrot_y  : Given displacement in rotation around Y-axis
    rrot_z  : Given displacement in rotation around Z-axis
    ..... (1 to npfix) .....
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) .....
node   dis-x   dis-y   dis-z   rot-x   rot-y   rot-z
    node  : Node number
    dis-x : Displacement in X-direction
    dis-y : Displacement in Y-direction
    dis-z : Displacement in Z-direction
    rot-x : Rotation around X-axis
    rot-y : Rotation around Y-axis
    rot-z : Rotation around Z-axis
    ..... (1 to npoin) .....
elem nodei   N_i   Sy_i   Sz_i   Mx_i   My_i   Mz_i
elem nodej   N_j   Sy_j   Sz_j   Mx_j   My_j   Mz_j
    elem  : Element number
    nodei : First node nymber
    nodej : Second node number
    N_i   : Axial force of node-i
    Sy_i  : Shear force of node-i in y-direction
    Sz_i  : Shear force of node-i in z-direction
    Mx_i  : Moment of node-i around x-axis (torsional moment)
    My_i  : Moment of node-i around y-axis (bending moment)
    Mz_i  : Moment of node-i around z-axis (bending moment)
    N_j   : Axial force of node-j
    Sy_j  : Shear force of node-j in y-direction
    Sz_j  : Shear force of node-j in z-direction
    Mx_j  : Moment of node-j around x-axis (torsional moment)
    My_j  : Moment of node-j around y-axis (bending moment)
    Mz_j  : Moment of node-j around z-axis (bending moment)
    ..... (1 to nele) .....
n=(total degrees of freedom)  time=(calculation time)

項目	説明
`dis-x, dis-y, dis-z`	x方向変位，y方向変位，z方向変位
`rot-x, rot-y, rot-z`	x軸周り回転量，y軸周り回転量，z軸周り回転量
`N, Sy, Sz`	軸力，y方向せん断力，z方向せん断力
`Mx, My, Mz`	x軸周りモーメント（ねじり），y軸周り曲げモーメント，z軸周り曲げモーメント

プログラミング

プログラム本体を，「プログラムの構成」に沿って以下に示します．

1. モジュール読み込み

numpy：配列演算用
scipy.sparse.linalg.spsolve：疎行列連立一次方程式解法
scipy.sparse import csr_matrix：疎行列圧縮格納
sys：コマンドライン引数使用
time：計算時間計測用

# ======================
# 3D Frame Analysis
# ======================
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time

2. 関数：データ入力 (def inpdata_3dfrm)

基本数値読み込み（npoin, nele, nsec, npfix, nlod）
配列宣言（ae, node, x, deltaT, mpfix, rdis, fp）
材料特性データ読み込み（配列 ae）
要素構成節点番号および材料番号読み込み（配列 node）
節点座標および節点温度変化量読み込み（配列 x, deltaT）
境界条件読み込み（配列 mpfix, rdis）
節点荷重データ読み込み（配列 fp）

def inpdata_3dfrm(fnameR,nod,nfree):
    f=open(fnameR,'r')
    text=f.readline()
    text=text.strip()
    text=text.split()
    npoin=int(text[0]) # Number of nodes
    nele =int(text[1]) # Number of elements
    nsec =int(text[2]) # Number of sections
    npfix=int(text[3]) # Number of restricted nodes
    nlod =int(text[4]) # Number of loaded nodes
    # array declaration
    ae    =np.zeros((12,nsec),dtype=np.float64)      # Section characteristics
    node  =np.zeros((nod+1,nele),dtype=np.int)       # node-element relationship
    x     =np.zeros((3,npoin),dtype=np.float64)      # Coordinates of nodes
    deltaT=np.zeros(npoin,dtype=np.float64)          # Temperature change of nodes
    mpfix =np.zeros((nfree,npoin),dtype=np.int)      # Ristrict conditions
    rdis  =np.zeros((nfree,npoin),dtype=np.float64)  # Ristricted displacement
    fp    =np.zeros(nfree*npoin,dtype=np.float64)    # External force vector
    # section characteristics
    for i in range(0,nsec):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ae[0,i] =float(text[0])  # E     : elastic modulus
        ae[1,i] =float(text[1])  # po    : Poisson's ratio
        ae[2,i] =float(text[2])  # aa    : section area
        ae[3,i] =float(text[3])  # aix   : tortional constant
        ae[4,i] =float(text[4])  # aiy   : moment of inertia aroutd y-axis
        ae[5,i] =float(text[5])  # aiz   : moment of inertia around z-axis
        ae[6,i] =float(text[6])  # theta : chord angle
        ae[7,i] =float(text[7])  # alpha : thermal expansion coefficient
        ae[8,i] =float(text[8])  # gamma : unit weight of material
        ae[9,i] =float(text[9])  # gkX   : acceleration in X-direction
        ae[10,i]=float(text[10]) # gkY   : acceleration in Y-direction
        ae[11,i]=float(text[11]) # gkZ   : acceleration in Z-direction
    # element-node
    for i in range(0,nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        node[0,i]=int(text[0]) #node_1
        node[1,i]=int(text[1]) #node_2
        node[2,i]=int(text[2]) #section characteristic number
    # node coordinates
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        x[0,i]=float(text[0])    # x-coordinate
        x[1,i]=float(text[1])    # y-coordinate
        x[2,i]=float(text[2])    # z-coordinate
        deltaT[i]=float(text[3]) # Temperature change
    # boundary conditions (0:free, 1:restricted)
    for i in range(0,npfix):
        text=f.readline()
        text=text.strip()
        text=text.split()
        lp=int(text[0])              #fixed node
        mpfix[0,lp-1]=int(text[1])   #fixed in x-direction
        mpfix[1,lp-1]=int(text[2])   #fixed in y-direction
        mpfix[2,lp-1]=int(text[3])   #fixed in z-direction
        mpfix[3,lp-1]=int(text[4])   #fixed in rotation around x-axis
        mpfix[4,lp-1]=int(text[5])   #fixed in rotation around y-axis
        mpfix[5,lp-1]=int(text[6])   #fixed in rotation around z-axis
        rdis[0,lp-1]=float(text[7])  #fixed displacement in x-direction
        rdis[1,lp-1]=float(text[8])  #fixed displacement in y-direction
        rdis[2,lp-1]=float(text[9])  #fixed displacement in z-direction
        rdis[3,lp-1]=float(text[10]) #fixed rotation around x-axis
        rdis[4,lp-1]=float(text[11]) #fixed rotation around y-axis
        rdis[5,lp-1]=float(text[12]) #fixed rotation around z-axis
    # load
    if 0&lt;nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[6*lp-6]=float(text[1]) #load in x-direction
            fp[6*lp-5]=float(text[2]) #load in y-direction
            fp[6*lp-4]=float(text[3]) #load in z-direction
            fp[6*lp-3]=float(text[4]) #moment around x-axis
            fp[6*lp-2]=float(text[5]) #moment around y-axis
            fp[6*lp-1]=float(text[6]) #moment around z-axis
    f.close()
    return npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp

3. 関数：入力データ書き出し (def prinp_3dfrm)

def prinp_3dfrm(fnameW,npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s} {4:&gt;5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout)
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
    .format('sec','E','po','A','J','Iy','Iz','theta'),file=fout)
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s}'
    .format('sec','alpha','gamma','gkX','gkY','gkZ'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i],ae[6,i]),file=fout)
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e}'
        .format(i+1,ae[7,i],ae[8,i],ae[9,i],ae[10,i],ae[11,i]),file=fout)
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s} {8:&gt;15s} {9:&gt;15s} {10:&gt;15s}'
    .format('node','x','y','z','fx','fy','fz','mx','my','mz','deltaT'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e} {8:15.7e} {9:15.7e} {10:15.7e}'
        .format(lp,x[0,i],x[1,i],x[2,i],fp[6*lp-6],fp[6*lp-5],fp[6*lp-4],fp[6*lp-3],fp[6*lp-2],fp[6*lp-1],deltaT[i]),file=fout)
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s} {4:&gt;5s} {5:&gt;5s} {6:&gt;5s} {7:&gt;15s} {8:&gt;15s} {9:&gt;15s} {10:&gt;15s} {11:&gt;15s} {12:&gt;15s}'
    .format('node','kox','koy','koz','kmx','kmy','kmz','rdis_x','rdis_y','rdis_z','rrot_x','rrot_y','rrot_z'),file=fout)
    for i in range(0,npoin):
        if 0&lt;mpfix[0,i]+mpfix[1,i]+mpfix[2,i]+mpfix[3,i]+mpfix[4,i]+mpfix[5,i]:
            lp=i+1
            print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d} {6:5d} {7:15.7e} {8:15.7e} {9:15.7e} {10:15.7e} {11:15.7e} {12:15.7e}'
            .format(lp,mpfix[0,i],mpfix[1,i],mpfix[2,i],mpfix[3,i],mpfix[4,i],mpfix[5,i],rdis[0,i],rdis[1,i],rdis[2,i],rdis[3,i],rdis[4,i],rdis[5,i]),file=fout)
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s}'.format('elem','i','j','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout)
    fout.close()

4. 関数：計算結果書き出し (def prout_3dfrm)

def prout_3dfrm(fnameW,npoin,nele,node,disg,fsec):
    fout=open(fnameW,'a')
    # displacement
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s}'
          .format('node','dis-x','dis-y','dis-z','rot-x','rot-y','rot-z'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
              .format(lp,disg[6*lp-6],disg[6*lp-5],disg[6*lp-4],disg[6*lp-3],disg[6*lp-2],disg[6*lp-1]),file=fout)
    # section force
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
    .format('elem','nodei','N_i','Sy_i','Sz_i','Mx_i','My_i','Mz_i'),file=fout)
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
    .format('elem','nodej','N_j','Sy_j','Sz_j','Mx_j','My_j','Mz_j'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
        .format(ne+1,node[0,ne],fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
        print('{0:5d} {1:5d} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
        .format(ne+1,node[1,ne],fsec[6,ne],fsec[7,ne],fsec[8,ne],fsec[9,ne],fsec[10,ne],fsec[11,ne]),file=fout)
    fout.close()

5. 関数：要素剛性マトリックス作成 (def sm_3dfrm)

\begin{equation}
\boldsymbol{[k]}=
\begin{bmatrix}
\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} & 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} \\
0 & 0 & \cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 & 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} \\
-\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} & 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} \\
0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 & 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} \\
\end{bmatrix}
\end{equation}

def sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2):
    ek=np.zeros((12,12),dtype=np.float64) # local stiffness matrix
    xx=x2-x1
    yy=y2-y1
    zz=z2-z1
    el=np.sqrt(xx**2+yy**2+zz**2)
    ek[ 0, 0]= EA/el
    ek[ 0, 6]=-EA/el
    ek[ 1, 1]= 12*EIz/el**3
    ek[ 1, 5]=  6*EIz/el**2
    ek[ 1, 7]=-12*EIz/el**3
    ek[ 1,11]=  6*EIz/el**2
    ek[ 2, 2]= 12*EIy/el**3
    ek[ 2, 4]= -6*EIy/el**2
    ek[ 2, 8]=-12*EIy/el**3
    ek[ 2,10]= -6*EIy/el**2
    ek[ 3, 3]= GJ/el
    ek[ 3, 9]=-GJ/el
    ek[ 4, 2]= -6*EIy/el**2
    ek[ 4, 4]=  4*EIy/el
    ek[ 4, 8]=  6*EIy/el**2
    ek[ 4,10]=  2*EIy/el
    ek[ 5, 1]=  6*EIz/el**2
    ek[ 5, 5]=  4*EIz/el
    ek[ 5, 7]= -6*EIz/el**2
    ek[ 5,11]=  2*EIz/el
    ek[ 6, 0]=-EA/el
    ek[ 6, 6]= EA/el
    ek[ 7, 1]=-12*EIz/el**3
    ek[ 7, 5]= -6*EIz/el**2
    ek[ 7, 7]= 12*EIz/el**3
    ek[ 7,11]= -6*EIz/el**2
    ek[ 8, 2]=-12*EIy/el**3
    ek[ 8, 4]=  6*EIy/el**2
    ek[ 8, 8]= 12*EIy/el**3
    ek[ 8,10]=  6*EIy/el**2
    ek[ 9, 3]=-GJ/el
    ek[ 9, 9]= GJ/el
    ek[10, 2]= -6*EIy/el**2
    ek[10, 4]=  2*EIy/el
    ek[10, 8]=  6*EIy/el**2
    ek[10,10]=  4*EIy/el
    ek[11, 1]=  6*EIz/el**2
    ek[11, 5]=  2*EIz/el
    ek[11, 7]= -6*EIz/el**2
    ek[11,11]=  4*EIz/el
    return ek

6. 関数：座標変換マトリックス作成 (def tm_3dfrm)

\begin{gather}
\ell=\sqrt{(x_j-x_i)^2+(y_j-y_i)^2+(z_j-z_i)^2}\\
l=\cfrac{x_j-x_i}{\ell} \qquad m=\cfrac{y_j-y_i}{\ell} \qquad n=\cfrac{z_j-z_i}{\ell}
\end{gather}

\begin{equation}
\boldsymbol{[t_1]}=
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos\theta & \sin\theta \\
0 & -\sin\theta & \cos\theta
\end{bmatrix}
\qquad \text{($\theta$ : chord angle)}
\end{equation}

\begin{align}
&\boldsymbol{[t_2]}=
\begin{bmatrix}
0 & 0 & n \\
n & 0 & 0 \\
0 & 1 & 0
\end{bmatrix}
& &\text{部材長軸（部材座標系x軸）が全体座標系Z軸と平行な場合}\\
&\boldsymbol{[t_2]}=
\begin{bmatrix}
l & m & n \\
-\cfrac{m}{\sqrt{l^2+m^2}} & \cfrac{l}{\sqrt{l^2+m^2}} & 0 \\
-\cfrac{l \cdot n}{\sqrt{l^2+m^2}} & -\cfrac{m \cdot n}{\sqrt{l^2+m^2}} & \sqrt{l^2+m^2}
\end{bmatrix}
& &\text{部材長軸（部材座標系x軸）が全体座標系Z軸と平行でない場合}
\end{align}

\begin{equation}
\boldsymbol{[t]}=\boldsymbol{[t_1]}\cdot\boldsymbol{[t_2]}
\end{equation}

\begin{equation}
\boldsymbol{[T]}=
\begin{bmatrix}
\boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0} \\
\boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0} \\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0} \\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0}   &\boldsymbol{[t]}
\end{bmatrix}
\end{equation}

def tm_3dfrm(theta,x1,y1,z1,x2,y2,z2):
    tt=np.zeros((12,12),dtype=np.float64) # transformation matrix
    t1=np.zeros((3,3),dtype=np.float64)
    t2=np.zeros((3,3),dtype=np.float64)
    xx=x2-x1
    yy=y2-y1
    zz=z2-z1
    el=np.sqrt(xx**2+yy**2+zz**2)
    t1[0,0]=1
    t1[1,1]= np.cos(theta)
    t1[1,2]= np.sin(theta)
    t1[2,1]=-np.sin(theta)
    t1[2,2]= np.cos(theta)
    ll=(x2-x1)/el
    mm=(y2-y1)/el
    nn=(z2-z1)/el
    if x2-x1==0.0 and y2-y1==0.0:
        t2[0,2]=nn
        t2[1,0]=nn
        t2[2,1]=1.0
    else:
        qq=np.sqrt(ll**2+mm**2)
        t2[0,0]=ll
        t2[0,1]=mm
        t2[0,2]=nn
        t2[1,0]=-mm/qq
        t2[1,1]= ll/qq
        t2[2,0]=-ll*nn/qq
        t2[2,1]=-mm*nn/qq
        t2[2,2]=qq
    t3=np.dot(t1,t2)
    tt[0:3,0:3]  =t3[0:3,0:3]
    tt[3:6,3:6]  =t3[0:3,0:3]
    tt[6:9,6:9]  =t3[0:3,0:3]
    tt[9:12,9:12]=t3[0:3,0:3]
    return tt

7. 関数：節点慣性力ベクトル作成 (def bfvec_3dfrm)

\begin{equation*}
\{\boldsymbol{F_{b}}\}=\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot A \cdot L}{2}
\begin{Bmatrix}
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0   \\
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0
\end{Bmatrix}
\end{equation*}

def bfvec_3dfrm(A,gamma,gkX,gkY,gkZ,x1,y1,z1,x2,y2,z2):
    # Body force vector in global coordinate system
    bfe_g=np.zeros(12,dtype=np.float64)
    xx=x2-x1
    yy=y2-y1
    zz=z2-z1
    el=np.sqrt(xx**2+yy**2+zz**2)
    bfe_g[0]=0.5*gamma*A*el*gkX
    bfe_g[1]=0.5*gamma*A*el*gkY
    bfe_g[2]=0.5*gamma*A*el*gkZ
    bfe_g[6]=0.5*gamma*A*el*gkX
    bfe_g[7]=0.5*gamma*A*el*gkY
    bfe_g[8]=0.5*gamma*A*el*gkZ
    return bfe_g

8. 関数：節点温度荷重ベクトル作成 (def tfvec_3dfrm)

\begin{equation*}
\{\boldsymbol{f_{t}}\}=EA\cdot\alpha\cdot\Delta T
\begin{Bmatrix}
-1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0
\end{Bmatrix}
\end{equation*}

def tfvec_3dfrm(EA,alpha,tem):
    # Thermal load vector  in local coordinate system
    tfe_l=np.zeros(12,dtype=np.float64)
    tfe_l[0]=-EA*alpha*tem
    tfe_l[6]= EA*alpha*tem
    return tfe_l

9. 関数：要素断面力計算 (def calsecf_3dfrm)

$N_1'$ および $N_2'$ は，温度変化による軸力補正である．

\begin{align}
&\{\boldsymbol{f_{sec}}\}=[\boldsymbol{k}]\cdot ([\boldsymbol{T}]\{\boldsymbol{U}\})\\
&N_1'=N_1+EA\cdot\alpha\cdot\Delta T \\
&N_2'=N_2-EA\cdot\alpha\cdot\Delta T
\end{align}

def calsecf_3dfrm(EA,GJ,EIy,EIz,theta,alpha,tem,dis,x1,y1,z1,x2,y2,z2):
    ek=sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2) # Stiffness matrix in local coordinate
    tt=tm_3dfrm(theta,x1,y1,z1,x2,y2,z2) # Transformation matrix
    secf=np.dot(ek,np.dot(tt,dis))
    secf[0]=secf[0]+EA*alpha*tem
    secf[6]=secf[6]-EA*alpha*tem
    return secf

10. 関数：メインルーチン

メインルーチン処理手順

コマンドラインより入出力ファイル名取得：
コマンドライン引数として入出力データファイル名を取得．
入力データファイル名：fnameR，出力データファイル名：fnameW
データ入力：
関数 inpdata_3dfrm により入力データ取得
入力データ書き出し：
関数 prinp_3dfrm により入力データ書き出し
剛性マトリックス組み立て：
各要素毎に，剛性マトリックス ck （座標変換後），温度荷重ベクトル tfe （座標変換後），慣性力ベクトル bfe を計算し，全体剛性方程式に組み込んでいく．
ここで剛性マトリックス，温度荷重ベクトル，慣性力ベクトル計算の関数の引数は，全てスカラーとした．
要素番号が指定されれば，要素を構成する節点座標や要素物性値は特定できるため，他の要素の値を含んで収納する配列よりは，スカラーを送り込んだほうがプログラム的にスッキリするように思う．
全体剛性マトリックスは変数 gk に，荷重ベクトルは変数 fp に加算されていく．
要素剛性マトリックスを全体剛性マトリックスに組み込むための位置情報は，配列 ir[] に格納されている．
境界条件処理（変位既知自由度の処理）：
全体剛性マトリックスの大きさは変えず，変位既知節点の処理を行う．

項目	説明
`npoin`	総節点数
`nfree`	１節点自由度（ここでは 6）
`mpfix`	各節点毎の拘束条件を格納した配列（=1：拘束，=0：自由）
`fp`	節点荷重ベクトル
`rdis`	変位拘束する節点の変位量
`gk`	全体剛性マトリックス
`disg`	全体剛性方程式の解（変位）

剛性方程式を解く（変位計算）：
全体剛性マトリックス gk を CSR 変換後，scipy.sparse.linalg.spsolve() により連立一次方程式を解く．
境界条件処理（既知変位の組み込み）：
連立一次方程式を解いた後，既知変位に相当する自由度の項に，既知変位を入れ直すことを忘れないようにする．
要素断面力計算：
関数 calsecf_3dfrm による要素毎の断面力計算実施
計算結果書き出し：
関数 prout_3dfrm による計算結果出力
計算時間等情報表示：
以下により，計算完了時刻とスタート時刻の差分を処理時間として出力する．

import time

start=time.time()
......
dtime=time.time()-start

メインルーチンプログラム

def main_3dfrm():
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=2              # Number of nodes per element
    nfree=6            # Degree of freedom per node
    # data input
    npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_3dfrm(fnameR,nod,nfree)
    # print out of input data
    prinp_3dfrm(fnameW,npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp)
    # array declaration
    ir=np.zeros(nod*nfree,dtype=np.int)   # Work vector for matrix assembly
    gk=np.zeros((nfree*npoin,nfree*npoin),dtype=np.float64)   # Global stiffness matrix
    # assembly of stiffness matrix & load vectors
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]; y1=x[1,i]; z1=x[2,i]
        x2=x[0,j]; y2=x[1,j]; z2=x[2,j]
        ee   =ae[0,m]  # elastic modulus
        po   =ae[1,m]  # Poisson's ratio
        aa   =ae[2,m]  # section area
        aix  =ae[3,m] # tortional constant
        aiy  =ae[4,m] # moment of inertia around y-axis
        aiz  =ae[5,m] # moment of inertia around z-axis
        theta=np.radians(ae[6,m]) # chord angle
        alpha=ae[7,m] # thermal expansion coefficient
        gamma=ae[8,m] # unit weight of material
        gkX  =ae[9,m]   # seismic coefficient in X-direction
        gkY  =ae[10,m]  # seismic coefficient in Y-direction
        gkZ  =ae[11,m]  # seismic coefficient in Z-direction
        A=aa  # section area
        EA=ee*aa
        GJ=ee/2/(1+po)*aix
        EIy=ee*aiy
        EIz=ee*aiz
        tem=0.5*(deltaT[i]+deltaT[j])                # average temperature change
        tt   =tm_3dfrm(theta,x1,y1,z1,x2,y2,z2)         # Transformation matrix
        ek   =sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2) # Stiffness matrix in local coordinate
        ck   =np.dot(np.dot(tt.T,ek),tt)             # Stiffness matrix in global coordinate
        tfe_l=tfvec_3dfrm(EA,alpha,tem)              # Thermal load vector in local coordinate
        tfe  =np.dot(tt.T,tfe_l)                     # Thermal load vector in global coordinate
        bfe  =bfvec_3dfrm(A,gamma,gkX,gkY,gkZ,x1,y1,z1,x2,y2,z2) # Body force vector in global coordinate
        ir[11]=6*j+5; ir[10]=ir[11]-1; ir[9]=ir[10]-1; ir[8]=ir[9]-1; ir[7]=ir[8]-1; ir[6]=ir[7]-1
        ir[5] =6*i+5; ir[4] =ir[5]-1 ; ir[3]=ir[4]-1 ; ir[2]=ir[3]-1; ir[1]=ir[2]-1; ir[0]=ir[1]-1
        for i in range(0,nod*nfree):
            it=ir[i]
            fp[it]=fp[it]+tfe[i]+bfe[i]
            for j in range(0,nod*nfree):
                jt=ir[j]
                gk[it,jt]=gk[it,jt]+ck[i,j]
    # treatment of boundary conditions
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp[iz]=0.0
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp=fp-rdis[j,i]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0
    # solution of simultaneous linear equations
    #disg = np.linalg.solve(gk, fp)
    gk = csr_matrix(gk)
    disg = spsolve(gk, fp, use_umfpack=True)
    # recovery of restricted displacements
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                disg[iz]=rdis[j,i]
    # calculation of section force
    dis=np.zeros(12,dtype=np.float64)
    fsec  =np.zeros((nod*nfree,nele),dtype=np.float64)  # Section force vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]; y1=x[1,i]; z1=x[2,i]
        x2=x[0,j]; y2=x[1,j]; z2=x[2,j]
        ee   =ae[0,m]  # elastic modulus
        po   =ae[1,m]  # Poisson's ratio
        aa   =ae[2,m]  # section area
        aix  =ae[3,m] # tortional constant
        aiy  =ae[4,m] # moment of inertia around y-axis
        aiz  =ae[5,m] # moment of inertia around z-axis
        theta=np.radians(ae[6,m]) # chord angle
        alpha=ae[7,m] # thermal expansion coefficient
        EA=ee*aa
        GJ=ee/2/(1+po)*aix
        EIy=ee*aiy
        EIz=ee*aiz
        tem=0.5*(deltaT[i]+deltaT[j]) # average temperature change
        dis[0]=disg[6*i]  ; dis[1] =disg[6*i+1]; dis[2]= disg[6*i+2]
        dis[3]=disg[6*i+3]; dis[4] =disg[6*i+4]; dis[5]= disg[6*i+5]
        dis[6]=disg[6*j]  ; dis[7] =disg[6*j+1]; dis[8]= disg[6*j+2]
        dis[9]=disg[6*j+3]; dis[10]=disg[6*j+4]; dis[11]=disg[6*j+5]
        fsec[:,ne]=calsecf_3dfrm(EA,GJ,EIy,EIz,theta,alpha,tem,dis,x1,y1,z1,x2,y2,z2)
    # print out of result
    prout_3dfrm(fnameW,npoin,nele,node,disg,fsec)
    # information
    dtime=time.time()-start
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
    fout=open(fnameW,'a')
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
    fout.close()

11. メインルーチン実行

#==============
# Execution
#==============
if __name__ == '__main__': main_3dfrm()

以上

Pythonで立体骨組構造解析