Question 3 of 3: The bilinear Q4 element and a reinforced plate

Question 3: The bilinear Q4 element and a reinforced plate (50 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.

Given.

Figure 3(a) fixes the local numbering: node 1 at the origin, node 2 at $(L,0)$, node 3 at $(L,L)$ and node 4 at $(0,L)$, counter-clockwise. Figure 3(b) is a square plate $2L$ wide and $2L$ high, divided into four $L\times L$ elements. Its top and bottom edges run on rollers (vertical movement prevented, horizontal movement free) and each of those two edges is reinforced by a bar of area $A$ — these are the "truss bar elements" labelled on the figure. Edge loads $P$, $2P$, $P$ act outward on the left edge and outward on the right edge at the top, mid-height and bottom nodes respectively, so the plate is stretched horizontally by a self-equilibrated system totalling $4P$ each side. Nodes 1, 2, 3 and 4 named in 3.4 are the four corners of the top-right element: node 1 the centre of the plate, node 2 the mid-height point of the right edge, node 3 the top-right corner, node 4 the mid-point of the top edge. Point C lies inside that element at local coordinates $(L/4, L/4)$ from node 1.

Find. (3.1) the four bilinear shape functions; (3.2) the 3×8 strain-displacement matrix; (3.3) the three stiffness coefficients of the $\{F_2,F_3\}$–$\{u_2,u_3\}$ relation; (3.4) $u_2$, $u_3$ and the strain components at C.

[Figure not reproduced: Fig. 3(a) redrawn: the square bilinear Q4 element, side L, with local axes originating at node 1 and the nodes numbered counter-clockwise. See the official exam paper.]

Approach. Build the shape functions as products of one-dimensional linear Lagrange functions, differentiate them to get $[B]$, integrate $[k]=t\int\int [B]^{T}[D][B]\,dx\,dy$ for the two entries the question names, then use the double symmetry of Figure 3(b) to reduce the nine-node model to the same two unknowns.

1. Part 3.1 — construct the shape functions as Lagrange products. The Q4 displacement field is bilinear, $u(x,y) = a_1 + a_2x + a_3y + a_4xy$, four constants matched to four nodal values, and the same for $v$. Rather than inverting a 4×4, note that the element is a rectangle whose sides are parallel to the axes, so the interpolation separates: in $x$ the one-dimensional linear Lagrange functions are $\left(1-\frac{x}{L}\right)$ and $\frac{x}{L}$, and in $y$ they are $\left(1-\frac{y}{L}\right)$ and $\frac{y}{L}$. The shape function of a node is the product of the two that equal one at that node: $$\begin{aligned} N_1 &= \left(1-\tfrac{x}{L}\right)\left(1-\tfrac{y}{L}\right) &&\text{node 1 at } (0,0), \\ N_2 &= \tfrac{x}{L}\left(1-\tfrac{y}{L}\right) &&\text{node 2 at } (L,0), \\ N_3 &= \tfrac{x}{L}\cdot\tfrac{y}{L} = \tfrac{xy}{L^{2}} &&\text{node 3 at } (L,L), \\ N_4 &= \left(1-\tfrac{x}{L}\right)\tfrac{y}{L} &&\text{node 4 at } (0,L), \end{aligned}$$ which is exactly the set the question asks us to show.
2. Verify the three properties that make them admissible. First the Kronecker-delta property: each $N_i$ takes the value 1 at its own node and 0 at the other three, because at any other node at least one of its two factors vanishes. Second, partition of unity, $$\sum_{i=1}^{4}N_i = \left[\left(1-\tfrac{x}{L}\right)+\tfrac{x}{L}\right] \left[\left(1-\tfrac{y}{L}\right)+\tfrac{y}{L}\right] = 1 ,$$ so a rigid-body translation is reproduced exactly. Third, linear completeness: $\sum N_i x_i = x$ and $\sum N_i y_i = y$, so a constant strain state is reproduced exactly and the element passes the patch test. Along any edge one coordinate is fixed and the field degenerates to a linear function of the other, so adjacent elements sharing an edge stay in contact — the element is $C^{0}$ conforming.
3. Part 3.2 — differentiate to obtain $[B]$. The small-strain definitions are $\varepsilon_x = \partial u/\partial x$, $\varepsilon_y = \partial v/\partial y$ and $\gamma_{xy} = \partial u/\partial y + \partial v/\partial x$. Differentiating each shape function, $$\begin{aligned} \frac{\partial N_1}{\partial x} &= -\frac{1}{L}\left(1-\frac{y}{L}\right), & \frac{\partial N_1}{\partial y} &= -\frac{1}{L}\left(1-\frac{x}{L}\right), \\ \frac{\partial N_2}{\partial x} &= +\frac{1}{L}\left(1-\frac{y}{L}\right), & \frac{\partial N_2}{\partial y} &= -\frac{x}{L^{2}}, \\ \frac{\partial N_3}{\partial x} &= +\frac{y}{L^{2}}, & \frac{\partial N_3}{\partial y} &= +\frac{x}{L^{2}}, \\ \frac{\partial N_4}{\partial x} &= -\frac{y}{L^{2}}, & \frac{\partial N_4}{\partial y} &= +\frac{1}{L}\left(1-\frac{x}{L}\right). \end{aligned}$$ Because $\{d\}$ interleaves the two components node by node, $[B]$ is assembled from four 3×2 blocks, one per node, in the order $[B]=\left[\,[B_1]\ [B_2]\ [B_3]\ [B_4]\,\right]$.
4. Write the block and read off what it means. Each block is $$[B_i] = \begin{bmatrix} \dfrac{\partial N_i}{\partial x} & 0 \\[6pt] 0 & \dfrac{\partial N_i}{\partial y} \\[6pt] \dfrac{\partial N_i}{\partial y} & \dfrac{\partial N_i}{\partial x} \end{bmatrix},$$ so, for example, $$[B_2] = \begin{bmatrix} \dfrac{1}{L}\left(1-\dfrac{y}{L}\right) & 0 \\[6pt] 0 & -\dfrac{x}{L^{2}} \\[6pt] -\dfrac{x}{L^{2}} & \dfrac{1}{L}\left(1-\dfrac{y}{L}\right)\end{bmatrix}, \qquad [B_3] = \begin{bmatrix} \dfrac{y}{L^{2}} & 0 \\[6pt] 0 & \dfrac{x}{L^{2}} \\[6pt] \dfrac{x}{L^{2}} & \dfrac{y}{L^{2}}\end{bmatrix}.$$ Two features matter later. The entries are linear in $x$ and $y$, not constant, so unlike the constant-strain triangle the Q4 carries a varying strain field — $\varepsilon_x$ varies with $y$ only and $\varepsilon_y$ with $x$ only, while $\gamma_{xy}$ varies with both. And every row of $[B]$ sums to zero across the eight columns, which is the statement that a rigid-body motion produces no strain.
5. Part 3.3 — integrate the two coefficients the question names. For plane stress $$[D] = \frac{E}{1-\nu^{2}}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \dfrac{1-\nu}{2}\end{bmatrix}, \qquad [k] = t\int_{0}^{L}\!\!\int_{0}^{L}[B]^{T}[D][B]\,dx\,dy .$$ The entry coupling $u_2$ with itself picks up only the $D_{11}$ and $D_{33}$ terms, because the $u_2$ column of $[B]$ has zero in the $\varepsilon_y$ row: $$k_{11}=k[u_2,u_2] = t\!\int\!\!\int\!\left[D_{11}\left(\frac{\partial N_2}{\partial x}\right)^{2} + D_{33}\left(\frac{\partial N_2}{\partial y}\right)^{2}\right]dA .$$ The two integrals are elementary: $\int\!\!\int (\partial N_2/\partial x)^{2}dA = \frac{1}{L^{2}}\cdot L\cdot\frac{L}{3} = \frac13$ and $\int\!\!\int (\partial N_2/\partial y)^{2}dA = \frac{1}{L^{4}}\cdot\frac{L^{3}}{3}\cdot L = \frac13$, both independent of $L$. Hence $$k_{11} = \frac{Et}{1-\nu^{2}}\left[\frac13 + \frac{1-\nu}{2}\cdot\frac13\right] = \frac{Et}{1-\nu^{2}}\left[\frac13 + \frac{1-\nu}{6}\right],$$ exactly as printed. The same two integrals with $N_3$ in place of $N_2$ give $\frac13$ and $\frac13$ again, so $k_{22}=k[u_3,u_3]=k_{11}$, which is the symmetry the question asserts.
6. Integrate the cross term, and note the sign. The coupling entry is $$k_{12}=k[u_2,u_3] = t\!\int\!\!\int\!\left[D_{11}\frac{\partial N_2}{\partial x} \frac{\partial N_3}{\partial x} + D_{33}\frac{\partial N_2}{\partial y} \frac{\partial N_3}{\partial y}\right]dA .$$ The direct-stress part gives $\int\!\!\int \frac{1}{L}\left(1-\frac{y}{L}\right)\frac{y}{L^{2}}\,dA = \frac{1}{L^{2}}\left(\frac{L^{2}}{2}-\frac{L^{2}}{3}\right) = \frac16$, but the shear part carries a negative integrand, because $\partial N_2/\partial y = -x/L^{2}$ while $\partial N_3/\partial y = +x/L^{2}$: $$\int\!\!\int \left(-\frac{x}{L^{2}}\right)\left(\frac{x}{L^{2}}\right)dA = -\frac{1}{L^{4}}\cdot\frac{L^{3}}{3}\cdot L = -\frac13 .$$ Therefore $$\boxed{\ k_{12} = \frac{Et}{1-\nu^{2}}\left[\frac16 - \frac{1-\nu}{6}\right] = \frac{Et\,\nu}{6\left(1-\nu^{2}\right)}\ }$$ Assembling the three coefficients into the 2×2 relation $\{F_2,F_3\}^{T} = [k_{11}, k_{12}; k_{12}, k_{22}]\{u_2,u_3\}^{T}$ completes 3.3 — the relation holds whenever the element's other six degrees of freedom are held at zero, which is exactly the condition symmetry will create in 3.4.

[Figure not reproduced: Fig. 3(b) redrawn: the 2L × 2L plate meshed with four L × L Q4 elements. Rollers on the top and bottom edges prevent vertical movement; the heavy lines on those two edges are the reinforcing truss bars of area A. The labelled nodes 1, 2, 3, 4 are the corners of the top-right element, and. See the official exam paper.]

7. Part 3.4 — use the double symmetry to reduce the model. Place the origin at the centre node (node 1). The structure — plate, bars and rollers — is symmetric about both centrelines. The loading is symmetric about the horizontal centreline (P, 2P, P repeat above and below) and antisymmetric about the vertical one (outward on both sides). Antisymmetry of the horizontal load forces $u(-x,y) = -u(x,y)$, hence $u = 0$ along the whole centre column. Symmetry of the load about $y=0$ forces $v(x,-y) = -v(x,y)$, hence $v = 0$ along the mid-row; and the rollers already set $v = 0$ on the top and bottom edges. Every vertical displacement in the nine-node mesh therefore vanishes, and the only unknowns left are the horizontal movements of the right edge, $u_2$ at mid-height and $u_3$ at the top (equal by symmetry to the one at the bottom).
8. Assemble the reduced system. Inside the top-right element every degree of freedom except $u_2$ and $u_3$ is now zero, so its contribution is precisely the 2×2 matrix of 3.3. Node 3 belongs to one element only and carries the applied $P$, plus the top reinforcing bar, which spans from node 4 to node 3 with $u_4 = 0$ and so adds $\frac{EA}{L}u_3$. Node 2 is shared by the top-right and bottom-right elements and carries $2P$; the two elements contribute identically, so that equation may be divided through by two. The result is $$\begin{bmatrix} k_{11} & k_{12} \\ k_{12} & k_{11}+\dfrac{EA}{L}\end{bmatrix} \begin{Bmatrix} u_2 \\ u_3\end{Bmatrix} = \begin{Bmatrix} P \\ P \end{Bmatrix} .$$ Only the node-3 equation feels the reinforcement, because the bars lie on the top and bottom edges and node 2 sits between them.
9. Evaluate the coefficients. With $Et = (200\,000)(5) = 1.000\times10^{6}$ N/mm and $\nu = 0.25$, so that $1-\nu^{2}=0.9375$, $$\begin{aligned} k_{11}=k_{22} &= \frac{1.000\times10^{6}}{0.9375}\left[\frac13+\frac{0.75}{6}\right] = 4.8889\times10^{5}\ \text{N/mm}, \\ k_{12} &= \frac{(1.000\times10^{6})(0.25)}{6(0.9375)} = 4.4444\times10^{4}\ \text{N/mm}, \\ \frac{EA}{L} &= \frac{Et}{2} = 5.000\times10^{5}\ \text{N/mm}. \end{aligned}$$ The reinforcement is worth slightly more than the plate itself at that node, which is why it changes the answer so much.
10. Solve for the two displacements. Substituting and taking $P = 5000$ N, $$\begin{bmatrix} 4.8889 & 0.4444 \\ 0.4444 & 9.8889 \end{bmatrix} \times 10^{5}\begin{Bmatrix} u_2 \\ u_3\end{Bmatrix} = \begin{Bmatrix} 5000 \\ 5000 \end{Bmatrix}\ \text{N},$$ whose determinant is $4.8146\times10^{11}$, giving $$\boxed{\ u_2 = 9.808\times10^{-3}\ \text{mm}, \qquad u_3 = 4.615\times10^{-3}\ \text{mm}\ }$$ Both are outward, as the loading demands, and node 2 moves 2.1 times as far as node 3 even though it carries only twice the force spread over twice as much element — the edge bar is what holds node 3 back. Of the 5.00 kN applied at node 3, the bar takes $\frac{EA}{L}u_3 = 2.31$ kN and the plate the remaining 2.69 kN.
11. Compute the strains at C. Point C lies inside the top-right element at $x = y = L/4 = 5$ mm from node 1. With $v \equiv 0$ and $u_1 = u_4 = 0$, only the $u_2$ and $u_3$ columns of $[B]$ survive: $$\begin{aligned} \varepsilon_x &= \frac{\partial N_2}{\partial x}u_2 + \frac{\partial N_3}{\partial x}u_3 = \frac{1}{L}\left(1-\frac{y_c}{L}\right)u_2 + \frac{y_c}{L^{2}}u_3 = \frac{0.75\,u_2 + 0.25\,u_3}{L}, \\ \varepsilon_y &= \frac{\partial N_2}{\partial y}v_2 + \frac{\partial N_3}{\partial y}v_3 = 0 \quad\text{(all vertical displacements vanish)}, \\ \gamma_{xy} &= \frac{\partial N_2}{\partial y}u_2 + \frac{\partial N_3}{\partial y}u_3 = \frac{x_c}{L^{2}}\left(u_3-u_2\right) = \frac{0.25\left(u_3-u_2\right)}{L} . \end{aligned}$$ Substituting $u_2$, $u_3$ and $L = 20$ mm, $$\boxed{\ \varepsilon_x = 4.255\times10^{-4}, \qquad \varepsilon_y = 0, \qquad \gamma_{xy} = -6.490\times10^{-5}\ }$$
12. Interpret, and convert to stress as a cross-check. That $\varepsilon_y$ is exactly zero everywhere is a consequence of the rollers, not an accident of the point chosen: the plate is stretched in $x$ while being prevented from contracting in $y$, so it is in a state of confined extension. The shear strain is negative and non-zero even though the loading is purely horizontal, because $u$ varies with $y$ across the element — the top of the right edge lags the middle. Applying the plane-stress law, $$\begin{aligned} \sigma_x &= \frac{E}{1-\nu^{2}}\left(\varepsilon_x+\nu\varepsilon_y\right) = 90.8\ \text{MPa}, \ \sigma_y &= \frac{E\nu}{1-\nu^{2}}\,\varepsilon_x = 22.7\ \text{MPa}, \ \tau_{xy} &= \frac{E}{2(1+\nu)}\,\gamma_{xy} = -5.19\ \text{MPa}. \end{aligned}$$ The transverse stress $\sigma_y$ is exactly $\nu\sigma_x$, the reaction the rollers must supply, and it is a real design quantity: a plate detailed for the 90.8 MPa longitudinal stress alone would be missing a quarter of that again across the grain of the load.