Systematic Fitting and Comparison of Hyperelastic Continuum Models for Elastomers

Hyperelasticity is a common modeling approach to reproduce the nonlinear mechanical behavior of rubber materials at finite deformations. It is not only employed for stand-alone, purely elastic models but also within more sophisticated frameworks like viscoelasticity or Mullins-type softening. The choice of an appropriate strain energy function and identification of its parameters is of particular importance for reliable simulations of rubber products. The present manuscript provides an overview of suitable hyperelastic models to reproduce the isochoric as well as volumetric behavior of nine widely used rubber compounds. This necessitates firstly a discussion on the careful preparation of the experimental data. More specific, procedures are proposed to properly treat the preload in tensile and compression tests as well as to proof the consistency of experimental data from multiple experiments. Moreover, feasible formulations of the cost function for the parameter identification in terms of the stress measure, error type as well as order of the residual norm are studied and their effect on the fitting results is illustrated. After these preliminaries, invariant-based strain energy functions with decoupled dependencies on all three principal invariants are employed to identify promising models for each compound. Especially, appropriate parameter constraints are discussed and the role of the second invariant is analyzed. Thus, this contribution may serve as a guideline for the process of experimental characterization, data processing, model selection and parameter identification for existing as well as new materials.


Introduction
Elasticity describes an idealized mechanical material behavior whose stress response depends only on the current deformation state. Once the external load is removed, the material is expected to recover its initial, undeformed configuration and to release the entire work done by the external load, i.e., no hysteretic effects occur. The concept of hyperelasticity ensures such a behavior by defining the material response in terms of a potential (commonly called strain energy density function W) from which the constitutive equations are derived. Thus, a path independent mechanical work is guaranteed and the relation between the stress and deformation tensor is provided by a scalar function that is conveniently to handle. Furthermore, the laws of thermodynamics can readily be evaluated and a variational functional for the balance of the linear momentum can be formulated. As a consequence of these desirable properties, hyperelasticity is a common modeling approach in solid mechanics. Especially materials that can undergo large deformations, like elastomers or soft tissues, are typically described by a strain energy function.
In the literature, a vast number of strain energy functions have been published, particularly for rubber materials at finite strains. However, their validity is rarely tested for different rubber types or compounds such that an extensive database for a well-founded model selection is missing. Moreover, most publications focus on modeling either the isochoric response or, less often, the volumetric behavior. Also review articles which compile lists and benchmark tests of strain energy functions for rubbers are limited to a few compounds and often stick to the classical Treloar [95] or Kawabata data [56] stemming from sulfur crosslinked natural and isoprene rubber, respectively. An overview of these 1 3 articles is given in Tables 1 and 2 including the number of considered models, the type of rubber used for fitting as well as the objective of the work.
The present manuscript provides a systematic fitting of both, isochoric and volumetric hyperelastic models. They are fitted to the experimental data of nine different rubber compounds and, thus, an extensive database is established. In addition and in contrast to the existing literature, the dependence of the strain energy functions on the first and second principal invariant is investigated. Particularly, the role and importance of the second principal invariant for rubber models is studied in detail. Also, the effects of possible cost function formulations on the fitting results are demonstrated, Table 1 Review articles on isochoric hyperelastic models for rubber materials (number of models n includes only models for rubbers; CB: carbon black, HBNR: hydrogenated nitrile butadiene rubber, NR: natural rubber, S: sulfur, Si: silicone) Literature n Experimental data Objective [37] 10 None Deriving a new model [34] 8 4 S-crosslinked NR from literature Visual comparison, giving attention to small strain behavior, considering applicability of Valanis-Landel hypothesis [103] 5 S-crosslinked NR Discussing the models' capability and drawbacks, considering also swollen rubber, proposing a model extension to improve the fitting at small strain [13] 8 Treloar data Visual comparison, discussing predictivity to other deformation modes, addressing the effects of compressibility [88] 6 lowly CB-filled HNBR Visual comparison, considering conditioned and unconditioned state, discussing predictivity to other deformation modes, studying of finite element simulations [96] 45 none Providing list of references [46] 5 none Discussing behavior at large strains, presenting extensions to thermoelasticity, anisotropy and compressibility [70] 7 Si-rubber Visual comparison and in terms of the correlation coefficient, considering also pig muscular tissue [68] 20 Treloar & Kawabata data Compiling rankings [47] 25 Treloar data & NR & Si-rubber from [71] Grading of the model performance, deriving a new model [91] & [48] 14 & 11 Treloar data Providing and testing consistent linearizations, discussing performance and predictivity to other deformation modes [11] 25 Treloar data & S-crosslinked NR from [104] Presenting building-strategy for new models [49] 8 Treloar data Considering micromechanically motivated chain and network models, visual comparison, studying predictivity to other deformation modes [33] 8 Treloar data Developing an application for parameter identification and verification [74] 12 Treloar data Modeling of nonlinear elasticity using deformation dependent parameters [22] 44 Treloar & Kawabata data Compiling rankings, discussing alternative cost functions with weights for each experiment, presenting a generic fitting algorithm, discussing predictivity to other deformation modes [30] 7 2 CB-filled, S-crosslinked NRs Proposing a new model, analyzing temperature dependencies [41] 75 Treloar data & CB-filled HNBR from [45] Compiling top ten rankings, comparing run time of finite element simulations Table 2 Review articles on volumetric hyperelastic models (n is the number of models; CB: carbon black, EPDM: ethylene propylene diene rubber, S: sulfur) Literature n Experimental data Objective [66] 5 None Discussing constraints on volumetric strain energy functions [26] 5 None Discussing constraints on volumetric strain energy functions, proposing new models [38] 10 None Discussing polyconvexity (also for isochoric models), proposing new models [83] 4 CB-filled, S-crosslinked EPDM Illustrating the effect of CB content, crosslink content and curing time on the bulk modulus reasonable constraints on the parameter bounds are recapped and crucial model properties are highlighted. Another difference between the present work and existing contributions is that the treatment of the experimental data is highlighted. Since an accurate parameter identification requires a careful preparation of the raw test data, some common issues regarding data processing are addressed. On the one hand, a physically based approach is outlined to consider the preload in biaxial tensile tests as well as compression tests. On the other hand, a straightforward method is proposed to check the consistency between data stemming from multiple experiments.
The manuscript is organized as follows. First, the theoretical foundations are laid in Sect. 2 including continuum mechanics, thermodynamics, the experimental setup and parameter identification. Next in Sect. 3, the treatment of the experimental data is explained. Finally in Sects. 4 and 5, the fitting procedure for the isochoric and volumetric strain energy functions is outlined, conducted and discussed.

Kinematics
To describe the deformation of a solid body, the motion of a material point is tracked. Denoting its position in a reference configuration (typically the undeformed state) and the current configuration by X and x(X, t) , respectively, the displacement u = x − X is introduced. Then, the deformation gradient is computed which is an invertible, two-field tensor. I denotes the identity tensor. Rotational and distortional contributions to the deformation gradient are separated by the polar decomposition with the proper orthogonal rotation R and the symmetric right stretch tensor U . Moreover, U can be split into a volume-preserving deformation and a pure dilatation represented by the isochoric (i.e., unimodular) right stretch tensor Ū = J −1∕3 U and the volumetric contribution J 1∕3 I , respectively. Here, J = det(U) denotes the determinant of U and measures the local volume change. Note that det(F) = det(U) holds true due to the multiplicativity of the determinant and det(R) = 1 . The multiplicative decomposition of the deformation gradient yields The square of U is called right Cauchy-Green tensor whose principal invariants are with the trace tr( ⋅ ) , adjugate adj( ⋅ ) = det( ⋅ )( ⋅ ) −1 and Frobenius norm ‖ ⋅ ‖ . Accordingly, the square of Ū will be referred to as isochoric right Cauchy-Green tensor and the isochoric principal invariants are defined as The eigenvalues of U are called principal stretches and denoted by 1 , 2 , 3 whereas the eigenvalues of Ū will be referred to as isochoric stretches ̄k = J −1∕3 k , k = 1, 2, 3 . Both stretch tensors (and their squares) share the same eigenvectors G k , viz.,

Stresses
The Cauchy stress (also known as true stress) is defined as force per unit current area. More precisely, the stress tensor relates the normal n of an imaginary cut surface da in the current configuration and the resultant force df acting on the cutting surface such that In contrast, the 1st Piola-Kirchhoff stress P (also known as engineering stress) describes the force per unit reference area, viz., where N and dA are the corresponding normal vector and cutting surface in the reference configuration. The 1st Piola-Kirchhoff stress is a two-field tensor, work-conjugated to the material time derivative of the deformation gradient Ḟ and linked to the Cauchy stress by with the transposed inverse ( ⋅ ) −T . The third stress measure used in this manuscript is the 2nd Piola-Kirchhoff stress S which is work-conjugated to 1∕2Ċ . It is based in the reference configuration, lacks a physical meaning and can be computed by The Cauchy stress can be additively decomposed into an isochoric (i.e., traceless) and a volumetric contribution iso and vol , respectively, commonly called deviatoric and hydrostatic stress. They are given by where p is called hydrostatic pressure. Applying Eq. (12) to the deviatoric and hydrostatic Cauchy stress leads to The stress power per unit reference volume is given by In analogous manner to the stress, it splits additively into an isochoric and a volumetric part, viz., where the identity ̇J = J∕2 C −1 ∶Ċ is used.

Thermodynamics
The second law of thermodynamics for isothermal conditions can be represented via the Clausius-Planck inequality per unit reference volume with the mechanical dissipation rate D m and the free energy per unit reference volume Ψ . For hyperelasticity, the free energy is given by a strain energy density function W. Assuming that W is defined in terms of the right Cauchy-Green tensor, the chain rule Ẇ = W∕ C ∶Ċ and Eq. (15) can be applied resulting in Thus, the Clausius-Planck inequality is always fulfilled by defining the stress-deformation relation (16) P iso = P − P vol and P vol = −ṗJ cf. Eq. (12), implying a material behavior free of dissipation, i.e., D m = 0.

Strain Energy Density Functions and Imposed Restrictions
The strain energy functions discussed in this manuscript are assumed to be composed of an isochoric part W iso in terms of the isochoric invariants Ī 1 , Ī 2 and a volumetric part W vol depending on the volume change J, viz., Such a separate modeling of the deviatoric and the hydrostatic response is a reasonable and widely used approach for nearly incompressible materials like elastomers 1 . Employing Eq. (19) 1 with the isochoric-volumetric split of the strain energy according to Eq. (20) leads automatically to a decoupled stress response, cf. Eq. (14). With the derivatives An alternative form of a decoupled strain energy is obtained by defining the isochoric part of Eq. (20) in terms of the isochoric stretches ̄k . To satisfy a priori symmetry requirements due to isotropy, the Valanis-Landel assumption [97] is employed leading to the form ̄k such that all three stretches are treated uniformly and separately by a scalar function , cf. [53,98,55] for the validity and limitations of this assumption regarding experimental findings. The corresponding 2nd Piola-Kirchhoff stress is obtained as see for instance [80]. The forms W iso Ī 1 ,Ī 2 and W iso ̄1 ,̄2,̄3 will be referred to as invariant and stretch formulation, respectively.
The second law of thermodynamics, cf. Sect. 2.3, does not impose any restrictions on the construction of the strain energy function. Therefore, many authors discussed requirements for an objective, physically plausible and numerically desirable material behavior. Objectivity is a very fundamental demand for material frame indifference, see for instance [39]. It is automatically satisfied by an appropriate choice of arguments of the strain energy function, e.g., the invariants in Eq. (20) or the eigenvalues in Eq. (23). Physically motivated restrictions stem from experimental or empirical observations, e.g., • Greater stress occurs always in the direction of the greater strain, Baker-Ericksen inequalities [5] • Incremental mechanical work done by external loads must be positive, Drucker-Hill postulate [28,44] • Tangential shear and bulk modulus must be positive • Sound speed in a material must be a non-complex number These wordings must then be turned into verifiable mathematical formulations within the finite strain theory and typically lead to constraints on the set of feasible parameters. Also constraints when approaching certain limits or deformation states, e.g.,  [6,7,69,77] for an overview and implications between these conditions. Since some of the numerical restrictions have a direct physical consequence, the boundary between both categories is blurred. Furthermore, a few restrictions do not even affect the constitutive equations but are reasonable for convenience, e.g., normalization conditions like W(F = I) = 0.
The choice of a suitable restriction depends on the application and is often a trade-off between numerical stability and the implied limitations. For instance, [55,9] showed experimentally that W∕Ī 2 < 0 for small biaxial deformations of natural rubber. Too strong restrictions can easily rule out the capability of a model to reproduce this behavior. Moreover, the demand for polyconvexity would preclude a priori several W iso,2 -functions from the study in this manuscript as [38] showed that Ī m 2 is not a polyconvex term for m < 3∕2 . In addition, strict convexity of W(F) is violated in the case of buckling or phase transitions. Sects. 4.1.1 and 5 discuss the restrictions made for each function W iso,1 Ī 1 , W iso,2 Ī 2 and W vol (J).

Incompressibility
Nearly incompressible materials like rubbers and fluids are characterized by a stress response being much larger for volumetric deformations than for isochoric deformations. Thus, the material volume tends to remain nearly constant, viz., J ≈ 1 and their Poisson ratio is close to 0.5. For analytical considerations, nearly incompressible materials are often assumed to be perfectly incompressible, i.e., J = 1 always holds true and = 0.5 . In this case, the stretch tensor U and the isochoric stretch tensor Ū as well as their invariants and eigenvalues coincide. Furthermore, the hydrostatic pressure p in Eq. (14) is not given by a constitutive equation as in Eq. (21) but is determined by the static equilibrium.

Deformation Modes for Material Characterization
The deformation and stress states of five frequently used experiments for the characterization of the mechanical material behavior of rubbers are considered in this work: uniaxial tension, pure shear (also known as planar tension), equibiaxial tension, simple shear and confined compression (also known as volumetric compression). These idealized deformation modes are illustrated in App. B. The uniaxial tension, pure shear and equibiaxial tension state (denoted by ux, ps and bx) can be derived as special cases of the general biaxial deformation mode where 1 = denotes the stretch in load direction x. The z-direction is considered as stress-free and the corresponding stretch 3 = 1∕( 1 2 ) follows from the assumption of incompressibility. See Table 3 for the lateral stretch 2 in and R = I y-direction, the experimental stress obtained from the measured force F x and the hyperelastic model response in each deformation mode.
Next, the simple shear test (ss) is introduced with the shear strain s, the deformation gradient and J = 1.
Furthermore, the deformation gradient of the confined compression test (cc) reads In case of nearly incompressible materials, xx ≈ yy = zz can be observed, i.e., the deviatoric stress contribution iso is negligible such that the approximation ≈ vol = −pI is valid. Thus, contrary to the biaxial tensile and simple shear tests which provide information about the isochoric response, the confined compression test reveals the volumetric behavior.

Experiments
Experimental data of nine unfilled rubbers from quasistatic, uniaxial tensile tests as well as confined compression tests are employed, see Table 4, Fig. 1 and [84]. The uniaxial test data are used for the parameter identification of (27) Pure shear (ps)/ incompressible planar tension the Ī 1 -dependent strain energy functions W iso,1 in Sect. 4.2, whereas the volumetric strain energy functions W vol are assessed with the data from confined compression tests in Sect. 5. All experiments were conducted with a Zwick 1445 universal test machine, cf. Supplementary Fig. 24(a). The tensile test procedure is based on the DIN 53504 standard (preload: 0.1 N ; displacement rate: 200 mm∕min ; flat, shouldered specimen with 20 mm gauge length; optical strain measurement; loading until failure). For the confined compression tests, a tight-fitting, cylindrical specimen (diameter 7.5 mm ; height 10 mm ) is placed in the hole of a thickwalled metal cylinder. Thus, the lateral strain is suppressed when a compressive force is applied and a volume change is enforced. The data of both experiments are processed according to the procedures in Sect. 3.1 and 3.3, respectively, and afterwards resampled with equidistant strain increments of 0.5 % (i.e., Δ = 0.005 ) to ensure an equal weight of all strain regions in the fitting procedure. Furthermore, the classical Treloar data set [95] comprised of uniaxial tension, pure shear and equibiaxial tension data of unfilled, sulfur-crosslinked natural rubber, cf. Fig. 2, is used for the discussions on Ī 2 -dependent strain energy functions W iso,2 in Sect. 4.3. It is processed according to Sect. 3.1 and its consistency is checked in Sect. 3.2. The choice of the Treloar data is explained by means of Supplementary  Table 11 which lists frequently used experimental data sets for benchmarking constitutive models for rubbers (filled, thermoplastic and foamed elastomers are excluded since they cannot generally be idealized as volume preserving or free of dissipation and, hence, conflict with the model assumptions of perfect elasticity and incompressibility). First, it can be noticed that a lot of publications in the table provide general biaxial tension data sets rather than separate  Table 4 (data were processed according to procedures in Sect. 3.1 and 3.3, respectively; legend in (b) also applies to (a)) uniaxial, planar and equibiaxial tension data, cf. [52,53,56,85,98]. Such data require further discussions which are beyond the scope of this work, e.g., how to properly treat the second stress value P yy in the cost function, see Sect. 2.8, or how to apply the data preparation and consistency check from Sect. 3. Moreover, the listed general biaxial tensile tests do not always cover the full range of biaxial deformation from uniaxial ( 2 = 1∕ √ 1 ) to equibiaxial tension ( 2 = 1 ) at all 1 -levels, cf. [53,85,98], or the maximum stretch is comparatively small. The remaining references often lack in detailed information about the experimental procedure or present incomplete data sets, i.e., the uniaxial, planar or equibiaxial tensile test is missing, cf. [1,3,104]. In addition, as nearly all data sets in the table are based on similar natural rubber compounds with comparable properties, the evaluation of several of these sets would lead to redundant findings. Due to the outlined considerations, the probably most frequently used and well established Treloar data are considered here, see for example [11,13,16,36,57,58,63,67,73] or Table 1. Although his measurement methods are outdated and he did not table the test data, Treloar [95] explained in detail his experimental setups and findings. Moreover, he presented interesting discussions on the systematic error in pure shear experiments, the equivalence of equibiaxial tension and uniaxial compression and the need of sample preconditioning by comparing the loading and unloading curves. Nevertheless, extensive, complete and consistent data sets for several rubber compounds from a modern test machine and up to large stretches with an analysis of error sources are still missing in the literature.

Parameter Fitting
Parameter fitting (also known as model calibration or parameter identification) is an optimization procedure which searches for the best parameter set of a material model such that the difference between the model response and given experimental data is minimal in a certain sense. Here, the p-norm of the residual array is minimized leading to the cost function with the residual of the i-th data point (of experiment k with m k data points) where j , j = 1, … , n are the material parameters. exp,i and P exp,xx,i are the measured stretch and 1st Piola-Kirchhoff stress in direction of tension/compression and s exp,i and P exp,xy,i are the measured shear strain and stress. The model response P mod,xx,i or P mod,xy,i is computed according to Table 3. Common choices for the norm order are p = 1 (taxicab norm), p = 2 (Euclidean norm) or p = ∞ (maximum norm), i.e., the mean error, the root mean square error (RMSE) or the maximum error are minimized. When multiple experiments are fitted simultaneously, the cost function is given by F = ∑ k F k . In this case, the factor 1∕m k in Eq. (28) eliminates a weighting of the experiments due to different numbers of data points.
The residual in Eq. (29) is defined in terms of the absolute error of the 1st Piola-Kirchhoff stress. Other reasonable error types are the relative error or the normalized error where the indices xx or xy are omitted for the sake of a more general notation. These two residuals scale the error from each experiment in order to eliminate a weighting when multiple experiments with different stress ranges are fitted, see for example the Treloar data in Fig. 2 with a maximum uniaxial stress of 6.6 MPa and a maximum equibiaxial stress of 2.5 MPa . In addition, the relative error in residual (30) gives a higher weight to the range of low stresses than the 2 Treloar data [95] with close-up of the small strain regime (data were processed according to procedure in Sect. 3.1; ux: uniaxial tension, ps: pure shear, bx: equibiaxial tension) absolute and normalized error but may generate dubious r i when dividing by noisy stress values P exp,i close to zero. An alternative formulation of the cost function is obtained by replacing the 1st Piola-Kirchhoff stress P by the Cauchy stress . For the biaxial tensile tests (ux, ps and bx), this is equivalent to multiplying the residual by exp,i > 1 , cf. Table 3. Thus, it amplifies deviations from large stresses. Vice versa, replacing P by the 2nd Piola-Kirchhoff stress S leads to a division by exp,i for these deformation modes and, hence, errors at small stress are more penalized. For simple shear, the 1st Piola-Kirchhoff stress is identical to the Cauchy stress P xy = xy . Whereas the 2nd Piola-Kirchhoff stress obtained from S xy = P xy − s P yy is experimentally not accessible since the stress component P yy is usually not measured 2 . For confined compression, the 1st Piola-Kirchhoff stress is equal to the Cauchy stress, too, whereas the 2nd Piola-Kirchhoff stress scales P by 1∕ exp,i . In Sect. 4.1.3, the effects of the norm order, stress measure and error type on the fitting result are illustrated so that a reasonable choice for the subsequent studies and model comparisons can be made.
For the parameter identification, the SciPy optimize library is used. The implementation of the hyperelastic models is done with the Mathematica Add-On AceGen, cf. [62]. More precisely, the strain energy function is programmed in the symbolic Wolfram language such that the stresses as well as the exact Jacobian, i.e., derivative of the stress w.r.t. the material parameters are obtained via automatic differentiation. Then, AceGen generates numerically efficient Fortran or C subroutines which are callable within the SciPy algorithms, e.g., via the F2PY tool of the NumPy library.

Preload Correction in Tensile Tests
For tensile tests with flat specimens, a small preload is typically applied to ensure a plane, non-buckling reference configuration. Moreover, it avoids difficulties for the controller of the test machine when encountering such a possibly unstable load state. The preload is ideally chosen as small as possible while it concurrently guarantees a straightened specimen. This balance may depend on the material's stiffness, the test machine settings and the specimen clamping, e.g., for different loading modes. Thus, no universal, best choice exists. Possible preload correction procedures are discussed below.
Since the preload leads to a slight stress offset of the first data point, many test machines provide an option to zero the load cell after applying the preload. As a consequence, the whole stress-strain curve is shifted along the stress-axis such that the first data point provides zero stress at zero strain. This procedure (hereinafter referred to as initial stress correction or vertical shifting) is not ideal since the preload is actually present and should not be ignored. To overcome this drawback, an alternative to the stress correction is provided. Let us assume that the testing machine measures optically the distance in loading direction between two reference points on the specimen rather than the displacement of the traverse. The measured length in the current state at time t i and in the preloaded reference state are denoted by i and * 0 , respectively, leading to the erroneous stretch * i = i ∕ * 0 . Hence, the initial length * 0 must be corrected because the elongation due to the preload, i.e., the prestretch pre = * 0 ∕ 0 with the true initial length 0 is not captured by the testing machine. This correction does not lead to a shifting of the stress-axis but to a scaling of the stretch-axis by pre , viz., the true stretch is obtained by Therefore, this procedure is referred to as initial stretch correction or horizontal scaling. In the following, a procedure to find a good approximation of pre is proposed: 1. Extract the data at small strains. The exact range must be checked for each material and deformation mode. For the materials depicted in Fig. 1(a), 1 < < 1.1 is a good overall choice. 2. Choose an appropriate strain energy function which is able to reproduce the small strain behavior well. For the unfilled rubbers in Fig. 1(a), a Neo-Hooke approach W iso = c 10 Ī 1 − 3 is sufficient 3 . 3. Fit the chosen strain energy function to the extracted data with a modified stretch input according to Eq. (32). That is, i = * i pre where * i are the measured stretches and pre is the prestrech which is treated as an additional fitting parameter. If data sets from multiple deformation modes are available, this fitting should be done for each experiment separately. 4. Scale all measured stretches * i by the fitted prestretch pre to obtain the corrected stretches: i = * i pre .
2 Note that a simple shear deformation state F = I + s e x ⊗ e y is assumed here which generally does not lead to a simple shear stress state = xy e x ⊗ e y + e y ⊗ e x , i.e., xx = yy = zz = 0 for finite strains, see for instance [93].

3
The proposed procedure is applied to all compounds in Fig. 1(a), see Table 6 for the fitted shear moduli G 0 = 2 c 10 . It is illustrated for the CR+S compound in Fig. 3(a). The fitted value pre = 1.051 states that the specimen of this compound was prestretched by 5.1 % due to the preload. The effect of the initial stretch correction is compared to the initial stress correction in Fig. 3(b). The difference between both approaches is pronounced particularly at large strains. Furthermore, Fig. 4 depicts the effect of the initial stretch correction on the stress-stretch curves of uniaxial tensile tests with the same material but different preloads from −0.25 N (i.e., a buckled specimen with an underestimated initial length) to 4 N (i.e., a strongly prestretched specimen with an overestimated initial length). It can be observed that the corrected curves nearly coincide over a large range of preloads. Only the stress responses with the negative and highest preload deviate slightly at large strains. In contrast, the stress correction is not able to reasonably compensate the preload dependency as it would lead to a shift of all data points by only −0.03 … 0.5 MPa ( −0.25 … 4 N on a cross section area of approx. 2 × 4 mm 2 ). The effect would be barely visible in the stress-stretch plot. Thus, the presented initial stretch correction should be preferred. However, the choice of the applied preload is still to be made carefully and depends on several factors as outlined at the beginning of this section. In general, negative and very large preloads should be avoided.

Consistency Check of Tensile and Shear Tests
Experimental data obtained from uniaxial tensile tests are very reliable and reproducible. In contrast, pure shear, equibiaxial tension and simple shear measurements require more difficult, error-prone setups which do not yield exactly the desired deformation states. For instance, the pure shear state The additional fitting parameter pre is used to scale the stretch data (orange); For illustration, the conventional initial stress correction which shifts the data vertically, viz., zeroing the load cell after applying the preload is depicted as well (green). (b) Effect of the initial stretch as well as initial stress correction on the whole data range (the fitted parameters are pre = 1.051 , c 10 = 0.225) Fig. 4 Proposed initial stretch correction applied to tensile tests with the same material but different preloads: The raw data are plotted as dashed lines and, for the sake of readability, with an offset of 10 MPa , cf. right axis, whereas the corrected data are given as solid lines is commonly obtained from a tensile test using a specimen with a high aspect ratio (width to length), cf. Supplementary  Fig. 24(b). The high aspect ratio minimizes the lateral contraction but, nevertheless, leads to a non-homogeneous deformation at the edges. Moreover, this setup requires highly parallel and concentric lower and upper clamping as well as specimens with a constant thickness over the specimen width. In case of equibiaxial tensile tests with a clamping frame as shown in Supplementary Fig. 24(c), the friction between the clamps and the rails causes an overstiff response. Note also that the simple shear deformation given in Eq. (26) is just an idealization of the real distortion, especially for large strains.
Due to the variety of errors, experimental data can readily be inconsistent, i.e., the measured response from different deformation modes do not comply with theoretical considerations. For invariant-or stretch-based hyperelastic models under the assumption of incompressibility, the ratios between the initial tangent moduli are fixed and independent of the chosen strain energy function: If the experimental data do not exhibit these ratios, no model would be able to reproduce the behavior properly. Therefore, it is advisable to check the consistency before running the model calibration: 1. Run the initial stretch correction from Sect. 3.1 for each experiment. If an initial stretch correction is not desired or needed, run the procedure anyway but keep the prestretch constant pre = 1 , i.e., the strain energy function is fitted to the unmodified data.
Use the identified parameters of the strain energy function to compute the initial moduli according to Eq. (33).
For illustration, this procedure is applied to the Treloar data in Fig. 2 considering the Neo-Hooke model and the strain range 1 < < 1.35 . The fitted initial moduli and prestretches are given in Table 5. Both, the initial modulus of the pure shear data as well as the equibiaxial data are slightly too stiff but still in good agreement with the uniaxial Young's modulus. Thus, reasonable results can be expected when fitting hyperelastic models under the assumption of incompressibility.

Data Processing of Confined Compression Tests
Confined compression tests typically show a too compliant behavior at small deformations because the sample diameter is slightly less than the bore diameter of the specimen holder. When the pressure is applied, the material can expand a bit in lateral direction such that initially the response is predominantly isochoric. Once the nearly incompressible rubber comes into contact with the inner surface of the specimen holder, the stress response increases notably due to the much stiffer volumetric behavior. Unfortunately, this stress upturn is smooth rather than abrupt. Hence, the challenge is to find a point that represents this transition in order to approximate the specimen height h 0 when the volumetric deformation begins. A suitable approach is proposed below and exemplified in Fig. 5(a) using the IR+S compound. The crucial differences to the preload correction in Sect. 3.1 are: it is presumed that the applied preload is negligibly small, the height of the undeformed specimen h * 0 is known and the displacement is measured rather than the strain.
1. Choose a lower pressure bound (e.g., 15 MPa in Fig. 5(a)) above which the response can be definitely assumed to be volumetric and truncate the data below.   Table 6 Initial shear moduli G 0 and bulk moduli K 0 according to the procedures in Sect. 3.1 and 3.3, respectively, as well as Poisson ratios obtained by 0 = 3 K 0 ∕G 0 − 2 ∕ 6 K 0 ∕G 0 + 2 and measured mass densities 0 of the rubbers in Fig. 1  where u * i and u i are the measured and corrected displacement at time t i . Note that the displacement is defined to be positive in compression direction. u 0 is the applied displacement, i.e., the height reduction when the specimen comes into contact with the inner surface of the specimen holder and the volumetric deformation begins. It is considered as an additional fitting parameter besides K 0 . 4. Use the identified u 0 to derive the volume change from the measured displacement: The preload correction is applied to all confined compression tests from Fig. 1(b), see Table 6 for the fitted bulk moduli K 0 . Moreover, the K 0 values are plotted against the shear moduli G 0 from Sect. 3.1 in Fig. 5(b). No correlation between the bulk modulus, shear modulus, Poisson ratio or mass density can be observed. However, the provided initial shear and bulk moduli may serve as a data source for future numerical work lacking in experimental data and, hence, realistic parameter values.

Fitting Procedure, Considered Models and Parameter Limits
The fitting procedure for the isochoric strain energy is divided into two parts. At first, the uniaxial tension data of the nine rubbers depicted in Fig. 1(a) are used to identify promising Ī 1 -dependent strain energy functions for each compound. Thereafter, the best performing Ī 1 -based models are combined with several Ī 2 -terms and fitted simultaneously to the uniaxial tension, pure shear and equibiaxial tension data of the natural rubber from Treloar [95]. This separated procedure is followed to systematically investigate the role of both terms.   19 Boyce & Arruda [13] Eq. (27)  33 Fu [30] Eqs. (10)-(11) The considered strain energy functions of this study are collected from the literature and compiled in Tables 7, 8 and 9. Note that for models depending on both invariants Ī 1 and Ī 2 one of them is formally set equal to 3 to obtain separated formulations (in a similar manner, possible stretch-based terms are eliminated by setting all ̄k = 1 ). That is, some models differ from their original formulation (indicated by an asterisk * in the tables) and approaches with coupled terms are beyond the scope of this work. All strain energy functions are proven to satisfy the following requirements: • For Ī 1 -and Ī 2 -based strain energy functions, the maximum number of parameters is five and two, respectively. • A closed-form representation of the strain energy in terms of elementary functions must be available. In addition, models requiring numerical integration, internal iterations or other computationally costly aspects are excluded as well. • Only models that are particularly developed for rubbers are considered. For instance, some publications consider also models for isotropic soft tissues. However, due to the characteristic tension-compression asymmetry of tissues (i.e., the stress response is more pronounced under compression than under tension), different modeling approaches are necessary, see also [75,15].
All remaining material parameters can assume any real number. Indeed, these requirements are very generous and may lead to numerically undesirable or physically implausible behavior for arbitrary deformation states. However, finding suitable parameter bounds for each strain energy function is not reasonable for a large number of models. It depends on the imposed restrictions, cf. Sect. 2.4, and often leads to an optimization problem subjected to additional constraints. Moreover, promising models can easily be overlooked due to too strong constraints. Especially, when combining Ī 1 -and Ī 2 -models, parameter limits are often too restrictive, e.g., demanding convexity for both W iso,1 and W iso,2 . Therefore, the results are proven after the fitting to fulfill monotonicity conditions, viz., the Cauchy stress as well as the 1st Piola-Kirchhoff stress increase monotonically in the experimental strain range. A stability check beyond the maximum strain of the considered data is not taken into account since the specimens are loaded until failure. Any discussion on a physical plausible behavior beyond the elongation at break is pointless and, for numerical stability, an extrapolation with a constant modulus can be implemented.  Demanding monotonicity of the Cauchy stress is a fundamental, physically motivated requirement and equivalent to satisfying the Baker-Ericksen inequalities Whereas monotonicity of the 1st Piola-Kirchhoff stress is a necessary condition to ensure a stable material behavior in the sense of Hill and Drucker [44,28]. That is, the incremental work defined by ΔF ⋅ Δu must be positive where Δu is the displacement increment due to the additional external force ΔF . For the biaxial deformation states, it is equivalent to ΔP xx Δ > 0 and for simple shear to ΔP xy Δs > 0 . Note that monotonicity of the 1st Piola-Kirchhoff stress implies monotonicity of the Cauchy stress.

Role of the Second Principal Invariant
Many hyperelastic material models are defined in terms of the first principal invariant Ī 1 or in terms of the first and second principal invariant Ī 1 and Ī 2 , see for instance [20] for an overview. The additional dependency on Ī 2 , its physical meaning and its role within the fitting procedure are discussed by many authors. A comprehensive review on experimental findings was given by Kawabata & Kawai [55]. Assuming incompressibility, they discussed the stress dependency (more precisely, the derivatives W∕Ī 1 and W∕Ī 2 ) for biaxial deformation modes within the small as well as large strain regime and at different temperatures. Also the relaxation behavior, the effect of filler and different polymers were investigated. In context of construction and calibration of strain energy functions, two important outcomes are: W∕Ī 1 > W∕Ī 2 holds true for all deformations and W∕Ī 2 can take negative values at small strains. Under the assumption of a simple network, Kawabata & Kawai [55] derived a micromechanical interpretation of their observations stating that the Ī 1 -dependency is ,,related primarily to intramolecular forces" while the Ī 2 -dependency is "a manifestation of intermolecular interactions". However, it can be noticed that physically motivated modeling approaches rarely result in a strain energy function in terms of the second principal invariant. For instance, so-called tube models are typically comprised of a crosslink part defined in terms of Ī 1 and a stretch-based part accounting additionally for chain entanglements, cf. [29,42].
Destrade et al. [24] discussed the Ī 2 -dependency based on observations from systematic fittings. They stated that "the Ī 2 -dependence is precisely the missing ingredient to obtain excellent agreement in the small-to-moderate regime" when fitting the Mooney-Rivlin model to uniaxial tension of the Treloar data. Indeed, the stress-strain curve is well fitted in the considered strain regime and deformation mode. However, their parameters lead to c 10 < c 01 , i.e., W∕Ī 1 < W∕Ī 2 what is contrary to the experimental findings explained above. Moreover, these parameters result in an overstiff response under equibiaxial tension.
The authors of the present manuscript are of the opinion that the ratio between the Ī 1 -and Ī 2 -dependency, i.e., W∕Ī 1 and W∕Ī 2 is essential to balance the model response in different deformation modes. This interpretation is motivated by the dependence of the stress response on Ī 2 for each (red lines) the simplified tube model without the Ī 2 -term is fitted simultaneously to ux and bx data; (green lines) the model without the Ī 2 -term is fitted only to the ux data and predicts the response to the bx deformation deformation mode in Table 3. On the one hand, compared to the W∕Ī 1 -contribution, the W∕Ī 2 -term in the constitutive equation is scaled by 1∕ < 1 for uniaxial tension and by 2 > 1 for equibiaxial tension, whereas the W∕Ī 1and W∕Ī 2 -term are weighted equally for simple and pure shear. On the other hand, it should be noted that for these four deformation modes Ī 1 and Ī 2 are not independent of each other, see App. C for the Ī 2 Ī 1 -dependencies. In particular, the relations Ī 1 >Ī 2 for uniaxial tension, Ī 1 <Ī 2 for equibiaxial tension and Ī 1 =Ī 2 for simple as well as pure shear can be shown. That is, uniaxial tension is primarily driven by Ī 1 and equibiaxial tension by Ī 2 whereas the shear deformations are equally affected by both invariants. As a consequence of this predetermined weighting between the Ī 1 -and Ī 2 -contribution in different deformation modes, test data from only one experiment is insufficient for the calibration of models with Ī 1 -and Ī 2 -dependency because a plausible extrapolation to other deformation modes cannot be guaranteed.
This conclusion should not be taken as a drawback of invariant-based models compared with stretch-based models. In fact, leaving out the Ī 2 -dependency when fitting to only one deformation mode is a possibility to avoid unforeseen behavior in the other deformation modes. In contrast, there is no straightforward method for stretch-based models to ensure in advance a reasonable behavior for other deformation modes. This is due to the fact that both terms � (̄1) and � (̄3) always appear in the constitutive equations, cf. Table 3.
To illustrate the considerations regarding the Ī 2 -dependency, the simplified tube model by [81] is employed with G c , G e and n being the crosslink modulus, entanglement modulus and elastically active chain length, respectively. It depends on both principal invariants and is calibrated on the Treloar data (with an Euclidean norm, normalized error, 1st Piola-Kirchhoff stress) in four different ways, see Fig. 6. First, the model is fitted to the uniaxial as well as equibiaxial data. Then, only the uniaxial data are considered for the parameter identification and the calibrated model is used to predict the equibiaxial response. Afterwards, this procedure is repeated with G e = 0 such that the Ī 2 -dependency is turned off. It can be observed that the model with Ī 2 -dependency is generally able to fit the uniaxial and equibiaxial data very well. When the equibiaxial data are not considered for the model calibration, the uniaxial fit at small strains ( < 4 ) improves. However, the equibiaxial response is drastically overestimated. In contrast, the model without Ī 2 -dependency cannot capture both deformation modes simultaneously since the ratio between the model's response under uniaxial tension and equibiaxial tension is fixed and does not match with the experimental data. Thereby, the equibiaxial prediction is less random and still acceptable compared to the prediction with both Ī 1 -and Ī 2 -dependency. A similar finding was described by Seibert & Schöche [88] who attest "pure Ī 1 -dependent formulations [...] the best predictive properties of multiaxial deformation states when the material parameters were determined on the basis of uniaxial data" only. Also Boyce & Arruda [13] studied the predictive capability of hyperelastic models and observed that "strain energy expressions which contain the second invariant [...], Ī 2 , should be used with caution; forms such as the Mooney-Rivlin model are found to be overly stiff in certain types of deformation".

Effect of the Stress Measure, Norm Order and Error Type in the Cost Function
As mentioned in Sect. 2.8, there are several ways to define the cost function for the fitting procedure. To illustrate the effect of different stress measures and norm orders, the 8-chain model by [3] with the inverse Langevin function approximation by [19] is fitted to the data of the CR+S compound. Furthermore, to illustrate the influence of the error types, the Treloar data [95] are employed. Fig. 7 shows the influence of the stress measure on the fitting result. As explained in Sect. 2.8, different stress measures emphasize different sections of the stress-strain curve. The optimization w.r.t. the Cauchy stress leads to the smallest deviations at large stretches ( > 7.5 ) whereas the 2nd Piola-Kirchhoff stress provides the smallest errors in the range 1.75 < < 6.5 . The 1st Piola-Kirchhoff stress is a tradeoff with a stress-strain curve in between. Note that this and the following effects can be more or less pronounced depending on the experimental data and the model capability, e.g., a perfectly fitting model would be unaffected by the definition of the cost function.
In Fig. 8, the order of the p-norm in the cost function Eq. (28) is varied. The taxicab norm with p = 1 optimizes the mean absolute error, i.e., high deviations at large strains can be compensated by a good fit at low strains. In contrast, the maximum norm ( p = ∞ ) takes only the maximum error of the whole strain range into account. As a consequence, the deviations of the stress response at very larges stretches ( > 9.5 ) are much lower in exchange for a worsened fit at < 7 . The result with the Euclidean norm ( p = 2 ) is a Fig. 9 Effect of the error type in the cost function on the fitting result (the black solid lines are the experimental data from Treloar [95]) Fig. 10 Fitting results for the BIIR+S compound: (a) ranking sorted by RMSE, (b) stress-stretch curve of the best performing models with five, four, three and two parameters (experimental data are shown as grey, solid line; number of parameters is given in brackets) Fig. 11 Fitting results for the BR+S compound: (a) ranking sorted by RMSE, (b) stress-stretch curve of the best performing models with five, four, three and two parameters (experimental data are shown as grey, solid line; number of parameters is given in brackets)

Fig. 12
Fitting results for the CR+P compound: (a) ranking sorted by RMSE, (b) stress-stretch curve of the best performing models with five, four, three and two parameters (experimental data are shown as grey, solid line; number of parameters is given in brackets) Fig. 13 Fitting results for the CR+S compound: (a) ranking sorted by RMSE, (b) stress-stretch curve of the best performing models with five, four, three and two parameters (experimental data are shown as grey, solid line; number of parameters is given in brackets)

Fig. 14
Fitting results for the IR+S compound: (a) ranking sorted by RMSE, (b) stress-stretch curve of the best performing models with five, four, three and two parameters (experimental data are shown as grey, solid line; number of parameters is given in brackets) Fig. 15 Fitting results for the NBR+S compound: (a) ranking sorted by RMSE, (b) stress-stretch curve of the best performing models with five, four, three and two parameters (experimental data are shown as grey, solid line; number of parameters is given in brackets)

Fig. 16
Fitting results for the NR+S compound: (a) ranking sorted by RMSE, (b) stress-stretch curve of the best performing models with five, four, three and two parameters (experimental data are shown as grey, solid line; number of parameters is given in brackets) Fig. 17 Fitting results for the SBR+P compound: (a) ranking sorted by RMSE, (b) stress-stretch curve of the best performing models with five, four, three and two parameters (experimental data are shown as grey, solid line; number of parameters is given in brackets) compromise with a slight bias towards the result obtained with p = 1.
Finally, the effect of the error type in the residual is shown in Fig. 9. Here, the 8-chain model is simultaneously fitted to all three deformation modes of the Treloar data using an absolute, normalized and relative error, cf. Eqs. (29)-(31). In general, it can be observed that the chosen model is not able to capture the stress response of the three experiments simultaneously, see also Sect. 4.1.2. More specific, the uniaxial response tends to be overestimated and the equibiaxial response underestimated. Using an absolute error, the fit to the uniaxial tensile test is slightly better than for the other error types, especially for > 2 , but at the expense of the goodness of fit to the equibiaxial data. Vice versa, the relative error leads to better results for the equibiaxial experiment but worse results for the uniaxial one, see explanation in Sect. 2.8. Moreover, it improves the results for small strains ( < 1.5 ) of all deformation modes. The result obtained with the normalized error lies in between the other results. Note that, if one compares the effect of the error type by means of only one experiment, the absolute error and normalized error would coincide.
In summary, the formulation of the cost function has an essential influence on the rankings in the subsequent sections. Therefore, this choice must be made carefully. In the present manuscript, the • Euclidean norm • Normalized error • 1st Piola-Kirchhoff stress are preferred and employed in all parameter optimizations and numerical investigations due to the following reasons. The Euclidean norm leads to a least square problem for which many specialized, efficient algorithms exist. Moreover, many statistical measures are based on the squared error, e.g., the coefficient of determination. The normalized error gives an equal weight to all deformation modes and does not overemphasize small stresses from very low and possibly noisy forces. The 1st Piola-Kirchhoff stress scales linearly w.r.t. the measured forces and does not need a conversion based on the strain measurement or assumption of perfect incompressibility. Thus, it allows a direct interpretation when testing or simulating the forces on the final rubber product. Furthermore, in contrast to the 2nd Piola-Kirchhoff stress, it is readily accessible for all deformation modes and does not lack of a physical meaning. Unlike the Cauchy stress, it does not give too much weight to the large strain regime which is of lower importance for most engineering applications.

Fitting of ̄I 1 -Dependent Strain Energy Functions
The fitting results of the strain energy functions depending only on Ī 1 , cf. Tables 7 and 8, are given as bar charts in  18(b). These models are chosen for plotting rather than the overall top three models since the stressstretch curves of the best ranked models nearly coincide for the most compounds and are hardly distinguishable.
Comparing the most promising candidates for each compound, it can be observed that there is not a best model which is generally applicable to all rubbers. However, counting the top three appearances, there are some strain energy functions with a high potential: the network averaging tube [58], Gregory [35], Beda [10], extended Yeoh [103] and Hoss and Marczak [47] model. The first one is a physically motivated model with only three parameters. It ranks first for two materials (BIIR+S, IR+S) and provides only for the CR+S compound unsatisfying results. In contrast, the others are phenomenological models with five parameters showing acceptable results for all compounds. However, for each material, there exist four-or three-parameter strain energy functions with a fitting quality nearly as good as the aforementioned ones. In general, models with a lower number of parameters are preferable to avoid overparametrization. That is, they are prone to non-sensitive, non-unique and numerically undesirable parameters with a strong dependence on the initial guess, high correlations and a large number of iterations within in the fitting procedure.
Even models with two parameters are acceptable for some compounds. For example, they are suitable for materials that can stand comparatively low strains and exhibit no upturn, e.g., the peroxide-crosslinked rubbers. For the CR+P compound in Fig. 12, yet the Neo-Hooke model can be reasonable. Furthermore, the response of some sulfur-crosslinked compounds with large strains and a pronounced upturn can be captured acceptably with only two parameters, e.g., the BIIR+S, CR+S, NBR+S and SBR+S compound. Notably, in the authors' opinion, the stress-strain curves of these materials do not have any obvious property in common.
Looking into more detail for the sulfur-crosslinked compounds, two frequent shortcomings even of top-ranked models can be noticed. On the one hand, models typically fail to fit accurately the concave shape at small to moderate strains leading to an underestimation of the stress response in this region, see for instance Fig. 15(b). On the other hand, the material behavior close to the maximum strain is not reproduced properly by numerous models, e.g., the Swanson model in Fig. 14(b). Hence, the maximum error usually lies either in that small to moderate strain region or close to the failure strain. Since rubber parts are rarely designed for large deformations close to failure and the sample variance is high at very large strains, the former is considered as more critical. As discussed in Sect. 4.1.3, a residual with a 2nd Piola-Kirchhoff stress or a relative error can ease the shortcoming in the small strain regime to some extent.
The "a posteriori" check of monotonicity conditions, as discussed in Sect. 4.1.1, is passed by all strain energy functions except the extended Yeoh model. It violates monotonicity of the 1st Piola-Kirchhoff stress when calibrated for the IR+S compound. This issue is due to its exponential term, cf. Table 7, which is responsible for the good fit of the aforementioned concave shape. However, the exponential reduces faster than the polynomial terms increase such that the stress-stretch curve temporarily shows a slightly negative slope when changing from the concave to convex curvature.

Fitting of ̄I 2 -Dependent Strain Energy Functions
For the identification of promising Ī 2 -dependent strain energy functions, the natural rubber from Treloar [95] is used since experimental data from multiple biaxial tensile tests are needed, cf. Sect. 4.1.2. First, all Ī 1 -models from Tables 7 and 8 were fitted only to the uniaxial Treloar data to identify the best performing two-, three-and four-parameter Ī 1 -model. These are the strain energy functions by Xiang et al. [100], Khiêm & Itskov [58] and Yeoh & Fleming [104] which show a very similar goodness of Fig. 19 Fitting of the best performing two-, three-and four-parameter Ī 1 -model to the uniaxial tension data by Treloar [95] and their predictions for pure shear and equibiaxial tension (the stresses of the predictions are shifted by 1 and 2 MPa , respectively, for the sake of readability) fit, cf. Fig. 19 for the resulting stress-strain curves. Each of the three Ī 1 -models is then combined with all Ī 2 -models from Table 9.
The results are depicted in Fig. 20. The additional Ī 2 -dependency improves the fitting quality significantly for all three considered Ī 1 -models. However, the degree of improvement depends on the combination of the Ī 1 -and Ī 2parts. For example, the Ī 1 -model by Xiang et al. [100] almost always achieves the smallest improvement whereas the Network averaging tube model [58] benefits strongly from all Ī 2 -parts. Only when combined with the Alexander model, the Xiang model shows the largest improvement. Since such interdependencies between the Ī 1 -and Ī 2 -parts can hardly be foreseen, general conclusions for all Ī 1 -Ī 2 -combinations must be made carefully. More precisely, it is also possible that an Ī 1 -model which performs poorly on its own may show great results with highly different parameter values when combined with a proper Ī 2 -part. Fig. 20 a Ranking of the Ī 2 -dependent strain energy functions for the Treloar data [95] (number of parameters is given in brackets); fitting of the best performing one-and two-parameter Ī 2 -model as well as the Ī 2 -independent approach combined with the best ranked Ī 1 -models with b two, c three and c four parameters Another desirable feature of the Ī 2 -models is the improved fit of the concave shape at low to moderate strains ( < 2 ) of the uniaxial fits, cf. Fig. 19 and Sect. 4.2. This shortcoming is largely compensated when a Ī 2 -part is added, cf. green and orange lines in Figs. 20(b-d). In case of the Yeoh & Fleming Ī 1 -model, this effect can lead to anomalous parameters in the exponential term (i.e., negative B or large A, cf. Table 7) since this term is already designed to reproduce the behavior in this strain regime and, hence, loses its original purpose.
For all Ī 1 -models, the Ī 2 -part by Chevalier & Marco [18] provides the best improvement, followed by the Klingbeil & Shield approach [60]. Both strain energy functions require two additional parameters. For the latter, the exponent m determines whether W iso,2 is a concave ( m ≤ 1 ) or convex ( m ≥ 1 ) function. This parameter is fitted to m = 0.46 (Xiang), m = 0.60 (Network averaging tube) and m = 0.48 (Yeoh & Fleming), i.e., W iso,2 is always concave. The two best performing one-parameter Ī 2 -models are special cases of this approach with m = 1∕2 [16] and m = 1∕3 [21] and, hence, concave functions, too. Furthermore, the aforementioned, top-ranked model by Chevalier & Marco [18] is concave as it is the sum of the two special cases m = 1 (Mooney-Rivlin, linear) and m → 0 (Gent & Thomas, concave). Interestingly, these two cases on their own show only moderate results. Hartmann & Neff [38] introduced the exponent m = 3∕2 to ensure polyconvexity for all deformation states but at the expense of fitting quality. In case of the Mansouri & Darijani model [67], the parameter a was fitted to small negative values (order of magnitude −10 −3 ) implying slight concavity of W iso,2 . Also, c 02 of the Haines & Wilson model [36] was optimized to small negative values ( −10 −6 ) resulting in a negligibly concave Ī 2 -part. The exponent of the Lambert-Diani & Rey model [63] tends to the lower bound m → 1 , i.e., it reduces from a strictly convex function to the linear Mooney-Rivlin approach. Apparently, a concave Ī 2 -dependency is essential for a good fit to several deformation modes. To prove this conclusion, the model , with a symmetric dependency on Ī 1 and Ī 2 was fitted to all three deformation modes of the Treloar data. The optimized parameters fulfill c 10 > 0 , m 1 > 1 , c 01 > 0 , 0 < m 2 < 1 , i.e., a convex Ī 1 -and a concave Ī 2 -dependency. The section closes with a comparison of the invariantbased models in Figs. 20(b-d) to the widely used, stretchbased Ogden model [78], viz., cf. Eq. (23). Implementations with two, three and four terms ( n = 2, 3, 4 ), i.e., four, six and eight parameters are considered. The fitting results in Fig. 21 reveal that at least three terms are needed for an accurate fit. Comparing the threeand four-term implementation, both exhibit similar stress responses for equibiaxial tension and pure shear. Only for the uniaxial load case at moderate strains ( 2 < < 4 ) and at very large strains ( > 7.5 ), the model with four terms is superior. However, the four-term implementation with eight parameters shows signs of overparametrization, e.g., non-reproducible results from different initial guesses. In contrast, all model combinations in Figs. 20(b-d) with both Ī 1 -and Ī 2 -dependency can reproduce the experiments over the full strain range in all three deformation modes with less and unique material parameters. For instance, the Yeoh & Fleming model combined with the Chevalier & Marco Ī 2 -term requires six parameters and the combination of the Xiang with the Carroll model leads to only three parameters in total.

Fitting of Volumetric Strain Energy Functions
The fitting of volumetric strain energy functions W vol is more straightforward than the calibration of the isochoric models in the previous section. Due to the nearly linear material behavior, i.e., slightly convex pressure-volume curves, cf. Fig. 1(b), most of the models in the literature have one or two parameters. Moreover, W vol is a function of only one argument J. That allows, in contrast to the isochoric strain

Fig. 21
Fitting of the Ogden model with two, three and four terms to the Treloar data [95] (number of parameters is given in brackets) energy functions, a simple "a priori" imposition of reasonable parameter limits. In total, 23 models were found in the literature, cf. [84], but only nine of them, cf. Table 10, fulfill the following requirements: • Strict convexity, i.e., a positive bulk modulus (even for J → ∞ ): K = 2 W vol ∕ J 2 > 0 • Non-finite response when approaching zero volume: J → 0 ⇒ W vol → ∞ • Stress-free state at undeformed configuration: J = 1 ⇒ W vol ∕ J = 0 (implies that if W vol is strict convex, then J = 1 is its global minimum) • Maximum two parameters Although it violates the second requirement, the widely used strain energy with a constant bulk modulus is taken as reference. For the two-parameter models, a lower bound for the second parameter , cf. Table 10, is prescribed to ensure that the aforementioned restrictions are always fulfilled. It is worth mentioning that the first three requirements are numerically motivated rather than physically. More specific, the validity of the considered models is bounded on the one hand by cavitation damage which can be observed for J > 1 , see for instance [27]. On the other hand, hydrostatic compression ( J < 1 ) can lead to pressure-induced glass transition [17]. Both phenomena cause a significant change of the material properties and cannot be captured by the models. Hence, a discussion for J → ∞ or J → 0 is actually meaningless from a physical point of view. Moreover, if cavitation  Table 10 Volumetric strain energy functions (n is the number of material parameters with restrictions on the second parameter ; the initial bulk modulus K 0 must be positive for all models; comments on specific models: 6) both citations originally present three-parameter strain energy functions which are reduced here to the same two-parameter model, see also [83,87] damage is to be modeled, the strict convexity requirement can be violated and has to be abandoned. The results are very similar for all compounds. Therefore, they are exemplified by means of the IR+S compound in Fig. 22 and those for the remaining materials are given in App. E. For all materials, the two-parameter models by Doll & Schweizerhof [26]/Schröder & Neff [86] and Neff et al. [77] show the best results. However, for both models, the second parameter takes large values ( > 18 for Doll/ Schröder and > 57 for Neff, cf. Table 10). Since this parameter appears in the exponent and adjusts the curvature of the pressure-volume curve, it can readily affect the stability of finite element simulations. Hence, specialized elements are required, see for instance [87]. The third twoparameter model by Ogden [79] reduces for all materials to the Simo & Taylor II model with the limiting case → 2 .
Here, greater -values would generate a less convex or even concave curvature of the pressure-volume curve.
The one-parameter models fail to reproduce the curvature of the experimental data and show moderate fitting quality. For instance, the top-ranked Ansys model [2] and Hartmann & Neff model [38] are special cases of the Doll/Schröder model by fixing the second parameter. However, the fixed -values are too small to properly reproduce the experimental data (fixed = 4 and = 5 vs. fitted = 18 … 42 ). In general, the gain in fitting quality of the one-parameter models in comparison to the reference approach with a constant bulk modulus is comparatively small. Hence, in practice, the possibly little improvement of the fit must be weighted up against the numerical robustness of a constant bulk modulus. Moreover, the confined compression test is always subjected to friction, cf. [83,50], such that the modulus and the curvature of the experimental pressure-volume curve tend to be overestimated.

Conclusion
The present manuscript provides a guideline for the process of model selection and calibration of hyperelastic models for rubber materials. For this purpose, the paper begins with a discussion on crucial aspects of the preparation of the experimental data. That includes the proposition of physically motivated procedures for the preload correction in tensile and compression tests. Moreover, a consistency check for multiple experiments in different deformation modes is derived. It was shown that the presented preload correction is a suitable approach to compensate an improper preload. The consistency check was applied to the widely used Treloar data [95] proving a good agreement between its uniaxial tensile, pure shear and equibiaxial tensile test.
Then, the importance of the Ī 2 -contribution to the strain energy function for an accurate and plausible model behavior in several deformation modes is pointed out and illustrated. However, models depending on both Ī 1 and Ī 2 show a limited predictivity to other deformation modes when fitted to only one experiment. Furthermore, the effect of different cost function formulations (i.e., stress measure, error type and order of the residual norm) on the fitting result is demonstrated. This choice can strongly affect the fitting results and cause different rankings. It allows to design tailor-made cost functions which emphasize particular strain regimes or deformation modes.
Based on this preliminary work, 33 strain energy functions in terms of Ī 1 were fitted to uniaxial tension data of nine different rubber compounds. The models with the largest number of material parameters are able to reproduce the stress response of most of the materials, e.g., the extended Yeoh [103], Beda [10], Gregory [35] or Hoss & Marczak [47] model with five parameters each. However, they often tend to numerically undesirable parameters and a slow optimization procedure. In contrast, strain energy functions with three or four parameters, like the network average tube model [58], show satisfying results for specific compounds but are less generally applicable. Two-parameter models provide acceptable fittings for the two peroxide cross-linked rubbers, which exhibit the smallest elongation at break, and for a few of the sulfur cross-linked compounds.
For the identification of promising Ī 2 -based strain energy functions, the Treloar data were considered. The bestperforming one-, two and three-parameter Ī 1 -model were combined with twelve Ī 2 -approaches. It was observed that concave Ī 2 -dependent strain energy functions with a moderate curvature are well suited to improve the goodness of fit, e.g., the Chevalier & Marco [18] or Carroll model [16] with two and one parameter, respectively. More precisely, the additional Ī 2 -dependency is able to balance the stress response of different deformation modes and can compensate shortcomings of a pure Ī 1 -dependency. However, since the test framework considers those Ī 1 -models which perform already well on their own, these findings do not allow general statements about preferable Ī 2 -approaches.
Finally, modeling approaches for the volumetric behavior were analyzed. For this purpose, nine strain energy functions are fitted to the confined compression data of the nine abovementioned rubber compounds. In contrast to the stress-strain curves of the uniaxial tensile test, the pressure-volume curves of all materials are very similar. Hence, suitable volumetric strain energy functions are ranked similarly, too. Clearly, the twoparameter models Doll/Schröder [26,86] and Neff [77] provide the best results. However, their pronounced non-linearity may lead to additional numerical challenges.
For future work, an expansion of the presented database is reasonable. More precisely, up-to-date data sets comprised of uniaxial tensile, pure shear and equibiaxial tensile tests on several compounds are desirable. This would allow well-founded investigations of further Ī 1 -Ī 2 model combinations, of coupled Ī 1 -Ī 2 -terms and of stretch-based models. Alternatively, general biaxial tension data could be employed for these purposes. In that case, discussions on the preparation and on the consistency check of the data as well as on the treatment of two stress values in the cost function are necessary. Furthermore, the accuracy and error sources of the biaxial experimental setup should be assessed.