Finite element method beginning (one-dimensional element stiffness matrix)

About this article

This article is written in Java with reference to The first finite element method by Java published by Corona Publishing Co., Ltd. It is a program made with python. In addition, this time it is up to the program to create the element stiffness matrix.

Outline of the finite element method

The finite element method is used for simulating the stress and strain of solid mechanics. In solid mechanics, an object is divided into many small regions (elements), and the state of each element is calculated to determine the state of deformation of the entire substance. The state of an object is expressed by three differential equations (stress parallel equation, strain-displacement relational equation, stress-strain relational equation (construction equation)), but in the finite element method, these differential equations are expressed by equivalent integral expressions ( By substituting with the virtual work principle), the displacement and strain can be calculated approximately. The integral representation replaces the differential equations with the stiff equations, which are simultaneous linear equations. In addition, the calculation flow of the finite element method is to divide the object into each element, establish a stiffness equation (element stiffness matrix) for each element, and then integrate the stiffness equations of all the elements to obtain the overall stiffness equation (overall). It's like assembling a stiffness matrix) and calculating an approximate solution.

Basic differential equation

For simplicity, one endpoint of an object with a cross-sectional area of $ A \ \ rm {[L ^ {2}]} $ and a length of $ l \ \ rm {[L]} $ is fixed and the other endpoint. Consider the case where a one-dimensional and uniform tensile load $ F \ \ rm {[F]} $ works on. (Here, let us assume that $ \ rm {[L]} $ represents the dimension of length and $ \ rm [F] $ represents the dimension of force. Using such a notation, for example, pressure is a unit area. Since it is a hitting force, it can be written as the dimension of force divided by the square of the dimension of length, that is, $ \ rm {[FL ^ {-2}]} $. The dimension of the parameter is clarified. Doing makes it easier to see what that parameter represents and how it relates to other parameters.) In addition, this substance is made into an elastic body (a substance that deforms when a force is applied and returns to its original shape when the force is stopped), and its Young's modulus is $ E \ \ rm {[FL ^ {-2}]}. Suppose it is $. Young's modulus is the proportional coefficient of strain $ \ varepsilon $ [dimensionless] of Hooke's law that "stress is proportional to strain" that holds in elastic bodies. (Similar to the relationship between load (or stress) and displacement when a force is applied to a spring, but be aware that the dimensions are different.)

Assuming these, the balance of forces in this object can be written as follows, assuming that the stress generated in the object is $ \ sigma \ \ rm {[FL ^ {-2}]} $.

$ F = \ sigma A $ (force balance, stress equilibrium equation) (1')

Furthermore, if the stress is $ \ sigma $, Hooke's law can be written as follows using Young's modulus $ E $ and strain $ \ varepsilon $.

$ \ sigma = E \ \ varepsilon $ (material properties, stress-strain equation) (2')

The strain $ \ varepsilon $ can be expressed as follows, where the elongation of the object under a tensile load is $ u \ rm {[L]} $.

$ \ Varepsilon = u / L $ (shape change, strain-displacement relation) (3')

These equations (1')-(3') are special cases of basic differential equations. Considering the more general situation, the following three basic differential equations are derived.

$ d (\ sigma A) / dx = 0 $ (force balance, stress equilibrium equation) (1) $ \ sigma = E \ \ varepsilon $ (material properties, stress-strain equation) (2) $ \ Varepsilon = du / dx $ (shape change, strain-displacement relational expression) (3)

When calculating the strain generated in an elastic body analytically, it is necessary to solve the differential equations of equations (1)-(3).

Virtual work principle

The principle of virtual work is explained as follows. Please refer to the references as the detailed explanation will be long.

When an object is deformed and in a balanced state, the virtual work done by the external force is equal to the virtual work done by the internal force.

The principle of virtual work can be expressed by a mathematical formula as follows.

Arbitrary virtual displacement and virtual strain with balanced force         \int_V \sigma \delta\varepsilon dV = F\delta u Is established.

Where $ \ delta \ varepsilon $ is the virtual strain and $ \ delta u $ is the virtual displacement.

Finite element method for one dimension (up to the fabrication of element stiffness matrix)

Division of elements, the principle of virtual work

Now, the finite element method is ready. In the finite element method, an object is first divided into small areas, "elements". Then, apply the virtual work principle to each element and create an element stiffness equation (element stiffness matrix) for each element. Let's consider the principle of virtual work for one element when a one-dimensional object is divided. Since an element has two nodes (joints between elements), one is node 1 (coordinate $ x_1 $) and the other is node 2 (coordinate $ x_2 $), and the internal force acting on each is $ f_1. Let the displacement of the points of $ and $ f_2 $ be $ u_1 $ and $ u_2 $, respectively. If the stress here is $ \ sigma $ and the virtual displacements of each node are $ \ delta u_1 $ and $ \ delta u_2 $, the principle of virtual work of this element is

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 $  (4)

Can be written.

Linear approximation of displacement

Now suppose that inside each element, the displacement $ u $ changes with a linear function of coordinates $ x $.

u(x) = ax+b  (5)

This displacement must be equal to $ u_1 $, $ u_2 $ at nodes 1 and 2, so

u(x_1) = ax_1+b =u_1  (6)         u(x_1) = ax_2+b =u_2  (7)

Summarizing this for $ a $ and $ b $,

a=\frac{u_2 - u_1}{x_2 - x_1} = \frac{u_2 - u_1}{l}  (8)

b=\frac{x_2 u_1 - x_1 u_2}{x_2 - x_1}=\frac{x_2 u_1 - x_1 u_2}{l}  (9)

It will be. Where $ l = x_2-x_1 $ is the length of the element. Using (8) and (9), (5) is

u(x) = ax+b = \frac{u_2 - u_1}{l}x+\frac{x_2 u_1 - x_1 u_2}{l}  (10)

Furthermore, by rewriting equation (10),

u(x) = \frac{x_2 - x}{l}u_1+\frac{x - x_1}{l} u_2 =N_1(x)u_1 + N_2(x)u_2  (11)

Therefore, we were able to rewrite the equation with the nodal displacements $ u_1 $ and $ u_2 $ as unknowns. Also, $ N_1 (x) $ and $ N_2 (x) $ are called shape functions. In addition, equation (11) can be expressed in a matrix,

$u(x) = {\bf N}^{{\rm T}} {\bf U} $  (12)

Can be written. Where $ {\ bf N} $ and $ {\ bf U} $ are


  {\bf N} = \left(
    \begin{array}{c}
      N_1(x) \\
      N_2(x) 
    \end{array}
  \right),

  {\bf U} = \left(
    \begin{array}{c}
      u_1 \\
      u_2
    \end{array}
  \right)  (13)

Is the vector.

Approximate display of strain and stress

Using the displacement equation (11) and the strain-displacement equation (3), the strain $ \ varepsilon $ is

\varepsilon = \frac{d}{dx}(N_1(x)u_1 + N_2(x)u_2) = \frac{dN_1(x)}{dx}u_1 + \frac{dN_2(x)}{dx}u_2 (14)

Can be written. This is also a matrix notation,

$\varepsilon = {\bf B}^{{\rm T}} {\bf U} $  (15)

Can be rewritten as. Here, $ {\ bf B} $ is called the $ {\ bf B} $ matrix, or strain-displacement matrix, and can be written as follows.

  {\bf B} = \left(
    \begin{array}{c}
      \frac{dN_1(x)}{dx} \\
      \frac{dN_2(x)}{dx} 
    \end{array}
  \right)
  =
   \left(
    \begin{array}{c}
      \frac{-1}{l} \\
      \frac{1}{l} 
    \end{array}
  \right)
 (16)

The stress $ \ sigma $ is calculated from the stress-strain relation equations (2) and (15).

$\sigma = E \varepsilon = E{\bf B}^{{\rm T}} {\bf U} $  (17)

Can be expressed as.

Virtual displacement and strain

From equations (12) and (15), the virtual displacement $ \ delta u $ and the virtual strain $ \ delta \ varepsilon $ are calculated by matrix notation.

$\delta u = {\bf N}^{{\rm T}} \delta{\bf U} $  (18)         $\delta \varepsilon = {\bf B}^{{\rm T}} \delta{\bf U} $  (19)

Can be written as.

Element stiffness equation (element stiffness matrix)

Now, let's consider equation (4) of the principle of virtual work.

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 $  (4)

First, rewriting the inside of the integral on the left side of equation (4) with matrix notation,

\sigma \delta \varepsilon = E{\bf B}^{{\rm T}} {\bf U}{\bf B}^{{\rm T}} \delta{\bf U} = E {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U}  (20)

Therefore, if this is applied to the left side of equation (4),

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = \int_{x_1}^{x_2} E {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U}A dx $  (21)

Considering that the cross-sectional area $ A $ and Young's modulus $ E $ are constant in the elements in Eq. (21),

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = E {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U}A \int_{x_1}^{x_2} dx = E V {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U} $  (22)

Where $ V $ is $ A \ int_ {x_1} ^ {x_2} dx = El = V $, which represents the volume of the element. Also, if you set $ {\ bf K} = E V {\ bf B} {\ bf B} ^ {{\ rm T}} $, you can rewrite it as follows.

$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = {\bf U}^{{\rm T}} {\bf K} \delta{\bf U} $  (23)

$ {\ bf K} $ is called the element stiffness matrix.

Next, if the right side of equation (4) is expressed in a matrix,

f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 = {\bf F}^{{\rm T}} \delta{\bf U}  (24)

Where $ {\ bf F} $ is:

  {\bf F} = \left(
    \begin{array}{c}
      f_1 \\
      f_2 
    \end{array}
  \right)
 (25)

In this way, equation (4) of the principle of virtual work can be rewritten as follows using equations (23) and (24).

{\bf U}^{{\rm T}} {\bf K} \delta{\bf U} = {\bf F}^{{\rm T}} \delta{\bf U}  (26)

Also, if the above formula is summarized by $ \ delta {\ bf U} $,

({\bf U}^{{\rm T}} {\bf K} - {\bf F}^{{\rm T}} )\delta{\bf U} = {\bf 0}  (27)

Here, the principle of virtual work holds for any virtual displacement $ \ delta {\ bf U} $, so

{\bf U}^{{\rm T}} {\bf K} = {\bf F}^{{\rm T}}  (28)

If you take this transpose,

$ {\bf K}{\bf U} = {\bf F}$    (29)

It will be. Equation (29) is an equation called the element stiffness equation. As mentioned at the beginning, after formulating this element stiffness equation for each element, the stiffness equations of all the elements are integrated to construct the overall stiffness equation (overall stiffness matrix) and calculate the approximate solution. It becomes.

Calculation program of element stiffness equation (element stiffness matrix) of one-dimensional element

Now, I made a program to calculate the one-dimensional element stiffness equation with python. The library used is numpy and the code is below. From the given conditions, the stress can be calculated for one element by $ {\ bf B} $ matrix and element stiffness matrix $ {\ bf K} $.


import numpy as np

class Element:
    #Number of contacts
    NumberOfNodeInElement = 2
    
    def __init__(self):
        #Element length
        self.length = 0
        #Cross-sectional area of the element
        self.area = 0
        #Young's modulus
        self.young = 0
        #Contact numbers that make up the element
        self.ix = np.zeros(self.NumberOfNodeInElement)
        #Coordinate values of the contacts that make up the element
        self.xe = np.zeros(self.NumberOfNodeInElement)
        #Contact displacement of elements
        self.ue = np.zeros(self.NumberOfNodeInElement)
        #Bmatrix
        self.Bmatrix = np.zeros(self.NumberOfNodeInElement)
        #Element stiffness matrix
        self.Kmatrix = np.zeros((self.NumberOfNodeInElement,self.NumberOfNodeInElement))
    
    #Bmatrix calculation
    def makeBmatrix(self):
        self.length = self.xe[1]-self.xe[0]
        self.Bmatrix = np.array([[-1/self.length, 1/self.length]])

    #Calculation of element stiffness matrix
    def makeElementStiffness(self):
        self.volume = self.area*self.length
        self.Kmatrix = self.young*self.volume*np.dot(self.Bmatrix.T,self.Bmatrix)

    #Stress calculation
    def stressCalculation(self):
        return self.young*np.vdot(self.Bmatrix,self.ue)

For example, the calculation is performed on an element whose coordinates of node 1 and node 2 are $ x_1 = 50 $ (mm), $ x_2 = 150 $ (mm), cross-sectional area 100 $ mm ^ 2 $, and Young's modulus 200 GPa. When

Ele = Element()
Ele.xe = np.array([50,150])
Ele.area = 100
Ele.young = 200000
Ele.ue = np.array([0.01,0.025])

Ele.makeBmatrix()
Ele.makeElementStiffness()
scal = Ele.stressCalculation()

The calculation result in this case is as follows.

#Bmatrix
[[-0.01  0.01]]
#Bmatrix
[[ 200000. -200000.]
 [-200000.  200000.]]
#stress
30.000000000000004

Summary

This time, we explained up to the one-dimensional element stiffness matrix of the finite element method and introduced the element stiffness matrix fabrication program. Next time, we will assemble and calculate the overall stiffness equation (overall stiffness matrix).

Reference book

Shigeru Nagagi (2010) "First finite element method by Java" Corona Publishing Co., Ltd.

Recommended Posts

Finite element method beginning (one-dimensional element stiffness matrix)
Finite element method Beginning 2 (Preparation of one-dimensional overall stiffness matrix)
Solving one-dimensional wave equation by finite difference method (python)
Use IvyFEM (Finite Element Method Library for .NET) from Python