Interactive Lecture
Understanding the Finite Element Method Through 1D Bar Elements
by Dr Arnaldo Delli Carri
The Finite Element Method (FEM) is a numerical technique for finding approximate solutions to boundary value problems for partial differential equations. It subdivides a large problem into smaller, simpler parts called finite elements, which are connected at points called nodes.
FEM transforms complex continuous problems into discrete algebraic equations that can be solved computationally. The key is to approximate the solution using simple functions (shape functions) over small regions (elements), then assemble these into a global system.
The simplest finite element is the 1D bar element, which models axial deformation under tension or compression. Consider a uniform bar of length \(L\), cross-sectional area \(A\), and elastic modulus \(E\).
Governing Differential Equation:
\[\frac{d}{dx}\left(EA\frac{du}{dx}\right) + f(x) = 0\]
where \(u(x)\) is the axial displacement and \(f(x)\) is the distributed body force per unit length.
For a linear bar element with two nodes (1 and 2), we approximate the displacement field using linear shape functions:
\[u(x) = N_1(x)u_1 + N_2(x)u_2\]
where the shape functions are:
\[N_1(x) = 1 - \frac{x}{L}, \quad N_2(x) = \frac{x}{L}\]
Properties of shape functions:
Using the principle of minimum potential energy or Galerkin's method, we derive the element stiffness matrix:
\[[k^e] = \frac{EA}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\]
The element equilibrium equation is:
\[[k^e]\{u^e\} = \{f^e\}\]
where \(\{u^e\} = \begin{bmatrix} u_1 \\ u_2 \end{bmatrix}\) and \(\{f^e\} = \begin{bmatrix} f_1 \\ f_2 \end{bmatrix}\)
For multiple elements, we assemble the global stiffness matrix by adding contributions from elements sharing common nodes. This is called the direct stiffness method.
Assembly Rule: The global stiffness matrix [K] is assembled by adding element stiffness matrices at positions corresponding to their global degrees of freedom. Elements sharing a node contribute to the same entry in [K], ensuring compatibility and equilibrium.
Global System:
\[[K]\{U\} = \{F\}\]
After applying boundary conditions and solving, we obtain nodal displacements \(\{U\}\).
Post-Processing: Element strains and stresses are calculated from:
\[\varepsilon = \frac{du}{dx} = \frac{u_2 - u_1}{L}\]
\[\sigma = E\varepsilon = \frac{E(u_2 - u_1)}{L}\]
As we increase the number of elements (mesh refinement), the FEM solution converges to the exact solution. This is guaranteed if:
This simulation compares the exact analytical solution with the FEM solution for a stepped shaft under axial loading. The shaft has two sections with different diameters.
The Finite Element Method is one of the most widely used numerical techniques in engineering:
| Concept | Equation | Description |
|---|---|---|
| Governing Equation (1D Bar) | \(\frac{d}{dx}\left(EA\frac{du}{dx}\right) + f(x) = 0\) | Equilibrium equation for axial bar |
| Shape Functions (Linear) | \(N_1 = 1 - \frac{x}{L}, \quad N_2 = \frac{x}{L}\) | Linear interpolation functions |
| Displacement Approximation | \(u(x) = N_1(x)u_1 + N_2(x)u_2\) | FEM approximation of displacement |
| Element Stiffness Matrix | \([k^e] = \frac{EA}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\) | Bar element stiffness |
| Global System | \([K]\{U\} = \{F\}\) | Assembled system of equations |
| Element Strain | \(\varepsilon = \frac{u_2 - u_1}{L}\) | Constant strain in element |
| Element Stress | \(\sigma = E\varepsilon = \frac{E(u_2 - u_1)}{L}\) | Constant stress in element |
A bar of total length 2 m is divided into two elements. Element 1 has \(E_1 = 200\) GPa, \(A_1 = 100\) mm², and length 1 m. Element 2 has \(E_2 = 100\) GPa, \(A_2 = 50\) mm², and length 1 m. A force of 10 kN is applied at the free end. Calculate the nodal displacements and element stresses.
Show Solution
Step 1: Element stiffnesses
\(k_1 = \frac{E_1 A_1}{L_1} = \frac{200 \times 10^9 \times 100 \times 10^{-6}}{1} = 20 \times 10^6\) N/m
\(k_2 = \frac{E_2 A_2}{L_2} = \frac{100 \times 10^9 \times 50 \times 10^{-6}}{1} = 5 \times 10^6\) N/m
Step 2: Assemble global stiffness matrix
\([K] = 10^6 \begin{bmatrix} 20 & -20 & 0 \\ -20 & 25 & -5 \\ 0 & -5 & 5 \end{bmatrix}\) N/m
Step 3: Apply boundary conditions
\(u_1 = 0\) (fixed), \(F_3 = 10000\) N
Step 4: Solve reduced system
\(10^6 \begin{bmatrix} 25 & -5 \\ -5 & 5 \end{bmatrix} \begin{bmatrix} u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 10000 \end{bmatrix}\)
Step 5: Nodal displacements and element stresses
\(u_2 = 0.5\) mm, \(u_3 = 2.5\) mm
\(\sigma_1 = \frac{E_1(u_2 - u_1)}{L_1} = \frac{200 \times 10^9 \times 0.5 \times 10^{-3}}{1} = 100\) MPa
\(\sigma_2 = \frac{E_2(u_3 - u_2)}{L_2} = \frac{100 \times 10^9 \times 2.0 \times 10^{-3}}{1} = 200\) MPa
Answer: \(u_2 = 0.5\) mm, \(u_3 = 2.5\) mm, \(\sigma_1 = 100\) MPa, \(\sigma_2 = 200\) MPa.
For a uniform bar with \(E = 200\) GPa, \(A = 100\) mm², \(L = 1\) m, and \(P = 50\) kN, the exact tip displacement is \(u_{exact} = PL/(EA) = 2.5\) mm. Calculate the FEM solution using 1, 2, and 4 elements and discuss convergence.
Show Solution
1 Element:
\(k_e = \frac{EA}{L} = \frac{200 \times 10^9 \times 100 \times 10^{-6}}{1} = 20 \times 10^6\) N/m
\(u_{tip} = \frac{P}{k_e} = \frac{50000}{20 \times 10^6} = 2.5\) mm
Error = 0% (exact for uniform bar!)
2 Elements:
\(k_e = \frac{EA}{L/2} = 40 \times 10^6\) N/m
\(u_{tip} = 2.5\) mm
Error = 0%
4 Elements:
\(k_e = \frac{EA}{L/4} = 80 \times 10^6\) N/m
\(u_{tip} = 2.5\) mm
Error = 0%
Answer: \(u_{tip} = 2.5\) mm with 0% error for 1, 2, and 4 elements.
Conclusion: For the uniform bar under constant axial load, the linear element gives the exact solution regardless of the number of elements because the exact solution is itself linear. For non-uniform bars or varying loads, more elements improve accuracy.
Verify that the linear shape functions \(N_1(x) = 1 - x/L\) and \(N_2(x) = x/L\) satisfy the essential properties: (a) Kronecker delta property, (b) partition of unity, and (c) can represent constant strain.
Show Solution
(a) Kronecker Delta Property:
At node 1 (\(x = 0\)): \(N_1(0) = 1\), \(N_2(0) = 0\)
At node 2 (\(x = L\)): \(N_1(L) = 0\), \(N_2(L) = 1\)
✓ Each shape function equals 1 at its own node and 0 at other nodes.
(b) Partition of Unity:
\(N_1(x) + N_2(x) = \left(1 - \frac{x}{L}\right) + \frac{x}{L} = 1\) for all \(x\)
✓ Shape functions sum to unity everywhere, ensuring rigid body motion can be represented.
(c) Constant Strain Representation:
If \(u_1 = 0\) and \(u_2 = \varepsilon_0 L\) (constant strain \(\varepsilon_0\)):
\(u(x) = N_1(x) \cdot 0 + N_2(x) \cdot \varepsilon_0 L = \frac{x}{L} \cdot \varepsilon_0 L = \varepsilon_0 x\)
\(\varepsilon = \frac{du}{dx} = \varepsilon_0\) (constant)
✓ Linear shape functions can exactly represent constant strain states.
Answer: All three properties hold — (a) Kronecker delta, (b) partition of unity, and (c) exact constant-strain representation.