Discovering Thermodynamically Admissible Dissipation Potentials via Grammar-Based Symbolic Regression

Grammar-Based Symbolic Regression Discovers Thermodynamically Consistent Material Laws

In the field of constitutive modeling—the science of defining how materials like rubber, metals, or biological tissues deform under stress—getting the math right is a matter of physical survival. If a model predicts a material will gain energy during a deformation cycle instead of dissipating it, the simulation is physically impossible.

Current data-driven approaches, such as black-box neural networks, often sacrifice this interpretability. Even when researchers use physics-encoded architectures, the resulting models remain opaque. They are difficult to extrapolate and cumbersome to integrate into standard finite element codes (numerical tools used to simulate how objects respond to physical forces). This paper proposes a way to bridge that gap. It uses symbolic regression (a method that searches for mathematical formulas) to "write" laws that are guaranteed to obey the laws of thermodynamics by construction.

The Problem

The status quo in material modeling is caught between two suboptimal poles. On one side are classical phenomenological models. These rely on expert intuition and ad hoc assumptions. They are highly interpretable and easy to deploy in simulation software. However, they often lack the expressivity needed to capture complex, nonlinear phenomena like the Payne effect (amplitude-dependent softening) in elastomers.

On the other side are modern machine learning approaches. Neural networks can approximate arbitrary constitutive maps with high precision. But they are essentially black boxes. Even "Physics-Informed Neural Networks" (PINNs) typically treat physical laws as soft penalty terms in a loss function. This means they offer no formal guarantee of thermodynamic admissibility (ensuring the model obeys energy conservation laws) outside the training range. For an engineer deploying a model in a safety-critical application, a model that might "invent" energy is a liability.

How It Works

The authors tackle this by shifting the search from weight optimization to searching the space of mathematical expressions. The core innovation is a composition-extended convexity-preserving grammar, denoted as $G_{cvx}^{comp}$.

The methodology follows these logical steps:

1. Constraint Encoding: Instead of penalizing violations, the authors encode the requirements of the Clausius–Duhem inequality directly into the grammar's production rules. Specifically, they enforce convexity and non-negativity. This ensures any generated expression is thermodynamically admissible by design.
2. Hierarchical Construction: The grammar builds candidate potentials using a structured hierarchy. It starts with "inner primitives" (basic building blocks like power laws). It then allows them to be composed with "outer functions" (convex, non-decreasing functions). As shown in, this creates an expression tree.

Every leaf and internal node respects the necessary mathematical properties. 3. Evolutionary Search: The system uses genetic programming to evolve a population of these expressions. They use term-level crossover and mutation. These operations swap or modify entire blocks of the equation. This ensures the "offspring" stay within the valid grammar. 4. Hybrid Optimization: For every candidate expression, the authors perform a two-mode local search. If the expression is smooth, they use L-BFGS-B (a quasi-Newton optimization algorithm). If it contains non-differentiable elements like the Macaulay bracket (a function used to model yield thresholds), they switch to the gradient-free Nelder–Mead simplex. 5. Trajectory Rollout: Fitness is not measured by simple curve fitting. Instead, the system performs a "forward rollout." It integrates the resulting evolution law through time. This checks if the predicted stress trajectories match the experimental data .

Numbers

The authors validate the framework against three synthetic benchmark families and real-world experimental data. The results suggest the symbolic approach is efficient at recovering ground-truth physics.

In the Newtonian viscosity benchmark (E1), the authors report a near-exact recovery rate of 97.8% under moderate noise. For the more complex power-law family (E2), they achieve a 93.3% exact recovery rate. Most impressively, for the Bingham viscoplastic model (E3), the framework achieved a 100% exact recovery rate [Table 9]. This means the system perfectly identified the yield threshold and the flow law.

In the experimental validation using a synthetic elastomer, the GP-discovered potential outperformed the industry-standard baseline. When compared against a calibrated linear Zener model (a standard linear viscoelastic model), the symbolic model reduced the multi-history stress error (nRMSE) from 0.1103 to 0.0822. This represents a ~25% improvement in predictive accuracy for the stress signal. Furthermore, as seen in, the discovered potential successfully captured the amplitude-dependent softening of the material's dynamic moduli. The linear Zener baseline is structurally incapable of representing this nonlinear behavior.

What's Missing

While the results are compelling, there are clear boundaries to the current implementation.

• Dimensionality Limits: The framework is currently restricted to a single internal variable ($N=1$). Real-world complex materials involve multiple coupled dissipative mechanisms. Extending this would turn identification into a much harder simultaneous state estimation problem.
• Tensor Generalization: The current derivation focuses on one-dimensional or simplified scalar settings. Moving to a full three-dimensional continuum requires generalizing the grammar. It must operate on invariants (scalar values derived from tensors) of the stress and strain tensors.
• Assumption of Known Elasticity: The methodology assumes the elastic free energy $\psi$ is known a priori. In many industrial scenarios, you cannot characterize the equilibrium state separately. You need to discover both the elastic and dissipative laws simultaneously.

Should You Prototype This

Yes, but with caveats. If you are working in a domain where "black-box" models hit a wall due to physical inconsistency, this is a high-value target. Because the output is a closed-form symbolic expression, you can deploy it directly. You can drop the discovered law into an ABAQUS or LS-DYNA subroutine without a custom neural network inference engine.

However, do not expect a "plug-and-play" experience for 3D tensors or multi-mechanism materials yet. Start by prototyping this on 1D or 2D simplified loading cases. Test if your material's nonlinearity is well-captured by the provided primitive library. If your goal is to replace a tuned phenomenological model with something more accurate and mathematically rigorous, this framework provides a viable path.

Figures from the paper