Q4 shape functions, plane-stress stiffness term, and a reinforced plate: 07-Str-B3, May 2017, Q3

Question 3: Q4 shape functions, plane-stress stiffness term, and a reinforced plate (equal value)

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.

3.1 Derivation of the bilinear shape functions

Given. The rectangular four-node element of Figure 3(a), of side $a$ in $x$ and $b$ in $y$, with node 1 at the origin, node 2 at $(a,0)$, node 3 at $(a,b)$ and node 4 at $(0,b)$, and two translational degrees of freedom per node.

Find. The four shape functions $N_i(x,y)$ that interpolate $u$ and $v$ from the nodal values.

Figure 3(a) — the four-node rectangular (Q4) element, nodes numbered counter-clockwise from the origin.

Approach. Assume the four-term bilinear polynomial that a four-node element can support, enforce the nodal values, and identify each shape function as the coefficient of its own nodal displacement; then verify the two properties (Kronecker delta and partition of unity) that any valid shape function must satisfy.

Choose the polynomial the element can support. Four nodes and one scalar field per node admit exactly four generalised coordinates, and the term set must be geometrically symmetric in $x$ and $y$ and complete to first order: $$u(x,y)=\alpha_1+\alpha_2x+\alpha_3y+\alpha_4xy$$ This is bilinear: linear along any line $x=\text{const}$ or $y=\text{const}$, and therefore linear along every edge, which is what guarantees inter-element compatibility.
Enforce the four nodal values. Substituting the four corner coordinates, $$u_1=\alpha_1,\quad u_2=\alpha_1+\alpha_2a,\quad u_3=\alpha_1+\alpha_2a+\alpha_3b+\alpha_4ab,\quad u_4=\alpha_1+\alpha_3b$$ Solving in turn: $\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$.
Substitute back and collect terms in each nodal value. Putting the four $\alpha$'s into the polynomial and grouping the coefficient of each $u_i$ gives $u=\sum N_i u_i$ with $$\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{xy}{ab},\quad N_4=\left(1-\frac{x}{a}\right)\frac{y}{b}\;}$$ which is the required result. The same four functions interpolate $v$, since the element uses one interpolation for both displacement components.
Recognise the product structure. Each $N_i$ is the product of a one-dimensional linear function in $x$ that is 1 at the node's own $x$ and 0 at the other, with the corresponding function in $y$. This is the tensor product of two two-node bar interpolations, and it is why the family generalises immediately to eight- and nine-node Lagrange elements.
Verify the Kronecker-delta property. Evaluating at each corner, $N_i(x_j,y_j)=\delta_{ij}$: for instance $N_2(a,0)=1\cdot1=1$ while $N_2(a,b)=1\cdot0=0$. This is what makes $u_i$ the physical displacement of node $i$ rather than an abstract coefficient.
Verify the partition of unity and first-order completeness. Expanding, $\sum N_i=1$ identically, so a rigid-body translation $u_i=c$ gives $u=c$ everywhere and generates no strain. Further, $\sum N_i x_i=x$ and $\sum N_i y_i=y$, so the element reproduces any linear displacement field exactly — it therefore passes the patch test and converges as the mesh is refined.

Final results — part 3.1
Shape function	Expression	Value at its own node
$N_1$ (node at $0,0$)	$(1-x/a)(1-y/b)$	1
$N_2$ (node at $a,0$)	$(x/a)(1-y/b)$	1
$N_3$ (node at $a,b$)	$xy/ab$	1
$N_4$ (node at $0,b$)	$(1-x/a)(y/b)$	1
Checks	$N_i(x_j,y_j)=\delta_{ij}$; $\sum N_i=1$; $\sum N_ix_i=x$, $\sum N_iy_i=y$

Topic. Generalised-coordinate derivation of the bilinear rectangle and the properties that make it a legitimate element.

Key relations. $u=\alpha_1+\alpha_2x+\alpha_3y+\alpha_4xy$; $\{\alpha\}=[A]^{-1}\{d\}$; $[N]=\lfloor 1\ x\ y\ xy\rfloor[A]^{-1}$; $N_i(x_j,y_j)=\delta_{ij}$; $\sum N_i=1$.

Why this works. Any interpolation that is complete to first order and continuous across element boundaries will converge, and the bilinear set achieves both cheaply. Along the edge $y=0$, $N_3$ and $N_4$ vanish and the remaining pair reduces to the linear bar interpolation, so two adjacent elements sharing that edge see the same displacement along it — the element is $C^{0}$ conforming. The $xy$ term is the price of having four nodes rather than three: it makes the strain field vary linearly rather than being constant, which is what gives the Q4 its advantage over the constant-strain triangle, and also what causes its shear-locking in bending-dominated meshes.

Common pitfall. Adding $x^{2}$ or $y^{2}$ to the polynomial to "improve" it. With only four nodes the system becomes unsolvable, and dropping $xy$ instead destroys the Kronecker-delta property and makes edges non-conforming. A second pitfall is renumbering the nodes clockwise, which flips the sign of the Jacobian and, in an isoparametric implementation, yields a negative element area.

Source. Logan, 6th ed., §10.1–10.2; Cook, Malkus, Plesha & Witt, 4th ed., §3.6 and §6.2; Bathe, Finite Element Procedures, 2nd ed., §4.2 and §5.3.

3.2 The plane-stress stiffness term $k_{33}$

Given. The same element, thickness $t$, modulus $E$, Poisson's ratio $\nu$, in plane stress; the degrees of freedom ordered $\{u_1,v_1,u_2,v_2,u_3,v_3,u_4,v_4\}$, so that degree of freedom 3 is $u_2$, the $x$-displacement of node 2.

Find. A closed-form expression for the diagonal stiffness term $k_{33}$.

Approach. Extract the third column of $[B]$, which involves only $N_2$; form the scalar $\left[B_3\right]^{T}[D]\left[B_3\right]$ for the plane-stress constitutive matrix; and integrate over the rectangle, where the double integral separates into two elementary one-dimensional integrals.

Write the strain–displacement column for degree of freedom 3. In plane elasticity $\{\varepsilon\}=\{\partial u/\partial x,\ \partial v/\partial y,\ \partial u/\partial y+\partial v/\partial x\}^{T}$. Degree of freedom 3 is $u_2$, which enters only $u$, so with $N_2=(x/a)(1-y/b)$, $$\left[B_3\right]=\begin{Bmatrix}\partial N_2/\partial x\\ 0\\ \partial N_2/\partial y\end{Bmatrix}=\begin{Bmatrix}\dfrac{1}{a}\left(1-\dfrac{y}{b}\right)\\[4pt] 0\\[4pt] -\dfrac{x}{ab}\end{Bmatrix}$$ Note that the direct strain from this degree of freedom varies with $y$ only, and the shear strain with $x$ only.
Write the plane-stress constitutive matrix. $$[D]=\frac{E}{1-\nu^{2}}\begin{bmatrix}1&\nu&0\\ \nu&1&0\\ 0&0&\dfrac{1-\nu}{2}\end{bmatrix}$$
Form the integrand. Because the second entry of $\left[B_3\right]$ is zero, the $\nu$ coupling drops out of a diagonal term and only two products survive: $$\left[B_3\right]^{T}[D]\left[B_3\right]=\frac{E}{1-\nu^{2}}\left[\left(\frac{\partial N_2}{\partial x}\right)^{2}+\frac{1-\nu}{2}\left(\frac{\partial N_2}{\partial y}\right)^{2}\right]$$ The first term is the direct-stretching contribution, the second the shearing contribution.
Integrate the direct-stretching term. $$\int_0^{b}\!\!\int_0^{a}\frac{1}{a^{2}}\left(1-\frac{y}{b}\right)^{2}dx\,dy=\frac{1}{a^{2}}\cdot a\cdot\frac{b}{3}=\frac{b}{3a}$$ using $\int_0^b(1-y/b)^{2}dy=b/3$.
Integrate the shearing term. $$\frac{1-\nu}{2}\int_0^{b}\!\!\int_0^{a}\frac{x^{2}}{a^{2}b^{2}}\,dx\,dy=\frac{1-\nu}{2}\cdot\frac{1}{a^{2}b^{2}}\cdot\frac{a^{3}}{3}\cdot b=\frac{\left(1-\nu\right)a}{6b}$$
Assemble. Multiplying by the thickness $t$ and the constitutive factor, $$\boxed{\;k_{33}=\left(\frac{b}{3a}+\frac{1-\nu}{6}\frac{a}{b}\right)\frac{Et}{1-\nu^{2}}\;}$$ which is the required expression. Two-by-two Gauss quadrature reproduces it exactly, because the integrand is at most quadratic in each variable.
Note the two companions needed for part 3.3. The identical procedure applied to degree of freedom 4 ($v_2$) and to degree of freedom 5 ($u_3$) gives $$k_{44}=\left(\frac{a}{3b}+\frac{1-\nu}{6}\frac{b}{a}\right)\frac{Et}{1-\nu^{2}},\qquad k_{55}=k_{33}$$ and the off-diagonal coupling terms are $k_{34}=-\dfrac{Et}{8\left(1-\nu\right)}$ and $k_{56}=+\dfrac{Et}{8\left(1-\nu\right)}$, equal in magnitude and opposite in sign. That sign reversal is used directly in the next part.

Final results — part 3.2
Term	Degree of freedom	Expression
$k_{33}$	$u_2$	$\left(\dfrac{b}{3a}+\dfrac{(1-\nu)a}{6b}\right)\dfrac{Et}{1-\nu^{2}}$
$k_{44}$	$v_2$	$\left(\dfrac{a}{3b}+\dfrac{(1-\nu)b}{6a}\right)\dfrac{Et}{1-\nu^{2}}$
$k_{55}$	$u_3$	equal to $k_{33}$
$k_{34}$, $k_{56}$	$u_2v_2$, $u_3v_3$	$\mp\,Et/\!\left[8(1-\nu)\right]$
Quadrature	—	$2\times2$ Gauss is exact for these terms

Topic. Explicit integration of a plane-stress element stiffness term, and what a single diagonal entry tells you about the element.

Key relations. $[k]=t\displaystyle\int_0^b\!\!\int_0^a [B]^{T}[D][B]\,dx\,dy$; plane stress $[D]$ as above; for a diagonal $u$-term only $\left(\partial N/\partial x\right)^{2}$ and $\frac{1-\nu}{2}\left(\partial N/\partial y\right)^{2}$ survive.

Why this works. On a rectangle the shape function derivatives are products of a constant and a linear function, so the integrand separates and the double integral collapses to a product of two elementary integrals — no Jacobian and no numerical quadrature are needed. Reading the result physically is worth the effort: the first term, $b/3a$, is a stretching stiffness that grows as the element becomes short and wide, while the second, $(1-\nu)a/6b$, is a shearing stiffness that grows as it becomes tall and thin. Their sum has a minimum near $a\approx b$, which is the analytical reason element aspect ratio matters and why highly distorted Q4 meshes become badly conditioned.

Common pitfall. Assuming the $\nu$ coupling contributes to a diagonal $u$-term. It cannot, because $[B_3]$ has a zero in the $\varepsilon_y$ row, so the $D_{12}$ entry is always multiplied by zero; candidates who carry a spurious $\nu$ term get $k_{33}$ wrong by exactly $\nu\,Et/(1-\nu^{2})$ times a geometric factor. A second trap is index confusion: $k_{33}$ is the third degree of freedom, which is $u_2$, not anything to do with node 3.

Source. Cook, Malkus, Plesha & Witt, 4th ed., §3.6 (explicit bilinear rectangle stiffness) and §6.5; Logan, 6th ed., §10.3; Chandrupatla & Belegundu, 4th ed., §7.2.

3.3 The bar-reinforced plate

Given. The plate of Figure 3(b): 8 cm wide and 12 cm high, thickness 0.1 cm, built in along its left, top and bottom edges and free on the right. It is meshed as two Q4 elements, each $a=L=8$ cm by $b=l=6$ cm, and reinforced along its horizontal centre line by a bar of area $A=0.5\ \text{cm}^{2}$ running from the built-in left edge to the free right edge. A horizontal force $P=6000$ N acts at node 4, the mid-height node of the free edge.

Given data
Quantity	Symbol	Value
Element width, height	$a$, $b$	8 cm, 6 cm
Plate thickness	$t$	0.1 cm
Reinforcing bar area, length	$A$, $L_b$	0.5 cm$^{2}$, 8 cm
Applied load at node 4	$P$	6000 N ($+x$)
Modulus of elasticity	$E$	$30\times10^{6}$ N/cm$^{2}$
Poisson's ratio	$\nu$	0.3
Point $C$ (from node 3)	$(x_C,y_C)$	$(L/4,\ 3l/4) = (2,\ 4.5)$ cm

Find. The horizontal displacement $u_4$ and the strains $\varepsilon_x$, $\varepsilon_y$, $\gamma_{xy}$ at point $C$.

Figure 3(b) — the reinforced plate. Nodes 1, 2 and 3 lie on built-in edges; node 4 on the free edge is the only unrestrained node, and it carries the load and the far end of the bar.

Approach. Identify the active degrees of freedom, assemble only the rows and columns that survive, show that symmetry about the centre line kills the coupling and the vertical displacement, solve one scalar equation for $u_4$, then push that displacement back through the element $[B]$ matrix at point $C$.

Count the active degrees of freedom. The support hatching runs along the left, top and bottom edges. Node 2 (top left) and node 1 (top right) sit on the built-in top edge, node 3 (mid left) on the built-in left edge, and both lower nodes of element 2 lie on the built-in left and bottom edges. Since a Q4 edge displacement is linear between its two end nodes, restraining the end nodes restrains the whole edge. Only node 4 is free, giving exactly two degrees of freedom, $u_4$ and $v_4$.
Identify where node 4 sits in each element. Taking the local origin of each element at its own lower-left corner, node 4 is the lower-right corner (local node 2, degrees of freedom 3 and 4) of the upper element and the upper-right corner (local node 3, degrees of freedom 5 and 6) of the lower element. Part 3.2 already supplies every term needed.
Assemble the $2\times2$ system. Adding the two element contributions, $$K_{uu}=k_{33}+k_{55}=2k_{33},\qquad K_{vv}=k_{44}+k_{66}=2k_{44},\qquad K_{uv}=k_{34}+k_{56}=0$$ The coupling cancels exactly, because the two contributions are equal and opposite (step 7 of part 3.2). That is the algebraic statement of the obvious physical fact: the plate is symmetric about the bar line and so is the load, so the loaded node cannot move vertically.
Evaluate the plate stiffness. With $a=8$, $b=6$, $t=0.1$, $\nu=0.3$, $$k_{33}=\left(\frac{6}{24}+\frac{0.7\times8}{36}\right)\frac{30\times10^{6}\times0.1}{1-0.09}=0.40556\times3.29670\times10^{6}=1.3370\times10^{6}\ \text{N/cm}$$ so the two plate elements together contribute $2k_{33}=2.6740\times10^{6}$ N/cm.
Add the reinforcing bar. The bar runs from the built-in node 3 to node 4 along $x$, so it is a spring acting on $u_4$ alone: $$k_{\text{bar}}=\frac{EA}{L_b}=\frac{30\times10^{6}\times0.5}{8}=1.8750\times10^{6}\ \text{N/cm}$$ It is not a small addition — a bar of 0.5 cm$^{2}$ is stiffer than a 0.1 cm thick plate 12 cm deep.
Solve for the displacements. The system is diagonal, so $$K_{uu}=2k_{33}+k_{\text{bar}}=4.5490\times10^{6}\ \text{N/cm}$$ $$\boxed{\;u_4=\frac{P}{K_{uu}}=\frac{6000}{4.5490\times10^{6}}=1.319\times10^{-3}\ \text{cm}=13.19\ \mu\text{m},\qquad v_4=0\;}$$ The bar carries $k_{\text{bar}}u_4=2473$ N, that is 41.2 % of the applied load, with the plate taking the remaining 3527 N.
Set up the strain recovery at $C$. Point $C$ lies in the upper element at $(x,y)=(2,\ 4.5)$ cm from that element's local origin at node 3. Of the eight element degrees of freedom only $u$ at local node 2 — that is, $u_4$ — is non-zero, so the whole displacement field of the element reduces to a single term: $$u(x,y)=N_2(x,y)\,u_4=\frac{x}{a}\left(1-\frac{y}{b}\right)u_4,\qquad v(x,y)=0$$
Differentiate for the strains. $$\varepsilon_x=\frac{\partial u}{\partial x}=\frac{1}{a}\left(1-\frac{y}{b}\right)u_4,\qquad \varepsilon_y=\frac{\partial v}{\partial y}=0,\qquad \gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}=-\frac{x}{ab}\,u_4$$ Substituting $x=2$, $y=4.5$, $a=8$, $b=6$ and $u_4=1.319\times10^{-3}$ cm, $$\boxed{\;\varepsilon_x=+4.122\times10^{-5},\qquad \varepsilon_y=0,\qquad \gamma_{xy}=-5.496\times10^{-5}\;}$$ $\varepsilon_y$ vanishes identically because no node of the element has a vertical displacement — not because the real plate has no vertical strain there, a distinction taken up in the concept note.
Convert to stresses as a check. With $[D]$ for plane stress, $$\sigma_x=1359\ \text{N/cm}^{2},\qquad \sigma_y=408\ \text{N/cm}^{2},\qquad \tau_{xy}=-634\ \text{N/cm}^{2}$$ $\sigma_y$ is non-zero even though $\varepsilon_y$ is zero: the Poisson coupling in $[D]$ turns the $\varepsilon_x$ stretching into a transverse stress once the element is prevented from contracting.

Final results — part 3.3
Quantity	Symbol	Value
Plate stiffness per element at node 4	$k_{33}$	$1.337\times10^{6}$ N/cm
Bar stiffness	$EA/L_b$	$1.875\times10^{6}$ N/cm
Assembled horizontal stiffness	$K_{uu}=2k_{33}+EA/L_b$	$4.549\times10^{6}$ N/cm
Horizontal displacement of node 4	$u_4$	$1.319\times10^{-3}$ cm (13.19 µm)
Vertical displacement of node 4	$v_4$	0 (by symmetry)
Load carried by the bar	$F_{\text{bar}}$	2473 N (41.2 %)
Load carried by the plate	$F_{\text{plate}}$	3527 N (58.8 %)
Direct strain at $C$	$\varepsilon_x$	$+4.122\times10^{-5}$
Transverse strain at $C$	$\varepsilon_y$	0
Shear strain at $C$	$\gamma_{xy}$	$-5.496\times10^{-5}$
Stresses at $C$	$\sigma_x,\sigma_y,\tau_{xy}$	1359, 408, −634 N/cm$^{2}$

Topic. Assembling a mixed plate-plus-bar model, exploiting symmetry to decouple the equations, and recovering strains from nodal displacements.

Key relations. $K_{uu}=\sum_e k_{uu}^{(e)}+EA/L_b$; symmetry ⇒ $K_{uv}=0$ and $v_4=0$; $u_4=P/K_{uu}$; $\{\varepsilon\}=[B]\{d\}$ with $\varepsilon_x=\frac{1}{a}(1-y/b)u_4$, $\gamma_{xy}=-\frac{x}{ab}u_4$; $\{\sigma\}=[D]\{\varepsilon\}$.

Why this works. A bar element and a plane-stress element assemble into the same global matrix simply because both are derived from the same variational statement — the bar contributes to translational degrees of freedom only, and adding $EA/L_b$ to the $u_4$ diagonal is the whole of its effect here. The symmetry argument then does the rest: because node 4 appears at the lower-right corner of one element and the upper-right corner of its mirror image, the two $u$–$v$ coupling terms are equal and opposite and cancel identically, so one hand-solvable scalar equation replaces a $2\times2$ solve. The result is worth reading as design information rather than just arithmetic: a bar of 0.5 cm$^{2}$ takes 41 % of the load from a 0.1 cm plate, which is why edge stiffeners are so effective in thin-plate structures.

Common pitfall. Treating the plate as though only the upper element existed, which halves the plate stiffness and inflates $u_4$ by about 30 %; and reporting $\varepsilon_y=0$ as a physical result. It is a modelling result: this two-element mesh simply has no degree of freedom capable of producing a vertical strain at $C$, and a refined mesh would show a small non-zero $\varepsilon_y$. A third pitfall is to measure $C$ from the wrong origin — the dimensions $L/4$ and $3l/4$ are taken from node 3, the mid-height point on the built-in left edge, not from a corner of the plate.

Source. Logan, 6th ed., §6.5 (plane-stress stress recovery) and Ch. 10 for the Q4; Cook, Malkus, Plesha & Witt, 4th ed., §2.7 (assembly of dissimilar elements) and §8.7 (symmetry); Chandrupatla & Belegundu, 4th ed., §7.3–7.4; Bathe, 2nd ed., §4.2.

Question 3 of 3: Q4 shape functions, plane-stress stiffness term, and a reinforced plate

Question 3: Q4 shape functions, plane-stress stiffness term, and a reinforced plate (equal value)

3.1 Derivation of the bilinear shape functions

3.2 The plane-stress stiffness term $k_{33}$

3.3 The bar-reinforced plate