⚡ electrical engineering

A Physically Constrained Energy-Based Excess Pore Pressure Model for Saturated Sands with an Ordinary Differential Equation Interpretation

This study reformulates the energy-based excess pore pressure model for saturated clean sands by interpreting it as a solution to a Cauchy–Euler ordinary differential equation, thereby ensuring physical boundedness, reducing calibration to a single parameter, and achieving high accuracy across various relative densities.

Original authors: Ferruh Gündoğan, Ayfer Erken

Published 2026-06-30

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Ferruh Gündoğan, Ayfer Erken

Original paper licensed under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

The Big Picture: Why Do Sandcastles Collapse?

Imagine you are at the beach building a sandcastle. If you just stand on it, it holds its shape. But if you start jumping up and down on it rhythmically, the sand grains start to jostle. In a saturated sandcastle (one soaked with water), this shaking forces the water trapped between the grains to push harder and harder against the sand.

Eventually, the water pressure becomes so strong that it pushes the sand grains apart completely. The sand stops acting like a solid and starts acting like a thick liquid (like a milkshake). This is called liquefaction, and it's what causes buildings to sink or tilt during earthquakes.

Engineers need a way to predict exactly when this happens. They use a concept called dissipated energy—think of it as the total amount of "jiggling" or work the sand has to do to rearrange itself. The more you shake it, the more energy it absorbs, and the higher the water pressure gets.

The Old Problem: A Broken Ruler

For decades, engineers have used a famous formula (created by Berrill and Davis in 1985) to predict this water pressure. It's like a simple ruler that says: "The more energy you put in, the higher the pressure goes."

However, this old ruler has two major flaws:

It has no ceiling: The formula keeps going up forever. In the real world, water pressure can't exceed the weight of the soil above it (it can't go above 100%). The old formula sometimes predicts impossible numbers, like 150% pressure, which breaks the physics.
It has tangled knobs: The formula has two settings (parameters) that are supposed to be independent, but they are actually stuck together. If you change one, the other changes too. This makes it very hard to calibrate the model for different types of sand.

The New Solution: A Smart, Bounded Model

The authors of this paper (Gündoğan and Erken) decided to fix this ruler. They looked at data from real lab tests on two types of sand: Ottawa F65 (a standard, clean sand used in many studies) and Akpınar sand (a local sand from Istanbul).

Instead of just drawing a straight line on a graph, they realized the relationship between shaking energy and water pressure is more like a three-stage journey:

The Start: A slow, gentle rise.
The Middle: A rapid, steep climb.
The End: A smooth curve that flattens out exactly at the limit (100% pressure).

To capture this, they treated the math like a recipe. They broke the formula down into two parts:

The Growth Engine: A basic term that handles the rapid rise.
The Mediator: A "traffic controller" that slows the growth down as it approaches the limit, ensuring it never breaks the 100% ceiling.

The "Synthetic" Secret: A Mathematical Story

Here is the clever part of their method. They noticed that the shape of their new curve looked exactly like the solutions to a specific type of physics equation called an Ordinary Differential Equation (ODE).

Think of an ODE as a story about how a system changes over time.

The Left Side of the equation represents the sand's internal behavior (how it resists).
The Right Side represents the external shaking (the earthquake).

By rewriting their model as this equation, they could apply strict "rules of the road" (physical constraints):

Rule 1: The pressure must always go up (you can't un-shake the sand).
Rule 2: The pressure must stop rising exactly when it hits 100% (liquefaction).
Rule 3: At the moment it hits 100%, the speed of the rise must be zero (it levels off smoothly, not a jagged spike).

By forcing the math to obey these rules, they were able to eliminate the "tangled knobs." They reduced the whole complex model down to just one single number (called b) that they need to find. The rest of the math is fixed by the rules of physics.

What They Found

They tested this new "One-Number Model" against the lab data:

Accuracy: It matched the real lab results incredibly well (99% accuracy).
Safety: It never predicted impossible numbers (it stays between 0% and 100%).
Simplicity: It is much easier to use than the old model because you only need to tune one setting.

The "Aha!" Moment: Two Stages of Failure

The math revealed something interesting about how sand fails.

Early Stage: The sand's own internal structure (how tightly packed the grains are) controls the pressure rise. This is the "Homogeneous" part of the equation.
Late Stage: Once the sand gets very close to liquefying (around 75% pressure), it loses its "memory." The internal structure breaks down, and the pressure rise is driven almost entirely by the external shaking energy. This is the "Particular" part of the equation.

Summary

The authors took a messy, unbounded, and complicated formula for predicting sand liquefaction and rebuilt it. They used a mathematical "story" (an equation) to force the model to behave like real physics: rising smoothly and stopping exactly at the limit. The result is a simpler, safer, and more accurate tool for engineers to understand when sand will turn to liquid during an earthquake.

Technical Summary: A Physically Constrained Energy-Based Excess Pore Pressure Model for Saturated Sands

Problem Statement
Soil liquefaction remains a critical geotechnical hazard, traditionally assessed through stress-based, strain-based, pore-pressure-based, and energy-based methods. Among these, Energy-Based Methods (EBM) offer a distinct advantage by relating excess pore pressure generation to cumulative dissipated energy, thereby accommodating irregular seismic loading without the uncertainty of equivalent-cycle conversions. The most influential EBM formulation is the Berrill–Davis (1985) power law ( $r_u = \alpha W^\beta$ ). However, this study identifies two fundamental limitations in the Berrill–Davis model:

Unbounded Response: The power law is mathematically unbounded, allowing the excess pore pressure ratio ( $r_u$ ) to exceed unity, which is physically impossible for saturated sands where $r_u$ should asymptotically approach 1.0 at liquefaction.
Parameter Interdependence: The calibration parameters $\alpha$ and $\beta$ are often treated as independent, yet empirical studies show they are mutually correlated. This dependence complicates calibration, reduces physical interpretability, and hinders extrapolation.

Methodology
The authors propose a reformulated model for uniform, clean, saturated sands that addresses these limitations through a multi-step derivation process:

Data Analysis and Decomposition: Using cyclic triaxial data from Ottawa F65 sand (LEAP-2017) and Akpınar sand, the relationship between the normalized pore pressure ratio ( $r_u$ ) and normalized dissipated energy ( $W_N$ ) is analyzed on a log–log scale. The authors observe that the data does not follow a single straight line but rather consists of at least three linear segments.
Polynomial Approximation: The data is approximated by a sixth-order polynomial without a constant term. This function is decomposed into a primary growth term ( $W_N^b$ ) and a mediation function $g(W_N)$ .
Taylor Expansion and Simplification: The mediation function is further analyzed. A Taylor expansion of the exponential component reveals that a third-order polynomial adequately represents the data trends. This leads to a candidate model form: $r_u = W_N^b (a_0 + a_1 W_N + a_2 W_N^2 + a_3 W_N^3)$ .
Ordinary Differential Equation (ODE) Interpretation: Recognizing the mathematical structure, the model is reinterpreted as the solution to a non-homogeneous Cauchy–Euler ordinary differential equation. The left-hand side represents the system's internal response, while the right-hand side represents external loading conditions.
Physical Constraints: To eliminate unphysical fluctuations and parameter interdependence, the model is subjected to specific boundary conditions at the liquefaction point ( $W_N = 1, r_u = 1$ $W_{N} = 1, r_{u} = 1$ ):
- $r_u = 1$ at $W_N = 1$ .
- A horizontal tangent ( $r_u' = 0$ ) at $W_N = 1$ , ensuring the curve flattens as it reaches the limit.
- Monotonic increase ( $r_u' > 0$ ) throughout the loading process.
- A specific relationship between coefficients ( $a_2 = -3a_3$ ) derived from the vanishing of the non-homogeneous term at liquefaction.

These constraints analytically determine the coefficients of the mediation polynomial in terms of the single parameter $b$ , reducing the model to a form governed solely by $b$ .

Key Contributions

Physically Constrained Formulation: The proposed model is constructed to inherently satisfy $0 \le r_u \le 1$ , preventing the unphysical overprediction of pore pressure common in the Berrill–Davis model.
Single Calibration Parameter: By analytically linking the mediation coefficients to the growth exponent $b$ , the model eliminates the need for independent calibration of multiple parameters, resolving the interdependence issue found in previous formulations.
ODE Framework: The study provides a novel interpretation of the empirical pore pressure relationship as the solution to a non-homogeneous Cauchy–Euler ODE. This framework separates the response into a homogeneous part (governing internal material characteristics, dominant in early stages) and a particular part (governing external loading conditions).
Validation: The model is validated against high-quality cyclic triaxial data for Ottawa F65 and Akpınar sands across a range of relative densities (from loose to very dense), achieving $R^2 \ge 0.99$ .

Results

The simplified model (Eq. 20) successfully reproduces measured pore-pressure histories for both Ottawa F65 and Akpınar sands.
The parameter $b$ was found to vary with relative density (e.g., $b \approx 0.36$ for 51% relative density, $b \approx 0.43$ for 80–91% relative density).
The model captures the non-linear evolution of pore pressure, including the initial rapid rise and the subsequent plateauing as liquefaction is approached.
Comparison with the original Berrill–Davis model demonstrates that the proposed formulation offers a more consistent description of pore pressure evolution, particularly near the terminal condition of complete liquefaction, where the original model exhibits significant scatter and unbounded growth.
The application of the simplified model to the remaining LEAP-2017 dataset confirms its capability to simulate pore pressure generation across various relative densities and cyclic stress ratios, provided $b \ge 0.27$ to avoid local maxima.

Significance
The paper claims that this work offers a viable, compact, and physically consistent alternative to the classical Berrill–Davis power law for liquefaction studies involving uniform, clean sands. By bounding the response and removing parameter interdependence through an ODE-based derivation, the model simplifies calibration while maintaining high accuracy. The authors emphasize that while the formulation is empirical, its derivation from regression under rigorous physical constraints allows it to capture the underlying mechanism of liquefaction: an initial phase governed by internal granular characteristics (homogeneous solution) followed by a phase where the material loses its mechanical memory and requires substantial energy to reach full liquefaction (particular solution). The scope is explicitly restricted to uniform, clean, saturated sands.