2020/06/15投稿

はじめに

微小変形弾性理論に基づく、Pythonによる軸対称応力解析プログラムを紹介します。
プログラムの特徴は以下の通りです。

軸対称要素として１節点２自由度の四角形アイソパラメトリック要素を用いています。１要素の回転方向角は１ラジアンとしています。
当初は「回転軸を水平右向きを正、半径方向軸を鉛直上向きを正」として、四角形要素を用いた２次元応力解析プログラム（https://qiita.com/damyarou/items/320bad2052bb5fccd75f）を書き換えて使っていましたが、「回転軸を鉛直上向きを正、半径方向軸を水平右向きを正」とした場合、変位及び応力が逆向きに求まってしまうため、鉛直軸の向きにより計算を補正するしくみを導入しました。
ある程度大規模な解析ができるよう、連立方程式は疎行列処理を行って解いています。

要素剛性方程式

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

\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}

要素剛性方程式の基本形

上図に示すような直交極座標系 $(z-r)$ 平面上の四角形要素を考える．ここに $ｚ$ は回転軸とし、水平右向きを正、$r$ は半径方向軸とし、鉛直上向きを正としておく。また、１要素の回転方向角は単位角度、すなわち１ラジアンとし、定式化・積分計算などの簡略化を図る。
軸対称解析用四角形要素は２次元応力解析用と同様に、１要素４節点とし，１節点２自由度（$z$方向変位および$r$方向変位）を有するものとする．よって基本的な要素剛性方程式は，以下の形をとる．ここに，四角形要素を構成する節点番号を $(z-r)$ 平面上で、反時計回りに，$i$，$j$，$k$，$l$ とする．

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

\begin{equation}
\begin{Bmatrix} f_{z,i} \\ f_{r,i} \\ f_{z,j} \\ f_{r,j} \\ f_{z,k} \\ f_{r,k} \\ f_{z,l} \\ f_{r,l} \end{Bmatrix}
=\begin{bmatrix}
k_{11} & k_{12} & k_{13} & k_{14} & k_{15} & k_{16} & k_{17} & k_{18} \\
k_{21} & k_{22} & k_{23} & k_{24} & k_{25} & k_{26} & k_{27} & k_{28} \\
k_{31} & k_{32} & k_{33} & k_{34} & k_{35} & k_{36} & k_{37} & k_{38} \\
k_{41} & k_{42} & k_{43} & k_{44} & k_{45} & k_{46} & k_{47} & k_{48} \\
k_{51} & k_{52} & k_{53} & k_{54} & k_{55} & k_{56} & k_{57} & k_{58} \\
k_{61} & k_{62} & k_{63} & k_{64} & k_{65} & k_{66} & k_{67} & k_{68} \\
k_{71} & k_{72} & k_{73} & k_{74} & k_{75} & k_{76} & k_{77} & k_{78} \\
k_{81} & k_{82} & k_{83} & k_{84} & k_{85} & k_{86} & k_{87} & k_{88} 
\end{bmatrix}
\begin{Bmatrix} w_i \\ u_i \\ w_j \\ u_j \\ w_k \\ u_k \\ w_l \\ u_l \end{Bmatrix}
\end{equation}

\begin{align}
&f_{z,i} & &\text{：節点$i$の$z$方向荷重} & & w_i & &\text{：節点$i$の$z$方向変位}\\
&f_{r,i} & &\text{：節点$i$の$r$方向荷重} & & u_i & &\text{：節点$i$の$r$方向変位}\\
&f_{z,j} & &\text{：節点$j$の$z$方向荷重} & & w_i & &\text{：節点$j$の$z$方向変位}\\
&f_{r,j} & &\text{：節点$j$の$r$方向荷重} & & u_i & &\text{：節点$j$の$r$方向変位}\\
&f_{z,k} & &\text{：節点$k$の$z$方向荷重} & & w_i & &\text{：節点$k$の$z$方向変位}\\
&f_{r,k} & &\text{：節点$k$の$r$方向荷重} & & u_i & &\text{：節点$k$の$r$方向変位}\\
&f_{z,l} & &\text{：節点$l$の$z$方向荷重} & & w_l & &\text{：節点$l$の$z$方向変位}\\
&f_{r,l} & &\text{：節点$l$の$r$方向荷重} & & u_l & &\text{：節点$l$の$r$方向変位}
\end{align}

線形弾性論による基本式

軸対称問題における応力-ひずみ関係とひずみ-変位関係

回転軸方向の応力・ひずみを添字 $z$，半径方向の応力・ひずみを添字 $r$，円周方向の応力・ひずみを添字 $\theta$ により表すとすれば，軸対称問題における応力-ひずみ関係式は，以下の通りとなる．

\begin{equation}
\begin{cases}
\epsilon_z-\alpha T=\cfrac{1}{E}\{\sigma_z-\nu(\sigma_r+\sigma_{\theta})\} \\
\epsilon_r-\alpha T=\cfrac{1}{E}\{\sigma_r-\nu(\sigma_z+\sigma_{\theta})\} \\
\epsilon_{\theta}-\alpha T=\cfrac{1}{E}\{\sigma_{\theta}-\nu(\sigma_z+\sigma_r)\} \\
\gamma_{zr}=\cfrac{2(1+\nu)}{E}\cdot\tau_{zr}
\end{cases}
\end{equation}

ここに，$\nu$ はポアソン比，$\alpha$ は熱膨張率，$T$ は温度変化量である．

上式を変形することにより

\begin{equation}
\begin{cases}
\sigma_z=\cfrac{E}{(1+\nu)(1-2\nu)}\{(1-\nu)\epsilon_z+\nu\epsilon_r+\nu\epsilon_{\theta}-(1+\nu)\alpha T\} \\
\sigma_r=\cfrac{E}{(1+\nu)(1-2\nu)}\{\nu\epsilon_z+(1-\nu)\epsilon_r+\nu\epsilon_{\theta}-(1+\nu)\alpha T\} \\
\sigma_{\theta}=\cfrac{E}{(1+\nu)(1-2\nu)}\{\nu\epsilon_z+\nu\epsilon_r+(1-\nu)\epsilon_{\theta}-(1+\nu)\alpha T\} \\
\tau_{zr}=\cfrac{E}{2(1+\nu)}\cdot\gamma_{zr}
\end{cases}
\end{equation}

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

の形で表せば，

\begin{equation}
\{\boldsymbol{\sigma}\}=\begin{Bmatrix} \sigma_z \\ \sigma_r \\ \sigma_{\theta} \\ \tau_{zr} \end{Bmatrix} \qquad
\{\boldsymbol{\epsilon}\}=\begin{Bmatrix} \epsilon_z \\ \epsilon_r \\ \epsilon_{\theta} \\ \gamma_{zr} \end{Bmatrix} \qquad
\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha T \\ \alpha T \\ \alpha T \\ 0 \end{Bmatrix}
\end{equation}

\begin{equation}
[\boldsymbol{D_e}]=\cfrac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix}
1-\nu & \nu & \nu & 0 \\
\nu & 1-\nu & \nu & 0 \\
\nu & \nu & 1-\nu & 0 \\
0 & 0 & 0 & \cfrac{1-2\nu}{2}
\end{bmatrix}
\end{equation}

ひずみー変位関係式は、回転軸方向の変位 $w$ および半径方向の変位 $u$ を用いて以下の通り表される。

\begin{equation}
\{\boldsymbol{\epsilon}\}
=\begin{Bmatrix} \epsilon_z \\ \epsilon_r \\ \epsilon_{\theta} \\ \gamma_{zr} \end{Bmatrix}
=\begin{Bmatrix} \cfrac{\partial w}{\partial z} \\ \cfrac{\partial u}{\partial r} \\ \cfrac{u}{r} \\ \cfrac{\partial w}{\partial r}+\cfrac{\partial u}{\partial z} \end{Bmatrix}
\end{equation}

四角形要素を用いた定式化

ひずみ-変位関係行列の導入

軸対称要素内任意点の回転軸方向変位 $w$ および半径方向変位 $u$ を以下のように仮定する．ここに座標 ($a$,$b$) は [$-1$,$1$] の範囲を有する正規化パラメータ座標，$w_{i,j,k,l}$ および $u_{i,j,k,l}$ は要素を構成する節点の変位を示す．

\begin{align}
w=&N_1(a,b)\cdot w_i+N_2(a,b)\cdot w_j+N_3(a,b)\cdot w_k+N_4(a,b)\cdot w_l \\
u=&N_1(a,b)\cdot u_i+N_2(a,b)\cdot u_j+N_3(a,b)\cdot u_k+N_4(a,b)\cdot u_l
\end{align}

要素内任意点の変位 $\{w, u\}$ は，節点変位 $\{\boldsymbol{u_{nd}}\}$ を用いて以下の通り表される。，

\begin{equation}
\begin{Bmatrix} w \\ u \end{Bmatrix}
=\begin{bmatrix}
N_1 & 0   & N_2 & 0   & N_3 & 0   & N_4 & 0   \\
0   & N_1 & 0   & N_2 & 0   & N_3 & 0   & N_4
\end{bmatrix}
\begin{Bmatrix}
w_i \\ u_i \\ w_j \\ u_j \\ w_k \\ u_k \\ w_l \\ u_l
\end{Bmatrix}
=[\boldsymbol{N}]\{\boldsymbol{u_{nd}}\}
\end{equation}

\begin{align}
N_1=&\frac{1}{4}(1-a)(1-b) \qquad N_2=\frac{1}{4}(1+a)(1-b) \\
N_3=&\frac{1}{4}(1+a)(1+b) \qquad N_4=\frac{1}{4}(1-a)(1+b)
\end{align}

\begin{align}
&\cfrac{\partial N_1}{\partial a}=-\cfrac{1}{4}(1-b)  &\cfrac{\partial N_2}{\partial a}=+\cfrac{1}{4}(1-b) & &\cfrac{\partial N_3}{\partial a}=+\cfrac{1}{4}(1+b) & &\cfrac{\partial N_4}{\partial a}=-\cfrac{1}{4}(1+b)\\
&\cfrac{\partial N_1}{\partial b}=-\cfrac{1}{4}(1-a)  &\cfrac{\partial N_2}{\partial b}=-\cfrac{1}{4}(1+a) & &\cfrac{\partial N_3}{\partial b}=+\cfrac{1}{4}(1+a) & &\cfrac{\partial N_4}{\partial b}=+\cfrac{1}{4}(1-a)
\end{align}

要素内任意点のひずみ $\{\boldsymbol{\epsilon}\}$ は，節点変位 $\{\boldsymbol{u_{nd}}\}$ を用いて，

\begin{equation}
\{\boldsymbol{\epsilon}\}=\begin{Bmatrix} \cfrac{\partial w}{\partial z} \\ \cfrac{\partial u}{\partial r} \\ \cfrac{u}{r} \\ \cfrac{\partial w}{\partial r}+\cfrac{\partial u}{\partial z} \end{Bmatrix}
=\begin{bmatrix}
\cfrac{\partial N_1}{\partial z} & 0 & \cfrac{\partial N_2}{\partial z} & 0 & \cfrac{\partial N_3}{\partial z} & 0 & \cfrac{\partial N_4}{\partial z} & 0 \\
0 & \cfrac{\partial N_1}{\partial r} & 0 & \cfrac{\partial N_2}{\partial r} & 0 & \cfrac{\partial N_3}{\partial r} & 0 & \cfrac{\partial N_4}{\partial r} \\
0 & \cfrac{N_1}{r} & 0 & \cfrac{N_2}{r} & 0 & \cfrac{N_3}{r} & 0 & \cfrac{N_4}{r} \\
\cfrac{\partial N_1}{\partial r} & \cfrac{\partial N_1}{\partial z} & \cfrac{\partial N_2}{\partial r} & \cfrac{\partial N_2}{\partial z} & \cfrac{\partial N_3}{\partial r} & \cfrac{\partial N_3}{\partial z} & \cfrac{\partial N_4}{\partial r} & \cfrac{\partial N_4}{\partial z}
\end{bmatrix}
\begin{Bmatrix}
w_i \\ u_i \\ w_j \\ u_j \\ w_k \\ u_k \\ w_l \\ u_l
\end{Bmatrix}
=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\}
\end{equation}

ここで，内挿関数 $N_i$ に関する微分は，Jacobi行列 $[J]$ を用いて以下の通り表せる．

\begin{equation}
\begin{Bmatrix} \cfrac{\partial N_i}{\partial a} \\ \cfrac{\partial N_i}{\partial b} \end{Bmatrix}
=[J] \begin{Bmatrix} \cfrac{\partial N_i}{\partial z} \\ \cfrac{\partial N_i}{\partial r} \end{Bmatrix} \qquad
\begin{Bmatrix} \cfrac{\partial N_i}{\partial z} \\ \cfrac{\partial N_i}{\partial r} \end{Bmatrix}
=[J]^{-1} \begin{Bmatrix} \cfrac{\partial N_i}{\partial a} \\ \cfrac{\partial N_i}{\partial b} \end{Bmatrix}
\end{equation}

\begin{equation}
[J]=\begin{bmatrix}
\cfrac{\partial z}{\partial a} & \cfrac{\partial r}{\partial a} \\
\cfrac{\partial z}{\partial b} & \cfrac{\partial r}{\partial b}
\end{bmatrix}
=\begin{bmatrix}
\displaystyle \sum_{i=1}^4\left(\cfrac{\partial N_i}{\partial a}z_i\right) & \displaystyle \sum_{i=1}^4\left(\cfrac{\partial N_i}{\partial a}r_i\right) \\
\displaystyle \sum_{i=1}^4\left(\cfrac{\partial N_i}{\partial b}z_i\right) & \displaystyle \sum_{i=1}^4\left(\cfrac{\partial N_i}{\partial b}r_i\right)
\end{bmatrix}
=\begin{bmatrix} J_{11} & J_{12} \\ J_{21} & J_{22} \end{bmatrix}
\end{equation}

\begin{equation}
[J]^{-1}=\frac{1}{\det(J)}\begin{bmatrix} J_{22} & -J_{12} \\ -J_{21} & J_{11} \end{bmatrix} \qquad
\det(J)=J_{11}\cdot J_{22}-J_{12}\cdot J_{21}
\end{equation}

以上より、ひずみ-変位関係行列$[\boldsymbol{B}]$の要素は、以下の通り表せる。

\begin{equation}
\cfrac{\partial N_i}{\partial z}=\cfrac{1}{\det(J)}\left\{J_{22}\cfrac{\partial N_i}{\partial a}-J_{12}\cfrac{\partial N_i}{\partial b}\right\} \qquad
\cfrac{\partial N_i}{\partial r}=\cfrac{1}{\det(J)}\left\{-J_{21}\cfrac{\partial N_i}{\partial a}+J_{11}\cfrac{\partial N_i}{\partial b}\right\}
\end{equation}

なお，この段階では$[\boldsymbol{B}]$は変数$a$および$b$の関数として定義されていることに注意する．

以下に$[J]$の要素を示しておく．ここに$z_{i,j,k,l}$および$r_{i,j,k,l}$は要素を構成する節点座標である．

\begin{align*}
&J_{11}=\frac{\partial N_1}{\partial a} z_i+\frac{\partial N_2}{\partial a} z_j+\frac{\partial N_3}{\partial a} z_k+\frac{\partial N_4}{\partial a} z_l \\
&J_{12}=\frac{\partial N_1}{\partial a} r_i+\frac{\partial N_2}{\partial a} r_j+\frac{\partial N_3}{\partial a} r_k+\frac{\partial N_4}{\partial a} r_l \\
&J_{21}=\frac{\partial N_1}{\partial b} z_i+\frac{\partial N_2}{\partial b} z_j+\frac{\partial N_3}{\partial b} z_k+\frac{\partial N_4}{\partial b} z_l \\
&J_{22}=\frac{\partial N_1}{\partial b} r_i+\frac{\partial N_2}{\partial b} r_j+\frac{\partial N_3}{\partial b} r_k+\frac{\partial N_4}{\partial b} r_l
\end{align*}

また，要素内任意点の$r$座標値は，以下により計算する．

\begin{equation}
r=N_1\cdot r_i+N_2\cdot r_j+N_3\cdot r_k+N_4\cdot r_l
\end{equation}

ここに，$r_i$，$r_j$，$r_k$，$r_l$は要素を構成する節点の$r$座標値である．

ガウス積分

以降の要素剛性マトリックスや荷重ベクトル算定で必要となる要素に対する面積積分は，以下に示すガウス積分により行う．

積分点を 4 点 ($n=2$) とすれば，$a$，$b$ および重み $H$ の値は下表の通りであり，$a$，$b$ を変化させた 4 回の値の総和が積分の近似値となる．また，重みが $H=1$ であることから，計算は更に単純になる．

\begin{equation}
\int_{-1}^{1}\int_{-1}^{1} f(a,b)\cdot da\cdot db
\doteqdot \sum_{i=1}^n\sum_{j=1}^n H_i\cdot H_j\cdot f(a_i,b_j)
= \sum_{kk=1}^4 f(a_{kk}, b_{kk})
\end{equation}

\begin{align}
& i  & j & &  a              & &  b              & & H              & & kk \\
& 1  & 1 & & -0.57735 02692  & & -0.57735 02692  & & 1.00000 00000  & & 1  \\\
& 1  & 2 & & +0.57735 02692  & & -0.57735 02692  & & 1.00000 00000  & & 2  \\\
& 2  & 1 & & +0.57735 02692  & & +0.57735 02692  & & 1.00000 00000  & & 3  \\\
& 2  & 2 & & -0.57735 02692  & & +0.57735 02692  & & 1.00000 00000  & & 4
\end{align}

ここで、軸対称要素の場合，回転中心より半径 $r$ の位置にある微小要素の体積は，$r\times d\theta\times dr\times dz$ となる．積分値や荷重ベクトルなどを，円周方向は1ラジアンで評価することとし，Gauss積分を用いるための積分変数変換を行えば，

\begin{equation}
r\times d\theta\times dr\times dz=r\times dr\times dz=r\times \det(J)\times da\times db
\end{equation}

\begin{equation}
r=N_1\cdot r_i+N_2\cdot r_j+N_3\cdot r_k+N_4\cdot r_l
\end{equation}

ここに，$r_i$，$r_j$，$r_k$，$r_l$は要素を構成する節点の$r$座標値である．

要素剛性マトリックス

4節点パラメトリック要素の剛性行列$[\boldsymbol{k}]$は，応力-ひずみ関係行列を$[\boldsymbol{D_e}]$とすれば以下の通りとなる．
なお、１要素の回転方向角度は、１ラジアンとしているため、積分計算上は現れない。

\begin{align}
[\boldsymbol{k}]=&\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dV \\
=&\int_{-1}^{1}\int_{-1}^{1}[\boldsymbol{B}]^{T}[\boldsymbol{D_e}][\boldsymbol{B}]\cdot r\cdot\det(J)\cdot da\cdot db \\
=&\sum_{kk=1}^4 \left\{[\boldsymbol{B}]^{T}[\boldsymbol{D_e}][\boldsymbol{B}]\cdot r\cdot\det(J)\right\}_{kk}
\end{align}

節点温度荷重ベクトル

初期ひずみとして温度ひずみを考える．
温度変化によるひずみ${\boldsymbol{\epsilon_0}}$による節点力ベクトル${\boldsymbol{f_t}}$は以下の通りとなる．

\begin{align}
\{\boldsymbol{f_t}\}=&\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV \\
=&\int_{-1}^{1}\int_{-1}^{1}[\boldsymbol{B}]^{T}[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}\cdot r\cdot\det(J)\cdot da\cdot db \\
=&\sum_{kk=1}^4 \left\{[\boldsymbol{B}]^{T}[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}\cdot r\cdot\det(J)\right\}_{kk}
\end{align}

温度変化によるひずみの成分は以下のとおりである．

\begin{equation}
T_p=N_1\cdot\phi_i+N_2\cdot\phi_j+N_3\cdot\phi_k+N_4\cdot\phi_l
\end{equation}

\begin{equation}
\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_p \\ \alpha\cdot T_p \\ \alpha\cdot T_p \\ 0 \end{Bmatrix}
\end{equation}

ここに，$\alpha$は線膨張係数，$T_p$は要素任意点の温度変化量，${\phi_i~~\phi_j~~\phi_k~~\phi_l}^T$は節点温度変化量，$[N_1~~N_2~~N_3~~N_4]$は$[\boldsymbol{N}]$の要素である．また，温度変化量の符号は，温度上昇を正とする．

節点慣性力ベクトル

\begin{align}
\{\boldsymbol{f_b}\}=&\gamma \cdot \left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} \\
=&\gamma\cdot\int_{A}[\boldsymbol{N}]^T[\boldsymbol{N}]dA\cdot\{\boldsymbol{w_{nd}}\} \\
=&\gamma\cdot\sum_{kk=1}^4 \left\{[\boldsymbol{N}]^T[\boldsymbol{N}]\cdot r\cdot\det(J)\right\}_ {kk}\cdot\{\boldsymbol{w_{nd}}\}
\end{align}

\begin{equation}
\{\boldsymbol{w_{nd}}\}=\begin{Bmatrix}k_z & 0 & k_z & 0 & k_z & 0 & k_z & 0 \end{Bmatrix}^T
\end{equation}

ここに $\gamma$ は要素の単位体積重量（質量$\times$重力加速度），$k_z$ は回転軸方向加速度（重力加速度に対する比率）である．なお、半径方向慣性力を作用させることは実用上ないと思われるため、ここではゼロとして扱っている。

節点外力ベクトル

節点外力ベクトルについては，ここで紹介するプログラムでは，直接，等価節点外力を入力データファイルから入力することとしており，特別な計算処理は行っていない．

なお、軸対称解析では、荷重の入力方法には注意を要する。入力荷重の考え方の一例を以下に示す。

(入力荷重の考え方の一例)

半径r=2000mm，軸方向長さz=500mmの円筒1要素にp=1N/mm$^2$(=1MPa)の内圧が作用する場合，要素の内壁に作用する全荷重は，以下の通りとなる．

\begin{equation*}
p\times r\times 1(rad)\times z=1(N/mm^2)\times 2,000(mm)\times 1(rad)\times 500(mm)=1,000,000(N)
\end{equation*}

このため，等価節点力の考え方により，(z,r)=(0,2,000)の節点に500,000(N)，(z,r)=(500,2,000)の節点に500,000(N)を半径正方向に載荷すればよい．

要素応力ベクトル

要素応力ベクトル ${\boldsymbol{\sigma}}$ は，以下のように計算できる．

\begin{align}
&\{\boldsymbol{\epsilon}\}=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\} & &\text{（節点変位から求まるひずみ）}\\
&T_p=N_1\cdot\phi_i+N_2\cdot\phi_j+N_3\cdot\phi_k+N_4\cdot\phi_l & &\text{（要素内任意点の温度変化量）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_p \\ \alpha\cdot T_p \\ \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{（温度ひずみ）}\\
&\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\} & &\text{（要素応力ベクトル）}
\end{align}

ここに，$N_1 \sim N_4$ は内挿関数，$\phi_i$，$\phi_j$，$\phi_k$，$\phi_l$ は要素節点 $i$，$j$，$k$，$l$ の温度変化量である．

プログラムでの座標系の扱い

以下に座標系と要素のイメージ図を示す。
(Case-1)および(Case-2)では、回転軸をｚ軸、半径方向軸をｒ軸としている。

(Reference)

(Reference) で示された図は通常の２次元応力解析用の座標である。ｘ軸を水平右向き、ｙ軸を鉛直上向きに設定している。
この場合、通常、要素番号は、節点番号を反時計回りに指定して要素を定義する。

(Case-1)

(Case-1)は、ｚ軸（回転軸）を水平右向き、ｒ軸（半径方向）を鉛直上向きに設定している。
この場合も、(Reference)の２次元応力解析と同様に、要素番号は、節点番号を反時計回りに指定して要素を定義する。（そうするようにプログラムが組まれている）

(Case-2)

(Case-2)は、ｚ軸（回転軸）を鉛直上向き、ｒ軸（半径方向）を水平右向きに設定している。
この場合、(Reference) の２次元応力解析と同様に反時計回りの節点番号で要素を定義するすると、(Case-1)と比べてｚ軸とｒ軸が入れ替わっているため、要素の節点番号が逆周りで定義されたことになり、剛性マトリックスや物体力ベクトル、温度荷重ベクトルなど要素座標を用いて計算する行列やベクトルが通常と逆符号で計算されることとなる。

対応策

このため、今回整備した軸対称応力解析プログラムでは、変数 nzdir を導入し、座標の方向に応じて、剛性マトリックス、物体力ベクトル、温度荷重ベクトルの方向を逆転させることにした。
すなわち、上図(Case-1)では、nzdir=1、上図(Case-2)ではnzdir=-1とし、メインプログラム上で、求まった剛性マトリックス、物体力ベクトル、温度荷重ベクトルの符号を調整するようにしている。

        sm=sm*float(nzdir) # 剛性マトリックス
        tfe=tfe*float(nzdir) # 温度荷重ベクトル
        bfe=bfe*float(nzdir) # 物体力（慣性力）ベクトル

入力データ・フォーマット

npoin  nele  nsec  npfix  nlod  nzdir # Basic values for analysis
E  po  alpha  gamma  gkz              # Material properties
    ..... (1 to nsec) .....           #
node1  node2  node3  node4  isec      # Element connectivity, material set number
    ..... (1 to nele) .....           #   (counter-clockwise order of node number)
z  r  deltaT                          # Coordinate, Temperature change of node
    ..... (1 to npoin) .....          #
node  koz  kor  rdisz  rdisr          # Restricted node and restrict conditions
    ..... (1 to npfix) .....          #
node  fz  fr                          # Loaded node and loading conditions
    ..... (1 to nlod) .....           #   (omit data input if nlod=0)

変数	説明
npoin, nele, nsec	節点数，要素数，材料特性数
npfix, nlod	拘束節点数，載荷節点数
nzdir	ｚ軸方向（水平右向き: １，鉛直上向き: −１）
E, po, alpha	弾性係数，ポアソン比，線膨張係数
gamma, gkz	単位体積重量，ｚ方向加速度（gの比）
node1, node2, node3, node4, isec	節点1，節点2，節点3，節点4，材料特性番号
z, r, deltaT	節点z座標，節点r座標，節点温度変化
node, koz, kor	拘束節点番号，z及びr方向拘束の有無（拘束: 1，自由: 0）
rdisz, rdisr	z及びr方向変位（無拘束でも0を入力）
node, fz, fr	載荷重節点番号，z方向荷重，r方向荷重

（注）回転軸ｚの方向が水平（nzdir=1）であっても鉛直（nzdir=-1）であっても、座標入力は (z , r) の順に行う。

プログラム全文

################################
# Axisymmetric Stress Analysis
################################
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time


def inpdata_ax4(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
    nzdir =int(text[5]) # direction of z-axis (=1 or =-1)
    # array declanation
    ae    =np.zeros([5,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)         # Temperature change of node
    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]) #alpha: Thermal expansion coefficient
        ae[3,i]=float(text[3]) #gamma: Unit weight of material
        ae[4,i]=float(text[4]) #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]) #node_3
        node[3,i]=int(text[3]) #node_4
        node[4,i]=int(text[4]) #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]) # z-coordinate
        x[1,i]   =float(text[1]) # r-coordinate
        deltaT[i]=float(text[2]) # Temperature change of node
    # boundary conditions (0:free, 1:restricted)
    if 0<npfix:
        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 z-direction
            mpfix[1,lp-1]=int(text[2])  #fixed in r-direction
            rdis[0,lp-1]=float(text[3]) #fixed displacement in z-direction
            rdis[1,lp-1]=float(text[4]) #fixed displacement in r-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[2*lp-2]=float(text[1]) #load in z-direction
            fp[2*lp-1]=float(text[2]) #load in r-direction
    f.close()
    return npoin,nele,nsec,npfix,nlod,nzdir,ae,node,x,deltaT,mpfix,rdis,fp


def prinp_ax4(fnameW,npoin,nele,nsec,npfix,nlod,nzdir,ae,x,fp,deltaT,mpfix,rdis,node):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}'
    .format('npoin','nele','nsec','npfix','nlod','nzdir'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}'
    .format(npoin,nele,nsec,npfix,nlod,nzdir),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s}'
    .format('sec','E','po','alpha','gamma','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}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s}'
    .format('node','z','r','fz','fr','deltaT','koz','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:5d} {7:5d}'
        .format(lp,x[0,i],x[1,i],fp[2*i],fp[2*i+1],deltaT[i],mpfix[0,i],mpfix[1,i]),file=fout)
    if 0<npfix:
        print('{0:>5s} {1:>5s} {2:>5s} {3:>15s} {4:>15s}'
        .format('node','koz','kor','rdis_z','rdis_r'),file=fout)
        for i in range(0,npoin):
            if 0<mpfix[0,i]+mpfix[1,i]:
                lp=i+1
                print('{0:5d} {1:5d} {2:5d} {3:15.7e} {4:15.7e}'
                .format(lp,mpfix[0,i],mpfix[1,i],rdis[0,i],rdis[1,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}'
    .format('elem','i','j','k','l','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}'
        .format(ne+1,node[0,ne],node[1,ne],node[2,ne],node[3,ne],node[4,ne]),file=fout)
    fout.close()


def prout_ax4(fnameW,npoin,nele,disg,strs):
    fout=open(fnameW,'a')
    # displacement
    print('{0:>5s} {1:>15s} {2:>15s}'.format('node','dis-z','dis-r'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e}'.format(lp,disg[2*lp-2],disg[2*lp-1]),file=fout)
    # section force
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'
    .format('elem','sig_z','sig_r','sig_t','tau_zr','p1','p2','ang'),file=fout)
    for ne in range(0,nele):
        sigz =strs[0,ne]
        sigr =strs[1,ne]
        sigt =strs[2,ne]
        tauzr=strs[3,ne]
        ps1,ps2,ang=pst_ax4(sigz,sigr,tauzr)
        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(ne+1,sigz,sigr,sigt,tauzr,ps1,ps2,ang),file=fout)
    fout.close()


def dmat_ax4(E,po):
    d=np.zeros([4,4],dtype=np.float64)
    d[0,0]=1.0-po; d[0,1]=po    ; d[0,2]=po    ; d[0,3]=0.0
    d[1,0]=po    ; d[1,1]=1.0-po; d[1,2]=po    ; d[1,3]=0.0
    d[2,0]=po    ; d[2,1]=po    ; d[2,2]=1.0-po; d[2,3]=0.0
    d[3,0]=0.0   ; d[3,1]=0.0   ; d[3,2]=0.0   ; d[3,3]=0.5*(1.0-2.0*po)
    d=d*E/(1.0+po)/(1.0-2.0*po)
    return d


def bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4):
    bm=np.zeros([4,8],dtype=np.float64)
    #[N]
    nm1=0.25*(1.0-a)*(1.0-b)
    nm2=0.25*(1.0+a)*(1.0-b)
    nm3=0.25*(1.0+a)*(1.0+b)
    nm4=0.25*(1.0-a)*(1.0+b)
    rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4
    #dN/da,dN/db
    dn1a=-0.25*(1.0-b)
    dn2a= 0.25*(1.0-b)
    dn3a= 0.25*(1.0+b)
    dn4a=-0.25*(1.0+b)
    dn1b=-0.25*(1.0-a)
    dn2b=-0.25*(1.0+a)
    dn3b= 0.25*(1.0+a)
    dn4b= 0.25*(1.0-a)
    #Jacobi matrix and det(J)
    J11=dn1a*z1 + dn2a*z2 + dn3a*z3 + dn4a*z4
    J12=dn1a*r1 + dn2a*r2 + dn3a*r3 + dn4a*r4
    J21=dn1b*z1 + dn2b*z2 + dn3b*z3 + dn4b*z4
    J22=dn1b*r1 + dn2b*r2 + dn3b*r3 + dn4b*r4
    detJ=J11*J22-J12*J21
    #[B]=[dN/dx][dN/dy]
    bm[0,0]= J22*dn1a-J12*dn1b
    bm[0,2]= J22*dn2a-J12*dn2b
    bm[0,4]= J22*dn3a-J12*dn3b
    bm[0,6]= J22*dn4a-J12*dn4b
    bm[1,1]=-J21*dn1a+J11*dn1b
    bm[1,3]=-J21*dn2a+J11*dn2b
    bm[1,5]=-J21*dn3a+J11*dn3b
    bm[1,7]=-J21*dn4a+J11*dn4b
    bm[2,1]=nm1/rr*detJ
    bm[2,3]=nm2/rr*detJ
    bm[2,5]=nm3/rr*detJ
    bm[2,7]=nm4/rr*detJ
    bm[3,0]=-J21*dn1a+J11*dn1b
    bm[3,1]= J22*dn1a-J12*dn1b
    bm[3,2]=-J21*dn2a+J11*dn2b
    bm[3,3]= J22*dn2a-J12*dn2b
    bm[3,4]=-J21*dn3a+J11*dn3b
    bm[3,5]= J22*dn3a-J12*dn3b
    bm[3,6]=-J21*dn4a+J11*dn4b
    bm[3,7]= J22*dn4a-J12*dn4b
    bm=bm/detJ
    return bm,detJ


def ab_ax4(kk):
    if kk==0:
        a=-0.5773502692
        b=-0.5773502692
    if kk==1:
        a= 0.5773502692
        b=-0.5773502692
    if kk==2:
        a= 0.5773502692
        b= 0.5773502692
    if kk==3:
        a=-0.5773502692
        b= 0.5773502692
    return a,b


def sm_ax4(E,po,z1,r1,z2,r2,z3,r3,z4,r4):
    sm=np.zeros([8,8],dtype=np.float64)
    #Stiffness matrix [sm]=[B]T[D][B]*t*det(J)
    d=dmat_ax4(E,po)
    for kk in range(0,4):
        a,b=ab_ax4(kk)
        bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4)
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4
        sm=sm+np.dot(bm.T,np.dot(d,bm))*rr*detJ
    return sm


def calst_ax4(E,po,alpha,dt1,dt2,dt3,dt4,wd,z1,r1,z2,r2,z3,r3,z4,r4):
    eps0  =np.zeros(4,dtype=np.float64)
    stress=np.zeros(4,dtype=np.float64)
    #stress vector {stress}=[D][B]{u}
    d=dmat_ax4(E,po)
    for kk in range(0,4):
        a,b=ab_ax4(kk)
        bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4)
        eps=np.dot(bm,wd)
        # Strain due to temperature change
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        tem=nm1*dt1 + nm2*dt2 + nm3*dt3 + nm4*dt4
        #Thermal strain
        eps0[0]=alpha*tem
        eps0[1]=alpha*tem
        eps0[2]=alpha*tem
        eps0[3]=0.0
        stress=stress+np.dot(d,(eps-eps0))
    stress=0.25*stress
    return stress


def pst_ax4(sigz,sigr,tauzr):
    ps1=0.5*(sigz+sigr)+np.sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr)
    ps2=0.5*(sigz+sigr)-np.sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr)
    if sigz==sigr:
        if tauzr >0.0: ang= 45.0
        if tauzr <0.0: ang=135.0
        if tauzr==0.0: ang=  0.0
    else:
        ang=0.5*np.arctan(2.0*tauzr/(sigz-sigr))
        ang=180.0/np.pi*ang
        if sigz>sigr and tauzr>=0.0: ang=ang
        if sigz>sigr and tauzr <0.0: ang=ang+180.0
        if sigz<sigr:                ang=ang+90.0
    return ps1,ps2,ang


def tfvec_ax4(E,po,alpha,dt1,dt2,dt3,dt4,z1,r1,z2,r2,z3,r3,z4,r4):
    eps0=np.zeros(4,dtype=np.float64)
    tfe =np.zeros(8,dtype=np.float64)
    # {tfe=[B]T[D]{eps0}
    d=dmat_ax4(E,po)
    for kk in range(0,4):
        a,b=ab_ax4(kk)
        bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4)
        # Strain due to temperature change
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        tem=nm1*dt1 + nm2*dt2 + nm3*dt3 + nm4*dt4
        #Thermal strain
        eps0[0]=alpha*tem
        eps0[1]=alpha*tem
        eps0[2]=alpha*tem
        eps0[3]=0.0
        rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4
        tfe=tfe+np.dot(bm.T,np.dot(d,eps0))*rr*detJ
    return tfe


def bfvec_ax4(gamma,gkz,z1,r1,z2,r2,z3,r3,z4,r4):
    ntn=np.zeros([8,8],dtype=np.float64)
    nm =np.zeros([2,8],dtype=np.float64)
    bfe=np.zeros(8,dtype=np.float64)
    for kk in range(0,4):
        a,b=ab_ax4(kk)
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4
        dn1a=-0.25*(1.0-b)
        dn2a= 0.25*(1.0-b)
        dn3a= 0.25*(1.0+b)
        dn4a=-0.25*(1.0+b)
        dn1b=-0.25*(1.0-a)
        dn2b=-0.25*(1.0+a)
        dn3b= 0.25*(1.0+a)
        dn4b= 0.25*(1.0-a)
        #Jacobi matrix and det(J)
        J11=dn1a*z1 + dn2a*z2 + dn3a*z3 + dn4a*z4
        J12=dn1a*r1 + dn2a*r2 + dn3a*r3 + dn4a*r4
        J21=dn1b*z1 + dn2b*z2 + dn3b*z3 + dn4b*z4
        J22=dn1b*r1 + dn2b*r2 + dn3b*r3 + dn4b*r4
        detJ=J11*J22-J12*J21
        nm[0,0]=nm1
        nm[0,2]=nm2
        nm[0,4]=nm3
        nm[0,6]=nm4
        nm[1,1]=nm[0,0]
        nm[1,3]=nm[0,2]
        nm[1,5]=nm[0,4]
        nm[1,7]=nm[0,6]
        ntn=ntn+np.dot(nm.T,nm)*gamma*rr*detJ
    w=np.array([gkz,0,gkz,0,gkz,0,gkz,0],dtype=np.float64)
    bfe=np.dot(ntn,w)
    return bfe


def main_ax4():
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=4          # Number of nodes per element
    nfree=2        # Degree of freedom per node
    # data input
    npoin,nele,nsec,npfix,nlod,nzdir,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_ax4(fnameR,nod,nfree)
    # print out of input data
    prinp_ax4(fnameW,npoin,nele,nsec,npfix,nlod,nzdir,ae,x,fp,deltaT,mpfix,rdis,node)
    # array declanation
    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 vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        l=node[3,ne]-1
        m=node[4,ne]-1
        z1=x[0,i]; r1=x[1,i]
        z2=x[0,j]; r2=x[1,j]
        z3=x[0,k]; r3=x[1,k]
        z4=x[0,l]; r4=x[1,l]
        E     =ae[0,m] #E    : Elastic modulus
        po    =ae[1,m] #po   : Poisson's ratio
        alpha =ae[2,m] #alpha: Thermal expansion coefficient
        gamma =ae[3,m] #gamma: Unit weight of material
        gkz   =ae[4,m] #gkz  : Acceleration in z-direction
        dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l] # temperature change
        sm=sm_ax4(E,po,z1,r1,z2,r2,z3,r3,z4,r4) # element stiffness matrix
        tfe=tfvec_ax4(E,po,alpha,dt1,dt2,dt3,dt4,z1,r1,z2,r2,z3,r3,z4,r4) # termal load vector
        bfe=bfvec_ax4(gamma,gkz,z1,r1,z2,r2,z3,r3,z4,r4) # body force vector
        sm=sm*float(nzdir)
        tfe=tfe*float(nzdir)
        bfe=bfe*float(nzdir)
        ir[7]=2*l+1; ir[6]=ir[7]-1
        ir[5]=2*k+1; ir[4]=ir[5]-1
        ir[3]=2*j+1; ir[2]=ir[3]-1
        ir[1]=2*i+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]+sm[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 element stress
    wd    =np.zeros(8,dtype=np.float64)
    strs  =np.zeros([4,nele],dtype=np.float64)      # Section force vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        l=node[3,ne]-1
        m=node[4,ne]-1
        z1=x[0,i]; r1=x[1,i]
        z2=x[0,j]; r2=x[1,j]
        z3=x[0,k]; r3=x[1,k]
        z4=x[0,l]; r4=x[1,l]
        wd[0]=disg[2*i]; wd[1]=disg[2*i+1]
        wd[2]=disg[2*j]; wd[3]=disg[2*j+1]
        wd[4]=disg[2*k]; wd[5]=disg[2*k+1]
        wd[6]=disg[2*l]; wd[7]=disg[2*l+1]
        E    =ae[0,m]
        po   =ae[1,m]
        alpha=ae[2,m]
        dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l]
        strs[:,ne]=calst_ax4(E,po,alpha,dt1,dt2,dt3,dt4,wd,z1,r1,z2,r2,z3,r3,z4,r4)
    # print out of result
    prout_ax4(fnameW,npoin,nele,disg,strs)
    # 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()


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

以　上