Previous article was up to the program for creating the element stiffness matrix, so this time I will explain the program for creating the overall stiffness matrix from the element stiffness matrix. .. Similar to the previous time, it is based on The first finite element method by Java published by Corona Publishing.
Last time, I made an element stiffness matrix and an element stiffness equation for one element of an object. We will formulate this for each element in the object. If an element i has nodes j, k, then the element stiffness matrix and the element stiffness equations are:
{\bf K}^{(i)} \left(
\begin{array}{c}
u_j \\
u_k
\end{array}
\right)
=
\left(
\begin{array}{c}
K_{11}^{(i)}\ K_{12}^{(i)} \\
K_{21}^{(i)}\ K_{22}^{(i)}
\end{array}
\right)
\left(
\begin{array}{c}
u_j \\
u_k
\end{array}
\right)
=
\left(
\begin{array}{c}
f^{(i)}_j \\
f^{(i)}_k
\end{array}
\right)
(1)
Can be written as. For example, if there are three elements and the i-th element has two nodes (i, i + 1), the element stiffness equation is as follows.
\left(
\begin{array}{c}
K_{11}^{(1)}\ K_{12}^{(1)} \\
K_{21}^{(1)}\ K_{22}^{(1)}
\end{array}
\right)
\left(
\begin{array}{c}
u_1 \\
u_2
\end{array}
\right)
=
\left(
\begin{array}{c}
f^{(1)}_1 \\
f^{(1)}_2
\end{array}
\right)
(2)
\left(
\begin{array}{c}
K_{11}^{(2)}\ K_{12}^{(2)} \\
K_{21}^{(2)}\ K_{22}^{(2)}
\end{array}
\right)
\left(
\begin{array}{c}
u_2 \\
u_3
\end{array}
\right)
=
\left(
\begin{array}{c}
f^{(2)}_2 \\
f^{(2)}_3
\end{array}
\right)
(3)
\left(
\begin{array}{c}
K_{11}^{(3)}\ K_{12}^{(3)} \\
K_{21}^{(3)}\ K_{22}^{(3)}
\end{array}
\right)
\left(
\begin{array}{c}
u_3 \\
u_4
\end{array}
\right)
=
\left(
\begin{array}{c}
f^{(3)}_3 \\
f^{(3)}_4
\end{array}
\right)
(4)
Considering that each element shares each node, these equations can be combined into one,
\left(
\begin{array}{cccc}
K_{11}^{(1)}& K_{12}^{(1)} & 0 & 0 \\
K_{21}^{(1)}& K_{22}^{(1)} + K_{11}^{(2)} & K_{12}^{(2)} & 0 \\
0 & K_{21}^{(2)} & K_{22}^{(2)} + K_{11}^{(3)} & K_{12}^{(3)} \\
0 & 0 & K_{21}^{(3)} & K_{22}^{(3)} \\
\end{array}
\right)
\left(
\begin{array}{c}
u_1 \\
u_2 \\
u_3 \\
u_4
\end{array}
\right)
=
\left(
\begin{array}{c}
f^{(1)}_1\\
f^{(1)}_2+f^{(2)}_2 \\
f^{(2)}_3+f^{(3)}_3 \\
f^{(3)}_4
\end{array}
\right)
(5)
It will be. The sum of the element stiffness equations in this way is called the overall stiffness equation, and the coefficient matrix on the left side is called the overall stiffness matrix. The total stress and strain are calculated by solving the overall stiffness equation in consideration of external forces and boundary conditions.
When the node coordinates, element area / Young's modulus, force applied to the node, and node displacement constraint condition of each node are input, the program that outputs the corresponding overall stiffness equation is as follows. (The force applied to the nodes and the node displacement constraints of each node are not used this time, but they will be necessary parameters when solving the next overall stiffness equation.) In the program, the Element class created last time is used. I will.
import numpy as np
#Data entry and calculation preparation
#list_x :Node coordinates
#list_area :Area of each element
#list_young :Young's modulus of each element
#list_force :Force applied to each node
#list_nodeCondition :Nodal displacement constraint condition of each node
def __init__(self,list_x,list_area,list_young,list_force,list_nodeCondition):
#Setting the number of nodes NumberOfNode, nnode
self.NumberOfNode = int(list_x.shape[0])
nnode = self.NumberOfNode
#Setting the number of elements NumberOfElement, nelm
self.NumberOfElement = nnode-1
nelm = self.NumberOfElement
#Create lists and arrays for data storage, calculation, and output
self.x = np.zeros(nnode) #Nodal coordinates
self.force = np.zeros(nnode) #Force on the node
self.nodeCondition = np.zeros((nnode,2)) #Nodal displacement constraint condition
self.totalStiffness = np.zeros((nnode,nnode))#Overall stiffness matrix
self.disp = np.zeros(nnode) #Nodal displacement
self.stress = np.zeros(nelm) #stress
#List for storing Element classes
self.elem = np.array([])
#Initial condition storage of node coordinates, force applied to nodes, and node displacement conditions
self.x = list_x
self.force = list_force
self.nodeCondition = list_nodeCondition
#Create Element class and store in elem list
#Then store the data about the element. Implemented for the number of elements.
for i in range(nelm):
self.elem = np.append(self.elem, Element())
self.elem[i].young=list_young[i]
self.elem[i].area=list_area[i]
for j in range(2):
#Stores the contact numbers ix and coordinates xe that make up each element
self.elem[i].ix[j]=i+j
self.elem[i].xe[j]=list_x[i+j]
#Fabrication of total stiffness matrix
def assembleMatrix(self):
#Bmatrix / element stiffness matrix creation for each element
print("\nAssemble Stiffness Matrix : \n")
nelm = self.NumberOfElement
for i in range(nelm):
self.elem[i].makeBmatrix()
self.elem[i].makeElementStiffness()
#Create an overall stiffness matrix from each element stiffness matrix
for k in range(nelm):
print(" Element",k)
for i in range(2):
for j in range(2):
ii = int(self.elem[k].ix[i])
jj = int(self.elem[k].ix[j])
self.totalStiffness[ii,jj] += self.elem[k].Kmatrix[i,j]
The above operation example is as follows. Here, it is assumed that there are four nodes, each coordinate is (50, 150, 250, 350) (mm), the cross-sectional area of each element is 100 mm $ ^ 2 $, and Young's modulus is 200 GPa. Also, although not used this time, list_nodeCondition specifies that the node coordinates are fixed at 50 mm, and list_force specifies that a force of 1000 N can be applied to the node coordinates of 350 mm.
input
#node parameters
list_x = np.array([50,150,250,350])
list_force = np.array([0.0,0.0,0.0,1000.])
list_nodeCondition = np.array([[1,0],[0,0],[0,0],[0,0]])
#Element parameters
list_area = np.array([100,100,100])
list_young = np.array([200000,200000,200000])
Fem = Fem1dim(list_x,list_area,list_young,list_force,list_nodeCondition)
Fem.assembleMatrix()
print(Fem.totalStiffness)
output
Assemble Stiffness Matrix :
Element 0
Element 1
Element 2
[[ 200000. -200000. 0. 0.]
[-200000. 400000. -200000. 0.]
[ 0. -200000. 400000. -200000.]
[ 0. 0. -200000. 200000.]]
This time, we explained the formula of the overall stiffness equation from the element stiffness equation and introduced the overall stiffness matrix fabrication program. Next time, I will introduce a program that applies displacement boundary conditions and calculates the overall stiffness equation (overall stiffness matrix).
Shigeru Nagagi (2010) "First finite element method by Java" Corona Publishing Co., Ltd.