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.
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.
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).
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.
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.
Now suppose that inside each element, the displacement $ u $ changes with a linear function of coordinates $ x $.
This displacement must be equal to $ u_1 $, $ u_2 $ at nodes 1 and 2, so
Summarizing this for $ a $ and $ b $,
It will be. Where $ l = x_2-x_1 $ is the length of the element. Using (8) and (9), (5) is
Furthermore, by rewriting equation (10),
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.
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.
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.
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,
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,
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).
Also, if the above formula is summarized by $ \ delta {\ bf U} $,
Here, the principle of virtual work holds for any virtual displacement $ \ delta {\ bf U} $, so
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.
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
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).
Shigeru Nagagi (2010) "First finite element method by Java" Corona Publishing Co., Ltd.
Recommended Posts