Finite Element Method (Simplex, Linear Elastic)

Finite Element method is a numerical method, rather than analytic method of solving PDEs.

Interpolation function

  • Defines the type of polynomials used approximate the real solution
  • Higher order polynomials have some flexibilities to approximate the actual solution more accurately
  • However it requires more information and operations to compute a higher order element. For example, in 1D, a 2nd order function requires information from 3 nodes to determine the 3 constants:
    • \phi = \alpha_1 + \alpha_2 x  2 nodes
    • \phi = \alpha_1 + \alpha_2 x + \alpha_3 x^2  3 nodes
FEM_1DElementOrder

Upper: 2-node element; Lower: 3-node element

  • In 2D, a simple linear equation would have requried 3 nodes to inform the 3 parameters:
    • \phi = \alpha_1 + \alpha_2 x + \alpha_3 y  3 nodes
    • \phi = \alpha_1 + \alpha_2 x + \alpha_3 y + \alpha_4 xy + \alpha_5 x^2 + \alpha_6 y^2  6 nodes
FEM_2DElementOrder

Upper: 3-node 2D element for linear interpolation function; Lower: 6-node 2D element for quadratic interpolation function

Example problem: Plane Stress/Strain and Isotropic Elasticity

We try to minimise the potential energy \Pi of the system:

\Pi = \Lambda - W

where \Lambda is the strain energy of the structure and W is the work done by external force.

Strain Energy

To calculate the strain energy of the structure \Lambda, we have \Lambda = \int_{V} \frac{1}{2} \varepsilon^T \sigma dV.

  • For 1D linear elasticity, we have \sigma = E\varepsilon.
  • For 2D linear elasticity, we have \mathbf{\sigma} = \mathbf{D\varepsilon}.
    • Inside this equation:
    • \mathbf{\sigma}^T = \begin{Bmatrix} \sigma_x & \sigma_y & \tau_{xy} \end{Bmatrix}
    • For Plane stress: \mathbf{D} = \frac{E}{(1-v^2)} \begin{bmatrix}  1 & v & 0 \\ v & 1 & 0 \\ 0 & 0 & \frac{1-v}{2} \end{bmatrix}  
    • For Plane strain: \mathbf{D} = \frac{E}{(1+v)(1-2v)}  \begin{bmatrix}  1-v & v & 0 \\ v & 1-v & 0 \\ 0 & 0 & \frac{1-2v}{2} \end{bmatrix}  
    • \mathbf{\varepsilon}^T = \begin{Bmatrix} \varepsilon_x & \varepsilon_y & \gamma_{xy}\end{Bmatrix}

At each node: \mathbf{u = Nu_e}

  • \mathbf{u_e} = \begin{Bmatrix} u_i & v_i & u_j & v_j & ... & v_n \end{Bmatrix}
  • Remember the following:
    • \varepsilon_x = \frac{\partial u}{\partial x}
    • \varepsilon_y = \frac{\partial v}{\partial y}
    • \gamma_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}
  • After differentiation, \mathbf{\varepsilon = Bu_e}
    • ie. for a 2-D triangular element: u = N_iu_i+N_ju_j+N_ku_k
    • \varepsilon_x = \frac{\partial u}{\partial x} = \frac{\partial N_i}{\partial x}u_i+\frac{\partial N_j}{\partial x}u_j+\frac{\partial N_k}{\partial x}u_k

For simplex elements, \frac{\partial \mathbf{N}}{\partial x} return constants.

Then it sums up to the element stiffness matrix \mathbf{K_e}:

\Lambda = \frac{1}{2} \int_{V_e} \mathbf{u^T_e} \mathbf{B^T} \mathbf{DBu_e} dV_e = \frac{1}{2}\mathbf{u}^T_e \mathbf{u}_e \mathbf{K_e}

where \mathbf{K_e} = \int_{V_e} \mathbf{B^T DB} dV_e.

Work done by external force

For simplicity, we assume a nodal load: W = \mathbf{u^T_ef_e}

Assembling the equations of the entire system: \Pi = \sum^{E}_{e=1} \frac{1}{2} \mathbf{u^T_eK_eu_e} - \mathbf{u^T_ef_e}

To get:

\frac{\partial \Pi}{\partial U} = \sum^{E}_{e=1} \mathbf{K_eU}-\mathbf{F} = 0

\mathbf{K} = \sum^{E}_{e=1} \mathbf{K_e} is called the global stiffness matrix.

2D Simplex elements

FEM_2DSimplex

2D Simplex element definition

Nodal information: \begin{Bmatrix}u \\ v \end{Bmatrix} = \begin{bmatrix} N_i & 0 & N_j & 0 & N_k & 0 \\ 0 & N_i & 0 & N_j & 0 & N_k \end{bmatrix} \begin{Bmatrix} u_{2i-1} \\ u_{2i} \\u_{2j-1} \\ u_{2j} \\u_{2k-1} \\ u_{2k} \end{Bmatrix} = \mathbf{Nu_e}

Shape function: N_\beta = \frac{1}{2A}(a_\beta+b_\beta x+c_\beta y) for \beta \in \{i,j,k\}

Nodal Strain:

  • \varepsilon_x = \begin{bmatrix} \frac{\partial N_i}{\partial x} & 0 & \frac{\partial N_j}{\partial x} & 0 & \frac{\partial N_k}{\partial x} & 0 \end{bmatrix} \mathbf{u_e} = \frac{1}{2A} \begin{bmatrix} b_i & 0 & b_j & 0 & b_k & 0 \end{bmatrix} \mathbf{u_e}
  • \varepsilon_y = \begin{bmatrix} 0 & \frac{\partial N_i}{\partial y} & 0 & \frac{\partial N_j}{\partial y} & 0 & \frac{\partial N_k}{\partial y} \end{bmatrix} \mathbf{u_e} = \frac{1}{2A} \begin{bmatrix} 0 & c_i & 0 & c_j & 0 & c_k \end{bmatrix} \mathbf{u_e}
  • \gamma_{xy} = \begin{bmatrix} \frac{\partial N_i}{\partial y} & \frac{\partial N_i}{\partial x} & \frac{\partial N_j}{\partial y} & \frac{\partial N_i}{\partial x} & \frac{\partial N_k}{\partial y} & \frac{\partial N_i}{\partial x} \end{bmatrix} \mathbf{u_e} = \frac{1}{2A} \begin{bmatrix} c_i & b_i & c_j & b_j & c_k & b_k \end{bmatrix} \mathbf{u_e}
  • Boom!
  • Assemble these equations to give \mathbf{B}:
  • \begin{Bmatrix} \varepsilon_x \\ \varepsilon_y \\ \gamma_{xy} \end{Bmatrix} = \frac{1}{2A} \begin{bmatrix} b_i & 0 & b_j & 0 & b_k & 0 \\ 0 & c_i & 0 & c_j & 0 & c_k \\ c_i & b_i & c_j & b_j & c_k & b_k \end{bmatrix} \mathbf{u_e} = \mathbf{Bu_e}

Get \mathbf{K} :

\mathbf{K_e} = \int_V \mathbf{B^T DB} dV = \mathbf{B^T DB} \int_V dV = \mathbf{B^T DB} tA

where t is the thickness, and A is the element area.

Now we get \mathbf{K_e}, how do we construct $latex \mathbf{K} = \sum^{E}_{e=1} \mathbf{K_e}$?

We do a numerical integration using Gaussian quadrature:

\int^1_{-1} f(x) dx = \sum^n_{i=1} w_i f(x_i)

which means that we numerically integrate function f(x) by summing up the function output f(x_i) at sample points x_i multiplied by the corresponding weights w_i. n is the number of sample points.

The weight w_i is calculated by w_i = \frac{2}{(1-x_i^2)(P'_n (x_i))^2}  with P_n(x) being the Legendre polynomial. Or (a look-up table):

No. of points n Points x_i Weights w_i
1 0 2
2 \pm \sqrt{\frac{1}{3}} 1
3 0 \frac{8}{9}
\pm \sqrt{\frac{3}{5}} \frac{5}{9}
4 \pm \sqrt{\frac{3}{7}-\frac{2}{7}\sqrt{\frac{6}{5}}} \frac{18+\sqrt{30}}{36}
\pm \sqrt{\frac{3}{7}-\frac{2}{7}\sqrt{\frac{6}{5}}} \frac{18-\sqrt{30}}{36}
5 0 \frac{128}{225}
\pm \frac{1}{3} \sqrt{5-2\sqrt{\frac{10}{7}}} \frac{322+13\sqrt{70}}{900}
\pm \frac{1}{3} \sqrt{5+2\sqrt{\frac{10}{7}}} \frac{322-13\sqrt{70}}{900}

 

This sets up the global stiffness matrix. Together with boundary conditions, initial conditions and the finite element mesh, we can solve the linear/non-linear matrix to obtain nodal results. This current example is a linear model. We will soon explore a non-linear mechanical model.

Sources:

One thought on “Finite Element Method (Simplex, Linear Elastic)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s