Unsupervised discovery of interpretable hyperelastic constitutive laws

Moritz Flaschel; Siddhant Kumar; Laura De Lorenzis

Source-linked AI summary

Unsupervised discovery of interpretable hyperelastic constitutive laws

Moritz Flaschel, Siddhant Kumar, Laura De Lorenzis

arXiv:2010.13496v3Oct 26, 2020cs.CE

TL;DR

The paper addresses automated constitutive discovery when stress data are unavailable and interpretable models are desired. It proposes one-shot unsupervised sparse regression constrained by bulk and boundary equilibrium, and numerically demonstrates accurate discovery across five hyperelastic models plus accurate response prediction when a feature is missing.

Problem
Existing constitutive-discovery approaches commonly use stress data or produce uninterpretable model-free or black-box representations.
Method
The method fits a sparse invariant-based hyperelastic model using displacement and global force data while enforcing bulk and boundary momentum balance and physics-based constraints.
Results
The approach accurately discovers five hyperelastic models of different complexity and can surrogate a missing library feature while accurately reproducing strain-energy and stress responses.
Takeaways & Limitations
Interpretable sparse constitutive laws can be discovered from experimentally measurable displacement and force data without stress-field measurements.
Takeaways & Limitations
The initial study assumes homogeneous, isotropic, hyperelastic materials and uses empirical physical-admissibility checks rather than pursuing properties such as polyconvexity or coercivity.

Abstract

from arXiv · show

We propose a new approach for data-driven automated discovery of isotropic hyperelastic constitutive laws. The approach is unsupervised, i.e., it requires no stress data but only displacement and global force data, which are realistically available through mechanical testing and digital image correlation techniques; it delivers interpretable models, i.e., models that are embodied by parsimonious mathematical expressions discovered through sparse regression of a large catalogue of candidate functions; it is one-shot, i.e., discovery only needs one experiment - but can use more if available. The problem of unsupervised discovery is solved by enforcing equilibrium constraints in the bulk and at the loaded boundary of the domain. Sparsity of the solution is achieved by l_p regularization combined with thresholding, which calls for a non-linear optimization scheme. The ensuing fully automated algorithm leverages physics-based constraints for the automatic determination of the penalty parameter in the regularization term. Using numerically generated data including artificial noise, we demonstrate the ability of the approach to accurately discover five hyperelastic models of different complexity. We also show that, if a "true" feature is missing in the function library, the proposed approach is able to surrogate it in such a way that the actual response is still accurately predicted.

1. Introduction

Existing data-driven approaches either bypass constitutive modeling or encode it with black-box models, often requiring stress data. This paper instead proposes one-shot, unsupervised discovery of interpretable isotropic hyperelastic laws using displacement and global force measurements with physics-based constraints.

Model-free and neural-network approaches replace or encode constitutive relations, limiting material-model interpretability.
Most existing approaches use supervised input-output data and therefore rely on stress measurements, which are difficult to obtain in general mechanical tests.
The proposed method automatically discovers isotropic hyperelastic constitutive models rather than bypassing or surrogating them.
Discovery uses full-field displacements and global forces, enforces physics constraints, and produces parsimonious mathematical expressions through sparse regression.
The method is one-shot: it needs one experiment for discovery, although additional experiments can be used.

2.1. Problem setting

The problem setting uses a quasi-static displacement-controlled test with full-field displacement measurements and global reaction forces. The initial study assumes homogeneous, isotropic, hyperelastic materials and seeks discovery from one sufficiently complex experiment.

The test is quasi-static, with Dirichlet and Neumann boundary conditions applied on complementary boundary portions.
Under displacement control, Neumann conditions are homogeneous and machine-measured forces correspond to reaction forces on Dirichlet boundaries.
The available data are full-field nodal displacements from DIC and net reaction forces measured on loaded boundary portions.
One-shot discovery requires an experiment complex enough to activate diverse deformation modes; the study uses biaxial tension of a plate with a hole.
The initial scope assumes homogeneous, isotropic, hyperelastic material behavior in a two-dimensional setting.

2.2. Displacement field approximation via finite element mesh

Measured nodal displacements are associated with a finite element mesh to construct a continuous displacement field. This field is differentiated with respect to reference coordinates to obtain deformation quantities needed by the discovery procedure.

The measured displacement grid is represented by a finite element mesh sharing the known nodal set.
The displacement field is approximated using nodal displacements and associated shape functions.
The deformation gradient is obtained by differentiating the finite element displacement approximation.
The gradient operator is taken with respect to the reference coordinates.

2.3. Material model library

The constitutive law is represented through an isotropic strain-energy density expressed as a linear combination of invariant-based candidate features. The library combines generalized Mooney–Rivlin, volumetric, and logarithmic features while preserving physical structure and interpretability.

For hyperelasticity, the target is the strain energy density W, from which the first Piola–Kirchhoff stress is derived.
Objectivity and isotropy reduce W to a function of the right Cauchy–Green tensor and its invariants.
The library represents W as a linear combination of nonlinear functions of the invariants I1, I2, and I3 with unknown parameters θ.
With N = M = 7, the catalogue contains 43 features spanning generalized Mooney–Rivlin, volumetric, and logarithmic terms.
The logarithmic feature can be deliberately excluded to test whether the method can surrogate missing constitutive features.
The feature-library design supports physical interpretability and can encode objectivity, isotropy, and a stress-free reference configuration automatically.

2.4. Equilibrium constraints for unsupervised discovery

The method identifies constitutive parameters by enforcing linear momentum balance in the material bulk and reaction-force balance on loaded boundaries, using displacement and global force measurements without stress-field data.

Bulk and boundary equilibrium: Equilibrium constraints substitute for stress data, enabling unsupervised discovery from interpolated displacement measurements and a material-model library.The constraints enforce linear momentum balance in the bulk and at the loaded boundary.
Bulk and boundary equilibrium: The weak formulation is preferred because it avoids the double spatial derivatives of the strong formulation and is less sensitive to displacement noise.The formulation uses admissible test functions that vanish on the Dirichlet boundary.
Bulk and boundary equilibrium: Free-degree-of-freedom equations are assembled into a linear system A_freeθ = b_free under homogeneous body-force and material-homogeneity assumptions.The material parameters are spatially constant and angular momentum balance is automatically fulfilled in the described setting.
Bulk and boundary equilibrium: Measured reaction-force sums over boundary degree-of-freedom subsets provide additional linear equations A_fixθ = b_fix for the loaded Dirichlet boundary.Each subset corresponds to one loaded side and one displacement direction, matching experimentally available global force measurements.
Bulk and boundary equilibrium: The overdetermined bulk-and-boundary system is solved collectively in least squares, with λ_r weighting reaction-force balance against interior balance.For the chosen discretization, λ_r = 100 compensates for the smaller number of boundary equations.
Multiple load steps: Multiple load steps can be incorporated by additively combining their least-squares losses, although the derivation is shown for one load step.The superscript l denotes evaluation at load step l = 1, . . . , L.

2.5. Parsimony and sparsity promotion

The method promotes parsimonious constitutive laws by regularizing feature coefficients, balancing interpretability against computational tractability and sensitivity to noise.

Motivation for sparsity: The unregularized inverse problem is highly ill-posed because displacement-dependent feature matrices are nonlinear, poorly conditioned, and sensitive to measurement noise.The lack of stress-field measurements contributes to the unsupervised problem's ill-posedness.
Motivation for sparsity: Dense coefficient vectors reduce interpretability, extrapolation quality, and physical admissibility, motivating models with as few active features as possible.This motivation is based on prior supervised studies of material models with many parameters.
Regularization choice: ℓ_p regularization promotes sparsity by controlling which feature-library terms remain active through the penalty parameters λ_p > 0 and p ∈ [0, 1].The ℓ_0 penalty counts active features, while ℓ_1 is a computationally tractable convex relaxation.
Regularization choice: Smaller p more aggressively promotes sparsity and interpretability but increases computational complexity through greater non-convexity; p = 1/4 was used in the numerical experiments.The choice p = 1/4 was reported as reasonable for the scope of this work.

2.6. Numerical optimization strategy

The optimization uses fixed-point iterations for non-convex ℓ_p regularization, thresholding to remove inactive features, and physics-based admissibility checks to select the penalty parameter.

Fixed-point optimization: For 0 < p < 1, ℓ_p regularization makes the objective non-convex and non-smooth, so conventional gradient-based and Newton-like methods perform poorly.The paper generalizes a weighted ℓ_2 fixed-point strategy to the ℓ_p-regularized problem.
Fixed-point optimization: The nonlinear optimality system is solved iteratively by replacing nonlinear coefficient terms with weights computed from the previous iterate, producing a linear system at each step.The resulting sequence of coefficient estimates is advanced through a fixed-point scheme.
Fixed-point optimization: The algorithm runs parallel fixed-point schemes from different random initial guesses because non-convexity can lead to different local minima.Iterations continue until convergence or a prescribed maximum is reached.
Penalty selection: The penalty parameter is increased multiplicatively from a small initial value until the resulting model passes physical admissibility checks.The strategy avoids models that are either too rich or too sparse and simple.
Physical admissibility: Admissibility checks require nonnegative strain energy on observed deformation gradients and monotonically increasing energy along selected tension, compression, and shear paths.The path condition is tested on finitely many samples up to γ_max, and violations reject the model.
Physical admissibility: The checks are empirical rather than rigorous, while the discovered expressions permit later semi-analytical verification of properties such as polyconvexity or coercivity.Those more rigorous checks were not pursued in this work.
Feature elimination: Thresholding permanently removes coefficients below a prescribed magnitude, after which regularization is deactivated and the reduced problem is solved unregularized.This eliminates inactive features from the initial library.

3. Numerical benchmarks

Numerical benchmarks test the unsupervised discovery algorithm on five hyperelastic solids using simulated plate-with-hole data, including representative displacement noise. Exact or accurate recovery depends on noise level and deformation coverage.

3.1. Data generation: A quadrant of a square plate with a hole under displacement-controlled asymmetric biaxial tension generates the benchmark data.The simulations use plane-strain conditions, symmetry boundaries, nodal displacements, and global reaction forces.
3.1. Data generation: The benchmark set comprises NH2, NH4, IH, HW, and GT hyperelastic models spanning quadratic, biquadratic, cubic, and logarithmic energy terms.The exact strain-energy densities are used as hidden models to be recovered from displacement and reaction-force data.
3.1. Data generation: Displacement noise is modeled as zero-mean Gaussian noise, with σ = 10^-4 and σ = 10^-3 representing low and high normalized noise levels.Noisy displacements are denoised by kernel ridge regression before entering the discovery algorithm.
3.2. Results: At σ = 0, all five models are exactly discovered, while at σ = 10^-4 all correct features are identified with only small parameter differences.The negligible differences in strain-energy and selected stress-path plots corroborate the accuracy of these discoveries.
3.2. Results: At σ = 10^-3, NH2 and NH4 agree well with the hidden models, whereas IH and HW deviate especially under large pure shear; GT still agrees well despite a different functional form.The high-noise IH and HW models therefore approximate the true response mainly in the small-deformation regime.
3.2. Results: When the logarithmic feature is excluded, GT⋆ replaces it with the linear feature (Ī2 − 3) while still reproducing strain energy and stress responses accurately.The same surrogate behavior occurs at high noise even when the logarithmic feature is available.

4. Conclusions and outlook

The paper presents an unsupervised sparse-regression method for discovering interpretable hyperelastic constitutive laws from measurable displacement and global-force data. Numerical experiments recover five models from one loading experiment and show accurate response prediction when a library feature is missing.

4. Conclusions and outlook: The method enforces interior and boundary momentum balance while using ℓp regularization, thresholding, and physics-based automatic penalty selection.It discovers parsimonious constitutive expressions from a large candidate-function library.
4. Conclusions and outlook: Numerical finite-element experiments show that one loading experiment can correctly identify five chosen hyperelastic models from displacement and global-force data.The data include cases with and without artificial noise.
4. Conclusions and outlook: When a true feature is absent from the candidate library, the approach can surrogate it while accurately predicting the actual material response.This extends the demonstrated capability beyond direct feature identification.
4. Conclusions and outlook: Future developments include experimental two- and three-dimensional applications and extensions to non-homogeneous, anisotropic, and irreversible materials.These directions mark the scope beyond the present work.

Appendix A. Denoising via kernel ridge regression (KRR)

Appendix A describes kernel ridge regression for denoising noisy displacement components using a high-dimensional feature representation, while avoiding direct inversion in feature space through the kernel trick.

KRR formulation: KRR represents the denoised scalar field as a linear combination in a potentially infinite-dimensional feature space.The feature dimension d is much larger than the number of measurements nn, and may be infinite.
KRR formulation: The coefficients are fitted by minimizing measurement discrepancies with ℓ2 regularization controlled by ξ > 0.The regularization balances fidelity to noisy measurements against the coefficient norm.
Kernel reformulation: Directly solving the feature-space system requires inverting a d × d matrix, which is intractable when d is infinite.This motivates the equivalent kernel formulation.
Kernel reformulation: The kernel reformulation instead requires inversion of an nn × nn matrix, making computation tractable when d ≫ nn.The resulting denoised field is obtained by substituting the coefficient solution into the field ansatz.
Kernel choice: The kernel trick eliminates the need to define the feature map explicitly by using a valid positive-definite kernel.The chosen RBF kernel has lengthscale χ > 0, and the resulting regression is non-parametric.
Implementation: Each displacement component is denoised separately, with ξ and χ optimized as the method’s hyperparameters.The appendix also notes that KRR is deterministic and was observed to be faster than Gaussian-process regression for moderately sized problems.

Appendix B. Numerical protocol and computational costs

The numerical protocol varies the number of load steps with model complexity and reports approximately ten minutes of computation for one fixed-point initialization.

Numerical protocol: NH2 and NH4 use L = 4 loadsteps, whereas IH, HW, and GT use L = 8.The higher-order models use more load steps to reach strain levels comparable to those of the simpler models.
Computational costs: One execution with a single initial guess requires approximately ten minutes on an average modern processor.Runs using different initial guesses can be distributed in an embarrassingly parallel manner across a cluster.