- 2017.12.27 投稿
- 2018.08.05 修正（プログラム改良）
（情報）初期のプログラムから以下の点を変更しています．
- 関数名を小文字にした
- データ入力部を関数にした
- 計算用関数（データ入力・プリントアウト部を除く）の引数を全てスカラーとした
- 無意味な配列初期化の部分をカットした
- 配列初期化時の形状をタプルで指定するようにした（何故かこれまでリストで渡していた）
- メインルーチンも関数化した．これによりグローバル変数を関数の引数で読み込むなどのおかしな点はなくなったと思います．
はじめに
２次元応力解析にひき続き，平面骨組弾性解析プログラムの解説を作ってみました．
基本となる剛性方程式，節点の変位拘束の処理方法などは２次元応力解析と同じですが，座標変換マトリックスを導入する必要があるところが大きな違いです．
基本的な事項は，以下の記事にまとめてあります．
梁要素の場合，曲げに関する剛性マトリックスの計算は４ｘ４の計算となり，めんどうなのでSympyによる計算チェックを行ってみました．記号のまま逆行列を求めたり，微分・積分が実行でき，なかなか便利です（最後の章参照）．
理論展開およびプログラムには，以下の条件が入っています．
- 線形弾性材料の微小変形を仮定している
- 物体力として慣性力を考慮する
- 初期ひずみとして温度変化によるひずみを考慮する
このプログラムでは，モデルが大規模になることは無い（自由度1000程度以下）との想定で，疎行列処理は導入していません．
このシリーズの投稿のリンク
- Python：FEMでの疎行列計算利用による高速化（３次元骨組構造解析での事例）
- Pythonによる２次元応力解析プログラム（四角形要素）の理論を少し詳しく解説する（改訂版）
- Pythonによる平面骨組構造解析プログラムの理論を少し詳しく解説する（改訂版）
- Pythonによる２次元飽和ー不飽和定常浸透流解析プログラム（三角形要素）の理論を少し詳しく解説する（改訂版）
- Pythonによる２次元飽和ー不飽和定常浸透流解析プログラム（四角形要素）の理論を少し詳しく解説する
平面骨組構造解析の理論
要素剛性方程式
要素剛性方程式の一般形は，２次元応力解析におけるものと同一である．
変位を求めるための方程式
\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}
軸力と曲げを受ける棒部材の剛性マトリックス
基本事項
平面骨組構造解析のための棒要素を考える．
棒の長手軸方向を$x$軸とし右向きを正，長手軸に直行する方向を$y$軸とし上向きを正とする．
ここで，せん断ひずみを無視すれば，発生ひずみは骨組みの軸方向ひずみ $\epsilon_x$ のみであるが，曲げを受ける場合，たわみ $v$ により $y$ 方向に $\epsilon_x$ が線形分布することになる．
この関係をひずみ-変位関係で表すと以下の通りとなる．
\begin{equation}
\epsilon_x=\frac{du}{dx}-y\frac{d^2 v}{dx^2}
\end{equation}
ここに $u$ は $x$ 方向変位，$v$ は $y$ 方向変位（たわみ）であり，$y=0$ の位置は棒部材の中立軸とする．
また，棒要素の断面は一定であり，荷重は部材端部の節点のみに作用するものとする．
ここでは，線形理論であるため，軸方向変形と曲げ変形を分離してそれぞれの剛性マトリックスを求め，あとで重ね合わせることを考える．
棒要素では，応力も長手方向応力 $\sigma_x$ のみであるため，部材の弾性係数を $E$ として
\begin{equation*}
\{\boldsymbol{\sigma}\}=\sigma_x \qquad [\boldsymbol{D_e}]=E \qquad \{\boldsymbol{\epsilon}\}=\epsilon_x
\end{equation*}
であり，軸方向変形および曲げ変形に関する剛性方程式の形式は以下の通りとなる．
\begin{align}
[\boldsymbol{k_t}]\{\boldsymbol{u_{tnd}}\}=&\{\boldsymbol{f_{tnd}}\} &\text{軸方向の剛性方程式}\\
&\{\boldsymbol{u_{tnd}}\}=\begin{Bmatrix} u_i & u_j \end{Bmatrix}^T &\text{（節点の軸方向変位）}\\
&\{\boldsymbol{f_{tnd}}\}=\begin{Bmatrix} N_i & N_j \end{Bmatrix}^T &\text{（節点の軸力）}\\
[\boldsymbol{k_b}]\{\boldsymbol{u_{bnd}}\}=&\{\boldsymbol{f_{bnd}}\} &\text{曲げに関する剛性方程式}\\
&\{\boldsymbol{u_{tnd}}\}=\begin{Bmatrix} v_i & \theta_i & v_j & \theta_j \end{Bmatrix}^T &\text{（節点の軸直交方向変位・回転角）}\\
&\{\boldsymbol{f_{tnd}}\}=\begin{Bmatrix} S_i & M_i & S_j & M_j \end{Bmatrix}^T &\text{（節点のせん断力・曲げモーメント）}\\
[\boldsymbol{k}]\{\boldsymbol{u_{nd}}\}=&\{\boldsymbol{f_{nd}}\} &\text{軸力と曲げを受ける部材の剛性方程式}\\
&\{\boldsymbol{u_{tnd}}\}=\begin{Bmatrix} u_i & v_i & \theta_i & u_j & v_j & \theta_j \end{Bmatrix}^T &\text{（節点変位）}\\
&\{\boldsymbol{f_{tnd}}\}=\begin{Bmatrix} N_i & S_i & M_i & N_j & S_j & M_j \end{Bmatrix}^T &\text{（節点力）}
\end{align}
軸方向の剛性マトリックス
曲げを考えず軸の伸びのみを考える場合，棒要素は１節点１自由度，２節点のため２自由度を有する．
このため，２個の未定係数 $\alpha_{1\sim 2}$ を用い，以下の変位を仮定する．
\begin{equation}
u=\alpha_1+\alpha_2 x
=\begin{bmatrix} 1 & x \end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \end{Bmatrix}
\end{equation}
ここで，棒の節点を $i$ および $j$ とし，$x$軸の原点を$i$端として，それぞれの節点変位を棒の長さ $L$ と ${\boldsymbol{\alpha}}$ で表せば以下の通りとなる．
\begin{equation}
\begin{Bmatrix} u_i \\ u_j \end{Bmatrix}
=\begin{bmatrix} 1 & 0 \\ 1 & L \end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \end{Bmatrix}
\end{equation}
これを ${\boldsymbol{\alpha}}$ について解けば
\begin{equation}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \end{Bmatrix}
=\begin{bmatrix} 1 & 0 \\ -1/L & 1/L \end{bmatrix}
\begin{Bmatrix} u_i \\ u_j \end{Bmatrix}
\end{equation}
よって要素任意点の変位 $u$ および要素任意点の軸方向一様ひずみ $\epsilon_t$ は以下の通り表される．
\begin{equation}
u=
\begin{bmatrix} 1 & x \end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \end{Bmatrix}
=\begin{bmatrix} 1 & x \end{bmatrix}
\begin{bmatrix} 1 & 0 \\ -1/L & 1/L \end{bmatrix}
\begin{Bmatrix} u_i \\ u_j \end{Bmatrix}
=\begin{bmatrix} 1-\cfrac{1}{L}x & \cfrac{1}{L}x \end{bmatrix}
\begin{Bmatrix} u_i \\ u_j \end{Bmatrix}
\end{equation}
上式より明らかなように，内挿関数 $[\boldsymbol{N_t}]$ は以下の通り定義できる．
\begin{equation}
u=[\boldsymbol{N_t}]\{\boldsymbol{u_{tnd}}\}
\qquad\qquad
[\boldsymbol{N_t}]=
\begin{bmatrix} 1-\cfrac{1}{L}x & \cfrac{1}{L}x \end{bmatrix}
\end{equation}
よって，軸ひずみは，
\begin{equation}
\epsilon_t=\frac{du}{dx}
=\begin{bmatrix}\cfrac{d N_{t1}}{d x} & \cfrac{d N_{t2}}{d x}\end{bmatrix}\{\boldsymbol{u_{tnd}}\}
=\begin{bmatrix} -\cfrac{1}{L} & \cfrac{1}{L} \end{bmatrix}
\begin{Bmatrix} u_i \\ u_j \end{Bmatrix}
\end{equation}
上式から明らかなように，ひずみー変位関係マトリックス $[\boldsymbol{B_t}]$ は
\begin{equation}
\epsilon_t=[\boldsymbol{B_t}]\{\boldsymbol{u_{tnd}}\}
\qquad\qquad
[\boldsymbol{B_t}]=\begin{bmatrix} -\cfrac{1}{L} & \cfrac{1}{L} \end{bmatrix}
\end{equation}
以上により，軸剛性マトリックス $[\boldsymbol{k_t}]$ は以下の通り定義できる．
\begin{align}
[\boldsymbol{k_t}]=&E \cdot \int_V [\boldsymbol{B_t}]^T [\boldsymbol{B_t}] dV=E \cdot \int_A dA \times \int_0^{L}[\boldsymbol{B_t}]^T [\boldsymbol{B_t}] dx \\
=&\frac{EA}{L}
\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \qquad\qquad A=\int_A dA
\end{align}
ここに，$E$ は棒の弾性係数，$A$ は棒の断面積である．
曲げに対する剛性マトリックス
軸の伸びを考えず曲げのみを考える場合，棒要素は「たわみ」と「たわみ角=回転角」の２自由度を有する節点が２つ存在するため，４自由度を有する．
このため，４個の未定係数 $\alpha_{1\sim 4}$ を用い，以下の変位を仮定する．
\begin{align}
&v=\alpha_1+\alpha_2 x+\alpha_3 x^2+\alpha_4 x^3 \qquad \text{（たわみ）}\\
&\theta=\frac{dv}{dx}=\alpha_2+2\alpha_3 x+3\alpha_4 x^2 \qquad \text{（回転角）}
\end{align}
上式を行列表示すれば，
\begin{equation}
\begin{Bmatrix} v \\ \theta \end{Bmatrix}
=\begin{bmatrix} 1 & x & x^2 & x^3 \\ 0 & 1 & 2x & 3x^2 \end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4\end{Bmatrix}
\end{equation}
ここで，棒の節点を $i$ および $j$ としてそれぞれの節点変位を棒の長さ $L$ と ${\boldsymbol{\alpha}}$ で表せば以下の通りとなる．
なお $x$ 軸は部材軸方向右向きが正（原点は $i$ 端），$y$ 軸は部材軸直行方向上向きが正（原点は部材中立軸），回転角 $\theta$ はいずれの節点でも反時計回りを正とする．
\begin{equation}
\begin{Bmatrix} v_i \\ \theta_i \\ v_j \\ \theta_j \end{Bmatrix}
=\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
1 & L & L^2 & L^3 \\
0 & 1 & 2L & 3L^2
\end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{Bmatrix}
\end{equation}
これを ${\boldsymbol{\alpha}}$ について解けば
\begin{equation}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{Bmatrix}
=\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
-3/L^2 & -2/L & 3/L^2 & -1/L \\
2/L^3 & 1/L^2 & -2/L^3 & 1/L^2
\end{bmatrix}
\begin{Bmatrix} v_i \\ \theta_i \\ v_j \\ \theta_j \end{Bmatrix}
\end{equation}
よって要素任意点の変位$v$は以下の通り表される．
\begin{equation}
v=
\begin{bmatrix} 1 & x & x^2 & x^3 \end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{Bmatrix}
=[\boldsymbol{N_b}]\{\boldsymbol{u_{bnd}}\}
\end{equation}
\begin{equation}
[\boldsymbol{N_b}]
=\left[\begin{matrix}1 - \frac{3 x^{2}}{L^{2}} + \frac{2 x^{3}}{L^{3}} & x - \frac{2 x^{2}}{L} + \frac{x^{3}}{L^{2}} & \frac{3 x^{2}}{L^{2}} - \frac{2 x^{3}}{L^{3}} & - \frac{x^{2}}{L} + \frac{x^{3}}{L^{2}}\end{matrix}\right]
\end{equation}
ここで，曲げひずみ $\epsilon_b$ は
\begin{equation}
\epsilon_b=-y\frac{d^2 v}{dx^2}
\end{equation}
で表せるため，これを行列表示すると以下の通りとなる．
\begin{equation}
\epsilon_b=[\boldsymbol{B_b}]\{\boldsymbol{u_{bnd}}\}
\end{equation}
\begin{equation}
[\boldsymbol{B_b}]=-y \left[\cfrac{d^2 \boldsymbol{N_b}}{dx^2}\right]
=\left[\begin{matrix}- y \left(- \frac{6}{L^{2}} + \frac{12 x}{L^{3}}\right) & - y \left(- \frac{4}{L} + \frac{6 x}{L^{2}}\right) & - y \left(\frac{6}{L^{2}} - \frac{12 x}{L^{3}}\right) & - y \left(- \frac{2}{L} + \frac{6 x}{L^{2}}\right)\end{matrix}\right]
\end{equation}
よって曲げ剛性行列 $[\boldsymbol{k_b}]$ は以下の通り定義できる．
\begin{align}
[\boldsymbol{k_b}]=&E \cdot \int_V [\boldsymbol{B_b}]^T[\boldsymbol{B_b}] dV
=E \cdot \int_A \int_0^L [\boldsymbol{B_b}]^T[\boldsymbol{B_b}] dx dA \\
=&E \cdot \int_A
\left[\begin{matrix}\frac{12 y^{2}}{L^{3}} & \frac{6 y^{2}}{L^{2}} & - \frac{12 y^{2}}{L^{3}} & \frac{6 y^{2}}{L^{2}}\\ \frac{6 y^{2}}{L^{2}} & \frac{4 y^{2}}{L} & - \frac{6 y^{2}}{L^{2}} & \frac{2 y^{2}}{L}\\- \frac{12 y^{2}}{L^{3}} & - \frac{6 y^{2}}{L^{2}} & \frac{12 y^{2}}{L^{3}} & - \frac{6 y^{2}}{L^{2}}\\ \frac{6 y^{2}}{L^{2}} & \frac{2 y^{2}}{L} & - \frac{6 y^{2}}{L^{2}} & \frac{4 y^{2}}{L}\end{matrix}\right] dA
=E \cdot \int_A y^2 dA \cdot
\left[\begin{matrix}\frac{12}{L^{3}} & \frac{6}{L^{2}} & - \frac{12}{L^{3}} & \frac{6}{L^{2}}\\ \frac{6}{L^{2}} & \frac{4}{L} & - \frac{6}{L^{2}} & \frac{2}{L}\\ - \frac{12}{L^{3}} & - \frac{6}{L^{2}} & \frac{12}{L^{3}} & - \frac{6}{L^{2}}\\ \frac{6}{L^{2}} & \frac{2}{L} & - \frac{6}{L^{2}} & \frac{4}{L}\end{matrix}\right] \\
=&EI \cdot
\left[\begin{matrix}\frac{12}{L^{3}} & \frac{6}{L^{2}} & - \frac{12}{L^{3}} & \frac{6}{L^{2}}\\ \frac{6}{L^{2}} & \frac{4}{L} & - \frac{6}{L^{2}} & \frac{2}{L}\\ - \frac{12}{L^{3}} & - \frac{6}{L^{2}} & \frac{12}{L^{3}} & - \frac{6}{L^{2}}\\ \frac{6}{L^{2}} & \frac{2}{L} & - \frac{6}{L^{2}} & \frac{4}{L}\end{matrix}\right]
\qquad\qquad
I=\int_A y^2 dA
\end{align}
ここに，$E$ は棒の弾性係数，$I$ は棒の断面２次モーメントである．
軸力と曲げを受ける平面骨組要素の剛性マトリックス
以上により，軸剛性行列と曲げ剛性行列を，それぞれの自由度の行・列で足し合わせることにより，軸力と曲げを受ける棒要素の要素剛性方程式が以下の通り得られる．
\begin{equation}
\begin{Bmatrix} N_i \\ S_i \\ M_i \\ N_j \\ S_j \\ M_j \end{Bmatrix}
=\begin{bmatrix}
EA/L & 0 & 0 & -EA/L & 0 & 0 \\
0 & 12EI/L^3 & 6EI/L^2 & 0 & -12EI/L^3 & 6EI/L^2 \\
0 & 6EI/L^2 & 4EI/L & 0 & -6EI/L^2 & 2EI/L \\
-EA/L & 0 & 0 & EA/L & 0 & 0 \\
0 & -12EI/L^3 & -6EI/L^2 & 0 & 12EI/L^3 & -6EI/L^2 \\
0 & 6EI/L^2 & 2EI/L & 0 & -6EI/L^2 & 4EI/L
\end{bmatrix}
\begin{Bmatrix} u_i \\ v_i \\ \theta_i \\ u_j \\ v_j \\ \theta_j \end{Bmatrix}
\end{equation}
平面骨組要素の節点慣性力ベクトル
物体力として慣性力を考え，節点ベクトルとして「加速度」を与える．
ここでは，最も単純な方法として，要素に作用する慣性力を単純に節点力に配分する．
要素に作用する全慣性力は，全体座標系での値として以下のとおりとなる．
\begin{align}
F_{m(x)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_h \cdot g) \\
F_{m(y)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_v \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_h \\
k_v \\
0 \\
k_h \\
k_v \\
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 \\
1 \\
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}]$ の転置マトリックス（この場合は逆行列に等しい）である．
平面骨組要素の節点外力ベクトル
節点外力ベクトルについては，ここで紹介するプログラムでは，直接，等価節点外力を入力データファイルから入力することとしており，特別な計算処理は行っていない．また荷重は，全体座標系における水平方向力，鉛直方向力，モーメントとしている．
座標変換マトリックス
これまでに求めた要素剛性方程式は，要素に固定された「要素座標系」でのものである．
実際の構造物では要素は色々な方向を向いているので，統一的に定義した「全体座標系」との関係を求めておく必要がある．
上図に示すように，全体座標系を $(X,Y)$，要素座標系を $(x,y)$ とする．
ここで要素座標系での変位を $(u,v)$ とすれば，全体座標系での変位 $(U,V)$ は以下の通りとなる．
\begin{align}
&U=u\cos\phi-v\sin\phi \\
&V=u\sin\phi+v\cos\phi
\end{align}
よって
\begin{equation}
\begin{Bmatrix} U \\ V \end{Bmatrix}=
\begin{bmatrix} \cos\phi & -\sin\phi \\ \sin\phi & \cos\phi \end{bmatrix}
\begin{Bmatrix} u \\ v \end{Bmatrix}
\end{equation}
したがって，全体座標系の変位を要素座標系の変位に変換するマトリックスは以下の通りとなる．
\begin{equation}
\begin{Bmatrix} u \\ v \end{Bmatrix}=
\begin{bmatrix} \cos\phi & \sin\phi \\ -\sin\phi & \cos\phi \end{bmatrix}
\begin{Bmatrix} U \\ V \end{Bmatrix}
\end{equation}
これを，回転角は全体座標系・要素座標系とも同じ（$\theta=\Theta$）ことに注意して，要素全自由度についての座標変換マトリックスとして表すと以下の通りとなる．
\begin{equation}
\begin{Bmatrix} u_i \\ v_i \\ \theta_i \\ u_j \\ v_j \\ \theta_j \end{Bmatrix}
=\begin{bmatrix}
\cos\phi & \sin\phi & 0 & 0 & 0 & 0 \\
-\sin\phi & \cos\phi & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & \cos\phi & \sin\phi & 0 \\
0 & 0 & 0 & -\sin\phi & \cos\phi & 0 \\
0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}
\begin{Bmatrix} U_i \\ V_i \\ \Theta_i \\ U_j \\ V_j \\ \Theta_j \end{Bmatrix}
\end{equation}
ここで，上で求めた全体座標系変位を要素座標系変位に変換する座標変換マトリックスを $[\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}]$ は，全体座標系変位・節点力を要素座標系変位・節点力に変換する座標変換マトリックスである．
\begin{equation}
[\boldsymbol{T}]
=\begin{bmatrix}
\cos\phi & \sin\phi & 0 & 0 & 0 & 0 \\
-\sin\phi & \cos\phi & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & \cos\phi & \sin\phi & 0 \\
0 & 0 & 0 & -\sin\phi & \cos\phi & 0 \\
0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}
\end{equation}
全体剛性方程式
以上より，全体座標系での剛性方程式は，以下のとおり表現される．
\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 軸正方向とする．
- 荷重として，節点外力，節点強制変位（ゼロ変位含む），節点温度変化，要素への慣性力を扱える
- 外力は，節点に作用するもののみ扱える．要素中間への荷重は扱えない．
- 連立一次方程式は，その大きさ（１節点自由度ｘ総節点数）を変えずに，節点強制変位導入による処理を行い解いている．
- 連立一次方程式は，numpy.linalg.solve(A,b) で解く
- 出力は，節点変位および要素断面力とする
- 入力データファイル，出力データファイルの文字区切りは空白とする
- 入力データファイル名，出力データファイル名は，コマンドライン引数として入力する
扱える境界条件の説明
| 項目 | 説明 |
|---|---|
| 節点外力 | 載荷節点数，積荷節点番号，荷重値を指定 |
| 節点温度変化 | 全節点の温度変化量を指定．温度上昇を正値とする |
| 要素慣性力 | 材料特性データの中で水平方向および鉛直方向の加速度を重力加速度に対する比で指定 |
| 節点強制変位 | 強制変位を与える節点数，節点番号，変位値を指定 |
プログラムの構成（プログラム名：py_fem_frame.py）
- モジュール読み込み
- 関数：データ入力 (def inpdata_frm)
(1) 基本数値読み込み
(2) 配列宣言
(3) 材料特性データ読み込み
(4) 要素構成節点番号および材料番号読み込み
(5) 節点座標および節点温度変化量読み込み
(6) 境界条件読み込み
(7) 節点荷重データ読み込み
- 関数：入力データ書き出し (def prinp_frm)
- 関数：計算結果書き出し (def prout_frm)
- 関数：要素剛性マトリックス作成 (def sm_frm)
- 関数：座標変換マトリックス作成 (def tm_frm)
- 関数：節点慣性力ベクトル作成 (def bfvec_frm)
- 関数：節点温度荷重ベクトル作成 (def tfvec_frm)
- 関数：要素断面力計算 (def calsecf_frm)
- 関数：メインルーチン
(1) コマンドラインより入出力ファイル名取得
(2) データ入力
(3) 入力データ書き出し
(4) 剛性マトリックス組み立て
(5) 境界条件処理（変位既知自由度の処理）
(6) 剛性方程式を解く（変位計算）
(7) 境界条件処理（既知変位の組み込み）
(8) 要素断面力計算
(9) 計算結果書き出し
(10) 計算時間等情報表示
- メインルーチン実行
実行コマンドと入出力データ
解析実行コマンド
python3 py_fem_frame.py inp.txt out.txt
| 変数 | 説明 |
|---|---|
| py_fem_frame.py | Python による FEM プログラム名 |
| inp.txt | 入力データファイル名（空白区切りデータ） |
| out.txt | 出力データファイル名（空白区切りデータ） |
入力データファイルフォーマット
npoin nele nsec npfix nlod # Basic values for analysis
E A I alpha gamma gkh gkv # Material properties
..... (1 to nsec) ..... #
node_1 node_2 isec # Element connectivity, material set number
..... (1 to nele) ..... #
x y deltaT # Node coordinate, temperature change of node
..... (1 to npoin) ..... #
lp fix_x fix_y fix_r rdis_x rdis_y rdis_r # Restricted node and restrict conditions
..... (1 to npfix) ..... #
lp fp_x fp_y fp_r # Loaded node and loading conditions
..... (1 to nlod) ..... # (omit data input if nlod=0)
| 変数 | 説明 |
|---|---|
| npoin, nele, nsec | 節点数，要素数，断面特性数 |
| npfix, nlod | 拘束節点数，載荷節点数 |
| E, A, I, alpha | 要素弾性係数，要素断面積，要素断面２次モーメント，線膨張係数 |
| gamma, gkh, gkv | 単位体積重量，水平及び鉛直方向加速度（gの比） |
| node_1, node_2, isec | 節点１，節点２，断面特性番号 |
| x, y, deltaT | 節点x座標，節点y座標，節点温度変化 |
| lp, fix_x, fix_y, fix_r | 拘束節点番号，x及びy方向拘束の有無（拘束: 1，自由: 0） |
| r_disx, r_disy, rdis_r | x及びy方向変位（無拘束でも0を入力） |
| lp, fp_x, fp_y, fp_r | 載荷重節点番号，x方向荷重，y方向荷重 |
出力データファイルフォーマット
npoin nele nsec npfix nlod
(Each value of above)
sec E A I alpha gamma gkh gkv
sec : Material number
E : Elastic modulus
A : Section area
I : Moment of inertia
alpha : Thermal expansion coefficient
gamma : Unit weight
gkh : Acceleration in horizontal direction (ratio to g)
gkv : Acceleration in vertical direction (ratio to g)
..... (1 to nsec) .....
node x y fx fy fr deltaT 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
deltaT : Temperature change of node
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) .....
node kox koy kor rdis_x rdis_y rdis_r
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)
kor : Index of restriction in rotation (0: free, 1: fixed)
rdis_x : Given displacement in x-direction
rdis_y : Given displacement in y-direction
rdis_r : Given displacement in rotation
..... (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-r
node : Node number
dis-x : Displacement in x-direction
dis-y : Displacement in y-direction
dis-r : Displacement 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) .....
n=(total degrees of freedom) time=(calculation time)
| 変数 | 説明 |
|---|---|
| dis-x, dis-y, dis-r | x方向変位，y方向変位，回転変位 |
| N, S, M | 軸力，せん断力，曲げモーメント |
プログラミング
プログラム本体を，「プログラムの構成」に沿って以下に示します．
1. モジュール読み込み
#======================
# 2D Frame Analysis
#======================
import numpy as np
import sys
import time
2. 関数：データ入力 (def inpdata_frm)
(1) 基本数値読み込み（npoin, nele, nsec, npfix, nlod）
(2) 配列宣言（ae, node, x, deltaT, mpfix, rdis, fp）
(3) 材料特性データ読み込み（配列 ae）
(4) 要素構成節点番号および材料番号読み込み（配列 node）
(5) 節点座標および節点温度変化量読み込み（配列 x, deltaT）
(6) 境界条件読み込み（配列 mpfix, rdis）
(7) 節点荷重データ読み込み（配列 fp）
def inpdata_frm(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((7,nsec),dtype=np.float64) # Section characteristics
node =np.zeros((nod+1,nele),dtype=np.int) # node-element relationship
x =np.zeros((2,npoin),dtype=np.float64) # Coordinates of nodes
deltaT=np.zeros(npoin,dtype=np.float64) # temparature 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]) # AA : Section area
ae[2,i]=float(text[2]) # AI : Moment of inertia
ae[3,i]=float(text[3]) # alpha : Thermal expansion coefficient
ae[4,i]=float(text[4]) # gamma : Unit weight
ae[5,i]=float(text[5]) # gkh : Acceleration in horizontal direction
ae[6,i]=float(text[6]) # gkv : Acceleration in vertical 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
deltaT[i]=float(text[2]) # temparature 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 rotation
rdis[0,lp-1]=float(text[4]) #fixed displacement in x-direction
rdis[1,lp-1]=float(text[5]) #fixed displacement in y-direction
rdis[2,lp-1]=float(text[6]) #fixed displacement in z-direction
# 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[3*lp-3]=float(text[1]) #load in x-direction
fp[3*lp-2]=float(text[2]) #load in y-direction
fp[3*lp-1]=float(text[3]) #load in rotation
f.close()
return npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp
3. 関数：入力データ書き出し (def prinp_frm)
def prinp_frm(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','A','I','alpha','gamma','gkh','gkv'),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:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>5s} {8:>5s} {9:>5s}'
.format('node','x','y','fx','fy','fr','deltaT','kox','koy','kor'),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:5d} {8:5d} {9:5d}'
.format(lp,x[0,i],x[1,i],fp[3*lp-3],fp[3*lp-2],fp[3*lp-1],deltaT[i],mpfix[0,i],mpfix[1,i],mpfix[2,i]),file=fout)
print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>15s} {5:>15s} {6:>15s}'
.format('node','kox','koy','kor','rdis_x','rdis_y','rdis_r'),file=fout)
for i in range(0,npoin):
if 0<mpfix[0,i]+mpfix[1,i]+mpfix[2,i]:
lp=i+1
print('{0:5d} {1:5d} {2:5d} {3:5d} {4:15.7e} {5:15.7e} {6:15.7e}'
.format(lp,mpfix[0,i],mpfix[1,i],mpfix[2,i],rdis[0,i],rdis[1,i],rdis[2,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_frm)
def prout_frm(fnameW,npoin,nele,disg,fsec):
fout=open(fnameW,'a')
# displacement
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'.format('node','dis-x','dis-y','dis-r'),file=fout)
for i in range(0,npoin):
lp=i+1
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'.format(lp,disg[3*lp-3],disg[3*lp-2],disg[3*lp-1]),file=fout)
# section force
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
.format('elem','N_i','S_i','M_i','N_j','S_j','M_j'),file=fout)
for ne in range(0,nele):
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
.format(ne+1,fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
fout.close()
5. 関数：要素剛性マトリックス作成 (def sm_frm)
\begin{equation}
[\boldsymbol{k}]
=\begin{bmatrix}
EA/L & 0 & 0 & -EA/L & 0 & 0 \\
0 & 12EI/L^3 & 6EI/L^2 & 0 & -12EI/L^3 & 6EI/L^2 \\
0 & 6EI/L^2 & 4EI/L & 0 & -6EI/L^2 & 2EI/L \\
-EA/L & 0 & 0 & EA/L & 0 & 0 \\
0 & -12EI/L^3 & -6EI/L^2 & 0 & 12EI/L^3 & -6EI/L^2 \\
0 & 6EI/L^2 & 2EI/L & 0 & -6EI/L^2 & 4EI/L
\end{bmatrix}
\end{equation}
def sm_frm(EA,EI,x1,y1,x2,y2):
ek=np.zeros((6,6),dtype=np.float64) # local stiffness matrix
xx=x2-x1
yy=y2-y1
el=np.sqrt(xx**2+yy**2)
ek[0,0]= EA/el; ek[0,1]= 0.0; ek[0,2]= 0.0; ek[0,3]=-EA/el; ek[0,4]= 0.0; ek[0,5]= 0.0
ek[1,0]= 0.0; ek[1,1]= 12*EI/el**3; ek[1,2]= 6*EI/el**2; ek[1,3]= 0.0; ek[1,4]=-12*EI/el**3; ek[1,5]= 6*EI/el**2
ek[2,0]= 0.0; ek[2,1]= 6*EI/el**2; ek[2,2]= 4*EI/el ; ek[2,3]= 0.0; ek[2,4]= -6*EI/el**2; ek[2,5]= 2*EI/el
ek[3,0]=-EA/el; ek[3,1]= 0.0; ek[3,2]= 0.0; ek[3,3]= EA/el; ek[3,4]= 0.0; ek[3,5]= 0.0
ek[4,0]= 0.0; ek[4,1]=-12*EI/el**3; ek[4,2]=-6*EI/el**2; ek[4,3]= 0.0; ek[4,4]= 12*EI/el**3; ek[4,5]=-6*EI/el**2
ek[5,0]= 0.0; ek[5,1]= 6*EI/el**2; ek[5,2]= 2*EI/el ; ek[5,3]= 0.0; ek[5,4]= -6*EI/el**2; ek[5,5]= 4*EI/el
return ek
6. 関数：座標変換マトリックス作成 (def tm_frm)
\begin{equation}
[\boldsymbol{T}]
=\begin{bmatrix}
\cos\phi & \sin\phi & 0 & 0 & 0 & 0 \\
-\sin\phi & \cos\phi & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & \cos\phi & \sin\phi & 0 \\
0 & 0 & 0 & -\sin\phi & \cos\phi & 0 \\
0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}
\end{equation}
def tm_frm(x1,y1,x2,y2):
tt=np.zeros((6,6),dtype=np.float64) # transformation matrix
xx=x2-x1
yy=y2-y1
el=np.sqrt(xx**2+yy**2)
s=yy/el
c=xx/el
tt[0,0]= c; tt[0,1]= s; tt[0,2]=0.0; tt[0,3]=0.0; tt[0,4]=0.0; tt[0,5]=0.0
tt[1,0]= -s; tt[1,1]= c; tt[1,2]=0.0; tt[1,3]=0.0; tt[1,4]=0.0; tt[1,5]=0.0
tt[2,0]=0.0; tt[2,1]=0.0; tt[2,2]=1.0; tt[2,3]=0.0; tt[2,4]=0.0; tt[2,5]=0.0
tt[3,0]=0.0; tt[3,1]=0.0; tt[3,2]=0.0; tt[3,3]= c; tt[3,4]= s; tt[3,5]=0.0
tt[4,0]=0.0; tt[4,1]=0.0; tt[4,2]=0.0; tt[4,3]= -s; tt[4,4]= c; tt[4,5]=0.0
tt[5,0]=0.0; tt[5,1]=0.0; tt[5,2]=0.0; tt[5,3]=0.0; tt[5,4]=0.0; tt[5,5]=1.0
return tt
7. 関数：節点慣性力ベクトル作成 (def bfvec_frm)
\begin{equation*}
\{\boldsymbol{F_{b}}\}=\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot A \cdot L}{2}
\begin{Bmatrix}
k_h \\
k_v \\
0 \\
k_h \\
k_v \\
0
\end{Bmatrix}
\end{equation*}
def bfvec_frm(gamma,A,gkh,gkv,x1,y1,x2,y2):
# Body force vector in global coordinate system
bfe_g=np.zeros(6,dtype=np.float64)
xx=x2-x1
yy=y2-y1
el=np.sqrt(xx**2+yy**2)
bfe_g[0]=0.5*gamma*A*el*gkh
bfe_g[1]=0.5*gamma*A*el*gkv
bfe_g[2]=0.0
bfe_g[3]=0.5*gamma*A*el*gkh
bfe_g[4]=0.5*gamma*A*el*gkv
bfe_g[5]=0.0
return bfe_g
8. 関数：節点温度荷重ベクトル作成 (def tfvec_frm)
\begin{equation*}
\{\boldsymbol{f_{t}}\}=EA\cdot\alpha\cdot\Delta T
\begin{Bmatrix}
-1 \\
0 \\
0 \\
1 \\
0 \\
0
\end{Bmatrix}
\end{equation*}
def tfvec_frm(EA,alpha,tem):
# Thermal load vector in local coordinate system
tfe_l=np.zeros(6,dtype=np.float64)
tfe_l[0]=-EA*alpha*tem
tfe_l[1]=0.0
tfe_l[2]=0.0
tfe_l[3]= EA*alpha*tem
tfe_l[4]=0.0
tfe_l[5]=0.0
return tfe_l
9. 関数：要素断面力計算 (def calsecf_frm)
$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_frm(EA,EI,alpha,tem,dis,x1,y1,x2,y2):
ek=sm_frm(EA,EI,x1,y1,x2,y2) # Stiffness matrix in local coordinate
tt=tm_frm(x1,y1,x2,y2) # Transformation matrix
secf=np.dot(ek,np.dot(tt,dis))
secf[0]=secf[0]+EA*alpha*tem
secf[3]=secf[3]-EA*alpha*tem
return secf
10. 関数：メインルーチン
メインルーチン処理手順
(1) コマンドラインより入出力ファイル名取得
コマンドライン引数として入出力データファイル名を取得．
入力データファイル名：fnameR，出力データファイル名：fnameW
(2) データ入力
関数 inpdata_frm により入力データ取得
(3) 入力データ書き出し
関数 prinp_frm により入力データ書き出し
(4) 剛性マトリックス組み立て
各要素毎に，剛性マトリックス ck （座標変換後），温度荷重ベクトル tfe （座標変換後），慣性力ベクトル bfe を計算し，全体剛性方程式に組み込んでいく．ここで剛性マトリックス，温度荷重ベクトル，慣性力ベクトル計算の関数の引数は，全てスカラーとした．要素番号が指定されれば，要素を構成する節点座標や要素物性値は特定できるため，他の要素の値を含んで収納する配列よりは，スカラーを送り込んだほうがプログラム的にスッキリするように思う．
全体剛性マトリックスは変数 gk に，荷重ベクトルは変数 fp に加算されていく．
要素剛性マトリックスを全体剛性マトリックスに組み込むための位置情報は，配列 ir[] に格納されている．
6 x 6 （6 は１節点自由度３と１要素節点数２の積）の次元の要素剛性マトリックスを NｘN （N は１節点自由度３と総節点数の積）の次元の全体剛性マトリックスにどのように配置するかを示したのが下図である．
碁盤目に区切られた大きな正方形は全体剛性マトリックス（N x N）を示し，$i$，$j$ は要素を構成する節点，$k_{11}\sim k_{66}$ は要素剛性マトリックス（6 x 6）の成分を示す．
(5) 境界条件処理（変位既知自由度の処理）
全体剛性マトリックスの大きさは変えず，変位既知節点の処理を行う．
| 変数 | 説明 |
|---|---|
| npoin | 総節点数 |
| nfree | １節点自由度（ここでは 3） |
| mpfix | 各節点毎の拘束条件を格納した配列（=1：拘束，=0：自由） |
| fp | 節点荷重ベクトル |
| rdis | 変位拘束する節点の変位量 |
| gk | 全体剛性マトリックス |
| disg | 全体剛性方程式の解（変位） |
(6) 剛性方程式を解く（変位計算）
numpy.linalg.solve() により連立一次方程式を解く．
(7) 境界条件処理（既知変位の組み込み）
連立一次方程式を解いた後，既知変位に相当する自由度の項に，既知変位を入れ直すことを忘れないようにする．
(8) 要素断面力計算
関数 calsecf_frm による要素毎の断面力計算実施
(9) 計算結果書き出し
関数 prout_frm による計算結果出力
(10) 計算時間等情報表示
以下により，計算完了時刻とスタート時刻の差分を処理時間として出力する．
import time
start=time.time()
......
dtime=time.time()-start
メインルーチンプログラム
def main_frm():
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=3 # Degree of freedom per node
# data input
npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_frm(fnameR,nod,nfree)
# print out of input data
prinp_frm(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]
x2=x[0,j]; y2=x[1,j]
ee =ae[0,m]
aa =ae[1,m]
ai =ae[2,m]
alpha=ae[3,m]
gamma=ae[4,m]
gkh =ae[5,m]
gkv =ae[6,m]
A =aa
EA=ee*aa
EI=ee*ai
tem=0.5*(deltaT[i]+deltaT[j])
tt=tm_frm(x1,y1,x2,y2) # Transformation matrix
ek=sm_frm(EA,EI,x1,y1,x2,y2) # Stiffness matrix in local coordinate
ck =np.dot(np.dot(tt.T,ek),tt) # Stiffness matrix in global coordinate
tfe_l=tfvec_frm(EA,alpha,tem) # Thermal load vector in local coordinate
tfe =np.dot(tt.T,tfe_l) # Thermal load vector in global coordinate
bfe=bfvec_frm(gamma,A,gkh,gkv,x1,y1,x2,y2) # Bpdy force vector in global coordinate
ir[5]=3*j+2; ir[4]=ir[5]-1; ir[3]=ir[4]-1
ir[2]=3*i+2; 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)
# 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(6,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]
x2=x[0,j]; y2=x[1,j]
ee =ae[0,m]
aa =ae[1,m]
ai =ae[2,m]
alpha=ae[3,m]
EA=ee*aa
EI=ee*ai
tem=0.5*(deltaT[i]+deltaT[j])
dis[0]=disg[3*i]; dis[1]=disg[3*i+1]; dis[2]=disg[3*i+2]
dis[3]=disg[3*j]; dis[4]=disg[3*j+1]; dis[5]=disg[3*j+2]
fsec[:,ne]=calsecf_frm(EA,EI,alpha,tem,dis,x1,y1,x2,y2)
# print out of result
prout_frm(fnameW,npoin,nele,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_frm()
（付録）Sympy による行列演算の確認
軸剛性の計算は，２ｘ２の行列なので，手計算でも簡単に行えるが，曲げ剛性の計算は４ｘ４の行列となるので，手計算では厳しい．
そこで，Sympyを使って曲げ成分剛性マトリックス演算を確認してみる．
(1) まず，次の行列の逆行列を求めてみる．
\begin{equation}
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
1 & L & L^2 & L^3 \\
0 & 1 & 2L & 3L^2
\end{bmatrix}
\end{equation}
from sympy import *
init_printing()
L, x, y = symbols('L x y')
Xnd= Matrix([
[1,0,0,0],
[0,1,0,0],
[1,L,L**2,L**3],
[0,1,2*L,3*L**2]
])
C = simplify(Xnd.inv())
C
上記コードを実行すると以下の結果が得られる．
\begin{equation}
\left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\- \frac{3}{L^{2}} & - \frac{2}{L} & \frac{3}{L^{2}} & - \frac{1}{L}\\ \frac{2}{L^{3}} & \frac{1}{L^{2}} & - \frac{2}{L^{3}} & \frac{1}{L^{2}}\end{matrix}\right]
\end{equation}
(2) 続いて，以下のコードを実行．これは内挿関数を求めるもの．
XX= Matrix([[1,x,x**2,x**3]])
N=XX*C
N
上記コードを実行すると以下の結果が得られる．内挿関数を計算できた．
\begin{equation}
\left[\begin{matrix}1 - \frac{3 x^{2}}{L^{2}} + \frac{2 x^{3}}{L^{3}} & x - \frac{2 x^{2}}{L} + \frac{x^{3}}{L^{2}} & \frac{3 x^{2}}{L^{2}} - \frac{2 x^{3}}{L^{3}} & - \frac{x^{2}}{L} + \frac{x^{3}}{L^{2}}\end{matrix}\right]
\end{equation}
(3) 続いて次を実行．ひずみー変位関係マトリックスを求めるもの．上で求めた内挿関数を２階微分して，$-y$ を乗じている．
dv=diff(N,x)
ddv=diff(dv,x)
eb=-y*ddv
eb
上記コードを実行すると以下の結果が得られる．ひずみー変位関係マトリックスを計算できた．
\begin{equation}
\left[\begin{matrix}- y \left(- \frac{6}{L^{2}} + \frac{12 x}{L^{3}}\right) & - y \left(- \frac{4}{L} + \frac{6 x}{L^{2}}\right) & - y \left(\frac{6}{L^{2}} - \frac{12 x}{L^{3}}\right) & - y \left(- \frac{2}{L} + \frac{6 x}{L^{2}}\right)\end{matrix}\right]
\end{equation}
(4) 次に以下を実行．ひずみー変位関係マトリックスの積を軸方向（$x$ 方向）に$0$から$L$まで積分している．
A=eb.T*eb
BB=integrate(A,(x,0,L))
BB
上記コードを実行すると以下の結果が得られる．$x$ 方向（要素軸方向）に積分しているため，全要素に $y^2$ の項が入っているが，これをくくりだして断面２次モーメントに置き換えれば，要素剛性マトリックスを構成するマトリックスが得られる．
\begin{equation}
\left[\begin{matrix}\frac{12 y^{2}}{L^{3}} & \frac{6 y^{2}}{L^{2}} & - \frac{12 y^{2}}{L^{3}} & \frac{6 y^{2}}{L^{2}}\\ \frac{6 y^{2}}{L^{2}} & \frac{4 y^{2}}{L} & - \frac{6 y^{2}}{L^{2}} & \frac{2 y^{2}}{L}\\ - \frac{12 y^{2}}{L^{3}} & - \frac{6 y^{2}}{L^{2}} & \frac{12 y^{2}}{L^{3}} & - \frac{6 y^{2}}{L^{2}}\\ \frac{6 y^{2}}{L^{2}} & \frac{2 y^{2}}{L} & - \frac{6 y^{2}}{L^{2}} & \frac{4 y^{2}}{L}\end{matrix}\right]
\end{equation}
以 上