Question 3 of 3: Shape functions, a stiffness coefficient and consistent nodal loads for the four-node rectangle (Problem 3)

Question 3: Shape functions, a stiffness coefficient and consistent nodal loads for the four-node rectangle (Problem 3)

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 four-node rectangular plane-stress element of width $a$ along $x$, height $b$ along $y$ and uniform thickness $t$, with the origin of the coordinate system at node 1 and the nodes numbered counter-clockwise: node 1 at $(0,0)$, node 2 at $(a,0)$, node 3 at $(a,b)$ and node 4 at $(0,b)$. The material is linear elastic and isotropic with modulus $E$ and Poisson's ratio $\nu$, and the plane-stress elasticity matrix $[C]$ is as printed on the paper. Each node carries two degrees of freedom, and the element displacement vector is ordered $\mathbf{d} = \{u_1,\ v_1,\ u_2,\ v_2,\ u_3,\ v_3,\ u_4,\ v_4\}$, so that degree of freedom 3 is the horizontal displacement $u_2$ of node 2. In part 3.3 a horizontal force $Q$ acts at the interior point $(3a/4,\ b/2)$.

Find. (3.1) The four shape functions of the element. (3.2) A derivation of the stiffness coefficient $k_{33}$ in the form printed on the paper. (3.3) The work-equivalent nodal forces produced by the interior point load $Q$.

[Figure not reproduced: Figure 3(a) redrawn: the four-node rectangle, origin at node 1, nodes numbered counter-clockwise. See the official exam paper.]

Approach. Build the interpolation from the complete bilinear polynomial and identify the shape functions as products of one-dimensional linear functions; differentiate them to obtain the strain-displacement matrix; integrate the appropriate entry of $\mathbf{B}^{\mathsf T}\mathbf{C}\,\mathbf{B}$ over the rectangle; and finally evaluate the shape functions at the load point to lump the interior force onto the nodes.

1. Part 3.1 — choose the displacement field. The element has four nodes and one displacement component per direction to interpolate, so exactly four constants are available. Take the bilinear polynomial $$u(x,y) = \alpha_1 + \alpha_2 x + \alpha_3 y + \alpha_4 xy ,$$ and the same form with constants $\alpha_5 \dots \alpha_8$ for $v(x,y)$. It contains the complete linear polynomial, which secures the rigid body and constant strain states, plus the $xy$ term needed to make the field vary linearly along each of the four edges.
2. Impose the nodal values and invert. Writing $u(0,0) = u_1$, $u(a,0) = u_2$, $u(a,b) = u_3$, $u(0,b) = u_4$ and solving the four equations for $\alpha_1 \dots \alpha_4$ gives $\alpha_1 = u_1$, $\alpha_2 = (u_2-u_1)/a$, $\alpha_3 = (u_4-u_1)/b$ and $\alpha_4 = (u_1-u_2+u_3-u_4)/(ab)$. Collecting the coefficient of each nodal value produces the interpolation $u = \sum_i N_i u_i$ with
3. Read off the shape functions. The result is the product of two one-dimensional linear Lagrange functions, one in each direction: $$\boxed{\; N_1 = \left(1-\frac{x}{a}\right)\left(1-\frac{y}{b}\right), \quad N_2 = \frac{x}{a}\left(1-\frac{y}{b}\right), \quad N_3 = \frac{x}{a}\,\frac{y}{b}, \quad N_4 = \left(1-\frac{x}{a}\right)\frac{y}{b}\;}$$ Each is unity at its own node and zero at the other three, each varies linearly along every edge, and their sum is $\left[(1-x/a) + x/a\right]\left[(1-y/b) + y/b\right] = 1$ everywhere. The partition of unity is what guarantees the rigid body translation discussed in Problem 1, and the linear edge variation is what guarantees inter-element compatibility when two such rectangles are placed side by side.
4. Part 3.2 — form the strain-displacement matrix. With $\varepsilon_x = \partial u/\partial x$, $\varepsilon_y = \partial v/\partial y$ and $\gamma_{xy} = \partial u/\partial y + \partial v/\partial x$, the columns of $\mathbf{B}$ belonging to node $i$ are $$\mathbf{B}_i = \begin{bmatrix} \partial N_i/\partial x & 0 \\ 0 & \partial N_i/\partial y \\ \partial N_i/\partial y & \partial N_i/\partial x \end{bmatrix}, \qquad \mathbf{k} = t\int_{0}^{b}\!\!\int_{0}^{a} \mathbf{B}^{\mathsf T}\,\mathbf{C}\,\mathbf{B}\; dx\, dy .$$ The thickness is constant, so it comes outside the integral.
5. Isolate the entry that is wanted. Degree of freedom 3 is $u_2$, so $k_{33}$ involves only the column of $\mathbf{B}$ for that degree of freedom, whose two non-zero entries are $\partial N_2/\partial x$ in row 1 and $\partial N_2/\partial y$ in row 3. Because $C_{13} = C_{31} = 0$ for an isotropic material, the two contributions do not mix and $$k_{33} = t \int_{0}^{b}\!\!\int_{0}^{a} \left[ C_{11}\left(\frac{\partial N_2}{\partial x}\right)^{2} + C_{33}\left(\frac{\partial N_2}{\partial y}\right)^{2}\right] dx\, dy ,$$ with $C_{11} = E/(1-\nu^{2})$ the direct term and $C_{33} = E(1-\nu)/\left[2(1-\nu^{2})\right]$ the shear term. In words: pushing node 2 horizontally both stretches the element in $x$ and shears it in $y$, and $k_{33}$ is the sum of the two resistances.
6. Differentiate and integrate the direct term. From $N_2 = (x/a)(1 - y/b)$, $$\frac{\partial N_2}{\partial x} = \frac{1}{a}\left(1-\frac{y}{b}\right) \;\Longrightarrow\; \int_{0}^{b}\!\!\int_{0}^{a}\left(\frac{\partial N_2}{\partial x}\right)^{2} dx\,dy = \frac{1}{a^{2}}\,(a)\int_{0}^{b}\left(1-\frac{y}{b}\right)^{2} dy = \frac{1}{a}\cdot\frac{b}{3} = \frac{b}{3a} .$$
7. Differentiate and integrate the shear term. Similarly $$\frac{\partial N_2}{\partial y} = -\frac{1}{b}\,\frac{x}{a} \;\Longrightarrow\; \int_{0}^{b}\!\!\int_{0}^{a}\left(\frac{\partial N_2}{\partial y}\right)^{2} dx\,dy = \frac{1}{b^{2}}\,(b)\int_{0}^{a}\frac{x^{2}}{a^{2}}\, dx = \frac{1}{b}\cdot\frac{a}{3} = \frac{a}{3b} .$$
8. Assemble the two contributions. Substituting both integrals together with the two elasticity constants, $$k_{33} = t\left[\frac{E}{1-\nu^{2}}\cdot\frac{b}{3a} + \frac{E(1-\nu)}{2(1-\nu^{2})}\cdot\frac{a}{3b}\right] ,$$ and taking the common factor $Et/(1-\nu^{2})$ outside gives $$\boxed{\;k_{33} = \left(\frac{b}{3a} + \frac{(1-\nu)}{6}\,\frac{a}{b}\right)\frac{Et}{1-\nu^{2}}\;}$$ which is the expression printed on the paper. The same integrals with $N_1$ in place of $N_2$ give $k_{11} = k_{33}$, as symmetry of the rectangle requires: both are horizontal degrees of freedom at a corner.
9. Sanity-check the result numerically. For a steel element with $a = 200\ \text{mm}$, $b = 150\ \text{mm}$, $t = 10\ \text{mm}$, $E = 200\ \text{GPa}$ and $\nu = 0.30$ the two bracketed terms are $b/(3a) = 0.2500$ and $(1-\nu)a/(6b) = 0.15556$, while $Et/(1-\nu^{2}) = 2.1978 \times 10^{6}\ \text{N/mm}$, so $$k_{33} = 0.40556 \times 2.1978 \times 10^{6} = 8.912 \times 10^{5}\ \text{N/mm} = 891.2\ \text{kN/mm} .$$ Integrating $\mathbf{B}^{\mathsf T}\mathbf{C}\,\mathbf{B}$ numerically over the same rectangle reproduces this to machine precision, which is the check to run whenever a closed-form stiffness term is quoted from memory.
10. Part 3.3 — lump the interior point load onto the nodes. The work-equivalent nodal load vector for any body or point force is $\mathbf{f} = \int \mathbf{N}^{\mathsf T} \mathbf{b}\; dV$, and for a concentrated force at a single interior point $(x_0,y_0)$ the integral collapses to an evaluation: $$f_{xi} = Q\,N_i(x_0,\, y_0), \qquad f_{yi} = 0 ,$$ because the force has no vertical component. This is not an arbitrary lumping rule — it is the statement that the nodal forces do the same virtual work as $Q$ for every displacement field the element can represent.
11. Evaluate at the load point. With $x_0/a = 3/4$ and $y_0/b = 1/2$, $$N_1 = \tfrac{1}{4}\cdot\tfrac{1}{2} = \tfrac{1}{8}, \quad N_2 = \tfrac{3}{4}\cdot\tfrac{1}{2} = \tfrac{3}{8}, \quad N_3 = \tfrac{3}{4}\cdot\tfrac{1}{2} = \tfrac{3}{8}, \quad N_4 = \tfrac{1}{4}\cdot\tfrac{1}{2} = \tfrac{1}{8} ,$$ so the equivalent nodal forces, all horizontal and all in the direction of $Q$, are $$\boxed{\;f_{x1} = \frac{Q}{8}, \qquad f_{x2} = \frac{3Q}{8}, \qquad f_{x3} = \frac{3Q}{8}, \qquad f_{x4} = \frac{Q}{8}, \qquad f_{y1} = f_{y2} = f_{y3} = f_{y4} = 0 \;}$$
12. Check the lumping. The four forces sum to $Q(1/8 + 3/8 + 3/8 + 1/8) = Q$, as the partition of unity guarantees, and their resultant acts at $x = a\left(0 \cdot \tfrac{1}{8} + 1 \cdot \tfrac{3}{8} + 1 \cdot \tfrac{3}{8} + 0 \cdot \tfrac{1}{8}\right) = \tfrac{3a}{4}$, the correct line of action. The load divides three to one between the near and far edges in $x$, matching the position of the load, and equally between the bottom and top nodes because the load sits at mid-height. For $Q = 40\ \text{kN}$ the four values are 5, 15, 15 and 5 kN.

[Figure not reproduced: Figure 3(b) redrawn with the answer to 3.3: the horizontal force Q at (3a/4, b/2) and the four work-equivalent nodal forces, all horizontal. See the official exam paper.]

It is worth noticing what the lumping does not reproduce. The four nodal forces carry the correct resultant and the correct line of action in $x$, but they exert no couple about the $x$ axis, whereas the real load at mid-height exerts none either — here the two agree. Had the load been placed off mid-height, the equal split between bottom and top would have changed accordingly, and the consistent rule would still have been the shape function evaluation, not a hand-assigned share.