The weak form of linear elasticity, on one page

The finite element method does not solve the equilibrium equations directly. It solves their weak form. Here is the path for small strain linear elasticity.

Strong form

Find the displacement field \(\mathbf{u}\) on \(\Omega\) such that

\[ \nabla\cdot\boldsymbol{\sigma} + \mathbf{b} = \mathbf{0} \ \text{in}\ \Omega,\qquad \mathbf{u} = \bar{\mathbf{u}} \ \text{on}\ \Gamma_D,\qquad \boldsymbol{\sigma}\,\mathbf{n} = \bar{\mathbf{t}} \ \text{on}\ \Gamma_N, \]

with \(\boldsymbol{\sigma} = \mathbb{C} : \boldsymbol{\varepsilon}\) and \(\boldsymbol{\varepsilon}(\mathbf{u}) = \tfrac{1}{2}\left(\nabla\mathbf{u} + \nabla\mathbf{u}^{\mathsf{T}}\right)\).

Weak form

Multiply by a test function \(\mathbf{v}\) that vanishes on \(\Gamma_D\), integrate over \(\Omega\), and apply the divergence theorem. Because \(\boldsymbol{\sigma}\) is symmetric, \(\nabla\mathbf{v} : \boldsymbol{\sigma} = \boldsymbol{\varepsilon}(\mathbf{v}) : \boldsymbol{\sigma}\), which gives

\[ \int_\Omega \boldsymbol{\varepsilon}(\mathbf{v}) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{u})\,\mathrm{d}\Omega = \int_\Omega \mathbf{v}\cdot\mathbf{b}\,\mathrm{d}\Omega + \int_{\Gamma_N} \mathbf{v}\cdot\bar{\mathbf{t}}\,\mathrm{d}\Gamma \quad \forall\, \mathbf{v} \in \mathcal{V}. \]

The Neumann condition is now natural: it enters through the load term instead of being imposed on the solution space. The Dirichlet condition stays essential.

Discretisation

With shape functions collected in \(\mathbf{N}\) and their derivatives in the strain displacement matrix \(\mathbf{B}\), the Galerkin method gives \(\mathbf{K}\mathbf{u} = \mathbf{f}\) with element stiffness

\[ \mathbf{K}_e = \int_{\Omega_e} \mathbf{B}^{\mathsf{T}} \mathbf{D}\, \mathbf{B}\,\mathrm{d}\Omega . \]

A 1D check

A bar of length \(L\), clamped at \(x = 0\), free at \(x = L\), under a uniform axial line load \(q\). The exact solution of \(EA\,u'' = -q\) is

\[ u(x) = \frac{q}{EA}\left(Lx - \frac{x^2}{2}\right), \qquad u(L) = \frac{qL^2}{2EA}. \]

import numpy as np

E, A, L, q = 210e9, 1e-4, 1.0, 1e3   # [Pa], [m^2], [m], [N/m]
n = 8                                 # linear elements
h = L / n

K = np.zeros((n + 1, n + 1))
f = np.zeros(n + 1)
ke = E * A / h * np.array([[1.0, -1.0], [-1.0, 1.0]])
fe = q * h / 2 * np.array([1.0, 1.0])  # consistent load vector
for e in range(n):
    idx = [e, e + 1]
    K[np.ix_(idx, idx)] += ke
    f[idx] += fe

u = np.zeros(n + 1)
u[1:] = np.linalg.solve(K[1:, 1:], f[1:])   # u(0) = 0 imposed by elimination

print(u[-1], q * L**2 / (2 * E * A))        # both 2.381e-05 m

The nodal values match the exact solution for any \(n\). This is a known property of 1D problems with constant \(EA\) and a consistent load vector: the Green's function lies in the finite element space, so nodal values are exact. Between the nodes the linear interpolant still differs from the quadratic exact solution, and the element strains are constant approximations of a linear field.

References

T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover, 2000. Chapters 1 and 2.
O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, 7th ed., Butterworth Heinemann, 2013.