Structural Engineering Dynamics

CHAPTER 4: FINITE ELEMENT METHOD

4.7 INTERPOLATION FUNCTIONS (SHAPE FUNCTIONS)

PAST EXAM QUESTIONS ADDRESSED IN THIS SECTION:

Define shape function in finite analysis and write down the properties of shape function. [3] (2080 Chaitra, Q9)
Derive the expressions for Hermite shape function used for the interpolation in beam elements. [7-8] (2080 Chaitra, Q9; 2079 Jestha, Q10/12; 2077 Chaitra, Q4)

Definition and Properties

In Finite Element Analysis, a Shape Function (or interpolation function) is a continuous mathematical polynomial used to approximate the exact solution of a field variable (such as displacement, temperature, or stress) at any point inside an element, based solely on the known discrete values of that variable at the element's nodes.

Mathematically, the displacement $ u(x) $ within an element is expressed as $ u(x) = \sum N_i u_i $, where $ N_i $ is the shape function associated with node $ i $, and $ u_i $ is the nodal displacement.

Fundamental Properties of Shape Functions:

Kronecker Delta Property: A shape function $ N_i $ has a value of unity (1.0) at its home node $ i $, and a value of zero (0.0) at all other nodes. mathematically: $ N_i(x_j) = \delta_{ij} $.
Partition of Unity: The sum of all shape functions evaluated at any arbitrary point within the element is exactly equal to one. $ \sum N_i = 1 $. This ensures that a rigid body translation is accurately modeled without inducing fictitious strains.
Polynomial Order: The degree of the chosen polynomial strictly depends on the number of nodes and degrees of freedom (DOF) assigned to the element.

Derivation of Hermite Shape Functions for Beam Elements

Unlike a simple truss element which only displaces axially, a Euler-Bernoulli beam element undergoes both transverse deflection ($ v $) and rotation ($ \theta = \frac{dv}{dx} $). A 2-node beam element therefore has 4 degrees of freedom: $ v_1, \theta_1 $ at node 1, and $ v_2, \theta_2 $ at node 2.

Figure 4.1: A 2-node beam element showing 4 degrees of freedom.

Because there are 4 unknowns, we must assume a cubic polynomial displacement function. To simplify the algebra, we introduce a non-dimensional natural coordinate $ \xi = \frac{x}{L} $. The element spans from $ \xi = 0 $ to $ \xi = 1 $.

$$ v(\xi) = a_0 + a_1 \xi + a_2 \xi^2 + a_3 \xi^3 $$

The slope (rotation) is the derivative with respect to $ x $. Using the chain rule, $ \theta = \frac{dv}{dx} = \frac{dv}{d\xi} \frac{d\xi}{dx} = \frac{1}{L} \frac{dv}{d\xi} $. Therefore:

$$ \frac{dv}{d\xi} = L \cdot \theta(\xi) = a_1 + 2a_2 \xi + 3a_3 \xi^2 $$

We apply the four physical boundary conditions to solve for constants $ a_0, a_1, a_2, a_3 $:

At Node 1 ($ \xi = 0 $): $ v(0) = v_1 \implies a_0 = v_1 $
At Node 1 ($ \xi = 0 $): $ \frac{dv}{d\xi}(0) = L \theta_1 \implies a_1 = L \theta_1 $
At Node 2 ($ \xi = 1 $): $ v(1) = v_2 \implies v_1 + L\theta_1 + a_2 + a_3 = v_2 $
At Node 2 ($ \xi = 1 $): $ \frac{dv}{d\xi}(1) = L \theta_2 \implies L\theta_1 + 2a_2 + 3a_3 = L\theta_2 $

Solving equations (3) and (4) simultaneously for $ a_2 $ and $ a_3 $:

$$ \begin{aligned} a_2 &= 3(v_2 - v_1) - L(2\theta_1 + \theta_2) \\ a_3 &= 2(v_1 - v_2) + L(\theta_1 + \theta_2) \end{aligned} $$

Substituting $ a_0, a_1, a_2, a_3 $ back into the original cubic equation $ v(\xi) $, and factoring out the nodal variables $ v_1, \theta_1, v_2, \theta_2 $, we get:

$$ \begin{aligned} v(\xi) &= (1 - 3\xi^2 + 2\xi^3)v_1 \\ &\quad + L(\xi - 2\xi^2 + \xi^3)\theta_1 \\ &\quad + (3\xi^2 - 2\xi^3)v_2 \\ &\quad + L(-\xi^2 + \xi^3)\theta_2 \end{aligned} $$

This perfectly defines the displacement in terms of the four nodal degrees of freedom: $ v(\xi) = N_1 v_1 + N_2 \theta_1 + N_3 v_2 + N_4 \theta_2 $. The resulting Hermite Shape Functions are:

$$ \begin{aligned} N_1 &= 1 - 3\xi^2 + 2\xi^3 \\ N_2 &= L(\xi - 2\xi^2 + \xi^3) \\ N_3 &= 3\xi^2 - 2\xi^3 \\ N_4 &= L(-\xi^2 + \xi^3) \end{aligned} $$

(Note: These functions guarantee the continuity of both displacement and slope across element boundaries, which is a strictly required condition (C1 continuity) for beam bending problems.)

4.8 CONSTANT STRAIN TRIANGLE (CST)

PAST EXAM QUESTIONS ADDRESSED IN THIS SECTION:

Derive the relation of strain-Displacement [B] matrix for constant strain triangle. [6] (2073 Magh, Q9)

Derivation of the Strain-Displacement [B] Matrix

The Constant Strain Triangle (CST) is the simplest 2D solid element. It is a 3-node triangle where each node has two degrees of freedom (horizontal displacement $ u $ and vertical displacement $ v $). Thus, there are 6 total DOF for the element.

Figure 4.2: Constant Strain Triangle showing nodal coordinates and displacements.

The linear displacement field within the triangle is assumed as:

$$ \begin{aligned} u(x,y) &= N_1 u_1 + N_2 u_2 + N_3 u_3 \\ v(x,y) &= N_1 v_1 + N_2 v_2 + N_3 v_3 \end{aligned} $$

Where the linear shape functions $ N_i $ are defined geometrically in terms of the element's area ($ A $):

$$ \begin{aligned} N_1 &= \frac{1}{2A}(a_1 + b_1 x + c_1 y) \\ N_2 &= \frac{1}{2A}(a_2 + b_2 x + c_2 y) \\ N_3 &= \frac{1}{2A}(a_3 + b_3 x + c_3 y) \end{aligned} $$

The constants are derived from the coordinate geometry, where $ b_1 = y_2 - y_3 $, $ c_1 = x_3 - x_2 $, with cyclic permutations for nodes 2 and 3. Notice that differentiating these shape functions yields pure constants: $ \frac{\partial N_i}{\partial x} = \frac{b_i}{2A} $ and $ \frac{\partial N_i}{\partial y} = \frac{c_i}{2A} $.

We compute the 2D strains (normal strains and shear strain) using the derivatives of the displacement field:

$$ \begin{aligned} \{\epsilon\} &= \begin{Bmatrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{Bmatrix} \\ &= \begin{Bmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{Bmatrix} \end{aligned} $$

Substituting the shape functions into the strain definition:

$$ \begin{aligned} \epsilon_x &= \frac{\partial}{\partial x}(N_1 u_1 + N_2 u_2 + N_3 u_3) \\ &= \frac{1}{2A}(b_1 u_1 + b_2 u_2 + b_3 u_3) \\ \epsilon_y &= \frac{\partial}{\partial y}(N_1 v_1 + N_2 v_2 + N_3 v_3) \\ &= \frac{1}{2A}(c_1 v_1 + c_2 v_2 + c_3 v_3) \\ \gamma_{xy} &= \left[ \frac{\partial}{\partial y}(N_1 u_1 + N_2 u_2 + N_3 u_3) \right] + \left[ \frac{\partial}{\partial x}(N_1 v_1 + N_2 v_2 + N_3 v_3) \right] \\ &= \frac{1}{2A}(c_1 u_1 + c_2 u_2 + c_3 u_3) + \frac{1}{2A}(b_1 v_1 + b_2 v_2 + b_3 v_3) \end{aligned} $$

Because these derivatives yield constants (no $x$ or $y$ variables remain), the strain is perfectly uniform everywhere inside the element—hence the name "Constant Strain Triangle". Factoring out the nodal displacement vector $ \{d\} $, we construct the Strain-Displacement matrix $ [B] $:

$$ \begin{aligned} \{\epsilon\} &= [B]\{d\} \\ \begin{Bmatrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{Bmatrix} &= \frac{1}{2A} \begin{bmatrix} b_1 & 0 & b_2 & 0 & b_3 & 0 \\ 0 & c_1 & 0 & c_2 & 0 & c_3 \\ c_1 & b_1 & c_2 & b_2 & c_3 & b_3 \end{bmatrix} \begin{Bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \end{Bmatrix} \end{aligned} $$

4.9 ISOPARAMETRIC FORMULATION

PAST EXAM QUESTIONS ADDRESSED IN THIS SECTION:

What is Jacobian matrix? [2] (2079 Chaitra, Q10)
Define isoparametric, super parametric, and sub parametric elements. [3] (2079 Jestha, Q10/12)
Explain / What is Iso-parametric formulation / element? [1-2] (2076 Bhadra, Q9; 2074 Bhadra, Q9; 2072 Magh, Q5)

Isoparametric Formulation

In the analysis of complex real-world structures, meshes must accommodate curved and highly irregular boundaries. Performing complex calculus (integration) over skewed, distorted shapes is mathematically tedious and computationally inefficient.

Isoparametric Formulation solves this by using the exact same shape functions ($ N_i $) to map both the geometric coordinates of the element and the field variables (displacements) within the element. By doing so, any irregularly shaped element in the global coordinate system $ (x, y) $ is mathematically transformed (mapped) into a perfectly normalized, geometrically simple "master element" in a natural coordinate system $ (\xi, \eta) $ spanning from -1 to 1. All numerical integration is performed rapidly on this simple master element.

Figure 4.3: Isoparametric mapping of a distorted physical quad into a perfect computational square.

Element Classifications

Isoparametric Element: The number of nodes used to define the element's geometry is exactly equal to the number of nodes used to define the displacement field. (e.g., A 4-node quad where geometry and displacement are both linear).
Super-parametric Element: More nodes are used to define the geometry than are used for displacements. (Used to accurately model highly curved boundaries, even if the assumed strain state inside is kept simple/linear).
Sub-parametric Element: Fewer nodes are used to define the geometry than are used for displacements. (Used when the physical boundaries are perfectly straight, but a complex higher-order quadratic strain field is required inside the element).

The Jacobian Matrix [J]

The Jacobian matrix is a vital mathematical tensor of partial derivatives that links the spatial derivatives of the global coordinate system $ (x, y) $ to the spatial derivatives of the natural, normalized coordinate system $ (\xi, \eta) $.

By the chain rule of differentiation, derivatives with respect to natural coordinates are tied to global coordinates via:

$$ \begin{aligned} \begin{Bmatrix} \frac{\partial N_i}{\partial \xi} \\ \frac{\partial N_i}{\partial \eta} \end{Bmatrix} &= \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{bmatrix} \begin{Bmatrix} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \end{Bmatrix} \\ &= [J] \begin{Bmatrix} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \end{Bmatrix} \end{aligned} $$

The matrix $ [J] $ is the Jacobian matrix. Its determinant, $ |J| $, is fundamentally necessary to map elemental areas from the physical space into the computational space during the numerical stiffness integration (Gauss Quadrature), such that:

$$ dx \cdot dy = |J| \cdot d\xi \cdot d\eta $$

If $ |J| \le 0 $ anywhere inside the element, it physically means the element mesh is invalid (e.g., folded over on itself, or overly distorted), and the computer algorithm will fail.