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.

Let \mathbf{L} be the strain-displacement differential operator for converting \mathbf{u} to \mathbf{\epsilon}.

Let \mathbf{E} be the 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. For convenience, we will apply the substitution E' = E/(1-\nu^2).

Newton’s second law at any point in the material may be written as follows:

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.

Let \mathbf{f}_b be the sum of other forces per unit volume, such as body forces like gravity and magnetism, or reference-frame forces like centrifugal force and coriolis force.

Basics of Finite Element

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.

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 1

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.

Example 2

Lets calculate \mathbf{B} for the trail functions from example 1 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.

Now lets rewrite Newton's second law from above in terms of the approximate solution.

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.

We will give each term its own variable and integrate over the domain of the element.

Where \mathbf{f}_{s} is the force due to the stress, \mathbf{K} is the stiffness, and M is the mass.

The term inside the integral for the stiffness matrix evaluate to the following matrix. \mathbf{B}^T_j\mathbf{E}\mathbf{B}_i =
\begin{bmatrix} E' \frac{\partial N_i}{\partial x} \frac{\partial N_j}{\partial x} + G \frac{\partial N_i}{\partial y} \frac{\partial N_j}{\partial y} & E'\nu \frac{\partial N_i}{\partial y} \frac{\partial N_j}{\partial x} + G \frac{\partial N_i}{\partial x} \frac{\partial N_j}{\partial y} \\ E'\nu \frac{\partial N_i}{\partial x} \frac{\partial N_j}{\partial y}+ G \frac{\partial N_i}{\partial y} \frac{\partial N_j}{\partial x} & E' \frac{\partial N_i}{\partial y} \frac{\partial N_j}{\partial y} + G \frac{\partial N_i}{\partial x} \frac{\partial N_j}{\partial x} \end{bmatrix}

Example 3

Now lets calculate mass values and stiffness matrices using the trails functions from 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_{21} = \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 matrices.

\mathbf K_11 = \frac{1}{2} \begin{bmatrix} E'+G & E'\nu+G \\ E'\nu+G & E'+G \end{bmatrix}
\mathbf K_21 = \mathbf K_12^T = \frac 1 2 \begin{bmatrix} -E' & -E'\nu \\ -G & -G \end{bmatrix}
\mathbf K_31 = \mathbf K_13^T = \frac 1 2 \begin{bmatrix} -G & -G \\ -E'\nu & -E' \end{bmatrix}
\mathbf K_32 = \mathbf K_23^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.

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.

Where \hat{\mathbf{N}} is the normal unit vector, \hat{\mathbf{T}} is the tangent unit vector, and ds is the differential length.

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 4

Lets calculate \mathbf f_s. We will consider a bearing stress \sigma_b on 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.

With these constraints, Newton's second law on the element simplifies as follows.

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)

First solve for (\mathbf{u}_u) in the top equation. Then solve for (\mathbf{f}_u) using the bottom equation.

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 5

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}_{12} \mathbf{u}_2 + \mathbf{K}_{13} \mathbf{u}_3 = \mathbf{f}_{s1}
\mathbf{K}_{21} \mathbf{u}_1 + \mathbf{K}_{22} \mathbf{u}_2 + \mathbf{K}_{23} \mathbf{u}_3 = \mathbf{f}_{s2}
\mathbf{K}_{31} \mathbf{u}_1 + \mathbf{K}_{32} \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}_{13} \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}_{23} \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.

Elements Working Together

Using multiple elements is not that much different from using a single element. The elements interact through their mutual nodes. Each element needs to have its properties calculated. As a review, the trail functions are defined for each node on the domain of the element, then the stiffness matrix is calculated for each node of the element. You may notice that the stiffness matrices for some of the nodes have already been calculated from previous elements; however, they do not apply to the domain of any new element you add. To get the full stiffness matrix you need to sum the matrices over the domains of all the elements.

Where, D_k is the domain of the k^{th} element and D is the union of all the element domains. The same applies for the mass.

Example 6

Expanding from the previous example 1 of a single triangular linear element, we will add a second element to create a square. The new element will share Nodes 2 and 3 and introduce Node 4 at (1, 1). The domain of the old element will be denoted D_1 and this new element will be denoted D_2.

D_2 = \{ (x, y) : x<1,\,y<1,\,1<x+y\}

N_2 = 1-y
N_3 = 1-x
N_4 = x+y-1

Below is our new system of equations.

\mathbf{K}_{11} \mathbf{u}_1 + \mathbf{K}_{12} \mathbf{u}_2 + \mathbf{K}_{13} \mathbf{u}_3 + \mathbf{K}_{14} \mathbf{u}_4 = \mathbf{f}_{s1}
\mathbf{K}_{21} \mathbf{u}_1 + \mathbf{K}_{22} \mathbf{u}_2 + \mathbf{K}_{23} \mathbf{u}_3 + \mathbf{K}_{24} \mathbf{u}_4 = \mathbf{f}_{s2}
\mathbf{K}_{31} \mathbf{u}_1 + \mathbf{K}_{32} \mathbf{u}_2 + \mathbf{K}_{33} \mathbf{u}_3 + \mathbf{K}_{34} \mathbf{u}_4 = \mathbf{f}_{s3}
\mathbf{K}_{41} \mathbf{u}_1 + \mathbf{K}_{42} \mathbf{u}_2 + \mathbf{K}_{43} \mathbf{u}_3 + \mathbf{K}_{44} \mathbf{u}_4 = \mathbf{f}_{s4}

The only unknown Nodes are \mathbf {u}_3 and \mathbf {u}_4, while \mathbf {u}_1 = \mathbf {u}_2 = \mathbf {0}.

\mathbf{K}_{33} \mathbf{u}_3 + \mathbf{K}_{34} \mathbf{u}_4 = \mathbf{f}_{s3}
\mathbf{K}_{43} \mathbf{u}_3 + \mathbf{K}_{44} \mathbf{u}_4 = \mathbf{f}_{s4}

Now we calculate the four stiffness matrices we need. Note, we already calculated the stiffness matrix \mathbf{K}_33; however, we need to add the new domain to it.

\mathbf{K}_33 = \frac 1 2 \begin{bmatrix} G & 0 \\ 0 & E' \end{bmatrix} + \frac 1 2 \begin{bmatrix} E' & 0 \\ 0 & G \end{bmatrix} = \frac 1 2 \begin{bmatrix} E' + G & 0 \\ 0 & E' + G \end{bmatrix}
\mathbf{K}_34 = \mathbf{K}^T_43 = \frac 1 2 \begin{bmatrix} -E' & -E'\nu \\ -G & -G \end{bmatrix}
\mathbf{K}_44 = \frac 1 2 \begin{bmatrix} E' + G & E'\nu + G \\ E'\nu + G & E' + G \end{bmatrix}

Now we solve for our unknown displacements.

\frac 1 2 \begin{bmatrix} E'+ G & 0 & -E' & -E'\nu \\ 0 & E'+G & -G & -G \\ -E' & -G & E'+G & E'\nu+G \\ -E'\nu & -G & E'\nu+G & E'+G \end{bmatrix} \begin{pmatrix} u_{x3} \\ u_{y3} \\ u_{x4} \\ u_{y4} \end{pmatrix} = \begin{pmatrix} f_{sx3} \\ f_{sy3} \\ f_{sx4} \\ f_{sy4} \end{pmatrix}
\begin{bmatrix} 20000 & 0 & -15000 & -4500 \\ 0 & 20000 & -5000 & -5000 \\ -15000 & -5000 & 20000 & 9500 \\ -4500 & -5000 & 9500 & 20000 \end{bmatrix} \begin{pmatrix} u_{x3} \\ u_{y3} \\ u_{x4} \\ u_{y4} \end{pmatrix} = \begin{pmatrix} 15 \\ 0 \\ 0 \\ 0 \end{pmatrix}
\begin{pmatrix} u_{x3} \\ u_{y3} \\ u_{x4} \\ u_{y4} \end{pmatrix} = \begin{pmatrix} 0.00195 \\ 0.00035 \\ 0.00168 \\ -0.00027 \end{pmatrix}