Engineering-Oriented Symbolic Regression: LLMs as Physics Agents for Discovery of Simulation-Ready Constitutive Laws

Abstract

The discovery of constitutive laws for complex materials has historically faced a dichotomy between high-fidelity data-driven approaches, which demand prohibitive full-field experimental data, and traditional engineering fitting, which often yields numerically unstable models outside calibration regimes. In this work, we propose an Engineering-Oriented Symbolic Regression (EO-SR) framework that bridges this gap by leveraging Large Language Models (LLMs) as "Physics-Informed Agents." Unlike unconstrained symbolic regression, our framework utilizes an LLM Agent to zero-shot synthesize executable physical constraints -- specifically thermodynamic consistency and frame indifference -- transforming the search process from mathematical curve-fitting into a physics-governed discovery engine. We validate this approach on the hyperelastic modeling of rubber-like materials using standard Treloar datasets. The framework autonomously identifies a novel hybrid constitutive law that combines a Mooney-Rivlin linear base with a rational locking term. This discovered model not only achieves high predictive accuracy across multi-axial deformation modes (including zero-shot prediction of pure shear) but also guarantees unconditional convexity. Finite element validation demonstrates that while industry-standard models (e.g., Ogden N=3) fail due to numerical singularities under severe transverse compression, the EO-SR-discovered model maintains robust convergence. This study establishes a generalized, low-barrier pathway for discovering simulation-ready constitutive closures that satisfy both data accuracy and rigorous physical laws.

Figures (7)

• Figure 1: Overview of the Engineering-Oriented Symbolic Regression (EO-SR) framework. The workflow is designed to bridge physical rigor with engineering practicality through three integrated stages. Stage 1 (Skill Injection): A Foundation LLM acts as a domain agent, translating natural language problem descriptions into a structured constraint set $\mathcal{S}=\{\mathcal{T}, \mathcal{O}, \mathcal{C}\}$, defining coordinate transformations, operator whitelists, and physical inequalities (e.g., Drucker stability). Stage 2 (Constrained Discovery): The symbolic regression engine performs a genetic search governed by a physics-informed fitness function, balancing data fidelity with the penalties derived from $\mathcal{S}$. Stage 3 (FEM Verification): Discovered models undergo rigorous finite element validation. Unlike classical empirical models which risk numerical divergence (Path A) or idealized models that lack flexibility (Path B), the EO-SR approach ensures the identification of constitutive laws that are both empirically accurate and numerically robust (Optimal Path).
• Figure 2: Pareto frontier analysis of discovered constitutive models. The plots illustrate the trade-off between mathematical complexity (number of nodes) and accuracy. (Left) Fitting MSE on training data (Uniaxial + Equibiaxial). (Right) Generalization MSE on unseen Pure Shear data. The proposed model (Eq16, blue triangle) identifies an optimal "elbow" point, achieving a balance between sparsity and precision. Note that the unconstrained search found a lower-error solution (Sqrt-Eq, red cross), but it was rejected due to violation of the convexity (Drucker stability) constraint.
• Figure 3: Performance comparison of the discovered model against industry benchmarks. The stress-stretch responses are evaluated under: (Left) Uniaxial Tension (Training), (Center) Equibiaxial Tension (Training), and (Right) Pure Shear (Zero-shot Prediction). The proposed SR model (Eq16, solid red line) demonstrates superior multi-axial capability compared to the Yeoh model (N=3, dashed blue), which fails to capture the equibiaxial response. Crucially, in the pure shear regime (not used for calibration), the SR model accurately predicts the experimental data ($MSE \approx 0.0048$), verifying its physical generalization capability compared to polynomial overfitting.
• Figure 4: Visualization of thermodynamic consistency via energy landscape curvature.(Left) The proposed convex model (Eq16) exhibits a positive-definite Hessian (smooth bowl shape) across the entire invariant space, guaranteeing unconditional stability. (Right) The unconstrained SR solution (Sqrt-Eq) reveals regions of negative curvature (concavity), indicating physical instability despite high fitting accuracy. This comparison highlights the critical role of physics-informed constraints in filtering out mathematical artifacts.
• Figure 5: Analysis of Drucker stability and finite extensibility limit. The plot shows the evolution of tangent stiffness ($d^2W/d\lambda^2$) under uniaxial tension. Inset: The proposed model maintains strict positivity (min stiffness $\approx 0.30$ MPa), satisfying the stability prerequisite. Main Plot: At large deformations, the Ogden model (dashed green) predicts unbounded polynomial growth. In contrast, the proposed SR model (solid blue) successfully identifies the molecular chain rupture limit, exhibiting an asymptotic stiffness surge at $\lambda \approx 8.77$. This physical "locking" behavior prevents non-physical infinite stretching in simulation.