Axialsymmetrisches Spannungsanalyseprogramm von Python (quadratisches Element) [überarbeitete Ausgabe]

Gepostet am 2020/06/15

Einführung

Einführung eines axialsymmetrischen Spannungsanalyseprogramms mit Python basierend auf der Theorie der Mikrodeformationselastizität. Die Funktionen des Programms sind wie folgt.

• Ein quadratisches isoparametrisches Element mit 1 Knoten und 2 Freiheitsgraden wird als axialsymmetrisches Element verwendet. Der Drehrichtungswinkel eines Elements beträgt ein Bogenmaß.
• Zunächst ein zweidimensionales Spannungsanalyseprogramm mit quadratischen Elementen ([https://qiita.com/damyarou/items/320bad2052bb5fccd75f]), bei dem die Rotationsachse horizontal nach rechts positiv und die Radialachse vertikal nach oben positiv eingestellt ist. Früher habe ich https://qiita.com/damyarou/items/320bad2052bb5fccd75f)) umgeschrieben, aber wenn "die Rotationsachse vertikal nach oben positiv und die Radialachse horizontal nach rechts positiv ist", sind Verschiebung und Spannung gleich Da es in die entgegengesetzte Richtung berechnet wird, haben wir einen Mechanismus eingeführt, um die Berechnung entsprechend der Richtung der vertikalen Achse zu korrigieren.
• Gleichzeitige Gleichungen werden durch spärliche Matrixverarbeitung gelöst, so dass bis zu einem gewissen Grad eine Analyse in großem Maßstab durchgeführt werden kann.

Elementsteifigkeitsgleichung

Gleichung zum Finden der Verschiebung

\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{(Elementsteifigkeitsmatrix)}\\
&\{\boldsymbol{f}\}=\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA  & & \text{(Nodaler externer Kraftvektor)}\\
&\{\boldsymbol{f_t}\}=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV  & & \text{(Anfänglicher Knotenverzerrungsvektor)}\\
&\{\boldsymbol{f_b}\}=\gamma\left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} & & \text{(Knotenträgheitskraftvektor)}
\end{align}

Berechnungsformel für die interne Spannung des Elements

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

\begin{align}
&\{\boldsymbol{\epsilon}\}  & &\text{: Dehnungsvektor des Elements, erhalten aus der Verschiebung, die durch Lösen der Steifigkeitsgleichung erhalten wurde}\\
&\{\boldsymbol{\epsilon_0}\} & &\text{: Anfangsbelastung wie Temperaturbelastung}\\
&[\boldsymbol{D_e}]          & &\text{: Spannungs-Dehnungs-Beziehungsmatrix}
\end{align}

Grundform der Elementsteifigkeitsgleichung

Betrachten Sie ein rechteckiges Element in der Ebene des orthogonalen Polarkoordinatensystems $ (z-r) $, wie in der obigen Abbildung gezeigt. Hier ist $ z $ die Rotationsachse, die horizontale rechte Richtung ist positiv, $ r $ ist die radiale Achse und die vertikale Aufwärtsrichtung ist positiv. Ferner wird der Drehrichtungswinkel eines Elements auf einen Einheitswinkel, dh einen Bogenmaß, eingestellt, um die Formulierung und die Integralberechnung zu vereinfachen. Ähnlich wie bei der zweidimensionalen Spannungsanalyse muss das quadratische Element für die axiale Symmetrieanalyse 4 Knoten pro Element und 2 Freiheitsgrade ($ z $ Richtungsverschiebung und $ r $ Richtungsverschiebung) aufweisen. Daher hat die Grundelementsteifigkeitsgleichung die folgende Form. Hier seien die Knotennummern, aus denen das viereckige Element besteht, $ i $, $ j $, $ k $ und $ l $ gegen den Uhrzeigersinn in der $ (z-r) $ -Ebene.

\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{:Knoten$i$von$z$Richtungslast} & & w_i & &\text{:Knoten$i$von$z$Richtungsverschiebung}\\
&f_{r,i} & &\text{:Knoten$i$von$r$Richtungslast} & & u_i & &\text{:Knoten$i$von$r$Richtungsverschiebung}\\
&f_{z,j} & &\text{:Knoten$j$von$z$Richtungslast} & & w_i & &\text{:Knoten$j$von$z$Richtungsverschiebung}\\
&f_{r,j} & &\text{:Knoten$j$von$r$Richtungslast} & & u_i & &\text{:Knoten$j$von$r$Richtungsverschiebung}\\
&f_{z,k} & &\text{:Knoten$k$von$z$Richtungslast} & & w_i & &\text{:Knoten$k$von$z$Richtungsverschiebung}\\
&f_{r,k} & &\text{:Knoten$k$von$r$Richtungslast} & & u_i & &\text{:Knoten$k$von$r$Richtungsverschiebung}\\
&f_{z,l} & &\text{:Knoten$l$von$z$Richtungslast} & & w_l & &\text{:Knoten$l$von$z$Richtungsverschiebung}\\
&f_{r,l} & &\text{:Knoten$l$von$r$Richtungslast} & & u_l & &\text{:Knoten$l$von$r$Richtungsverschiebung}
\end{align}

Grundformel nach linearer Elastizitätstheorie

Spannungs-Dehnungs-Beziehung und Dehnungs-Verschiebungs-Beziehung bei axialer Symmetrie

Wenn die Spannung / Dehnung in der Rotationsachse durch den Index $ z $ dargestellt wird, wird die Spannung / Dehnung in radialer Richtung durch den Index $ r $ dargestellt, und die Spannung / Dehnung in Umfangsrichtung wird durch den Index $ \ theta $ dargestellt, die Spannung in der axialen Symmetrie problem- Der Stamm relationale Expression ist wie folgt.

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

Hier ist $ \ nu $ das Poisson-Verhältnis, $ \ alpha $ die Wärmeausdehnungsrate und $ T $ das Ausmaß der Temperaturänderung.

Durch Transformation der obigen Gleichung

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

Ausgedrückt in Form von

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

Der relationale Ausdruck von Dehnung und Verschiebung wird wie folgt ausgedrückt, wobei die Verschiebung $ w $ in der Rotationsachse und die Verschiebung $ u $ in radialer Richtung verwendet werden.

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

Formulierung mit viereckigen Elementen

Einführung der Dehnungs-Verschiebungs-Beziehungsmatrix

Die Rotationsachsenverschiebung $ w $ und die radiale Verschiebung $ u $ eines beliebigen Punktes in dem axialsymmetrischen Element werden wie folgt angenommen. Wobei die Koordinaten ($ a $, $ b ) normalisierte Parameterkoordinaten mit einem Bereich von [ -1 $, $ 1 $], $ w_ {i, j, k, l} $ und $ u_ {i, j sind , k, l} $ gibt die Verschiebung der Knoten an, aus denen das Element besteht.

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

Die Verschiebung eines beliebigen Punktes im Element $ \ {w, u \} $ wird unter Verwendung der Knotenverschiebung $ \ {\ boldsymbol {u_ {nd}} \} $ wie folgt ausgedrückt. ，

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

Die Verzerrung $ \ {\ boldsymbol {\ epsilon} \} $ an einem beliebigen Punkt im Element verwendet die Knotenverschiebung $ \ {\ 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}

Hier kann die Differenzierung in Bezug auf die Interpolationsfunktion $ N_i $ unter Verwendung der Jacobi-Matrix $ [J] $ wie folgt ausgedrückt werden.

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

Aus dem Obigen können die Elemente der Dehnungs-Verschiebungs-Beziehungsmatrix $ [\ boldsymbol {B}] $ wie folgt ausgedrückt werden.

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

Beachten Sie, dass zu diesem Zeitpunkt $ [\ boldsymbol {B}] $ als Funktion der Variablen $ a $ und $ b $ definiert ist.

Die Elemente von $ [J] $ sind unten gezeigt. Hier sind $ z_ {i, j, k, l} $ und $ r_ {i, j, k, l} $ die Knotenkoordinaten, aus denen das Element besteht.

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

Der $ r $ -Koordinatenwert eines beliebigen Punkts im Element wird wie folgt berechnet.

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

Hier sind $ r_i $, $ r_j $, $ r_k $ und $ r_l $ die $ r $ -Koordinatenwerte der Knoten, aus denen das Element besteht.

Gaußsche Integration

Die Flächenintegration für die Elemente, die für die nachfolgende Berechnung der Elementsteifigkeitsmatrix und des Lastvektors erforderlich sind, wird durch die unten gezeigte Gaußsche Integration durchgeführt.

Unter der Annahme, dass die Integrationspunkte 4 Punkte sind ($ n = 2 $), sind die Werte von $ a $, $ b $ und das Gewicht $ H $ wie in der folgenden Tabelle gezeigt, und $ a $ und $ b $ werden geändert 4 Die Summe der Zeitwerte ist der ungefähre Wert der Integration. Da das Gewicht $ H = 1 $ ist, wird die Berechnung noch einfacher.

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

Hier beträgt im Fall eines axialsymmetrischen Elements das Volumen eines winzigen Elements, das sich in einem Radius von $ r $ vom Rotationszentrum befindet, $ r \ mal d \ theta \ mal dr \ mal dz $. Wenn der integrierte Wert, der Lastvektor usw. mit 1 Bogenmaß in Umfangsrichtung ausgewertet werden und die Umwandlung der Integralvariablen zur Verwendung des Gauß-Integrals durchgeführt wird,

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

Hier sind $ r_i $, $ r_j $, $ r_k $ und $ r_l $ die $ r $ -Koordinatenwerte der Knoten, aus denen das Element besteht.

Elementsteifigkeitsmatrix

Die Steifheitsmatrix $ [\ boldsymbol {k}] $ des parametrischen 4-Knoten-Elements ist wie folgt, wobei die Spannungs-Dehnungs-Beziehungsmatrix $ [\ boldsymbol {D_e}] $ ist. Da der Drehrichtungswinkel eines Elements ein Bogenmaß ist, erscheint er nicht in der Integralberechnung.

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

Knotentemperaturlastvektor

Betrachten Sie die Temperaturbelastung als anfängliche Belastung. Der Knotenkraftvektor $ \ {\ boldsymbol {f_t} } $ aufgrund der Dehnung $ \ {\ boldsymbol {\ epsilon_0} } $ aufgrund der Temperaturänderung ist wie folgt.

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

Die Komponenten der Dehnung aufgrund von Temperaturänderungen sind wie folgt.

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

Hier ist $ \ alpha $ der lineare Expansionskoeffizient, $ T_p $ ist der Betrag der Temperaturänderung an einem beliebigen Elementpunkt und $ \ {\ phi_i ~~ \ phi_j ~~ \ phi_k ~~ \ phi_l } ^ T $ ist die Änderung der Knotentemperatur. Die Menge $ [N_1 ~~ N_2 ~~ N_3 ~~ N_4] $ ist ein Element von $ [\ boldsymbol {N}] $. Das Vorzeichen für das Ausmaß der Temperaturänderung ist, dass der Temperaturanstieg positiv ist.

Trägheitskraftvektor des Knotens

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

Hier ist $ \ gamma $ das Einheitsvolumengewicht des Elements (Masse $ \ mal $ Schwerkraftbeschleunigung) und $ k_z $ ist die Beschleunigung der Rotationsachse (Verhältnis zur Schwerkraftbeschleunigung). Da es anscheinend nicht praktikabel ist, eine radiale Trägheitskraft aufzubringen, wird sie hier als Null behandelt.

Knoten externer Kraftvektor

In Bezug auf den externen Kraftvektor des Knotens wird in dem hier vorgestellten Programm die äquivalente externe Kraft des Knotens direkt aus der Eingabedatendatei eingegeben, und es wird keine spezielle Berechnungsverarbeitung durchgeführt.

Bei der Axialsymmetrieanalyse ist bei der Lasteingabemethode Vorsicht geboten. Ein Beispiel für das Konzept der Eingangslast ist unten dargestellt.

(Beispiel für das Konzept der Eingangslast)

Wenn ein Innendruck von p = 1 N / mm $ ^ 2 $ (= 1 MPa) auf ein zylindrisches Element mit einem Radius von r = 2000 mm und einer axialen Länge von z = 500 mm wirkt, ist die auf die Innenwand des Elements wirkende Gesamtlast wie folgt. Es wird.

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

Basierend auf dem Konzept der äquivalenten Knotenkraft werden daher 500.000 (N) am Knoten von (z, r) = (0,2.000) und 500.000 (N) am Knoten von (z, r) = (500,2.000) festgelegt. Es sollte in die Richtung geladen werden.

Elementspannungsvektor

Der Elementspannungsvektor $ \ {\ boldsymbol {\ sigma} } $ kann wie folgt berechnet werden.

\begin{align}
&\{\boldsymbol{\epsilon}\}=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\} & &\text{(Dehnung aus Knotenverschiebung)}\\
&T_p=N_1\cdot\phi_i+N_2\cdot\phi_j+N_3\cdot\phi_k+N_4\cdot\phi_l & &\text{(Ausmaß der Temperaturänderung an einem beliebigen Punkt im Element)}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_p \\ \alpha\cdot T_p \\ \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{(Temperaturbelastung)}\\
&\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\} & &\text{(Elementspannungsvektor)}
\end{align}

Hier ist $ N_1 \ sim N_4 $ eine Einfügefunktion, $ \ phi_i $, $ \ phi_j $, $ \ phi_k $, $ \ phi_l $ sind Elementknoten $ i $, $ j $, $ k $, $ l Dies ist der Betrag der Temperaturänderung von $.

Handhabung des Koordinatensystems im Programm

Ein Bilddiagramm des Koordinatensystems und der Elemente ist unten gezeigt. In (Fall-1) und (Fall-2) ist die Rotationsachse die Z-Achse und die Radialachse die R-Achse.

(Reference) Die in (Referenz) gezeigte Abbildung ist die Koordinate für die normale zweidimensionale Spannungsanalyse. Die x-Achse ist horizontal nach rechts und die y-Achse vertikal nach oben eingestellt. In diesem Fall definiert die Elementnummer normalerweise das Element, indem die Knotennummer gegen den Uhrzeigersinn angegeben wird.

(Case-1) In (Fall 1) ist die z-Achse (Rotationsachse) horizontal nach rechts und die r-Achse (radiale Richtung) vertikal nach oben eingestellt. Auch in diesem Fall wie in der zweidimensionalen Spannungsanalyse von (Referenz) definiert die Elementnummer das Element, indem die Knotennummer gegen den Uhrzeigersinn angegeben wird. (Das Programm ist dafür ausgelegt)

(Case-2) In (Fall 2) ist die z-Achse (Rotationsachse) vertikal nach oben und die r-Achse (radiale Richtung) horizontal nach rechts eingestellt. In diesem Fall werden, wenn das Element durch die Knotennummer gegen den Uhrzeigersinn wie in der zweidimensionalen Spannungsanalyse von (Referenz) definiert ist, die z-Achse und die r-Achse im Vergleich zu (Fall-1) vertauscht, also die Knotennummer des Elements. Wird in entgegengesetzter Richtung definiert, und die Matrix oder der Vektor wird unter Verwendung der Elementkoordinaten wie der Steifheitsmatrix, des Objektkraftvektors und des Temperaturlastvektors mit dem entgegengesetzten Vorzeichen von dem üblichen berechnet.

Gegenmaßnahme

Daher wird in dem diesmal erstellten axialsymmetrischen Spannungsanalyseprogramm die Variable nzdir </ code> eingeführt, um die Richtungen der Steifigkeitsmatrix, des Objektkraftvektors und des Temperaturlastvektors gemäß der Richtung der Koordinaten umzukehren. tat. Das heißt, in der obigen Abbildung (Fall 1) ist nzdir = 1 </ code> und in der obigen Abbildung (Fall 2) nzdir = -1 </ code>, die im Hauptprogramm erhalten wurde. Die Vorzeichen der Steifigkeitsmatrix, des Objektkraftvektors und des Temperaturlastvektors werden angepasst.

        sm=sm*float(nzdir) #Steifheitsmatrix
        tfe=tfe*float(nzdir) #Temperaturlastvektor
        bfe=bfe*float(nzdir) #Objektkraftvektor (Trägheitskraft)

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)

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