3D skeleton structure analysis with Python

Big fix

Overview

I made a 3D skeleton structure analysis program in Python that I've always wanted. In the plane frame analysis program created earlier, the rigidity matrix and coordinate transformation matrix are mainly changed to support beam elements with 1 node and 6 degrees of freedom.

~~ The test is not enough yet, but ~~ I was inspired by the article below and posted it. http://qiita.com/sasaco/items/917429a497c6e9a08401

~~ Since the program is created on Jupyter, this article is also described according to the description of my Jupyter. The analysis program by Python is divided into three parts. ~~ I'm always worried, but the file export part is too cluttered. Can't you write better? .. ..

~~ The calculation is going around. .. .. ~~

The programming environment is Machine: MacBook Pro (Retina, 13-inch, Mid 2014) OS: macOS Sierra Python version: Python 3.6.1

Theory of three-dimensional skeleton structure analysis

Basic matters

The theory development and program include the following conditions.

In this program, assuming that the model will be large-scale, sparse matrix processing is introduced when solving simultaneous linear equations.

Element stiffness equation

The general form of the element stiffness equation is the same as in the two-dimensional stress analysis.

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}

Element stiffness matrix of three-dimensional skeleton structure

The element stiffness matrix of the three-dimensional skeleton structure is a 12 x 12 square matrix. This can be derived by extending the one for the plane skeleton, except to consider the twist. The derivation process is omitted, and only the results are described here.

Element stiffness equation (simplest expression)

\begin{equation}
\boldsymbol{\{f\}}=\boldsymbol{[k]}\boldsymbol{\{u\}}
\end{equation}

Nodal displacement vector

\begin{equation}
\boldsymbol{\{u\}}=
\begin{Bmatrix}
u_i & v_i & w_i & \theta_{xi} & \theta_{yi} & \theta_{zi} & u_j & v_j & w_j & \theta_{xj} & \theta_{yj} & \theta_{zj}
\end{Bmatrix}^T
\end{equation}

\begin{align}
&u & &\text{: X-axis displacement} & &\theta_x& &\text{: X-axis rotation amount (torsion angle)}\\
&v & &\text{: Y-axis displacement} & &\theta_y& &\text{: Y-axis rotation amount}\\
&w & &\text{: Z-axis displacement} & &\theta_z& &\text{: Z-axis rotation amount}
\end{align}

Nodal force vector

\begin{equation}
\boldsymbol{\{f\}}=
\begin{Bmatrix}
N_i & S_{yi} & S_{zi} & M_{xi} & M_{yi} & M_{zi} & N_j & S_{yj} & S_{zj} & M_{xj} & M_{yj} & M_{zj}
\end{Bmatrix}^T
\end{equation}

\begin{align}
&N  & &\text{: Axial force (x-axis directional force)} & &M_x& &\text{: X-axis moment (torsion moment)}\\
&S_y& &\text{: Y-axis shear force}   & &M_y& &\text{: Y-axis bending moment}\\
&S_z& &\text{: Z-axis shear force}   & &M_z& &\text{: Z-axis bending moment}
\end{align}

Element stiffness matrix

\begin{equation}
\boldsymbol{[k]}=
\begin{bmatrix}
\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} & 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} \\
0 & 0 & \cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 & 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} \\
-\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} & 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} \\
0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 & 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} \\
\end{bmatrix}
\end{equation}

\begin{align}
&L& &\text{: Element length}        & &J& &\text{: Torsion constant}\\
&E& &\text{: Elastic modulus}      & &I_y& &\text{: Second moment of inertia of area around y-axis}\\
&G& &\text{: Shear modulus}& &I_z& &\text{: Z-axis geometrical moment of inertia}\\
&A& &\text{: Cross-sectional area}        & &
\end{align}

Nodal inertial force vector of plane skeleton element

Consider the inertial force as the physical force, and give "acceleration" as the node vector. Here, the simplest method is to simply distribute the inertial force acting on the element to the nodal force. The total inertial force acting on the element is as follows as the value in the whole coordinate system.

\begin{align}
F_{m(x)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_x \cdot g) \\
F_{m(y)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_y \cdot g) \\
F_{m(z)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_z \cdot g) \\
\end{align}

Therefore, the nodal force $ \ {\ boldsymbol {f_ {b}} } $ due to the inertial force can be determined as follows.

\begin{equation*}
\{\boldsymbol{F_{b}}\}=\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot A \cdot L}{2}
\begin{Bmatrix}
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0   \\
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0
\end{Bmatrix}
\end{equation*}

Here, $ \ gamma $ is the unit volume weight of the element, $ A $ is the element cross-sectional area, $ L $ is the element length, and $ g $ is the gravitational acceleration.

Nodal temperature load vector of plane skeleton elements

It is assumed that the temperature change affects only the axial force of the member. Assuming that the temperature rise is a positive temperature change, the nodal force vector in the member coordinate system due to the temperature change is as follows.

\begin{equation*}
\{\boldsymbol{f_{t}}\}=EA\cdot\alpha\cdot\Delta T
\begin{Bmatrix}
-1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0
\end{Bmatrix}
\end{equation*}

Here, $ EA $ is the axial rigidity of the member, $ \ alpha $ is the coefficient of linear expansion of the member, and $ \ Delta T $ is the amount of temperature change of the member.

By transforming this coordinates, the nodal force vector in the overall coordinate system is determined.

\begin{equation*}
\{\boldsymbol{F_t}\}=[\boldsymbol{T}]^T\{\boldsymbol{f_t}\}
\end{equation*}

Here, $ [\ boldsymbol {T}] ^ T $ is the transposed matrix of the coordinate transformation matrix $ [\ boldsymbol {T}] $ (in this case, equal to the inverse matrix).

Nodal external force vector of plane skeleton element

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. The load is the horizontal force, vertical force, and moment in the overall coordinate system.

Coordinate transformation matrix

The element stiffness equations obtained so far are in the "element coordinate system" fixed to the element. In an actual structure, the elements are oriented in various directions, so it is necessary to find the relationship with the unified "whole coordinate system".

Here, if the coordinate transformation matrix that converts the total coordinate system displacement to the element coordinate system displacement is $ [\ boldsymbol {T}] $, The stiffness matrix $ [\ boldsymbol {K}] $ in the global coordinate system and the stiffness matrix $ [\ boldsymbol {k}] $ in the element coordinate system have the following relationship.

\begin{equation}
[\boldsymbol{K}]=[\boldsymbol{T}]^T[\boldsymbol{k}][\boldsymbol{T}]
\end{equation}

This relationship can be understood from the relationship between the element coordinate system displacement $ \ {\ boldsymbol {u} } $ and the global coordinate system displacement $ \ {\ boldsymbol {U} } $ as follows. ..

\begin{equation}
\{\boldsymbol{u}\}=[\boldsymbol{T}]\{\boldsymbol{U}\} \qquad\qquad \{\boldsymbol{u}\}^T=\{\boldsymbol{U}\}^T [\boldsymbol{T}]^T
\end{equation}

Here, $ [\ boldsymbol {T}] $ is a coordinate conversion matrix that converts the displacement / nodal force of the entire coordinate system into the displacement / nodal force of the element coordinate system.

Specific representation of the coordinate transformation matrix

Direction cosine of member axis (x axis)

Let the nodes that make up the element be $ i $ and $ j $, and let their coordinates be $ (x_i, y_i, z_i) $, $ (x_j, y_j, z_j) $.

\begin{gather}
l=\cfrac{x_j-x_i}{\ell} \qquad m=\cfrac{y_j-y_i}{\ell} \qquad n=\cfrac{z_j-z_i}{\ell} \\
\ell=\sqrt{(x_j-x_i)^2+(y_j-y_i)^2+(z_j-z_i)^2}
\end{gather}

Coordinate transformation submatrix

$ \ boldsymbol {[t]} $ is a submatrix (3 x 3) of the coordinate transformation matrix (12 x 12) that transforms the entire coordinate system into the member coordinate system.

\begin{equation}
\boldsymbol{\{x\}}=\boldsymbol{[t]}\boldsymbol{\{X\}}
\end{equation}

\begin{align}
&\boldsymbol{\{x\}}=\{x,y,z\}^T  & &\text{: Local coordinates}\\
&\boldsymbol{\{X\}}＝\{X,Y,Z\}^T & &\text{: Global coordinate}
\end{align}

\begin{equation}
\boldsymbol{[t]}=
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos\theta & \sin\theta \\
0 & -\sin\theta & \cos\theta
\end{bmatrix}
\cdot
\begin{bmatrix}
l & m & n \\
-\cfrac{m}{\sqrt{l^2+m^2}} & \cfrac{l}{\sqrt{l^2+m^2}} & 0 \\
-\cfrac{l \cdot n}{\sqrt{l^2+m^2}} & -\cfrac{m \cdot n}{\sqrt{l^2+m^2}} & \sqrt{l^2+m^2}
\end{bmatrix}
\qquad \text{($\theta$ : chord angle)}
\end{equation}

\begin{equation}
\boldsymbol{[t]}=
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos\theta & \sin\theta \\
0 & -\sin\theta & \cos\theta
\end{bmatrix}
\cdot
\begin{bmatrix}
0 & 0 & n \\
n & 0 & 0 \\
0 & 1 & 0
\end{bmatrix}
\qquad \text{($\theta$ : chord angle)}
\end{equation}

Coordinate transformation matrix

\begin{equation}
\boldsymbol{[T]}=
\begin{bmatrix}
\boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0} \\\
\boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0} \\\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0} \\\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0}   &\boldsymbol{[t]}
\end{bmatrix}
\end{equation}

Relationship between the member coordinate system (x-y-z) and the overall coordinate system (X-Y-Z) where the code angle is 0

Overall stiffness equation

From the above, the stiffness equation in the overall coordinate system is expressed as follows.

\begin{equation}
[\boldsymbol{K}]\{\boldsymbol{U}\}=\{\boldsymbol{F}\}+\{\boldsymbol{F_b}\}+\{\boldsymbol{F_t}\}
\end{equation}

\begin{align}
&\{\boldsymbol{U}\}                          & & \text{Nodal displacement vector in global coordinate system}\\
&\{\boldsymbol{F}\}                          & & \text{Node external force vector in the global coordinate system}\\
&\{\boldsymbol{F_b}\}=\{\boldsymbol{f_b}\}   & & \text{Nodal inertial force vector in global coordinate system}\\
&\{\boldsymbol{F_t}\}=[\boldsymbol{T}]^T\{\boldsymbol{f_t}\}  & & \text{Nodal temperature load vector in global coordinate system}\\
&[\boldsymbol{K}]=[\boldsymbol{T}]^T[\boldsymbol{k}][\boldsymbol{T}]   & & \text{Rigidity matrix in the global coordinate system}\\
&[\boldsymbol{k}]                            & & \text{Rigidity matrix in element coordinate system}\\
&\{\boldsymbol{f_b}\}                          & & \text{Nodal inertial force vector in element coordinate system}\\
&\{\boldsymbol{f_t}\}                          & & \text{Nodal temperature load vector in element coordinate system}\\
&[\boldsymbol{T}]                            & & \text{Coordinate transformation matrix}
\end{align}

By solving the above equation, the node displacement $ \ {\ boldsymbol {U} } $ in the global coordinate system can be obtained.

Element cross-sectional force

For the element section force, the total coordinate system displacement $ \ {\ boldsymbol {U} } $ obtained above is converted to the element coordinate system displacement by the coordinate conversion matrix $ [\ boldsymbol {T}] $. Calculated by multiplying the element stiffness matrix $ [\ boldsymbol {k}] $. That is, the element section force $ \ {\ boldsymbol {f_ {sec}} } $ is calculated by the following.

\begin{equation}
\{\boldsymbol{f_{sec}}\}=[\boldsymbol{k}]\cdot ([\boldsymbol{T}]\{\boldsymbol{U}\})
\end{equation}

Note that the calculated axial force calculated above does not include the effect of temperature change, so it is necessary to correct the temperature effect of the axial force as follows.

\begin{align}
N_1'=&N_1+EA\cdot\alpha\cdot\Delta T \\
N_2'=&N_2-EA\cdot\alpha\cdot\Delta T
\end{align}

Here, $ N_1'$ and $ N_2'$ are the axial forces at nodes 1 and 2 considering the temperature change, and $ N_1 $ and $ N_2 $ are the axial forces calculated from the displacement as the solution of the simultaneous equations. $ EA $ is the axial stiffness of the member, $ \ alpha $ is the coefficient of linear expansion of the member, and $ \ Delta T $ is the amount of temperature change of the member.

Three-dimensional frame structure analysis program

Program overview

Explanation of boundary conditions that can be handled

Item Description
Nodal external force Specify the number of loading nodes, loading node number, and load value
Nodal temperature change Specify the amount of temperature change at all nodes. Let the temperature rise be the positive value
Element inertial force In the material property data, specify the acceleration in the x / y / z direction as the ratio to the gravitational acceleration
Nodal forced displacement Specify the number of nodes, nodal number, and displacement value to give forced displacement

Program configuration (program name: py_fem_3dfrmx.py)

   1.Module loading
   2.Function: data entry(def inpdata_3dfrm)
       (1)Basic numerical reading
       (2)Array declaration
       (3)Material property data reading
       (4)Read element composition node number and material number
       (5)Read node coordinates and node temperature change
       (6)Boundary condition reading
       (7)Read node load data
   3.Function: Write input data(def prinp_3dfrm)
   4.Function: Export calculation result(def prout_3dfrm)
   5.Function: Create element stiffness matrix(def sm_3dfrm)
   6.Function: Create coordinate conversion matrix(def tm_3dfrm)
   7.Function: Node inertial force vector creation(def bfvec_3dfrm)
   8.Function: Node temperature load vector creation(def tfvec_3dfrm)
   9.Function: Element section force calculation(def calsecf_3dfrm)
  10.Function: Main routine(def main_3dfrm)
       (1)Get I / O file name from command line
       (2)Data entry
       (3)Input data export
       (4)Rigid matrix assembly
       (5)Boundary condition processing (processing of known displacement degree)
       (6)Solving the stiffness equation (displacement calculation)
       (7)Boundary condition processing (incorporation of known displacement)
       (8)Element section force calculation
       (9)Export calculation results
       (10)Information display such as calculation time
  11.Main routine execution

Execution command and I / O data

Analysis execution command

python3 py_fem_3dfrmx.py inp.txt out.txt
Item Description
py_fem_3dfrmx.py Python FEM program name
inp.txt Input data file name (blank delimited data)
out.txt Output data file name (blank separated data)

Input data file format

npoin  nele  nsec  npfix  nlod                   # Basic values for analysis
E  po  A  Ix  Iy  Iz  theta  alpha  gamma  gkX  gkY  gkZ
    ..... (1 to nsec) .....                      # Material properties
node_1  node_2  isec                             # Element connectivity, material set number
    ..... (1 to nele) .....                      #
x  y  z  deltaT                                  # Node coordinate, temperature change of node
    ..... (1 to npoin) .....                     #
lp  kox  koy  koz  kmx  kmy  kmz  rdis_x  rdis_y  rdis_z  rrot_x  rrot_y  rrot_z
    ..... (1 to npfix) .....                     # Restricted node and restrict conditions
lp  fp_x  fp_y  fp_z  mp_x  mp_y  mp_z           # Loaded node and loading conditions
    ..... (1 to nlod) .....                      #     (omit data input if nlod=0)
Item Description
npoin, nele, nsec Number of nodes, number of elements, number of cross-sectional characteristics
npfix, nlod Number of constrained nodes, number of loaded nodes
E, po, A Element elastic modulus, element Poisson's ratio, element cross-sectional area
Ix, Iy, Iz Torsion constant, moment of inertia of area around y-axis, moment of inertia of area around z-axis >
theta, alpha Code angle, coefficient of linear expansion
gamma, gkX, gkY, gkZ Unit volume weight, acceleration in x / y / z direction (ratio of g)
node_1, node_2, isec Node 1, Node 2, Sectional characteristic number
x, y, z, deltaT Node x coordinate, node y coordinate, node z coordinate, node temperature change
lp, kox, koy, koz Constraint node number, x / y / z Directional displacement With or without constraint (constraint: 1, freedom: 0)
kmx, kmy, kmz x / y / z With or without axial rotation constraint (constraint: 1, freedom: 0) < / tr>
r_disx, r_disy, rdis_z x / y / z direction known displacement (enter 0 even if unconstrained)
rrot_x, rrot_y, rrot_z x / y / z Known rotation amount around axis (enter 0 even if unconstrained)
lp, fp_x, fp_y, fp_z Load node number, x-direction load, y-direction load, z-direction load >
mp_x, mp_y, mp_z x-axis moment (torsion), y-axis moment, z-axis moment

Output data file format

npoin  nele  nsec npfix  nlod
    (Each value of above)
sec  E      po     A    J    Iy    Iz    theta
sec  alpha  gamma  gkX  gkY  gkZ
    sec   : Material number
    E     : Elastic modulus
    po    : Poisson's ratio
    A     : Section area
    J     ; Torsional constant
    Iy    : Moment of inertia around y-axis
    Iz    : Moment of inertia around z-axis
    theta : Chord angle
    alpha : Thermal expansion coefficient
    gamma : Unit weight
    gkX   : Acceleration in X-direction (ratio to g)
    gkY   : Acceleration in Y-direction (ratio to g)
    gkZ   : Acceleration in Z-direction (ratio to g)
    ..... (1 to nsec) .....
node   x   y   z   fx   fy   fz   mx   my   mz  deltaT
    node   : Node number
    x      : x-coordinate
    y      : y-coordinate
    fx     : Load in X-direction
    fy     : Load in Y-direction
    fz     : Load in Z-direction
    mx     : Moment load around X-axis
    my     : Moment load around Y-axis
    mz     : Moment load around Z-axis
    deltaT : Temperature change of node
    ..... (1 to npoin) .....
node   kox   koy   koz   kmx   kmy   kmz   rdis_x   rdis_y   rdis_z   rrot_x   rrot_y   rrot_z
    node    : Node number
    kox     : Index of restriction in X-direction (0: free, 1: fixed)
    koy     : Index of restriction in Y-direction (0: free, 1: fixed)
    koz     : Index of restriction in Z-direction (0: free, 1: fixed)
    kmx     : Index of restriction in rotation around X-axis  (0: free, 1: fixed)
    kmy     : Index of restriction in rotation around Y-axis  (0: free, 1: fixed)
    kmz     : Index of restriction in rotation around Z-axis  (0: free, 1: fixed)
    rdis_x  : Given displacement in X-direction
    rdis_y  : Given displacement in Y-direction
    rdis_z  : Given displacement in Z-direction
    rrot_x  : Given displacement in rotation around X-axis
    rrot_y  : Given displacement in rotation around Y-axis
    rrot_z  : Given displacement in rotation around Z-axis
    ..... (1 to npfix) .....
elem     i     j   sec
    elem : Element number
    i    : Node number of start point
    j    : Node number of end point
    sec  : Material number
    ..... (1 to nele) .....
node   dis-x   dis-y   dis-z   rot-x   rot-y   rot-z
    node  : Node number
    dis-x : Displacement in X-direction
    dis-y : Displacement in Y-direction
    dis-z : Displacement in Z-direction
    rot-x : Rotation around X-axis
    rot-y : Rotation around Y-axis
    rot-z : Rotation around Z-axis
    ..... (1 to npoin) .....
elem nodei   N_i   Sy_i   Sz_i   Mx_i   My_i   Mz_i
elem nodej   N_j   Sy_j   Sz_j   Mx_j   My_j   Mz_j
    elem  : Element number
    nodei : First node nymber
    nodej : Second node number
    N_i   : Axial force of node-i
    Sy_i  : Shear force of node-i in y-direction
    Sz_i  : Shear force of node-i in z-direction
    Mx_i  : Moment of node-i around x-axis (torsional moment)
    My_i  : Moment of node-i around y-axis (bending moment)
    Mz_i  : Moment of node-i around z-axis (bending moment)
    N_j   : Axial force of node-j
    Sy_j  : Shear force of node-j in y-direction
    Sz_j  : Shear force of node-j in z-direction
    Mx_j  : Moment of node-j around x-axis (torsional moment)
    My_j  : Moment of node-j around y-axis (bending moment)
    Mz_j  : Moment of node-j around z-axis (bending moment)
    ..... (1 to nele) .....
n=(total degrees of freedom)  time=(calculation time)
Item Description
dis-x, dis-y, dis-z x-direction displacement, y-direction displacement, z-direction displacement
rot-x, rot-y, rot-z x-axis rotation amount, y-axis rotation amount, z-axis rotation amount td>
N, Sy, Sz Axial force, y-direction shear force, z-direction shear force
Mx, My, Mz x-axis moment (torsion), y-axis bending moment, z-axis bending moment

programming

The main body of the program is shown below according to the "program structure".

1. Module loading

# ======================
# 3D Frame Analysis
# ======================
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time

2. Function: Data entry (def inpdata_3dfrm)

def inpdata_3dfrm(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
    # array declaration
    ae    =np.zeros((12,nsec),dtype=np.float64)      # Section characteristics
    node  =np.zeros((nod+1,nele),dtype=np.int)       # node-element relationship
    x     =np.zeros((3,npoin),dtype=np.float64)      # Coordinates of nodes
    deltaT=np.zeros(npoin,dtype=np.float64)          # Temperature change of nodes
    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])  # aa    : section area
        ae[3,i] =float(text[3])  # aix   : tortional constant
        ae[4,i] =float(text[4])  # aiy   : moment of inertia aroutd y-axis
        ae[5,i] =float(text[5])  # aiz   : moment of inertia around z-axis
        ae[6,i] =float(text[6])  # theta : chord angle
        ae[7,i] =float(text[7])  # alpha : thermal expansion coefficient
        ae[8,i] =float(text[8])  # gamma : unit weight of material
        ae[9,i] =float(text[9])  # gkX   : acceleration in X-direction
        ae[10,i]=float(text[10]) # gkY   : acceleration in Y-direction
        ae[11,i]=float(text[11]) # 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]) #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])    # x-coordinate
        x[1,i]=float(text[1])    # y-coordinate
        x[2,i]=float(text[2])    # z-coordinate
        deltaT[i]=float(text[3]) # Temperature change
    # boundary conditions (0:free, 1:restricted)
    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 x-direction
        mpfix[1,lp-1]=int(text[2])   #fixed in y-direction
        mpfix[2,lp-1]=int(text[3])   #fixed in z-direction
        mpfix[3,lp-1]=int(text[4])   #fixed in rotation around x-axis
        mpfix[4,lp-1]=int(text[5])   #fixed in rotation around y-axis
        mpfix[5,lp-1]=int(text[6])   #fixed in rotation around z-axis
        rdis[0,lp-1]=float(text[7])  #fixed displacement in x-direction
        rdis[1,lp-1]=float(text[8])  #fixed displacement in y-direction
        rdis[2,lp-1]=float(text[9])  #fixed displacement in z-direction
        rdis[3,lp-1]=float(text[10]) #fixed rotation around x-axis
        rdis[4,lp-1]=float(text[11]) #fixed rotation around y-axis
        rdis[5,lp-1]=float(text[12]) #fixed rotation around z-axis
    # load
    if 0&lt;nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[6*lp-6]=float(text[1]) #load in x-direction
            fp[6*lp-5]=float(text[2]) #load in y-direction
            fp[6*lp-4]=float(text[3]) #load in z-direction
            fp[6*lp-3]=float(text[4]) #moment around x-axis
            fp[6*lp-2]=float(text[5]) #moment around y-axis
            fp[6*lp-1]=float(text[6]) #moment around z-axis
    f.close()
    return npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp

3. Function: Write input data (def prinp_3dfrm)

def prinp_3dfrm(fnameW,npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s} {4:&gt;5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout)
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
    .format('sec','E','po','A','J','Iy','Iz','theta'),file=fout)
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s}'
    .format('sec','alpha','gamma','gkX','gkY','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} {6:15.7e} {7:15.7e}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i],ae[6,i]),file=fout)
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e}'
        .format(i+1,ae[7,i],ae[8,i],ae[9,i],ae[10,i],ae[11,i]),file=fout)
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s} {8:&gt;15s} {9:&gt;15s} {10:&gt;15s}'
    .format('node','x','y','z','fx','fy','fz','mx','my','mz','deltaT'),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:15.7e} {7:15.7e} {8:15.7e} {9:15.7e} {10:15.7e}'
        .format(lp,x[0,i],x[1,i],x[2,i],fp[6*lp-6],fp[6*lp-5],fp[6*lp-4],fp[6*lp-3],fp[6*lp-2],fp[6*lp-1],deltaT[i]),file=fout)
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s} {4:&gt;5s} {5:&gt;5s} {6:&gt;5s} {7:&gt;15s} {8:&gt;15s} {9:&gt;15s} {10:&gt;15s} {11:&gt;15s} {12:&gt;15s}'
    .format('node','kox','koy','koz','kmx','kmy','kmz','rdis_x','rdis_y','rdis_z','rrot_x','rrot_y','rrot_z'),file=fout)
    for i in range(0,npoin):
        if 0&lt;mpfix[0,i]+mpfix[1,i]+mpfix[2,i]+mpfix[3,i]+mpfix[4,i]+mpfix[5,i]:
            lp=i+1
            print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d} {6:5d} {7:15.7e} {8:15.7e} {9:15.7e} {10:15.7e} {11:15.7e} {12:15.7e}'
            .format(lp,mpfix[0,i],mpfix[1,i],mpfix[2,i],mpfix[3,i],mpfix[4,i],mpfix[5,i],rdis[0,i],rdis[1,i],rdis[2,i],rdis[3,i],rdis[4,i],rdis[5,i]),file=fout)
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s}'.format('elem','i','j','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout)
    fout.close()

4. Function: Export calculation result (def prout_3dfrm)

def prout_3dfrm(fnameW,npoin,nele,node,disg,fsec):
    fout=open(fnameW,'a')
    # displacement
    print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s}'
          .format('node','dis-x','dis-y','dis-z','rot-x','rot-y','rot-z'),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:15.7e}'
              .format(lp,disg[6*lp-6],disg[6*lp-5],disg[6*lp-4],disg[6*lp-3],disg[6*lp-2],disg[6*lp-1]),file=fout)
    # section force
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
    .format('elem','nodei','N_i','Sy_i','Sz_i','Mx_i','My_i','Mz_i'),file=fout)
    print('{0:&gt;5s} {1:&gt;5s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
    .format('elem','nodej','N_j','Sy_j','Sz_j','Mx_j','My_j','Mz_j'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
        .format(ne+1,node[0,ne],fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
        print('{0:5d} {1:5d} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
        .format(ne+1,node[1,ne],fsec[6,ne],fsec[7,ne],fsec[8,ne],fsec[9,ne],fsec[10,ne],fsec[11,ne]),file=fout)
    fout.close()

5. Function: Create element stiffness matrix (def sm_3dfrm)

\begin{equation}
\boldsymbol{[k]}=
\begin{bmatrix}
\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} & 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} \\
0 & 0 & \cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 & 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} \\
-\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\
0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} & 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} \\
0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 & 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 \\
0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 \\
0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 \\
0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} \\
\end{bmatrix}
\end{equation}

def sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2):
    ek=np.zeros((12,12),dtype=np.float64) # local stiffness matrix
    xx=x2-x1
    yy=y2-y1
    zz=z2-z1
    el=np.sqrt(xx**2+yy**2+zz**2)
    ek[ 0, 0]= EA/el
    ek[ 0, 6]=-EA/el
    ek[ 1, 1]= 12*EIz/el**3
    ek[ 1, 5]=  6*EIz/el**2
    ek[ 1, 7]=-12*EIz/el**3
    ek[ 1,11]=  6*EIz/el**2
    ek[ 2, 2]= 12*EIy/el**3
    ek[ 2, 4]= -6*EIy/el**2
    ek[ 2, 8]=-12*EIy/el**3
    ek[ 2,10]= -6*EIy/el**2
    ek[ 3, 3]= GJ/el
    ek[ 3, 9]=-GJ/el
    ek[ 4, 2]= -6*EIy/el**2
    ek[ 4, 4]=  4*EIy/el
    ek[ 4, 8]=  6*EIy/el**2
    ek[ 4,10]=  2*EIy/el
    ek[ 5, 1]=  6*EIz/el**2
    ek[ 5, 5]=  4*EIz/el
    ek[ 5, 7]= -6*EIz/el**2
    ek[ 5,11]=  2*EIz/el
    ek[ 6, 0]=-EA/el
    ek[ 6, 6]= EA/el
    ek[ 7, 1]=-12*EIz/el**3
    ek[ 7, 5]= -6*EIz/el**2
    ek[ 7, 7]= 12*EIz/el**3
    ek[ 7,11]= -6*EIz/el**2
    ek[ 8, 2]=-12*EIy/el**3
    ek[ 8, 4]=  6*EIy/el**2
    ek[ 8, 8]= 12*EIy/el**3
    ek[ 8,10]=  6*EIy/el**2
    ek[ 9, 3]=-GJ/el
    ek[ 9, 9]= GJ/el
    ek[10, 2]= -6*EIy/el**2
    ek[10, 4]=  2*EIy/el
    ek[10, 8]=  6*EIy/el**2
    ek[10,10]=  4*EIy/el
    ek[11, 1]=  6*EIz/el**2
    ek[11, 5]=  2*EIz/el
    ek[11, 7]= -6*EIz/el**2
    ek[11,11]=  4*EIz/el
    return ek

6. Function: Create coordinate transformation matrix (def tm_3dfrm)

\begin{gather}
\ell=\sqrt{(x_j-x_i)^2+(y_j-y_i)^2+(z_j-z_i)^2}\\
l=\cfrac{x_j-x_i}{\ell} \qquad m=\cfrac{y_j-y_i}{\ell} \qquad n=\cfrac{z_j-z_i}{\ell}
\end{gather}

\begin{equation}
\boldsymbol{[t_1]}=
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos\theta & \sin\theta \\
0 & -\sin\theta & \cos\theta
\end{bmatrix}
\qquad \text{($\theta$ : chord angle)}
\end{equation}

\begin{align}
&\boldsymbol{[t_2]}=
\begin{bmatrix}
0 & 0 & n \\
n & 0 & 0 \\
0 & 1 & 0
\end{bmatrix}
& &\text{When the member long axis (member coordinate system x-axis) is parallel to the overall coordinate system Z-axis}\\
&\boldsymbol{[t_2]}=
\begin{bmatrix}
l & m & n \\
-\cfrac{m}{\sqrt{l^2+m^2}} & \cfrac{l}{\sqrt{l^2+m^2}} & 0 \\
-\cfrac{l \cdot n}{\sqrt{l^2+m^2}} & -\cfrac{m \cdot n}{\sqrt{l^2+m^2}} & \sqrt{l^2+m^2}
\end{bmatrix}
& &\text{When the member long axis (member coordinate system x-axis) is not parallel to the overall coordinate system Z-axis}
\end{align}

\begin{equation}
\boldsymbol{[t]}=\boldsymbol{[t_1]}\cdot\boldsymbol{[t_2]}
\end{equation}

\begin{equation}
\boldsymbol{[T]}=
\begin{bmatrix}
\boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0} \\
\boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0}   & \boldsymbol{0} \\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{[t]} & \boldsymbol{0} \\
\boldsymbol{0}   & \boldsymbol{0}   & \boldsymbol{0}   &\boldsymbol{[t]}
\end{bmatrix}
\end{equation}

def tm_3dfrm(theta,x1,y1,z1,x2,y2,z2):
    tt=np.zeros((12,12),dtype=np.float64) # transformation matrix
    t1=np.zeros((3,3),dtype=np.float64)
    t2=np.zeros((3,3),dtype=np.float64)
    xx=x2-x1
    yy=y2-y1
    zz=z2-z1
    el=np.sqrt(xx**2+yy**2+zz**2)
    t1[0,0]=1
    t1[1,1]= np.cos(theta)
    t1[1,2]= np.sin(theta)
    t1[2,1]=-np.sin(theta)
    t1[2,2]= np.cos(theta)
    ll=(x2-x1)/el
    mm=(y2-y1)/el
    nn=(z2-z1)/el
    if x2-x1==0.0 and y2-y1==0.0:
        t2[0,2]=nn
        t2[1,0]=nn
        t2[2,1]=1.0
    else:
        qq=np.sqrt(ll**2+mm**2)
        t2[0,0]=ll
        t2[0,1]=mm
        t2[0,2]=nn
        t2[1,0]=-mm/qq
        t2[1,1]= ll/qq
        t2[2,0]=-ll*nn/qq
        t2[2,1]=-mm*nn/qq
        t2[2,2]=qq
    t3=np.dot(t1,t2)
    tt[0:3,0:3]  =t3[0:3,0:3]
    tt[3:6,3:6]  =t3[0:3,0:3]
    tt[6:9,6:9]  =t3[0:3,0:3]
    tt[9:12,9:12]=t3[0:3,0:3]
    return tt

7. Function: Node inertial force vector creation (def bfvec_3dfrm)

\begin{equation*}
\{\boldsymbol{F_{b}}\}=\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot A \cdot L}{2}
\begin{Bmatrix}
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0   \\
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0
\end{Bmatrix}
\end{equation*}

def bfvec_3dfrm(A,gamma,gkX,gkY,gkZ,x1,y1,z1,x2,y2,z2):
    # Body force vector in global coordinate system
    bfe_g=np.zeros(12,dtype=np.float64)
    xx=x2-x1
    yy=y2-y1
    zz=z2-z1
    el=np.sqrt(xx**2+yy**2+zz**2)
    bfe_g[0]=0.5*gamma*A*el*gkX
    bfe_g[1]=0.5*gamma*A*el*gkY
    bfe_g[2]=0.5*gamma*A*el*gkZ
    bfe_g[6]=0.5*gamma*A*el*gkX
    bfe_g[7]=0.5*gamma*A*el*gkY
    bfe_g[8]=0.5*gamma*A*el*gkZ
    return bfe_g

8. Function: Node temperature load vector creation (def tfvec_3dfrm)

\begin{equation*}
\{\boldsymbol{f_{t}}\}=EA\cdot\alpha\cdot\Delta T
\begin{Bmatrix}
-1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 1 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0
\end{Bmatrix}
\end{equation*}

def tfvec_3dfrm(EA,alpha,tem):
    # Thermal load vector  in local coordinate system
    tfe_l=np.zeros(12,dtype=np.float64)
    tfe_l[0]=-EA*alpha*tem
    tfe_l[6]= EA*alpha*tem
    return tfe_l

9. Function: Element section force calculation (def calsecf_3dfrm)

$ N_1'$ and $ N_2'$ are axial force corrections due to temperature changes.

\begin{align}
&\{\boldsymbol{f_{sec}}\}=[\boldsymbol{k}]\cdot ([\boldsymbol{T}]\{\boldsymbol{U}\})\\
&N_1'=N_1+EA\cdot\alpha\cdot\Delta T \\
&N_2'=N_2-EA\cdot\alpha\cdot\Delta T
\end{align}

def calsecf_3dfrm(EA,GJ,EIy,EIz,theta,alpha,tem,dis,x1,y1,z1,x2,y2,z2):
    ek=sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2) # Stiffness matrix in local coordinate
    tt=tm_3dfrm(theta,x1,y1,z1,x2,y2,z2) # Transformation matrix
    secf=np.dot(ek,np.dot(tt,dis))
    secf[0]=secf[0]+EA*alpha*tem
    secf[6]=secf[6]-EA*alpha*tem
    return secf

10. Function: Main routine

Main routine processing procedure

Item Description
npoin Total number of nodes less
mpfix Array containing constraint conditions for each node (= 1: constraint, = 0: free)
fp Nodal load vector
rdis Displacement amount of the node to be displaced
gk Overall stiffness matrix
disg Solution (displacement) of the overall stiffness equation

  • Solving the stiffness equation (displacement calculation) </ b>: After CSR conversion of the overall stiffness matrix gk </ var>, solve simultaneous linear equations with scipy.sparse.linalg.spsolve () </ var>.

  • Boundary condition processing (incorporation of known displacement) </ b>: After solving the simultaneous linear equations, remember to re-enter the known displacement in the term of degrees of freedom corresponding to the known displacement.

  • Element section force calculation </ b>: Calculation of section force for each element by the function calsecf_3dfrm </ em>

  • Export calculation results </ b>: Calculation result output by function prout_3dfrm </ em>

  • Information display such as calculation time </ b>: By the following, the difference between the calculation completion time and the start time is output as the processing time.

import time

start=time.time()
......
dtime=time.time()-start

Main routine program

def main_3dfrm():
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=2              # Number of nodes per element
    nfree=6            # Degree of freedom per node
    # data input
    npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_3dfrm(fnameR,nod,nfree)
    # print out of input data
    prinp_3dfrm(fnameW,npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp)
    # array declaration
    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 vectors
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]; y1=x[1,i]; z1=x[2,i]
        x2=x[0,j]; y2=x[1,j]; z2=x[2,j]
        ee   =ae[0,m]  # elastic modulus
        po   =ae[1,m]  # Poisson's ratio
        aa   =ae[2,m]  # section area
        aix  =ae[3,m] # tortional constant
        aiy  =ae[4,m] # moment of inertia around y-axis
        aiz  =ae[5,m] # moment of inertia around z-axis
        theta=np.radians(ae[6,m]) # chord angle
        alpha=ae[7,m] # thermal expansion coefficient
        gamma=ae[8,m] # unit weight of material
        gkX  =ae[9,m]   # seismic coefficient in X-direction
        gkY  =ae[10,m]  # seismic coefficient in Y-direction
        gkZ  =ae[11,m]  # seismic coefficient in Z-direction
        A=aa  # section area
        EA=ee*aa
        GJ=ee/2/(1+po)*aix
        EIy=ee*aiy
        EIz=ee*aiz
        tem=0.5*(deltaT[i]+deltaT[j])                # average temperature change
        tt   =tm_3dfrm(theta,x1,y1,z1,x2,y2,z2)         # Transformation matrix
        ek   =sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2) # Stiffness matrix in local coordinate
        ck   =np.dot(np.dot(tt.T,ek),tt)             # Stiffness matrix in global coordinate
        tfe_l=tfvec_3dfrm(EA,alpha,tem)              # Thermal load vector in local coordinate
        tfe  =np.dot(tt.T,tfe_l)                     # Thermal load vector in global coordinate
        bfe  =bfvec_3dfrm(A,gamma,gkX,gkY,gkZ,x1,y1,z1,x2,y2,z2) # Body force vector in global coordinate
        ir[11]=6*j+5; ir[10]=ir[11]-1; ir[9]=ir[10]-1; ir[8]=ir[9]-1; ir[7]=ir[8]-1; ir[6]=ir[7]-1
        ir[5] =6*i+5; ir[4] =ir[5]-1 ; ir[3]=ir[4]-1 ; ir[2]=ir[3]-1; ir[1]=ir[2]-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]+ck[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 section force
    dis=np.zeros(12,dtype=np.float64)
    fsec  =np.zeros((nod*nfree,nele),dtype=np.float64)  # Section force vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]; y1=x[1,i]; z1=x[2,i]
        x2=x[0,j]; y2=x[1,j]; z2=x[2,j]
        ee   =ae[0,m]  # elastic modulus
        po   =ae[1,m]  # Poisson's ratio
        aa   =ae[2,m]  # section area
        aix  =ae[3,m] # tortional constant
        aiy  =ae[4,m] # moment of inertia around y-axis
        aiz  =ae[5,m] # moment of inertia around z-axis
        theta=np.radians(ae[6,m]) # chord angle
        alpha=ae[7,m] # thermal expansion coefficient
        EA=ee*aa
        GJ=ee/2/(1+po)*aix
        EIy=ee*aiy
        EIz=ee*aiz
        tem=0.5*(deltaT[i]+deltaT[j]) # average temperature change
        dis[0]=disg[6*i]  ; dis[1] =disg[6*i+1]; dis[2]= disg[6*i+2]
        dis[3]=disg[6*i+3]; dis[4] =disg[6*i+4]; dis[5]= disg[6*i+5]
        dis[6]=disg[6*j]  ; dis[7] =disg[6*j+1]; dis[8]= disg[6*j+2]
        dis[9]=disg[6*j+3]; dis[10]=disg[6*j+4]; dis[11]=disg[6*j+5]
        fsec[:,ne]=calsecf_3dfrm(EA,GJ,EIy,EIz,theta,alpha,tem,dis,x1,y1,z1,x2,y2,z2)
    # print out of result
    prout_3dfrm(fnameW,npoin,nele,node,disg,fsec)
    # 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()

11. Main routine execution

#==============
# Execution
#==============
if __name__ == '__main__': main_3dfrm()

that's all

Recommended Posts