Question 4 of 7: Basis and shape functions, mesh convergence, Galerkin–Ritz equivalence, and Timoshenko beam locking

Question 4: Basis and shape functions, mesh convergence, Galerkin–Ritz equivalence, and Timoshenko beam locking (20 marks)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

(a) Basis function versus shape function [3 marks]

A basis function is any member of a set of linearly independent functions that spans the approximation space in which the unknown field is sought. Writing the trial solution as a generalized expansion $u(x)=\sum_i c_i\phi_i(x)$, the $\phi_i$ are the basis and the $c_i$ are generalized coordinates — pure algebraic multipliers with no physical identity. The monomials $1, x, x^{2},\ldots$ and the trigonometric set $\sin(n\pi x/L)$ are classical examples.

A shape function (interpolation function) is the particular basis obtained after those generalized coordinates have been eliminated in favour of the nodal values of the field, so that the approximation reads $u(x)=\sum_{i=1}^{n} N_i(x)\, u_i$ with $u_i$ the value of $u$ at node $i$. The elimination forces two properties an arbitrary basis need not possess: the cardinal (Kronecker-delta) property $N_i(\mathbf{x}_j)=\delta_{ij}$, and partition of unity $\sum_i N_i = 1$ over the element. In short, every shape function is a basis function of the element space, but not every basis function is a shape function — the shape functions are the nodal basis of that same space, and moving from $\phi_i$ to $N_i$ is a change of basis, not a change of the approximation.

(b) Why finite element solutions improve as the number of elements increases [3 marks]

Refining the mesh enlarges the trial space: with the polynomial order fixed, each additional element adds degrees of freedom, and the space of piecewise polynomials on the finer mesh contains the space on the coarser mesh. Since the Galerkin solution is the projection of the exact solution onto the trial space in the energy norm, enlarging a nested space cannot worsen that projection and generally improves it — formally $\|e\|_E = O(h^{p})$ in the energy norm for a smooth solution, where $h$ is the element size and $p$ the complete polynomial order.

The mechanical reading is equally simple: each element can only represent a polynomial variation of the field, so within one element the strain (or temperature gradient, or flux) is constrained to a low-order form. Shrinking the elements shrinks the interval over which that crude local representation must serve, so a piecewise-polynomial curve tracks the true curved solution ever more closely. Two caveats belong in the answer: convergence requires the elements to satisfy the completeness and compatibility (conformity) conditions, and it is monotone in the energy norm only for a properly formulated displacement-based element — a locking element converges so slowly that refinement alone is not a practical remedy, which is exactly the situation in part (d).

(c) Why Bubnov–Galerkin and Ritz give the identical discretization [4 marks]

They coincide because Galerkin takes its weight functions from the trial space itself. Writing the weighted-residual statement with $w_j=N_j$, $$\int_\Omega N_j\,\mathcal{R}(u_h)\, d\Omega = 0,\qquad j=1,\ldots,n$$ and comparing with the stationarity conditions $\partial\Pi/\partial u_j=0$ of the Ritz functional $\Pi$, the two sets of algebraic equations are term-for-term the same whenever a quadratic functional $\Pi$ exists whose Euler–Lagrange equation is the governing differential equation. That requires the differential operator $\mathcal{L}$ to be linear, self-adjoint and positive definite, and requires both methods to use the same approximation functions, satisfying the same essential boundary conditions, with the natural boundary conditions carried by the weak form.

For structural and thermal problems those conditions hold — the operator is self-adjoint and $\Pi$ is the total potential energy — so “Galerkin finite elements” and “Ritz finite elements” describe the same matrices, and the resulting stiffness matrix is symmetric. The equivalence fails as soon as self-adjointness fails: the convective term in an advection–diffusion problem admits no such functional, so Ritz is simply unavailable while Galerkin (or a Petrov–Galerkin variant with weights drawn from a different space) remains applicable and yields an unsymmetric system.

(d) Resolving the slow convergence of the Timoshenko beam [10 marks]

The phenomenon being observed is shear locking (parasitic shear). In a Timoshenko beam the transverse displacement $w$ and the cross-section rotation $\theta$ are independent fields, and the transverse shear strain is the difference between a derivative of one and the other field itself: $$\gamma_{xz}=\frac{dw}{dx}-\theta,\qquad \Pi = \frac12\int_0^L EI\left(\frac{d\theta}{dx}\right)^{2}dx + \frac12\int_0^L k_sGA\left(\frac{dw}{dx}-\theta\right)^{2}dx$$ As the beam becomes slender the dimensionless ratio $k_sGAL^{2}/EI$ grows without bound, so the shear term behaves as a numerical penalty that drives the discrete $\gamma_{xz}$ towards zero at every integration point. If the polynomial space of $dw/dx$ does not coincide with the space of $\theta$, that penalty cannot be satisfied by the shear mode alone and must also annihilate part of the bending mode. The element then stores spurious shear energy in what should be pure bending, becomes far too stiff, and the computed deflection creeps toward the correct value only under extreme refinement — the reported “very slow convergence”.

The requirement being violated is field consistency: $dw/dx$ and $\theta$ must be of the same polynomial order, so $w$ must be interpolated one order higher than $\theta$. Any of the following resolves the problem, ranked by how commonly they appear in production codes.

1. Reduced or selective integration of the shear term. Integrate the bending energy exactly but under-integrate $\int k_sGA\,\gamma_{xz}^{2}\, dx$ with one Gauss point fewer than exactness requires — one point for a linear element, two for a quadratic. Sampling $\gamma_{xz}$ only where it is field-consistent removes the spurious constraints; a stabilization or hourglass-control term is then added to suppress any zero-energy mode that under-integration introduces.
2. Consistent or linked interpolation. Use the consistent-interpolation element — $w$ quadratic and $\theta$ linear, with the interior $w$ degree of freedom condensed against the rotation field — or Reddy's interdependent-interpolation element with $w$ cubic and $\theta$ quadratic, which reproduces the Euler–Bernoulli solution exactly for constant $EI$.
3. Mixed or assumed-strain formulations. Interpolate $\gamma_{xz}$ (or the shear force) as an independent field in a Hellinger–Reissner two-field element, or tie the shear strain to sampling points as in the MITC family. These are the locking-free beam and plate elements used in commercial software.
4. $p$-refinement. Raising the element order dilutes the ratio of constraints to degrees of freedom, so high-order elements lock far less. This mitigates rather than cures, and is the fallback when the element library cannot be changed.