Nonlinear Finite Element Analysis of Concrete Strip Foundations on Elastic Support

Introduction

Strip foundations are fundamental structural elements that transfer building loads to the underlying soil. Traditional analysis methods assume linear elastic behavior, which fails to capture the actual response of concrete foundations that undergo cracking, creep, and nonlinear soil interaction. This article describes a computational framework that addresses these limitations through:

1. Beam-on-elastic-foundation theory with Winkler spring model
2. Nonlinear material modeling via moment-curvature relationships
3. Time-dependent effects through creep and aging
4. Advanced reinforcement including steel bars and structural fibers
5. Iterative FEA solution with adaptive mesh refinement

Theoretical Foundation

Governing Differential Equation

A beam on elastic foundation is governed by the fourth-order differential equation:

where:

• $w(x)$ = vertical deflection [m]
• $EI$ = flexural rigidity [MN·m²]
• $k(x)$ = modulus of subgrade reaction [kN/m/m]
• $q(x)$ = applied load per unit length [kN/m]

For nonlinear analysis, the stiffness $EI$ becomes a function of curvature:

where $\kappa$ is the curvature [m⁻¹] and $M$ is the corresponding moment [kNm].

Characteristic Length

The characteristic length defines the zone of influence for concentrated loads and governs mesh refinement:

where $W$ is the foundation width [m]. This parameter determines:

• The distance over which load effects propagate
• Required analysis length for bounded solutions
• Optimal element sizing for FEA mesh

Modulus of Subgrade Reaction

The soil stiffness is characterized by the modulus of subgrade reaction, modified by Terzaghi's width correction:

where $B$ is the foundation width [m] and $k_0$ is the reference modulus at $B = 0.3$ m.

For layered soil profiles, the code implements depth-dependent stiffness based on Janbu's tangent modulus:

where $m$ and $j$ are soil parameters, $\sigma_v$ is vertical stress, and $z$ is depth.

Material Modeling

Concrete Constitutive Behavior

The concrete stress-strain relationship follows Eurocode 2 with parabolic-rectangular distribution:

where:

• $f_{cd} = \alpha_{cc} f_{ck} / \gamma_c$ = design compressive strength
• $\varepsilon_{c2}$ = strain at peak stress
• $n$ = parabola exponent (typically 2)

The effective modulus accounts for creep:

where $\varphi(\infty, t_0)$ is the creep coefficient from EN 1992-1-1.

Creep Modeling

Long-term deflections require accounting for creep. The code implements the age-adjusted effective modulus method:

where:

• $\varphi_0$ = notional creep coefficient depending on humidity, member size, and concrete strength
• $\beta_c(t - t_0)$ = time development function
• $t_0$ = age at loading
• $t$ = age at analysis

The notional size is:

where $A_c$ is the cross-sectional area and $u$ is the perimeter exposed to drying.

Reinforcement

Steel Reinforcement

Bilinear stress-strain with or without strain hardening:

where $f_{yd} = f_{yk}/\gamma_s$ and $\varepsilon_{yd} = f_{yd}/E_s$.

Fiber Reinforcement

The code implements the 𝜎-𝜀 design method from EN 1992-1-1:2023 Annex L.5.5.2 with residual strength parameters

Post-Cracking Stress-Strain Relationship

The bilinear constitutive law for structural analysis (EN 1992-1-1 L.5.5.2) is:

Design Residual Strengths

The effective residual tensile strengths are calculated from characteristic values:

Where:

• $f_{R1k}$, $f_{R3k}$ = Characteristic residual flexural tensile strengths at CMOD = 0.5 mm and 2.5 mm
• $k_o$ = Orientation factor (accounts for fiber alignment, typically 0.5 to 2.0)
• $k_G$ = Global safety factor for static indeterminacy

Strain Limits

Where $s_{rm,cal,F}$ is the mean crack spacing (see EN 1992-1-1 L.9.3, Formula L.26). As a simplification, a constant value of $l_{cs} = 125$ mm may be used unless explicitly precluded by the National Annex (EN 1992-1-1 L.5.5.2 Note 1).

Moment-Curvature Analysis

For each section, the moment-curvature relationship is computed by:

1. Assuming a curvature $\kappa$
2. Determining neutral axis position $x_{\text{NA}}$ such that:

where $\varepsilon(y) = \kappa \cdot (y - x_{\text{NA}})$

3. Computing moment:

4. Repeating for range of curvatures to build $M(\kappa)$ curve

The resulting relationship provides the variable stiffness function:

Finite Element Formulation

Element Stiffness Matrix

The foundation is discretized into Euler-Bernoulli beam elements. For a single element with length $L_e$ and constant $EI$:

For vertical deflection-only analysis (single DOF per node), the reduced stiffness is:

Spring Elements

The Winkler foundation is modeled as discrete springs at each node:

where $L_{\text{trib},i}$ is the tributary length for node $i$:

The spring stiffness is added directly to the diagonal of the global stiffness matrix.

Load Vector

The nodal forces account for distributed and concentrated loads. For a uniformly distributed load $q$ over element length $L_e$:

For concentrated loads, the code distributes the load over its specified length using tributary areas:

where $w_j(x_i)$ is the weight function for load $j$ at node $i$.

Global System

Assembling all elements yields the global system:

where $[K_{\text{global}}]$ includes both beam and spring contributions:

Nonlinear Solution Algorithm

Iterative Procedure

For nonlinear analysis where $EI = EI(\kappa)$, the code implements a fixed-point iteration:

The element update function computes the average curvature at element midpoint:

Spring Update for Tension

Soil cannot sustain tension. When node displacement is positive (upward), the spring stiffness is reduced:

This creates an additional source of nonlinearity handled by the iteration.

Convergence Criteria

The iteration converges when:

If convergence fails, the utilization ratio is set to 9.99 to flag the issue.

Adaptive Mesh Refinement

Mesh Generation Strategy

The code implements adaptive mesh refinement around load positions:

Refinement zones:

• Fine mesh: $|x - x_{\text{load}}| < 0.3 l_c$, element size $\Delta x = 0.3h$
• Transition: $0.3 l_c < |x - x_{\text{load}}| < l_c$, linearly interpolated
• Coarse mesh: $|x - x_{\text{load}}| > l_c$, element size $\Delta x = 2h$

For multiple loads, the minimum distance to any load determines the local element size:

Analysis Length Determination

The analysis length must be sufficient to ensure boundary effects do not influence the region of interest. For a semi-infinite beam on elastic foundation, load effects decay exponentially with distance characterized by the characteristic length $l_c$.

For concentrated loads:

For uniform loads:

where $h$ is the foundation height and $L_{\text{span}}$ is the extent of loading. These expressions ensure the analysis domain extends at least 2-3 characteristic lengths beyond load application points, allowing load effects to decay to negligible levels before reaching the free-end boundaries.

If the analysis length is reduced below these recommendations, the free-end boundary conditions will increasingly influence the solution in the loaded region, potentially underestimating moments and deflections for what would otherwise behave as a longer beam.

Section Force Recovery

Moment Calculation

Once displacements are known, moments at element nodes are recovered from curvature:

For linear analysis:

For nonlinear analysis:

Shear Force Calculation

Shear forces are computed from element equilibrium. For a beam element with end moments $M_1, M_2$ and distributed load $q$:

Point springs create discontinuities, so shear is averaged:

Crack Width Estimation

For serviceability limit state, crack widths are estimated from the moment or curvature:

Linear analysis:

Nonlinear analysis:

where $f_{cw}$ is the crack width function derived from the detailed section analysis.

Verification and Utilization Ratios

Ultimate Limit State (ULS)

For ULS verification, utilization ratios are computed:

Moment capacity:

where:

• Linear: $M_{Rd}$ from section capacity calculation
• Nonlinear: $UR_M = \max\left( \frac{\kappa_{Ed}}{\kappa_{Rd,\text{pos}}}, \frac{\kappa_{Ed}}{\kappa_{Rd,\text{neg}}} \right)$

Shear capacity:

where $V_{Rd}$ is calculated per EN 1992-1-1 Section 6.2, accounting for concrete contribution and stirrups.

Serviceability Limit State (SLS)

Crack width:

where $w_{\text{lim}}$ is the allowable crack width.

Settlement:

Implementation Details

Load Application for Multiple Concentrated Loads

When multiple concentrated loads are specified, each load $j$ contributes to node $i$ according to:

where $w_{ij}$ is the tributary length of node $i$ within load $j$'s distribution length.

The weight function accounts for load spread:

Stiffness Variation

For differential settlement analysis (uniform load case), soil stiffness varies:

where $\delta$ is the transition length.

For concentrated load cases, uniform minimum stiffness is used:

Numerical Integration

The section analysis requires numerical integration over the cross-section. The code uses fiber discretization:

where each fiber has area $A_i$ at position $\gamma_i$ with stress $\sigma_i(\varepsilon_i)$.

Validation and Limitations

Validation Cases

The implementation has been validated against:

1. Hetényi's analytical solutions for beams on elastic foundation
2. Finite element programs

Key Assumptions

• Plane sections remain plane (Bernoulli hypothesis)
• Winkler foundation model (no interaction between springs)
• Small deflections (linear geometry)
• No bond slip between concrete and reinforcement
• Uniform soil properties in transverse direction

Limitations

• Large deflections
• Soil-structure interaction beyond Winkler model
• Three-dimensional effects
• Dynamic loading
• Thermal effects

Conclusion

This implementation provides a comprehensive framework for analyzing concrete strip foundations with:

1. Nonlinear material behavior capturing cracking and creep
2. Variable stiffness finite elements for accurate load distribution
3. Adaptive meshing for computational efficiency
4. Multiple load configurations including concentrated and distributed
5. Full code compliance with Eurocode 2 and Eurocode 7

The method enables engineers to:

• Quantify the actual safety margins beyond empirical reinforcement rules
• Optimize reinforcement for specific loading conditions
• Account for time-dependent effects in long-term performance
• Validate or challenge the Danish SBi-guidance practice of using steel reinforcement equivalent to 0.2% of the concrete cross section

References

This theoretical framework forms the basis of the greenStrip analysis engine, implementing state-of-the-art computational methods for foundation design while maintaining computational efficiency suitable for practical engineering applications.