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.
- A quadrilateral isoparametric element with 1 node and 2 degrees of freedom is used as an axisymmetric element. The rotation direction angle of one element is one radian.
- Initially, a two-dimensional stress analysis program using quadrangular elements ([https://qiita.com/damyarou/items/320bad2052bb5fccd75f]) with the rotation axis set horizontally to the right and the radial axis set vertically upward. I used to rewrite https://qiita.com/damyarou/items/320bad2052bb5fccd75f)), but when "the rotation axis is positive vertically upward and the radial axis is horizontal right", the displacement and stress are Since it is calculated in the opposite direction, we introduced a mechanism to correct the calculation according to the direction of the vertical axis.
- Simultaneous equations are solved by sparse matrix processing so that large-scale analysis can be performed to some extent.
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
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
\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.
(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
Recommended Posts