The linear triangular element — stiffness block, plate stresses and local equilibrium: 16-Civ-B9 The Finite Element Method, May 2015, Q3

Problem 3: The linear triangular element — stiffness block, plate stresses and local equilibrium (one of three equal-value problems)

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.

Part 3.1 — deriving the block $[k_{33}]$

Given. The constant strain-displacement matrix $[B]$ printed above, the plane-stress constitutive matrix, and an element of thickness $t$ and area $A_{T}$.

Find. A closed form for the lower-right $2\times2$ block of $[k]=t\,A_{T}[B]^{T}[D][B]$, and a demonstration that it equals the printed expression.

Approach. Partition $[B]$ by node, so the block $[k_{ij}]$ is $t\,A_{T}[B_{i}]^{T}[D][B_{j}]$; then multiply out the single case $i=j=3$ with the plane-stress $[D]$.

Partition the strain-displacement matrix by node. Reading the printed $[B]$ column-pair by column-pair, $$[B]=\frac{1}{2A_{T}}\big[\,[\hat{B}_{1}]\;\big|\;[\hat{B}_{2}]\;\big|\;[\hat{B}_{3}]\,\big], \qquad [\hat{B}_{i}]=\begin{bmatrix}\beta_{i}&0\\0&\gamma_{i}\\\gamma_{i}&\beta_{i}\end{bmatrix}.$$ Because $[B]$ is constant over the element — a consequence of the linear displacement field — the volume integral in the stiffness definition is just a multiplication by the volume $t A_{T}$: $$[k]=\int_{V}[B]^{T}[D][B]\,dV=t\,A_{T}\,[B]^{T}[D][B].$$
Write the general block. Substituting the partition and taking the $1/(2A_{T})$ factors out of both $[B]^{T}$ and $[B]$, $$[k_{ij}]=t\,A_{T}\frac{1}{(2A_{T})^{2}}[\hat{B}_{i}]^{T}[D][\hat{B}_{j}] =\frac{t}{4A_{T}}[\hat{B}_{i}]^{T}[D][\hat{B}_{j}].$$ Setting $i=j=3$ leaves only $\beta_{3}$ and $\gamma_{3}$ in play, which is why the printed result contains no other coefficient.
Insert the plane-stress constitutive matrix. 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},$$ so $$[D][\hat{B}_{3}]=\frac{E}{1-\nu^{2}} \begin{bmatrix}\beta_{3}&\nu\gamma_{3}\\ \nu\beta_{3}&\gamma_{3}\\ \dfrac{1-\nu}{2}\gamma_{3}&\dfrac{1-\nu}{2}\beta_{3}\end{bmatrix}.$$
Premultiply by $[\hat{B}_{3}]^{T}$ and collect terms. With $[\hat{B}_{3}]^{T}=\begin{bmatrix}\beta_{3}&0&\gamma_{3}\\0&\gamma_{3}&\beta_{3}\end{bmatrix}$, the $(1,1)$ entry is $\beta_{3}^{2}+\frac{1-\nu}{2}\gamma_{3}^{2}$ and the $(2,2)$ entry is $\gamma_{3}^{2}+\frac{1-\nu}{2}\beta_{3}^{2}$. The off-diagonal entry needs one line of algebra: $$\nu\beta_{3}\gamma_{3}+\frac{1-\nu}{2}\gamma_{3}\beta_{3} =\gamma_{3}\beta_{3}\left(\nu+\frac{1-\nu}{2}\right) =\gamma_{3}\beta_{3}\left(\frac{2\nu+1-\nu}{2}\right) =\left(\frac{1+\nu}{2}\right)\gamma_{3}\beta_{3},$$ and the same value appears in the $(2,1)$ position, so the block is symmetric as required.
Assemble the result. Restoring the factor $t/(4A_{T})$ and the $E/(1-\nu^{2})$ from $[D]$, $$\boxed{\;[k_{33}]=\frac{Et}{4A_{T}(1-\nu^{2})}\begin{bmatrix} \beta_{3}^{2}+\dfrac{1-\nu}{2}\gamma_{3}^{2} & \left(\dfrac{1+\nu}{2}\right)\gamma_{3}\beta_{3}\\[6pt] \left(\dfrac{1+\nu}{2}\right)\gamma_{3}\beta_{3} & \gamma_{3}^{2}+\dfrac{1-\nu}{2}\beta_{3}^{2}\end{bmatrix}\;}$$ which is the printed expression.

Part 3.2 — stresses in the plate

Given.

Problem 3.2 — data from the question and Figure 3-b
Quantity	Symbol	Value
Base (horizontal leg)	—	30 mm
Height (vertical leg, built in)	—	20 mm
Thickness	$t$	10 mm
Modulus of elasticity	$E$	70 GPa $=70\,000$ MPa
Poisson ratio	$\nu$	0.3
Pressure on the hypotenuse, normal and inward	$p$	100 MPa
Node 1 (top of the built-in edge)	$(x_{1},y_{1})$	$(0,\,20)$ mm, fully fixed
Node 2 (bottom of the built-in edge)	$(x_{2},y_{2})$	$(0,\,0)$ mm, fully fixed
Node 3 (free corner)	$(x_{3},y_{3})$	$(30,\,0)$ mm, free

Find. The constant stress components $\sigma_{x}$, $\sigma_{y}$ and $\tau_{xy}$ that a single constant-strain triangle predicts for the plate.

Figure 3.1 — The plate of Figure 3-b as a single linear triangular element. Numbering the nodes counter-clockwise starting at the top of the fixed edge puts the one free node at position 3, so the block derived in 3.1 is exactly the one needed.

Approach. Choose the node numbering so that the free corner is node 3; compute $\beta_{3}$, $\gamma_{3}$ and $A_{T}$; assemble the $2\times2$ system $[k_{33}]\{u_{3}\,v_{3}\}^{T}=\{f_{3}\}$ with $\{f_{3}\}$ the half of the edge pressure that the consistent load vector sends to node 3; solve, then run the displacements back through $[B]$ and $[D]$.

Number the nodes and confirm the area. Taking node 1 at $(0,20)$, node 2 at $(0,0)$ and node 3 at $(30,0)$ traverses the triangle counter-clockwise, so the signed area is positive: $$A_{T}=\tfrac{1}{2}\big[(x_{2}-x_{1})(y_{3}-y_{1})-(x_{3}-x_{1})(y_{2}-y_{1})\big] =\tfrac{1}{2}\big[(0)(-20)-(30)(-20)\big]=300\ \text{mm}^{2}.$$ The whole of the left edge is built in, so nodes 1 and 2 are fully restrained and node 3 carries the only two active degrees of freedom in the model.
Evaluate the geometric coefficients for node 3. From the definitions printed in the question, $$\beta_{3}=y_{1}-y_{2}=20-0=20\ \text{mm},\qquad \gamma_{3}=x_{2}-x_{1}=0-0=0.$$ That $\gamma_{3}$ vanishes is a gift: it uncouples the two equations and makes $[k_{33}]$ diagonal.
Form the stiffness block. The scalar in front is $$\frac{Et}{4A_{T}(1-\nu^{2})}=\frac{(70\,000)(10)}{4(300)(1-0.3^{2})} =\frac{700\,000}{1092}=641.03\ \text{N}/\text{mm}^{3},$$ and with $\gamma_{3}=0$ the printed form collapses to $$[k_{33}]=641.03\begin{bmatrix}\beta_{3}^{2}&0\\0&\dfrac{1-\nu}{2}\beta_{3}^{2}\end{bmatrix} =641.03\begin{bmatrix}400&0\\0&140\end{bmatrix} =\begin{bmatrix}256\,410&0\\0&89\,744\end{bmatrix}\ \text{N}/\text{mm}.$$
Resolve the pressure into a nodal load vector. The loaded edge runs from node 1 to node 3, so its direction is $(30,-20)$ mm and its length is $$L_{13}=\sqrt{30^{2}+20^{2}}=\sqrt{1300}=36.056\ \text{mm}.$$ The outward unit normal to that edge is $\{n_{x},n_{y}\}=\{20,30\}/36.056=\{0.5547,\,0.8321\}$, and a pressure pushes against the outward normal, so the traction is $\{t\}=-p\{n\}$. Multiplying by the loaded area $L_{13}t$ gives the total edge force $$\{F\}_{\text{edge}}=-p\,\{n\}\,L_{13}\,t=\{-20\,000,\;-30\,000\}\ \text{N}=\{-20.0,\;-30.0\}\ \text{kN}.$$ Note the tidy shortcut: $p\,n_{x}L_{13}t=p(20)t$ and $p\,n_{y}L_{13}t=p(30)t$, because the normal components times the hypotenuse length are just the projected leg lengths.
Send half of it to the free node. For a linear triangle the shape functions vary linearly along an edge, so a uniform edge traction splits equally between the two nodes on that edge. Node 1 is restrained and its half goes straight into the reactions; node 3 receives $$\{f_{3}\}=\tfrac{1}{2}\{F\}_{\text{edge}}=\{-10\,000,\;-15\,000\}\ \text{N}.$$
Solve for the displacement of node 3. Because $[k_{33}]$ is diagonal the two equations separate: $$u_{3}=\frac{-10\,000}{256\,410}=-0.0390\ \text{mm},\qquad v_{3}=\frac{-15\,000}{89\,744}=-0.16714\ \text{mm},$$ $$\boxed{\;u_{3}=-0.0390\ \text{mm},\qquad v_{3}=-0.1671\ \text{mm}\;}$$ Both are negative: the corner moves left and down, into the plate and away from the load, as the picture demands.
Recover the strains. With every other nodal displacement zero, the constant strains follow from the first, second and third rows of $[B]$ using only the node-3 columns: $$\varepsilon_{x}=\frac{\beta_{3}u_{3}}{2A_{T}}=\frac{20(-0.0390)}{600}=-1.300\times10^{-3},$$ $$\varepsilon_{y}=\frac{\gamma_{3}v_{3}}{2A_{T}}=0,\qquad \gamma_{xy}=\frac{\gamma_{3}u_{3}+\beta_{3}v_{3}}{2A_{T}}=\frac{20(-0.16714)}{600}=-5.5714\times10^{-3}.$$ The exactly zero $\varepsilon_{y}$ is a direct consequence of $\gamma_{3}=0$, not a rounding artefact: with the whole left edge clamped, a single element simply has no mechanism for vertical stretching.
Convert to stresses. Applying the plane-stress law $\{\sigma\}=[D]\{\varepsilon\}$ with $E/(1-\nu^{2})=70\,000/0.91=76\,923$ MPa, $$\sigma_{x}=76\,923\,(\varepsilon_{x}+\nu\varepsilon_{y})=76\,923(-1.300\times10^{-3})=-100.0\ \text{MPa},$$ $$\sigma_{y}=76\,923\,(\nu\varepsilon_{x}+\varepsilon_{y})=76\,923(-0.390\times10^{-3})=-30.0\ \text{MPa},$$ $$\tau_{xy}=G\gamma_{xy}=\frac{E}{2(1+\nu)}\gamma_{xy}=26\,923(-5.5714\times10^{-3})=-150.0\ \text{MPa},$$ $$\boxed{\;\sigma_{x}=-100.0\ \text{MPa},\qquad \sigma_{y}=-30.0\ \text{MPa},\qquad \tau_{xy}=-150.0\ \text{MPa}\;}$$ These three numbers are the whole answer: a constant-strain triangle carries one stress state across the entire element, so this is the predicted stress distribution.
Interpret the state. The principal stresses follow from Mohr's circle with centre $(\sigma_{x}+\sigma_{y})/2=-65.0$ MPa and radius $\sqrt{35^{2}+150^{2}}=154.03$ MPa: $$\sigma_{1}=+89.0\ \text{MPa},\qquad \sigma_{2}=-219.0\ \text{MPa},\qquad \tau_{\max}=154.0\ \text{MPa},$$ on planes at $\theta_{p}=-51.6^{\circ}$ to the $x$ axis, with a von Mises equivalent stress of 274.6 MPa. The dominant action is a large shear — unsurprising for a cantilevered wedge pushed on its sloping face — and the element is telling us the plate is shear-critical near the built-in edge.

Check: two readings are taken from Figure 3-b rather than from the printed text. First, the hatching runs the full height of the left edge, so both nodes on that edge are taken as fully fixed; a support at only one node would leave four active degrees of freedom and different numbers. Second, the load arrows are drawn perpendicular to the hypotenuse, so $p$ is treated as a pressure normal to that face rather than as a vertical or horizontal traction. Both readings are the ones that make the problem solvable with a single element, which is what the question asks for.

Part 3.3 — local equilibrium of the approximated stress field

The question is a trap worth walking into deliberately, because the honest answer has two halves that point in opposite directions.

Inside the element, the differential equations of equilibrium are satisfied identically — and trivially. The two-dimensional equilibrium equations with no body force are

$$\frac{\partial\sigma_{x}}{\partial x}+\frac{\partial\tau_{xy}}{\partial y}=0,\qquad \frac{\partial\tau_{xy}}{\partial x}+\frac{\partial\sigma_{y}}{\partial y}=0,$$

and the constant-strain triangle delivers a stress field in which every one of those derivatives is zero. So the element passes the internal test, but it passes it for no good reason: it would pass whatever numbers came out of the analysis. A constant field can never detect an equilibrium violation, which is exactly why single-element CST results must be treated with suspicion.

On the boundaries, equilibrium fails badly. The real test is whether the computed stress reproduces the applied traction. On the hypotenuse, with outward normal $\{n\}=\{0.5547,\,0.8321\}$, Cauchy's relation gives

$$t_{x}=\sigma_{x}n_{x}+\tau_{xy}n_{y}=-100(0.5547)-150(0.8321)=-180.3\ \text{MPa},$$ $$t_{y}=\tau_{xy}n_{x}+\sigma_{y}n_{y}=-150(0.5547)-30(0.8321)=-108.2\ \text{MPa},$$

whereas the applied traction is $-p\{n\}=\{-55.5,\,-83.2\}$ MPa. The residual is $\{-124.8,\,-25.0\}$ MPa, of magnitude 127.3 MPa — larger than the applied pressure itself. The free bottom edge is worse in relative terms: with outward normal $\{0,-1\}$ the computed traction is $\{-\tau_{xy},\,-\sigma_{y}\}=\{+150,\,+30\}$ MPa on a face that should be entirely traction-free.

What is satisfied is the weak, integrated statement. Summing the computed tractions around the closed boundary of the element gives exactly zero resultant force, and the nodal forces recovered as $t A_{T}[B]^{T}\{\sigma\}$ reproduce the applied load at node 3 and the reactions at nodes 1 and 2 to machine precision. In other words the element is in global equilibrium, as the assembled equations $[K]\{d\}=\{F\}$ guarantee, while being pointwise wrong everywhere on its surface. This is precisely the phenomenon described in Part 1 of Problem 1, made numerical.

The practical conclusion is that a one-element idealisation of this plate is a demonstration, not an analysis. The true stress field in a wedge loaded on one face varies strongly with position and is singular at the re-entrant corner where the loaded edge meets the clamped edge; no constant field can represent it. Refining to a mesh of many CSTs, or switching to a linear-strain triangle (six-node LST) or a bilinear quadrilateral, reduces the inter-element traction jumps and the boundary residual towards zero. Monitoring exactly that residual — the discontinuity in recovered stress between neighbouring elements — is what error estimators do, and it is the standard trigger for adaptive mesh refinement.

Problem 3 — results
Quantity	Symbol	Value
Triangle area	$A_{T}$	300 mm$^{2}$
Geometric coefficients at node 3	$\beta_{3}$, $\gamma_{3}$	20 mm, 0
Stiffness block	$[k_{33}]$	diag(256 410, 89 744) N/mm
Consistent nodal load at node 3	$\{f_{3}\}$	$(-10.0,\,-15.0)$ kN
Displacement of node 3	$u_{3}$, $v_{3}$	$-0.0390$ mm, $-0.1671$ mm
Strains (constant)	$\varepsilon_{x}$, $\varepsilon_{y}$, $\gamma_{xy}$	$-1.300\times10^{-3}$, 0, $-5.5714\times10^{-3}$
Direct stress, $x$	$\sigma_{x}$	$-100.0$ MPa
Direct stress, $y$	$\sigma_{y}$	$-30.0$ MPa
Shear stress	$\tau_{xy}$	$-150.0$ MPa
Principal stresses	$\sigma_{1}$, $\sigma_{2}$	$+89.0$ MPa, $-219.0$ MPa at $-51.6^{\circ}$
Maximum in-plane shear / von Mises	$\tau_{\max}$, $\sigma_{vM}$	154.0 MPa, 274.6 MPa
Traction residual on the loaded edge	—	127.3 MPa (applied 100 MPa)
Internal differential equilibrium	—	satisfied identically (constant field)
Boundary traction conditions	—	violated on both the loaded and the free edge

Topic. The constant-strain triangle in plane stress — from the geometric coefficients through the stiffness block and the consistent edge load to the element stress state, and what that state does and does not satisfy.

Key relations. $[k]=tA_{T}[B]^{T}[D][B]$ with $[B]$ constant; $[k_{ij}]=\dfrac{t}{4A_{T}}[\hat{B}_{i}]^{T}[D][\hat{B}_{j}]$; plane stress $[D]=\dfrac{E}{1-\nu^{2}}\begin{bmatrix}1&\nu&0\\\nu&1&0\\0&0&\frac{1-\nu}{2}\end{bmatrix}$; consistent edge load $\{f\}=\frac{1}{2}(-p\{n\})L_{e}t$ split equally between the two edge nodes; Cauchy traction $t_{i}=\sigma_{ij}n_{j}$; equilibrium $\partial\sigma_{ij}/\partial x_{j}+b_{i}=0$.

Why this works. A linear displacement field over a triangle differentiates to a constant strain, which is what makes $[B]$ constant and lets the volume integral in the stiffness definition become a plain multiplication by $tA_{T}$ — the reason the CST is the one two-dimensional element whose stiffness can be written out by hand. The same constancy is what makes the block structure of part 3.1 so clean: each $2\times2$ block sees only its own node's $\beta$ and $\gamma$. In part 3.2 the geometry does the rest of the work. Numbering counter-clockwise with the free corner last makes the required block the one already derived; the vertical clamped edge makes $\gamma_{3}=0$, which diagonalises $[k_{33}]$ and kills $\varepsilon_{y}$; and the equal split of the edge traction follows directly from the linear variation of the shape functions along an edge. The answer to part 3.3 then falls out of the same constancy from the other side: a field with no gradient cannot violate the interior equations, and cannot satisfy any traction condition except by coincidence.

Common pitfall. Numbering the nodes clockwise, which makes $A_{T}$ negative and flips the sign of every stress; always check the signed area is positive before going on. The second is resolving the pressure into components by eye — the total edge force is $p\,L_{e}t$ along the normal, so its $x$ component is $p\,t\times$(the vertical projection, 20 mm) and its $y$ component is $p\,t\times$(the horizontal projection, 30 mm), not $p\,t\times$(the hypotenuse) in each direction. A third is answering 3.3 with a flat “equilibrium is not satisfied”; the marks are in distinguishing the interior equations, which hold trivially, from the boundary conditions, which do not.

Source. Logan, A First Course in the Finite Element Method, 6th ed., Ch. 6 (the constant-strain triangle, plane stress and plane strain, equivalent nodal forces on element edges); Cook, Malkus, Plesha and Witt, Concepts and Applications of Finite Element Analysis, 4th ed., Ch. 3 and Ch. 4 (patch test, stress recovery and error estimation); Zienkiewicz, Taylor and Zhu, The Finite Element Method: Its Basis and Fundamentals, 7th ed., Ch. 4 and Ch. 14.

Question 3 of 3: The linear triangular element — stiffness block, plate stresses and local equilibrium

Problem 3: The linear triangular element — stiffness block, plate stresses and local equilibrium (one of three equal-value problems)

Part 3.1 — deriving the block $[k_{33}]$

Part 3.2 — stresses in the plate

Part 3.3 — local equilibrium of the approximated stress field