大きな修正
- (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<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:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>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:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'
.format('sec','E','po','A','J','Iy','Iz','theta'),file=fout)
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>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:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s} {8:>15s} {9:>15s} {10:>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:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s} {6:>5s} {7:>15s} {8:>15s} {9:>15s} {10:>15s} {11:>15s} {12:>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<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:>5s} {1:>5s} {2:>5s} {3:>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:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>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:>5s} {1:>5s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'
.format('elem','nodei','N_i','Sy_i','Sz_i','Mx_i','My_i','Mz_i'),file=fout)
print('{0:>5s} {1:>5s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>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()
以上