2D FEM stress analysis program with Python

I think the post below is better for the finish as a program, so please refer to it.

Explanation of the theory of 2D stress analysis program (quadrilateral element) by Python (revised edition)
Finite element method and Python program (two-dimensional stress analysis program using triangular primary elements)

Overview

I made a linear elastic two-dimensional FEM stress analysis program by Python.

The element is an isoparametric element with 4 nodes per element, and each element has 4 Gauss integral points.

The loads that can be handled are as follows.

Nodal external force
Nodal temperature change
Self-weight and inertial force during an earthquake

Theoretical formula

Nodal displacement calculation formula

[\boldsymbol{k}]\\{\boldsymbol{u}\\}=\\{\boldsymbol{f}\\}+\\{\boldsymbol{f_t}\\}+\\{\boldsymbol{f_b}\\}

\begin{align} &[\boldsymbol{k}]=t\cdot \int_A [\boldsymbol{B}]^T[\boldsymbol{D}][\boldsymbol{B}] dA \\\ &\\{\boldsymbol{f_t}\\}=t\cdot \int_A [\boldsymbol{B}]^T[\boldsymbol{D}]\\{\boldsymbol{\epsilon_0}\\} dA \\\ &\\{\boldsymbol{f_b}\\}=t\cdot\gamma\cdot \int_A [\boldsymbol{N}]^T [\boldsymbol{N}] dA \cdot \\{\boldsymbol{w}\\} \end{align}

$ [\ boldsymbol {k}] $: Element stiffness matrix
$ \ {\ boldsymbol {u} \} $: Nodal displacement vector
$ \ {\ boldsymbol {f} \} $: Nodal external force vector
$ \ {\ boldsymbol {f_t} \} $: Nodal force vector due to temperature change
$ \ {\ boldsymbol {f_b} \} $: Node vector due to own weight and seismic acceleration
$ [\ boldsymbol {D}] $: Stress-strain relationship matrix
$ t $, $ \ gamma $: Element thickness, element unit volume weight
$ \ {\ boldsymbol {w} \} $: Nodal acceleration vector (input as a ratio to gravitational acceleration)

Element stress calculation formula

\\{\boldsymbol{\sigma}\\}=[\boldsymbol{D}]\\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\\}

$ \ {\ boldsymbol {\ epsilon_0} \} $ is distortion due to temperature change

Nodal displacement-arbitrary position strain relational expression

\\{\boldsymbol{\epsilon}\\}=[\boldsymbol{B}]\\{\boldsymbol{u}\\}

$ \ {\ boldsymbol {\ epsilon} \} $ is the distortion of any point in the element

Nodal displacement-Arbitrary position strain calculation formula

\\{\boldsymbol{v}\\}=[\boldsymbol{N}]\\{\boldsymbol{u}\\}

$ \ {\ boldsymbol {v} \} $ is the displacement of any point in the element

Can it be used?

Since Python is an interpreter, I was worried whether it would be useful in terms of calculation speed, but the calculation time for the model with the number of nodes: 3417, the number of elements: 3257, and the degree of freedom: 6834 is 7.9 seconds, which is a usable level. It was confirmed that.

program

The program is long, so I put a link to what I pasted in Gist. In addition, contour drawing only fills the elements with different colors according to the stress level. For my purpose, it is not necessary to insert it separately and draw a contour line.

Input data format

npoin  nele  nsec  npfix  nlod  NSTR #Basic amount
t  E  po  alpha  gamma  gkh  gkv     #Material property
..... (1~nsec) .....
node1  node2  node3  node4  isec     #element-Nodal relationship, material property number
..... (1~nele) .....
x  y  deltaT                         #Nodal coordinates, nodal temperature changes
..... (1~npoin) .....
node  kox  koy  rdisx  rdisy         #Displacement constraint conditions
..... (1~npfix) .....
node  fx  fy                         #External force
..... (1~nlod) .....

npoin, nele, nsec	Number of nodes, number of elements, number of material characteristics
npfix, nlod	Number of constrained nodes, number of loaded nodes
NSTR	Stress state (plane stress: 0, plane stress: 1)
t, E, po, alpha	Plate thickness, elastic modulus, Poisson's ratio, coefficient of linear expansion
gamma, gkh, gkv	Unit volume weight, horizontal and vertical acceleration (ratio of g)
x, y, deltaT	Node x coordinate, node y coordinate, node temperature change
node, kox, koy	Constrain node number, presence / absence of x and y direction constraints (constraint: 1, freedom: 0)
rdisx, rdisy	Displacement in x and y directions (enter 0 even if unconstrained)
node, fx, fy	Load node number, x-direction load, y-direction load

Output data format

npoin  nele  nsec npfix  nlod   NSTR
    4     1     1     2     2      1
  sec               t               E              po           alpha           gamma        gkh        gkv
    1   1.0000000e+00   1.0000000e+03   0.0000000e+00   1.0000000e-05   2.3000000e+00      0.000      0.000
 node               x               y              fx              fy          deltaT   kox   koy
    1   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     1     1
    2   1.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     0     1
    3   1.0000000e+00   1.0000000e+00   0.0000000e+00   1.0000000e+01   0.0000000e+00     0     0
    4   0.0000000e+00   1.0000000e+00   0.0000000e+00   1.0000000e+01   0.0000000e+00     0     0
 node   kox   koy          rdis_x          rdis_y
    1     1     1   0.0000000e+00   0.0000000e+00
    2     0     1   0.0000000e+00   0.0000000e+00
 elem     i     j     k     l   sec
    1     1     2     3     4     1
 node           dis-x           dis-y
    1   0.0000000e+00   0.0000000e+00
    2   8.8817842e-19   0.0000000e+00
    3   1.7763568e-18   2.0000000e-02
    4   4.1448326e-18   2.0000000e-02
 elem           sig_x           sig_y          tau_xy              p1              p2             ang
    1  -7.4014868e-16   2.0000000e+01   1.1333527e-16   2.0000000e+01   0.0000000e+00   9.0000000e+01
n=8  time=0.044 sec

node, dis-x, dis-y	Node number, displacement in x direction, displacement in y direction
elme, sig_x, sig_y, tau_xy	Element number, x-direction direct stress, y-direction direct stress, shear stress
p1, p2, ang	First principal stress, second principal stress, direction of first principal stress
n, time	total degrees of freedom, calculation time

Output example

that's all