A Convex Route to Thermoelasticity: Learning Internal Energy and Dissipation

Hagen Holthusen; Paul Steinmann; Ellen Kuhl

cs.CEMar 30, 2026

A Convex Route to Thermoelasticity: Learning Internal Energy and Dissipation

Authors: Hagen Holthusen, Paul Steinmann, Ellen Kuhl

Organizations: Institute of Applied Mechanics, University of Erlangen-Nuremberg, Egerlandstraße 5, 91058 Erlangen, Germany · Department of Mechanical Engineering, Stanford University, United States

Abstract

We present a physics-based neural network framework for the discovery of constitutive models in fully coupled thermomechanics. In contrast to classical formulations based on the Helmholtz energy, we adopt the internal energy and a dissipation potential as primary constitutive functions, expressed in terms of deformation and entropy. This choice avoids the need to enforce mixed convexity--concavity conditions and facilitates a consistent incorporation of thermodynamic principles. In this contribution, we focus on materials without preferred directions or internal variables. While the formulation is posed in terms of entropy, the temperature is treated as the independent observable, and the entropy is inferred internally through the constitutive relation, enabling thermodynamically consistent modeling without requiring entropy data. Thermodynamic admissibility of the networks is guaranteed by construction. The internal energy and dissipation potential are represented by input convex neural networks, ensuring convexity and compliance with the second law. Objectivity, material symmetry, and normalization are embedded directly into the architecture through invariant-based representations and zero-anchored formulations. We demonstrate the performance of the proposed framework on synthetic and experimental datasets, including purely thermal problems and fully coupled thermomechanical responses of soft tissues and filled rubbers. The results show that the learned models accurately capture the underlying constitutive behavior. All code, data, and trained models are made publicly available via https://doi.org/10.5281/zenodo.19248596.

Figures & tables

Figure 1: Overall network architecture consisting of the subnetworks $\mathrm{MLP}_{s}$ , $\mathrm{FICNN}_{s}$ , $\mathrm{FICNN}_{\mbox{\boldmath{\scriptsize{$F$}}}{}}$ , $\mathrm{FICNN}_{\mbox{\boldmath{\scriptsize{$F$}}}{},s}$ , and $\mathrm{PICNN}_{\mbox{\boldmath{\scriptsize{$g$}}}{}}$ , together with their interactions. In a first step, the entropy $s$ is predicted by $\mathrm{MLP}_{s}$ . Subsequently, $s$ and the deformation gradient $F$ are provided as inputs to the respective subnetworks, from which the first Piola–Kirchhoff stress $P$ is obtained via differentiation. In parallel, the referential heat gradient $g$ is processed by $\mathrm{PICNN}_{\mbox{\boldmath{\scriptsize{$g$}}}{}}$ , which is additionally parameterized by $(\mbox{\boldmath$F$}{},s)$ , yielding the referential heat flux $q$ through differentiation. — Figure 1: Overall network architecture consisting of the subnetworks $\mathrm{MLP}_{s}$ , $\mathrm{FICNN}_{s}$ , $\mathrm{FICNN}_{\mbox{\boldmath{\scriptsize{$ F $}}}{}}$ , $\mathrm{FICNN}_{\mbox{\boldmath{\scriptsize{$ F $}}}{},s}$ , and $\mathrm{PICNN}_{\mbox{\boldmath{\scriptsize{$ g $}}}{}}$ , together with their interactions. In a first step, the entropy $s$ is predicted by $\mathrm{MLP}_{s}$ . Subsequently, $s$ and the deformation gradient $F$ are provided as inputs to the respective subnetworks, from which the first Piola–Kirchhoff stress $P$ is obtained via differentiation. In parallel, the referential heat gradient $g$ is processed by $\mathrm{PICNN}_{\mbox{\boldmath{\scriptsize{$ g $}}}{}}$ , which is additionally parameterized by $(\mbox{\boldmath$ F $}{},s)$ , yielding the referential heat flux $q$ through differentiation.

Panel of Figure 2 — Figure 2: Boundary value problem for the transient heat diffusion example. A cubic specimen with edge length $1\,\mathrm{mm}$ is subjected to a time-dependent temperature at the exterior boundary and a constant temperature at the interior boundary. The prescribed loading is given by $T(t)=\left[300-293.15\right]\sin\!\left(\frac{\pi}{4}t\right)+293.15,$ while the initial temperature is $T_{\mathrm{init}}=293.15\,\mathrm{K}$ . The domain is discretized using $4\times 4\times 4$ elements with edge length $0.25\,\mathrm{mm}$ .

Panel of Figure 3 — Figure 3: Temperature distribution at two exemplary time steps along the wall thickness direction for the transient heat diffusion problem. Although the temperatures at the interior and exterior boundaries are identical, a nonlinear temperature profile emerges within the wall as a result of transient effects. The edge length of one finite element is $0.25$ mm.

Panel of Figure 4 — Figure 4: Left: Training history of the physics-based neural network for the transient heat diffusion problem. Shown are the total loss as well as its individual contributions, including the physics-based residual losses, the loss of the auxiliary network, and the regularization term, over the course of training iterations. Right: Relative activity according to Eq. 48 of the individual subnetworks for the transient heat diffusion problem. The diagram shows the normalized contributions of internal energy, dissipation, and auxiliary MLP to the overall activity. It can be observed that the thermal response is dominated by dissipative effects, while the remaining contributions are of significantly smaller magnitude.

Panel of Figure 5 — Figure 5: Comparison of the heat flux component $q_{1}$ for the transient heat-conduction problem. Left: Temporal evolution of the heat flux at selected inner and outer locations, comparing the reference solution with the predictions of the auxiliary MLP and the Newton solver. Right: Parity plot of predicted versus reference heat flux values.

$T=310.15$ K		$T=318.15$ K		$T=323.15$ K		$T=333.15$ K		$T=343.15$ K		$T=353.15$ K
$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]
1.0	0.0000	1.0	0.0000	1.0	0.0000	1.0	0.0000	1.0	0.0000	1.0	0.0000
1.2	0.0217	1.2	0.0146	1.2	0.0146	1.2	0.0108	1.2	0.0108	1.2	0.0108
1.4	0.0436	1.4	0.0314	1.4	0.0314	1.4	0.0221	1.4	0.0221	1.4	0.0221
1.5	0.0787	1.5	0.0527	1.5	0.0527	1.5	0.0320	1.5	0.0320	1.5	0.0320
1.6	0.1756	1.6	0.1069	1.6	0.0875	1.6	0.0469	1.6	0.0469	1.6	0.0469

Table 1: Piola stress

P_{11}

as a function of the axial stretch

F_{11}

for incompressible uniaxial loading at different temperatures for the porcine tissue [ 64 ] . The table summarizes the thermo-mechanical response used to generate the training data, illustrating the temperature dependence of the constitutive behavior.

Panel of Figure 6 — Figure 6: Left: Training history of the physics-based neural network for the temperature dependent porcine tissue dataset [ 64 ] . Shown are the total loss as well as its individual contributions, including the physics-based residual losses and the regularization term, over the course of training iterations. Right: Relative activity according to Eq. 48 of the individual subnetworks for the porcine tissue. The diagram shows the normalized contributions of internal energy. No transient thermal effects are present. Training is performed using an iterative solver instead of the auxiliary network.

Panel of Figure 7 — Figure 7: Comparison of the discovered model and the reference data for porcine tissue under incompressible uniaxial loading at different temperatures. The Piola stress $P_{11}$ is plotted as a function of the axial stretch $F_{11}$ . The reference data are shown as discrete points, while the continuous curves represent the response of the discovered model. Left: Enlarged view near $F_{11}=1$ . Right: Full stretch range, demonstrating that the discovered model accurately captures the strongly nonlinear and temperature-dependent constitutive response.

$T=283$ K		$T=303$ K		$T=323$ K		$T=343$ K		$T=363$ K		$T=373$ K		$T=383$ K		$T=393$ K
$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]	$F_{11}$ [-]	$P_{11}$ [MPa]
1.00000	0.000000	1.00000	0.000000	1.00000	0.000000	1.00000	0.000000	1.00000	0.000000	1.00000	0.000000	1.00000	0.000000	1.00000	0.000000
1.06605	0.450927	1.07206	0.398578	1.06455	0.351218	1.05925	0.302547	1.06152	0.285128	1.07430	0.397746	1.04399	0.238756	1.06077	0.318900
1.11608	0.638094	1.12397	0.601937	1.13816	0.617895	1.11824	0.581138	1.11104	0.480744	1.13038	0.598524	1.10160	0.509104	1.10742	0.530547
1.17084	0.836345	1.19056	0.852403	1.21525	0.849735	1.17495	0.765233	1.15913	0.622078	1.18448	0.781625	1.15635	0.746744	1.15566	0.724084
1.22446	1.003285	1.23686	0.971640	1.27205	1.014738	1.22774	0.899099	1.19393	0.738158	1.22118	0.874767	1.19527	0.863940	1.19929	0.861042

Table 2: Piola stress

P_{11}

as a function of the axial stretch

F_{11}=\lambda

for incompressible uniaxial loading at different temperatures for the carbon-filled black rubber [ 65 ] . The table summarizes the thermo-mechanical response used to generate the training data, illustrating the temperature dependence of the constitutive behavior.

Panel of Figure 8 — Figure 8: Left: Training history of the physics-based neural network for the temperature dependent carbon-filled black rubber dataset [ 65 ] . Shown are the total loss as well as its individual contributions, including the physics-based residual losses and the regularization term, over the course of training iterations. Right: Relative activity according to Eq. 48 of the individual subnetworks for the carbon-filled black rubber. The diagram shows the normalized contributions of internal energy. No transient thermal effects are present. Training is performed using an iterative solver instead of the auxiliary network.

Panel of Figure 9 — Figure 9: Comparison of the discovered model and the reference data for carbon-filled black rubber under incompressible uniaxial loading at different temperatures. The Piola stress $P_{11}$ is plotted as a function of the axial stretch $F_{11}$ . The reference data are shown as discrete points, while the continuous curves represent the response of the discovered model. Left: Full stretch range. Right: Enlarged view, demonstrating that the discovered model is unable to explain the thermal stiffening effect in high temperature ranges.

Figure 10: Setup of the boundary value problem and training scenarios for the plate with an elliptic hole. Top: Geometry and dimensions of the specimen [ 70 ] , the thickness is $t=1$ mm. Middle: Three distinct training scenarios with different combinations of mechanical loading $u_{y}(t)$ and thermal loading $T(t)$ applied at the top boundary, while the bottom boundary is fixed and the initial temperature is $T_{\mathrm{init}}=293.15\,\mathrm{K}$ . Bottom: Corresponding temporal evolution of the prescribed displacement and temperature for each training case. Solid black lines denote the prescribed displacement, whereas the temperature evolution is represented by colored lines transitioning from red to blue.

Panel of Figure 11 — Figure 11: Left: Training history of the physics-based neural network for the plate with an elliptic hole problem. Shown are the total loss as well as its individual contributions, including the physics-based residual losses, the loss of the auxiliary network, and the regularization term, over the course of training iterations. Right: Relative activity according to Eq. 48 of the individual subnetworks for the plate with an elliptic hole problem. The diagram shows the normalized contributions of internal energy, dissipation, and auxiliary MLP to the overall activity.

Figure 12: Boundary value problem and testing scenarios for the spring-like structure. Left: Geometry and boundary conditions, including prescribed displacement $u_{y}(t)$ and temperature $T(t)$ at the top boundary, while the bottom boundary is fixed and maintained at constant temperature $T=293.15\,\mathrm{K}$ . The initial temperature is set to $T_{\mathrm{init}}=293.15\,\mathrm{K}$ . Right: Two testing scenarios used for model evaluation. In Testing #1, coupled thermo-mechanical loading is applied via time-dependent displacement and temperature. In Testing #2, only mechanical loading is applied while the temperature is kept constant, such that thermal effects arise solely from the deformation. Solid black lines denote the prescribed displacement, whereas the temperature evolution is represented by colored lines transitioning from red to blue.

Figure 13: Results for Testing #1 of the spring-like structure under coupled thermo-mechanical loading. Shown are five representative time steps during the loading process. For each snapshot, the temperature field $T$ , the heat flow $Q$ , and the forces $F_{y}$ are depicted. The reference solution $(\bullet)_{\mathrm{ref}}$ is compared against the predictions of the discovered model. The results demonstrate that the learned model accurately captures the fully coupled thermo-mechanical response.

Panel of Figure 14 — Figure 14: Quantitative comparison of the predicted and reference responses for Testing #1 of the spring-like structure under coupled thermo-mechanical loading. Left: Temporal evolution of the reaction heat flow $Q$ (top) and the reaction forces $F_{y}$ (bottom) at the top and bottom boundaries. Right: Corresponding parity plots comparing predictions and reference values. The results show good agreement for both mechanical and thermal quantities, with only minor deviations in the heat flow at the bottom boundary.

Figure 15: Results for Testing #2 of the spring-like structure under purely mechanical loading with constant temperature boundary conditions. Shown are five representative time steps during the loading process. For each snapshot, the temperature field $T$ , the heat flow $Q$ , and the forces $F_{y}$ are depicted. The reference solution $(\bullet)_{\mathrm{ref}}$ is compared against the predictions of the discovered model. While the forces are accurately reproduced, noticeable discrepancies are observed in the heat flow, including non-physical fluxes within the interior of the domain.

Panel of Figure 16 — Figure 16: Quantitative comparison of the predicted and reference responses for Testing #2 of the spring-like structure under purely mechanical loading with constant temperature boundary conditions. Left: Temporal evolution of the reaction heat flow $Q$ (top) and the reaction forces $F_{y}$ (bottom) at the top and bottom boundaries. Right: Corresponding parity plots comparing predictions and reference values. While the reaction forces are accurately captured, significant deviations are observed in the heat flow, both in magnitude and temporal evolution, including spurious oscillations and increased scatter in the parity plot.

Appendix figures & tables2 assets

Supplementary material from the paper’s appendix.

Appendix

Category	Parameter	Value
Optimization	Number of epochs	-
	Learning rate	$1\times 10^{-3}$
	Gradient clipping	disabled
	Clip norm	$1.0$
	Optimizer	ADAM
Loss weighting	$\lambda_{\mathrm{D}}$ (Dirichlet force)	$1$

Appendix

Table 3: Training hyperparameters.

Network	Layers	Activations	Output	Output act.	Output bias
FICNN ${}_{\mbox{\boldmath{\scriptsize{$ F $}}}{},s}$ (joint energy)	[12, 12]	exp, exp	1	softplus	Yes
FICNN ${}_{\mbox{\boldmath{\scriptsize{$ F $}}}{}}$ (energy features)	[12, 12, 12]	exp, softplus, softplus	12	identity	Yes
FICNN s (entropy features)	[12, 12, 12]	exp, exp, exp	12	identity	Yes
PICNN ${}_{\mbox{\boldmath{\scriptsize{$ g $}}}{}}$ (dissipation potential)	[12, 12, 12]	softplus, softplus, softplus	1	ReLU	No
PICNN (coupling branch)	[6, 6, 6]	GELU, GELU, GELU	–	–	–
MLP s (auxiliary entropy)	[12, 12]	GELU, GELU	1	identity	Yes

Appendix

Table 4: Neural network architectures.

Materials Science

Materials science studies how composition, structure, processing, and defects determine material properties and how to design useful materials.

See all 301 papers

May 29, 2026

Discovering Thermodynamically Admissible Dissipation Potentials via Grammar-Based Symbolic Regression

Constitutive laws for inelastic materials must satisfy strict thermodynamic admissibility requirements, yet current data-driven approaches sacrifice interpretability, even when formal guarantees are provided by physics-encoded architectures. We propose a symbolic regression framework for the data-driven discovery of dissipation potentials governing the evolution of internal variables within the Generalized Standard Materials (GSM) formalism. Starting from the Clausius--Duhem inequality, we enforce the thermodynamic requirements, convexity and non-negativity, that the dual dissipation potential must satisfy to guarantee non-negative mechanical dissipation. These requirements are formulated in the general subdifferential setting, encompassing rate-dependent (viscoelastic) and viscoplastic dissipative mechanisms, including potentials with genuine elastic domains, within a unified framework. Candidate potentials are generated by a composition-extended convexity-preserving grammar that guarantees thermodynamic admissibility by construction. The framework is validated on synthetic datasets spanning Newtonian, power-law, and Bingham viscoplastic ground truths under process and measurement noise, and on experimental oscillatory shear measurements of a synthetic elastomer across multiple strain amplitudes and frequencies, where the discovered potentials reproduce the amplitude-dependent softening of the dynamic moduli and outperform a calibrated linear Zener baseline.

Federico Califano, Jacopo Ciambella

Materials Science Computational Materials Science Symbolic Regression

Ask PDF

Feb 19, 2026

Inelastic Constitutive Kolmogorov-Arnold Networks: A generalized framework for automated discovery of interpretable inelastic material models

A key problem of solid mechanics is the identification of the constitutive law of a material, that is, the relation between strain history and stress. Machine learning has lead to considerable advances in this field lately. Here we introduce inelastic Constitutive Kolmogorov-Arnold Networks (iCKANs). This novel artificial neural network architecture can discover in an automated manner symbolic constitutive laws describing both the elastic and inelastic behavior of materials. That is, it can translate data from material testing into corresponding elastic and inelastic potential functions in closed mathematical form. We demonstrate the advantages of iCKANs using both synthetic data and experimental data of the viscoelastic polymer materials VHB 4910 and VHB 4905. The results demonstrate that iCKANs accurately capture complex viscoelastic behavior while preserving physical interpretability. It is a particular strength of iCKANs that they can process not only mechanical data but also arbitrary additional information available about a material (e.g., about temperature-dependent behavior). This makes iCKANs a powerful tool to discover in the future also how specific processing or service conditions affect the properties of materials.

Chenyi Ji, Kian P. Abdolazizi, Hagen Holthusen +2

Computational Mechanics in Medicine, Applied Medical Engineering, RWTH Aachen University Pauwelsstraße 20, 52074 Aachen, Germany · Institute for Continuum and Material Mechanics, Hamburg University of Technology Eißendorfer Straße 42, 21073 Hamburg, Germany · Institute of Applied Mechanics, University of Erlangen-Nuremberg Egerlandstrasse 5, 91058 Erlangen, Germany +1

Materials Science Neural Network Interpretability Kolmogorov-Arnold Networks Computational Materials Science Symbolic Regression

Ask PDF

Mar 27, 2026

Hyperelastic constitutive model discovery with differentiable finite elements and structure-preserving neural networks

The discovery of constitutive laws from experimentally accessible measurements is a central problem in nonlinear computational mechanics. Many data-driven constitutive identification approaches rely either on paired strain-stress data or on full-field displacement measurements, both of which are difficult to obtain in realistic three-dimensional settings. We present a differentiable finite element framework for the discovery of hyperelastic material laws from partial observations, including boundary-only displacement measurements and global reaction forces. The method embeds the nonlinear finite element equilibrium problem directly into the learning loop, so that candidate strain-energy densities are assessed through the deformation fields and reactions they induce. This formulation enforces mechanical equilibrium as a constraint and allows the loss function to be evaluated only at observed locations. To ensure physical admissibility and promote numerical solvability throughout training, the constitutive response is represented by Hyperelastic Neural Networks, a structure-preserving neural class that enforces residual energy and stress-free conditions, frame indifference, isotropic material symmetry, polyconvexity, coercivity, and controlled volumetric growth by construction. The resulting PDE-constrained learning problem is solved using a quasi-Newton strategy combined with continuation and solver-aware backtracking. Numerical experiments in two- and three-dimensional finite elasticity demonstrate accurate recovery of hyperelastic isotropic responses from boundary-only data, robustness to measurement noise, and generalization across geometries, loading conditions, and boundary conditions.

Francesco Regazzoni

MOX, Department of Mathematics, Politecnico di Milano, Milan, Italy

Structure-Preserving Neural Networks Hyperelasticity Physics-Informed ML Differentiable Physics

Ask PDF

Abstract

Materials Science

Explore similar work

Discovering Thermodynamically Admissible Dissipation Potentials via Grammar-Based Symbolic Regression

Inelastic Constitutive Kolmogorov-Arnold Networks: A generalized framework for automated discovery of interpretable inelastic material models

Hyperelastic constitutive model discovery with differentiable finite elements and structure-preserving neural networks