Finite Element Analysis of Civil Structures

Introduction and context

Finite Element Analysis (FEA) is a numerical framework for predicting the response of structures whose geometry, loading, material behaviour, or boundary conditions make closed-form analysis impractical. Civil engineering applications include reinforced-concrete buildings, steel frames, bridges, shells, foundations, retaining structures, dams, tunnels, pavements, and soil–structure systems.

The central idea is to replace a continuous structure by a finite collection of interconnected subdomains called elements. Within each element, the unknown field—usually displacement—is approximated by interpolation functions. The governing differential equations are then converted into a finite system of algebraic equations. The solution provides nodal displacements, from which strains, stresses, reactions, and other response quantities are calculated.

FEA is not a substitute for structural reasoning. It is a model-based approximation whose reliability depends on idealisation, mesh quality, material laws, loading, boundary conditions, numerical conditioning, and engineering interpretation. A sophisticated analysis of an inappropriate model can be less reliable than a simple analysis based on sound mechanics.

Learning objectives

After completing this chapter, the learner should be able to:

1. Explain the displacement-based finite element concept and its relation to continuum mechanics.
2. Derive the basic element equations using the principle of virtual work.
3. Distinguish truss, beam, frame, plane-stress, plane-strain, axisymmetric, solid, shell, and plate elements.
4. Assemble and solve the global equilibrium equations while applying restraints and loads correctly.
5. Evaluate mesh convergence, numerical stability, stresses, reactions, and modelling assumptions.
6. Apply FEA to a civil-structure problem and critically assess the results.
7. Recognise nonlinear, dynamic, contact, uncertainty, and verification issues in professional practice.

1. Fundamental formulation

1.1 Continuum equilibrium

For a deformable body occupying volume Ω, static equilibrium in the absence of inertia is expressed by

where σ is the Cauchy stress tensor in force per unit area, b is the body-force vector per unit volume, and ∇·σ denotes the divergence of stress. On a traction boundary Γ_t,

where n is the outward unit normal and \bar{t} is prescribed traction. On a displacement boundary Γ_u, displacement is prescribed as u = \bar{u}. The boundary is partitioned so that Γ_u and Γ_t do not overlap and together form the complete boundary.

For small deformation, the strain–displacement relationship is

where ε is the infinitesimal strain tensor and u is the displacement vector. For a linear elastic material,

where D is the constitutive matrix. These equations assume small strains, small rotations, linear material behaviour, and a known, deformation-independent geometry. They are not adequate by themselves for yielding, cracking, large rotations, or changing contact conditions.

1.2 Weak form and virtual work

The finite element method commonly uses the weak form. Let δu be an admissible virtual displacement that is zero on Γ_u. Multiplication of equilibrium by δu and integration over Ω gives the principle of virtual work:

Here δε is virtual strain and P represents concentrated nodal or boundary loads. The left side is internal virtual work; the right side is external virtual work. The weak form reduces the differentiability requirements on the displacement field and is therefore suitable for piecewise interpolation.

1.3 Displacement approximation

Within element e, the displacement field is approximated by

where N is the shape-function matrix and dₑ is the vector of element nodal degrees of freedom. Consequently,

with B containing the spatial derivatives of shape functions in the appropriate strain arrangement. Substitution into virtual work gives the element stiffness equation

and

where kₑ is the element stiffness matrix and fₑ contains consistent equivalent nodal loads, body-force loads, surface loads, and any element-level prescribed effects. For linear elasticity, kₑ is symmetric. It is positive definite when rigid-body modes have been removed and the material stiffness is physically positive.

The dimensions are important: D has units of force/area, B has units of inverse length, and integration over volume produces kₑ with units of force/length for translational degrees of freedom.

2. Element families used in civil structures

Element selection should follow the structural idealisation, not software convenience.

A truss element has axial stiffness. For a two-node bar of length L and area A, aligned with its local axis,

where E is Young’s modulus. In a two-dimensional global system, the local matrix is transformed using direction cosines c = cos θ and s = sin θ. The resulting global matrix is

This expression applies only to a linear axial member with small displacement and pin-jointed force transfer.

For Euler–Bernoulli beam theory, plane sections remain plane and normal to the neutral axis; transverse shear deformation is neglected. This is often satisfactory for slender beams. Timoshenko beam theory includes shear deformation and is preferable for deep beams, short members, sandwich sections, and some dynamic problems. Shear locking can occur when an inappropriate low-order formulation is used for a slender beam or thin shell.

For isotropic plane stress, with engineering shear strain γ_xy, the constitutive relation is

where E is Young’s modulus and ν is Poisson’s ratio. Plane stress is appropriate for thin regions with free surfaces normal to the thickness; it is not generally appropriate for thick three-dimensional restraint.

3. Assembly of the global system

Element equations are mapped into a global degree-of-freedom numbering system. The assembled linear system is

where K is the global stiffness matrix, d is the global displacement vector, and F is the global load vector. Assembly is an additive operation: contributions from all elements sharing a node accumulate in the corresponding rows and columns. This automatically enforces displacement compatibility at shared nodes.

The matrix K is commonly sparse because each element connects only nearby degrees of freedom. Sparse storage and factorisation are therefore essential for large models. In a linear static analysis, direct solvers such as Cholesky factorisation are effective when K is symmetric positive definite. Iterative solvers may be advantageous for very large systems but require suitable preconditioning.

Boundary conditions

Essential boundary conditions prescribe displacement or rotation. Natural boundary conditions prescribe traction, force, or moment. A restrained degree of freedom must be treated explicitly, usually by elimination, partitioning, penalty methods, or Lagrange multipliers. For example, after partitioning unknown and prescribed degrees of freedom,

If dᵦ is prescribed, the unknown part satisfies

Reactions are recovered from the restrained rows of the original system. A model with insufficient restraints has rigid-body modes and a singular stiffness matrix. A model with excessive or contradictory restraints may be artificially stiff and produce incorrect reactions.

Load representation

A distributed load should normally be converted to consistent nodal loads by

This preserves virtual work. Applying the entire load to one arbitrary node can create unrealistic local response, especially in a coarse mesh. Point loads applied directly to nodes are legitimate when the physical load is truly transferred at that node; otherwise load distribution or a loading plate/contact region is more realistic.

4. Numerical integration and isoparametric elements

Most practical elements use natural coordinates and isoparametric interpolation: the same shape functions describe geometry and displacement. The element integral is transformed to a parent domain using the Jacobian J,

and evaluated numerically, commonly by Gauss quadrature. The number and location of integration points affect accuracy, computational cost, locking, and hourglass modes. Reduced integration can alleviate locking but may introduce zero-energy modes; stabilization or an appropriate element formulation is then required.

The Jacobian must have a physically acceptable determinant throughout the element. A negative or nearly zero determinant indicates inverted or severely distorted geometry. Mesh quality should be checked before interpreting stresses.

5. Worked example: two-dimensional truss

Consider a two-member pin-jointed truss. Node 1 is fixed; node 2 is at coordinates (3 m, 0); node 3 is at (3 m, 4 m). Members 1–3 and 2–3 have area A = 400 mm² and modulus E = 200 GPa. A vertical downward load of 20 kN acts at node 3. The objective is to illustrate stiffness assembly and interpretation, not to replace a complete design check.

Use N and mm. Thus E = 200,000 N/mm², A = 400 mm², and the geometry is 3000–4000–5000 mm. Member 1–3 has c = 3/5 and s = 4/5; member 2–3 is vertical, with c = 0 and s = 1. The axial stiffness factor is

for the diagonal. For the vertical member, AE/L = 200000 × 400/4000 = 20000 N/mm. Construct each four-degree-of-freedom global matrix using the direction-cosine expression, place it into the global matrix, restrain the two translations at node 1 and the horizontal translation at node 2, and solve for the two unknown translations at node 3.

The vertical equilibrium and geometry lead to the same result more transparently. The diagonal member carries axial force N_d and the vertical member carries N_v. Horizontal equilibrium at node 3 requires N_d(3/5) = 0, so N_d = 0 because there is no horizontal external load and the vertical member is vertical. Vertical equilibrium then gives N_v = 20 kN in tension. This illustrates a critical modelling insight: a member may have zero force under a particular load case even though its stiffness contributes to the global system for other load cases.

The vertical elongation of member 2–3 is

The result assumes pin joints, linear elasticity, no self-weight, small displacement, and ideal axial load transfer. A real connection may introduce eccentricity, bending, slip, or local yielding; those effects require a refined model.

6. Verification, validation, and result interpretation

Verification asks whether the equations have been solved correctly. Validation asks whether the mathematical model represents the physical structure adequately. Both are necessary.

Mesh convergence

A mesh-refinement study compares a response quantity as element size decreases. For a reliable quantity Q,

where h is a characteristic element size. Global displacements and total reactions often converge smoothly. Point stresses near a sharp re-entrant corner, point load, crack tip, or idealised support may not converge to a finite value because of a mathematical singularity. In such regions, use stress resultants, averaged stresses away from the singularity, fracture parameters, or a physically detailed connection model.

A practical study should refine the mesh in regions of high gradient, compare displacement and reaction results, inspect energy error indicators where available, and document the convergence criterion. Uniform refinement is not always efficient; adaptive refinement is preferable when error estimates are reliable.

Equilibrium and compatibility checks

The sum of reactions must balance applied loads within numerical tolerance. Element deformation must be compatible at shared nodes. The strain energy for a linear static model is

when the only loads are conservative and prescribed displacement effects are treated consistently. This identity is a useful audit check.

Stress recovery

Displacements are primary finite element variables in standard displacement formulations. Strains are obtained by differentiation and stresses from the constitutive law. Therefore, stresses are generally less accurate than displacements and may vary significantly between integration points and nodes. Nodal stress averaging can make contours visually smooth while concealing discontinuities; averaged values should not automatically be treated as more accurate.

For reinforced concrete, elastic principal stresses from an uncracked continuum model are not directly equivalent to design actions after cracking. Cracking, tension stiffening, reinforcement orientation, bond, creep, shrinkage, and construction sequence may need explicit treatment or separate design procedures.

7. Nonlinear and dynamic analysis

Material nonlinearity

When the constitutive relation depends on loading history,

the tangent stiffness changes during the solution. Plasticity, concrete cracking and crushing, steel yielding, viscoelasticity, and damage are examples. Incremental–iterative methods such as Newton–Raphson solve the residual equation

using

where K_t is the tangent stiffness and Δd is an incremental correction. Convergence must be assessed using force residuals, displacement corrections, and energy measures—not merely by reaching a prescribed iteration count.

Geometric nonlinearity

Large displacement or rotation changes the equilibrium configuration and may cause P–Δ effects, buckling, or post-buckling response. The tangent stiffness can be written conceptually as

where the geometric component depends on current stress and configuration. Initial imperfections and load-path control can strongly influence the computed critical and post-critical response.

Contact and interfaces

Contact introduces inequalities: bodies may transmit compression while separating in tension, and friction limits tangential traction. Penalty, augmented-Lagrangian, and constraint-based formulations are common. Contact results are sensitive to interface stiffness, friction coefficient, mesh compatibility, and contact-search tolerances. A very large penalty stiffness can cause ill-conditioning; a small value permits excessive penetration.

Dynamic analysis

For structural dynamics,

where M is the mass matrix, C is damping, K is stiffness, d is displacement, and F(t) is time-dependent loading. A consistent mass matrix follows the same interpolation principle as stiffness; a lumped mass matrix places mass at selected degrees of freedom. Modal analysis solves

where ω is circular natural frequency and φ is a mode shape. Frequencies are in radians per second; the cyclic frequency is f = ω/(2π). Damping assumptions should be stated because damping is not uniquely determined by static material properties.

8. Applied case: finite element assessment of a reinforced-concrete transfer floor

The following is an illustrative case, not a project-specific design.

Background and problem

A twelve-storey reinforced-concrete building contains a transfer floor where upper-storey columns do not align with lower-storey supports. The design team must assess slab deflection, load redistribution, punching-shear demand, and sensitivity to column stiffness. A simple beam model cannot represent two-way slab action and local stress concentration adequately.

Analysis approach

The slab and transfer beams are represented by shell elements, while columns are represented by beam-column elements with appropriate section properties. The model includes self-weight, superimposed dead load, factored gravity combinations, and selected lateral load combinations. Column–slab interfaces are checked for realistic connectivity rather than assuming all nodes share identical rotations. A local mesh refinement is introduced around transfer columns, openings, and beam intersections.

Material properties are assigned consistently with the design code and documented separately for service and ultimate combinations. The model is first solved elastically to check reactions, equilibrium, symmetry where applicable, and deflection plausibility. A second model includes cracked-section stiffness assumptions for service analysis. Sensitivity runs vary effective slab stiffness, column stiffness, and support idealisation.

Findings and interpretation

The coarse model predicts acceptable global deflection but produces a highly mesh-dependent peak moment at a column corner. Refinement changes the peak value substantially while the integrated support reaction and average strip moment stabilise. The correct engineering response is not to accept the largest contour value uncritically. Instead, design actions are extracted using rational column strips, support regions, or code-prescribed averaging, while local punching and anchorage are checked with an appropriate detailed model.

Sensitivity analysis shows that reducing cracked slab stiffness increases transfer-floor deflection and redistributes load toward stiffer beams. This result demonstrates that stiffness assumptions can be as influential as mesh density. The final design decision includes additional top reinforcement over supports, local shear reinforcement where required, construction-stage review, and monitoring of deflection during early occupancy.

Learning points

1. Global equilibrium is necessary but does not guarantee local accuracy.
2. A shell model requires careful interpretation of drilling degrees of freedom, thickness, transverse shear, and reinforcement representation.
3. Mesh convergence should be assessed for design actions, not only for colourful stress contours.
4. Sensitivity studies expose uncertainty in cracked stiffness and connectivity.
5. FEA output must be converted into code-compatible design quantities by a qualified engineer.

9. Professional practice, governance, and emerging developments

A finite element model should have a documented purpose, geometry, units, material sources, element types, mesh strategy, boundary conditions, load combinations, solver settings, convergence criteria, and post-processing rules. Version control is important because changing software versions, element libraries, or default settings can alter results.

Engineers must distinguish computed precision from physical certainty. Reporting six significant figures for a displacement based on uncertain stiffness, cracking, or support conditions is misleading. Independent hand calculations, benchmark problems, alternative element formulations, and peer review provide important safeguards.

AI-assisted meshing, surrogate models, automated parameter studies, and digital twins are increasingly used. They can improve efficiency but do not remove the need for mechanics-based verification. Training-data bias, opaque model decisions, data security, software traceability, and responsibility for safety-critical decisions must be addressed. Automated optimisation must include constructability, robustness, fire, durability, fatigue, inspection, and code requirements rather than optimising a single numerical objective.

Common failure modes include unit inconsistency, unintended releases, double-counted self-weight, inadequate restraints, unrealistic rigid links, poor aspect ratios, excessive mesh distortion, singular point loads, inappropriate linearisation, ignoring construction stages, and interpreting nodal peaks as physical stresses. A disciplined pre-processing checklist and independent post-processing review are among the most effective controls.

Chapter Summary

Finite Element Analysis converts the continuum equations of structural mechanics into a finite set of algebraic equations by subdividing a structure into elements and approximating its displacement field with shape functions. The principle of virtual work leads to the element stiffness matrix kₑ = ∫ BᵀDB dΩ, and assembly produces the global system Kd = F. Correct element selection, coordinate transformation, load representation, restraint treatment, numerical integration, and unit consistency are essential.

FEA results must be verified through equilibrium, compatibility, energy checks, solver diagnostics, and mesh convergence. Validation requires judgement about material behaviour, geometry, connections, construction stages, cracking, contact, and boundary conditions. Nonlinear and dynamic analyses extend the basic formulation but introduce path dependence, tangent stiffness, convergence, mass, damping, and stability issues. Professional use requires transparent documentation, independent checks, sensitivity studies, and responsible interpretation of uncertainty.

Key Terms

Finite element: A subdomain used to discretise a continuous structure.

Node: A point at which degrees of freedom and connectivity are defined.

Degree of freedom: An independent displacement or rotation variable.

Shape function: An interpolation function used to approximate a field within an element.

Element stiffness matrix: Matrix relating element nodal actions to element nodal displacements.

Global stiffness matrix: Assembled matrix representing the stiffness of the complete discretised structure.

Constitutive matrix: Matrix relating strains to stresses for a specified material model.

Virtual work: Work equation using admissible virtual displacements.

Essential boundary condition: Prescribed displacement or rotation.

Natural boundary condition: Prescribed force, moment, or traction.

Isoparametric element: Element using common interpolation functions for geometry and field variables.

Jacobian: Transformation matrix between natural and physical coordinates.

Mesh convergence: Stabilisation of a response quantity as the mesh is refined.

Locking: Artificial stiffness caused by an inadequate finite element formulation.

Rigid-body mode: A deformation mode involving motion without strain energy.

Tangent stiffness: Local stiffness matrix used to linearise a nonlinear problem.

Geometric nonlinearity: Nonlinearity caused by changes in configuration or equilibrium geometry.

Modal analysis: Determination of natural frequencies and mode shapes.

Verification: Demonstration that the numerical equations are solved correctly.

Validation: Assessment that the model represents the physical system adequately.

Review Questions

Short Answer Questions

1. What is the principal approximation made in the displacement-based finite element method?
2. Distinguish essential and natural boundary conditions.
3. Why can a stiffness matrix become singular?
4. State two situations in which plane stress is appropriate.
5. Why are stresses generally less reliable than displacements in a standard displacement formulation?

Long Answer Questions

1. Derive the linear elastic finite element equilibrium equation from the principle of virtual work and define all matrices involved.
2. Compare truss, beam, shell, and solid elements with reference to civil engineering applications and assumptions.
3. Explain the assembly and application of boundary conditions in the global equation Kd = F.
4. Discuss mesh convergence, stress singularities, and the interpretation of nodal stress contours.
5. Explain how material, geometric, and contact nonlinearities alter the solution procedure.

Analytical / Application Questions

1. Develop the global stiffness matrix for a two-dimensional two-member truss and identify the boundary conditions required to remove rigid-body modes.
2. Propose a finite element modelling strategy for a transfer slab with columns, openings, and concentrated loads. Justify element selection, mesh refinement, load application, and result extraction.
3. A nonlinear analysis fails to converge near the predicted peak load. Formulate a diagnostic plan considering load increments, tangent stiffness, material parameters, contact, mesh quality, and physical instability.

Multiple Choice Questions with Answers

Question 1

Which equation represents the assembled linear static finite element system?

A. Kd = F

B. K + d = F

C. K/F = d

D. K = Fd

Correct Answer: A. Kd = F

Question 2

For a two-dimensional truss member, which quantities determine its axial stiffness factor?

A. Poisson’s ratio and thickness only

B. Area, Young’s modulus, and length

C. Density, damping, and frequency

D. Bending moment and shear force

Correct Answer: B. Area, Young’s modulus, and length

Question 3

Which boundary condition prescribes displacement or rotation?

A. Natural boundary condition

B. Traction condition

C. Essential boundary condition

D. Equilibrium condition

Correct Answer: C. Essential boundary condition

Question 4

In a standard displacement formulation, strains are obtained primarily by

A. Integrating nodal forces

B. Differentiating the interpolated displacement field

C. Averaging support reactions

D. Solving the mass equation

Correct Answer: B. Differentiating the interpolated displacement field

Question 5

A model with insufficient restraints commonly produces

A. A singular stiffness matrix due to rigid-body modes

B. Automatically converged stresses

C. Zero material density

D. Exact local stress recovery

Correct Answer: A. A singular stiffness matrix due to rigid-body modes

Question 6

For a thin plate whose normal stresses through the thickness are negligible, the most suitable idealisation is generally

A. Plane stress

B. Plane strain

C. Axisymmetric solid only

D. One-dimensional diffusion

Correct Answer: A. Plane stress

Question 7

The purpose of a mesh-convergence study is to determine whether

A. The software licence is active

B. A response stabilises as element size is reduced

C. All stresses become exactly zero

D. The structure has no loads

Correct Answer: B. A response stabilises as element size is reduced

Question 8

In nonlinear analysis, the tangent stiffness is used to

A. Linearise the current equilibrium problem for an incremental correction

B. Replace all material properties with density

C. Eliminate the need for boundary conditions

D. Calculate only the undeformed geometry

Correct Answer: A. Linearise the current equilibrium problem for an incremental correction

Question 9

Which equation is the standard undamped-free modal eigenproblem when mass and stiffness matrices are used?

A. (K - ω²M)φ = 0

B. (K + ωM)φ = F

C. Mφ = K

D. ω = K/M

Correct Answer: A. (K - ω²M)φ = 0

Question 10

Why should a peak stress at an idealised point load or sharp re-entrant corner be interpreted cautiously?

A. It may be mesh-dependent or singular

B. It is always equal to the average stress

C. It cannot be affected by geometry

D. It represents only rigid-body motion

Correct Answer: A. It may be mesh-dependent or singular