Chapter 4: Finite Element Method

4.1 — Introduction to FEM
Collected & Prepared By: Er. Sanjay Kumar Sah (Lecturer, Department of Civil Engineering, Pulchowk Campus)

Element Shapes

(i) One Dimensional Element

(ii) Two Dimensional Elements

a) Constant Strain Triangle (CST)

→ Linear displacement triangle — 3-noded triangle.

b) Linear Strain Triangle (LST)

→ Quadratic displacement triangle — 6-noded triangle.

c) Quadratic Strain Triangle

→ Cubic displacement triangle — 10-noded triangle.

3-Noded (CST) 3 1 2 6-Noded (LST) 3 1 2 6 5 4 10-Noded (Cubic) 3 1 2 8 7 9 6 4 5 10 Corner node Mid-edge Internal
Fig: Visual representation of triangular elements — CST, LST, and Cubic

⇒ One can also use higher-order triangular elements like cubic strain triangle (15-noded triangle) and quartic strain triangle (21-noded triangle).

(d) Lagrange Elements

→ Elements in which nodes are in a grid pattern.

(e) Serendipity Elements

→ Elements in which nodes are provided in boundaries only.

(f) Curved 2-D Elements

→ Development of iso-parametric concepts has facilitated using curved elements also.

(iii) Three Dimensional Elements

a) Tetrahedron (4-noded solid)

4 1 3 2
Fig: Tetrahedron element (4-noded solid)

b) Brick Element (8-noded solid)

A standard 3D rectangular box/hexahedron with 8 corner nodes.

(iv) Axi-symmetric Elements

→ Generally 2D/3D elements having nodes all around about an axis.

Nodes

Definition

Nodes are specific finite points that are to be determined in finite element analysis. The basic unknowns at any points other than nodes are to be determined using interpolation functions in terms of nodal parameter values of elements.

Types of Nodes

(Example: In the 10-noded triangle shown above, nodes 1 through 9 are external boundary nodes, and node 10 is an internal node.)

Steps Involved in FEM

  1. Discretization of structure into the elements.
  2. Selection of displacement models.
  3. Derivation of stiffness matrix for individual element.
  4. Assembling stiffness matrix.
  5. Solution of unknown displacement.
  6. Computation of element stress & strain.

a) Discretization of Structure into the Elements

The given structure is divided into an equivalent system of finite elements. The equivalent system may consist of bar, truss, triangular area, quadrilateral area, tetrahedron solid, hexahedron solids etc. in 2D or 3D coordinate systems.

b) Selection of Displacement Models

It is generally not possible to select the displacement function that can represent exactly the actual variation of displacement in the element. Hence, basic approximation of the FEM is introduced at this stage. There are three inter-related factors which influence the selection of the displacement models:

  1. The type and degree of displacement model must be chosen.
  2. The particular displacement magnitude that describes the model must be selected. The displacements at nodal points are known after analysis, but they may also include the derivation of the displacement at some other points inside the element.
  3. Displacement must satisfy certain requirements, which ensures that the numerical result approaches the correct solutions.

c) Derivation of the Element Stiffness Matrix

The stiffness matrix consists of the coefficients of the equilibrium equation derived from the properties of material and element and obtained by any method. The stiffness relates the displacement at the nodal points to the applied force at the corresponding nodal points. The distributed forces applied to the structure are converted into equivalent concentrated forces at nodes. The equilibrium equation relates the stiffness matrix $[K]$, nodal force vector $\{f\}$, and the nodal displacement vector $\{u\}$.

1 2 F₁, V₁ M₁, θ₁ F₂, V₂ M₂, θ₂ L
Fig: Force and displacement acting at corresponding nodes (Beam Element)

Displacement matrix:

$$\{u\} = \begin{Bmatrix} V_1 \\ \theta_1 \\ V_2 \\ \theta_2 \end{Bmatrix}$$

Force matrix:

$$\{f\} = \begin{Bmatrix} F_1 \\ M_1 \\ F_2 \\ M_2 \end{Bmatrix}$$

Element Stiffness Matrix (for a beam element):

$$[K] = \frac{EI}{L^3} \begin{bmatrix} 12 & 6L & -12 & 6L \\ 6L & 4L^2 & -6L & 2L^2 \\ -12 & -6L & 12 & -6L \\ 6L & 2L^2 & -6L & 4L^2 \end{bmatrix}$$

d) Assembling Stiffness Matrix

The process includes the assemblage of the overall or global stiffness matrix for the entire body/structure from the individual element stiffness matrix. A typical assembled matrix for a structure having 2 beam elements overlaps at the common node.

$$\{F\} = [K]\{u\}$$

Visual concept of overlapping stiffness matrices:

$$[K_{\text{global}}] = \begin{bmatrix} [K_1] & & \\ & [K_1] + [K_2] & \\ & & [K_2] \end{bmatrix}$$

e) Solution of Unknown Displacement

The algebraic equations obtained in the above steps are solved for the unknown displacements. In a linear equilibrium equation, there is a straightforward application of matrix algebra techniques. However, for non-linear problems, the desired solutions are obtained by a sequence of steps, each step involving the modification of the stiffness matrix.

$$\text{Displacement at nodal points:} \quad \{u\} = [K]^{-1} \{F\}$$

f) Computation of the Element Stress and Strain

In some cases, the magnitude of the primary unknowns, that is the nodal displacement, will be all that is required for the engineering solutions. However, the other quantities derived from the primary unknowns like stress, strain must be computed.

Force acting on each element is:

$$F_i = K_i \cdot u_i$$

Coordinate System

(i) Global Coordinate System

The coordinate system used to define the points in the entire structure is known as the global coordinate system. For a randomly oriented bar element, $X$ and $Y$ are the global coordinate axes, and the global coordinates for nodes 1 & 2 are $(X_1, Y_1)$ and $(X_2, Y_2)$ respectively.

(ii) Local Coordinate System

For the convenience of deriving element properties in FEM, a separate coordinate system oriented with the element itself is used for each element, known as the local coordinate system. For a bar element, $x'$ and $y'$ are local coordinate axes. The local coordinates for nodes 1 & 2 are $(0,0)$ and $(l,0)$ respectively.

(iii) Natural Coordinate System

The coordinates of nodes of an element defined by a set of dimensionless numbers is known as the natural coordinate system. The value of these dimensionless numbers at a node never exceeds unity, and the sum of all coordinate values at any point is always unity. Hence, it can also be described by weightage of a point in terms of nodal coordinates/weightage.

⇒ Natural Coordinate system can be applied in two ways:

1. Natural Coordinates in the Form of $(L_1, L_2, \ldots)$

For the element with '$n$' nodes, '$n$' variables are used to define the location of a point such that the value of its term at node '$i$' is unity with all other terms' values as zero.

Total weightage (sum of all the terms) at any point within the element is always unity.

For a Two-Noded Bar Element

Location of a point $P$ is defined by coordinates $(L_1, L_2)$ where the value of $L_1$ is 1 at node 1 & zero at node 2. Similarly, the value of $L_2$ is zero at node 1 and 1 at node 2.

1 L₁=1, L₂=0 2 L₁=0, L₂=1 P (L₁, L₂)
Fig: Natural coordinates of a two-noded bar element

For a CST Element (Constant Strain Triangle)

With three nodes, three variables are used in the form of $(L_1, L_2, L_3)$ to locate a point $P$ within the element.

Node$L_1$$L_2$$L_3$
Node 1100
Node 2010
Node 3001

2. Natural Coordinates in the Form of $(\xi, \eta)$

In this method, for a 1-D element, the value of $\xi$ varies from $-1$ to $+1$ from one end to another end of the element.

For a 2-D element, two variables $\xi$ and $\eta$ are used; the value of both varies independently from $-1$ to $+1$ in 2-dimensions.

⇒ The number of variables in this system is independent of the number of nodes the element has, but dependent on whether the element is 1-D, 2-D, or 3-D.
1 ξ = −1 2 ξ = 0 3 ξ = +1 ξ
Fig: Natural coordinates of a 3-noded bar element in terms of $\xi$
η ξ (0,0) 1 (1,−1) 2 (1,1) 4 (−1,1) 5 (−1,−1) 3 (0,1) 6 (0,−1) 7 (1,0) 8 (−1,0)
Fig: Natural coordinates of an 8-noded serendipity element in form of $(\xi, \eta)$

Relation between Natural Coordinate and Cartesian Coordinate

i) For 1-D Beam Element with Natural Coordinate in $(L_1, L_2)$ Form

The total weightage of a point is unity:

$$L_1 + L_2 = 1 \quad \text{--- (i)}$$
$$L_1 x_1 + L_2 x_2 = x \quad \text{--- (ii)}$$

From similar triangles:

$$\frac{L_1}{x_2 - x} = \frac{1}{x_2 - x_1}$$
$$\frac{L_2}{x - x_1} = \frac{1}{x_2 - x_1}$$

Writing equations (i) & (ii) in matrix form:

$$\begin{bmatrix} 1 & 1 \\ x_1 & x_2 \end{bmatrix} \begin{Bmatrix} L_1 \\ L_2 \end{Bmatrix} = \begin{Bmatrix} 1 \\ x \end{Bmatrix}$$
$$\begin{Bmatrix} L_1 \\ L_2 \end{Bmatrix} = \begin{bmatrix} 1 & 1 \\ x_1 & x_2 \end{bmatrix}^{-1} \begin{Bmatrix} 1 \\ x \end{Bmatrix}$$
$$\begin{Bmatrix} L_1 \\ L_2 \end{Bmatrix} = \frac{1}{x_2 - x_1} \begin{bmatrix} x_2 & -1 \\ -x_1 & 1 \end{bmatrix} \begin{Bmatrix} 1 \\ x \end{Bmatrix}$$
$$\begin{Bmatrix} L_1 \\ L_2 \end{Bmatrix} = \frac{1}{l} \begin{Bmatrix} x_2 - x \\ x - x_1 \end{Bmatrix}$$
$$\therefore \quad L_1 = \frac{x_2 - x}{l}, \qquad L_2 = \frac{x - x_1}{l}$$

ii) For Two-Noded Beam with Natural Coordinate in $\xi$ Form ($-1$ to $+1$)

From similar mapping:

$$x_c = \frac{x_1 + x_2}{2}$$
$$x_2 - x_c = x_c - x_1 = \frac{l}{2}$$
$$\frac{\xi}{1} = \frac{x - x_c}{l/2}$$
$$\xi = \frac{2}{l} \left\{ x - \left(\frac{x_1 + x_2}{2}\right) \right\}$$
$$\xi = \frac{2}{l} \left\{ \frac{2x - (x_1 + x_2)}{2} \right\}$$
$$\xi = \frac{2}{l} (x - x_c)$$

Alternatively, to find $x$:

$$x - x_1 = \frac{l}{2}(1 + \xi)$$
$$x = x_1 + \frac{l}{2}(1 + \xi)$$

iii) For CST with Natural Coordinate in $(L_1, L_2, L_3)$ Form

From the definition of natural coordinate and displacement function:

$$L_1 + L_2 + L_3 = 1$$ $$L_1 x_1 + L_2 x_2 + L_3 x_3 = x$$ $$L_1 y_1 + L_2 y_2 + L_3 y_3 = y$$

In matrix form:

$$\begin{bmatrix} 1 & 1 & 1 \\ x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \end{bmatrix} \begin{Bmatrix} L_1 \\ L_2 \\ L_3 \end{Bmatrix} = \begin{Bmatrix} 1 \\ x \\ y \end{Bmatrix}$$

Or,

$$\begin{Bmatrix} L_1 \\ L_2 \\ L_3 \end{Bmatrix} = \begin{bmatrix} 1 & 1 & 1 \\ x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \end{bmatrix}^{-1} \begin{Bmatrix} 1 \\ x \\ y \end{Bmatrix}$$
$$\begin{Bmatrix} L_1 \\ L_2 \\ L_3 \end{Bmatrix} = \frac{1}{2A} \begin{bmatrix} x_2 y_3 - x_3 y_2 & y_2 - y_3 & x_3 - x_2 \\ x_3 y_1 - x_1 y_3 & y_3 - y_1 & x_1 - x_3 \\ x_1 y_2 - x_2 y_1 & y_1 - y_2 & x_2 - x_1 \end{bmatrix} \begin{Bmatrix} 1 \\ x \\ y \end{Bmatrix}$$

Where:

$$2A = \det \begin{vmatrix} 1 & 1 & 1 \\ x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \end{vmatrix} = 2 \times \text{Area of triangle}$$

Or it can be written as:

$$L_1 = \frac{a_1 + b_1 x + c_1 y}{2A}$$ $$L_2 = \frac{a_2 + b_2 x + c_2 y}{2A}$$ $$L_3 = \frac{a_3 + b_3 x + c_3 y}{2A}$$
$$\therefore \quad \begin{Bmatrix} L_1 \\ L_2 \\ L_3 \end{Bmatrix} = \frac{1}{2A} \begin{bmatrix} a_1 & b_1 & c_1 \\ a_2 & b_2 & c_2 \\ a_3 & b_3 & c_3 \end{bmatrix} \begin{Bmatrix} 1 \\ x \\ y \end{Bmatrix}$$

Where:

$$a_1 = x_2 y_3 - x_3 y_2; \quad b_1 = y_2 - y_3; \quad c_1 = x_3 - x_2$$ $$a_2 = x_3 y_1 - x_1 y_3; \quad b_2 = y_3 - y_1; \quad c_2 = x_1 - x_3$$ $$a_3 = x_1 y_2 - x_2 y_1; \quad b_3 = y_1 - y_2; \quad c_3 = x_2 - x_1$$