Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials

A novel data-driven constitutive modeling approach is proposed, which combines the physics-informed nature of modeling based on continuum thermodynamics with the benefits of machine learning. This approach is demonstrated on strain-rate-sensitive soft materials. This model is based on the viscous dissipation-based visco-hyperelasticity framework where the total stress is decomposed into volumetric, isochoric hyperelastic, and isochoric viscous overstress contributions. It is shown that each of these stress components can be written as linear combinations of the components of an irreducible integrity basis. Three Gaussian process regression-based surrogate models are trained (one per stress component) between principal invariants of strain and strain rate tensors and the corresponding coefficients of the integrity basis components. It is demonstrated that this type of model construction enforces key physics-based constraints on the predicted responses: the second law of thermodynamics, the principles of local action and determinism, objectivity, the balance of angular momentum, an assumed reference state, isotropy, and limited memory. The three surrogate models that constitute our constitutive model are evaluated by training them on small-size numerically generated data sets corresponding to a single deformation mode and then analyzing their predictions over a much wider testing regime comprising multiple deformation modes. Our physics-informed data-driven constitutive model predictions are compared with the corresponding predictions of classical continuum thermodynamics-based and purely data-driven models. It is shown that our surrogate models can reasonably capture the stress-strain-strain rate responses in both training and testing regimes, and provide improvements in terms of prediction accuracy, generalizability to multiple deformation modes, and compatibility with limited data.


Introduction
Constitutive models are equations (or sets of equations) that describe the response of a material to imposed loads, deformations, or temperature changes.The philosophy behind formulating constitutive models has evolved considerably over the years, developing through four stages which are illustrated in Fig. 1.The oldest stage can be traced back to Robert Hooke in the Latin anagram "ut tensio, sic vis" (as the extension, so the force) (Hooke (1678)), and involved writing simplified and often empirical equations to describe observations of a material's response within a restricted range of loading conditions (e.g., the ideal Hookean elastic solid and the ideal Newtonian viscous fluid).The emergence of the modern continuum theory of constitutive modeling in the 1950s introduced a different philosophy: it begins with a very general functional constitutive equation that conforms to a set of constraints imposed by certain physical laws and thermodynamic considerations, and simplifying assumptions to obtain the finalized model (based on the observed material response) are imposed as late and as little as possible (Coleman and Noll (1963); Malvern (1969); Truesdell and Noll (2004)).This approach especially benefited the formulation of constitutive models of materials that exhibit complex deformation features (e.g., large deformations) by mitigating the requirement of an extensive experimental exploration to propose constitutive equations.One such class of materials is soft materials, whose mechanical response involves large deformation, nonlinear stress-strain behavior, and strain rate (or time)-dependence.The present study focuses on the constitutive modeling of soft materials (e.g., elastomers, hydrogels, and soft tissues) with a rate-sensitive response.
The third stage in the evolution of constitutive modeling emerged with the popularity of machine learning (ML) tools in computer sciences.These "black-box" models create simple mappings between stress and deformation, and have been extensively applied in the past 30 years on many material classes including soft materials (Wu et al. (1990); Ghaboussi et al. (1998); Lefik and Schrefler (2003); Jung and Ghaboussi (2006); Huang et al. (2020); Fuhg et al. (2021a); Logarzo et al. (2021)).ML offers numerous advantages in regard to describing mechanical responses: (i) it can capture complex trends when given sufficient training data, (ii) it can directly utilize experimental data with no requirement of physical laws or expert knowledge of response trends, and (iii) it can offer potentially lower computational cost and high-throughput.Despite these advantages, the accuracy of these black-box models is limited due to their requirement of a large amount of training data, which is often not available in real experimental settings.In addition, recent studies have shown that because these models are not restricted by physical laws, they show very poor prediction accuracy in the deformation regimes (say, > 10% strain) that are not considered during model training (say, 10% strain) (Fuhg and Bouklas (2022)).
The latest stage in the constitutive modeling evolution attempts to combine the physics-informed nature of the continuum theory of constitutive modeling (i.e., the second stage) with the flexibility and efficiency of classical ML (i.e., the third stage) to formulate "physics-informed data-driven constitutive models".For example, Liu et al. (2020) recently proposed a physics-informed neural network material model (NNMat) for isotropic hyperelastic soft tissues.This model is based on a mapping between the invariants of the right Cauchy-Green deformation tensor and the derivatives of the strain energy density with respect to the invariants, and imposes convexity constraints in the NN loss function to ensure physically reasonable predictions.Unlike this study, Frankel et al. (2020) employed Gaussian process regression (GPR) to create a mapping between the invariants of the right Cauchy-Green deformation tensor and the coefficients of the irreducible integrity basis of the stress tensor (for isotropic hyperelastic solids).Fuhg and Bouklas (2022) showed that this type of generalized functional constitutive equation (i.e., stress written as a linear function of the irreducible integrity basis) automatically ensures a number of physical constraints on material behavior: material frame-indifference, material symmetry, the balance of angular momentum, and the second law of thermodynamics (as the Clausius-Planck inequality).Further, owing to the ability of GPR to exactly predict a training data point (called the exact inference property), the physical constraint of a stress-free reference state was achieved by simply including the zero stress-zero strain data point in the training data.
Although important and useful, these GPR-based models are limited in applicability to hyperelastic solids that assume no strain-rate-dependence in material behavior, and were trained using a large amount of artificially generated data.In practice, soft materials often show strain-rate dependence in material response.In addition, sparse data at only a few strain rate levels and in a limited range of strain values is usually available from mechanical experiments (Upadhyay et al. (2019a(Upadhyay et al. ( , 2021)); Luo et al. (2019)).Thus, the goal of the present study is to develop a physics-informed data-driven constitutive model that can capture the strain-rate-sensitive mechanical response of visco-hyperelastic soft materials from limited training data, consistent with current experimental practice.To this end, this work leverages the recent contributions by Frankel et al. (2020) and Fuhg and Bouklas (2022) for hyperelastic soft materials.Note that strain rate-dependence in visco-hyperelasticity can be classified into short-time (e.g., high strain rate deformation) and long-time (e.g., creep and relaxation) responses (Upadhyay et al. (2020b)).This work focuses only on the short-time response, which is of special interest in the field of injury biomechanics (e.g., simulations of crashes, blast, and ballistic impact) and in the design of protective equipment (Bracq et al. (2018); Harrigan et al. (2010); Upadhyay et al. (2022a); Payne et al. (2015)).
The paper is organized as follows: Section 2 formulates a generalized functional constitutive equation for viscohyperelastic soft materials, which forms the basis for the data-driven model proposed in this study.In Section 3, a data-driven mapping is defined for the generalized constitutive equation.The physics-based constraints on this data-driven mapping, which stem from both the generalized model construction and the mapping approach, are also summarized.Gaussian process regression is the primary supervised learning method utilized in this work, which allows the imposition of several additional physical constraints.These are described in Section 4. Section 5 demonstrates the fitting and prediction performance of our GPR-based physics-informed data-driven constitutive model on several numerical tests that consider multiple deformation modes and a wide range of loading rates.The performance of our model is also compared against both conventional visco-hyperelastic constitutive models (i.e, from the second stage of the constitutive modeling evolution) and the classical ML mapping models (i.e., from the third stage).Finally, Section 6 presents a summary of this work.

Generalized Visco-hyperelastic Constitutive Framework
In general, every admissible thermomechanical process (and so the choice of a constitutive model) must satisfy the Clausius-Duhem inequality, Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT where ρ 0 is the mass density, η 0 is the specific entropy, Q is the heat flux, R is the rate of internal heating per unit mass, and θ is the absolute temperature.Subscript 0 denotes the reference system configuration, and the operator (∇ 0 •) represents divergence.Further, considering the reference configuration local form of the conservation of energy, where u 0 is the specific internal energy, T 0 is the nominal stress tensor, and F is the deformation gradient tensor.The symbol •• represents the tensor scalar product, such that T 0 • • Ḟ = T 0 ij Ḟji in rectangular Cartesians.By substituting R from Eq. (2) in Eq. ( 1) and introducing the Helmholtz free energy function, Restricting the focus to isothermal deformations only, we particularize the second law of thermodynamics through the Clausius-Planck inequality as where Ξ int is the internal dissipation or local entropy production (Limbert and Middleton (2004)).Equation 4represents the primary physical constraint on the possible mathematical forms of our generalized constitutive equation.
In addition to the Clausius-Planck inequality, another common physical law that constrains the forms of constitutive equations is the principle of local action, which states that material response at a given point depends only on conditions in the close vicinity of that point (e.g., ψ = ψ(θ, F, ∇ 0 F, . . .)).Assuming that the Helmholtz free energy function depends only on θ and F (i.e., neglecting higher-order gradients) so that ψ = ψ(θ, F), Eq. ( 4) becomes In terms of the right Cauchy-Green deformation tensor, C = F T F, and the hyperelastic strain energy density function, W h , For an ideal hyperelastic material, the inequality in Eq. ( 6) reduces to equality under the assumption that the material undergoes zero dissipation during deformation.This leads to a rate-independent form of the constitutive model.While this assumption holds reasonably well in the case of quasi-static deformations, the response of soft materials under high strain rate deformation is associated with irreversible thermodynamic processes resulting from viscous dissipation effects stemming from the microscale level (Vogel et al. (2017); Pioletti (2006)).Note that "high strain rate" here is synonymous to the short-time response, when the deformation time scale is very small compared to the time it takes for the internal material microstructure to rearrange or relax under load.Under these conditions, the existence of a viscous dissipation potential (also called pseudo-potential) W v for accounting energy dissipation is usually assumed (Germain (1998); Pioletti et al. (1998); Pioletti and Rakotomanana (2000); Zhurov et al. (2007); Upadhyay et al. (2020b)), Rearranging Eq.( 7) and introducing the second Piola-Kirchhoff stress S (S = F −1 T 0 T ), From Eq. (8, the total stress in the material is additively decomposed into a hyperelastic stress component (S h ) and a viscous overstress component (S v ).Note that because both the hyperelastic strain energy density and the viscous dissipation potential are functions of symmetric tensors (i.e., C and Ċ), the stress tensor S is also symmetric.This ensures the physical constraint of the balance of angular momentum.Further, it can be readily confirmed that Eq. ( 8) also satisfies the principle of determinism or causality (i.e., material response at a given instant is only a function of past and present events) Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT Next, imposing another physical constraint of isotropic material symmetry, Eq. ( 8) can be written in terms of certain scalar invariants of the tensors C and Ċ (Malvern (1969); Upadhyay et al. (2020b)), where The generalized constitutive equation in Eq. ( 9) will now be converted into an equivalent form, without loss of generality, to be compatible with our data-driven mapping approach (described later).To this end, following the standard approach of modeling compressible materials (Holzapfel (2000)), the tensors F and C are multiplicatively decomposed into dilatational and isochoric components, where J is an invariant that is given as and the modified deformation gradient tensor F along with the modified right Cauchy-Green deformation tensor C are introduced.These are associated with the isochoric (volume-preserving) part of the deformation, such that, The decomposition of the deformation gradient in Eq. ( 11) allows an equivalent decomposition of the stress tensor in Eq. ( 9) into dilatational and isochoric stress components (see Zhurov et al. (2007)), where are the invariants.In Eq. ( 14), S vol is the volumetric stress component, S h,vol is the isochoric hyperelastic stress component, and S v,vol is the isochoric viscous overstress component.These three stress components are captured by the volumetric energy density function U (J), the modified hyperelastic strain energy density Wh , and the modified viscous dissipation potential Wv , respectively.For convenience, Wh and Wv will henceforth be simply referred to as the strain energy density and the viscous dissipation potential, respectively.
Equation ( 14) is the generalized functional form of constitutive equation considered in this work.This generalized framework (referring to both Eq. ( 9) and Eq. ( 14)) is called external state variable-driven viscous dissipation-based visco-hyperelasticity (Pioletti and Rakotomanana (2000); Upadhyay et al. (2020b)).Here, C and Ċ are the external thermodynamic state variables that are employed to relate the strain rate stiffening or softening of material response to viscous dissipation.This is in contrast to the internal state variable-driven framework in which the non-equilibrium part of the Helmholtz free energy function is described via a set of internal variables that a priori lack any physical meaning Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT (e.g., see Simo (1987); Reese and Govindjee (1998); Holzapfel and Gasser (2001); Garcia-Gonzalez et al. (2018)).Unlike this alternative class of visco-hyperelastic models, many of which result in a hereditary-integral-based equation for stress, the external state variable-driven viscous dissipation-based visco-hyperelastic constitutive framework of this work is bound by the constraint of limited memory (a subset of the principle of fading memory).Here, limited memory means that the viscous material behavior is dependent only on the instantaneous deformation rate (i.e., very recent history), and dependence on the entire previous loading history (often described via a hereditary-integral) is neglected.First proposed by Pioletti et al. (1998), constitutive models of this framework have been employed to successfully capture the rate-dependent response of numerous soft materials under rapid loading: human ligaments and tendons (Jiang et al. (2015);Zhurov et al. (2007); Ahsanizadeh and Li (2015)), skeletal muscles (Lu et al. (2010)), hydrogels and elastomers (Upadhyay et al. (2020b(Upadhyay et al. ( , 2022a))), tongue tissue (Yousefi et al. (2018)), and the brain and pericardium (Kulkarni et al. (2016); Upadhyay et al. (2021Upadhyay et al. ( , 2022b))), among others.These studies assume specific mathematical forms for U (J), Wh ( Ī1 , Ī2 ) and Wv ( Ī1 , Ī2 , J1 , J2 , J3 , J4 , J5 , J6 , J7 ) based on expert knowledge and experience.The physics-informed data-driven mapping approach of the present study (Section 3) aims to eliminate this need for expert intervention by employing ML to flexibly discover the mapping between stress-and strain-like variables directly from limited experimental data.
The partial derivatives representing the three stress components in Eq. ( 14) can be expanded via the chain rule, leading to the following generalized equations for the three stress components (see derivation in the Appendix A.): where Dev(•) = (•) − 1 3 ((•) : C)C −1 is the deviatoric operator in the Lagrangian description (Holzapfel (2000)).Notice that individual stress components in these equations are linear combinations of the components of an integrity basis G of the total stress S, defined as where . ., Φ 7 are coefficients of the integrity basis components in the stress equations (Eqs.(16-18)), and are functions of the invariants of the tensors C and Ċ.For several commonly used constitutive models in the literature (i.e., for a given choice of U , Wh , and Wv ), the explicit mathematical forms of these coefficients are provided in Tables A.1-A.3 of the Appendix A.. The physics-informed data-driven mapping approach of the present study discovers a mapping between invariants and coefficients of the integrity basis components directly from stress-strain-strain rate data.
3 Proposed Physics-Informed Data-Driven Mapping Approach Consider a dataset D consisting of strain and strain rate as input and stress components as output, Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT Here, constituent datasets D vol , D h,iso and D v,iso correspond to the material response under hydrostatic, isochoric quasi-static, and isochoric dynamic conditions, respectively.Superscripts i,j, and k denote a particular data point in the datasets D vol , D h,iso , and D v,iso , respectively.A dataset like D for a soft material can be obtained in practice by compiling its hydrostatic stress-strain data as D vol (Nie et al. (2013); Saraf et al. (2007)), quasi-static stress-strain data under uniaxial and/or shear deformations as D h,iso (Yang et al. (2000); Upadhyay et al. (2020a)), and viscous overstress (total stress minus stress under quasi-static loading)-strain data from high strain rate testing at multiple strain rate levels under uniaxial and/or shear deformations as D v,iso (Yang et al. (2000); Saraf et al. (2007); Bracq et al. (2017); Upadhyay et al. (2022aUpadhyay et al. ( , 2019a))).For every data point in the dataset D, the inputs C and Ċ can be used to compute the set of invariants (using Eqs. ( 12) and ( 15)) as well as the components of the integrity basis G (using Eq. ( 20)).The coefficients of the integrity basis in Eqs.(16-18) can then be obtained by solving the following systems of equations, respectively: where vec(•) denotes the Voigt form of a matrix (for an arbitrary second-order symmetric tensor Z, vec(Z) = Z 11 , Z 22 , Z 33 , Z 23 , Z 13 , Z 12 T ).Notice that the above systems of equations are identical to Eqs. (16)(17)(18), and that each of these systems is written in the form of b = A x .The vectors x containing the coefficients of the integrity basis can be obtained by solving the corresponding optimization problems-min x A x − b 2 2 -via, for example, QR-decomposition.
The above procedure of obtaining invariants and the coefficients of the integrity basis is carried out for every data point in D. Using this information, an alternative dataset is generated that comprises invariants and integrity basis coefficients corresponding to each of the stress component types, The constituent datasets in D * are used in the present study to train the three surrogate models that comprise our physics-informed data-driven constitutive model: Specifically, M vol captures the volumetric stress component and is a mapping between the random vector of the invariant J (i.e., J vol ) and the random vector of the coefficient ζ 1 (i.e., ζ).M h,iso captures the isochoric hyperelastic stress component and is a mapping between the random vector of the invariants [ Ī1 , Ī2 ] (i.e., I) and the random vector of the corresponding coefficients [Γ 1 , Γ 2 ] (i.e., Γ ).Finally, M v,iso captures the isochoric viscous overstress and is a mapping between the random vector of the invariants [ Ī1 , Ī2 , J1 , J4 , J6 ] (i.e., J ) and the random vector of the Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT corresponding coefficients [Φ 1 , . . ., Φ 7 ] (i.e., Φ).In a more explicit form, we have Note, the rationale behind not considering invariants J2 , J3 , J5 and J7 in the surrogate model M v,iso (Eq.(28/31)) is described in the next section.Briefly, this assumption of a negligible effect of certain invariants on the viscous overstress allows imposition of the physical constraint of a stress-free reference state; i.e., to ensure S v,iso (C = I) = 0 regardless of the applied rate of the right Cauchy-Green deformation tensor ( Ċ).
In this work, the Gaussian Process Regression (GPR) supervised learning method (Williams and Rasmussen (1995)) is chosen for surrogate modeling purposes (see details of GPR in Section 4), i.e., to learn the mappings in Eqs.(26-28) from a given input-output training dataset D and then predict integrity basis coefficients for a new set of input tensors C and Ċ.Given the linear relationship between individual stress components and their corresponding integrity basis coefficients (Eq.(16-18)), the trained surrogate models allow prediction of the three stress components as In Eqs.(32-34), the accent symbols • over stress and coefficients of the integrity basis denote that these are approximated values predicted via surrogate models.Note, for any given input strain and strain rate, the integrity basis components (G i/j/k ) can be obtained using Eq. ( 20)).Fuhg et al. (2022a) have recently proposed a data-driven method that allows discovering which of the integrity basis components play an active role in the stress response, i.e. discovery of both the type and orientation of the anisotropy.The entire process of the development of our physics-informed data-driven constitutive model from a stress-strain-strain rate dataset (i.e., D) is summarized in Algorithm 1.
As the proposed surrogate models in Eqs.(26-28) are based on the generalized functional forms of the visco-hyperelastic stress equations (i.e., Eq. (16-18)), our data-driven constitutive model and its predictions (Eqs.(32-34)) automatically satisfy a number of physics-based constraints: • Principle of local action: Because the integrity basis only captures spatially local deformation (via C and Ċ).
• Balance of angular momentum: Because the integrity basis is necessarily symmetric.
• Principle of determinism (or causality): Because the predicted stress is only a function of past and present events.• Principle of material frame-indifference (objectivity): Because S, C, and Ċ are all objective tensors associated with the reference configuration and thus remain unaffected under change of observer.• Isotropic material symmetry: Because the surrogate model stress equations ) remain unaffected by any transformation from the proper orthogonal material symmetry group SO(3) (i.e., under rigid body rotation).This is a direct consequence of the generalized model formulation in Eq. ( 14) in the form of isotropic invariants of tensors C and Ċ.
Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT Algorithm 1 Development of a physics-informed data-driven constitutive model for strain-rate-sensitive soft materials Result: Corresponding surrogate models M vol , M h,iso , and M v,iso Obtain the invariant J i from C i using Eq. ( 12) 3: Obtain the integity basis component G i 1 using Eq. ( 20). 4: Obtain the Voigt form of Obtain the Voigt form of S i vol : vec(S i vol ) ∈ R 6 6: Solve the system of equations in Eq. ( 22) to obtain Obtain the invariants J j , Īj 1 , and Īj 2 using Eqs.( 12) and ( 15) 12: Obtain the integrity basis components G j 2 and G j 3 using Eq. ( 20) 13: Obtain the Voigt forms and construct the matrix Obtain the Voigt form of S j h,iso : vec(S j h,iso ) ∈ R 6 15: Solve the system of equations in Eq. ( 23) to obtain Obtain the invariants J k , Īk 1 , Īk 2 , Jk 1 , Jk 4 , and Jk 5 using Eqs.( 12) and ( 15) 21: Obtain the integrity basis components , and G k 8 using Eq. ( 20) 22: Obtain the Voigt forms and construct the matrix Obtain the Voigt form of S k v,iso : vec(S k v,iso ) ∈ R 6 24: Solve the system of equations in Eq. ( 23) to obtain • Limited memory constraint (a subset of the principle of fading memory): Because the viscous overstress component (Eq.( 34)) depends only on the instantaneous deformation rate (i.e., very recent history).
Two additional physics-based constraints will be imposed via the employment of GPR as the chosen surrogate modeling technique: • the stress-free reference state, and • the second law of thermodynamics.These additional physical constraints along with the GPR technique will now be described.
Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT

Gaussian Process Regression and the Enforcement of Additional Physical Constraints
Gaussian process regression has recently seen increased interest as a tool for building surrogate models for describing material behaviors (Rocha et al. (2021); Fuhg et al. (2022b); Chen et al. ( 2023)).This is in part because GPR is derived from a convenient statistical background, offers strict convergence guarantees (unlike, for example, artificial neural networks), and has empirically shown excellent performance for out-of-sample model predictions (Park et al. (2014); Datta et al. (2016)).For an excellent overview of GPR, the reader is referred to Rasmussen (2003).In this section, the basic ideas behind GPR and its specific employment for regression in our surrogate models are discussed.Two distinct but related GPR formulations are introduced: the standard GPR and the constrained GPR.The latter allows enforcement of thermodynamic consistency in the viscous overstress component of our data-driven model.

Standard Gaussian Process Regression (GPR)
Consider a general dataset D g of N data points, where x ∈ R ni and y ∈ R no represent input and output data points of dimensions n i and n o , respectively.The sets of all input and output data points can be reformulated as matrices X ∈ R N ×ni and Y ∈ R N ×no , respectively.The objective of a GPR model that is trained using the dataset D g is to predict an output y for a new input data point x that is not contained in the training dataset.
In GPR, the notion of the similarity between input data points is critical.The assumption is that two input data points that are closer together are more likely to have closer target values than the two input data points that are farther away from each other.Typically, the similarity between two input data points x and x is modeled by a user-defined covariance function.The choice of this function is an important part of GPR.Based on the work of Laurent et al. (2019), which compared different functional forms of the covariance function for computer experiments when no prior knowledge is available, we restrict ourselves to the Matérn 3/2 kernel (Matérn, 1960), given as where (•) 2 represents the L2-norm, σ f is a scaling factor, and l is the characteristic length-scale of the covariance function.σ f and l can be combined into a vector θ = [σ f , l] T that collects the trainable parameters.The term α denotes a small positive value that helps to overcome numerical instabilities Peng and Wu (2014).δ x,x is the Kronecker delta function.Using the functional form of the covariance function in Eq. ( 36), a covariance matrix K(X, X ) can be constructed to define similarity between sets of points X and X .
The optimal set of parameters θ that best describes the training dataset (Eq.( 35)) is obtained using a maximum log-likelihood approach (Fuhg et al., 2021b), After finding the best parameters, the GPR regression model is fully defined.Given a new input data point x , the predicted output value y of the Gaussian process regressor reads As α → 0 in Eq. ( 36), the GPR predictor becomes an exact interpolator (see Marrel et al. (2008)), i.e. y(x i ) = y i for all points i = 1, . . ., N of the training dataset of Eq. ( 35).This is termed as the exact inference property of GPR.Since the prediction is probabilistic, the variance σ 2 of the regressor fit can also be obtained, In this study, standard GPR (as defined above) is employed to fit the surrogate models M vol (Eq.( 26) and M h,iso (Eq.( 27) using datasets of the form of D * vol and D * h,iso (Eq.( 25), respectively; remember, the datasets D * vol and D * h,iso are first derived from the initial training datasets D vol and D h,iso (Eq.( 21), respectively (see Algorithm 1).Further, the exact inference property of GPR is harnessed to enforce the physical constraint of a stress-free reference (i.e., Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT undeformed) state, also called the normalization condition (Holzapfel (2000)).For the two stress components S v,iso and S h,iso that are modeled using standard GPR, the normalization condition reads To enforce the above constraint, the stress-free reference states-{C = I, S vol = 0} and {C = I, S h,iso = 0}-are included in the training dataset.This guarantees that the normalization condition is achieved, provided the stabilization constant α of Eq. ( 36) is kept small enough.Note, in terms of invariants, C = I translates to Ī1 = Ī2 = 3.
Another physics-based constraint of interest in this work is the second law of thermodynamics.In light of the Clausius-Planck inequality, this constraint requires that the internal dissipation Ξ int in a visco-hyperelastic soft material as described in Eq. ( 7) is non-negative.Among the three stress components S vol , S h,iso and S v,iso , the first two are rate-independent and therefore do not contribute to dissipation.Thus, enforcing the second law of thermodynamics constraint on the surrogate models M vol and M h,iso that are based on standard GPR becomes trivial.In other words, the specific construction of the generalized constitutive equation in this work automatically ensures the second law constraint on the rate-independent portion of our data-driven constitutive model.

Constrained Gaussian Process Regression (C-GPR)
Unlike the rate-independent stress components S vol and S h,iso , the rate-dependent isochoric viscous overstress S v,iso results in viscous dissipation.In this regard, the Clausius-Planck inequality in Eq. ( 7) can be rewritten as (Pioletti and Rakotomanana (2000)) For arbitrary input strain and strain rates, the isochoric viscous overstress prediction of the surrogate model M v,iso (Eq.( 34)) must satisfy the above inequality.This type of inequality constraint on input-output mappings can be enforced by employing the constrained GPR (C-GPR) formulation recently proposed by Pensoneault et al. (2020).In standard GPR, the model output remains unconstrained.
The basic idea of C-GPR is to restrict the solution space in the hyperparameter optimization problem (i.e., Eq. ( 37)) such that the predicted output at a set of user-defined "constraint points" follows the desired inequality constraint.
Accordingly, in the present study, the following optimization problem is devised for training the GPR regressor of the surrogate model M v,iso (cf.Eq. ( 37)): Here, S v,iso (C k , Ċk ) denotes the isochoric viscous overstress prediction of the surrogate at the k th constraint point (i.e., for the {C k , Ċk } input).N c is the total number of constraint points.The constrained optimization of Eq. ( 42) imposes the Clausius-Planck inequality (Eq.( 41)) on the rate-dependent part of our physics-informed data-driven constitutive model.Note that the inequality condition in Eq. ( 42) is not a functional constraint; rather, it is applied at a finite set of input data points.Therefore, C-GPR only applies a "weak" constraint on the model output.
Next, the normalization condition constraint for the isochoric viscous overstress is written as The above condition requires that the isochoric viscous overstress vanishes in the undeformed state (i.e., C = I) regardless of the applied rate of deformation Ċ in that state.Owing to the particular choice of the strain rate invariants J1 , J4 , and J6 in the mapping of the M v,iso surrogate (Eq.( 31)), this condition can be simply enforced by including the stress-free reference state-{C = I, S h,iso = 0}-in the training data.This is possible because in the undeformed state, these three invariants hold a fixed value of zero regardless of the loading rate, i.e., For a derivation of the above result, see Appendix B.. From Eq. ( 44), C = I translates to a single set of input values in the mapping of M v,iso : [ Ī1 , Ī2 , J1 , J4 , J6 ] = [3, 3, 0, 0, 0], at which the exact inference property of GPR will always lead to a zero viscous overstress value.
Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT Unlike the invariants J1 , J4 , and J6 , the invariants J2 , J3 , J5 , and J7 in the undeformed state of C = I are functions of the tensor Ċ (see Appendix B.).Therefore, if these invariants were included in the mapping of M v,iso (Eq.( 31)), there would be infinitely many possible combinations of invariant values (i.e., input data points) that would correspond to the stress-free reference state, making it impossible to utilize the GPR exact inference property for ensuring the normalization condition.

Model Performance
In this section, the performance of our physics-informed data-driven constitutive model is demonstrated in terms of (i) the accuracy of fitting training datasets from various deformation modes (e.g., hydrostatic compression, and unconfined uniaxial tension and compression), (ii) the accuracy and physical plausibility of out-of-sample predictions outside the training regime, and (iii) the training dataset size requirement.As the three surrogate models that make up our constitutive model-M vol , M h,iso , and M v,iso -are independent of each other in terms of the specific dataset type required for training them (i.e., D vol , D vol , and D vol in Eq. ( 21), respectively), the following subsections will evaluate these surrogate models individually using multiple deformation modes for training and testing purposes, for a deeper insight.In practice, given multiple datasets spanning volumetric, isochoric quasi-static, and isochoric dynamic deformation modes, the individually trained M vol , M h,iso , and M v,iso surrogate models can be conveniently combined to form one constitutive model.
To assess the fitting accuracy of the surrogate models, the following percent relative error metric is employed: where (•) F represents the Frobenius norm, S is the stress tensor predicted by a surrogate model, and S is the true value of the stress tensor (i.e., the ground truth).The scalar error, err, is obtained with respect to individual data points.
For convenience, a percent mean error value will also be reported, where err i is the error value at the ith data point, and N is the total size of the dataset (training or testing).
To assess the accuracy and physical plausibility of model predictions outside the training regime, two types of out-ofsample "testing regions" will be considered: (i) Testing region that corresponds to the same deformation mode that was considered during training-for example, training a model using uniaxial tension stress-strain data in the [0, 0.25] strain range, and then predicting tensile stresses for the wider [0, 0.5] strain range; (ii) Testing region that corresponds to a different deformation mode from those considered during model training-for example, training a model using uniaxial tension stress-strain data and then using the trained model to predict material responses under uniaxial compression and simple shear deformation modes.Further, the effect of training dataset size on both fitting errors and out-of-sample predictions will be studied.
Finally, for each surrogate model case, the performance of our model will be compared to (i) a conventional viscohyperelastic constitutive model, and to (ii) a corresponding surrogate model that employs a classical ML black-box mapping between the Voigt forms of tensors C and Ċ (as strain and strain rate input) and S (as the second Piola-Kirchoff stress output) (Ghaboussi et al. (1998); Lefik and Schrefler (2003)), i.e., Classical mapping: Standard GPR will be utilized to learn the classical mappings.Of course, tensor Ċ will not be considered in the case of quasi-static loading conditions.
All the calculations presented in this work are performed in Python (Van Rossum and Drake ( 2009)).For machine learning tasks, the Scikit-learn python package (Pedregosa et  Consider the following deformation gradient tensor, Deformations of the form of Eq. ( 48) are applied in practice via confined compression experiments to study the hydrostatic bulk material response (Sasson et al. (2012); Kim et al. (2019); Nie et al. (2013)).These experiments provide hydrostatic stress / pressure versus volumetric strain (i.e., J) response of the material, which in turn can be used to calibrate a particular form of the volumetric energy density function U (J) (see Eq. ( 14)).In this study, the commonly used Simo-Miehe volumetric energy density (Simo and Miehe (1992)) is employed to generate the hydrostatic stresses for training the M vol surrogate model, where κ is the bulk modulus.We choose κ = 10.The corresponding volumetric stress component S SM vol (see Eq. ( 13)) is given by where C = C app = F T app F app .C app and S SM vol constitute the training dataset D vol of Eq. ( 21) (note, D = D vol in this case) that is utilized to train the M vol surrogate model (Eq.( 26)) by following the procedure detailed in Section 3. Standard GPR is used for mapping.As the stress-free undeformed state (i.e., F app = I) is included in the training data, the normalization condition is simply enforced by assigning a near-zero Gaussian noise of α = 10 −4 (Eq.( 36)) at all data points.
The confined compression deformation mode in Eq. ( 47) generates three non-zero volumetric stress components: S vol,11 , S vol,22 , and S vol,33 .Given the 11-loading direction, S vol,22 = S vol,33 , resulting in only two independent stress components.Figure 2(a) shows the evolution of S vol,11 and S vol,22 versus J from the training data and compares it with the corresponding S vol,11 -J and S vol,22 -J responses predicted by the trained surrogate model.26 data points were considered for training in this case (i.e., a strain increment of 0.01).The corresponding percent relative error (Eq.( 45)) versus J plot is shown in Fig. 2(b).From these plots, a very good agreement between training data and surrogate model predictions in the training regime is evident, which is expected from an ML-based model.Overall, the maximum and mean relative errors in the training regime are 1.12% and 0.12%, respectively, suggesting a very good fitting accuracy of our surrogate model.
Next, the performance of our surrogate model for a wider range of J (i.e., the testing regime) is analyzed and compared with the corresponding predictions of (i) a conventional volumetric strain energy density function of the neo-Hookean  47)).The neo-Hookean volumetric strain energy density and its corresponding volumetric stress tensor are given as The calibration of the volumetric neo-Hookean model against the training data yielded a κ NH value of 11.56, which is close to the ground truth of κ =10 considered in the Simo-Miehe model for generating the training data.
Figures 3(a-b) show the S vol,11 -J and S vol,22 -J responses as predicted by the surrogate model proposed in this work, the volumetric neo-Hookean model, and the classical mapping of Eq. ( 47), in a strain range (J ∈ [0.5, 1]) that outstretches the training regime in the compression deformation mode (i.e., the deformation mode considered in training).The ground truth (i.e., Simo-Miehe Model) in this extended strain regime is also shown in the same figure for comparison, with the corresponding percent relative error err (with respect to the ground truth) versus J plots for the three models shown in the left half of Fig. 3(e).From these figures, it is clear that within the training regime (i.e., J ∈ [0.75, 1]), all three models exhibit reasonable accuracy in predicting the stress-strain data.Among the three models, the two ML-based data-driven models -the surrogate model of this study and the classical mapping model -resulted in much smaller percent mean relative errors (err) of 0.12% and 0.27% compared to the volumetric neo-Hookan model, which resulted in a err of 4.90%.This was expected because unlike ML-based mappings, conventional constitutive models have limited fitting accuracy due to their fixed mathematical forms.Outside the training regime when J ∈ [0.5, 0.75], our surrogate model clearly outperforms the other two models by predicting stress components that are in excellent agreement with the ground truth.In terms of the percent mean relative error, our surrogate model yields a err of 1.07% in this regime, which is an order-of-magnitude smaller than the corresponding errors of the volumetric neo-Hookean model (i.e., 14.86%) and the classical mapping model (i.e., 12.72%).Figures 3(c-d) compare the predictions of the three models with the ground truth (Simo-Miehe model) in the tension (confined) deformation mode (J ∈ [1, 1.5]) that was not considered at all in the training data.Here, the classical mapping model results in physically implausible predictions in that the stresses at large volumetric strains start to decrease, violating thermodynamic stability for incremental deformations.In fact, for J > 1.6 (not shown in the plots), this model predicts compressive stresses for a confined tension loading, which is not possible.This type of erroneous model behavior is attributed to the purely black-box mapping in these conventional ML-based constitutive models, which do not respect any physical or mechanistic constraints on the response of the continuum.From Figs.
3(c-d), our physics-informed data-driven surrogate model does not suffer from this issue and makes physically-plausible and trustworthy predictions even in the deformation mode that was not considered in training.Not surprisingly, the predictions of the volumetric neo-Hookean model are also physically reasonable (as κ NH > 0).Further, our surrogate model outperforms the other two models in terms of agreement of the predicted stress components with the ground truth.
The percent mean relative error of our model in the tensile regime is 6.66%, which is significantly lower compared to the errors of the volumetric neo-Hookean model (i.e., 23.59 %) and the classical mapping model (i.e., 37.23%).
Lastly, the effect of training dataset size on the performance of the two data-driven models is analyzed in Figs.4( )), which includes both confined compression and tension deformation modes, the classical mapping model no longer shows improvement in the performance with an increase in training dataset size.In fact, for N vol > 50, err versus N vol increases monotonically, which coincides with highly physically implausible predictions in the confined tension regime (see an example of this in the supplementary material, Section SM1).Our surrogate model does not suffer from this limitation due to its physics-informed nature, and its predictions continue to improve with an increase in training dataset size in both the training and testing regimes.Further, as the inclusion of physics-based constraints confines the feature space of possible GPR hyperparameters, our surrogate model predictions are in very good agreement with the ground truth even at N vol = 26.Clearly, our physics-informed data-driven constitutive model combines the high fitting accuracy of ML-based models with the physical nature and small-data compatibility of conventional constitutive models.

5.2
The M h,iso Surrogate Model Under Quasi-static Uniaxial and Shear Loading Now, consider the following isochoric (volume-preserving) uniaxial deformation mode, where λ is called the uniaxial stretch.Under the quasi-static assumption, the rate of stretching λ is approximately zero.Deformations like Eq. ( 52) are routinely applied in the mechanical characterization of soft materials during where A 10 and A 01 are model parameters.We choose A 10 = 1 and A 01 = 0.5 for generating the training data.Also, J = J app = det(F app ), C = C app = F T app F app , and C = Capp = J −2/3 C app .Using C app and S MR h,iso as the training dataset D h,iso (in this case, D = D h,iso ), a standard GPR-based mapping as in Eq. ( 27) was trained.Like the hydrostatic loading case, the normalization condition, in this case, is automatically enforced owing to the presence of the stress-free reference state in the training data.This trained surrogate model can predict the isochoric hyperelastic stress component S h,iso for any arbitrary quasi-static deformation.
Figure 5 compares the stress-stretch response considered for training the surrogate model with the corresponding predictions of the trained surrogate model in this regime.Like the hydrostatic case (Section 5.1), there are two independent stress components resulting from the uniaxial loading condition, S h,iso,11 and S h,iso,22 , with S h,iso,22 = S h,iso,33 .26 data points were considered for training in this case.From Fig. 5(a), the model predictions are in excellent agreement with the training data, leading to relatively small errors as shown in Fig. 5(b).The mean percent relative error in this case is 4.75%.
With a good fitting accuracy established from Fig. 5, the performance of our surrogate model is now analyzed in a wider testing regime.The testing regime in this case consists of uniaxial deformation mode (Eq.( 52)) with λ ∈ [0.5, 1.5], and the simple shear deformation mode, where γ is the shear strain.The surrogate model performance is compared with that of the classical mapping model (Eq.( 47)) and another polynomial hyperelastic model, the 2-parameter Yeoh model (henceforth, referred to simply as the Yeoh model).The Yeoh model (Yeoh (1993)) is given by  In the overall testing regime, our surrogate model results in a markedly better performance than the other two models both in terms of the physical plausibility of predicted responses and the prediction accuracy.For example, the compressive stress-stretch responses predicted by the classical mapping model (Fig. 6(c-d)) exhibit unreasonable mechanical features: (i) the S h,iso,11 − λ response shows a softening response that is uncharacteristic of hyperelastic soft materials, which typically exhibit compression-tension asymmetry in the loading direction (i.e., the material exhibits a highly nonlinear and stiffer response in compression than in tension) (Budday et al. (2017); Upadhyay et al. (2020a);Treloar (1944)), and (ii) the S h,iso,22 − λ response starts to decrease following a maximum at approximately λ = 0.57, which violates the expected monotonicity in stress-stretch responses (prior to damage/failure).The latter behavior is also exhibited by the Yeoh model, for which the predicted S h,iso,22 versus λ response violates monotonicity at large compressive strains of λ < 0.61.Non-physical response predictions from conventional hyperelastic models have been observed in the literature for certain model parameter values and are attributed to the violation of the second law of thermodynamics by the model with those parameter values (Liu (2012); Upadhyay et al. (2019b); Fontenele et al. (2022)).Owing to its physics-informed formulation, our surrogate model prevents this issue and yields a physically-plausible mechanical response even at large deformations.Further, it results in a mean percent relative error of 10.98% in the uniaxial testing regime, which is slightly lower compared to the err of 11.08% of the Yeoh model, and significantly lower than the err of 49.94% of the classical mapping model.
Figure 7(a) shows the predicted shear stress versus shear strain responses from the three models and the ground truth.
Here, the classical mapping model predicts zero stress regardless of the applied strain.This is attributed to the sole consideration of uniaxial deformation in the training dataset, in which non-diagonal components in the C and S h,iso tensors were always zero.As our surrogate model is based on the mapping of strain invariants with response coefficients, it does not suffer from this limitation and leads to accurate shear stress predictions as shown in Fig. 7(a).Note that even though the simple shear deformation mode generates multiple non-zero stress components (viz., S h,iso,11 , S h,iso,12 , S h,iso,22 , and S h,iso,33 ), only the dominant 12-component is shown for brevity.In terms of the percent relative error metric that includes all non-zero stress components and is plotted in Fig. 7(b), our surrogate model leads to a maximum err of 13.46% across the shear deformation mode, with a mean error of 10.09%.For comparison, the maximum err for the Yeoh model is 14.59%, while the mean err is 7.54%.Overall, it is seen that both the surrogate model of this study and the Yeoh model exhibit reasonable prediction accuracy across the entire testing regime (uniaxial + shear), but only our surrogate model results in stress-strain predictions that are physically plausible across all the investigated deformation modes.
Figure 8 shows the evolution of mean percent relative error err in the training and overall testing regimes, with the training dataset size N h,iso .Similar to the hydrostatic loading case, our surrogate model results in a monotonically decreasing err with N h,iso response with asymptotic behavior in both cases.On the other hand, the classical mapping model shows an asymptotic decrease in err only in the testing regime, while showing a slight increase in the training regime.Ultimately, while both the models show very good fitting accuracy for every investigated training dataset size, our surrogate model results in significantly more accurate predictions across multiple deformation modes.This ability to seamlessly transition between deformation modes without loss of accuracy is a distinctive feature of trustworthy constitutive models (Beda (2014)).Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT

The M v,iso Surrogate Model Under Dynamic Uniaxial and Shear Loading
Thus far, the performance of the physics-informed data-driven constitutive model of this study has been evaluated only under quasi-static loading conditions, when the applied strain rate is so small that viscous dissipation remains negligible.Now, consider an isochoric dynamic uniaxial loading at a finite applied loading rate, equivalently following where λ = dλ/dt is the applied uniaxial strain rate.At any given time t during loading, λ = λt.
Dynamic loading of the form of Eq. ( 56) is applied during high strain rate tension tests on soft materials such as the tensile Kolsky bar experiment (e.g., see Upadhyay et al. (2021); Yang et al. (2000)).The stress-strain (or stretch)-strain rate data from these tests combined with the stress-strain (or stretch) data from quasi-static uniaxial tests are then used to calibrate visco-hyperelastic constitutive models.As a number of visco-hyperelastic constitutive models (see Upadhyay et al. (2020b) for a review) assume an additive decomposition of total stress into hyperelastic (rate-independent) and viscous overstress components (e.g., see Eq. ( 14)), the quasi-static test data is exclusively used to calibrate the hyperelastic model ( Wh ) parameters, while the viscous overstress that is calculated by subtracting quasi-static stress-strain data from high strain rate stress-strain data is used to calibrate the model parameters of the viscoelastic component of the constitutive model (e.g., Wv in Eq. ( 14)).In this study, the 3-parameter Upadhyay-Subhash-Spearot (USS) viscous dissipation potential (Upadhyay et al. (2020b(Upadhyay et al. ( , 2022a))) is utilized to generate the viscous overstress training data for developing the M v,iso surrogate model, where k 11 and k 21 are the linear and nonlinear rate sensitivity control parameters, respectively, and c 21 is the rate sensitivity index.We use k 11 = 1, k 21 = 1, and c 21 = 0.75 for generating the training data.The expressions for invariants Ī1 , Ī2 , J2 , and J5 are given in Eq. ( 15).Just like the quasi-static loading case, J = J app = det(F app ), C = C app = F T app F app , and C = Capp = J −2/3 C app .Ċ = Ċapp is the time rate of change of Capp .
The training dataset D v,iso = D = {C app , Ċapp , S USS v,iso } (see Eq. ( 21)) generated from the USS model is used to train the M v,iso surrogate model (Eq.( 28)).Similar to the previous quasi-static cases (involving S vol and S h,iso ), the normalization condition is implicitly enforced due to the consideration of the stress-free reference state in the training data.The second law of thermodynamics constraint, on the other hand, is enforced via C-GPR as the chosen regression method.Here, the inequality constraint on the hyperparameter optimization problem in Eq. ( 42) is enforced at all training input data points (i.e., the constraint points are the same as the training data points; N c = N v,iso ).The trained surrogate model can predict viscous overstress S v,iso for any given applied tensors C and Ċ.
Figure 9 investigates the fitting response of the data-driven surrogate model in the training regime of uniaxial tension loading (Eq.( 52)) that is composed of 5 uniformly distributed stretch rates in the range of λ ∈ [10, 100] (viz., λ = 10, 32.5, 55, 77.5, and 100), and 31 uniformly distributed stretch values in the range of λ ∈ [1, 1.5] for every stretch rate considered.The total number of training data points is thus 5 × 31 = 155.Both the non-zero stresscomponents under uniaxial tension, S v,iso,11 and S v,iso,22 , are plotted in Figs.9(a) and 9(b), respectively.An excellent agreement between model predictions and the stress-stretch-stretch rate training data is apparent from these plots.The corresponding scalar fitting error (err) versus stretch response is shown in Fig. 9(c); here, the line represents the mean err values averaged across all the investigated stretch rates, and the shaded region represents the vertical range of mean ± standard deviation.The err versus λ response starts at relatively high values at small stretches and then assumes consistently small values close to zero for λ > 1.1.The large err at small stretch values is attributed to the very small stresses in the small strain regime in the training data (i.e., the denominator in the formula for err), which causes spuriously high err even when the absolute differences between true and predicted stresses are very small.Overall, the mean percent relative error err across all the 155 training data points is 1.73%, suggesting a very good model fitting accuracy.
The accuracy and physical plausibility of our data-driven model predictions are now investigated in a wider testing regime that consists of tension, compression, and simple shear deformation modes.Here, the two uniaxial deformation modes (i.e., tension and compression) comprise stretch rates in the range of λ ∈ [−145, 145] and stretch values in the range of λ ∈ [0.5, 1.75] (cf.Eq. ( 56)).Specific stretch rate values considered in model testing are: λ = ±10, ±32.5, ±55, ±77.5, ±100, ±122.5, and ±145.Now, the dynamic simple shear deformation is given by where γ = dγ/dt is the applied shear strain rate.At any given time t during loading, γ = γt.For model testing, we consider shear strain rates in the range of γ ∈ [10, 145] (specifically, γ = 10, 32.5, 55, 77.5, 100, 122.5, and 145); for each strain rate, shear strains in the range of γ ∈ [0, 0.5] are considered.The overall testing regime outstretches the training regime both in terms of the maximum stretch/strain and stretch-/strain-rate magnitudes and spans multiple deformation modes.Like the previous two quasi-static cases, the present model predictions are compared with those from the classical mapping model (Eq.( 28)) and an existing constitutive model from the literature.Here, we choose the Pioletti model (Pioletti et al. (1998); Pioletti and Rakotomanana (2000)), The Pioletti model was calibrated using the tensile stress-stretch-stretch rate training dataset, resulting in η = 6.94.
Figures 10(a-b) compare the tensile stress versus stretch responses in the testing regime at three representative stretch rates (for clarity) as predicted by the present data-driven surrogate model, the classical mapping model, and the Pioletti model, with the ground truth (i.e., USS model predictions).Out of the three stretch rates shown in these figures, λ = 10 and λ = 55 belong to the training subset, while λ = 145 belongs to the testing dataset that was not considered during training.In the training regime (i.e., λ ∈ [1, 1.5] and λ ∈ [10, 100]), both the data-driven models exhibit excellent fitting performance as their predicted tensile stress-stretch responses nearly overlap the ground truth.The Pioletti model, on the other hand, results in a relatively poor fitting accuracy, which is attributed to its overly simple mathematical form with only one model parameter.In terms of the scalar percent error metric plotted in Fig. 11 (averaged across all investigated strain rates), the average prediction error of our surrogate model in the training regime is the lowest at 1.73%, followed by the classical mapping model regime in the tension deformation mode when λ > 1.5 or λ > 100 (see Figs. 10a-b), the Pioletti model remains the worst performing model among the three investigated models.This time, however, the stress predictions of the classical mapping model also differ considerably from the ground truth, especially at high stretch rates.Overall, in the tensile deformation mode, the mean percent relative error, err, of our surrogate model is just 3.40%, which is considerably lower than that of the classical mapping model (12.30%) and the Pioletti model (41.59%).
In the compression deformation mode (Figs.10(c-f)), the stress predictions of all three models show disagreement with the ground truth to varying degrees.Nevertheless, the stress-stretch plots of our data-driven surrogate model and the Pioletti model are physically reasonable and exhibit the expected compressive nature (i.e., negative S v,iso,11 and S v,iso,22 ), monotonicity at every investigated stretch rate, and compression-tension asymmetry in the 11-loading direction (i.e., stiffer response in compression than in tension).The classical mapping model, on the other hand, predicts tensile stresses under compression (see Figs. 10 (e-f)), which is physically implausible.Owing to these unreasonable predictions, the mean percent relative error err of this model, 253.14%, is an order of magnitude higher than the corresponding mean errors of our surrogate model (39.69%) and the Pioletti model (45.59%).
The predictions of the two data-driven models, the present surrogate model and the classical mapping model, are now analyzed in light of their conformity to the second law of thermodynamics in the testing regime.In this regard, Fig. 12 plots the internal viscous dissipation Ξ int as a function of stretch at multiple stretch rates as predicted by these two models.In the tension deformation mode, both the models result in non-negative viscous dissipation, thus satisfying the second law of thermodynamics (see Eq. ( 41).However, in compression, the classical mapping model predictions show a negative internal dissipation at every investigated stretch rate, which violates the second law of thermodynamics.This causes the model to predict physically unreasonable stress-stretch-stretch rate predictions as were seen in Figs.10(e-f).
Owing to its physics-informed construction that enforces the non-negativity of viscous dissipation, our surrogate model prevents this issue and its predicted uniaxial mechanical response conforms to the second law of thermodynamics throughout the tensile testing regime.
The performance of our surrogate model in the simple shear deformation mode is analyzed and compared with the Pioletti and classical mapping models in Fig. 13.Like the quasi-static hyperelastic case of Section 5.2, the classical model predicts a zero shear stress component (i.e., S v,iso,12 ) regardless of the applied shear strain and strain rates (see Fig. 13(a)).This non-physical behavior is attributed to the absence of any non-zero non-diagonal components  , it is clear that among the two data-driven constitutive models under comparison, only the surrogate model of this work results in trustworthy predictions that maintain physical plausibility even under large strain and at high strain rates and more importantly at deformation modes completely unseen by the model, confirming its ability to generalize proficiently.Fig. 13(c) plots the internal viscous dissipation versus shear strain at multiple shear strain rates from these models.Unsurprisingly, the model predictions of our surrogate model show non-negative dissipation, thus obeying the second law of thermodynamics.Similar to the compression deformation mode, the classical mapping model in shear deformation exhibits negative dissipation that is inconsistent with the second law of thermodynamics.
Finally, the effect of training dataset size N v,iso on the mean percent relative error err of the two data-driven models is studied in Fig. 14.Here, N v,iso is the product of the number of distinct tensile stretch rates in the range of λ ∈ [10, 100] considered during training and the number of stretch values in the range of λ ∈ [1, 1.5] for every stretch rate case (e.g., in Fig. 9, N v,iso was equal to 5 × 31 = 155).We consider nine N v,iso values: 155 (5 (stretch rates) × 31 (stretch values)), 186 (6 × 31), 217 (7 × 31), 305 (5 × 61), 366 (6 × 61), 427 (7 × 61), 455 (5 × 91), 546 (6 × 91), and 637 (7 × 91).The effect of N v,iso on err values computed for the training regime (i.e., the fitting error) and for the overall testing regime (tension + compression + shear) are shown in Fig. 14(a) and Fig. 14(b), respectively.From these plots, the mean percent relative errors of both the data-driven models in the training regime are very small (maximum err value of less than 5%), demonstrating a consistently high fitting accuracy.Further, both the err versus N v,iso responses exhibit a generally decreasing trend.In the testing regime, the responses of the two data-driven models differ significantly.On one hand, the classical mapping model leads to extremely high mean percent relative errors that generally increase with the training dataset size.On the other hand, our data-driven model consistently results in reasonably small mean percent relative errors (< 20%) for all investigated training dataset sizes, revealing a good prediction performance even when data availability is limited.

Summary and Discussion
This work presents a framework for the development of physics-informed data-driven constitutive models to describe the short-time, strain-rate-dependent mechanical response of soft materials.A major motivation of this work is the limitations of the traditional continuum thermodynamics-based and the machine learning-based constitutive models: while the former have a limited fitting and prediction accuracy owing to their fixed mathematical form and require expert model selection based on experimental observations (e.g., choice of the hyperelastic model), the latter require Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT exorbitant amounts of training data (experiments that would be required to uniformly sample strain space are not physically plausible) and generally result in poor out-of-sample predictions (poor generalization performance).Our proposed framework takes a significant step toward eliminating these limitations by combining the physics-informed nature of continuum thermodynamics with the highly accurate and flexible regression capability of supervised ML.The result is a fully data-driven constitutive model that can capture complex material response features with high accuracy without any expert intervention, yields physically reasonable and accurate out-of-sample predictions, and can be trained with a small amount of training data that is achievable from simple contemporary experiments.As some of these points where previously exhibited for hyperelasticity and elastoplasticity this work outlines the extension of this general framework toward capturing the rate-dependent response of viscoelastic materials at finite deformations.
The formulation of our data-driven constitutive model is based on the generalized stress-strain-strain rate equations of the continuum thermodynamics-based framework of external state-variable driven viscous dissipation-based viscohyperelasticity.In these equations, the total stress is additively decomposed into volumetric, isochoric hyperelastic, and isochoric viscous overstress components.Each of these stress components is written as linear combinations of the components of an integrity basis.This type of linear relationship allowed us to propose three data-driven surrogate model mappings-M vol , M h,iso , and M v,iso -to capture each of the three stress components, respectively.These surrogate models map strain / strain rate invariants to the coefficients of the integrity basis that make up their corresponding stress components.It is shown that this type of model construction ensures key physics-based constraints on the predicted response: principles of local action, determinism, material frame-indifference, the balance of angular momentum, isotropic material symmetry, and limited memory.Further, owing to the exact inference property of the GPR supervised learning method (both standard and C-GPR) and the special inequality constraint capability of C-GPR, the proposed surrogate models also respect the normalization condition and the second law of thermodynamics.
The performance of each of the three surrogate models that form our constitutive model was evaluated by fitting them to a small numerically-generated training dataset, each obtained from one deformation mode-corresponding to common experimental protocols-and then applying the trained model to predict material responses in a significantly wider testing regime comprising multiple distinct deformation modes.In every case, our model's performance was compared with those of a traditional continuum thermodynamics-based constitutive model and a classical ML-based mapping model.The results showed that our models provide critical improvements in describing material responses compared to the other modeling frameworks.For example, from the results of the M vol model under hydrostatic loading, it was seen that our surrogate model predictions are in excellent agreement with the ground truth both in the training regime of confined compression (J ∈ [0.75, 1]) and in the overall testing regime comprising confined compression and tension (J ∈ [0.5, 1.5]).In comparison, the traditional volumetric neo-Hookean model resulted in large prediction errors in the testing regime, and the classical mapping model resulted in large errors as well as a physically unreasonable softening response in confined tension.
The pathology of predicting unreasonable / thermodynamically unstable physical responses in deformation modes that are not considered in training was consistently noted from the classical mapping approach throughout this study.For example, in the isochoric quasi-static loading case that was considered to evaluate the M h,iso surrogate model, the classical mapping model trained under uniaxial tension violated stress-stretch monotonicity and compression-tension asymmetry under uniaxial compression, and resulted in a zero shear stress regardless of the applied deformation under simple shear loading.In the isochoric dynamic loading case considered for analyzing the M v,iso surrogate model, the classical mapping model trained under rate-dependent uniaxial tension predicted tensile stresses under compression as well as no shear stresses under simple shear.Notably, most efforts in ML-enabled constitutive models in the literature follow the so-called classical approach of mapping strain components to stress components, which is evidently insufficient for performing tasks that traditional phenomenological and micromechanical modeling excelled at.
Note that in the isochoric quasi-static loading case, even the traditional Yeoh model resulted in a physically unreasonable softening response in compression.Physically unreasonable model predictions are known to be caused by thermodynamic consistency issues (Upadhyay et al. (2019b)).Unsurprisingly, in the isochoric dynamic loading case, it was seen that the classical mapping model predictions violated the second law of thermodynamics through negative internal dissipation in both the compression and shear deformation modes (the two deformation modes not considered in training).Unlike these alternative model types, our physics-informed data-driven surrogate models resulted in physically-reasonable mechanical responses across all the investigated deformation modes (training and testing) and their prediction errors were consistently at the lowest level (among the three model types).
The physics-informed nature of our constitutive model not only promotes accurate and physically reasonable predictions but also eliminates the requirement of large training datasets by restricting the solution space of possible surrogate model parameters based on physical laws (Karniadakis et al. (2021) seen by simple experiments.For each of the three surrogate models of this work, a reasonably good prediction accuracy over the full testing regime was obtained with relatively small-sized training datasets like the ones commonly seen in the experimental literature.In addition, the prediction accuracy generally improved when the training dataset size was increased (notably, still in the limited-data regime, utilizing synthetic data in deformation modes consistent with easily accessible experiments).On the other hand, classical ML-based mapping models resulted in very high prediction errors at low training dataset sizes, and these errors did not necessarily decrease with increasing data volume.In fact, in the hydrostatic and the isochoric dynamic loading cases, the prediction errors of these black-box models increased when large datasets were employed for training them, which coincided with highly non-physical predictions in the testing regime.Overall, even the lowest mean percent relative error of the classical mapping models, in any case, was several times the peak mean relative error of our surrogate models.
Recall that the proposed data-driven constitutive model is only applicable to short-time responses that are dominant under high strain rate deformation.Under sufficiently slow loading rates, long-time effects such as relaxation and creep become important.In addition, our current model enforces material isotropy as one of the physics-based constraints.While many soft materials (e.g., hydrogels, tissues, and elastomers) are considered isotropic, some are anisotropic (e.g., the white matter of the brain (Eskandari et al. (2021))).In future work, efforts to add long-time effects and anisotropy modeling capabilities in the current framework will be pursued, as well as the utilization of direct experimental data.
Appendix A. Derivation of the generalized stress equations of the visco-hyperelastic constitutive model Starting with Eq. ( 14), the total stress in a visco-hyperelastic soft material is written as The individual stress components in Eq. ( 62) can be simplified using the following results 2 , The following standard results are used to derive the second part of Eq. ( 63): Given a scalar φ and second order tensors Y and Z, we have ∂(αY) ∂Z = α ∂Y ∂Z + Y ⊗ ∂α ∂Z .Ċ = d(J −2/3 C) dt = J −2/3 Ċ − 2 3 J −5/3 JC, where J = J 2 tr(F −T ĊF −1 ) Physics-informed Data-driven Discovery of Constitutive Models with Application to Strain-Rate-sensitive Soft Materials A PREPRINT P is called the referential fourth order projection tensor (Holzapfel (2000)), , I being the fourth order unit tensor (in index notation, I ijkl = δ ik δ jl , where δ is the Kronecker delta symbol).The tensor P projects any second order tensor to its referential deviatoric component, i.e., P : Z = DevZ, where Z is an arbritrary second order tensor.Remember, DevZ = Z − 1 3 (Z : C)C −1 .
By substituting Eq. ( 63) in Eq. ( 62) and invoking the projection property of P, we have The derivatives of invariants Īj with respect to C and those of invariants Jk with respect to C are given as Note, the result C2 = Ī1 C − Ī2 I + C−1 is due to the Cayley-Hamilton theorem (Holzapfel (2000)).By substituting Eq. (66) in Eq. ( 65), we get S = S vol + S h,iso + S v,iso From Eq. ( 67), the three stress-components that constitute the total stress in a visco-hyperelastic soft material can be expressed as  Equations (68-70) constitute the alternative form of the generalized functional constitutive equation in Eq. ( 14) and is utilized in the physics-informed mapping approach of this study.Of course, traditional constitutive models have explicit expressions for coefficients ζ 1 , Γ 1 , Γ 2 , Φ 1 , . . ., Φ 7 that are functions of their model parameters.Table A1 lists out these coefficients for several commonly used visco-hyperelastic constitutive models.Our physics-informed data-driven model is not limited by the fitting ability of such explicit mathematical expressions and can learn the mapping between invariants and the coefficients directly from the experimental data.Table A.2: Coefficients of the integrity bases in commonly employed hyperelastic strain energy density functions ( Wh ): the neo-Hookean model (Rivlin (1948a);Treloar (1943)), the Mooney-Rivlin model and the generalized Rivlin models (Mooney (1940); Rivlin (1948b)), the Yoeh model (Yeoh (1993)), the Gent model (Gent (1996)), and the Gent-Gent model (Pucci and Saccomandi (2002))

Figure 1 :
Figure 1: Schematic diagram showing the different stages in the evolution of constitutive modeling.

Figure 2 :
Figure 2: (a) Comparison of the numerically generated volumetric stress components (S vol,11 and S vol,22 ) versus volumetric strain (J) training data under confined compression with the corresponding surrogate model predictions.Inset shows a schematic illustration of the confined compression deformation mode, highlighting the reference and deformed states and the 11-loading direction.(b) Evolution of the percent relative error (err) of surrogate model predictions versus J in the training regime.

Figure 3 :
Figure 3: Comparison of the numerically generated volumetric stress-volumetric strain data from the Simo-Miehe model (i.e., ground truth) in the overall testing regime (J ∈ [0.5, 1.5]) with the corresponding predictions of our surrogate model, the volumetric neo-Hookean model, and the classical ML-based mapping model: (a) S vol,11 versus J under confined compression, (b) S vol,22 versus J under confined compression, (c) S vol,11 versus J under confined tension, and (d) S vol,22 versus J under confined tension.(d) Comparison of the percent relative error (err) versus J responses for the present surrogate model, the volumetric neo-Hookean model, and the classical mapping model.

Figure 4 :
Figure 4: (a) Comparison of the evolution of mean percent relative error (err) in the predictions of the surrogate model and the classical mapping model in the training regime (J ∈ [0.75, 1]), as a function of the training dataset size (N vol ).(b) The corresponding err versus N vol responses of the present surrogate model and the classical mapping model, for their predicted responses in the overall testing regime (J ∈ [0.5, 1.5]).
a-b), which plot the mean percent relative errors err in the training regime (J ∈ [0.75, 1]) and the overall testing regime (J ∈ [0.5, 1.5]) as a function of the training dataset size N vol (i.e., the number of data points in D vol ).In the training regime (Fig.(4(a))), the errors of both the models decrease asymptotically with N vol .This is a typical behavior of ML-based models, which yield improved predictions in the training regime (i.e., interpolation) as larger volumes of data are used in model training.Overall, both the models show excellent fitting accuracy (err) even with a small training dataset of 26 stress-strain values.In the testing regime(Fig.(4(b)

Figure 5 :
Figure 5: (a) Comparison of the numerically generated isochoric hyperelastic stress components (S h,iso,11 and S h,iso,22 ) versus uniaxial stretch (λ) training data under uniaxial tension with the corresponding surrogate model predictions.Inset shows a schematic illustration of the unconfined uniaxial tension deformation mode, highlighting the reference and deformed states and the 11-loading direction.(b) Evolution of the percent relative error (err) of surrogate model predictions versus λ in the training regime.

Figure 6 :
Figure 6: Comparison of the numerically generated isochoric hyperelastic stress-uniaxial stretch data from the Mooney-Rivlin model (i.e., ground truth) in the uniaxial testing regime (λ ∈ [0.5, 1.5]) with the corresponding predictions of our surrogate model, the Yeoh model, and the classical mapping model: (a) S h,iso,11 versus λ under uniaxial tension, (b) S h,iso,22 versus λ under uniaxial tension, (c) S h,iso,11 versus λ under uniaxial compression, and (d) S h,iso,22 versus λ under uniaxial compression.(e) Comparison of the percent relative error (err) versus λ responses of the present surrogate model, the Yeoh model, and the classical mapping model.

Figure 6
Figure6compares the stress versus uniaxial stretch predictions of the present surrogate model, the classical mapping model and the Yeoh model, with the ground truth (i.e., from the Mooney-Rivlin model that was utilized for generating training data).In the training regime of λ ∈ [1, 1.25] (see Fig.6(a-b)), all three models show an excellent agreement with the ground truth.In terms of the percent relative error that is shown in Fig.6(e), the mean errors in the three model predictions in this regime are 4.75% (surrogate model), 1.26% (classical mapping model), and 0.50% (Yeoh model).

Figure 7 :
Figure 7: (a) Comparison of the numerically generated isochoric hyperelastic shear stress (S h,iso,12 )-shear strain (γ) data from the Mooney-Rivlin model (i.e., ground truth) in the shear testing regime (γ ∈ [0, 0.5]) with the corresponding predictions of our surrogate model, the Yeoh model, and the classical mapping model.(b) Comparison of the percent relative error (err) versus γ responses of the present surrogate model, the Yeoh model, and the classical mapping model.

Figure 8 :
Figure 8: (a) Comparison of the evolution of mean percent relative error (err) in the predictions of the surrogate model and the classical mapping model in the training regime (λ ∈ [1, 1.25]), as a function of the training dataset size (N h,iso ).(b) The corresponding err versus N h,iso responses of the present surrogate model and the classical mapping model, for their predicted responses in the overall testing regime (λ ∈ [0.5, 1.5] ∪ γ ∈ [0, 0.5]).

Figure 9 :
Figure 9: (a) Comparison of the numerically generated training data of isochoric viscous overstress component Sv,iso,11 versus uniaxial stretch (λ) at multiple stretch rates ( λ), with the corresponding surrogate model predictions.Inset shows a schematic illustration of the unconfined uniaxial tension deformation mode.(b) The corresponding Sv,iso,11 versus λ responses at fixed λ values from the training data and the surrogate model.(c) Evolution of the percent relative error (err) of surrogate model predictions versus λ in the training regime (dashed line: mean err across all stretch rates; shaded region: mean ± standard deviation).

Figure 10 :
Figures10(a-b) compare the tensile stress versus stretch responses in the testing regime at three representative stretch rates (for clarity) as predicted by the present data-driven surrogate model, the classical mapping model, and the Pioletti model, with the ground truth (i.e., USS model predictions).Out of the three stretch rates shown in these figures, λ = 10 and λ = 55 belong to the training subset, while λ = 145 belongs to the testing dataset that was not considered during training.In the training regime (i.e., λ ∈ [1, 1.5] and λ ∈ [10, 100]), both the data-driven models exhibit excellent fitting performance as their predicted tensile stress-stretch responses nearly overlap the ground truth.The Pioletti model, on the other hand, results in a relatively poor fitting accuracy, which is attributed to its overly simple mathematical form with only one model parameter.In terms of the scalar percent error metric plotted in Fig.11(averaged across all investigated strain rates), the average prediction error of our surrogate model in the training regime is the lowest at 1.73%, followed by the classical mapping model (4.67%) and the Pioletti model (40.59%).Outside the training

Figure 11 :Figure 12 :
Figure 11: Comparison of the percent relative error (err) versus λ responses of the present surrogate model, the Pioletti model, and the classical mapping model.The lines represent mean err across all investigated strain rates, and the shaded regions represent mean ± standard deviation.

Figure 13 :
Figure 13: (a) Comparison of the numerically generated isochoric shear viscous overstress (Sv,iso,12)-shear strain (γ) data from the USS model (i.e., ground truth) in the shear testing regime (γ ∈ [0, 0.5], γ ∈ [10, 145]) with the corresponding predictions of our surrogate model, the Pioletti model, and the classical mapping model.(b) The corresponding percent relative error (err) versus γ responses of these three models.The lines represent mean err across all investigated strain rates, and the shaded regions represent mean ± standard deviation.(c) Comparison of the internal dissipation (Ξint) versus γ responses of the two data-driven models: our surrogate model and the classical mapping model.Note: in (a) and (c), the background light gray lines are the data corresponding to the strain rates γ = ±32.5,±77.5, ±100, and ±122.5, which were also evaluated in model testing but are not shown for clarity.