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()
以 上