Strain mode-dependent weighting functions in hyperelasticity accounting for verification, validation, and stability of material parameters

Optimized material parameters obtained from parameter identification for verification wrt a certain loading scenario are amenable to two deficiencies: Firstly, they may lack a general validity for different loading scenarios. Secondly, they may be prone to instability, such that a small perturbation of experimental data may ensue a large perturbation for the material parameters. This paper presents a framework for extension of hyperelastic models for rubber-like materials accounting for both deficiencies. To this end, an additive decomposition of the strain energy function is assumed into a sum of weighted strain mode related quantities. We propose a practical guide for model development accounting for the criteria of verification, validation and stability by means of the strain mode-dependent weighting functions and techniques of model reduction. The approach is successfully applied for 13 hyperelastic models with regard to the classical experimental data on vulcanized rubber published by Treloar (Trans Faraday Soc 40:59–70, 1944), showing both excellent fitting capabilties and stable material parameters.


State of the art on hyperelasticity
Rubber-like materials or elastomers, respectively, consist of randomly oriented chain-like macromolecules with more or less closely connected entanglements or cross-links. Two main characteristics are their ability for large deformations subject to relatively small stresses and their retaining of the initial configuration after unloading without considerable permanent deformation. This behavior is attributed to the network entropy as the orientation of chains alters with deformation. Due to these special properties, the materials have numerous technical application such as for tires, structural bearings, medical devices and base isolations of buildings; see, e.g., [7]. Numerous phenomenological and micro-mechanically motivated models have been proposed in the literature in order to capture the elastic and nearly incompressible mechanical of rubber-like materials. Moreover, the former can be classified into invariant-based and principal-stretch-based formulations, cf. e.g. [47].
Phenomenological invariant-based models are based on the theory of invariants as elaborated extensively by Spencer in [46] for anisotropic materials. For the case of isotropy, an appropriate set of invariants dependent on the right Cauchy-Green tensor are selected, which are included as polynomials with sufficiently high orders into the strain energy function, known as Rivlin's expansion, [43]. Classical examples such as the Neo-Hooke material are of Mooney-Rivlin type, cf., e.g., [33,[40][41][42]. Since then a vast variety of models have been developed. For example, Yeoh [55] revealed that cubic terms of the first invariant are able to reproduce the highly nonlinear S-shaped uniaxial behavior of rubber, also at very large strains. Alternatively to polynomial formulations, logarithmic formulations are presented, e.g., in the Gent-Thomas model [15], the Gent model [14] or the Pucci Saccomandi model [39]. In Khajehsaeid et al. [22], the combination of polynomials, logarithmic, and exponential formulations is investigated. Another approach is given in the model of Carrol [9], which combines powers of the first invariant with the square root of the second invariant. A prominent example for a phenomenological principal-stretch-based model is presented by Ogden in [37].
Micro-mechanical models are derived from statistical mechanics arguments on networks of idealized chain molecules. In this way, they account for a lower scale insight to the physical/chemical microstructure to explain phenomenological macroscopic mechanical behavior, although this might render higher computational costs.
Well known examples for micro-mechanical models are the 3-chain, the 4-chain, and the 8-chain model as well as the unit sphere (21-chain) model, cf. Arruda and Boyce [3] and Miehe et al. [32].
Many of the above-mentioned models for rubber-like materials are well advanced from the mathematical point of view, [17], and the numerical point of view [31,45].

Criteria for a "good" parameter vector
A well accepted first step for parameter identification is based on a least-squares functional, in order to minimize the distance of simulated data and experimental data with respect to a chosen norm. However, this step leaves open the precise meaning of a "good" vector of material parameters. In this work, we follow the outline in [27] accounting for the three criteria of Verification is related to the fit quality of a material model as a consequence of parameter identification. This step of model development assesses if the simulated outputs of the underlying model are adequate when compared to measurement data which contribute to the underlying least-squares functional. In this way, all measurement data involved in this step are known beforehand, such that more precisely it could also be termed backwards-verification. Comparative studies on verification for models in hyperelasticity have been performed for example in Marckmann and Verron [30], Boyce and Arruda [6], Seibert and Schöche [44]. Parameter identification for Rivlin and Saunders type hyperelasticity [43] has been examined, e.g., by Hartmann et al. [16]. The identification from inhomogeneous experiments on polyurethane is investigated in [53], where particular focus has been set on error-controlled adaptive mesh refinement.
Some aspects of validation and stability are outlined next.

An aspect of validation: Conformity of data sets and constitutive model
In contrast to verification, validation of material models is related to the prediction capability. Alternatively, it could also be termed forecast-verification, since this step of model development assesses whether the simulated outputs of the underlying model are adequate when compared to additional measurement data which are not known or considered, respectively, in the backwards-verification step; see, e.g., [27].
In this work, we are particularly interested in possible remedies in case of non-validity for material models in hyperelasticity. In particular, we will resort to a typical situation of laboratory testing, where a certain number of data sets each related to a certain loading scenario is available. An example for such a set of experimental data is provided by Treloar [48] for the three loading scenarios of uniaxial, biaxial, and pure shear deformations. An extensive comparative study on the fit quality of 14 material models on these data is provided by Steinmann et al. in [47]. Most of the models meet the criteria of verification, however, merely wrt to one individual loading scenario with a certain vector of optimized material parameters (occasionally with perfect agreement). Then, the same parameter vector is not able to capture different loading scenarios (occasionally even showing disappointing agreement). This lack of validity motivates the following definition: Definition 1 (on conformity between experimental data and constitutive model) Available experimental data sets, all of them related to different loading scenarios (or stress modes, strain modes, respectively), and a constitutive model are conform, if a reduced experimental data set is sufficient to validate the material model for the remaining data sets.
Definition 1 applied to the extensive comparative study in Steinmann et al. in [47] renders, e.g., a poor conformity between the Treloar data and the Isihara model [20], while the situation improves significantly for the Carroll model [9].
It should be emphasized that Definition 1 on conformity is restricted to available data and, in this sense, constitutes only a necessary but not a sufficient condition for general validation, that is, predictive simulation for all loading cases that could be envisaged. Examples for further possible loading scenarios are hydrostatic tension, uniaxial compression, equibiaxial compression, or hydrostatic compression (although experimentally challenging) in order to get a more comprehensive (though in general still not complete) characterization of the material. 1 The experimental data sets in Definition 1 may have different origins. In practice, typically they may refer to different loading scenarios, which allows us to exploit the extensive literature for modeling asymmetric effects within the fields of plasticity and creep. Along this line constitutive equations have been formulated, e.g., in [2,5,25,26,51,[56][57][58], among others. It appears from the above mentioned references that so far no common approach exists concerning the best strategy for taking account of individual loading scenarios in the constitutive equations. A general agreement is the incorporation of odd power terms for odd invariants of stresses; see, e.g., [5]. Furthermore, a scalar variable, which is expressed in terms of the ratio of the second and third basic invariant of the deviatoric stress tensor, can be used as an indicator for detection of differences in the loading modes. This quantity, stress mode angle or Lode-factor, respectively, has been applied, e.g., in [1,12,57], [25]. Based on [26], extensive use of the stress mode angle has also been made by this author and co-workers for modeling asymmetric effects of experimental data in tension, compression and shear for metallic materials (such as AISI 52100) [10,11,29], with applications for cutting processes like hard turning [49].
Aims of present work on non-conformity A main purpose of this work is the extension of some well-known hyperelastic material models in case they fail to be conform to given data sets according to Definition 1. To this end, we make extensive use of the above-mentioned works on asymmetric effects. In particular, as a counterpart to the stress mode angle used in [1,12,57], [25,26] a strain mode angle is introduced in order to characterize the individual loading scenario. The key idea consists in an additive decomposition of the logarithmic isochoric Hencky strain tensor, where each of the related quantities incorporates a weighting function dependent on the strain mode angle. The advantage of this approach is, that certain (though not all) material parameters, such as Mooney-Rivlin-type constants, can be obtained individually from specific loading modes such as uniaxial tension, equibiaxial tension and pure shear, investigated experimentally in the laboratory.

An aspect of stable material parameters: I-criteria
The solution of the underlying least-squares problem for parameter-identification might render a satisfactory agreement between simulated and experimental data; however, it might be susceptible to instability, in the sense, that a small perturbation of experimental data may lead to a large deviation of the resulting parameter solution. Two possible reasons for this undesirable effect (or even non-uniqueness) have to be distinguished: 1. Deficiency of the material model, where (too many) functional terms and/or parameters may be lead to (almost) linear dependencies within the model (overparametrization) or 2. Deficiency of experimental data, where certain material effects intended by the model are not properly activated (incomplete data). see, e.g., [27].
For detection of possible instabilities several indicators are proposed in the literature. A well-known example is the correlation matrix, which is defined, e.g., in [38]. In addition, an indicator for perturbations of the measurements is given in [50]. Alternatively, statistical methods can be considered as discussed, e.g., in [27]. Since the generation of experimental data might become costly, this approach can be supported by stochastic simulation, cf., e.g., [35,36].
In the field of optimum experimental design a confidence matrix is introduced, which typically may be constructed in terms of the Jacobian of the underlying least-squares functional. Then, different indicator functions for evaluation of a stable solution vector are defined, known as A-,E-,D-and M-criterion and subsequently generally denoted as I-criteria; see, e.g., [4] and [24] for precise definitions. In the field of optimal control problems, the indicator function may be dependent on further design variables such that an optimality problem can be formulated, cf., e.g., [4] and [24].

Aims of present work on stable material parameters
Following [50], we perform a first-order perturbation analysis in order to investigate the influence of perturbed experimental data to the perturbation of material parameters. This analysis motivates the definition of four so-called I-criteria known from optimal control problems. In this work, these criteria will be applied, to give a stability assessment of material parameters related to the strain mode decomposition discussed above, which eventually can be exploited for model reduction.

Structure of this work
The structure of the paper is as follows: Based on the preliminaries in Sect. 2 for hyperelasticity in continuum mechanics, Sect. 3 presents a general framework for strain mode-dependent hyperelasticity. To this end, strain mode related weighting functions are introduced, which are incorporated into a new formulation of a general strain energy function. We propose a general framework for convenient implementation, including the consistent tangent modulus, for any free energy density in terms of the principal isochoric stretches, the eigenvalues and respectively the invariants of the isochoric right Cauchy-Green strain tensor. Section 4 considers aspects of stable parameter identification. We perform a perturbation analysis to motivate four different so-called I-criteria known from optimal control theory. In Sect. 5 we propose a practical guide for model development accounting for the criteria of verification, validation and stability by means of the strain mode-dependent weighting functions. In Sect. 6 the approach is applied for 13 hyperelastic models with regard to the classical experimental data on vulcanized rubber published by Treloar [48]. Detailed investigations with a focus on verification, validation and stability by means of the strain mode-dependent weighting functions and techniques of model reduction will be presented.

Notations
Square brackets [•] are used throughout the paper to denote 'function of' in order to distinguish from mathematical groupings with parenthesis (•).

Preliminaries of hyperelasticty in continuum mechanics
This section provides a brief review on the modeling of hyperelasticty in the finite strain regime of continuum mechanics. Particular interest is directed to basic kinematics, spectral decomposition, and the derivation of stress tensors and tangent operators associated to a given strain energy density.
We focus on hyperelastic properties for rubber-like materials which exhibit a decoupled response to volumetric and isochoric deformations. For this purpose, the following kinematical quantities are indispensable Here the material gradient F is defined as the partial derivative of the nonlinear deformation map ϕ with respect to the position vector X ∈ R 3 in the Euclidean space R 3 , J is its determinant, C is the right Cauchy-Green tensor and C is the isochoric right Cauchy-Green tensor, cf. e.g. [18,19]. The decomposition (1.4) can also be derived from the multiplicative split F = (J 1/3 1) ·F of the deformation gradient that goes back to Flory [13] and satisfies the incompressibility condition detF = 1. Here, also the second-order unity tensor 1 has been used.
In the subsequent exposition, extensive use will be made of the following relations; see, e.g., [45]: Here Eq. (2.1) represents the spectral decomposition of the right Cauchy-Green tensor C in Eq. (1.3) with associated eigenvalues a and eigenvectors N a , a = 1, 2, 3. Analogously, Eq. (2.3) represents the spectral decomposition of the isochoric right Cauchy-Green tensor C in Eq. (1.4) with associated eigenvalues¯ a and eigenvectors N a , and where the eigenvalues a and¯ a , a = 1, 2, 3 are related by Eq. (2.4). According to standard notation, see, e.g., [19], we introduce the three principal stretches in terms of the eigenvalues a and¯ a , respectively. Note, that λ a and λ a are the eigenvalues of F andF, respectively. Moreover, three invariants I 1 := tr[C], I 2 := 1 2 tr[C] 2 − tr[C 2 ] and I 3 := det(C) of the right Cauchy-Green strain tensor C, and analogously defined for its isochoric part C are written in terms of the eigenvalues a and¯ a , respectively, as In the last relation the condition det[F]) = λ 1 λ 2 λ 3 =J = 1 for the isochoric incompressible material behavior has been exploited. In view of Eq. (1.4), the above invariants are related as A well accepted starting point for modelling hyperelastic materials is an additive decomposition of the strain energy function into volumetric (shape preserving) and isochoric (volume preserving but shape changing) parts: Additionally, in Eq. (6) we have accounted for a vector of material parameters κ ∈ K ⊂ R n p , where K ⊂ R n p denotes the n p -dimensional space of admissible parameters, which will be investigated more deeply in the ensuing Sect. 4. The second Piola-Kirchhoff stress tensor S is given as the derivative of the free energy wrt C multiplied by 2, that is S = 2∂ψ/∂C, and a further derivative wrt C multiplied by a factor 2 renders the fourth-order elasticity tensor (or tangent operator) C = 2∂S/∂C = 4∂ 2 ψ/∂C∂C, cf. e.g. [45]. Exploiting the additive decomposition Eq. (6) yields the corresponding decoupled stress tensor and tangent operator as By use of the chain rule, the volumetric stress tensor can be expressed aŝ Here p is known as the hydrostatic pressure and the result ∂ J /∂C = 1/2C −1 has been used, cf. e.g. [45].

Strain mode-dependent energy function
Many existing hyperelastic models known from the literature are well advanced both from the mathematical and the numerical point of view and offer a convincing fit quality for a variety of nonlinear experimental stress-strain data. However, in many cases, optimized material parameters obtained for a certain strain mode (e.g., uniaxial tension) lack a general validity for different strain modes (e.g., equibiaxial tension or pure shear), and consequently are not conform according to Definition 1 in Sect. 1.3. This section presents a framework for extension of hyperelastic models for rubber-like materials exhibiting different behaviors in different loading scenarios, which can be examined individually in the laboratory. The starting point for the strain mode related approach is the following isochoric part of the strain energy ψ in Eq. (6): The above mathematical structure is identical to the one in [26] for creep simulation of asymmetric effects by use of stress mode-dependent weighting functions. Equation (9.1) represents an additive decomposition of the energy function into S strain mode related quantities. Each of them incorporates a strain mode-dependent energy function ψ i dependent on the isochoric right Cauchy-Green tensor C, and a vector of material parameters κ i associated to the ith strain mode investigated individually in the laboratory, such that κ = {κ i } 3 i=1 . Furthermore, in the above skeleton structure (9) a weighting function w i is associated to each mode i, which is dependent on a strain tensor E(C), and for which it is stipulated that where δ i j is the Kronecker-delta. Also, the strain tensors E(C j ), j = 1, 2, .., S refer to independent characteristic strain modes, which for example can be investigated experimentally in uniaxial tension, equibiaxial tension and pure shear, respectively. Note, that Eq. (10.1) can be regarded as a completeness condition, whereas Eq. (10.2) constitutes a normalization condition for the weighting functions. We remark also, that above Eq. (9) are restricted to isotropic materials. The aspect of anisotropy combined to strain mode-dependent material behavior is not within the scope of this paper.

Choices for strain mode related weighting functions
For illustrative purpose, we consider three independent strain modes for uniaxial tension (UT), equibiaxial tension (ET) and pure shear (PS), respectively, as published by Treloar [48] for experimental data on vulcanized rubber. For UT, only one out of three principal stretches λ a , a = 1, 2, 3 is prescribed, λ 1 = λ say, for ET there are two prescribed values λ 1 = λ 2 = λ say, and for PS we require λ 1 = λ and λ 2 = 1 say. From the condition on incompressibilityĪ 3 = λ 2 1 λ 2 2 λ 2 3 = 1, and the assumption of isotropy, the complementary principal stretches follow accordingly. In summary, the corresponding deformation gradients and the right Cauchy-Green tensors are, cf. e.g. [47]: The above tensors are formulated in terms of powers of λ. This motivates the following choice for the strain tensor E in Eq. (9) as the Hencky strain tensor in logarithmic form where M aa is the dyadic basis in Eq. (2.2), and where the relation (3.2) has been used. Application of Eq. (14) to the three loading scenarios in Eq. (11) to Eq. (13) renders  With these properties at hand, we can follow the same procedure as extensively outlined by [12] for deviatoric stress tensors. The three modes in Eq. (15) are represented in the octahedral plane of the associated logarithmic isochoric strain tensor. To this end the quantities are defined. Here θ shall be referred to as the strain mode angle dependent on the strain mode factor ξ . Furthermore, in above Eq. (16.3) J 2 and J 3 denote the second and third basic invariant of the logarithmic isochoric strain tensor E, respectively. Exploiting the spectral decomposition of the logarithmic Hencky strain tensor in Eq. (14) renders A graphical interpretation of the strain mode angle θ is given at the top of Fig. 1. In particular, it becomes apparent, that the independent strain modes of UT, ET, and PS, respectively, are characterized by the strain mode angles and for each of the strain modes the following periodicity angles are obtained. Based on these observations, in addition to the relations (10) the following is required for the weighting functions For the strain modes related to the loading scenarios of uniaxial tension, equibiaxial tension and pure shear, respectively, we set S = 3 and the requirements (20) are satisfied by the following weighting function and where n = 0, 1, 2, ... are integer values. A graphical representation of the weighting functions (21) is given in the middle graph of Fig. 1. For the case, that experimental data are available only for loading in uniaxial tension and equibiaxial tension, respectively, with S = 2, the following weighting functions can be used These functions, which are illustrated at the bottom of Fig. 1, do also satisfy the requirements (20). Upon using the definition (16.1), alternatively the functional relationships (21) and (22) can be rewritten in terms of the strain mode factor ξ as 1. UT: w (1) respectively, which is more convenient for numerical implementation of the strain energy ψ in Eq. (9).

Hyperelasticity based on principal stretches
A general formulation for the strain energy function for isotropy may be written as ψ = ψ[λ a , κ], which we assume as continuously differentiable with respect to the three principal stretches λ a , a = 1, 2, 3 of Eq. (3.1).
In order to account for incompressibility and strain mode-dependent material behavior, the free energy density for the isochoric part in Eq. (9) is now postulated in a mixed formulation as in terms of the eigenvalues¯ a introduced in Eq. (2) and the principal isochoric stretches λ a , a = 1, 2, 3 of Eq. (3.2). Closed-form expressions for the related stresses and tangent moduli of ψ i [λ a , κ i ] in terms of the reference configuration as well as the current configuration have been derived in [45]. A summary for the isochoric second Piola-Kirchhoff stress tensor S and the corresponding tangent C of Eq. (7) is provided, e.g., in [47]. A reformulation of the mixed relation (25) purely in terms of the eigenvalues a requires a reparametrization of the individual energy functions from a stretch based formulation ψ i [λ a , κ i ] to an eigenvalue based , which is easily achieved by means of the relations (3). Accordingly, straightforward differentiation renders by means of the chain rule the isochoric second Piola-Kirchhoff stress tensor in Eq. (7.1) as 1. 3. Setting S = 1 and w i = const in the above general structure (25) ensues W i a = 0 for the stress like coefficients in Eq. (26.5), such that the stress tensor S boils down to classical results as, e.g., in [45].
Similarly, the isochoric tangent operator in Eq. (7.2) is derived as where the fourth-order basis tensors 1.

Hyperelasticity based on invariants
Alternatively to the principal-stretch-based formulation (25) a strain energy function can be written depending on the invariants of the right Cauchy-Green tensor. A general formulation reads ψ = ψ[I 1 (C), I 2 (C), I 3 (C)], which we assume as continuously differentiable with respect to the three invariants I A , A = 1, 2, 2 in Eq. (4.1) of the right Cauchy-Green strain tensor. In order to account for incompressibility and strain modedependent material behavior, the free energy density for the isochoric part in Eq. (9) is now postulated in a mixed formulation as in terms of the eigenvalues¯ a introduced in Eq. (2) and the two invariantsĪ A , A = 1, 2 in Eq. (4.2) of the isochoric right Cauchy-Green strain tensor. Closed-form expressions for the related stresses and tangent moduli of ψ i [Ī A , κ i ] have been derived, e.g., in [21] for a formulation based on invariants, see also [31]. A summary for the second Piola-Kirchhoff stress tensor and the corresponding tangent of Eq. (7) is provided e.g. in [47]. A possible formulation related to the ith strain mode for the energy function in Eq. (25) is given as the Mooney-Rivlin/Saunders model [43] where the vector of material parameters associated with each strain mode is defined as A reformulation of the mixed relation (30) purely in terms of the eigenvalues a requires a reparametrization of the individual energy functions from an invariant based formulation The coefficients ∂ψ i /∂Ī A must be obtained for the individual material model. Examples for some classical material models are provided in Appendix B. Moreover, the following results are required for evaluation of Eq. (33): and where ∂¯ b /∂ a is given according to Eq. (27.3). Analogously, the results for the corresponding tangent in Eq. (28) are reparametrized from an invariant based formulation to an eigenvalue based formulation. The tangent coefficients C i ab in Eq. (28.4) become The coefficients ∂ 2 ψ i /∂Ī A ∂Ī B must be obtained for the individual material model. Examples for some classical material models are provided in Appendix B. The coefficients ∂Ī A /∂¯ a are given in Eq. (34). Moreover, the following results are required for evaluation of Eq. (35): 1 else (36) and where the coefficients ∂ 2Ī A /∂¯ a ∂¯ b are obtained based on (34). Observe the analogous structure of the results (29.1) and (36.1).

Least-squares minimization problem (inverse problem)
We denote byd ∈ D = R n d a measurement vector of experimental data obtained from laboratory testing. The corresponding vector of n d simulated data shall be termed d [κ], where as before κ ∈ K denotes the vector of material parameters of the underlying material model. In general n d ≥ n p , such that the distance between both data types must be minimized by means of a suitable least-squares function f : D × K → R, e.g., in the following simple form: Occasionally, problem (37) will also be referred as the inverse problem. The first-order necessary and the second-order sufficient optimality condition, respectively, are Variational formulations for determination of M are provided e.g. in [28,50]. Details on solution strategies for the minimization problem (37) are given elsewhere and shall not be considered here, cf. e.g. [27]. At this stage we only point out, that the Jacobian J and the Gauss-Newton matrix H G N play key roles in the performance of different solution strategies. The solution κ * of problem (37) might render a satisfactory agreement between simulated and experimental data; however, it might be susceptible to instability, in the sense, that a small perturbation of experimental data, δd say, may lead to a large deviation δκ for κ * . Possible reasons for this undesirable effect [or even nonuniqueness of problem (37)] are two-fold, [27]: 1. Deficiency of the material model: Functional terms and/or parameters may lead to (almost) linear dependencies within the model (overparametrization), or 2. Deficiency of experimental data: Certain material effects intended by the model are not properly "activated" (incomplete data).
In the sequel, two different indicators are formulated in order to quantify stability for the solution κ * of the minimization problem (37): 1. Stability indicator by correlation matrix 2. Stability indicators by I-criteria.

Stability indicator by correlation matrix
Following, e.g., [38] the coefficients of the correlation matrix are given as The correlation coefficient −1 ≤ C i j ≤ 1 is a quantitative measure for the correlation or respectively linear dependence between parameter κ i and κ j . Positive correlations occur for C i j ≥ 0 where an increase of κ i results in an increase of κ j . Negative correlations with C i j < 0 lead to decreasing κ j for increasing κ i and vice versa. There are no correlations when C i j = 0. If |C i j | = 1, i = j, the correlation is called perfect and ensues Regarding the results of parameter identification |C i j | 1 , i = j is desirable for a stable-or respectively a robust-solution vector κ * .

First-order perturbation analysis
In order to perform a first-order perturbation analysis, following [50], let us define a function F : K × D → K. The first-order necessary condition Eq. (38) for a given measurement vectord becomes The implicit function theorem guaranties the existence of a neighborhood D 0 ⊂ D and a continuously differentiable functionκ : Here, H and J are the Hessian and the Jacobian in Eq. (39), respectively, and from Eq. (42) With the result in (44.2) Eq. (45) becomes Several possibilities can be envisaged to obtain an estimate from Eq. (46); see, e.g., [50]. Neglecting the second-order term O ||δd|| 2 D and taking a selected norm on K, we obtain Here and where μ i are the eigenvalues of the matrix C. Eq. (48) reveals that the eigenvalue μ max has the interpretation of an amplification factor for the weighted perturbation of data J T δd to the perturbation of parameters δκ.

Stability indicators by I-criteria
In addition to the eigenvalue μ max defined in Eq. (49), several amplification factors have been introduced in the field of optimum experimental design, subsequently denoted as I-criteria with general structure; see, e.g., [4] and [24] Common examples of φ I are, see, e.g., [4] and [24] 1. A-Criterion: 3. D-Criterion: Thus, the E-criterion in Eq. (51.2) is identical to the criterion in Eq. (49). Figure 2 provides geometrical interpretations of all four I-criteria based on the confidence ellipsoid for n p = 2, cf., e.g., [23,52]. The Acriterion is proportional to the average length of the confidence ellipsoid, the D-criterion to its volume, the E-criterion to its the largest expansion, and the M-criterion to the largest side length of a box around the confidence ellipsoid. The A-criterion in Eq. (51.1) is related to the E-criterion in Eq. (49.2) as A relation between the D-criterion in Eq. (51.3) and the E-criterion in Eq. (49.2) is  1. In the field of optimal control problems, the indicator function φ I [C[κ]] may be dependent on further design variables such that an optimality problem can be formulated, cf. e.g. [4] and [24]. Then, in summary, two requirements are formulated: The solution vector κ * (a) should minimize the least squares function (representing a model error) according to (37) (b) should minimize the confidence criterion according to (50). 2. In this work, the generality of the methodology of optimal control problems will not be exploited. Instead, the different criteria in (51) will simply be used, to give an assessment of existing solutions κ * for different well known formulations of the energy density function ψ[C, κ]. The treatment of parameter identification as an optimal design problem dependent on additional design variables will constitute an aspect of future work. To the knowledge of the author, there is no similar work so far which applies A-,E-,D-, and M-criteria in the field of hyperelasticity. 3. A mathematical correct definition of stability is given e.g. in [50]. From there, we point out that local stability guarantees local uniqueness but not global uniqueness of the inverse problem (37).

A practical guide for model development
In order to account for the criteria of verification, validation and stability by means of the strain mode-dependent weighting functions of Sect. 3, a practical guide for model development is summarized in Table 1. As input for the resulting flowchart, we assume that measurement datad i ∈ D i = R n i d and an initial set of energy functions ψ i , both related to strain modes i = 1, ..., S, are given.
Step 1: Verification for each mode data set In this step of Table 1 each set of experimental datad i ∈ D i = R n i d , i = 1, ..., S is used separately to minimize mode data related least squares functionals of the form (37) to obtain Step 2: Stability for each mode data set In order to account for mode data related stability of all solution vectors κ i , based on the general structure in Eq. (50), in Step 2 of Table 1 we evaluate (at least one of) the following I-criteria: Here we use p according to Eq. (39.1) applied to all modes i = 1, ..., S.
Step 3: Validation for complete data set In order to check conformity according to Definition 1, in Step 3 of Table 1 each solutionκ i of Step 2 is used to simulate the remaining modes i = 1, ..., S, j = i. This defines the following complete data related least-squares functional Note, that f Step 4: Stability for complete data set In order to account for complete data related stability in Step 4 of Table 1 based on the general structure in Eq. (50), we evaluate (at least one of) the following I-criteria Here we use C = [J T J ] −1 with Jacobian J ∈ R n d ×n p according to Eq. (39.1), and where, contrary to Eq. (55), n d = S i=1 n i d refers to the complete set of available experimental data.
Step 5: Selection of final parameter vector In Step 5 of Table 1 the final parameter vector is selected as Consequently, if this step is reached, there is no need for the strain mode-dependent formulation Eq. (9).
Step 6: Weighted mode-related model extension

In case Step 3 on validation and/or
Step 4 on stability fail, Step 6 in Table 1 is activated, which makes use of the strain mode-dependent formulation Eq. (9) with verified and stabilized mode-related material parameterŝ κ i as a result of Step 1 and Step 2.  Table 1 render verified and stabilized material parametersκ i related to each mode data setd i , i = 1, ..., S.

Provided
Step 3 and Step 4 in Table 1

F[d, κ]}
Evaluate: to obtain the final parameter vector. 5. Since the strain mode-dependent formulation (9) in Step 6 in Table 1

Selection of experimental data and hyperelastic models
In this section, the experimental data on vulcanized rubber published by Treloar [48] are used to demonstrate the steps for model development according to Table 1. Three loading scenarios, that is uniaxial tension (UT), equibiaxial tension (ET) compression and pure torsion (PS), respectively, are considered. The data in [48] are given in pairs of principal stretches λ k and principal nominal stresses P T reloar [λ T reloar k ] for all three loading cases.
Regarding the hyperelastic models, we follow very closely the selection provided by Steinmann et al. in [47]. Here, both phenomenological and micro-mechanically motivated models have been investigated in order to capture the elastic and nearly incompressible effects of rubber-like materials investigated in [48]. The complete list of hyperelastic models investigated in this paper is as follows: Regarding the extensive literature on the above hyperelastic models we refer also to [47], and therefore shall not be repeated here. As mentioned before, the corresponding isochoric free energy densities related to each strain mode ψ i in Eq.

Analytical stress formulations for UT, ET, and PS
Assuming an isotropic and incompressible material, an analytical P i a [λ a ] formulation is provided, e.g., in [19,47,54] as Our implementation departs slightly from the one in [47] by exploiting a generality for UT, ET, and PS as follows: • Uniaxial tension (UT): For given deformation gradient according to Eq. (11), the stresses normal to the load directions are zero, that is P UT 2 = P UT 3 = 0. By setting Eq. (60) to zero, for example for a = 3, and inserting the resulting pressure into Eq. (60) for a = 1, renders • Equibiaxial tension (ET): For given deformation gradient according to Eq. (12), here the stresses in the load directions are equal, that is P ET 2 = P ET 3 , while the third direction is stress-free. Therefore, by setting Eq. (60) to zero for a = 3, and inserting the resulting pressure into Eq. (60) for a = 1 and a = 2, renders Observe the same rule in Eq. (61) and in Eq. (62), which is due to the above choice a = 3 for UT. • Pure shear (PS): For given deformation gradient according to Eq. (13) the stress in the third direction is zero. Therefore, by setting Eq. (60) to zero for a = 3, and inserting the resulting pressure into Eq. (60) for a = 1 and a = 2, renders for i = U T, E T, P S. The partial derivatives ∂ψ i /∂λ a in Eq. (60) for principle stretch-based formulations, and the partial derivatives ∂ψ i /∂Ī A in Eq. (64) for invariant-based formulations, respectively, are summarized for thirteen hyperelastic models in compact form in Appendix B of this work. In our implementation Eq. (60) is used for the final stress calculation for all three load-cases (UT, ET and PS), by accounting for the case-dependent pressure p. It applies throughout for all phenomenological (invariant-based and stretch-based formulation) and micro-mechanically motivated models of the above list with thirteen models. This unified procedure might be at the expense of some numerical effort, but avoids the implementation of load case-dependent and material-dependent formulations as provided in [47]. With this framework, the material-dependent information provided in Appendix B is sufficient for the subsequent discussions. Table 1 Subsequently, the steps in the flowchart for model development in Table 1  ], i = U T, E T, P S are provided e.g. in [47].

Model development according to
Step 1: Verification for each mode data set Each set of experimental data for UT, ET and PS is used separately to minimize the mode data related least squares functionals f [d i , κ] in Eq. (54). We point out, that in our simulations no attempt has been made to modify respectively improve the optimized material parameters obtained in [47]. That is, for evaluation of the above criterion (54), for κ UT , κ ET , κ PS we exploit directly the results from [47] and also summarized in Eq.  Tables 2, 3 and 4, respectively.
Step 2: Stability for each mode data set In order to account for stability of all three vectors κ i , i = U T, E T, P S taken from [47], I-criteria ϕ I [κ], I = A, E, D, M according to Eq. (55) are calculated. A summary of related results is also given in Tables 2, 3 and 4, respectively.
Step 3: Validation for complete data set As in [47], each set of optimal material parameters κ UT , κ ET and κ PS which minimize the least squares functionals (54), is used to simulate the remaining two deformation modes. As mentioned before, for evaluation of the criterion (56) we useκ i = κ i , i = U T, E T, P S, according to Eq. (B,•,4). This step renders results for the complete data related least-squares functionals (56) in Tables 5, 6 and 7.
Step 4: Stability for complete data set In order to account for stability, in addition to [47] and based on the general structure in Eq. (50) In view of the extensive discussion in [47] on all thirteen models, in this work we restrict a detailed discussion to three representative hyperelastic models, that is the Neo-Hooke model (No. 1 in the list of Sect. 6.1), the Isihara model (No. 3 in the list of Sect. 6.1) and the Caroll model (No. 10 in the list of Sect. 6.1) to Appendix C of this work. The comparative study in Fig. 5 reveals a poor verification quality of the Neo-Hooke model. A good verification quality of the Isihara model with however disappointing validation quality is illustrated in Fig. 6. Figure 7 visualizes a superior model performance of the Carroll model on Treloar's data, due to both, a perfect fit quality and remarkable predictive results for UT, ET, PS, respectively.
Step 5: Selection of final parameter vector Here, according to Remark 5.4 we simply select a final parameter set κ * according to Eq. (58), in case Step 4 is successful. In view of the excellent model performance on verification and validation in Fig. 7 as well as on stability in Tables 5, 6 and 7, respectively, this applies only for the Caroll model. Consequently, the parameter vector κ UT in Eq. (B.10.4) could be selected for κ * , which renders lowest value for A [κ UT ] according to Table 5. The remaining models would require the weighted mode related model extension according to Step 6. In order to obtain an overview on the performance of all models, still, results of a fourth optimization step are shown next, that computes material parameters κ AL under simultaneous consideration of the complete data set set, that is, UT, ET and PS. To this end, as in the previous Sect. 6.3 no attempt has been made to modify and, respectively, improve the optimized material parameters κ AL obtained in [47] and summarized in Eq. (B,•,4) in Appendix B of this work. Clearly, as a consequence no additional data are left for validation.
A comparative study of the verification quality in this step for all original hyperelastic models on Treloar's data is given in Fig. 3 Table 8 for each model in the first line.
Contrary to the individual fit capability for single measurement sets, some of the models do not capture the S-shape behavior anymore. This deficiency becomes apparent, e.g., for the Isihara model, if the results in Fig. 6 are compared to those in Fig. 3.3. This illustrative effect is confirmed by the relatively high value for F[d, κ AL ] in Table 8. The deficiency of the material model of non-uniqueness discussed in Sect. C.2 for the Isihara model has dramatic consequences for almost all design functions I [κ] in Table 8.
As mentioned before, from the macroscopic models the Carroll model [9] captured very well all three deformation patterns with the lowest least-squares value F[d, κ AL ] of all models considered here. This positive phenomenon is accompanied by relatively low values for the design functions I [κ AL ].
From the microscopic models the eight-parameter model is able to reproduce almost perfectly all deformations. When compared with the Carroll model, higher values for the design functions I [κ AL ] confirm the conclusion of [47] of "a questionable sensitivity wrt the initial values", such that "small changes may already yield a completely different set of values."  Table 8, "simult. fitting"

Fig. 3 continued
Step 6: Weighted mode-related model extension All models, except the Carroll model, require a model extension according to Step 6 in Table 1 by means of the strain mode-dependent formulation (9.1). Comparative results are provided in Table 8, where additionally to the results of simultaneously fitting in Step 5 two cases are distinguished for each model: • The second line in Table 8 of each model corresponds to the weighted model according to Eq. (9). We use material parameters κ * = {κ i } 3 i=1 , i = U T, E T, P S taken from [47] and summarized in compact form in Eq. (B,•,4) in Appendix B of this work. The corresponding diagrams are provided in Fig. 4. Apart from the Neo-Hooke model, the Mooney-Rivlin model and the Gent-Thomas model, almost all models show a satisfying up to excellent fitting quality.
• The third line in Table 8 of each model corresponds to Eq. (9.1) with new with material parameters κ * = {κ i } 3 i=1 , i = U T, E T, P S obtained by model reduction, such that in generaln i p = dim |κ i | ≤ dim |κ i | = n i p holds. The resulting material parametersκ i are summarized in compact form in Eq. (B,•,5) in Appendix B of this work. Note, that model reduction has been performed only for some selected models. Table 8 reveals a drastic reduction of the corresponding I-criteria I [κ i ] in Eq. (57) in some cases. E.g. for the Isihara model the total number of material parameters n p is reduced from 9 to 6, with the same least squares functional F[κ i ] = 7.79e-1 in Eq. (56) but a reduction for the E-criterion E [κ] in Eq. (57) from ∞ to 7.97e-3.  Table 8, "weighted model" Model  Model Model

Summary and conclusions
This contribution presents a practical guide for model development accounting for the criteria of verification, validation and stability. As a prerequisite, it requires a set of different data sets, which can be related to different strain modes investigated experimentally in the laboratory. In case, the validity check on non-conformity according to Definition 1 in this work is not satisfying, an extension of hyperelastic material models is proposed by means of so-called strain mode-dependent weighting functions, as a counterpart to the extension in [26] with stress mode-dependent weighting functions. To this end, an additive decomposition of the strain energy function is assumed into a sum of weighted strain mode related quantities. This approach can easily be combined with techniques of model reduction, in order to obtain more stable material parameters.
In order to account for incompressibility and strain mode-dependent material behavior, the free energy density for the isochoric part has been postulated in a mixed formulation as ψ in terms of the eigenvalues¯ a and the principal isochoric stretches λ a , a = 1, 2, 3, and alternatively , κ i ] in terms of the eigenvalues¯ a and the three invariants I A , A = 1, 2, 3 of the isochoric right Cauchy-Green strain tensor. For convenient implementation, including the corresponding tangent operators, a general framework has been presented, which only requires the partial derivatives ∂ψ i /∂λ b , ∂ 2 ψ i /∂λ a ∂λ b in the first case and respectively the coefficients ∂ψ Table 8 Simultaneous fit for all hyperelastic models on Treloar's data for UT, ET and PS: Least squares functionals and I-criteria. For each model, the first line presents results with parameters κ * = κ AL , n * p = dim |κ AL | taken from [47], see also (B,•,4). The second line corresponds to the proposed extended formulation (9) including weighting functions. Parameters κ * = {κ i } 3 i=1 , i = U T, ET, P S are taken from [47], see also (B,•,4). Note the relation n * p = dim |{κ i } 3 i=1 | = 3 dim |κ AL | before model reduction. Any third line corresponds to model reduction of the scheme (9). The new parameters κ *   in the second case. The detailed results for these coefficients are provided for thirteen well-known material models in Appendix B. The approach is successfully applied for all 13 hyperelastic models with regard to the classical experimental data on vulcanized rubber published by Treloar [48], showing quadratic convergence for the finite-element equilibrium iteration as well as excellent fitting capabilties and stable material parameters.
We remark that in the field of optimal control problems, the indicator function φ I [C [κ]] is used to formulate an optimization criterion, which apart from material parameters may be dependent on further design variables, cf e.g. [4] and [24]. Then, in summary, two requirements become essential: The solution vector 1. Should minimize the least-squares functional (representing a model error) according to (37) 2. Should minimize the confidence criterion according to (50).
In practice, the above two requirements on small model error and stable results might be contradictory and, therefore, have to be carefully balanced. In this work, the generality of the methodology of optimal control problems has not be exploited and therefore will constitute an aspect of future work.

1.
∂ξ Consequently, from the relation (23)   Consequently, from the relations (A.2) the second derivative of the weighting functions becomes 1. UT: and respectively from (A.3)

UT:
The relations (A.5) and (A.6), respectively, serve to compute the coefficients W i ab in Eq. (28)

B Survey on hyperelastic models
In this section we summarize basic relations needed for the comparative study in Sect. 6 for the isochoric free energy density, the partial derivatives in (33) Fig. 6. In the following, we will discuss on the three aspects of 1. verification, 2. validation and 3. stability.
• Verification: As already discussed in [47], the fit quality is very high for all deformation modes. In particular for ET and PS almost perfect results are obtained. Regarding UT-data, the high initial stiffness is not properly reflected, however (in contrast to the Neo-Hooke model and the Mooney Rivlin model) Isihara's model can capture the characteristic S-shape at large stretches. Related results for f [κ i ], i = U T, E T, P S are given in Tables 2, 3 and 4. To summarize, the fit capability of Isihara model [20] for Treloa's data is very good. • Validation: Regarding the prediction capability, the results are rather disappointing for all three cases. In particular, the material parameters κ ET obtained from ET-data in the middle graph of Fig. 6 are not able to capture the S-shape curve of UT in the left graph. The same holds e.g. for the material parameters κ UT obtained from UT-data in the left graph, applied to ET-data in the middle graph. These illustrative results are confirmed by the comparative large values for the least-squares functionals F[κ i ], i = U T, E T, P S defined in Eq. (56) and summarized in Tables 5, 6 and 7. As a consequence, with reference to Definition 1 in the introduction of work we summarize: The conformity between the Treloar data [48] and the Isihara model [20] is very poor.  Tables 5, 6 and 7. To summarize: The stability of the Isihara model [20] is rather poor. A modification of experimental information, in whatever way, cannot account for this undesired property. A model reduction is indispensable. Fig. 6 Performance of the Isihara model on Treloar's data. The fit quality is acceptable (UT) to perfect (ET, PS) which renders high verification quality. Contrary, the validation quality is very poor, in particular in the middle graph where unrealistic high stresses are predicted for ET for parameters κ UT and κ PS , respectively Carroll [9] proposed a phenomenological approach based on a successive extension of the free energy. Here the additional terms are chosen such that they reduce the errors that remain in the stress response of the previous terms, compared to measurements, see e.g. [9] on more details. The free energy function related to the i-th strain mode is given in Eq. (B.10.1). Each of the first three sets of material parameters in Eq. (B.10.4) is used to predict the two complementary deformation modes, and results are plotted in Fig. 7. In the following, we will discuss on the three aspects of 1. verification, 2. validation and 3. stability.
• Verification: As already discussed in [47], the fit quality of the parameter sets are very good to perfect.
Related results for f [κ i ], i = U T, E T, P S are given in Tables 2, 3 and 4. • Validation: The corresponding curves on prediction in Fig. 7 and reveal a remarkable model performance.
These illustrative results are confirmed by the comparative low values for the least-squares functionals     Tables 2, 3 and 4, as well as of I [κ i ], I = A, E, D, M of Eq. (57) and listed in Tables 5, 6 and 7 reveals remarkable stability properties.

D The correlation matrix of Isiharas model
This part of the Appendix provides results for the correlation matrix C i j [κ * ] defined in Eq. (40). For briefness of presentation, only results for the Isihara model after model reduction are shown. The comparatively high correlations in Table 10 between some of the material parameters (e.g. within the UT regime -0.9172 between c U T 10 and c U T 20 ) is an indicator for overparameterization, but, no attempt has been made on further reduction. However, the results in Table 10 also demonstrate a clear independence between UT, ET and PS related material parameters. This is regarded as a main advantage of the proposed weighted strain energy in Eq. (9).

E Local check for quadratic convergence
This section intends to verify the results for the isochoric tangent operator in Eq. (28), consistent with the isochoric second Piola-Kirchhoff stress tensor in Eq. (26). Thereby, we are only interested in the convergence performance, such that the checks are performed as one-element tests (rather than for a complex finite element geometry). are used. Solutions of the resulting equilibrium problems (not discussed in more detail here) are performed with a Newton method, which iterates on the corresponding finite-element residual of the element. Table 11 summarizes the final residual norms for all thirteen models in the list of Sect. 6.1. In almost all cases the number of correct digits is doubled in every iteration, that is quadratic convergence is achieved, which verifies correct derivation and implementation of the corresponding tangent operators.