Question 1 of 3: Ten short questions on finite element theory

Problem 1: Ten short questions on finite element theory (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.

Given. No numerical data: the ten parts examine the theory that sits underneath every displacement-based finite element code, from interelement continuity to integration order and mesh refinement.

Find. A short, defensible answer to each of the ten parts, with the two “draw” parts supported by a figure and the stress-field part settled by differentiation rather than by inspection.

Part 1 — a $C^{0}$ field and a $C^{1}$ field, and their first derivatives

The superscript counts how many derivatives survive the element boundary. A $C^{0}$ function is itself continuous but its first derivative is not: across the shared node the value matches, the slope does not. A $C^{1}$ function carries its first derivative through as well; the discontinuity is pushed up to the second derivative. Written as a limit at the shared node $x_{n}$ of two elements,

The distinction is not academic: it fixes which elements may be connected to which. A bar, a plane-stress membrane and a solid element interpolate displacement only, so $C^{0}$ is both necessary and sufficient — the strains they need are first derivatives, and a finite jump in a first derivative is square-integrable. A thin beam or plate element needs curvature, a second derivative, so its trial field must be $C^{1}$; that is why the plane beam element carries a rotation as well as a translation at each node. Elements that fail the requirement are called non-conforming and are not guaranteed to converge, although several deliberately non-conforming plate elements do pass the patch test in practice.

Part 2 — why $\sum N_{i}=1$ fails for the plane beam element

The partition-of-unity property is a completeness requirement in disguise. For a $C^{0}$ element every degree of freedom is a sample of the same field, $u=\sum N_{i}u_{i}$, so representing an arbitrary rigid-body value $u\equiv c$ demands $\sum N_{i}c=c$ for every $c$, which is exactly $\sum N_{i}=1$.

The plane beam element interpolates one field from two different kinds of degree of freedom. With $\xi=x/L$ the Hermite functions are

$N_{2}$ and $N_{4}$ multiply rotations, so they carry the dimension of length; adding them to the dimensionless $N_{1}$ and $N_{3}$ is not even a legitimate operation. Completeness instead splits into one condition per rigid-body mode. A rigid translation $v\equiv c$ has $v_{1}=v_{2}=c$ and $\theta_{1}=\theta_{2}=0$, and a rigid rotation $v=\theta_{0}x$ has $v_{1}=0$, $v_{2}=\theta_{0}L$, $\theta_{1}=\theta_{2}=\theta_{0}$, so the two requirements are

Both are satisfied identically by the functions above, whereas the naive sum evaluates to $1+L(\xi-3\xi^{2}+2\xi^{3})$, which equals unity only at the two nodes. The correct general statement is therefore that the shape functions belonging to a single field variable form a partition of unity, and that the remaining functions must reproduce the higher rigid-body modes instead.

Part 3 — what the compatibility condition means

In two dimensions three strain components are built from only two displacement components, so the strains cannot be independent: there must be one differential relation among them. Substituting $\varepsilon_{x}=u_{,x}$, $\varepsilon_{y}=v_{,y}$ and $\gamma_{xy}=u_{,y}+v_{,x}$ into the printed equation gives

an identity, which is the point: the condition is not extra physics but the statement that a proposed strain field is integrable. Physically, compatibility guarantees that the deformed body still fits together — no gaps open up, no material overlaps itself, and the displacement field recovered by integrating the strains is single-valued and continuous. Violate it and the strain field describes a body that has been cut into pieces which no longer assemble.

Two consequences matter for finite element work. First, in a displacement-based formulation compatibility is satisfied automatically and for free, because the strains are differentiated from an assumed continuous displacement field; that is precisely why the method never needs to enforce it. Second, in a force or stress-based (flexibility) formulation, and in any analytical solution posed in terms of stresses, compatibility must be imposed explicitly — which is what Part 4 tests.

Part 4 — is the given stress field a valid plane-elasticity solution?

A candidate stress field has two hurdles to clear: pointwise equilibrium, and compatibility expressed in terms of stress. Take equilibrium first, with zero body force.

Both equations are satisfied identically, so the field is statically admissible and a candidate would be tempted to stop there. Compatibility, however, reduces for a plane problem with zero body force to the Lévy condition $\nabla^{2}(\sigma_{x}+\sigma_{y})=0$, obtained by writing the strains in terms of stress and substituting into the relation of Part 3. Here the first stress invariant is $\sigma_{x}+\sigma_{y}=3a_{1}x^{2}y+a_{1}y^{3}$, so

The field is therefore not a valid solution unless $a_{1}=0$ or attention is restricted to the line $y=0$. Running the plane-stress version of Part 3 explicitly gives the same verdict from a different direction: $\varepsilon_{x,yy}+\varepsilon_{y,xx}= -12\nu a_{1}y/E$ while $\gamma_{xy,xy}=-12(1+\nu)a_{1}y/E$, and the two agree only if $a_{1}y=0$. Equivalently, the field derives from the Airy function $\phi=\tfrac{1}{4}a_{1}y^{4}+\tfrac{1}{2}a_{1}x^{2}y^{2}$ — which reproduces the stresses and hence guarantees equilibrium, but is not biharmonic. Equilibrium alone is never sufficient in elasticity.

Part 5 — shape functions of the two bar elements

A bar carries one degree of freedom per node, the axial displacement, so each shape function is the Lagrange polynomial that takes the value one at its own node and zero at every other node. Figure 1-a has two nodes, giving linear functions; Figure 1-b has three, giving quadratics — and because its interior node sits at $L/3$ rather than at midspan, they are not the symmetric textbook quadratics.

Both sets satisfy $\sum N_{i}=1$ everywhere, which is what lets either element carry a rigid-body translation. The off-centre node makes the second set markedly asymmetric: $N_{1}$ dips to $-1/3$ at $x=2L/3$, $N_{2}$ dips to $-1/24$ at $x=L/6$, and $N_{3}$ peaks at $1.125$ at midspan rather than reaching its maximum of one at its own node. Overshoot of that size is the practical argument for keeping mid-side nodes near the middle: it inflates the condition number of the element and, if the node is pushed past $L/4$, the Jacobian of the isoparametric map changes sign inside the element.

The strain follows from differentiating: element 1-a gives $\varepsilon=(u_{2}-u_{1})/L$, a constant, while element 1-b gives a strain that varies linearly along the bar. That is the whole reason for the extra node — one three-node bar reproduces a linearly varying axial force, such as that caused by self-weight, exactly.

Part 6 — the physical meaning of $\det[J]$ in a Q4

The Jacobian matrix maps the parent square $-1\le\xi,\eta\le 1$ onto the physical quadrilateral, and its determinant is the local area magnification of that map:

So $\det[J]$ evaluated at an integration point is the physical area that one unit of parent area occupies there — the weight by which that Gauss point contributes to $[k]=\int[B]^{T}[D][B]\,t\,\det[J]\,d\xi\,d\eta$. For a rectangle $a\times b$ it is the constant $ab/4$; for a distorted quadrilateral it varies from point to point, and the ratio of its largest to its smallest value over the Gauss points is the “Jacobian ratio” that mesh generators report as a distortion measure.

It must be strictly positive throughout the element. A zero value means the map has collapsed — two nodes coincide, or three are collinear — and the inverse $[J]^{-1}$ needed to form $[B]$ does not exist. A negative value means the element has turned inside out at that point, which for a Q4 happens as soon as an interior angle reaches $180^\circ$; the integration then subtracts stiffness and the assembled matrix can lose positive definiteness. Practical codes test the sign of $\det[J]$ at every Gauss point and reject the mesh rather than the analysis.

Part 7 — non-zero eigenvalues of a square bilinear element in plane strain

Count the degrees of freedom, subtract the rigid-body modes. A Q4 has four nodes and two translations each, so $[k]$ is $8\times 8$; a plane body has three rigid-body modes, two translations and one in-plane rotation, and each of them produces zero strain and hence zero strain energy. Therefore

Computing the eigenvalues of a square plane-strain Q4 confirms it exactly: three vanishing values and five positive ones, the largest belonging to the uniform-dilatation mode, which in plane strain is stiffened by the bulk term $\lambda+2\mu$. The five non-zero modes are the two uniform extensions, the uniform shear, and the two in-plane bending (flexural) modes.

The answer is tied to the integration rule, and Part 8 is the reason. Full $2\times 2$ Gauss integration supplies $4\times 3=12$ independent strain sampling conditions, more than the five needed, so the rank is the theoretical maximum of five. One-point reduced integration supplies only three, so the rank falls to three: two extra zero-energy modes appear, the notorious hourglass modes, in which the element deforms with zero strain at its centre. Those are spurious mechanisms, not rigid-body motions, and they must be suppressed by hourglass control.

Part 8 — advantages of reduced integration

Reduced integration means deliberately using a Gauss rule of lower order than the one that would integrate $[B]^{T}[D][B]\det[J]$ exactly — one point instead of $2\times 2$ for a Q4, $2\times 2$ instead of $3\times 3$ for a Q8. Four advantages are worth marks.

Cost. The stiffness matrix and the stress recovery are formed at fewer points, and the saving grows steeply with dimension and order: a 20-node brick drops from 27 evaluations to 8, roughly a threefold reduction in element formation time.

It cures locking. This is the real reason. A fully integrated low-order element imposes more strain constraints than it has degrees of freedom to satisfy, and the surplus constraints make it artificially stiff. Under-integrating removes exactly those surplus constraints. In bending-dominated problems it eliminates shear locking, the parasitic shear strain that makes a single fully integrated Q4 far too stiff in flexure; in nearly incompressible work, $\nu\to 0.5$ or fully developed plastic flow, it eliminates volumetric locking, where the constraint $\varepsilon_{v}\to 0$ at every Gauss point over-determines the mesh. Coarse-mesh displacements improve by an order of magnitude.

Better stresses. The reduced points coincide with the Barlow, or optimal, sampling points, where the finite element stress is superconvergent — typically one order more accurate than at the nodes. Extrapolating from them is how commercial codes report element stresses.

A softer, more forgiving element. Because the model is no longer bounded from below in energy, convergence is often faster from a coarse mesh, which matters in explicit dynamics where one-point Q4 and brick elements are standard.

The price, which a complete answer should state, is rank deficiency: as computed in Part 7 the one-point Q4 has only three non-zero eigenvalues, so two spurious hourglass mechanisms exist. In a patch of elements they may or may not be restrained by the neighbours, so production codes add an artificial hourglass stiffness or use selective (deviatoric-only) reduced integration, or an assumed-strain B-bar formulation, to get the benefit without the mechanism.

Part 9 — h-refinement for convergent stresses

In $h$-refinement the element type and polynomial order $p$ are held fixed and the characteristic element size $h$ is reduced — each element is subdivided, typically $h\to h/2$, and the problem is re-solved. For a smooth solution the error in the energy norm and in the stresses then decays as a power of $h$,

The reason a deliberate refinement study is needed at all is that a displacement-based finite element stress field is discontinuous between elements and does not satisfy the differential equations of equilibrium pointwise. Displacements converge from below and converge quickly; stresses, being derivatives of an approximation, converge one order more slowly and must be demonstrated rather than assumed.

The working procedure is a loop. Solve on a coarse mesh. Form an error indicator — in practice the difference between the raw discontinuous stresses and a smoothed nodal stress field, which is the Zienkiewicz–Zhu superconvergent patch recovery estimator, or simply the traction jump across interelement boundaries. Refine where the indicator is largest, which is adaptive $h$-refinement, or refine uniformly if the field is smooth. Re-solve, and compare the stress at the same physical point on successive meshes; stop when the change falls below the tolerance the design needs. Two or three levels usually suffice, and Richardson extrapolation on the sequence gives an estimate of the converged value.

Part 10 — two techniques for modelling embedded reinforcement

1. Discrete reinforcement. Each bar is modelled explicitly as a one-dimensional bar or beam element lying along its true line, with cross-sectional area $A_{s}$ and the steel constitutive law. If its nodes are shared with the surrounding concrete elements, perfect bond is implied and the bar simply adds axial stiffness along that line. To model slip instead, the bar nodes are duplicated and connected to the concrete nodes by bond-link or interface elements — short springs whose axial and transverse laws come from a pull-out test. The method reproduces individual bar forces, anchorage and dowel action directly, at the cost of forcing the concrete mesh to follow the bar layout.

2. Smeared or embedded reinforcement. Here the steel is not given its own nodes. In the smeared form the reinforcement is treated as a uniformly distributed stiffness inside the concrete element, added to the constitutive matrix in the bar direction in proportion to the reinforcement ratio $\rho=A_{s}/A_{c}$; perfect bond is assumed by construction. In the closely related embedded form a bar segment is allowed to run anywhere through a parent solid element and its degrees of freedom are eliminated by tying them to the parent element shape functions, so the mesh need not conform to the bar. Both are the practical choice for large models, walls and slabs with dense uniform reinforcement.

In Canadian practice the constitutive laws behind either choice come from CSA A23.3, and the tension-stiffening and bond models needed for a smeared analysis are those of the Modified Compression Field Theory that the code embodies.