Learn Finite Element Analysis
Intro to Mechanics of Materials
Starting with a 2D solid material. Let \mathbf{u} be the
displacement field of the material. Let \mathbf{\epsilon} be the normal
(\epsilon) and shear (\gamma) strains. Let
\mathbf{\sigma} be the normal (\sigma) and
shear (\tau) stresses.
\mathbf{u} =
\begin{pmatrix}
u_x \\ u_y
\end{pmatrix} \:\:\:\: displacement
\mathbf{\epsilon} =
\begin{pmatrix}
\epsilon_x \\ \epsilon_y \\ \gamma_{xy}
\end{pmatrix} =
\begin{pmatrix}
\frac{\partial u_x}{\partial x} \\
\frac{\partial u_y}{\partial y} \\
\frac{\partial u_x}{\partial y} +
\frac{\partial u_y}{\partial x}
\end{pmatrix} \:\:\:\: strain
\mathbf{\sigma} =
\begin{pmatrix}
\sigma_x \\ \sigma_y \\ \tau_{xy}
\end{pmatrix} \:\:\:\: stress
Let \mathbf{L} be a differential operator for converting
\mathbf{u} to \mathbf{\epsilon}.
\mathbf{L(u)}=
\begin{bmatrix}
\frac{\partial}{\partial x} & 0 \\
0 & \frac{\partial}{\partial y} \\
\frac{\partial}{\partial y} &
\frac{\partial}{\partial x}
\end{bmatrix}
\begin{pmatrix}
u_x \\ u_y
\end{pmatrix}
= \mathbf{\epsilon}
Let \mathbf{E} be an elasticity matrix to convert strain into
stress, where E is Young’s modulus, G is the
modulus of rigidity, and \nu is Poisson’s ratio. The matrix
shown is for the plane stress condition; that is \sigma_z = 0.
\mathbf{\sigma} =
\mathbf{E\epsilon} =
\begin{bmatrix}
\frac{E}{1-\nu^2} & \frac{E\nu}{1-\nu^2} & 0 \\
\frac{E\nu}{1-\nu^2} & \frac{E}{1-\nu^2} & 0 \\
0 & 0 & G
\end{bmatrix}
\begin{pmatrix}
\epsilon_x \\ \epsilon_y \\ \gamma_{xy}
\end{pmatrix}
Newton’s second law at any point in the material may be written as follows:
\mathbf{f} = \rho
\frac{\partial^2 \mathbf{u}}{\partial t^2}
Where \mathbf{f} is the sum of forces per unit volume and
\rho is the density. We are concerned with the
stress in the material, but other forces can be added, as needed. Let
\mathbf{f}_s be the net force due to stress.
\mathbf{f}_s=
\mathbf{L}^T(\mathbf{\sigma}) =
\begin{bmatrix}
\frac{\partial}{\partial x} & 0 & \frac{\partial}{\partial y} \\
0 & \frac{\partial}{\partial y} & \frac{\partial}{\partial x}
\end{bmatrix}
\begin{pmatrix}
\sigma_x \\
\sigma_y \\
\tau_{xy}
\end{pmatrix}
Let \mathbf{f}_b be the sum of other forces per unit volume,
such as body forces like gravity and magnetism, reference-frame forces like
centrifugal force and coriolis force.
\mathbf{f}_b + \mathbf{L}^T(\mathbf{\sigma}) =
\rho \frac{\partial^2 \mathbf{u}}{\partial t^2}
Basics of Finite Elements
Finite element method is used to solve partial differential equations on domains
with complex geometry. This is in contrast to finite difference method which usually
requires a domain to be easily divided into equal sized tiles or blocks. Finite
element method uses elements made up of nodes. Adjacent elements share the nodes
at their common boundary. The solution is approximated as the sum below.
\mathbf{u}(\mathbf{x},t)\approx\sum_i{N_i(\mathbf{x}) \mathbf{u}_i(t)}
Where N_i are trail functions and \mathbf u_i
are the values of the solution at each node. The trail functions are determined
by the domain and the type of elements used. A trial function for node i will have
a value of 1 at node i and 0 at all other nodes. In addition, the trail function
will be 0 outside of any elements it is a part of. Depending on the type of
element, the trail function may be linear, quadratic, or some other distribution.
When all the trail functions in an element are summed, the value across the entire
element is 1.
Example
Lets look at an example of
trail functions. The simplest 2D element is a linear triangular element. This
element is not accurate for structural finite element, so we will look at better
elements later. Let
D be the triangular domain of an element
between nodes at (0, 0), (1, 0), and (0, 1).
D = \{ (x, y) : 0<x,\,0<y,\,x+y<1\}
Each node has its own trail function. The trail function is 1 at its node
and zero at the other nodes.
N_1 = 1-x-y
N_2 = x
N_3 = y
Since the approximate solution is a product of a space-only dependent
function (N_i(\mathbf{x})) and a time-only dependent value
(\mathbf{u}_i(t)), the stress and strain only depend on the
derivatives of the trail functions and only the value at the node. For this
reason, we let \mathbf{B} be the derivative of the trail
function used to calculate stress and strain.
\mathbf{B} = \mathbf{L}(N) =
\begin{bmatrix}
\frac{\partial N}{\partial x} & 0 \\
0 & \frac{\partial N}{\partial y} \\
\frac{\partial N}{\partial y} &
\frac{\partial N}{\partial x}
\end{bmatrix}
Example
Lets calculate \mathbf{B} for the trail functions from the
example above.
\mathbf{B}_1 =
\begin{bmatrix}
-1 & 0 \\
0 & -1 \\
-1 & -1
\end{bmatrix} ,\:
\mathbf{B}_2 =
\begin{bmatrix}
1 & 0 \\
0 & 0 \\
0 & 1
\end{bmatrix},\:
\mathbf{B}_3 =
\begin{bmatrix}
0 & 0 \\
0 & 1 \\
1 & 0
\end{bmatrix}
The strain and stress can be approximated using the equations below.
\mathbf{\epsilon}\approx\sum_i{\mathbf{B}_i \mathbf{u}_i(t)}
,\:
\mathbf{\sigma}\approx\mathbf{E}\sum_i{\mathbf{B}_i \mathbf{u}_i(t)}
Now lets rewrite Newton's second law from above in terms of the approximate
solution.
\mathbf{f}_b + \mathbf{L}^T(\mathbf{\sigma}) =
\sum_i{\rho N_i \frac{\mathrm{d} \mathbf{u}_i}{\mathrm{d} t^2}}
We will multiply Newton's second law above by N_j, the nodes'
trail functions, to generate a system of ordinary differential equations that
can be solved.
N_j \mathbf{f}_b + N_j \mathbf{L}^T(\mathbf{\sigma}) =
\sum_i{\rho N_i N_j \frac{\mathrm{d} \mathbf{u}_i}{\mathrm{d} t^2}}
We apply the product rule to the stress term. This will allow us to break up
the term into an internal stress force and an external stress force applied at
the boundary.
N_j \mathbf{f}_b +
\mathbf{L}^T(N_j \mathbf{\sigma}) -
\mathbf{L}^T(N_j) \mathbf{\sigma} =
\sum_i{\rho N_i N_j \frac{\mathrm{d} \mathbf{u}_i}{\mathrm{d} t^2}}
We will give each term its own variable and integrate over the domain of
the element.
\mathbf{f}_{bj} + \mathbf{f}_{sj} - \sum_i \mathbf{K}_{ij}\mathbf{u}_i =
\sum_i M_{ij} \frac{\mathrm{d} \mathbf{u}_i}{\mathrm{d} t^2}
\mathbf{f}_{bj} = \iint_D N_j \mathbf{f}_b \, \mathrm{dA}
\mathbf{f}_{sj} = \iint_D \mathbf{L}^T(N_j \mathbf{\sigma}) \, \mathrm{dA}
\mathbf{K}_{ij}\mathbf{u}_i =
\iint_D \mathbf{B}^T_j \mathbf{\sigma} \, \mathrm{dA} =
\left( \iint_D \mathbf{B}^T_j\mathbf{E}\mathbf{B}_i \mathrm{dA} \right) \mathbf{u}_i
M_{ij} = \iint_D \rho N_i N_j \mathrm{dA}
Where \mathbf{f}_{s} is the force due to the stress,
\mathbf{K} is the stiffness, and
M is the mass.
Example
Now lets calculate mass values and stiffness matrices using the trails
functions from the example above. We will start with the mass value for i=1 and
j=1. For convenience, we will apply
the substitution v = 1-x and w = v-y.
M_{11} = \rho \iint_D {N_1}^2 \mathrm{dA}
= \rho \int_0^1 \int_0^{1-x} (1-x-y)^2 \text{ dy dx}
= \rho \int_0^1 \int_0^{v} w^2 \text{ dw dv}
= \frac{\rho}{3} \int_0^1 v^3 \text{ dv}
= \frac{\rho}{12}
The result is the same for M_{22} and
M_{33}. Now we look at i=1 and
j=2.
M_{12} = \rho \iint_D N_1 N_2 \text{ dA}
= \rho \int_0^1 \int_0^{1-x} (1-x-y)x \text{ dy dx}
= \rho \int_0^1 (1-v) \int_0^{v} w \text{ dw dv}
= \frac{\rho}{2} \int_0^1 (1-v)v^2 \text{ dv}
= \frac{\rho}{24}
The result is the same for M_{13}, M_{21},
M_{23}, M_{31}, and M_{32}.
Written as a matrix, the result looks like this.
\mathbf{M} = \frac{\rho}{24}
\begin{bmatrix}
2 & 1 & 1 \\
1 & 2 & 1 \\
1 & 1 & 2
\end{bmatrix}
Now we will look at the stiffness matrix for i=1 and
j=1. For convenience, we will apply
the substitution E' = E/(1-\nu^2).
\mathbf{K}_{11} =
\iint_D \mathbf{B}^T_1\mathbf{E}\mathbf{B}_1 \mathrm{dA}
= \mathbf{B}^T_1
\begin{bmatrix}
E' & E'\nu & 0 \\
E'\nu & E' & 0 \\
0 & 0 & G
\end{bmatrix}
\begin{bmatrix}
-1 & 0 \\
0 & -1 \\
-1 & -1
\end{bmatrix}
\iint_D \mathrm{dA}
= \frac{1}{2}
\begin{bmatrix}
-1 & 0 & -1 \\
0 & -1 & -1
\end{bmatrix}
\begin{bmatrix}
-E' & -E'\nu \\
-E'\nu & -E' \\
-G & -G
\end{bmatrix}
\mathbf K_11 = \frac{1}{2}
\begin{bmatrix}
E'+G & E'\nu+G \\
E'\nu+G & E'+G
\end{bmatrix}
And the other stiffness matrices...
\mathbf K_12 = \mathbf K_21^T =
\frac 1 2
\begin{bmatrix}
-E' & -E'\nu \\
-G & -G
\end{bmatrix}
\mathbf K_13 = \mathbf K_31^T =
\frac 1 2
\begin{bmatrix}
-G & -G \\
-E'\nu & -E'
\end{bmatrix}
\mathbf K_23 = \mathbf K_32^T =
\frac 1 2
\begin{bmatrix}
0 & G \\
E'\nu & 0
\end{bmatrix}
\mathbf K_22 =
\frac 1 2
\begin{bmatrix}
E' & 0 \\
0 & G
\end{bmatrix}, \:
\mathbf K_33 =
\frac 1 2
\begin{bmatrix}
G & 0 \\
0 & E'
\end{bmatrix}
It is relatively straight forward to calculate these values, except
\mathbf{f}_{sj}. Here we need to use Green's theorem. This shows
us that \mathbf{f}_{sj} is only a result of the stress at the
boundary. \partial D is the boundary of domain D.
\mathbf{f}_{sj} = \iint_D \mathbf{L}^T(N_j \mathbf{\sigma}) \, \mathrm{dA}
=\iint_D
\begin{pmatrix}
\frac{\partial}{\partial x}(N_j\sigma_x) +
\frac{\partial}{\partial y}(N_j\tau_{xy})\\
\frac{\partial}{\partial x}(N_j\tau_{xy}) +
\frac{\partial}{\partial y}(N_j\sigma_y)
\end{pmatrix}
\mathrm{dA}
=\oint_{\partial D}
N_j
\begin{pmatrix}
\sigma_x dy -
\tau_{xy} dx \\
\tau_{xy} dy -
\sigma_y dx
\end{pmatrix}
It is more useful to apply boundary conditions in terms of
the bearing stress \sigma_b and the shear friction
stress \tau_f on the surface.
\mathbf{f}_{sj} =
\oint_{\partial D}
N_j
\left(
\sigma_b \, \hat{\mathbf{N}} +
\tau_f \, \hat{\mathbf{T}}
\right)
ds
Where \hat{\mathbf{N}} is the normal unit vector,
\hat{\mathbf{T}} is the tangent unit vector, and
ds is the differential length.
\hat{\mathbf{N}} =
\begin{pmatrix}
\frac{dy}{ds} \\
-\frac{dx}{ds}
\end{pmatrix},\:
\hat{\mathbf{T}} =
\begin{pmatrix}
\frac{dx}{ds} \\
\frac{dy}{ds}
\end{pmatrix}
ds = \sqrt{dx^2 + dy^2}
Note \sigma_b is negative for a compressive stress.
\tau_f, \hat\mathbf T and \oint_{\partial D} go
counterclockwise around D. \hat\mathbf N
points outward.
Example
Lets calculate \mathbf f_s. We will consider a bearing stress
\sigma_bon the left side of the element.
let,
x = 0
y = 1-t
N_1 = 1-x-y = t , \: \tau_f = 0
\hat{\mathbf{N}} \, ds=
\begin{pmatrix}
\frac{dy}{dt} \\
-\frac{dx}{dt}
\end{pmatrix} \, dt=
\begin{pmatrix}
-1 \\ 0
\end{pmatrix} \, dt
\mathbf{f}_{s1} =
\int_{\partial D}
N_1
\left(
\sigma_b \, \hat{\mathbf{N}} +
\tau_f \, \hat{\mathbf{T}}
\right)
ds
= \sigma_b
\begin{pmatrix}
-1 \\ 0
\end{pmatrix}
\int_0^1 t \, dt
= \frac{\sigma_b}{2}
\begin{pmatrix}
-1 \\ 0
\end{pmatrix}
N_2 is zero on the entire left side.
N_2 = x = 0 , \:
\mathbf{f}_{s2} =
\begin{pmatrix}
0 \\ 0
\end{pmatrix}
N_3 = y = 1 - t
\mathbf{f}_{s3} =
\sigma_b
\begin{pmatrix}
-1 \\ 0
\end{pmatrix}
\int^{1}_{0} 1 - t \, dt
= \frac{\sigma_b}{2}
\begin{pmatrix}
-1 \\ 0
\end{pmatrix}
Finally, lets solve for the displacement, strain, and stress of element.
For now, we will not consider any body forces, just the external forces.
Additionally, we will only consider static equilibrium.
These two constraints can be summarized by the two equations below.
\mathbf{f}_{bj} =
\begin{pmatrix}
0 \\ 0
\end{pmatrix} ,\:
\frac{\mathrm{d} \mathbf{u}_i}{\mathrm{d} t^2} =
\begin{pmatrix}
0 \\ 0
\end{pmatrix}
With these constraints, Newton's second law on the element simplifies as
follows.
\sum_i \mathbf{K}_{ij}\mathbf{u}_i =
\mathbf{f}_{sj}
This is a set of linear equations, two equations for each node. Unfortunately
this system of equations cannot be solved as-is. Right now the system has infinite
solutions. If you imagine a force applied to the element, there is nothing
keeping it from flying off. Some of the nodes need to be restrained. The
combination of restrained nodes and applied loads provides our boundary conditions.
Nodes are restrained by setting their displacement. Now we will rearrange the
equations by the equations with
known forces (\mathbf{f}_k), and the equations with unknown forces
(\mathbf{f}_u). The terms of the equations are rearranged in terms of
known displacements (\mathbf{u}_k) and unknown displacements
(\mathbf{u}_u)
\mathbf{K}_{kk}\mathbf{u}_k +
\mathbf{K}_{ku}\mathbf{u}_u=
\mathbf{f}_k
\mathbf{K}_{uk}\mathbf{u}_k +
\mathbf{K}_{uu}\mathbf{u}_u=
\mathbf{f}_u
First solve for (\mathbf{u}_u) in the top equation.
Then solve for (\mathbf{f}_u) using the bottom equation.
\mathbf{K}_{ku}\mathbf{u}_u=
\mathbf{f}_k -
\mathbf{K}_{kk}\mathbf{u}_k
\mathbf{f}_u =
\mathbf{K}_{uk}\mathbf{u}_k +
\mathbf{K}_{uu}\mathbf{u}_u
The top equation can be solved directly using
Gaussian elimination or approximated using Jacobi or Gauss-Seidel iterative
methods. Gaussian elimination is slow, while Jacobi or Gauss-Seidel iterative
methods are faster and easier to implement but less accurate. Gauss-Seidel
method is faster and uses less memory than Jacobi method; however, Jacobi
method lends itself well to massively parallel computing that can be run on GPUs.
There is a risk that Gauss-Seidel or Jacobi method will not converge on a
solution; however, static structural FEA has a good chance of converging with
these methods.
Example
Using the element from the previous examples we will calculate the the displacement,
strain, and stress. The element has three nodes, for a total of six equations to be solved.
\mathbf{K}_{11} \mathbf{u}_1 +
\mathbf{K}_{21} \mathbf{u}_2 +
\mathbf{K}_{31} \mathbf{u}_3 =
\mathbf{f}_{s1}
\mathbf{K}_{12} \mathbf{u}_1 +
\mathbf{K}_{22} \mathbf{u}_2 +
\mathbf{K}_{32} \mathbf{u}_3 =
\mathbf{f}_{s2}
\mathbf{K}_{13} \mathbf{u}_1 +
\mathbf{K}_{23} \mathbf{u}_2 +
\mathbf{K}_{33} \mathbf{u}_3 =
\mathbf{f}_{s3}
We will use the following values to get a concrete solution. Nodes 1 and 2
will be restrained to have no displacement.
E' = 30000 , \:\:
\nu = 0.3 , \:\:
G = 10000 , \:\:
\sigma_b = -30
\mathbf{u}_1 =
\mathbf{u}_2 =
\begin{pmatrix}
0 \\
0
\end{pmatrix}
To solve for \mathbf{u}_3, we use the third equation from
above. The equation has been simplified based on the constraints on Nodes 1 and 2.
\mathbf{K}_{33} \mathbf{u}_3 =
\mathbf{f}_{s3}
\frac 1 2
\begin{bmatrix}
G & 0_{} \\
0_{} & E'
\end{bmatrix}
\begin{pmatrix}
u_{x3} \\
u_{y3}
\end{pmatrix} =
\begin{pmatrix}
f_{sx3} \\
f_{sy3}
\end{pmatrix}
\begin{bmatrix}
5000 & 0_{} \\
0_{} & 15000
\end{bmatrix}
\begin{pmatrix}
u_{x3} \\
u_{y3}
\end{pmatrix} =
\begin{pmatrix}
15_{} \\
0_{}
\end{pmatrix}
u_{x3} = 0.003, \:
u_{y3} = 0
Now we can calculate the forces at Nodes 1 and 2.
\mathbf{f}_{s1} =
\mathbf{K}_{31} \mathbf{u}_3
\begin{pmatrix}
f_{sx1} \\
f_{sy1}
\end{pmatrix} =
\frac 1 2
\begin{bmatrix}
-G & -E'\nu \\
-G & -E'
\end{bmatrix}
\begin{pmatrix}
u_{x3} \\
u_{y3}
\end{pmatrix}
\begin{pmatrix}
f_{sx1} \\
f_{sy1}
\end{pmatrix} =
\begin{bmatrix}
-5000 & -4500 \\
-5000 & -15000
\end{bmatrix}
\begin{pmatrix}
0.003 \\
0
\end{pmatrix}
f_{sx1} = -15, \:
f_{sy1} = -15
\mathbf{f}_{s2} =
\mathbf{K}_{32} \mathbf{u}_3
f_{sx2} = 0, \:
f_{sy2} = 15
Now strains and stresses.
\mathbf{\epsilon} = \sum_i{\mathbf{B}_i \mathbf{u}_i} =
\mathbf{B}_3 \mathbf{u}_3
\begin{pmatrix}
\epsilon_x \\ \epsilon_y \\ \gamma_{xy}
\end{pmatrix} =
\begin{bmatrix}
0 & 0 \\
0 & 1 \\
1 & 0
\end{bmatrix}
\begin{pmatrix}
0.003 \\
0
\end{pmatrix} =
\begin{pmatrix}
0 \\ 0 \\ 0.003
\end{pmatrix}
\mathbf{\sigma} =
\mathbf{E\epsilon}
\begin{pmatrix}
\sigma_x \\ \sigma_y \\ \tau_{xy}
\end{pmatrix} =
\begin{bmatrix}
E' & E'\nu & 0 \\
E'\nu & E' & 0 \\
0 & 0 & G
\end{bmatrix}
\begin{pmatrix}
\epsilon_x \\ \epsilon_y \\ \gamma_{xy}
\end{pmatrix}
=
\begin{bmatrix}
30000 & 9000 & 0 \\
9000 & 30000 & 0 \\
0 & 0 & 10000
\end{bmatrix}
\begin{pmatrix}
0 \\ 0 \\ 0.003
\end{pmatrix}
\begin{pmatrix}
\sigma_x \\ \sigma_y \\ \tau_{xy}
\end{pmatrix} =
\begin{pmatrix}
0 \\ 0 \\ 30
\end{pmatrix}
Notice these stresses and strains are uniform throughout the entire element.
This makes for a poor approximation of the actual stresses and strains. This is
why linear elements are bad for modeling material mechanics. An actual part
with these loads would have non-uniform stresses, including bending. A single
linear element cannot not model bending or transverse shear stress. To
properly model these stresses we would need many linear elements or a few more
sophisticated elements, like quadratic elements.