07-Str-B3 · May 2015
Nivaar worked solution (AI-drafted; not reviewed by a licensed engineer)
Paper format. Professional Engineers Ontario / Engineers Canada national examination 07-Str-B3 — Applications of the Finite Element Method, May 2015. Three hours; four pages. Closed book, with two 8½ × 11 in pages of handwritten notes permitted and one approved non-communicating calculator. Three problems, all of equal value; candidates are instructed to attempt all three. Problem 1 is a ten-part concept paper, Problems 2 and 3 are calculations.
07-Str-B3 is a finite element methods paper rather than a member-design paper: bar and beam elements, isoparametric quadrilaterals, the constant-strain triangle, numerical integration and structural dynamics. The reference list below is therefore the finite-element literature.
Reference texts.
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.
The displacement-based finite element method starts from an assumed displacement field $\mathbf{u}^{h} = \mathbf{N}\mathbf{d}$ built from element shape functions. That field is made to satisfy compatibility exactly — it is single-valued, continuous across element boundaries and satisfies the essential (displacement) boundary conditions. Equilibrium, by contrast, is never imposed point by point. It enters only through the principle of virtual work, which is the weak (Galerkin, weighted-residual) form of the equilibrium equations:
$$\int_{V}\delta\boldsymbol{\varepsilon}^{T}\boldsymbol{\sigma}\,dV=\int_{V}\delta\mathbf{u}^{T}\mathbf{b}\,dV+\int_{S_{t}}\delta\mathbf{u}^{T}\bar{\mathbf{t}}\,dS$$Substituting the trial field turns this into $\mathbf{K}\mathbf{d}=\mathbf{F}$. What that statement enforces is that the residual of the differential equilibrium equations is orthogonal to the trial space, not that it vanishes:
$$\int_{V}\mathbf{N}^{T}\bigl(\nabla\cdot\boldsymbol{\sigma}^{h}+\mathbf{b}\bigr)\,dV=\mathbf{0}\qquad\text{but}\qquad\nabla\cdot\boldsymbol{\sigma}^{h}+\mathbf{b}\neq\mathbf{0}\ \text{pointwise.}$$Three consequences follow. First, the stresses are obtained by differentiating an approximate displacement field, $\boldsymbol{\sigma}^{h} = \mathbf{D}\mathbf{B}\mathbf{d}$, which lowers the polynomial order by one; a constant-strain triangle therefore offers a constant stress field and a bilinear quadrilateral a stress field that is at best linear in one direction. Such a low-order field simply has too few degrees of freedom to balance a general body force at every point. Second, because only $C^{0}$ continuity of displacement is enforced, the tractions computed on the two sides of an inter-element boundary do not agree: the mesh has a jump $[\![\boldsymbol{\sigma}\cdot\mathbf{n}]\!]\neq\mathbf{0}$, so equilibrium fails across element faces as well as inside them. Third, the natural (traction) boundary conditions are themselves satisfied only in the work-equivalent sense — the applied pressure is replaced by nodal forces that do the same work, not by a stress field that reproduces it.
What is satisfied exactly is equally worth stating: nodal equilibrium of the assembled system, global equilibrium of the whole body (the reactions always balance the applied loads), and the equilibrium of every individual element under its own work-equivalent nodal loads. Local equilibrium is recovered only in the limit of mesh refinement ($h\to0$) or of increasing polynomial order ($p\to\infty$); the size of the residual is in fact the basis of the residual-based error estimators used for adaptive meshing. Problem 3.3 of this same paper is a concrete illustration of the effect.
The correct continuation is the third box: the strain energy is lower than the exact value.
The argument is the minimum-total-potential-energy theorem. For a linear elastic body under prescribed loads the exact displacement field minimises $\Pi = U - W$ over all kinematically admissible fields. The finite element trial space is a subspace of that admissible set, so the minimum found within it can never be lower than the true minimum:
$$\Pi_{\text{FE}}\ \ge\ \Pi_{\text{exact}}$$At equilibrium the external work is twice the strain energy, $W = 2U$, so $\Pi = U - 2U = -U$. Substituting gives $-U_{\text{FE}}\ge-U_{\text{exact}}$, that is
$$\boxed{\,U_{\text{FE}}\le U_{\text{exact}}\,}$$Physically, restricting the displacement to a finite set of shape functions is an extra kinematic constraint. Constraints stiffen a structure, so the model is stiffer than the real thing, the computed displacements are too small, and the energy stored is too small with them. This is why a displacement-based mesh converges to the exact answer from below and why refining the mesh always softens the model.
Check: the bound has conditions. It holds for a conforming, fully integrated, displacement-based element with a consistent load vector and load (natural) boundary conditions. Reduced integration, non-conforming or incompatible modes, and hybrid or mixed formulations all remove the guarantee, and a model driven by prescribed displacements reverses the inequality because then the finite element energy is an upper bound.
Figure 1-a — the 2-node bar. With two nodes the axial displacement can only be linear, $u(x) = a_{0} + a_{1}x$. Imposing $u(0)=u_{1}$ and $u(L)=u_{2}$ gives the familiar pair of straight lines, each equal to one at its own node and zero at the other:
$$N_{1}(x)=1-\frac{x}{L},\qquad N_{2}(x)=\frac{x}{L},\qquad N_{1}+N_{2}=1$$The strain $\varepsilon = du/dx = (u_{2}-u_{1})/L$ is constant, which is why this element is the constant-strain bar.
Figure 1-b — the 3-node bar with the interior node at $L/3$. Three nodes support a complete quadratic, $u(x)=a_{0}+a_{1}x+a_{2}x^{2}$. The Lagrange construction gives each shape function as the product of the factors that vanish at the other two nodes. Writing $s = x/L$ and placing node 3 at $s = 1/3$:
$$N_{1}(s)=\frac{(s-\tfrac{1}{3})(s-1)}{(0-\tfrac{1}{3})(0-1)}=3s^{2}-4s+1,\qquad N_{3}(s)=\frac{s(s-1)}{\tfrac{1}{3}(\tfrac{1}{3}-1)}=\tfrac{9}{2}s(1-s),\qquad N_{2}(s)=\frac{s(s-\tfrac{1}{3})}{1\cdot(1-\tfrac{1}{3})}=\tfrac{3}{2}s^{2}-\tfrac{1}{2}s$$These are the three parabolas drawn on the right of the figure. Because the interior node is not at mid-length, the curves are visibly skewed and two features must appear in a correct sketch. $N_{1}$ falls steeply, crosses zero at $s=1/3$ and reaches a minimum of $-1/3$ at $s=2/3$ before climbing back to zero at the far end; $N_{3}$ is pushed to the right of its own node and peaks at $1.125$ at mid-length, overshooting unity; $N_{2}$ dips slightly negative (to $-1/24$ at $s=1/6$) before rising to one at node 2. All three still satisfy the two properties every shape function set must have: $N_{i}(x_{j})=\delta_{ij}$, and $\sum N_{i}=1$ everywhere, which is what lets the element represent a rigid-body translation. The strain now varies linearly along the element, so this is the linear-strain bar.
A hybrid element is one derived from a multi-field variational principle in which an extra field is assumed independently of the displacements, and in which that extra field is confined to the element and eliminated before assembly. The archetype is Pian's hybrid-stress element: a self-equilibrating stress field is assumed inside the element,
$$\boldsymbol{\sigma}=\mathbf{P}\boldsymbol{\beta}\quad\text{with}\quad\nabla\cdot\boldsymbol{\sigma}=\mathbf{0},$$while a compatible displacement field $\mathbf{u}=\mathbf{N}\mathbf{q}$ is assumed only on the element boundary, so that continuity with the neighbours is preserved. Substituting both into the Hellinger–Reissner (complementary energy) functional and making it stationary with respect to $\boldsymbol{\beta}$ gives $\boldsymbol{\beta} = \mathbf{H}^{-1}\mathbf{G}\mathbf{q}$, and static condensation at element level leaves an ordinary stiffness matrix in the nodal displacements alone:
$$\mathbf{k}=\mathbf{G}^{T}\mathbf{H}^{-1}\mathbf{G},\qquad \mathbf{H}=\int_{V}\mathbf{P}^{T}\mathbf{D}^{-1}\mathbf{P}\,dV,\qquad\mathbf{G}=\int_{S}\mathbf{R}^{T}\mathbf{N}\,dS$$The pay-off is accuracy on coarse meshes: the element is far less sensitive to geometric distortion, it is free of shear and volumetric locking, it passes the patch test, and the stresses come out directly from $\boldsymbol{\beta}$ instead of by differentiating an approximate displacement field, which is the least accurate step of the ordinary formulation. The price is a more elaborate derivation and the need to choose $\mathbf{P}$ with exactly the right number of parameters — too many and the element locks again, too few and it develops spurious mechanisms.
The distinction from a mixed element is worth making explicit, because the two terms are often blurred. In a mixed element (for example a displacement–pressure element for near-incompressible media) both fields are retained as global unknowns and appear in the assembled system. In a hybrid element the extra field is internal and is condensed out, so the assembled system looks exactly like a conventional displacement model.
Reduced integration means evaluating the element stiffness integral
$$\mathbf{k}=\int_{-1}^{1}\!\!\int_{-1}^{1}\mathbf{B}^{T}\mathbf{D}\mathbf{B}\,t\,|\mathbf{J}|\,d\xi\,d\eta$$with a Gauss rule of lower order than the one needed to integrate the integrand exactly. For a bilinear quadrilateral the exact rule is $2\times2$; the reduced rule is a single point at the centroid. For the 8-node quadratic quadrilateral the exact rule is $3\times3$ and the reduced rule is $2\times2$.
The obvious motive is cost, but the real motive is accuracy. Under-integration deliberately softens the element by discarding the higher-order part of the strain energy — precisely the part that is spurious. A fully integrated linear element in bending develops a parasitic transverse shear strain that does not exist physically, and its stiffness error grows as $(L/h)^{2}$; this is shear locking. The same fully integrated element in a near-incompressible material ($\nu\to0.5$) over-constrains the volumetric strain at every Gauss point and produces volumetric locking. Sampling at fewer points relaxes both constraints and restores sensible behaviour. A related bonus is that the reduced Gauss points are the Barlow (superconvergent) points at which the sampled stresses are one order more accurate than anywhere else in the element, which is why commercial codes report stresses there and extrapolate outward.
The cost is rank deficiency. One-point integration of a bilinear quadrilateral leaves the stiffness matrix with more zero eigenvalues than the three rigid-body modes, and the extra ones are spurious zero-energy (hourglass) mechanisms that can propagate through the mesh and destroy the solution. The remedies in use are hourglass control (adding a small artificial stiffness that resists only the spurious pattern) and selective reduced integration, in which the integrand is split and each part gets its own rule: full integration on the bending or deviatoric part, reduced integration on the shear or volumetric part. Selective schemes are the standard cure because they avoid mechanisms altogether.
[Figure not reproduced: Figure 1-c (redrawn) — two 4-node elements (1, 2) abutting an 8-node element (3), which in turn abuts a further element (4) over only half of its right-hand edge. The two dashed red interfaces are where the connection is defective. See the official exam paper.]
Defect A — a linear edge forced to match a parabolic edge. Along the vertical interface the two 4-node elements 1 and 2 vary linearly in each of their own edges, so the displacement along the whole interface is piecewise linear with a kink at the shared corner node. The 8-node element 3 varies parabolically along the same interface. The three nodes coincide, so the two fields agree at the nodes, but nowhere in between: a gap opens over one half of the interface and the material overlaps over the other. Inter-element $C^{0}$ continuity, which the convergence proof relies on, is broken; the assembly is a non-conforming mesh, it fails the patch test, monotone convergence is lost, and a spurious stress concentration appears along the interface no matter how fine the mesh becomes.
Defect B — the mid-side node is a hanging node. The mid-side node on the left edge of element 3 happens to coincide with the corner shared by elements 1 and 2, so it is at least connected; but the mid-side nodes on the right edge of element 3 and on the top edge of element 4 have no counterpart in the neighbouring element. A node that lies on the interior of an adjacent element's edge but is not a node of that element carries degrees of freedom that the neighbour cannot see. It can therefore displace independently: the mesh contains an unintended internal slit through which forces are not transmitted, material can inter-penetrate, and the assembled stiffness under-represents the true connection.
Defect C — element 4 is offset. Element 4 covers only the lower half of the right-hand edge of element 3. Half of that edge therefore behaves as a free surface inside the body, and the corner node of element 4 sits at the mid-length of element 3's edge. This is the worst of the three: it is a topological disconnection, not merely an interpolation mismatch, and it will show up as a visible discontinuity in the deformed plot.
The accepted remedies, in order of preference: mesh the region with one element family so that the transition never occurs; use a genuine transition element (a 5-, 6- or 7-node quadrilateral formed by deleting the offending mid-side nodes from the 8-node element, its shape functions being corrected in the standard way so that the affected edge becomes linear); tie the free mid-side degrees of freedom to their corner neighbours with multi-point constraint equations, $u_{\text{mid}}=\tfrac{1}{2}(u_{a}+u_{b})$; or, for the offset in defect C, simply re-mesh so that every element edge is shared in its entirety by exactly one neighbour.
A square bilinear (Q4) element has four nodes and two degrees of freedom per node, hence an $8\times8$ stiffness matrix. Integrated exactly, with the $2\times2$ Gauss rule, its rank is 5 and it therefore has
$$\boxed{\,3\ \text{zero eigenvalues}\,}$$The count is the number of rigid-body modes available to a plane body: two translations ($u$ = const, $v$ = const) and one in-plane rotation about the normal. Their physical significance is that a rigid-body motion produces no strain and hence no strain energy, so $\mathbf{k}\mathbf{d}=\mathbf{0}$ for each of them. This number is a diagnostic that every new element must pass: exactly three in two dimensions, six in three dimensions. Fewer means the element cannot represent rigid-body motion, it will show self-straining and it will not converge; more means the element has a spurious mechanism.
Two riders complete the answer. First, the plane-strain condition does not change the count. Plane strain alters only the constitutive matrix $\mathbf{D}$; it does not touch the kinematics, so the null space of $\mathbf{B}$, and with it the number of zero eigenvalues, is the same as in plane stress. Second, the count does change with the integration rule. Integrated at one Gauss point at the centroid, the same element has five zero eigenvalues: the three rigid-body modes plus two spurious hourglass modes, the deformation patterns proportional to $\xi\eta$ whose strains happen to vanish at the centroid. Those two are mechanisms, not rigid-body motions; if the mesh cannot restrain them the solution is meaningless, which is why one-point elements are always shipped with hourglass control.
The consistent mass matrix is formed from the same shape functions as the stiffness matrix, which is where its name comes from:
$$\mathbf{m}_{c}=\int_{V}\rho\,\mathbf{N}^{T}\mathbf{N}\,dV$$It is full within the element, symmetric, positive definite, and after assembly it has the same bandwidth as $\mathbf{K}$. Its advantages follow from being variationally consistent: the eigenvalue problem $(\mathbf{K}-\omega^{2}\mathbf{M})\boldsymbol{\phi}=\mathbf{0}$ is then a proper Rayleigh–Ritz procedure, so every computed natural frequency is an upper bound on the corresponding exact frequency and convergence is monotone from above. It represents distributed inertia correctly, it handles rotational degrees of freedom (a beam element's $\theta$ terms) properly, and it is the more accurate choice for flexural vibration and for the lower modes generally. Its disadvantages are computational: it must be assembled, stored and factorised, and in an explicit time-integration scheme the need to solve $\mathbf{M}\ddot{\mathbf{d}}=\mathbf{R}$ at every step removes the entire reason for using explicit integration.
The diagonal (lumped) mass matrix concentrates the element mass at its nodes, obtained by direct particle lumping, by row-summing the consistent matrix, or by the HRZ scaling that keeps the diagonal proportions of $\mathbf{m}_{c}$ and rescales to the correct total mass. Its advantages are equally computational and equally decisive: it is trivial to form and store, and $\mathbf{M}^{-1}$ is obtained by reciprocating the diagonal, so a central-difference explicit step becomes uncoupled and matrix-free. That single property is why every impact, blast and crash code uses lumped mass. It also tends to give frequencies that err on the low side, so running a model with both matrices brackets the exact frequency from above and below, and for wave-propagation problems the lumped matrix is often the more accurate of the two because its dispersion error partly cancels that of the stiffness matrix. Its disadvantages are that it is not variationally consistent, so no error bound applies; that the lumping of rotational inertia for beam and plate elements is ambiguous and a naive scheme leaves zeros on the diagonal, making $\mathbf{M}$ singular; that the answer depends on which lumping rule was chosen; and that accuracy in flexural vibration and in the higher modes is noticeably worse.
The practical rule that follows: use the consistent matrix for implicit dynamics, for eigenvalue extraction and for response-spectrum work, and the lumped matrix for explicit transient analysis, for very large models where storage governs, and as the cheap first pass that brackets a consistent-mass result.
The governing feature of this problem is that the strengthened beam has three materials and, more importantly, two interfaces: the reinforcement–concrete bond and the adhesive layer between the concrete soffit and the steel plate. Whether the model is worth building at all depends on capturing the second of these, because plate-end debonding, not flexural yielding, usually governs the capacity of a plated beam.
Idealisation. For a rectangular or T-beam the economical model is a two-dimensional plane-stress slice through the web with the flange width smeared into the thickness; where torsion, a skewed support or a wide flange matters, go to three-dimensional solid elements. A layered nonlinear frame element is a legitimate faster alternative if only the load–deflection curve is wanted, but it cannot show debonding.
Materials. Concrete needs a nonlinear model — a smeared-crack or concrete damaged-plasticity formulation with the compressive curve defined by $f'_{c}$, a tensile strength $f_{ct}\approx0.6\lambda\sqrt{f'_{c}}$ in the CSA A23.3 form, and a tension-stiffening branch, without which the model sheds load far too quickly after first cracking. Reinforcing steel and the plate are elasto-plastic with strain hardening; the adhesive is characterised by its shear strength and fracture energy.
Reinforcement. Either discrete bar (truss) elements whose nodes are shared with the concrete mesh, which assumes perfect bond; embedded reinforcement superimposed on the solid mesh, which frees the bar from the concrete mesh lines; or a smeared reinforcement ratio folded into the concrete constitutive matrix, which suits heavily reinforced regions and slabs.
The strengthening plate and its interface. Model the plate with its own plane-stress or solid elements at the soffit. Do not merge its nodes with the concrete unless debonding is known to be prevented: merged nodes impose perfect bond and will over-predict the strengthened capacity, sometimes by a large margin. Instead insert interface (cohesive) elements carrying a traction–separation law calibrated to a bond–slip curve, so that interfacial shear and normal peeling stresses at the plate end are computed and the debonding limit state can actually occur. Match the plate mesh to the concrete mesh node for node, or tie them with constraint equations, and keep the plate element aspect ratio under control — the plate is thin and it is easy to produce badly shaped elements.
Solution strategy. The analysis is materially nonlinear and should be run incrementally with Newton–Raphson iteration, using displacement or arc-length control so that the softening branch past the peak load can be traced. Staging matters: a real retrofit is applied to a beam that is already carrying its service load, so apply that load to the unstrengthened model first, activate the plate elements, and only then continue to failure. Skipping the staging credits the plate with strain it never experiences and inflates the answer.
Verification and reporting. Run a mesh-convergence study, check the cracked-section capacity against a hand strain-compatibility calculation to CSA A23.3, check the plate anchorage length separately, and compare against test data where available. Report the load–deflection curve, the crack pattern, the plate stress profile and the interfacial shear at the plate end.
Euler–Bernoulli. The kinematic assumption is that plane sections remain plane and normal to the deformed axis. Transverse shear strain is therefore zero and the section rotation is not independent: $\theta = dv/dx$. Only one field, $v(x)$, is interpolated, and because the strain energy contains $v''$ the interpolation must be $C^{1}$ continuous — hence the Hermite cubics, with a deflection and a rotation as the two degrees of freedom at each node. Since the cubics solve $EI\,v''''=0$ exactly, the element is exact for concentrated loads and gives nodally exact displacements for any distributed load once work-equivalent nodal loads are used. Problem 2 of this paper relies on that property.
Timoshenko. Plane sections remain plane but no longer normal, so a constant shear rotation is admitted:
$$\gamma=\frac{dv}{dx}-\theta,\qquad V=\kappa GA\,\gamma$$with $\kappa$ the shear correction factor (5/6 for a rectangle, $A_{w}/A$ for a wide-flange section). Now $v$ and $\theta$ are independent fields, only $C^{0}$ continuity is required, and the simplest element interpolates both linearly with one deflection and one rotation per node. In dynamics the theory also adds rotary inertia, which changes the dispersion relation at high frequency.
The locking trap. With equal-order interpolation and full integration the linear Timoshenko element cannot represent pure bending: $dv/dx$ is constant while $\theta$ varies linearly, so a spurious shear strain appears and the element becomes far too stiff, with the error growing as $(L/h)^{2}$. The cures are the ones named in question 1.5 — one-point (reduced or selective) integration of the shear term, assumed-strain or MITC formulations, or linked interpolation.
When to use each. Use Euler–Bernoulli for slender members, in practice $L/h\gtrsim10$, which covers ordinary building frames, bridge girders in global analysis, and any case in which the shear strain energy is a negligible fraction of the total. Use Timoshenko when $L/h\lesssim10$ — deep beams, coupling beams in shear walls, corbels, short spans; when the section has a low shear stiffness relative to its flexural stiffness, as in sandwich, laminated composite and open thin-walled sections; in wave-propagation or high-frequency dynamic analysis where shear and rotary inertia change the answer; and whenever the beam must connect to a two- or three-dimensional mesh, because the $C^{0}$ Timoshenko element couples to solid elements far more naturally than a $C^{1}$ Hermite element does.
