Chapter 4: Finite Element Method

4.5 — Beam Element
Collected & Prepared By: Er. Sanjay Kumar Sah (Lecturer, Department of Civil Engineering, Pulchowk Campus)

Beam Element Overview

This chapter develops the two-noded Euler-Bernoulli beam element from unit displacement actions, then extends it to the 6-DOF plane frame element. It also collects the fixed-end load formulas, Hermite shape functions, and representative numerical examples for deflection, slope, reactions, shear force, and bending moment.

4 DOFs: $[v_1,\theta_1,v_2,\theta_2]$Symmetric stiffness matrixHermite cubic interpolationEquivalent nodal loads
12$v_1$$v_2$$\theta_1$$\theta_2$Beam element length $L$

Section Map

1. Unit displacement forces2. Beam stiffness matrix3. 6-DOF frame element4. Equivalent nodal loads5. Hermite polynomials6. Shape function examples7. Numerical examples8. Summary and procedure
4
Degrees of freedom for the bending beam element
6
Degrees of freedom for the plane frame element
C1
Continuity achieved by cubic Hermite interpolation

1. Force Development Due to Unit Displacement

The stiffness coefficients of the beam element are found by applying a unit value at one degree of freedom while restraining the remaining three. The resulting shears and end moments at both nodes give one column of the element stiffness matrix.

Euler-Bernoulli Beam Element

A standard two-noded beam element has four degrees of freedom in the order:

$$\{\delta\} = \begin{Bmatrix} v_1 \\ \theta_1 \\ v_2 \\ \theta_2 \end{Bmatrix}$$

The nodal force vector follows the same order as vertical force, moment, vertical force, and moment.

1.1 Unit Transverse Displacement at Node 1

Apply $v_1=1$ upward, keep $\theta_1=v_2=\theta_2=0$.

$\delta=1$$12EI/L^3$$12EI/L^3$$6EI/L^2$$6EI/L^2$
$$K_{11}=\frac{12EI}{L^3},\quad K_{21}=\frac{6EI}{L^2},\quad K_{31}=-\frac{12EI}{L^3},\quad K_{41}=\frac{6EI}{L^2}$$

1.2 Unit Transverse Displacement at Node 2

Apply $v_2=1$ upward, keep $v_1=\theta_1=\theta_2=0$.

$\delta=1$$12EI/L^3$$12EI/L^3$$6EI/L^2$$6EI/L^2$
$$K_{13}=-\frac{12EI}{L^3},\quad K_{23}=-\frac{6EI}{L^2},\quad K_{33}=\frac{12EI}{L^3},\quad K_{43}=-\frac{6EI}{L^2}$$

1.3 Unit Rotation at Node 1

Apply $\theta_1=1$, keep $v_1=v_2=\theta_2=0$.

$6EI/L^2$$6EI/L^2$$4EI/L$$2EI/L$
$$K_{12}=\frac{6EI}{L^2},\quad K_{22}=\frac{4EI}{L},\quad K_{32}=-\frac{6EI}{L^2},\quad K_{42}=\frac{2EI}{L}$$

1.4 Unit Rotation at Node 2

Apply $\theta_2=1$, keep $v_1=\theta_1=v_2=0$.

$6EI/L^2$$6EI/L^2$$2EI/L$$4EI/L$
$$K_{14}=\frac{6EI}{L^2},\quad K_{24}=\frac{2EI}{L},\quad K_{34}=-\frac{6EI}{L^2},\quad K_{44}=\frac{4EI}{L}$$

2. Stiffness Matrix of Beam Element

Placing the four force columns together gives the complete element stiffness matrix for beam bending.

DOFDOF 1 $(v_1)$DOF 2 $(\theta_1)$DOF 3 $(v_2)$DOF 4 $(\theta_2)$
DOF 1 $(v_1)$$12EI/L^3$$6EI/L^2$$-12EI/L^3$$6EI/L^2$
DOF 2 $(\theta_1)$$6EI/L^2$$4EI/L$$-6EI/L^2$$2EI/L$
DOF 3 $(v_2)$$-12EI/L^3$$-6EI/L^2$$12EI/L^3$$-6EI/L^2$
DOF 4 $(\theta_2)$$6EI/L^2$$2EI/L$$-6EI/L^2$$4EI/L$
All diagonal terms are positive. The matrix is symmetric because $K_{ij}=K_{ji}$.
$$[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}$$

Sign structure

The coupling between the two nodes appears through opposite signs in the off-diagonal force terms, reflecting action-reaction behavior of the element.

Physical interpretation

The $12EI/L^3$ terms govern transverse shear stiffness, while the $4EI/L$ and $2EI/L$ terms govern rotational stiffness at the two ends.

3. Beam Element Including Axial Force and Displacement (Frame Element)

When axial deformation is also included, each node carries three degrees of freedom: axial displacement $u$, transverse displacement $v$, and rotation $\theta$. The result is a 6-DOF plane frame element.

12DOF 1: $u_1$DOF 4: $u_2$DOF 2: $v_1$DOF 5: $v_2$DOF 3: $\theta_1$DOF 6: $\theta_2$
$$[K]=\begin{bmatrix}EA/L & 0 & 0 & -EA/L & 0 & 0 \\ 0 & 12EI/L^3 & 6EI/L^2 & 0 & -12EI/L^3 & 6EI/L^2 \\ 0 & 6EI/L^2 & 4EI/L & 0 & -6EI/L^2 & 2EI/L \\ -EA/L & 0 & 0 & EA/L & 0 & 0 \\ 0 & -12EI/L^3 & -6EI/L^2 & 0 & 12EI/L^3 & -6EI/L^2 \\ 0 & 6EI/L^2 & 2EI/L & 0 & -6EI/L^2 & 4EI/L\end{bmatrix}$$

4. Equivalent Nodal Loads

Distributed loads and intermediate concentrated loads are converted into statically equivalent nodal shears and nodal moments before assembly. In beam-element formulation, the equivalent nodal load vector is the negative of the fixed-end force vector.

4.1 Central Point Load $P$

$P$AB
$$F_A=F_B=P/2,\qquad M_A=M_B=PL/8$$

4.2 Off-Centre Point Load

$P$$a$$b$
$$F_A=\frac{Pb^2(L+2a)}{L^3},\quad F_B=\frac{Pa^2(L+2b)}{L^3}$$ $$M_A=\frac{Pab^2}{L^2},\quad M_B=\frac{Pa^2b}{L^2}$$

4.3 Uniformly Distributed Load

$\omega$
$$F_A=F_B=\omega L/2,\qquad M_A=M_B=\omega L^2/12$$

4.4 Triangular Load, Left Maximum

$$F_A=\frac{7\omega L}{20},\quad M_A=\frac{\omega L^2}{20}$$ $$F_B=\frac{3\omega L}{20},\quad M_B=\frac{\omega L^2}{30}$$

4.5 Symmetric Triangular Load

$$F_A=F_B=\omega L/4,\qquad M_A=M_B=5\omega L^2/96$$

4.6 Applied Moment at an Interior Point

$M$$a$$b$
$$F_A=\frac{6Mab}{L^3}\;\text{up},\qquad F_B=\frac{6Mab}{L^3}\;\text{down}$$ $$M_A=\frac{Mb(2a-b)}{L^2},\qquad M_B=\frac{Ma(2b-a)}{L^2}$$
Loading$F_{left}$$M_{left}$$F_{right}$$M_{right}$
UDL $\omega$ over full span$\omega L/2$$\omega L^2/12$$\omega L/2$$\omega L^2/12$
$P$ at midspan$P/2$$PL/8$$P/2$$PL/8$
$P$ at $a,b$$Pb^2(L+2a)/L^3$$Pab^2/L^2$$Pa^2(L+2b)/L^3$$Pa^2b/L^2$
Triangular load, left maximum$7\omega L/20$$\omega L^2/20$$3\omega L/20$$\omega L^2/30$
Symmetric triangular load$\omega L/4$$5\omega L^2/96$$\omega L/4$$5\omega L^2/96$
Applied moment $M$ at $a,b$$6Mab/L^3$ upward$Mb(2a-b)/L^2$$6Mab/L^3$ downward$Ma(2b-a)/L^2$

5. Hermite Polynomials as Shape Functions

5.1 Introduction

Hermite polynomials are natural beam interpolation functions because the beam element has to interpolate both transverse displacement and slope. A one-dimensional Hermite polynomial of order $n$ has order $2n+1$.

Their key property is that the polynomial and its required derivatives take values 0 or 1 at the end points, which matches the beam nodal requirements exactly.

5.2 Two-Noded Beam Element Setup

12$P(x)$$l=x_2-x_1$$s=0$$s=1$$(v_1,\theta_1)$$(v_2,\theta_2)$

Using the natural coordinate $s=(x-x_1)/l$, the first-order Hermite functions for the beam element are:

$$H^1_{01}(s)=1-3s^2+2s^3$$ $$H^1_{11}(s)=ls(s-1)^2$$ $$H^1_{02}(s)=s^2(3-2s)$$ $$H^1_{12}(s)=ls^2(s-1)$$

5.3 Displacement Model

$$v=[N]\{\delta\}=\begin{bmatrix}N_1 & N_2 & N_3 & N_4\end{bmatrix}\begin{Bmatrix}v_1 \\ \theta_1 \\ v_2 \\ \theta_2\end{Bmatrix}$$ $$\theta=\frac{\partial v}{\partial x}=\begin{bmatrix}\partial N_1/\partial x & \partial N_2/\partial x & \partial N_3/\partial x & \partial N_4/\partial x\end{bmatrix}\begin{Bmatrix}v_1 \\ \theta_1 \\ v_2 \\ \theta_2\end{Bmatrix}$$

Shape function behavior

01$N_1$$N_3$$N_2$$N_4$

$N_1$ and $N_3$ control nodal displacement. $N_2$ and $N_4$ control nodal slope contributions.

5.4 Verification at nodes

NodeConditionValue requirement
Node 1$s=0$$N_1=1$, $N_2=N_3=N_4=0$ and $\partial N_2/\partial x=1$
Node 2$s=1$$N_3=1$, $N_1=N_2=N_4=0$ and $\partial N_4/\partial x=1$

The element therefore satisfies the exact interpolation requirements for both displacement and rotation degrees of freedom.

6. Shape Function Examples

Example 1: Shape Functions in Cartesian Coordinate System
Find the shape functions of a beam element in Cartesian coordinates with $x_1=0$ and $x_2=l$.
Solution
$$v=\alpha_1+\alpha_2x+\alpha_3x^2+\alpha_4x^3$$ $$\theta=\frac{\partial v}{\partial x}=\alpha_2+2\alpha_3x+3\alpha_4x^2$$
$$v_1=\alpha_1,\qquad \theta_1=\alpha_2$$ $$v_2=\alpha_1+\alpha_2l+\alpha_3l^2+\alpha_4l^3$$ $$\theta_2=\alpha_2+2\alpha_3l+3\alpha_4l^2$$

After inverting the nodal-constant relation and substituting back:

$$N_1=1-\frac{3x^2}{l^2}+\frac{2x^3}{l^3}$$ $$N_2=x-\frac{2x^2}{l}+\frac{x^3}{l^2}$$ $$N_3=\frac{3x^2}{l^2}-\frac{2x^3}{l^3}$$ $$N_4=-\frac{x^2}{l}+\frac{x^3}{l^2}$$
Example 2: Shape Functions in Natural Coordinates $s\in[0,1]$
Solution

Using $s=x/l$:

$$N_1=1-3s^2+2s^3$$ $$N_2=ls(s-1)^2$$ $$N_3=s^2(3-2s)$$ $$N_4=ls^2(s-1)$$
$$\frac{\partial N}{\partial x}=\frac{1}{l}\frac{\partial N}{\partial s}$$
Example 3: Shape Functions in Natural Coordinates $\xi\in[-1,1]$
Solution

Using the coordinate mapping $\xi=\frac{2}{l}(x-l/2)$ so that $x/l=(\xi+1)/2$:

$$N_1=\frac{(1-\xi)^2(2+\xi)}{4}$$ $$N_2=\frac{l(1-\xi)^2(\xi+1)}{8}$$ $$N_3=\frac{(1+\xi)^2(2-\xi)}{4}$$ $$N_4=\frac{l(1+\xi)^2(\xi-1)}{8}$$
$$N_i=\frac{(1+\xi\xi_i)^2(2-\xi\xi_i)}{4},\qquad H_i=\frac{l}{8}(1+\xi\xi_i)^2(\xi-\xi_i)$$ $$\frac{\partial N}{\partial x}=\frac{2}{l}\frac{\partial N}{\partial \xi}$$

7. Numerical Examples

Numerical Example 1: Slope and Deflection at Midpoint of Cantilever Beam
Cantilever beam with UDL $5\,\text{kN/m}$, $L=2\,\text{m}$, and $EI=2\times10^{12}\,\text{N mm}^2 = 2\times10^3\,\text{kN m}^2$.
$\omega=5\,\text{kN/m}$$2\,\text{m}$
Step 1: Equivalent nodal loads
$$F=\omega L/2=5\,\text{kN},\qquad M=\omega L^2/12=1.67\,\text{kN m}$$ $$\{F\}=\begin{Bmatrix}R_1-5 \\ M_1-1.67 \\ -5 \\ 1.67\end{Bmatrix}$$
Step 2: Displacement vector
$$\{u\}=\begin{Bmatrix}0 \\ 0 \\ v_2 \\ \theta_2\end{Bmatrix}$$
Step 3: Reduced system
$$\begin{Bmatrix}-5 \\ 1.67\end{Bmatrix}=1000\begin{bmatrix}3 & -3 \\ -3 & 4\end{bmatrix}\begin{Bmatrix}v_2 \\ \theta_2\end{Bmatrix}$$ $$v_2=-0.005\,\text{m},\qquad \theta_2=-0.00333\,\text{rad}$$
Step 4: Reactions
$$R_1=10\,\text{kN},\qquad M_1=10\,\text{kN m}$$
Step 5: Midpoint deflection and slope
$$v(1)=N_3(1)v_2+N_4(1)\theta_2=-0.001666\,\text{m}\approx -1.67\,\text{mm}$$ $$\theta(1)=\frac{\partial N_3}{\partial x}v_2+\frac{\partial N_4}{\partial x}\theta_2=-0.002916\,\text{rad}\approx -0.167^\circ$$

Free-end deflection is 5 mm downward. Midpoint deflection is 1.67 mm downward.

Numerical Example 2: Propped Cantilever, Two-Element Model
UDL $15\,\text{kN/m}$ on a 10 m beam divided into a 6 m element and a 4 m element. $EI=10000\,\text{kN m}^2$.
$$\{F\}=\begin{Bmatrix}R_A-45 \\ M_A-45 \\ -75 \\ 25 \\ R_C-30 \\ 20\end{Bmatrix},\qquad \{u\}=\begin{Bmatrix}0 \\ 0 \\ \delta_B \\ \theta_B \\ 0 \\ \theta_C\end{Bmatrix}$$
$$\begin{Bmatrix}-75 \\ 25 \\ 20\end{Bmatrix}=\frac{10000}{64}\begin{bmatrix}15.556 & 13.333 & 24 \\ 13.333 & 106.667 & 32 \\ 24 & 32 & 64\end{bmatrix}\begin{Bmatrix}\delta_B \\ \theta_B \\ \theta_C\end{Bmatrix}$$
$$\delta_B=-0.0810\,\text{m},\qquad \theta_B=0.00225\,\text{rad},\qquad \theta_C=0.03125\,\text{rad}$$ $$R_A=93.75\,\text{kN},\qquad M_A=187.5\,\text{kN m},\qquad R_C=56.25\,\text{kN}$$
A single element over the full 10 m span reproduces $\theta_C$, but the two-element model gives better internal-point results.
Numerical Example 3: Continuous Beam, Reactions, SFD, BMD, and Deflections
Beam A-B-C with a 30 kN point load on AB, a 40 kN m moment at B, and a 20 kN/m UDL on BC. $EI=10000\,\text{kN m}^2$.
$$\{F\}=\begin{Bmatrix}R_A-20.508 \\ M_A-35.156 \\ R_B-69.492 \\ -78.096 \\ R_C-60 \\ 60\end{Bmatrix},\qquad \{u\}=\begin{Bmatrix}0 \\ 0 \\ 0 \\ \theta_B \\ 0 \\ \theta_C\end{Bmatrix}$$
$$\begin{Bmatrix}-78.096 \\ 60\end{Bmatrix}=\frac{1250}{9}\begin{bmatrix}84 & 24 \\ 24 & 48\end{bmatrix}\begin{Bmatrix}\theta_B \\ \theta_C\end{Bmatrix}$$ $$\theta_B=-0.0108906\,\text{rad},\qquad \theta_C=0.0144453\,\text{rad}$$
$$R_A=10.298\,\text{kN},\quad M_A=7.930\,\text{kN m},\quad R_B=85.626\,\text{kN},\quad R_C=54.076\,\text{kN}$$ $$R_{B1}=19.702\,\text{kN},\quad M_{B1}=75.547\,\text{kN m},\quad R_{B2}=65.924\,\text{kN},\quad M_{B2}=35.547\,\text{kN m}$$

Mid-span deflections from the Hermite interpolation are $+10.89$ mm upward in span AB and approximately $-19$ mm downward in span BC.

Numerical Example 4: Deflection and Slope Under a Point Load
Two-element beam with a 20 kN point load in element 1 and a 10 kN/m UDL in element 2. $EI=10^3\,\text{kN m}^2$ and each element is 4 m long.
$$\{F\}=\begin{Bmatrix}R_1-10 \\ M_1-10 \\ R_2-30 \\ -8.333 \\ R_3-20 \\ 13.333\end{Bmatrix},\qquad \{u\}=\begin{Bmatrix}0 \\ 0 \\ 0 \\ \theta_2 \\ 0 \\ \theta_3\end{Bmatrix}$$
$$\begin{Bmatrix}-8.333 \\ 13.333\end{Bmatrix}=62.5\begin{bmatrix}32 & 8 \\ 8 & 16\end{bmatrix}\begin{Bmatrix}\theta_2 \\ \theta_3\end{Bmatrix}$$ $$\theta_2=-0.005714\,\text{rad},\qquad \theta_3=0.016190\,\text{rad}$$
$$v(2)=N_4(2)\theta_2=0.002857\,\text{m}=2.857\,\text{mm upward}$$ $$\theta(2)=\left.\frac{\partial N_4}{\partial x}\right|_{x=2}\theta_2=0.001429\,\text{rad}\approx0.0819^\circ$$

8. Summary of Key Formulas and Procedures

8.1 Beam Element Stiffness Matrix

$$[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}$$

8.2 Shape Functions in Three Coordinate Systems

FunctionCartesian $x$Natural $s$Natural $\xi$
$N_1$$1-3x^2/l^2+2x^3/l^3$$1-3s^2+2s^3$$(1-\xi)^2(2+\xi)/4$
$N_2$$x-2x^2/l+x^3/l^2$$ls(s-1)^2$$l(1-\xi)^2(\xi+1)/8$
$N_3$$3x^2/l^2-2x^3/l^3$$s^2(3-2s)$$(1+\xi)^2(2-\xi)/4$
$N_4$$-x^2/l+x^3/l^2$$ls^2(s-1)$$l(1+\xi)^2(\xi-1)/8$

8.3 Key Properties of Stiffness Matrix

8.4 Step-by-Step FEM Procedure for Beam Analysis

  1. Discretize the beam into elements and label nodes and degrees of freedom.
  2. Convert distributed or intermediate loads into equivalent nodal loads.
  3. Form the element stiffness matrices $[k]^e$.
  4. Assemble the global stiffness matrix $[K]$ and the global load vector $\{F\}$.
  5. Apply essential boundary conditions by eliminating restrained displacements.
  6. Solve the reduced system for unknown nodal displacements and rotations.
  7. Back-substitute to obtain reactions.
  8. Use element equations to obtain end shears and end moments.
  9. Use $v(x)=[N]\{\delta\}^e$ and $\theta(x)=[\partial N/\partial x]\{\delta\}^e$ for any interior point.

8.5 Accuracy and Convergence

End of beam element notes. Content organized from the supplied lecture-note material for a readable HTML textbook format.
← Previous: 4.4 Shape Functions