# Learning a hyperelastic constitutive model from 3D experimental data

Mathys Bourdyot, Marius Compans, Raphaël Langlois, Benjamin Smaniotto,
Emmanuel Baranger, Clément Jailin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mathys Bourdyot, Marius Compans, Raphaël Langlois, Benjamin Smaniotto, Emmanuel Baranger, et al.. Learning a hyperelastic constitutive model from 3D experimental data. Computer Methods in Applied Mechanics and Engineering, 2025, 450, pp.118592. 10.1016/j.cma.2025.118592. hal-05390159v2

HAL Id: hal-05390159

https://hal.science/hal-05390159v2

Submitted on 15 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Learning a hyperelastic constitutive model from 3D experimental data

M. Bourdyot<sup>a</sup>, M. Compans<sup>ID a</sup>, R. Langlois<sup>ID b</sup>, B. Smaniotto<sup>ID a,b</sup>, E. Baranger<sup>ID b</sup>, C. Jailin<sup>ID b,*</sup>

## ARTICLE INFO

**Keywords:**

Machine learning

Physics augmented neural network

Digital volume correlation

Constitutive behavior

Tomography

## ABSTRACT

A hyperelastic constitutive law is experimentally learned and validated from three-dimensional full-field kinematic data. 3D-printed thermoplastic polyurethane specimens were subjected to monotonic uniaxial tension within a laboratory micro-computed tomography system. Global Digital Volume Correlation yields volumetric displacement fields, reaching up to 35 % axial strain. A Physics-Augmented Neural Network (PANN) architecture, embedding polyconvexity constraints through an input-convex neural network, is trained in an unsupervised manner by enforcing mechanical equilibrium via the EUCLID loss, using only measured displacements and global reaction forces. This three-dimensional formulation eliminates the need for depth-related assumptions inherent in 2D measurements and provides improved representation of boundary conditions, which is critical for equilibrium-based training. For comparison, Neo-Hookean and Saint-Venant-Kirchhoff models are identified in parallel. The PANN achieves superior equilibrium satisfaction, outperforming both linear and Neo-Hookean benchmarks. In addition to the original specimen, a second experiment on a different sample was performed under comparable loading conditions. This configuration provides an independent validation dataset acquired under slightly varying experimental conditions, thereby testing the robustness and transferability of the learned constitutive model. This work presents the first experimental demonstration of unsupervised PANN model training and establishes a practical protocol for data-based hyperelastic characterization from volumetric measurements.

## 1. Introduction

The reliable prediction of material behavior remains a cornerstone of solid and structural mechanics. Classical constitutive laws, obtained within a thermodynamic framework and expressed through fixed analytical potentials, can be calibrated with a handful of physical parameters but struggle to reproduce complex behaviors such as strong non-linearities, micro-structural heterogeneity, and anisotropy. Over the last decade, a rapid dissemination of machine-learning (ML) techniques within mechanics has opened new possibilities for flexible material behavior descriptions [1]. Different methods have emerged with the development of machine-

<footnote><sup>a</sup> ENS Paris-Saclay, DGM - Department of Mechanical Engineering, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS - Laboratoire de Mécanique Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: mathys.bourdyot@ens-paris-saclay.fr (M. Bourdyot), marius.compans@ens-paris-saclay.fr (M. Compans), raphael.langlois@ens-paris-saclay.fr (R. Langlois), benjamin.smaniotto@ens-paris-saclay.fr (B. Smaniotto), emmanuel.baranger@ens-paris-saclay.fr (E. Baranger), clement.jailin@ens-paris-saclay.fr (C. Jailin).</footnote>
learning-based identification strategies [2]. These approaches aim at relaxing potential shape, not only adjusting the parameters but the shape of the function itself.

Recent progress in artificial intelligence shows that machine learning techniques, and neural networks in particular, can infer constitutive laws directly from experimental observations [2,3]. Their high-dimensional approximation capability makes it possible to capture complex non-linear material responses. Difficulties linked to physical inconsistencies and limited access to stress data can be resolved by embedding physics within the learning process [4]. One strategy consists of incorporating physical terms into the loss function [5,6]; the Thermodynamics-based Artificial Neural Network (TANN) framework is a representative example [7,8]. A different route, shared by several related frameworks such as Constitutive Artificial Neural Networks (CANNs), Mechanics-Informed, and Physics-Augmented Neural Networks (PANNs), integrates the mechanical constraints directly into the neural architecture itself [9–13]. Those architectures often rely on Input-Convex Neural Networks (ICNN) to enforce the polyconvexity of the free-energy potential [14]. It should also be acknowledged that these advances resonate with a wider research effort in computational mechanics, where multiple teams are actively exploring data-based and physics-augmented strategies for constitutive modeling.

Experimentally, stress cannot be directly measured and must be inferred indirectly from global reaction forces measured by load cells. Consequently, the so-called supervised learning approaches, which require stress labels at every point, are difficult to apply in practice. To overcome these obstacles, so-called “unsupervised” learning techniques have been developed. The frameworks EUCLID, Efficient Unsupervised Constitutive Law Identification and Discovery [15,16], a generalized implementation of the Equilibrium Gap method [17], are “unsupervised” trainings using only full-field displacement and a small set of global reaction forces (and not stresses). By enforcing equilibrium in its weak form, minimizing the residual between internal forces (predicted by the network) and external loads on a finite-element mesh, the network learns to map measured deformations to balanced stress fields without requiring direct stress measurements. This EUCLID loss has been recently extended with the virtual field method [18]. Several related methods employ a similar equilibrium-gap principle [19] or Equilibrium-Based Convolutional Neural Networks (ECNN) [20], differ in their network architectures and training algorithms. Recently, a modified constitutive relation error (mCRE) formulation has been proposed to improve robustness to measurement noise when training Physics-Augmented Neural Networks [21,22].

Experimental mechanics has evolved with more and more high-resolution imaging and quantitative image processing techniques [23]. In particular, Digital Image Correlation (DIC) on specimen surfaces [24,25] measures full-field displacement fields with high accuracy by registering a reference image with sequences of images captured under load in insitu experiments. Its extension in 3D, Digital Volume Correlation (DVC) [26,27] allows for the measurement of displacement fields from 3D volumetric images acquired, for example, thanks to laboratory X-ray tomography. The use of DVC enables the collection of a rich amount of deformation fields on a single test. These large datasets are an asset to train the PANN and data-based methods.

Despite the improvements of (1) ML model and (2) full-field measurements, the application of PANN to actual experimental remains limited both in quantity and in experimental configuration used: i.e., to one-dimensional learning [28–31] or two-dimensional configurations with unsupervised training [18,20,32,33]. Surface-based DIC studies require having strong depth assumptions, such as plane-stress, plane-strain, and incompressibility, that contradict the objective of reducing modelling hypotheses. Using 3D data fields eliminates the need for additional depth hypotheses and therefore constitutes an ideal foundation for unsupervised PANN training. In addition, it also provides richer information and more measured free surfaces that are precious equilibrium information for unsupervised training.

The aim of the present study is to learn an isotropic hyperelastic constitutive law directly from three-dimensional experimental data, without any stress labels. Polyurethane specimens produced by additive manufacturing were loaded monotonically in tension while being imaged in a lab computed-tomograph. This work proposes an experimental protocol that combines conventional micro-computed tomography with mechanically regularised DVC. The PANN was trained using only the measured displacement fields and the global reaction force. The predictive capability of the identified model is assessed against Neo-Hookean and linear-elastic references and is further evaluated on an unseen loading path. From the best of the authors' knowledge, this study is the first experimental demonstration of constitutive law learning in three dimensions.

## 2. Methods

Vectors and second-order tensors are denoted in bold upright font (e.g., **u**, **F**, **P**, **I**), while scalars remain italic. The reference gradient is **∇**. We use tr(·), cof(·), and det(·) for the trace, cofactor, and determinant operators, respectively. All quantities are expressed in the reference configuration.

### 2.1. PANN Architecture model

The neural network model used in this application is a PANN designed for isotropic hyperelastic applications. Its architecture is based on a neural network composed of dense layers that take as inputs the deformation tensor **F** and outputs the first Piola-Kirchhoff stresses **P**.

The strain-energy density is expressed as a function of the right Cauchy-Green deformation tensor $C = F^T F$, which ensures material frame indifference (objectivity) of the formulation. Following the study of [12], which has become a standard formulation in the PANN literature, the strain energy density is written as

$$W = W_{ICNN} + W_{\text{energy}} + W_{\text{stress}} + W_{\text{growth}} \quad (1)$$
The first term corresponds to the energy output of an input convex neural network $W_{ICNN}(I)$ [14], with I the input set of deformation invariants. The following two terms correspond to energy and stress normalization, and the final term to the energy growth. Various physical constraints are implemented as presented in pioneer papers in the literature [10,12]. The network input is the deformation gradient F, from which a dedicated hard-coded layer computes the invariants

$$I_1 = \text{tr}(C), \quad I_2 = \text{tr}(\text{cot}(C)), \quad I_3 = \text{det}(C) = J^2. \qquad (2)$$

As C is convex in F, $I_1$ is convex in F, $I_2$ is convex in $\text{cot}(F)$, and J is convex in $\text{det}(F)$. The extra invariant $I_1^* = -2J$ is included to allow the network to represent compressive (negative) stresses while preserving polyconvexity [10,12]. The scalar output of the network is the strain-energy density W (F). Thermodynamic consistency is enforced by obtaining the first Piola-Kirchhoff stress through

$$\mathbf{P}(\mathbf{F}) = \frac{\partial W(\mathbf{F})}{\partial \mathbf{F}}. \qquad (3)$$

To guarantee stability and positivity of the energy potential, the hidden layers are constructed as an Input-Convex Neural Network (ICNN) [14]: each layer uses the softplus activation $\Lambda(y) = \ln(1 + e^y)$ and enforces non-negative weights, $\omega^{(i)} \ge 0$. Writing $\mathbf{z}^{(i)}$ for the vector of pre-activation outputs at layer i, the layer-to-layer mapping is

$$\mathbf{z}^{(i)} = \Lambda(\omega^{(i)} \mathbf{z}^{(i-1)} + \mathbf{b}^{(i)}), \qquad (4)$$

with $\mathbf{b}^{(i)}$ the bias vector. By enforcing convexity in each of the invariants $I = (I_1, I_2, I_3, I_1^*)$, the network preserves polyconvexity of W with respect to F, $\text{cot}(F)$ and $\text{det}(F)$. Normalization conditions

$$W_{\text{energy}} = -W_{\text{ICNN}}|_{F=I}, \quad W_{\text{stress}}(J) = -\eta(J-1), \qquad (5)$$

are imposed by subtracting the energy and stress evaluated at the reference configuration, with

$$\eta = 2 \left( \frac{\partial W_{\text{ICNN}}}{\partial I_1} + 2 \frac{\partial W_{\text{ICNN}}}{\partial I_2} + \frac{\partial W_{\text{ICNN}}}{\partial I_3} + \frac{\partial W_{\text{ICNN}}}{\partial I_1^*} \frac{\partial I_1^*}{\partial I_3} \right) \Big|_{F=I} \qquad (6)$$

Finally, a growth term is added to follow state-of-the-art literature on PANN architecture:

$$W_{\text{growth}}(J) = \left(J + \frac{1}{J} - 2\right)^2. \qquad (7)$$

Unlike the work of Linka et al. [13] introducing an incompressibility constraint, here, compressible behavior is studied.

## 2.2. Traditional constitutive models

In order to have a comparison with the PANN model, two additional traditional models are used [34,35]. Because two models are used, it is clear that two sets of elastic parameters are to be defined, each corresponding to the optimal (in the EUCLID loss sense).

**Saint-Venant-Kirchhoff (SVK).** As a first benchmark, we employ the Saint-Venant-Kirchhoff model, which is a simple linear elastic model expressed using the Green-Lagrange strain

$$\mathbf{E} = \frac{1}{2}(\mathbf{F}^{\top}\mathbf{F} - \mathbf{I}). \qquad (8)$$

The quadratic energy density is then

$$W_{\text{SVK}} = \frac{\lambda_{\text{SVK}}}{2} [\text{tr}(\mathbf{E})]^2 + \mu_{\text{SVK}} \text{tr}(\mathbf{E}^2), \qquad (9)$$

where $\mu_{\text{SVK}}$ is the shear modulus and $\lambda_{\text{SVK}}$ the first Lamé parameter for this model. The SVK model provides a simple baseline against which to assess more advanced, nonlinear constitutive descriptions.

**Neo-Hookean law (NH).** To evaluate the model with a more complex non-linear law, the Neo-Hookean hyperelastic model is used. The strain energy density function W for the NH material is expressed in terms of the Lamé parameters and deformation invariants as

$$W_{\text{NH}} = \frac{\mu_{\text{NH}}}{2}(I_1 - 3) - \mu_{\text{NH}} \log(J) + \frac{\lambda_{\text{NH}}}{2} \log^2(J), \qquad (10)$$

where $\mu_{\text{NH}}$ and $\lambda_{\text{NH}}$ are the Lamé coefficients of this NH law. The NH model is widely used due to its relative simplicity and ability to capture essential non-linear elastic responses under moderate deformations. In this work, such a model is appropriate since the test material exhibits a richer non-linear elastic behavior.
## 2.3. Unsupervised constitutive training

In this framework, the training loss function quantifies the mechanical equilibrium error at the nodal scale. A loss value close to zero means that the predicted stress field satisfies the static equilibrium conditions, both inside the sample (in internal nodes) and at its boundaries. Physically, this ensures that the internal forces induced by the material response exactly compensate the applied external loads, without requiring direct stress measurements.

Consider a displacement field $\mathbf{u}(\mathbf{x})$ defined on the reference positions $\mathbf{x}$. This displacement field may come from DIC/DVC measurement procedures and is supposed to be known [27]. The displacement field is expressed using a finite element mesh, composed of $N_n$ nodes, with finite element shape functions $\psi(\mathbf{x})$ such that:

$$ \mathbf{u}(\mathbf{x}) = \sum_{i=1}^{N_n} \mathbf{u}_i \psi_i(\mathbf{x}), \qquad (11) $$

with $\mathbf{u}_i$ being the nodal displacement amplitudes. The deformation gradient field is then approximated as follows:

$$ \mathbf{F}(\mathbf{x}) = \mathbf{I} + \sum_{i=1}^{N_n} \mathbf{u}_i \otimes \nabla \psi_i(\mathbf{x}), \qquad (12) $$

where $\mathbf{I}$ is the identity matrix and $(\cdot)$ the gradient operator.

Considering a constitutive model $\mathcal{M}$ mapping deformation $\mathbf{F}$ onto first Piola-Kirchhoff stress $\mathbf{P}(\mathbf{F}) = \mathcal{M}[\mathbf{F}]$. In the finite strain regime, the first Piola-Kirchhoff stress tensor $\mathbf{P}$ and the deformation gradient $\mathbf{F}$ form a pair of thermodynamically conjugate variables. Consequently, any constitutive law of the form $\mathbf{P} = \mathcal{M}(\mathbf{F})$, whether derived from a traditional potential or learned through data, can be directly embedded into a finite element solve.

Let a sample be loaded with $N_R$ reaction forces $R_i$ with $i \in \{1, \dots, N_R\}$, obtained with external force sensors. The EUCLID framework [15] consists of minimizing a loss function based on the force balance residuals given by a first term $\mathcal{L}_{\text{int}}$ at all free nodes and a second term $\mathcal{L}_{\text{BC}}$ at the controlled boundary conditions. The $N_n$ nodes of the finite element mesh can hence be divided into an inner part $N_{\text{int}}$ and the controlled boundary condition part $N_{\text{BC}}$, with $N_n = N_{\text{int}} + N_{\text{BC}}$.

In case of negligible body forces, the nodal residual can be written by integrating the stresses onto the reference domain $\Omega$, considering tetrahedron elements with linear shape functions:

$$ \mathbf{f}_i = \sum_{e=1}^{N_{el}} (\mathbf{P}^e \nabla \psi_i) V_e, \qquad (13) $$

where $\mathbf{P}^e$ is the first Piola-Kirchhoff stress evaluated in element $e$, $\nabla \psi_i$ is the gradient of the shape function $\psi_a$ for node $a$ within that element computed in the reference configuration, and $V_e$ is the volume of the tetrahedral element in the reference configuration. For linear T4 elements used in this work, the Jacobian determinant is constant within each element and directly equal to $V_e$, so the above discrete formulation is equivalent to the general integral form.

The EUCLID framework proposes a first loss term minimizing the Euclidean norm of the nodal force residuals, thereby enforcing local equilibrium within the specimen. In the present experimental context, each DVC measurement used for the model training corresponds to a discrete loading step acquired at a specific time (or equivalently, at a given actuator displacement). To account for all available load levels chosen for training, the loss is written as a function of the temporal index $t$, which enumerates the different loading steps accessible in the dataset

$$ \mathcal{L}_{\text{int}}(t) = \sum_{i \in \text{int}} \|\mathbf{f}_i(t)\|_2^2, \qquad (14) $$

However, acquisition noise assumed to be white and Gaussian impacts the displacement measurement (thus significantly the deformation) and may induce a bias when using the square of each nodal force residual. The second term expects to balance the reaction forces

$$ \mathcal{L}_{\text{BC}}(t) = \sum_{j=1}^{N_R} \left\| \mathbf{R}_j(t) - \sum_{i \in \text{BC}_j} \mathbf{f}_i(t) \right\|_2^2. \qquad (15) $$

The complete loss consists of a weighted sum of the two terms described above, with $\gamma$ a weighting factor [15]. The loss also sums all $N_t$ temporal contributions such that:

$$ \mathcal{L}_t = \frac{1}{N_t} \sum_t^{N_t} [\mathcal{L}_{\text{int}}(t) + \gamma \mathcal{L}_{\text{BC}}(t)]. \qquad (16) $$

The EUCLID method has been extended with the learning of neural-network constitutive models and is referred to as NN-EUCLID [16]. In the latter, a simple weight factor was chosen: $\gamma = 1$.

The reference PANN model is composed of 2 hidden layers and 8 neurons each. An Adam optimizer [36] with a scheduled decreasing learning rate was used to optimize the model weights over 5000 epochs. The learning rate evolved in steps, starting with a rate of 0.02 over the first 1000 epochs and then changing to 0.01 until 2500 epochs, and finally 0.005. The training and validation datasets were not split into different batches. The code is written with the Tensorflow library and is accessible on a GitHub repository (https://github.com/cjailin/Experimental_PANN-learning).
## 2.4. Digital volume correlation

From the tomographic volumes acquired during the mechanical test, high-resolution kinematic measurements can be extracted. Digital Volume Correlation (DVC) [26,37] is a non-rigid image-registration technique that estimates the full three-dimensional displacement field between a reference scan and one or more deformed scans. DVC has been extensively studied and applied across many materials, including metals, polymers, ceramics, and granular media, and even in low-texture or biomedical contexts [27].

In DVC, the displacement field $\mathbf{u}(x)$ is obtained to minimize the summed squared difference between the reference volume $f$ (deformed by $\mathbf{u}$) and the target volume $g$. For every voxel $x$ in the region of interest $\Omega$, and assuming Gaussian noise with variance $\gamma_f^2$, the cost functional is

$$ \Gamma_{\text{DVC}}(\mathbf{u}; a; b) = \frac{1}{2\gamma_f^2} \sum_{x \in \Omega} \|(1 + a(x))g(\mathbf{x} + \mathbf{u}(x)) + b(x) - f(x)\|^2. \quad (17) $$

Here $a(x)$ and $b(x)$ are respectively local contrast and brightness corrections used to match the intensity change of the volume between the two scans. These changes may be due to X-ray beam variations, intensity changes with the deformation (volume thus density changes), artifacts, etc. The displacement field and the correction fields are then obtained with

$$ \mathbf{u}, a, b = \arg \min_{v, \tilde{a}, \tilde{b}} \Gamma_{\text{DVC}}(v; \tilde{a}; \tilde{b}). \quad (18) $$

Because there are far more unknown displacements than available image-intensity constraints, the problem is ill-posed and requires regularization. In global DVC [38], the displacement $\mathbf{u}(x)$ is projected onto a reduced kinematic basis, typically the shape-function basis of a finite-element mesh, so that

$$ \mathbf{u}(x) = \sum_{i=1}^{N_a} \mathbf{u}_i \psi_i(x), \quad (19) $$

where the $\psi_i$ are nodal shape functions and the $\mathbf{u}_i$ are the unknown coefficients.

This parametrization drastically reduces the number of unknowns, enforces continuity of the displacement field, and provides a simple interface with conventional finite-element models for subsequent mechanical analysis. Note that the same mesh and elements are chosen in the EUCLID approach to compute the equilibrium. For an additional remark, the same mesh is used to regularize $a$ and $b$ the brightness and contrast correction fields.

Minimization of $\Gamma_{\text{DVC}}$ is performed by iterative gradient-based methods, involving successive linearizations and updates to the reference image. Convergence is often accelerated by multiscale strategies: one applies Gaussian smoothing to create a pyramid of images from coarse to fine scales, solving first on the coarsest level before refining. Finally, the residual field

$$ \rho(\mathbf{x}; \tilde{\mathbf{u}}; \tilde{a}; \tilde{b}) = (1 + \tilde{a}(\mathbf{x}))g(\mathbf{x} + \tilde{\mathbf{u}}(\mathbf{x})) + \tilde{b}(\mathbf{x}) - f(\mathbf{x}), \quad (20) $$

provides critical diagnostic information. Patterns in this residual field reveal image noise, artifacts, model violations (e.g., discontinuities across cracks), or convergence issues, and thus serve to validate the quality of the DVC measurement. The computations were performed with the in-house global DIC/DVC software Correli [39].

Studying a hyperelastic 3D-printed material with important displacement presents two main challenges for DVC [27]: (i) the magnitude of the displacements can greatly exceed the characteristic length of the internal texture, and (ii) the periodic patterning (due to fused deposition) causes numerous local minima in the correlation functional. To ensure robust convergence, we therefore employed a two-fold multiscale strategy in both the image domain and the kinematic model.

**Image multiscale pyramid.** We constructed a four-level Gaussian pyramid by sequentially down-sampling the reconstructed volumes by factors of $l = \{8, 4, 2, 1\}$ [27]. At each level, the volume was first smoothed with a 3D Gaussian filter (with a standard deviation of $l$) and then selecting one voxel every $l$ in each direction. DVC was performed at the coarsest resolution level until convergence. The resulting displacement field was then rescaled and used to initialize the next level in the pyramid. Brightness and contrast adjustments, which do not require scaling, were similarly propagated. This coarse-to-fine scheme avoids entrapment in local minima by first capturing the large, low-frequency displacement modes and then progressively refining the high-frequency details.

**Mechanical regularization.** In parallel, we imposed a “mechanical smoothing” prior on the displacement field by minimizing a regularization functional based on the equilibrium gap of a locally homogeneous elastic material [40,41]. The regularization length scale was initially set to a large value, approximately one fourth of the smallest dimension of the volume to enforce a smooth initial displacement field. As the algorithm progressed, this length scale was gradually reduced to zero, progressively giving full weight to the data-based DVC term. This iterative strategy prevents the solution from falling to an nonphysical minimum, especially due to periodic printing patterns, while ultimately recovering the fine-scale kinematics once the global displacement trend is well established. It was chosen to completely remove the regularization not to interfere with the ML model learning.

**Intensity correction regularization.** Brightness and contrast correction fields were regularized during the DVC procedure through a squared Laplacian penalty. This approach penalized strong spatial variations in the brightness and contrast fields, promoting the required spatial smoothness. The regularization length scale was set to 100 voxels, corresponding to twice the element size.
## 2.5. Metrics

The accuracy of the model is evaluated with its loss computation, which must be as close to zero as possible. Plots of the nodal forces also helped to judge the quality of the prediction and evaluate where the errors are located. It is expected to see little nodal forces for the internal nodes and a concentration of forces on the boundary conditions in opposite to the imposed effort.

Different metrics will be employed to better appreciate the results: the Mean Internal Error (MIE) on internal nodes and the Mean Boundary Condition Error (MBCE) on boundary condition nodes. The MIE evaluates the average magnitude of residual forces across all free nodes, reflecting the overall equilibrium error within the specimen:

$$ \text{MIE}(t) = \frac{1}{N_{\text{int}}} \sum_{i \in \Omega_{\text{int}}} \| \mathbf{f}_i(t) \|_2, \qquad (21) $$

where $N_{\text{int}}$ is the number of internal nodes in $\Omega_{\text{int}}$ and $\mathbf{f}_i$ is the nodal force residual at node i.

The MBCE quantifies the average deviation of predicted global tension on the boundary from the actual applied and measured boundary conditions in the z direction, indicating how closely the model enforces the known loading and support constraints:

$$ \text{MBCE}(t) = \frac{1}{N_R} \sum_{k=1}^{N_R} \left| R_k(t) - \sum_{i \in \Omega_{BC,k}} \mathbf{f}_i(t) \cdot \mathbf{z}_i \right|, \qquad (22) $$

where $N_R$ is the number of distinct reaction forces, $R_k$ is the true applied load on the k-th boundary region, and $\Omega_{BC,k}$ is the set of nodes carrying that boundary condition.

# 3. 3D Experimental data

## 3.1. Experimental procedure

Two uniaxial tensile experiments have been performed on two 3D printed thermoplastic polyurethane (FDP TPU 92A) specimens (Figure 2). Both specimens had overall dimensions of 35 × 20×110 mm³ but differed in geometry: Experiment-1 contained multiple side holes designed to break symmetry and generate heterogeneous strain fields, whereas Experiment-2 featured a single large central hole to localize deformation within a simpler yet more severe geometry. These holes also intend to create a rich strain field, with a great diversity of strain for each load step. This printing technique may exhibit anisotropic behavior, the samples are massive and composed of 40 layers in the depth.

The samples are first manufactured without holes, which are then drilled in the bulk material using a column drilling machine. Both specimens were printed using identical process parameters and filament batch, with approximately 75–80% infill porosity, but Experiment-2 was produced and tested seven months after the first campaign. Small material and environmental variations (e.g., humidity, filament aging, machine calibration) are therefore expected. Due to the printing process, the resulting material may exhibit an anisotropic behavior which will not be explicitly accounted for in the various models used in this work.

For Experiment-1, the tensile test was controlled in displacement at 50 µms⁻¹. Four loading steps with a controlled stroke amplitude of 7 mm were done. A 25 kN tension/torsion load cell measured the tension applied force of {0, 280, 590, 863, 1110} N, with an uncertainty standard deviation (estimated during a constant loading phase) at ≈2.5 N. The second test was controlled in force to reach similar loadings. Registered forces on Experiment-2 are {237, 427, 687, 871} N. The two experiments are referred to as Experiment-1 and Experiment-2.

For the two tests, the torsion sensor registered a negligible torque throughout the whole experiments (standard deviation of 0.019 N.m).

At each loading step, the sample was allowed to relax for several minutes until the applied force reached a plateau (as seen with the two plots of Fig. 1), the acquisition was conducted using X-ray tomography. In the following analysis, the associated viscous behavior is neglected.

Each of the 10 acquired scan was composed of 900 projections (200 kV, 150 µA, step (10 ms) and shoot procedure with an average of 5 frames) for 360° rotation, acquired using a binning 2×2 leading to a definition of 1536×1944 pixels. This acquisition allows acquiring the required projections in 12 minutes. The reconstruction was carried out using the North Star Imaging – NSI laboratory CT software, using geometric parameters extracted from a pre-scanned calibration target. For Experiment-1, the reconstructed volumes have a definition of 769×1346×1089 voxels with a voxel size of 30 µm. For Experiment-2, the reconstructed volumes have a definition of 839×1412×935 voxels with a voxel size of 29.5 µm

Because of the conical geometry of the X-ray beam, the presence of radio-opaque top and bottom clamps, and the localized concentration of deformation, the analysis was restricted to a small sub-volume within the specimen where the deformation was prominent and beam quality was highest. The scanned zone is a rectangular parallelepiped about 2 cm height in the center of the sample (see Fig. 2(a-b)).

Fig. 3 presents slices of the reconstruction. The pattern of the fused filaments is clearly visible and creates a strong gray level texture providing strong image gradients, useful for DVC. However, the periodic pattern complicates the convergence, leading to multiple local minima. Various printing defects are visible and highlighted with red circles. Those defects can be seen in the two configurations
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** (a) Force-stroke curve of *Experiment-1* with the 5 acquired scans (red dots). The load drop corresponding to relaxation before each scan acquisition is visible. (b) Force-time curve of *Experiment-1* highlighting the relaxation parts. The gray zones correspond to tomography scanning. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** (a-b) *Respectively* *Experiment-1* and 2 - 3D printed sample with the scanned zone in the red rectangle (top and bottom holes are printed used for the clamping), and (c) picture of the lab-CT experiment, composed of (1) an X-ray scintillator, (2) the *insitu* tensile machine, with the specimen position shown with the dashed rectangle and (3) the X-ray source. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 3.2. DVC Measurement fields

For each experiment, five 3D scans were acquired using X-ray tomography: one reference (unloaded) scan and four scans at successive loading steps, yielding ten scans in total across the two experiments. First, a mesh was built on the reference volume using *Thermo Scientific™ Avizo* meshing software. The mesh was positioned based on the gray-level information such that it fits the edges of the sample. This is mandatory for the computation of the equilibrium in the EUCLID framework. This led to adapted finite-element meshes for each experiment:

* *Experiment-1*, a mesh of 4537 nodes and 20,801 first order tetrahedral elements, with an average element size of 49 voxels (≈1.5 mm).

* *Experiment-2*, a mesh of 2608 nodes and 11,138 first order tetrahedral elements, corresponding to an average element size of 57 voxels (≈1.7 mm).
<figure><img src="image_6.png" /><figcaption>Fig. 3. Slices of the reconstructed volumes, before and after the loading, for the two experiments. Printing defects are shown with the red circles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Experiment-1 - slice (x, y) of the residual field ρ(x) of the DVC procedure for step 4. (left) before registration and (right) after registration.</figcaption></figure>

In both cases, the element size is several times larger than the printed filament diameter (0.5 mm), ensuring that the displacement fields describe a homogenized mechanical response at a mesoscopic scale. The aim of this study is not to focus on the anisotropic behavior at a small scale, but rather on the larger isotropic behavior that results. This level of regularization smooths out anisotropy induced by the layer-by-layer printing process, supporting the isotropic assumption adopted for the constitutive modeling.

Using DVC, each of the deformed scans was registered to its reference configuration, yielding four three-dimensional displacement fields that capture the internal deformation of the specimen at successive loading stages. Fig. 5 and 6 display the four loading steps obtained after DVC for each experiment. In Experiment-1, displacement magnitudes reach 160 voxels in the step 4 which already represents a challenge in DVC techniques. At convergence, the residual fields reached low values, below 5% of the gray level dynamic range at the highest load, indicating effective convergence of the algorithm. A representative (x, y) slice of the residual field is shown in Fig. 4. Initially, the residual field exhibited significant positive and negative structures, reflecting strong mismatches due to uncorrected displacements. After convergence, these patterns were substantially diminished. Localized blue features, aligned with the printed layers, remained visible and may indicate residual inconsistencies between the corrected deformed image and the reference. These artifacts are hypothesized to originate from compressed porosities, i.e., gaps between printed layers, that are not adequately resolved by the relatively large DVC elements used in the analysis.

### 3.3. Data analysis

#### 3.3.0.1. Non-linearity in the measured fields.
To assess the non-linearity of the material response, a singular value decomposition (SVD) was applied to the space-time displacement gradient of Experiment-1 $\nabla u(x, t)$. For this, each tensor component was discretized and reshaped into a matrix $M \in \mathbb{R}^{4N_{el} \times N_t}$, where each row corresponds to a spatial position x and each column to a time step t. The decomposition reads:

$$M(x, t) = \sum_{m} A_m \Xi_m(x) \Theta_m^T(t), \qquad (23)$$

where $A_m$ are the singular values, $\Xi_m$ are spatial components (space singular vectors), and $\Theta_m$ are temporal components (time singular vectors). These satisfy the orthonormality conditions:

$$\Xi_m^T \Xi_n = \delta_{mn}, \quad \Theta_m^T \Theta_n = \delta_{mn}, \qquad (24)$$

with $\delta_{mn}$ the Kronecker delta, and singular values ordering: $A_m \ge A_{m+1}$. Fig. 7 presents the three most significant spatial components, out of the four available, associated with singular values A: {80.71, 3.35, 1.80}.

The dominant mode ($A_1 = 80.71$) captures the primary tensile loading. The presence of additional modes, despite their much smaller singular values, signals non-linear (with respect to (u) contributions to the overall response. In a purely linear elastic law in
<figure><img src="image_8.png" /><figcaption>Fig. 5. *Experiment-1* - Visualization of the 4 displacement fields (with the rigid body motion removed).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. *Experiment-2* - Visualization of the 4 displacement fields (with the rigid body motion removed). The colormap gives a visual indication of the displacement magnitude variation.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. *Experiment-1* - The first three deformation components of the singular value decomposition Ξ₁, Ξ₂, Ξ₃. For the display, only the *zz* deformation is represented (+ 1 to visually match with the deformation tensors).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 8. *Experiment-1* - K₁, K₂ and K₃ invariants for the 4 datasets. Note that loading 4 is used as training data and the three other loading as evaluation/test data.</figcaption></figure>

small strain, only a single mode would carry significant weight. Here, the non-negligible A₂ and A₃ confirm that the material behavior escapes from simple first gradient linearity.

All nodal coordinates, displacements, contrast, and brightness nodal corrections (u, a, b), force measurement data and connectivity matrix are available on a GitHub repository (https://github.com/cjailin/Experimental_PANN-learning).

3.3.0.2. *Phase space analysis.* Following an interesting data representation from [42], which characterizes the heterogeneous kinematic field, the Hencky logarithmic strain tensor, also known as the true strain tensor, is defined as: **H** = 1/2 log (FFT). Criscione et al. [43] introduced three relevant invariants of **H**, denoted as K₁, K₂, and K₃. They are defined as follows:

$$K_1 = \text{tr}(\mathbf{H}), \quad K_1 \in \mathbb{R}, \qquad (25)$$

$$K_2 = \sqrt{\text{dev}(\mathbf{H}) : \text{dev}(\mathbf{H})}, \quad K_2 > 0, \qquad (26)$$

$$K_3 = \frac{3\sqrt{6}}{K_2^3} \det(\text{dev}(\mathbf{H})), \quad K_3 \in [-1, 1], \qquad (27)$$

These invariants provide a characterization of different aspects of the deformation field: the amount of dilatation, the magnitude of distortion, and the mode of distortion, respectively. K₂ is a real positive number that increases with the magnitude of strain, while K₃ is a real number that describes the nature of the deformation state. The particular values of K₃ and their corresponding deformation states, as described in [42], are K₃ = 1 for uniaxial tension or equibiaxial compression; K₃ = 0 for planar tension (pure shear); and K₃ = -1 for uniaxial compression or equibiaxial tension.

The analysis of the invariants obtained from the DVC fields is presented in Fig. 8 and 9, respectively for the first and second experiment. The two datasets are primarily characterized by uniaxial tension or equibiaxial compression, as expected in tensile tests. However, a non-negligible number of elements also exhibit pure shear (with K₃ ≈ 0) and uniaxial compression or equibiaxial tension (with K₃ ≈ -1). This limited but notable complexity in the data arises from the presence of holes, which are intended to induce multiaxial deformation states.

The comparison between the two datasets shows that *Experiment-1* explores a richer phase space of deformation modes, while *Experiment-2* remains dominated by uniaxial tension around the central hole. This difference is deliberate: *Experiment-1* was designed for training diversity, whereas *Experiment-2* was designed for validation robustness under new geometric and experimental conditions.
<figure><img src="image_12.png" /><figcaption>Fig. 9. *Experiment-2* - K₁, K₂ and K₃ invariants for the 4 datasets. Note that all loadings of *Experiment-2* are used for the evaluation/test data.</figcaption></figure>

## 4. Results

The overall dataset therefore comprises eight distinct loading steps: four from *Experiment-1* and four from *Experiment-2*. The fourth loading step of *Experiment-1* was used for training, while the remaining seven loading steps serve as unseen validation or test configurations. In particular, *Experiment-2* provides a fully independent validation campaign, conducted months apart and involving a different geometry. Because *Experiment-2* explores a lower-dimensional deformation manifold than *Experiment-1*, generalization could be anticipated as the mechanics excited in *Experiment-2* remains included in the training manifold learned from *Experiment-1*. This setting therefore constitutes a meaningful yet tractable validation of model transfer under changes of geometry and acquisition context.

Before analyzing the quantitative results, the main modeling assumptions are summarized for clarity. The 3D-printed TPU specimen, although manufactured by fused deposition, is treated as macroscopically isotropic and homogeneous at the scale of the DVC finite elements. The potential anisotropy induced by the layered printing process and filament orientation is therefore neglected, as the relatively coarse volumetric mesh (element size ≈1.5 mm) implicitly performs spatial averaging.

The material is assumed to behave in a purely elastic and rate-independent manner. Minor time-dependent effects observed during the load-hold phases (see Fig. 1b) are interpreted as small viscoelastic relaxations and are neglected in the present formulation. Temperature variations during scanning are considered negligible.

Finally, all boundary conditions are directly imposed from the measured displacements and global reaction forces, without additional symmetry or regularization constraints. These simplifying assumptions define the framework within which the equilibrium-based learning (EUCLID loss) and subsequent model comparisons are conducted.

### 4.1. Training phase

The finite-element mesh was partitioned into boundary nodes, where traction or displacement conditions are applied, and free (internal) nodes. In particular, nodes with minimal z-coordinate define the bottom surface (Fig. 10). Boundary nodes are shown in red, while free nodes appear in black.

The constitutive response is first approximated using the SVK and NH models. For each model, a grid search over the two elastic parameters was performed, and the EUCLID loss on the fourth load step was evaluated at every grid point. The resulting loss-landscapes, expressed in log decimal scale, are shown in Fig. 11.

The two heatmaps reveal a high sensitivity to the Young's modulus and a relatively low sensitivity to the Poisson's ratio, as indicated by the steep and steady gradients of the loss surfaces with respect to these parameters, respectively. This disparity in
<figure><img src="image_13.png" /><figcaption>Fig. 10. Boundary nodes (red for the top nodes, green for the bottom nodes) and free (internal) nodes (black) for Experiment-1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. EUCLID loss as a function of material parameters for (left) Saint-Venant-Kirchhoff and (right) Neo-Hookean models on the fourth load step.</figcaption></figure>

sensitivity is a well-known characteristic of DVC-based identification procedures [44]. The NH model yielded a lower minimum error compared to the linear model. The optimal model parameters were identified by selecting the parameter set corresponding to the minimum value on the loss map. For the SVK model, the minimum loss of $L_{SVK}$ = 11070 N² occurs at

$$E_{SVK} = 4.86 \text{ MPa}, \quad v_{SVK} = 0.014. \qquad (28)$$

The Neo-Hookean model achieves a lower minimum, $L_{NH}$ = 4216 N², at

$$E_{NH} = 10.29 \text{ MPa}, \quad v_{NH} = 0.286. \qquad (29)$$

The PANN was trained for 5000 epochs on the fourth load step using the same EUCLID loss. The convergence curve is shown in black on Fig. 12.

The final PANN training loss, represented by the black plain line is $L_{train}^{PANN}$ = 2526 N², substantially lower than for SVK and NH models. The PANN reaches low values very quickly in approximately 200 epochs while being initialized randomly. Then the evolution is slow and decreases.

Fig. 13 displays the nodal forces after learning on the training step with the three models. These vector plots are to be read in two steps.

First, a concentration of vertical forces is visible on the top and bottom surfaces. The sum of those vectors on each of those surfaces is expected to compensate for the measured load. The mean absolute error modeling the applied force MBCE is, respectively, for the SVK, NH, and PANN model, {51.2, 8.4, 4.8} N.

Then, the non-zero internal nodal forces illustrate the internal EUCLID Loss value. In the center of the sample, non-zero internal nodal forces are visible near the holes and on the free surfaces. These zones are presenting the most important strain and probably the main non-linearity. More important internal forces are visible for the Saint-Venant-Kirchhoff prediction than for the two other models. The MIE metric is, respectively, for the SVK, NH, and PANN models, {1.16, 1.02, 0.80}.

Both the NH and PANN models show similar results in the force plot, although a difference appears at the left and right edges, where the NH model exhibits stronger forces on the free surface. The summary of the losses and metrics in the training step is proposed in the Table 1.
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Convergence curves during the PANN training, monitored with the three validation steps of *Experiment-1*. In the legend, NH-i refers to the Neo-Hookean results for each scan.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Nodal forces magnitudes on the fourth load step (R = 1110 N), used for training, for (left) SVK model, (center) NH model, and (right) the 3D-PANN using the identified parameters.</figcaption></figure>

<table><caption>**Table 1**

EUCLID losses for each constitutive model evaluated on the training set.</caption><thead><tr><th>Model</th><th><i>E</i> [MPa]</th><th><i>ν</i></th><th>Loss [N<sup>2</sup>]</th><th>Loss/R<sup>2</sup> [%]</th><th>MBCE [N]</th><th>MIE [N]</th></tr></thead><tbody><tr><td>SVK</td><td>4.86</td><td>0.014</td><td>11070</td><td>8.98</td><td>51.2</td><td>1.16</td></tr><tr><td>NH</td><td>10.29</td><td>0.286</td><td>4216</td><td>3.42</td><td>8.9</td><td>1.02</td></tr><tr><td>3D-PANN</td><td>-</td><td>-</td><td>2526</td><td>2.05</td><td>4.8</td><td>0.80</td></tr></tbody></table>

## 4.2. Generalization to unseen test steps

To assess predictive capability, which is the true objective of the method, the learned models were applied to the first three loads of data.

The nodal force distributions are presented for *Experiment-1* (Fig. 14), each time for the three identified behaviors.

In all three loading steps, the internal errors are lower for the PANN model. The SVK model has higher internal forces, especially around the edges. The summary of all the losses and metrics in the test steps is proposed in the Table 2.

Throughout the training process based on the data of the fourth loading step, inference was performed at each epoch for the three validation datasets. This results and the losses corresponding to the optimal NH model are presented on Fig. 12 and allows to assess how the identified PANN model behaves with respect to the validation datasets and the classical parametric models. Here, the average validation loss decreases overall, indicating that the model is learning relevant features during training. However, Validation 2 and 3 exhibit relatively low loss values (same order of magnitude as the training loss) at early epochs, followed by a gradual increase, suggesting potential overfitting or a shift in generalization behavior. Validation 1, surprisingly, reaches a very low error within the
<figure><img src="image_17.png" /><figcaption>Fig. 14. Nodal forces magnitudes for (left) SVK model, (center) NH model, and (right) the 3D-PANN using the identified parameters.</figcaption></figure>

<caption>**Table 2**</caption>

Generalization results on the three test steps of *Experiment-1*, for the three materials models.

<table><thead><tr><th>Load</th><th>Model</th><th>Loss [N<sup>2</sup>]</th><th>Loss/R<sup>2</sup> [%o]</th><th>MBCE [N]</th><th>MIE [N]</th></tr></thead><tbody><tr><td rowspan="3">Load 3</td><td>SVK</td><td>34231</td><td>45.9</td><td>126.9</td><td>0.71</td></tr><tr><td>NH</td><td>3101</td><td>4.16</td><td>18.4</td><td>0.76</td></tr><tr><td>3D-PANN</td><td>2743</td><td>3.68</td><td>21.0</td><td>0.62</td></tr><tr><td rowspan="3">Load 2</td><td>SVK</td><td>45464</td><td>130</td><td>149.6</td><td>0.43</td></tr><tr><td>NH</td><td>1432</td><td>4.11</td><td>12.8</td><td>0.52</td></tr><tr><td>3D-PANN</td><td>1898</td><td>5.45</td><td>20.9</td><td>0.46</td></tr><tr><td rowspan="3">Load 1</td><td>SVK</td><td>10077</td><td>128</td><td>70.1</td><td>0.21</td></tr><tr><td>NH</td><td>10550</td><td>135</td><td>70.6</td><td>0.28</td></tr><tr><td>3D-PANN</td><td>4241</td><td>54.1</td><td>42.8</td><td>0.28</td></tr></tbody></table>

first few epochs, but then shows a significant increase, up to 6000 before decreasing again toward the end of training. This transient instability in Validation 1 is also observed in the baseline models, both NH and SVK, indicating that this particular loading step may be more difficult to generalize.

For Validation 1 and the Training set, the identified PANN model provides solution with lower losses compared to the optimal NH model after few epochs (below 100). It is not the case for the Validation 2 and Validation 3 sets where the NH loss is close the PANN during training and at convergence Table 2 at. This suggests that the model may have learned step-specific patterns more efficiently for step 1, while still matching baseline performance for the others.

The comparison of predicted top and bottom reaction forces across different models and simulation steps, as shown in Fig. 15, highlights the stronger performance of the proposed PANN model.
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Experiment-1 - Absolute predicted top and bottom forces for all models across four simulation steps, with ground truth forces shown as dashed black lines.</figcaption></figure>

<caption>**Table 3**</caption>

Generalization results on the four test steps of *Experiment-2*, for the three materials models.

<table><thead><tr><th>Load</th><th>Model</th><th>Loss [N²]</th><th>Loss/R² [%o]</th><th>MBCE [N]</th><th>MIE [N]</th></tr></thead><tbody><tr><td rowspan="3">Step 4</td><td>SVK</td><td>255546</td><td>336.8</td><td>356.56</td><td>0.790</td></tr><tr><td>NH</td><td>6006</td><td>7.9</td><td>46.42</td><td>0.920</td></tr><tr><td>PANN</td><td>3111</td><td>4.1</td><td>28.40</td><td>0.850</td></tr><tr><td rowspan="3">Step 3</td><td>SVK</td><td>167591</td><td>355.1</td><td>288.85</td><td>0.590</td></tr><tr><td>NH</td><td>1526</td><td>3.2</td><td>16.60</td><td>0.690</td></tr><tr><td>PANN</td><td>991</td><td>2.1</td><td>3.69</td><td>0.690</td></tr><tr><td rowspan="3">Step 2</td><td>SVK</td><td>69683</td><td>382.2</td><td>186.31</td><td>0.350</td></tr><tr><td>NH</td><td>622</td><td>3.4</td><td>11.28</td><td>0.420</td></tr><tr><td>PANN</td><td>1286</td><td>7.1</td><td>20.72</td><td>0.450</td></tr><tr><td rowspan="3">Step 1</td><td>SVK</td><td>21277</td><td>378.8</td><td>102.98</td><td>0.180</td></tr><tr><td>NH</td><td>150</td><td>2.7</td><td>5.28</td><td>0.210</td></tr><tr><td>PANN</td><td>675</td><td>12.0</td><td>16.57</td><td>0.240</td></tr></tbody></table>

Among all evaluated methods, PANN consistently delivers predictions closest to the ground truth values at all levels of loading. The model also preserves symmetry between top and bottom forces, indicative of a physically consistent response. In contrast, the SVK model shows significant deviations from the expected values, particularly at higher load levels. The NH model performs better than SVK but is less accurate than PANN, with moderate discrepancies in force magnitude. These results confirm the effectiveness of the PANN approach in capturing complex mechanical behaviors with higher fidelity and stability, particularly in data-driven simulations of nonlinear elasticity.

Beyond intra-experiment generalization, the model ability to extrapolate across experiments is evaluated using the independent datasets of *Experiment-2*. This test is more demanding, as it combines differences in geometry, local boundary conditions, micro-structure, and acquisition environment. Such cross-experiment validation represents a realistic assessment of the reproducibility of data-driven (and also traditional) constitutive models. The obtained results are presented in Table 3.

The residual-force distributions are presented for completely unseen second experiment (Fig. 16), each time for the three identified behaviors.

Despite being trained solely on *Experiment-1*, the PANN maintained satisfactory equilibrium on the second specimen, outperforming or remaining comparable to the Neo-Hookean baseline. Interestingly, the internal equilibrium error (MIE) remains relatively low across all load steps, indicating that the overall stress-strain relationship and the shape of the constitutive behavior have been properly captured. Such deviations could arise from differences in the effective stiffness between the two experiments, possibly caused by minor material aging or humidity absorption in the TPU [45]. Additional factors, such as small boundary alignment differences or voxel-size variations in the DVC reconstructions, could further contribute to this amplitude mismatch. These observations highlight
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** *Experiment-2* - Nodal forces magnitudes for (left) SVK model, (center) NH model, and (right) the 3D-PANN using the identified parameters from the fourth load of *Experiment-1*.</figcaption></figure>

both the capacity of the model to generalize the constitutive form and the sensitivity of the absolute stress scaling to physical and experimental variations.

### 4.3. Impact of the PANN architecture

The number of layers and the number of neurons in the hidden layers of the PANN were tested. The computations were performed using the same datasets and similar optimization parameters. The tested number of layers included 1, 2, 3, 4 and 5, while the number of neurons per layer was assessed at 2, 4, 8, and 16. In total, 20 computations were evaluated for all combinations. The extracted metric is the sum of the three losses from the three test loading steps, without including step 4, which would highly bias the results as step 4 may reduce as the number of parameters increases.

### 4.4. Impact of 3D data measurement

To assess the importance of incorporating full 3D data into learning procedures, the impact of data dimensionality on model performance is investigated. To emulate a surface Digital Image Correlation (DIC) measurement setup, only the in-plane displacement components (x, y, z) on the front surface of the 3D dataset (with surface normal along y) are used. At this stage, a plane strain assumption
<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Results of the training with various model sizes. The first number *p* corresponds to the number of trainable parameters in the network. The second value corresponds to the average loss on the three (unseen during training) loading steps of *Experiment-1*.</figcaption></figure>

<table><caption>**Table 4**
Results of the 2D-PANN training in *Experiment-1*.</caption><thead><tr><th>Load</th><th>Model</th><th>Loss [N²]</th><th>MBCE [N]</th><th>MIE [N]</th></tr></thead><tbody><tr><td>Load 4</td><td>2D-PANN</td><td>40405</td><td>126.3</td><td>2.52</td></tr><tr><td>Load 3</td><td>2D-PANN</td><td>31510</td><td>111.7</td><td>2.39</td></tr><tr><td>Load 2</td><td>2D-PANN</td><td>14961</td><td>76.9</td><td>1.65</td></tr><tr><td>Load 1</td><td>2D-PANN</td><td>2207</td><td>25.0</td><td>0.62</td></tr></tbody></table>

is chosen, which is common for thick samples where through-thickness deformation is negligible. Some assumption on the displacement field across the thickness is necessary in 2D, with common options including plane strain, plane stress, or incompressibility.

It is important to note that the plane strain assumption is applied only during the training (step 4) and the trained model is referred to as 2D-PANN. During evaluation, the full 3D ground truth displacement field is used. The model training process is conducted with a loss function considering only the *x* and *z* force components, deliberately excluding the *y*-direction forces used in the fully 3D approach. This selective approach mitigates potential errors introduced by the plane strain simplification. Final model performance evaluation uses the full loss function, consistent with prior analyses, utilizing all available displacement and force data.

Training follows an identical methodology, with outcomes presented on the first experiment, presented in Table 4.

# 5. Discussion

## 5.1. Training phase

The almost 0-value identified Poisson's ratio for the linear-elastic (SVK) model is non-physical for such materials. The one identified with the NH model is slightly below to the expected value for printed TPU, which ranges from $v \approx 0.35$ to $v \approx 0.4$ [46,47]. This underestimation may come from two sources: (i) the mismatch between the model assumption and the true hyperelastic response of the pore-filled specimen, and (ii) the important 20–25% porosity introduced by the additive-manufacturing process. The Poisson's ratio estimated with the SVK model is completely incorrect.

The SVK behavior is not well adapted to the hyperelastic behavior of the TPU, yet the minimization tries to find a compromise. Because it is trained on a specimen that is already highly loaded, it underestimates the Young's modulus compared to the NH model. In contrast, the NH model, with its nonlinear strain-energy density, better captures the hyperelastic behavior and thus attains a significantly lower EUCLID loss.

The PANN further reduces unbalanced residual forces, achieving the closest approach to mechanical equilibrium among all three models. The mean reaction force error is only 4.8 N (i.e., corresponding to 0.43% of error and twice the sensor standard deviation uncertainty). This superior performance is consistent with its higher expressivity: the PANN contains 121 trainable parameters versus only two for the SVK and NH models. Consequently, it can adapt to subtle deviations from idealized constitutive laws. Ultimately, however, the true test of each model lies in its generalization to unseen loading states. In the following section, we compare their performance on the validation load step to assess predictive capability.
<figure><img src="image_21.png" /><figcaption>Fig. 18. Uniaxial response of the three models.</figcaption></figure>

The responses of the models for a simple uniaxial synthetic loading are presented in Fig. 18. For this computation, a deformation tensor $F_{\text{uniaxial}}$, composed of $\alpha$ and $\beta$ amplitudes, was optimized with gradient descent so that the radial components of the stresses vanish.

$$F_{\text{uniaxial}} = \begin{bmatrix} \alpha & 0 & 0 \\ 0 & \alpha & 0 \\ 0 & 0 & \beta \end{bmatrix} \qquad (30)$$

The optimization was performed for each loading $\beta$ and each model.

The linear model underestimates the Young's modulus required to capture the non-linear response, in contrast to the NH and PANN models. In this uniaxial synthetic test, the NH and PANN models deliver close response, with a slight underestimation of the PANN model. Most of the differences between the NH and PANN models appear at large deformations, consistently with the trends observed in the quantitative metrics. This indicates that the PANN more accurately captures the material response in the high-strain regime.

## 5.2. Generalization to unseen loading steps

The capacity of each constitutive model to predict material response under unseen loading conditions was evaluated by inferring on the first three load steps. In every case, the PANN demonstrated superior generalization performance compared to both the SVK and NH benchmarks. Quantitatively, the average equilibrium error (MIE) and boundary-condition deviation (MBCE) increased only modestly for the PANN when moving from the training step (step 4) to the validation steps (steps 1-3), whereas the classical models exhibited substantially larger degradations (Table 2).

Visually, the residual-force maps on steps 3, 2 and 1 reveal that the PANN maintains a uniformly low level of internal nodal forces across the perforated specimen, even in regions of high strain concentration around the holes (Figs. 14). In contrast, the SVK model produces larger force imbalances at the free surfaces and near geometric singularities, indicating a failure to extrapolate its small-strain assumption to intermediate load levels. Although the NH model captures some nonlinearity and reduces these local residuals relative to SVK, it still accumulates larger errors than the PANN, particularly at lower applied loads where the tension-strain relationship deviates from the simple Neo-Hookean form.

This robustness of the PANN arises from several factors. First, the network's expressivity, composed of 121 trainable parameters, permits the learning of subtle constitutive features that extend beyond the monotonic stress-strain curve obtained during training. Second, the thermodynamic constraints ensures a smooth learning and prevent the network from converging to nonphysical solution. Third, the inclusion of multiple invariants at the input layer ensures that the network can represent a wide range of deformation modes (dilatation, distortion and their couplings), which naturally enhances its ability to interpolate and mildly extrapolate across the observed phase space. Although at small scales the material may exhibit anisotropy due to the printing process and mesostructure, the chosen DVC element size implicitly introduces a homogenization effect. As such, the measured behavior may be interpreted as isotropic at the scale of the mesh.

Furthermore, the PANN's ability to locally adapt its constitutive response allows it to handle measurement noise and data irregularities in a fundamentally different manner from classical models. During DVC, localized artefacts and Gaussian noise are inevitably present in the reconstructed displacement fields. Traditional analytical models, whose constitutive mappings are globally defined and smooth, tend to propagate such perturbations uniformly across the domain, which limits their sensitivity to local outliers but also constrains their flexibility. In contrast, the PANN can modify its response through learning, allowing it to accommodate small fluctuations in the measured fields without creating large force imbalances. However, this same adaptability can also make it more sensitive to structured or correlated noise, which may bias the identified constitutive behavior if not properly regularized. In that sense, traditional models are robust but restrictive, while PANNs are flexible but potentially sensitive to data quality. The impact of
<figure><img src="image_22.png" /><figcaption>Fig. 19. Internal residual force differences between the NH and PANN model.</figcaption></figure>

experimental noise on unsupervised constitutive learning therefore remains an important topic for further investigation, and advanced formulations explicitly accounting for measurement uncertainty [48] could improve robustness and interpretability.

The PANN's generalization performance, evidenced by lower MIE and MBCE across unseen load steps, demonstrates that machine learning models can learn a truly transferable representation of hyperelastic behavior from limited data. The difference between the PANN and NH residual forces at step 3 is shown in Fig. 19. The majority of the force differences are identified on the free surfaces. This is due to the loss that minimizes internal squared residual forces and finds minimal solution when edge nodes have residuals due to less neighbouring nodes. Note that the forces at the BC are very similar between PANN and NH.

Furthermore, it seems that the mesh could have been refined near the holes. Indeed, the stress computed around geometrical singularities changes significantly from one element to another. Thus, it suggests that a more refined mesh could enable for a better prediction. This can also be an explanation for the high force values near the holes.

### 5.3. Architecture of the PANN model

In this application, the architecture (specifically, the number of neurons and the number of hidden layers) does not appear to have a significant impact (see Fig. 17). With a very small number of trainable parameters (e.g., with only 2 neurons per layer), the model lacks the capacity to learn effectively. However, with parameter increases, the model acquires the expressiveness needed to capture the nonlinear behaviors in the training configuration. No strong differences are evident among the models. This might be attributed to the simplicity of the experiments, which does not show strong and complex nonlinear patterns and is hence not complex to model as it does not probe a large deformation space.

### 5.4. 3D Data impact

Using only surface measurements under the plane strain hypothesis notably reduces data richness and severely compromises model performance. Large increases in both training and testing losses are evident when compared to previous fully 3D analyses (see Table 4). The important 2D-trained loss underscores the inadequacy of capturing complex internal displacement fields through simplified plane strain assumptions. Given the complex internal architecture of the specimen, characterized by holes and significant thickness, the plane strain approximation proves inadequate, emphasizing the necessity for comprehensive 3D data in such scenarios.

To characterize the anisotropy of the strain state, a singular value decomposition was performed on the six independent components of the Green-Lagrange strain tensor *E* across all nodes and loading steps (see Fig. 20). Note that compared to Section 3.3.0.1, the SVD is here performed on the space of {nodal information, 6 components of *E*} and not on the space {nodal information, time}. The first singular vector accounted for over 96% of the total variance, indicating that the deformation is mostly uniaxial, as expected from the experimental setup. This dominant mode was primarily composed of the axial strain component *E*<sub>22</sub>, consistent with the direction of applied loading. Additional components show that the deformations are 3D (although associated with small variance) and really require a 3D data measurement technique. Two transverse components, *E*<sub>00</sub> and *E*<sub>11</sub>, contributed unequally to this mode, with *E*<sub>00</sub> exhibiting a noticeably higher weight than *E*<sub>11</sub>. This asymmetry in the transverse response suggests a slight deviation from ideal isotropy, potentially arising from material heterogeneity, asymmetric boundary conditions, or geometrical effects. This anisotropy is not currently taken into account in the present formulation of the model.

### 5.5. Limitations and perspectives

Several limitations of the present study must be acknowledged, as they outline promising avenues for future work.

#### 5.5.0.1. Dataset coverage and generalization.
The training dataset was limited to a single monotonic tensile loading path composed of only four load steps. Although the PANN demonstrated satisfactory generalization to lower-magnitude loading conditions, its behavior under more complex scenarios, such as multiaxial or non-monotonic loading remains unexplored. Demonstrating generalization capability is inherently dependent on the choice of test data [49]. The design of experiments involving combined tension, compression, and shear through multi-axial loading paths would enable the PANN to learn a more complete constitutive response manifold. Nevertheless, even such efforts may yield a limited phase space in terms of variability and richness. In addition, an interesting opportunity lies in learning the variability between *Experiment-1* and *Experiment-2*. While both tests follow the same general constitutive
<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Results of the SVD decomposition of the {Green Lagrange component, nodes} space. (left) Singular vectors and (right) explained variance ratio. The first component represents 96% of the total variance.</figcaption></figure>

behavior, they exhibit small yet physically meaningful differences associated with factors such as temperature, humidity, filament aging, and effective stiffness. Capturing these variations within a unified learning framework [50] could enable the identification of a global model with shared constitutive structure and experiment-specific corrections, thereby improving both interpretability and robustness.

In contrast to 2D-DIC, which benefits from high temporal resolution, DVC, due to the long acquisition procedure, is typically constrained to a small number of scans, often between five to ten. This limitation restricts the temporal richness of the training data. Employing faster acquisition strategies, such as projection-based DVC [51], and it's coupling with DVC, could significantly enhance the temporal mechanical coverage of the dataset, thereby improving learning efficacy and potentially identify the viscous relaxation effects.

### 5.5.0.2. Spatial resolution and microstructural effects.
The finite element mesh used for both DVC regularization and EUCLID equilibrium enforcement was of moderate resolution and approximately uniform. Consequently, the DVC measurements were conducted at a macroscopic scale, without resolving the meso- or microstructural features of the printed layers. This led to a homogenized displacement field, which may obscure anisotropic effects or localized phenomena. While a finer DVC resolution could reveal these behaviors, doing so requires advanced correlation methods, as the filament microstructure often exhibits low signal-to-noise ratio. In addition, coarse meshes may fail to capture sharp strain gradients near geometric singularities, such as holes or sharp boundaries resulting in localized residual forces. Addressing this limitation could involve adaptive meshing strategies or the integration of multiscale formulations. Incorporating transverse anisotropy into the PANN framework also presents a relevant perspective [12,18,52]. The use of rich 3D experimental data will facilitate these learning. In addition, strategies to learn and compensate for experimental biases, such as material behavior variability across experiment or acquisition artifacts [50] could be integrated to improve model robustness and predictive accuracy.

### 5.5.0.3. Loss formulation and optimization.
Finally, the current EUCLID training procedure employs a fixed weighting factor $\gamma = 1$ between the internal equilibrium term and the boundary traction term in the loss function (see Eq. 16). However, the relative importance of these two components may vary depending on the load level, material behavior, or noise characteristics of the data. Modifying the loss term [18,21,22], or incorporating adaptive weighting strategies that account for data quality, boundary uncertainty, and local noise characteristics [48,53] could further enhance convergence stability and improve the balance between enforcing interior equilibrium and boundary condition fidelity.

## 6. Conclusion

This study presents, to the best of the authors' knowledge, the first fully unsupervised learning of a hyperelastic constitutive law directly from three-dimensional experimental data using a Physics-Augmented Neural Network (PANN). A rich experimental protocol, combining mechanically regularized Digital Volume Correlation (DVC) and micro-computed tomography, enabled the learning of material behavior from volumetric displacement fields and global reaction forces, without any measured internal stress data. The proposed framework integrates polyconvexity constraints through an Input-Convex Neural Network and enforces equilibrium using the EUCLID loss. Compared to traditional constitutive models such as Saint-Venant-Kirchhoff and Neo-Hookean laws, the 3D-PANN
demonstrated better performance in terms of equilibrium satisfaction and generalization to unseen loading steps. This confirms the ability of the network to learn physically consistent and transferable representations of hyperelastic behavior from limited data.

The transition from surface-based 2D measurements to volumetric 3D fields is an important advance. It not only removes restrictive modeling assumptions such as plane-stress or incompressibility but also significantly improves equilibrium enforcement by leveraging richer measured boundary conditions. Free surfaces, volumetric gradients, and rich spatial measurements contribute to a better-posed inverse problem and more accurate learning outcomes. This approach offers a path toward more realistic and robust data-driven constitutive modeling, particularly in complex material/structure couplings. Future investigations may explore multiaxial loadings, and anisotropic formulations to further exploit the potential of volumetric experimental data in combination with machine learning.

## Code and data availability

The global DVC datasets and the scripts used to train and validate the SVK, NH, and PANN models for the current study are publicly available in the GitHub repository https://github.com/cjailin/Experimental_PANN-learning.

## Declaration of generative AI and AI-assisted technologies in the writing process

During the writing process of this work the authors used *Grammarly* and *ChatGPT 4o* in order to proofread sentences of the manuscript and correct English language errors. After using these tools, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

## CRediT authorship contribution statement

**M. Bourdyot:** Writing – original draft, Investigation, Formal analysis, Data curation; **M. Compans:** Writing – original draft, Investigation, Formal analysis, Data curation; **R. Langlois:** Writing – original draft, Methodology, Formal analysis; **B. Smaniotto:** Writing – original draft, Visualization, Investigation, Data curation; **E. Baranger:** Writing – original draft, Methodology, Formal analysis; **C. Jailin:** Writing – original draft, Supervision, Project administration, Methodology, Conceptualization.

## Data availability

Data and code are shared in a Github repository

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Jailin reports financial support was provided by French National Research Agency. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This project was made possible by the *Agence Nationale de la Recherche* (ANR) grant No. ANR-22-CPJ2-0046-01.

## References

[1] H. Jin, E. Zhang, H.D. Espinosa, Recent advances and applications of machine learning in experimental solid mechanics: a review, Appl. Mech. Rev. 75 (6) (2023) 061001.

[2] J.N. Fuhg, G. Anantha Padmanabha, N. Bouklas, B. Bahmani, W. Sun, N.N. Vlassis, M. Flaschel, P. Carrara, L. De Lorenzis, A review on data-driven constitutive laws for solids, Arch. Comput. Methods Eng. (2024) 1–43.

[3] A. Hussain, A.H. Sakhaei, M. Shafiee, Machine learning-based constitutive modelling for material non-linearity: a review, Mech. Adv. Mater. Struct. (2024) 1–19.

[4] L.A. Herrmann, S. Kollmannsberger, Deep learning in computational mechanics: a review, Comput. Mech. (2024).

[5] M. Raissi, P. Perdikaris, G.E. Karniadakis, Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378 (2019) 686–707.

[6] D.W. Abueidda, S. Koric, E. Guleryuz, N.A. Sobh, Enhanced physics-informed neural networks for hyperelasticity, Int. J. Numer. Methods Eng. 124 (7) (2023) 1585–1601.

[7] F. Masi, I. Stefanou, P. Vannucci, V. Maffi-Berthier, Thermodynamics-based artificial neural networks for constitutive modeling, J. Mech. Phys. Solids 147 (2021) 104277.

[8] F. Masi, I. Stefanou, Multiscale modeling of inelastic materials with thermodynamics-based artificial neural networks (TANN), Comput. Methods Appl. Mech. Eng. 398 (2022) 115190.

[9] M. Fernández, M. Jamshidian, T. Böhlke, K. Kersting, O. Weeger, Anisotropic hyperelastic constitutive models for finite deformations combining material theory and data-driven approaches with application to cubic lattice metamaterials, Comput. Mech. 67 (2) (2021) 653–677.

[10] D.K. Klein, M. Fernández, R.J. Martin, P. Neff, O. Weeger, Polyconvex anisotropic hyperelasticity with neural networks, J. Mech. Phys. Solids 159 (2022) 104703.

[11] F. As'ad, P. Avery, C. Farhat, A mechanics-informed artificial neural network approach in data-driven constitutive modeling, Int. J. Numer. Methods Eng. 123 (12) (2022) 2738–2759.

[12] L. Linden, D.K. Klein, K.A. Kalina, J. Brummund, O. Weeger, M. Kästner, Neural networks meet hyperelasticity: a guide to enforcing physics, J. Mech. Phys. Solids (2023) 105363.
[13] K. Linka, E. Kuhl, A new family of constitutive artificial neural networks towards automated model discovery, Comput. Methods Appl. Mech. Eng. 403 (2023) 115731.

[14] B. Amos, L. Xu, J.Z. Kolter, Input convex neural networks, in: International Conference on Machine Learning, PMLR, 2017, pp. 146–155.

[15] M. Flaschel, S. Kumar, L. De Lorenzis, Unsupervised discovery of interpretable hyperelastic constitutive laws, Comput. Methods Appl. Mech. Eng. 381 (2021) 113852.

[16] P. Thakolkaran, A. Joshi, Y. Zheng, M. Flaschel, L. De Lorenzis, S. Kumar, NN-EUCLID: Deep-learning hyperelasticity without stress data, J. Mech. Phys. Solids 169 (2022) 105076.

[17] D. Claire, F. Hild, S. Roux, A finite element formulation to identify damage fields: the equilibrium gap method, Int. J. Numer. Methods Eng. 61 (2) (2004) 189–208.

[18] S. Meng, A.A.K. Yousefi, S. Avril, Machine-learning-based virtual fields method: application to anisotropic hyperelasticity, Comput. Methods Appl. Mech. Eng. 434 (2025) 117580.

[19] R. Shi, H. Yang, J. Chen, K. Hackl, S. Avril, Y. He, Deep learning without stress data on the discovery of multi-regional hyperelastic properties, Comput. Mech. (2025) 1–30.

[20] L. Li, C.Q. Chen, Equilibrium-based convolution neural networks for constitutive modeling of hyperelastic materials, J. Mech. Phys. Solids 164 (2022) 104931.

[21] A. Benady, E. Baranger, L. Chamoin, Unsupervised learning of history-dependent constitutive material laws with thermodynamically-consistent neural networks in the modified constitutive relation error framework, Comput. Methods Appl. Mech. Eng. 425 (2024) 116967.

[22] A. Benady, E. Baranger, L. Chamoin, NN-Mcre: a modified constitutive relation error framework for unsupervised learning of nonlinear state laws with physics-augmented neural networks, Int. J. Numer. Methods Eng. 125 (8) (2024) e7439.

[23] J. Neggers, O. Allix, F. Hild, S. Roux, Big data in experimental mechanics and model order reduction: today's challenges and tomorrow's opportunities, Arch. Comput. Methods Eng. 25 (2018) 143–164.

[24] B.D. Lucas, T. Kanade, An iterative image registration technique with an application to stereo vision, in: IJCAI'81: 7th International Joint Conference on Artificial Intelligence, 2, 1981, pp. 674–679.

[25] M.A. Sutton, W.J. Wolters, W.H. Peters, W.F. Ranson, S.R. McNeill, Determination of displacements using an improved digital correlation method, Image Vis. Comput. 1 (3) (1983) 133–139.

[26] B.K. Bay, Methods and applications of digital volume correlation, J. Strain Anal. Eng. Des. 43 (8) (2008) 745–760.

[27] A. Buljac, C. Jailin, A. Mendoza, J. Neggers, T. Taillandier-Thomas, A. Bouterf, B. Smaniotto, F. Hild, S. Roux, Digital volume correlation: review of progress and challenges, Exp. Mech. 58 (2018) 661–708.

[28] B. Jordan, M.B. Gorji, D. Mohr, Neural network model describing the temperature-and rate-dependent stress-strain response of polypropylene, Int. J. Plast. 135 (2020) 102811.

[29] M. Diamantopoulou, N. Karathanasopoulos, D. Mohr, Stress-strain response of polymers made through two-photon lithography: micro-scale experiments and neural network modeling, Addit. Manuf. 47 (2021) 102266.

[30] S.R.S. Pierre, D. Rajasekharan, E.C. Darwin, K. Linka, M.E. Levenston, E. Kuhl, Discovering the mechanics of artificial and real meat, Comput. Methods Appl. Mech. Eng. 415 (2023) 116236.

[31] M. Flaschel, H. Yu, N. Reiter, J. Hinrichsen, S. Budday, P. Steinmann, S. Kumar, L. De Lorenzis, Automated discovery of interpretable hyperelastic material models for human brain tissue with EUCLID, J. Mech. Phys. Solids 180 (2023) 105404.

[32] C. Jailin, A. Benady, R. Legroux, E. Baranger, Experimental learning of a hyperelastic behavior with a physics-augmented neural network, Exp. Mech. 64 (9) (2024) 1465–1481. https://doi.org/10.1007/s11340-024-01106-5.

[33] M.A. Moreno-Mateos, S. Wiesheier, A. Esmaeili, M. Hossain, P. Steinmann, Biaxial characterization of soft elastomers: experiments and data-adaptive configurational forces for fracture, (2025). arXiv preprint arXiv:2505.20244.

[34] A.H. Muhr, Modeling the stress-strain behavior of rubber, Rubber Chem. Technol. 78 (3) (2005) 391–425.

[35] S.K. Melly, L. Liu, Y. Liu, J. Leng, A review on material models for isotropic hyperelasticity, Int. J. Mech. Syst. Dyn. 1 (1) (2021) 71–88.

[36] D.P. Kingma, J. Ba, Adam: A method for stochastic optimization, (2014). arXiv preprint arXiv:1412.6980.

[37] B.K. Bay, T.S. Smith, D.P. Fyhrie, M. Saad, Digital volume correlation: three-dimensional strain mapping using X-ray tomography, Exp. Mech. 39 (1999) 217–226.

[38] S. Roux, F. Hild, P. Viot, D. Bernard, Three-dimensional image correlation from X-ray computed tomography of solid foam, Composites Part A: Appl. Sci. Manuf. 39 (8) (2008) 1253–1265.

[39] H. Leclerc, J. Neggers, F. Mathieu, F. Hild, S. Roux, Correli 3.0 [IDDN. FR. 001.520008. 000. SP 2015.000. 31500], Agence pour la Protection des Programmes, Paris (France) (2015).

[40] H. Leclerc, J.-N. Périé, S. Roux, F. Hild, Voxel-scale digital volume correlation, Exp. Mech. 51 (2011) 479–490.

[41] A. Mendoza, J. Neggers, F. Hild, S. Roux, Complete mechanical regularization applied to digital image and volume correlation, Comput. Methods Appl. Mech. Eng. 355 (2019) 27–43.

[42] L. Costecalde, A. Leygue, M. Coret, E. Verron, Data-Driven identification of hyperelastic models by measuring the strain energy density field, Rubber Chem. Technol. 96 (4) (2023) 443–454.

[43] J.C. Criscione, J.D. Humphrey, A.S. Douglas, W.C. Hunter, An invariant basis for natural strain which yields orthogonal stress response terms in isotropic hyperelasticity, J. Mech. Phys. Solids 48 (12) (2000) 2445–2465.

[44] J. Neggers, F. Mathieu, F. Hild, S. Roux, N. Swiergiel, Improving full-field identification using progressive model enrichments, Int. J. Solids Struct. 118 (2017) 213–223.

[45] A. Boubakri, N. Haddar, K. Elleuch, Y. Bienvenu, Impact of aging conditions on mechanical properties of thermoplastic polyurethane, Mater. Des. 31 (9) (2010) 4194–4201.

[46] D.M. Haid, O. Duncan, J. Hart, L. Foster, Characterisation of thermoplastic polyurethane (TPU) for additive manufacturing, in: Tagungsband Zum 14. Symposium Der Sektion Sportinformatik Und Sporttechnologie Der Deutschen Vereinigung Für Sportwissenschaft, p. 85.

[47] H. Lee, R.-i. Eom, Y. Lee, Evaluation of the mechanical properties of porous thermoplastic polyurethane obtained by 3D printing for protective gear, Adv. Mater. Sci. Eng. 2019 (1) (2019) 5838361.

[48] S. Roux, F. Hild, Comprehensive full-field measurements via digital image correlation, Comprehens. Mech. Mater. (2024).

[49] T. Kirchoferer, M. Ortiz, Data-driven computational mechanics, Comput. Methods Appl. Mech. Eng. 304 (2016) 81–101.

[50] C. Jailin, E. Baranger, Material-embedding physics-Augmented neural networks: a first application to constitutive law parameterization, Comput. Methods Appl. Mech. Eng. 445 (2025) 118188.

[51] C. Jailin, A. Bouterf, M. Poncelet, S. Roux, In situ μ CT-scan mechanical tests: fast 4D mechanical identification, Exp. Mech. 57 (2017) 1327–1340.

[52] A.A. Jadoon, K.A. Kalina, M.K. Rausch, R. Jones, J.N. Fuhg, Inverse design of anisotropic microstructures using physics-augmented neural networks, J. Mech. Phys. Solids (2025) 106161.

[53] S. Roux, F. Hild, Optimal procedure for the identification of constitutive parameters from experimentally measured displacement fields, Int. J. Solids Struct. 184 (2020) 14–23.