Programme d'analyse des contraintes symétriques axiales par Python (élément carré) [édition révisée]

Publié le 15/06/2020

introduction

Présentation d'un programme d'analyse des contraintes à symétrie axiale utilisant Python basé sur la théorie de l'élasticité des microdéformations. Les caractéristiques du programme sont les suivantes.

Équation de rigidité des éléments

Equation pour trouver le déplacement

\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{(Matrice de rigidité des éléments)}\\
&\{\boldsymbol{f}\}=\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA  & & \text{(Vecteur de force externe nodale)}\\
&\{\boldsymbol{f_t}\}=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV  & & \text{(Vecteur de distorsion du nœud initial)}\\
&\{\boldsymbol{f_b}\}=\gamma\left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} & & \text{(Vecteur de force d'inertie nodale)}
\end{align}

Formule de calcul des contraintes internes de l'élément

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

\begin{align}
&\{\boldsymbol{\epsilon}\}  & &\text{: Vecteur de déformation de l'élément obtenu à partir du déplacement obtenu en résolvant l'équation de rigidité}\\
&\{\boldsymbol{\epsilon_0}\} & &\text{: Déformation initiale telle que déformation thermique}\\
&[\boldsymbol{D_e}]          & &\text{: Matrice de relation contrainte-déformation}
\end{align}

Forme de base de l'équation de rigidité des éléments

fig

Considérons un élément rectangulaire sur le plan du système de coordonnées polaires orthogonales $ (z-r) $ comme indiqué dans la figure ci-dessus. Ici, $ z $ est l'axe de rotation, la direction horizontale droite est positive, $ r $ est l'axe radial et la direction verticale ascendante est positive. De plus, l'angle de direction de rotation d'un élément est réglé sur un angle unitaire, c'est-à-dire un radian, afin de simplifier la formulation et le calcul d'intégration. Comme pour l'analyse bidimensionnelle des contraintes, l'élément carré pour l'analyse de symétrie axiale doit avoir 4 nœuds par élément et 2 degrés de liberté (déplacement de direction $ z $ et déplacement de direction $ r $). Par conséquent, l'équation de rigidité de l'élément de base prend la forme suivante. Ici, laissez les nombres de nœuds qui composent l'élément quadrangulaire être $ i $, $ j $, $ k $ et $ l $ dans le sens antihoraire sur le plan $ (z-r) $.

\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{:nœud$i$de$z$Charge directionnelle} & & w_i & &\text{:nœud$i$de$z$Déplacement directionnel}\\
&f_{r,i} & &\text{:nœud$i$de$r$Charge directionnelle} & & u_i & &\text{:nœud$i$de$r$Déplacement directionnel}\\
&f_{z,j} & &\text{:nœud$j$de$z$Charge directionnelle} & & w_i & &\text{:nœud$j$de$z$Déplacement directionnel}\\
&f_{r,j} & &\text{:nœud$j$de$r$Charge directionnelle} & & u_i & &\text{:nœud$j$de$r$Déplacement directionnel}\\
&f_{z,k} & &\text{:nœud$k$de$z$Charge directionnelle} & & w_i & &\text{:nœud$k$de$z$Déplacement directionnel}\\
&f_{r,k} & &\text{:nœud$k$de$r$Charge directionnelle} & & u_i & &\text{:nœud$k$de$r$Déplacement directionnel}\\
&f_{z,l} & &\text{:nœud$l$de$z$Charge directionnelle} & & w_l & &\text{:nœud$l$de$z$Déplacement directionnel}\\
&f_{r,l} & &\text{:nœud$l$de$r$Charge directionnelle} & & u_l & &\text{:nœud$l$de$r$Déplacement directionnel}
\end{align}

Formule de base par théorie de l'élasticité linéaire

Relation contrainte-déformation et relation déformation-déplacement dans un problème de symétrie axiale

Si la contrainte / déformation dans l'axe de rotation est représentée par l'indice $ z $, la contrainte / déformation dans la direction radiale est représentée par l'indice $ r $, et la contrainte / déformation dans la direction circonférentielle est représentée par l'indice $ \ theta $, la contrainte dans le problème de symétrie axiale- L'expression relationnelle de déformation est la suivante.

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

Ici, $ \ nu $ est le coefficient de Poisson, $ \ alpha $ est le taux de dilatation thermique et $ T $ est la quantité de changement de température.

En transformant l'équation ci-dessus

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

Exprimé sous la forme de

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

L'expression relationnelle déformation-déplacement est exprimée comme suit en utilisant le déplacement $ w $ dans l'axe de rotation et le déplacement $ u $ dans la direction radiale.

\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 à l'aide d'éléments quadrangulaires

Introduction de la matrice de relation déformation-déplacement

Le déplacement d'axe de rotation $ w $ et le déplacement radial $ u $ de tout point de l'élément axialement symétrique sont supposés comme suit. Où les coordonnées ($ a $, $ b ) sont des coordonnées de paramètres normalisés avec une plage de [ -1 $, $ 1 $], $ w_ {i, j, k, l} $ et $ u_ {i, j , k, l} $ indique le déplacement des nœuds qui composent l'élément.

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

Le déplacement de tout point de l'élément $ \ {w, u \} $ est exprimé comme suit en utilisant le déplacement de nœud $ \ {\ 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}

La distorsion $ \ {\ boldsymbol {\ epsilon} \} $ en tout point de l'élément utilise le déplacement de nœud $ \ {\ 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}

Ici, la différenciation liée à la fonction d'interpolation $ N_i $ peut être exprimée comme suit en utilisant la matrice de 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}

D'après ce qui précède, les éléments de la matrice de relation déformation-déplacement $ [\ boldsymbol {B}] $ peuvent être exprimés comme suit.

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

Notez qu'à ce stade, $ [\ boldsymbol {B}] $ est défini en fonction des variables $ a $ et $ b $.

Les éléments de $ [J] $ sont indiqués ci-dessous. Ici, $ z_ {i, j, k, l} $ et $ r_ {i, j, k, l} $ sont les coordonnées des nœuds qui composent l'élément.

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

La valeur de coordonnée $ r $ de tout point de l'élément est calculée comme suit.

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

Ici, $ r_i $, $ r_j $, $ r_k $ et $ r_l $ sont les valeurs de coordonnées $ r $ des nœuds qui composent l'élément.

Intégration gaussienne

L'intégration de surface pour les éléments requis pour la matrice de rigidité d'élément ultérieure et le calcul du vecteur de charge est effectuée par l'intégration gaussienne indiquée ci-dessous.

En supposant que les points d'intégration sont de 4 points ($ n = 2 $), les valeurs de $ a $, $ b $ et le poids $ H $ sont comme indiqué dans le tableau ci-dessous, et $ a $ et $ b $ sont modifiés 4 La somme des valeurs des temps est la valeur approximative de l'intégration. De plus, comme le poids est $ H = 1 $, le calcul devient encore plus simple.

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

Ici, dans le cas d'un élément axialement symétrique, le volume d'un élément minuscule situé à un rayon de $ r $ du centre de rotation est $ r \ fois d \ theta \ fois dr \ fois dz $. Si la valeur intégrée, le vecteur de charge, etc. sont évalués avec 1 radian dans la direction circonférentielle et que la conversion de la variable intégrale pour l'utilisation de l'intégrale de Gauss est effectuée,

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

Ici, $ r_i $, $ r_j $, $ r_k $ et $ r_l $ sont les valeurs de coordonnées $ r $ des nœuds qui composent l'élément.

Matrice de rigidité des éléments

La matrice de rigidité $ [\ boldsymbol {k}] $ de l'élément paramétrique à 4 nœuds est la suivante, où la matrice de relation contrainte-déformation est $ [\ boldsymbol {D_e}] $. Étant donné que l'angle de direction de rotation d'un élément est d'un radian, il n'apparaît pas dans le calcul de l'intégrale.

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

Vecteur de charge de température de nœud

Considérez la déformation thermique comme la déformation initiale. Le vecteur de force nodale $ \ {\ boldsymbol {f_t} } $ dû à la déformation $ \ {\ boldsymbol {\ epsilon_0} } $ due au changement de température est le suivant.

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

Les composants de la déformation due aux changements de température sont les suivants.

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

Ici, $ \ alpha $ est le coefficient de dilatation linéaire, $ T_p $ est la quantité de changement de température en tout point d'élément, et $ \ {\ phi_i ~~ \ phi_j ~~ \ phi_k ~~ \ phi_l } ^ T $ est le changement de température du nœud. La quantité, $ [N_1 ~~ N_2 ~~ N_3 ~~ N_4] $, est un élément de $ [\ boldsymbol {N}] $. Le signe de la quantité de changement de température est que l'élévation de température est positive.

Vecteur de force d'inertie de nœud

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

Ici, $ \ gamma $ est le poids volumique unitaire de l'élément (masse $ \ fois $ accélération de gravité), et $ k_z $ est l'accélération de l'axe de rotation (rapport à l'accélération de gravité). Puisqu'il semble qu'il n'est pas pratique d'appliquer une force d'inertie radiale, elle est ici traitée comme nulle.

Vecteur de force externe nodale

En ce qui concerne le vecteur de force externe du nœud, dans le programme présenté ici, la force externe du nœud équivalent est directement entrée à partir du fichier de données d'entrée, et aucun traitement de calcul spécial n'est effectué.

Dans l'analyse de symétrie axiale, il faut faire attention à la méthode d'entrée de la charge. Un exemple du concept de charge d'entrée est présenté ci-dessous.

(Exemple du concept de charge d'entrée)

Lorsqu'une pression interne de p = 1N / mm $ ^ 2 $ (= 1MPa) agit sur un élément cylindrique avec un rayon de r = 2000mm et une longueur axiale de z = 500mm, la charge totale agissant sur la paroi interne de l'élément est la suivante. Il devient.

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

Par conséquent, sur la base du concept de force nodale équivalente, 500 000 (N) sont définis au nœud de (z, r) = (0,2 000) et 500 000 (N) sont définis au nœud de (z, r) = (500,2 000). Il doit être chargé dans le sens.

Vecteur de contrainte d'élément

Le vecteur de contrainte d'élément $ \ {\ boldsymbol {\ sigma} } $ peut être calculé comme suit.

\begin{align}
&\{\boldsymbol{\epsilon}\}=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\} & &\text{(Déformation obtenue à partir du déplacement du nœud)}\\
&T_p=N_1\cdot\phi_i+N_2\cdot\phi_j+N_3\cdot\phi_k+N_4\cdot\phi_l & &\text{(Quantité de changement de température en tout point de l'élément)}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_p \\ \alpha\cdot T_p \\ \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{(Déformation de température)}\\
&\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\} & &\text{(Vecteur de contrainte d'élément)}
\end{align}

Ici, $ N_1 \ sim N_4 $ est une fonction d'insertion, $ \ phi_i $, $ \ phi_j $, $ \ phi_k $, $ \ phi_l $ sont des nœuds d'élément $ i $, $ j $, $ k $, $ l C'est la quantité de changement de température de $.

Gestion du système de coordonnées dans le programme

Un diagramme en image du système de coordonnées et des éléments est présenté ci-dessous. Dans (Cas-1) et (Cas-2), l'axe de rotation est l'axe z et l'axe radial est l'axe r.

fig_coord.jpg

(Reference) La figure montrée dans (Référence) est les coordonnées pour une analyse de contrainte bidimensionnelle normale. L'axe des x est défini horizontalement vers la droite et l'axe des y est défini verticalement vers le haut. Dans ce cas, le numéro d'élément définit généralement l'élément en désignant le numéro de nœud dans le sens antihoraire.

(Case-1) Dans (Cas-1), l'axe z (axe de rotation) est défini horizontalement vers la droite et l'axe r (direction radiale) est défini verticalement vers le haut. Dans ce cas également, comme dans l'analyse de contrainte bidimensionnelle de (Référence), le numéro d'élément définit l'élément en désignant le numéro de nœud dans le sens antihoraire. (Le programme est conçu pour ce faire)

(Case-2) Dans (Cas-2), l'axe z (axe de rotation) est défini verticalement vers le haut et l'axe r (direction radiale) est défini horizontalement vers la droite. Dans ce cas, si l'élément est défini par le numéro de nœud dans le sens antihoraire comme dans l'analyse de contrainte bidimensionnelle de (Référence), l'axe z et l'axe r sont interchangés par rapport à (Cas-1), donc le numéro de nœud de l'élément. Est défini dans la direction opposée, et la matrice ou le vecteur à calculer à l'aide des coordonnées d'élément telles que la matrice de rigidité, le vecteur de force d'objet et le vecteur de charge de température est calculé avec le signe opposé à celui habituel.

contre-mesure

Par conséquent, dans le programme d'analyse des contraintes à symétrie axiale préparé cette fois, la variable nzdir </ code> est introduite pour inverser les directions de la matrice de rigidité, du vecteur de force objet et du vecteur de charge de température en fonction de la direction des coordonnées. fait. C'est-à-dire, dans la figure ci-dessus (Cas-1), nzdir = 1 </ code>, et dans la figure ci-dessus (Cas-2), nzdir = -1 </ code>, qui a été obtenu sur le programme principal. Les signes de la matrice de rigidité, du vecteur de force de l'objet et du vecteur de charge de température sont ajustés.

        sm=sm*float(nzdir) #Matrice de rigidité
        tfe=tfe*float(nzdir) #Vecteur de charge de température
        bfe=bfe*float(nzdir) #Vecteur de force de l'objet (force d'inertie)

Format des données d'entrée

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 Nombre de nœuds, nombre d'éléments, nombre de propriétés du matériau
npfix, nlod Nombre de nœuds de contrainte, nombre de nœuds de chargement
nzdir direction de l'axe z (horizontal vers la droite: 1, vertical vers le haut: -1)
E, po, alpha Coefficient d'élasticité, coefficient de Poisson, coefficient de dilatation linéaire
gamma, gkz Poids volumique unitaire, accélération dans la direction z (rapport de g)
nœud1, nœud2, nœud3, nœud4, isec Nœud 1, Nœud 2, Nœud 3, Nœud 4, Numéro de propriété du matériau
z, r, delta T Coordonnée du nœud z, coordonnée du nœud r, changement de température du nœud
nœud, koz, kor Contrainte numéro de nœud, direction z et r Contrainte (contrainte: 1, liberté: 0)
rdisz, rdisr déplacement de direction z et r (entrez 0 même sans contrainte)
nœud, fz, fr Numéro de nœud de charge, charge z-direction, charge r-direction

(Note) Que la direction de l'axe de rotation z soit horizontale (nzdir = 1) ou verticale (nzdir = -1), la saisie des coordonnées est effectuée dans l'ordre de (z, r).

Programme complet

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

c'est tout