Question 1 of 3: Reactions by virtual work, and why the element must contain the rigid body modes (Problem 1)

Question 1: Reactions by virtual work, and why the element must contain the rigid body modes (Problem 1)

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.

Given. A prismatic beam of depth $h$ rests on two simple supports a distance $L = a + b$ apart and carries a single downward concentrated force $P$ at a distance $a$ from the left support (support 1) and $b$ from the right support (support 2). Usual beam theory applies, so the material is linear elastic and the only strain measure that does work is the axial strain associated with bending. To make the arithmetic concrete the worked numbers below use $P = 60\ \text{kN}$, $a = 4\ \text{m}$ and $b = 2\ \text{m}$, so that $L = 6\ \text{m}$; the symbolic result is what the question asks for.

Find. (1.1) The two support reactions $R_1$ and $R_2$, obtained from the principle of virtual work rather than from statics. (1.2) A demonstration that the same reaction can be recovered from a finite element model only if the element displacement functions are able to reproduce rigid body motion exactly.

[Figure not reproduced: Figure 1(a) redrawn: simply supported beam of depth h, span L = a + b, carrying the concentrated force P at distance a from support 1. See the official exam paper.]

Approach. Release one reaction at a time and give the beam the rigid body virtual displacement that the release permits; a rigid motion produces no virtual strain, so the internal virtual work vanishes identically and the external virtual work equation alone delivers the reaction.

1. Part 1.1 — state the principle of virtual work. For a body in equilibrium, any kinematically admissible virtual displacement field satisfies $$\delta W_{\text{ext}} = \delta U_{\text{int}}, \qquad \delta U_{\text{int}} = \int_V \sigma_{ij}\, \delta \varepsilon_{ij}\, dV .$$ Here $\sigma_{ij}$ is the real stress field carried by the beam and $\delta\varepsilon_{ij}$ is the strain derived from the virtual displacement. The two fields are independent: that independence is what makes the method useful.
2. Choose a virtual field that releases $R_1$. Remove the left support, keep the right one, and let the beam rotate as a rigid bar about support 2. Taking the upward virtual displacement at support 1 as $\delta_1$, the field is $$\delta v(x) = \delta_1\left(1 - \frac{x}{L}\right), \qquad 0 \le x \le L ,$$ measured from support 1. It is kinematically admissible for the released structure because it respects the one remaining support, $\delta v(L) = 0$.
3. Show the internal virtual work is zero. The virtual field is linear in $x$, so its curvature is $$\delta\kappa = -\frac{d^{2}(\delta v)}{dx^{2}} = 0 ,$$ and therefore $\delta\varepsilon_{x} = -y\,\delta\kappa = 0$ at every point of the cross-section. A rigid body motion stretches nothing, so $\delta U_{\text{int}} = \int_V \sigma_x \delta\varepsilon_x\, dV = 0$ whatever the real stress $\sigma_x$ happens to be. This is the whole trick of the method.
4. Write the external virtual work and solve for $R_1$. Only two forces move: the released reaction, through $\delta_1$, and the applied load, through the value of the virtual field at the load point, $\delta v(a) = \delta_1 b / L$. Because the load acts downward while the virtual displacement is upward, its contribution is negative: $$\delta W_{\text{ext}} = R_1 \delta_1 - P\,\frac{\delta_1 b}{L} = 0 .$$ Dividing through by the arbitrary amplitude $\delta_1$ gives $$\boxed{\,R_1 = \frac{P\,b}{L} = \frac{P\,b}{a+b}\,}$$ which for the illustrative numbers is $R_1 = 60 \times 2 / 6 = 20\ \text{kN}$.
5. Repeat with the other release for $R_2$. Now remove the right support and rotate the bar about support 1, so that $\delta v(x) = \delta_2\, x / L$ and $\delta v(a) = \delta_2 a / L$. The same argument gives $R_2 \delta_2 - P \delta_2 a / L = 0$, hence $$\boxed{\,R_2 = \frac{P\,a}{L} = \frac{P\,a}{a+b}\,}$$ that is $R_2 = 60 \times 4 / 6 = 40\ \text{kN}$.
6. Check. Adding the two results, $R_1 + R_2 = P(b + a)/L = P$, and taking moments about support 1 returns $R_2 L = P a$. Virtual work with rigid body fields has simply reproduced the equations of global equilibrium, which is exactly what it must do.

Two features of that derivation are worth naming before Part 1.2, because they are the hinges of the argument. First, the answer never used the beam's stiffness, its depth $h$ or its material: the reaction of a statically determinate structure is a matter of equilibrium alone. Second, the reason the stiffness dropped out is that the chosen virtual field was rigid, so the internal work integral collapsed. Any virtual field that bent the beam would have brought the unknown real stress distribution back into the equation and would have left $R_1$ undetermined.

Part 1.2 — why the element must contain the rigid body modes

Model the same beam with the four-node element of Figure 1(b). The displacement inside the element is interpolated from the nodal displacement vector $\mathbf{d}$ through the shape function matrix, $\mathbf{u} = \mathbf{N}\,\mathbf{d}$, the strains follow as $\boldsymbol{\varepsilon} = \mathbf{B}\,\mathbf{d}$ with $\mathbf{B} = \partial \mathbf{N}$, and the discrete equilibrium equations obtained by applying virtual work to the assembled model are $$\mathbf{K}\,\mathbf{D} = \mathbf{R}, \qquad \mathbf{K} = \int_V \mathbf{B}^{\mathsf T}\, \mathbf{C}\, \mathbf{B}\; dV .$$ The reaction at support 1 is recovered afterwards from the row of this system that belongs to the supported degree of freedom.

Now repeat the Part 1.1 argument inside the discrete model. Let $\mathbf{D}_{\text{rb}}$ be the nodal values of the rigid rotation $\delta_1(1 - x/L)$ about support 2, and use it as the virtual displacement vector. Pre-multiplying the discrete equilibrium equations by it gives $$\mathbf{D}_{\text{rb}}^{\mathsf T}\, \mathbf{K}\, \mathbf{D} = \mathbf{D}_{\text{rb}}^{\mathsf T}\, \mathbf{R} .$$ The right-hand side is precisely the external virtual work of Step 4, namely $R_1 \delta_1 - P \delta_1 b / L$. So the finite element model reproduces the exact reaction if and only if the left-hand side vanishes.

The left-hand side vanishes exactly when the interpolation can represent the rigid body mode. Suppose it can, so that some nodal vector $\mathbf{D}_{\text{rb}}$ generates the rigid field pointwise, $\mathbf{N}\,\mathbf{D}_{\text{rb}} = \delta_1(1 - x/L)$ throughout the element. Differentiating, the strains it produces are identically zero, $$\mathbf{B}\,\mathbf{D}_{\text{rb}} = \mathbf{0} \;\Longrightarrow\; \mathbf{K}\,\mathbf{D}_{\text{rb}} = \int_V \mathbf{B}^{\mathsf T} \mathbf{C}\, \left(\mathbf{B}\,\mathbf{D}_{\text{rb}}\right) dV = \mathbf{0} ,$$ and because $\mathbf{K}$ is symmetric, $\mathbf{D}_{\text{rb}}^{\mathsf T} \mathbf{K} \mathbf{D} = \left(\mathbf{K}\,\mathbf{D}_{\text{rb}}\right)^{\mathsf T} \mathbf{D} = 0$. What survives is $\mathbf{D}_{\text{rb}}^{\mathsf T} \mathbf{R} = 0$, which written out is $R_1 \delta_1 - P \delta_1 b / L = 0$ and therefore $R_1 = Pb/L$: the finite element model returns the Part 1.1 answer.

Suppose instead the interpolation cannot represent the rigid mode. Then $\mathbf{B}\,\mathbf{D}_{\text{rb}} \ne \mathbf{0}$ for the nodal vector closest to a rigid motion, the element stores strain energy while merely being carried through space, $\mathbf{K}\,\mathbf{D}_{\text{rb}} \ne \mathbf{0}$, and the identity above fails. The recovered reaction becomes $$R_1 = \frac{P\,b}{L} + \frac{1}{\delta_1}\, \mathbf{D}_{\text{rb}}^{\mathsf T}\, \mathbf{K}\, \mathbf{D} ,$$ that is, the exact value plus a spurious term that has nothing to do with the applied load. The model then violates global equilibrium: the computed reactions do not sum to $P$. This is not a discretisation error that a finer mesh will remove, because the defect is present in every element at every size, which is why the ability to represent rigid body modes is one of the two halves of the completeness requirement (the other being the ability to represent constant strain states, which governs convergence rather than equilibrium).

An equivalent and very quick statement of the same condition is the row-sum property. Give the whole model a rigid translation, $\mathbf{D} = \mathbf{1}$ in one direction. No force is needed to carry a free body through space, so $\mathbf{K}\,\mathbf{1} = \mathbf{0}$, which says that every row of $\mathbf{K}$ sums to zero. Since reactions are read from exactly those rows, a stiffness matrix whose rows do not sum to zero cannot give equilibrating reactions. It is worth checking numerically on any new element: assemble $\mathbf{K}$ before applying supports, sum each row, and confirm the result is zero to machine precision.

The four-node element in Figure 1(b) does satisfy the requirement. Its bilinear shape functions form a partition of unity, $\sum_i N_i = 1$, so the nodal set $u_i = c_1,\ v_i = c_2$ reproduces a rigid translation exactly; and because the interpolation contains the complete linear polynomial, the small rigid rotation $u = -c_3 y$, $v = c_3 x$ is also reproduced exactly. The element therefore has the three rigid body modes a plane element must have, and the virtual work argument above goes through unchanged.