Question 1 of 7: Basis vs. shape functions, Galerkin–Ritz equivalence, adaptivity, and Timoshenko beam convergence [20 marks]

Question 1: Basis vs. shape functions, Galerkin–Ritz equivalence, adaptivity, and Timoshenko beam convergence [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 [5 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 functions $\phi_i$ are the basis and the $c_i$ are generalized coordinates: they are pure algebraic multipliers with no physical identity. The monomials $1,\,x,\,x^2,\dots$ and the trigonometric set $\sin(n\pi x/L)$ are classic bases.

A shape function (interpolation function) is the particular basis obtained after the generalized coordinates have been eliminated in favour of the nodal values of the field variable, 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$. That elimination forces two properties that an arbitrary basis need not have: the Kronecker-delta (cardinal) property $N_i(\mathbf{x}_j)=\delta_{ij}$, and partition of unity $\sum_i N_i = 1$ over the element, which is what guarantees that a rigid-body (constant) field is reproduced exactly.

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 the change from $\phi_i$ to $N_i$ is a change of basis, not a change of the approximation itself.

(b) When Bubnov–Galerkin coincides with Ritz [4 marks]

The two discretizations produce the identical algebraic system when the governing differential operator $\mathcal{L}$ is linear, self-adjoint and positive definite, so that a quadratic functional (a variational principle, e.g. total potential energy) exists whose Euler–Lagrange equation is the governing equation, and the same approximation functions — satisfying the same essential boundary conditions, with the natural boundary conditions carried in the weak form — are used in both methods.

The reason is structural: in the Bubnov–Galerkin method the weight functions are taken from the trial space itself ($w_j=N_j$), so the weighted-residual statement $\int_\Omega N_j\,\mathcal{R}\,d\Omega = 0$ is exactly the stationarity condition $\partial \Pi/\partial u_j = 0$ of the Ritz functional. If the operator is not self-adjoint — the convection term in an advection–diffusion problem is the standard counterexample — no such functional exists, Ritz is simply unavailable, and Galerkin (or a Petrov–Galerkin variant with different weights) must be used on its own.

(c) h-version and p-version adaptivity [4 marks]

In the h-version the polynomial order of the element interpolation is held fixed and accuracy is bought by reducing the element size $h$ — refining the mesh, either uniformly or, in an adaptive scheme, only where an a-posteriori error indicator is large. Convergence is algebraic, with the discretization error behaving as $\|e\| = O(h^{p})$ in the energy norm for a smooth solution.

In the p-version the mesh is held fixed and accuracy is bought by raising the polynomial degree $p$ of the element interpolation, normally through hierarchic shape functions so that the previously assembled equations are retained and only new rows and columns are appended. For a smooth solution the p-version converges exponentially in $p$, which is far faster than h-refinement, but it degenerates to algebraic convergence near a re-entrant corner or other singularity. The combined hp-version — h-refinement toward singularities, p-enrichment in smooth regions — recovers exponential convergence for the whole problem and is the basis of most modern adaptive codes.

(d) The Timoshenko beam convergence problem [7 marks]

The phenomenon is shear locking (also called parasitic shear). In a Timoshenko beam the two independent fields are the transverse displacement $w$ and the cross-section rotation $\theta$, and the transverse shear strain is the difference of a derivative of one field and the other field itself:

As the beam becomes slender the ratio $k_sGAL^{2}/EI$ grows without bound, so the shear term behaves as a numerical penalty that drives the discrete $\gamma_{xz}$ to zero at every integration point. If the polynomial space of $dw/dx$ does not coincide with the polynomial space of $\theta$, that penalty cannot be satisfied by the shear mode alone — it must also annihilate part of the bending mode. The element then carries spurious shear energy in what should be pure bending, is far too stiff, and the computed deflection creeps toward the correct answer only as the mesh is refined enormously: exactly the “very slow convergence” being reported.

The requirement 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, and they are ranked here 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. Sampling $\gamma_{xz}$ only where it is field-consistent removes the spurious constraints; a hourglass-control or stabilization term is then added to suppress the zero-energy mode that reduced integration can introduce.
2. Consistent (linked) interpolation. Use the consistent-interpolation element — $w$ quadratic, $\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 / assumed-strain formulations. Interpolate $\gamma_{xz}$ (or the shear force) independently 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 plate and beam elements used in commercial software.
4. p-refinement. Raising the element order dilutes the locking constraint ratio, so high-order elements lock much less; this mitigates rather than cures, and is the fallback when the element library cannot be changed.