Axisymmetric stress analysis program by Python (square element) [revised edition]

Posted on 2020/06/15

Introduction

Introducing an axisymmetric stress analysis program using Python based on the theory of micro-deformation elasticity. The features of the program are as follows.

Element stiffness equation

Equation for finding displacement

\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{(Element stiffness matrix)}\\
&\{\boldsymbol{f}\}=\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA  & & \text{(Nodal external force vector)}\\
&\{\boldsymbol{f_t}\}=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV  & & \text{(Initial node distortion vector)}\\
&\{\boldsymbol{f_b}\}=\gamma\left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} & & \text{(Nodal inertial force vector)}
\end{align}

Element internal stress calculation formula

\begin{equation}
\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\}
\end{equation}
\begin{align}
&\{\boldsymbol{\epsilon}\}  & &\text{: Strain vector of the element obtained from the displacement obtained by solving the stiffness equation}\\
&\{\boldsymbol{\epsilon_0}\} & &\text{: Initial strain such as temperature strain}\\
&[\boldsymbol{D_e}]          & &\text{: Stress-strain relation matrix}
\end{align}

Basic form of element stiffness equation

fig

Consider a quadrilateral element on the orthogonal polar coordinate system $ (z-r) $ plane as shown in the above figure. Here, $ z $ is the axis of rotation, the horizontal right direction is positive, $ r $ is the radial axis, and the vertical upward direction is positive. Further, the rotation direction angle of one element is set to a unit angle, that is, one radian, and the formulation / integral calculation is simplified. As with the two-dimensional stress analysis, the quadrilateral element for axisymmetric analysis shall have four nodes per element and two degrees of freedom (displacement in the $ z $ direction and displacement in the $ r $ direction). Therefore, the basic element stiffness equation takes the following form. Here, let the node numbers that make up the quadrilateral element be $ i $, $ j $, $ k $, and $ l $ counterclockwise on the $ (z-r) $ plane.

\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{:node$i$of$z$Directional load} & & w_i & &\text{:node$i$of$z$Directional displacement}\\
&f_{r,i} & &\text{:node$i$of$r$Directional load} & & u_i & &\text{:node$i$of$r$Directional displacement}\\
&f_{z,j} & &\text{:node$j$of$z$Directional load} & & w_i & &\text{:node$j$of$z$Directional displacement}\\
&f_{r,j} & &\text{:node$j$of$r$Directional load} & & u_i & &\text{:node$j$of$r$Directional displacement}\\
&f_{z,k} & &\text{:node$k$of$z$Directional load} & & w_i & &\text{:node$k$of$z$Directional displacement}\\
&f_{r,k} & &\text{:node$k$of$r$Directional load} & & u_i & &\text{:node$k$of$r$Directional displacement}\\
&f_{z,l} & &\text{:node$l$of$z$Directional load} & & w_l & &\text{:node$l$of$z$Directional displacement}\\
&f_{r,l} & &\text{:node$l$of$r$Directional load} & & u_l & &\text{:node$l$of$r$Directional displacement}
\end{align}

Basic formula based on the theory of linear elasticity

Stress-strain relationship and strain-displacement relationship in axisymmetric problem

If the stress / strain in the rotation axis direction is represented by the subscript $ z $, the stress / strain in the radial direction is represented by the subscript $ r $, and the stress / strain in the circumferential direction is represented by the subscript $ \ theta $, the stress in the axial symmetry problem- The strain relation equation is as follows.

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

Here, $ \ nu $ is the Poisson's ratio, $ \ alpha $ is the coefficient of thermal expansion, and $ T $ is the amount of temperature change.

By transforming the above equation

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

Expressed in the form of

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

The strain-displacement relational expression is expressed as follows using the displacement $ w $ in the axis of rotation and the displacement $ u $ in the radial direction.

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

Formulation using quadrilateral elements

Introduction of strain-displacement relationship matrix

The rotational axial displacement $ w $ and the radial displacement $ u $ of any point in the axisymmetric element are assumed as follows. Here the coordinates ($ a $, $ b ) are the normalized parameter coordinates with the range [ -1 $, $ 1 $], $ w_ {i, j, k, l} $ and $ u_ {i, j , k, l} $ indicates the displacement of the nodes that make up the element.

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

The displacement of any point in the element $ \ {w, u \} $ is expressed as follows using the node displacement $ \ {\ 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}

Strain $ \ {\ boldsymbol {\ epsilon} \} $ at any point in the element uses node displacement $ \ {\ 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}

Here, the derivative for the interpolation function $ N_i $ can be expressed as follows using the Jacobian matrix $ [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}

From the above, the elements of the strain-displacement relationship matrix $ [\ boldsymbol {B}] $ can be expressed as follows.

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

Note that at this stage $ [\ boldsymbol {B}] $ is defined as a function of the variables $ a $ and $ b $.

The elements of $ [J] $ are shown below. Here, $ z_ {i, j, k, l} $ and $ r_ {i, j, k, l} $ are the nodal coordinates that make up the element.

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

The $ r $ coordinate value of any point in the element is calculated by the following.

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

Here, $ r_i $, $ r_j $, $ r_k $, and $ r_l $ are the $ r $ coordinate values of the nodes that make up the element.

Gauss integral

The surface integrals for the elements required for the subsequent element stiffness matrix and load vector calculation are performed by the Gaussian integrals shown below.

Assuming that the integration points are 4 points ($ n = 2 $), the values of $ a $, $ b $ and the weight $ H $ are as shown in the table below, and $ a $ and $ b $ are changed 4 The sum of the values of the times is an approximation of the integral. Also, since the weight is $ H = 1 $, the calculation becomes even simpler.

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

Here, in the case of an axisymmetric element, the volume of the minute element located at a radius of $ r $ from the center of rotation is $ r \ times d \ theta \ times dr \ times dz $. If the integral value, load vector, etc. are evaluated in 1 radian in the circumferential direction, and the integral variable transformation for using the Gauss integral is performed,

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

Here, $ r_i $, $ r_j $, $ r_k $, and $ r_l $ are the $ r $ coordinate values of the nodes that make up the element.

Element stiffness matrix

The stiffness matrix $ [\ boldsymbol {k}] $ of the 4-node parametric element is as follows, where the stress-strain relationship matrix is $ [\ boldsymbol {D_e}] $. Since the rotation direction angle of one element is one radian, it does not appear in the integral calculation.

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

Nodal temperature load vector

Consider temperature strain as the initial strain. The nodal force vector $ \ {\ boldsymbol {f_t} } $ due to the strain $ \ {\ boldsymbol {\ epsilon_0} } $ due to temperature change is as follows.

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

The components of strain due to temperature changes are as follows.

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

Here, $ \ alpha $ is the coefficient of linear expansion, $ T_p $ is the amount of temperature change at any element point, and $ \ {\ phi_i ~~ \ phi_j ~~ \ phi_k ~~ \ phi_l } ^ T $ is the node temperature change. The quantity, $ [N_1 ~~ N_2 ~~ N_3 ~~ N_4] $, is an element of $ [\ boldsymbol {N}] $. The sign of the amount of temperature change is that the temperature rise is positive.

Nodal inertial force vector

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

Here, $ \ gamma $ is the unit volume weight of the element (mass $ \ times $ gravitational acceleration), and $ k_z $ is the rotational axis acceleration (ratio to gravitational acceleration). Since it seems that it is not practical to apply a radial inertial force, it is treated as zero here.

Nodal external force vector

Regarding the nodal external force vector, in the program introduced here, the equivalent nodal external force is directly input from the input data file, and no special calculation processing is performed.

In the axisymmetric analysis, care must be taken in the load input method. An example of the concept of input load is shown below.

(Example of the concept of input load)

When an internal pressure of p = 1N / mm $ ^ 2 $ (= 1MPa) acts on one element of a cylinder with radius r = 2000mm and axial length z = 500mm, the total load acting on the inner wall of the element is as follows. It becomes.

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

Therefore, based on the concept of equivalent nodal force, 500,000 (N) is set at the node of (z, r) = (0,2,000) and 500,000 (N) is set at the node of (z, r) = (500,2,000). It should be loaded in the direction.

Element stress vector

The element stress vector $ \ {\ boldsymbol {\ sigma} } $ can be calculated as follows.

\begin{align}
&\{\boldsymbol{\epsilon}\}=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\} & &\text{(Strain obtained from node displacement)}\\
&T_p=N_1\cdot\phi_i+N_2\cdot\phi_j+N_3\cdot\phi_k+N_4\cdot\phi_l & &\text{(Amount of temperature change at any point in the element)}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_p \\ \alpha\cdot T_p \\ \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{(Temperature strain)}\\
&\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\} & &\text{(Element stress vector)}
\end{align}

Here, $ N_1 \ sim N_4 $ is the interpolation function, $ \ phi_i $, $ \ phi_j $, $ \ phi_k $, $ \ phi_l $ are the element nodes $ i $, $ j $, $ k $, $ l. The amount of temperature change in $.

Programmatic handling of coordinate system

An image diagram of the coordinate system and elements is shown below. In (Case-1) and (Case-2), the rotation axis is the z-axis and the radial axis is the r-axis.

fig_coord.jpg

(Reference) The figure shown in (Reference) is the coordinates for normal two-dimensional stress analysis. The x-axis is set horizontally to the right and the y-axis is set to vertically upward. In this case, the element number usually defines the element by designating the node number counterclockwise.

(Case-1) In (Case-1), the z-axis (rotational axis) is set horizontally to the right and the r-axis (radial direction) is set to vertically upward. In this case as well, as in the two-dimensional stress analysis of (Reference), the element number defines the element by designating the node number counterclockwise. (The program is designed to do so)

(Case-2) In (Case-2), the z-axis (rotational axis) is set vertically upward, and the r-axis (radial direction) is set horizontally to the right. In this case, if the element is defined by the counterclockwise node number as in the two-dimensional stress analysis of (Reference), the z-axis and r-axis are interchanged compared to (Case-1), so the node number of the element. Is defined in the opposite direction, and the matrix or vector calculated using the element coordinates such as the rigidity matrix, the object force vector, and the temperature load vector is calculated with the opposite sign from the usual one.

counter-measure

Therefore, in the axisymmetric stress analysis program prepared this time, the variable nzdir </ code> is introduced to reverse the directions of the stiffness matrix, object force vector, and temperature load vector according to the direction of the coordinates. did. That is, in the above figure (Case-1), nzdir = 1 </ code> is set, and in the above figure (Case-2), nzdir = -1 </ code> is set, and it is obtained on the main program. The signs of the stiffness matrix, physical strength vector, and temperature load vector are adjusted.

        sm=sm*float(nzdir) #Rigidity matrix
        tfe=tfe*float(nzdir) #Temperature load vector
        bfe=bfe*float(nzdir) #Object force (inertial force) vector

Input data format

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)
Variables Description
npoin, nele, nsec Number of nodes, number of elements, number of material characteristics
npfix, nlod Number of constrained nodes, number of loaded nodes
nzdir z-axis direction (horizontal right direction: 1, vertical upward direction: -1)
E, po, alpha Elastic modulus, Poisson's ratio, coefficient of linear expansion
gamma, gkz Unit volume weight, z-direction acceleration (ratio of g)
node1, node2, node3, node4, isec Node 1, node 2, node 3, node 4, material property number
z, r, delta T Nodal z coordinate, nodal r coordinate, nodal temperature change
node, koz, kor Constrain node number, z and r direction Constrained (constraint: 1, freedom: 0)
rdisz, rdisr z and r displacement (enter 0 even if unconstrained)
node, fz, fr Load node number, z-direction load, r-direction load

(Note) Coordinate input is performed in the order of (z, r) regardless of whether the direction of the rotation axis z is horizontal (nzdir = 1) or vertical (nzdir = -1).

Full program

################################
# 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()

that's all