A review on material models for isotropic hyperelasticity

Dozens of hyperelastic models have been formulated and have been extremely handy in understanding the complex mechanical behavior of materials that exhibit hyperelastic behavior (characterized by large nonlinear elastic deformations that are completely recoverable) such as elastomers, polymers, and even biological tissues. These models are indispensable in the design of complex engineering components such as engine mounts and structural bearings in the automotive and aerospace industries and vibration isolators and shock absorbers in mechanical systems. Particularly, the problem of vibration control in mechanical system dynamics is extremely important and, therefore, knowledge of accurate hyperelastic models facilitates optimum designs and the development of three‐dimensional finite element system dynamics for studying the large and nonlinear deformation behavior. This review work intends to enhance the knowledge of 15 of the most commonly used hyperelastic models and consequently help design engineers and scientists make informed decisions on the right ones to use. For each of the models, expressions for the strain‐energy function and the Cauchy stress for both arbitrary loading assuming compressibility and each of the three loading modes (uniaxial tension, equibiaxial tension, and pure shear) assuming incompressibility are provided. Furthermore, the stress–strain or stress–stretch plots of the model's predictions in each of the loading modes are compared with that of the classical experimental data of Treloar and the coefficient of determination is utilized as a measure of the model's predictive ability. Lastly, a ranking scheme is proposed based on the model's ability to predict each of the loading modes with minimum deviations and the overall coefficient of determination.


| INTRODUCTION
Elastomers (derived from elastic polymers) are polymeric materials that exhibit both viscous and elastic behavior. Characteristics of these materials, which are also known as rubber-like materials as they exhibit rubber-like properties, include high reversible deformations (up to 500% strains) without fracture, high damping properties, low thermal and electrical conductivity, durability, and hysteresis under cyclic loading. 1 The most important property is its large nonlinear elastic deformation leading to their categorization as hyperelastic/ green elastic materials. These properties are highly desirable in a myriad of engineering applications such as engine mountings, structural bearings, vibration absorbers, corrosion protection, tires, medical devices, shock isolators, and springs. 1 Both hyperelasticity and linear elasticity are reversible processes (no internal energy dissipation) meaning that the work done during the loading process is completely recovered when the load is removed.
The difference is that for hyperelastic materials, the deformation can be extremely large and the relationship between the stress and strain is nonlinear. Derivation of the stress-strain relationship for the hyperelastic materials is based on a function known as Helmholtz free energy per unit reference volume (Ψ). 4 It is also known as strain energy density, strain energy function, or the elastic potential and is a scalarvalued function that relates the strain energy to the state of deformation. The unique trait of hyperelastic materials is that the strain energy density is dependent only on the current strain and not on the loading history. 5 The formulation of hyperelastic material models begins with the development of a suitable strain energy density function.
Several assumptions are adhered to in deducing the strain energy function. These include that the material is isotropic, homogenous, free of hysteresis, strain-rate independent, and nearly or purely incompressible. 6 Furthermore, various restrictions on the strain energy density such as being nonnegative, having a zero value at undeformed state, and being invariant under coordinate transformations have been proposed by several texts. 7 Over the years, dozens of constitutive models for hyperelastic materials have been formulated by various research groups. Based on the approach with which the models are formulated, they can be categorized as either phenomenological or micromechanical. As the name suggests, phenomenological-based models arise from the observation of rubber-like materials under different conditions of homogenous deformation and thereafter fitting mathematical equations to the experimental data. 8 These equations result in polynomial formulations which further classify the phenomenological models into subcategories based on strain invariants, principal stretches, or a combination of both. On the other hand, the micromechanical models exploit the techniques of statistical mechanics to describe the behavior of hyperelastic materials at the microscopic and macroscopic levels. Even though the phenomenological models make up a larger percentage of the hyperelastic models in the literature, micromechanical-based models have got more attention thanks to their governing parameters that can relate the mechanical behavior with the physical or chemical structure of the material. 9 Apart from the synthesis of new hyperelastic materials that require new constitutive relations to predict their mechanical behavior, formulating models that accurately describe the complete behavior of the material remains the main motivating factor to researchers. By complete behavior, we mean the capability of the model to predict the behavior under uniaxial extension, biaxial extension, and pure shear. 10 Treloar in his paper, 11 provided the experimental data that has been the main reference for researchers in testing new models and their predictive abilities in all three modes of loading. The challenges involving the material models include their predictive ability and the number of parameters required. Most models may predict the behavior in one mode of loading such as the uniaxial extension to good accuracy and fail in the other two modes. Again, some models require many material parameters that may be difficult or timeconsuming to achieve through experimental means. It is desirable to have a model that can reproduce experimental data in all the loading modes and has the least number of parameters that are easy to obtain experimentally. Choosing a model for a specific application is determined by several factors such as its accuracy, available or achievable parameters, ease of implementation, material type, loading type, and so on.
It is worth noting that most of the new model developments are modified versions of the already established models to improve their predictive abilities. For instance, Arruda-Boyce's eight-chain (EC) model 12 which is micromechanically inspired and one of the most accurate models, is known to be relatively inaccurate in predicting the behavior during biaxial extension loading. As such, several modified versions of the EC model have been reported in the literature some of which were subjected to a comparative study by Hossain et al. 13 where they found that the bootstrapped EC model by Miroshnychenko  With dozens of models in the literature, it is imperative to the design engineers to understand the type of models to utilize in their simulations depending on the material type, accuracy, loading type, computational cost, availability of parameters, and so on. A review of the models provides the engineers and upcoming researchers the opportunity to hasten their understanding of the same. There have been remarkable reviews in the literature about constitutive models for rubber-like materials. One of the most referenced reviews on rubber elasticity models was authored by Boyce and Arruda 21 more than two decades ago. They reviewed a total of 11 rubber elasticity constitutive models that are formulated based on both the statistical mechanics and invariant or stretch-based continuum mechanics theories. Importantly, they included a discussion on two common approaches of modeling the compressibility of elastomers. Another frequently cited review work on hyperelastic models was authored by Marckmann and Verron 10 in the year 2006 in which they compared the performance of 20 hyperelastic models based on their ability to reproduce experimental data in different types of loading conditions for two sets of classical experimental data.
Moreover, they proposed a ranking scheme to list the models from the best to the worst based on a number of indicators namely the capability of the model to achieve complete behavior, the number of required material parameters, the ability of the model to reproduce both sets of experimental data with the same set of parameters, and whether the parameters are physically motivated. Models which could reproduce experimental behavior in different loading conditions, required few parameters, reproduced both sets of data without changing parameters, and whose parameters have physical meaning were highly ranked. In 2019, Dal et al. 22 presented a comparative study of 40 hyperelastic models while focusing on the parameter identification process. They proposed a novel parameter identification toolbox based on a multiobjective optimization technique that selects the best model and its parameters from an input of uniaxial tension, equibiaxial extension, and pure shear experimental data. The models were then ranked based on the quality of fit for simultaneous fits for two sets of experimental data.
Most recently, Dal et al. 23 extended the previously mentioned work into a comprehensive state of the art involving 44 hyperelastic models. The ranking of the models was not only based on the quality of fit for simultaneous fits but also the number of material parameters and the validity range.
Considering the growing interest in hyperelastic materials and the role played by numerical simulation in product design, this review work aims to bring about the knowledge of the most commonly utilized material models for hyperelastic materials concisely for the benefit of design engineers and researchers.
In this study, while presenting our discussions, we classify the models into the two mentioned main categories namely phenomenological and micromechanical where the former is further subdivided into those based on invariants or principal stretches.
The discussion on each model mainly focuses on the strain energy density expression, nominal stress-strain or stress-stretch expressions for arbitrary or specific loading conditions (uniaxial tension, equibiaxial extension, or pure shear), material parameters, and their predictive capabilities according to Treloar's data. Illustrations that compare the experimental and the predicted stress-strain or stress-stretch curves are presented and the accuracy is quantified using the coefficient of determination.
Importantly, this study proposes a novel ranking scheme that considers the behavior of the model in each loading mode, the overall behavior, and the deviation of the fitness coefficient for each loading mode in the experimental data set.

| PHENOMENOLOGICAL MODELS
Constitutive models emanating from the phenomenological approach are formulated by fitting mathematical equations to the experimentally observed behavior of the material. 24 The formulation considers the macroscopic nature of the material hence treating the problem from the continuum mechanics viewpoint. 25 There are two categories of phenomenological constitutive models; those that are based on the invariants of the Cauchy-Green deformation tensors as introduced by Rivlin 26  | 73 the real behavior of elastomeric materials. Practically, these materials undergo volume changes, especially when subjected to hydrostatic deformations. Consequently, the most accurate constitutive models must account for compressibility. Not only is it important for accuracy but also for avoiding numerical problems during finite element implementation of the model. To include a compressibility term, the Ψ expression is additively decomposed into volumetric and deviatoric parts which are responsible for volume and shape changes respectively as shown in Equation (1).
where the subscripts d and v represent the deviatoric and volumetric parts respectively. J is the volume ratio and is obtained by the determinant of the deformation gradient J F = det whereas I* 1 and I* 2 This study utilizes the most commonly used expression for the volumetric part according to 30 and is given by K is the bulk modulus of the material.

| Neo-Hookean (NH)
It was firstly formulated by the renowned rubber expert, Ronald . It is an extension of Hook's law for linear elasticity to include large deformations where the stress-strain relationship is nonlinear. It is the most commonly used and well-known hyperelastic model due to its simplicity as it requires only two material parameters that can be easily determined; the shear modulus μ ( ) and the bulk modulus κ ( ). Of the two material parameters, the former controls the deviatoric response whereas the latter controls the volumetric response. 31 If incompressibility is assumed, that is, no volume change during deformation (as is frequently done in theoretical calculations), the expression becomes even simpler as only the shear modulus parameter is required. As with all the hyperelastic constitutive equations, the NH model is derived from its corresponding expression for strain energy density function (Ψ) as shown in Equation (3).
It is worth noting that several authors have slightly varying expressions for Ψ. This study is based on the expressions by Bergström. 32 The terms in Equation (3) are as follows; μ = shear modulus, where F is the deformation gradient, I* 1 = tr C ( *), where J C C * = (−2/3) (distortional part of the right Cauchy-Green deformation tensor), k = bulk modulus. By assuming an incompressible deformation, J = 1, thus reducing the energy expression to I μ I Ψ( *) = /2( * − 3) 1 1 . From the energy expression, the NH model's Cauchy stress for an arbitrary loading mode is obtained as shown in Equation (4).
where J B B * = (−2/3) (distortional part of the left Cauchy-Green deformation tensor). When incompressibility of the material is assumed, the Cauchy stress for the three loading modes; uniaxial, planar, and biaxial are given in Equation (5).
The subscripts u p , , and b represent uniaxial, planar, and biaxial loading respectively whereas λ represent the stretch applied.
The accuracy of the NH model, as with any other hyperelastic model, can be approximated by comparing its stress-strain predictions to that of the classical experimental data by Treloar 11 as shown in Figure 1A. The coefficient of determination R ( ) 2 as shown in Equation (6) where n = number of data points, e i = experimental data at a point i, p = i prediction data at a point i, and e m = the mean of the experimental data.

| Yeoh
While experimenting on the behavior of carbon black reinforced rubber, Yeoh 33,34 observed that the shear modulus significantly dropped at low strains. The available hyperelastic models at that time could not predict this behavior accurately. He proposed the addition of exponentially decaying terms to the strain energy function resulting in an expression that is a polynomial of I* 1 without the dependence of the second invariant as shown in Equation (7). The work of Kawabata et al. 35 found that the Helmholtz free energy for most rubber-like materials is less dependent on I* 2 than on I* 1 . Furthermore, it is more strenuous to determine by experimental means how I* 2 influences the energy function. This motivated Yeoh to drop the dependence of I* 2 in his expression.

∑ ∑
The terms C i0 and k m in Equation (7) above are the material parameters to be determined whereas N is an input integer representing the order of the polynomial.
If three orders are considered in the strain energy function (hence N = 3), the Cauchy stress considering compressibility in arbitrary loading mode as given by Bergström 32 is shown in Equation (8).
On the other hand, Cauchy stress expression when incompressibility is assumed for each of the three loading modes is given in Equation (9).
It is worth noting that if only the first order is considered, that is,  36 and has a reputation for predicting the response of hyperelastic materials to a high level of accuracy, therefore, it is well known and most commonly preferred model. Such improved accuracy is a result of the inclusion of linear dependence on I* 2 in the strain energy function. This means that the deviatoric response is defined by both the first and the second invariant. One way of expressing the strain energy function is shown in Equation (10).
where C 10 , C 01 , and k are the material parameters required for this Considering compressibility during deformation, the Cauchy stress expression for arbitrary loading is obtained as shown in Equation (11).
The incompressible forms of Cauchy stress for uniaxial, planar, and biaxial loading are given in Equation (12).
The Mooney-Rivlin model works well from small to medium strains. The plots which are shown in Figure 1C compare the model's stress-strain predictions with the experimental results of Treloar. The model has improved accuracy relative to that of the NH model.

| Gent
The Gent model 37  The original work by Gent assumed incompressibility in formulating the strain energy density expression. The expression shown in Equation (13) has been modified by Bergström 32 to include a com- The Cauchy stress expression obtained for arbitrary loading and considering compressibility is given in Equation (14) σ From Equation (14), it is worth noting that as J m tends towards positive infinity, the Gent model's Cauchy stress expression becomes equivalent to that of NH. The incompressible forms of Cauchy stresses for the Gent model are given in Equation (15).
The Gent model has the advantage of being able to predict large strain loadings (up to 300%), higher accuracy, and also requires only three parameters. As shown in Figure 1D, the Gent model does relatively well in predicting the material behavior and its accuracy is close to that of the Yeoh model.

| Isihara
The model was formulated by Isihara et al. 38 and it is considered as a specialized form of the Mooney-Rivlin model since its strain energy density function (see Equation 16) includes the dependence on the second invariant.
The material parameters required by the incompressible version are C C , 10 20 , and C 01 . As per the original work, the uniaxial tension Cauchy stress expression for an incompressible deformation is shown in Equation (17).
The predictive capability of the Isihara model is shown in Figure 2A where it is compared with the classical experimental data. It can be observed that uniaxial tension loading predictions are accurate to moderate stretch values. The model poorly estimates the behavior in biaxial loading mode.

| Horgan and Saccomandi
The researchers Horgan and Saccomandi 40 aimed to improve the predictive ability of the Gent model by introducing the dependence on I* 2 in the strain energy density function. Their work resulted in a hyperelastic model in their names, Horgan and Saccomandi which we shall shorten to HS for convenience. The difference between the limiting chain stretches in Gent and HS models is that the former depends on the maximum value of the first invariant whereas the latter on the maximum allowable stretch thus more physical significance. 32 Strain energy expression that assumes compressibility for the HS model is given in Equation (18).
The parameter λ max represents the limiting chain stretch.
HS's model Cauchy stress expression derived from the strain energy density function is given in Equation (19).
As in Gent's model, the HS equates to the NH as the limiting chain stretch approaches infinity. Considering incompressibility, the Cauchy stress for uniaxial tension, planar shear, and equibiaxial tension loading modes are given in Equation (20).
One main advantage of the HS is its stability provided that the shear modulus and limiting chain stretch are greater than zero and one,  (21).
where A, B, and C are material parameters found by Carroll to be 0.15, 3.1 × 10 7 , and 0.095, respectively (units = MPa).
According to Steinmann 9, the expressions for the stress-stretch relationship for the three loading modes considering incompressibility are given in Equation (22).
where λ λ , 1 2 , and λ 3 are the principal extension ratios. As in the other hyperelastic materials, the expressions for stress can be obtained by taking the partial derivative of Ψ with respect to Consequently, the components of Cauchy stress in arbitrary loading are given as shown in Equation (24). 41 where i = 1, 2, 3.
The Valanis-Landel model requires only one material parameter for the incompressible form which is the shear modulus. This makes the model simple but reduced accuracy in predicting the material behavior as shown in Figure 2D.

| Ogden
Formulated by Ogden 28 The parameters are as follows: N = the order of the model (normally taken as 3), μ k and α k 2 are material constants, and D k is a parameter that indicates volume change. The expression for the principal stresses (σ i , = 1, 2, 3 i ) is shown in Equation (26).
[(( *) + ( *) ) + ( *) ] An interesting observation in Equation (26) is that when both N and α k are equal to 1, then the Ogden model equates to the NH model. Considering incompressibility, the stress expressions for the three loading modes (subscripts μ p , , and b representing uniaxial, planar, and biaxial loading, respectively) are given in Equation (27).
The predicted stress-stretch response of the Ogden model has been found to agree very well with the classical experimental results of Treloar as shown in Figure 4A. The model is highly suitable for predicting large deformation behavior and can describe well the sharp increase in stiffness at large strains. The main drawback of this model is that material parameters are specific to every deformation mode and, therefore, require different sets of parameters.

| Shariff
In his work, Shariff 42  and most importantly being linear in its material parameters as shown in Equation (28).
where i = 0, 1, 2, 3, …, n, α i = material parameters, α = 1 0 , φ i = smooth function that takes up to four expansion points as shown in Equation (29) according to Badienia. 39 To obtain the expressions for the Cauchy stress, the first derivative of Equation (28)  given as shown in Equation (30).
The Shariff model requires five material parameters which are E α α α , , , 1 2 3 , and α 4 . As shown in Figure 4B, the model accurately captures the observed experimental behavior, particularly, the S-shaped curve behavior at high strains. While extra effort is needed for parameter fitting due to the complexity of the constraint equations, less time is required compared to that of the Ogden model since its energy function is a function of linear parameters. 43

| MICROMECHANICAL MODELS
Micromechanical models arise from analyzing the deformation behavior of rubber-like materials from the microstructural point of view. As it is well known, the microstructure of rubber-like materials consists of randomly oriented long polymeric chains that are joined together into a network structure. 44 Furthermore, the polymeric chains are comprised of N rigid beams which are commonly known as Kuhn segments each of equal length l. 39 The ideal maximum length of the chain r max at full elongation is given by r Nl = max . 9 Statistical mechanics theory (the random walk theory 45 ) is utilized to describe the average end-to-end distance of stress-free undeformed chain which is found as r N l = . The distribution of the end-to-end distance r of a polymeric chain is given by the Gaussian probability distribution function shown in  Langevin function is utilized to account for the finite chain extensibility and the resulting force-extension expression is given in The parameters are as follows; k = Boltzmann's constant, .
To formulate a micromechanical-based constitutive model, a representative network structure is necessary to relate the chain stretch of the individual polymeric chains to the applied deformation.
Various micromechanical-inspired hyperelastic models with their corresponding network structures are presented in this section.

| Three-chain model
The representative network structure of the 3-chain model 47,48 has three chains (3-chain) with each located along the axes of the undeformed cubic cell as shown in Figure 3A. When the deformation is applied, the chains will deform in an affine manner with the cubic structure with the stretch on each chain corresponding to the principal stretch value. 21 The non-Gaussian distribution is utilized in deriving the strain energy function expression for the 3-chain model as shown in Equation (34) where λ i = principle stretch in the i th axis, , and N = the total number of chains in the three principal directions.
The stress-stretch expressions for the three loading modes when incompressibility is considered according to Steinmann et al. 9 are given in Equation (35).
The subscripts UT ET , , and PS represents the uniaxial tension, equibiaxial tension, and pure shear, respectively.

The 3-chain model requires two material parameters (μ and N)
for each of the loading modes. The model's predictions compared to Treloar's data are plotted in Figure 4C. It can be observed that the model fits excellently to the uniaxial tension data but fails to capture accurately the biaxial and planar shear behavior.

| Eight-chain model
The stored energy function for the 8-chain model is dependent on the effective distortional stretch and is expressed as shown in Equation (37).
The arbitrary loading Cauchy stress expression for the 8-chain model considering compressibility according to Bergström 32 is shown in Equation (38). (38) The parameter λ lock = the fully extended stretch of the chain (limiting chain stretch). Three material parameters are required which are μ λ , lock and k. Assuming that the material is incompressible, only two material properties will be required (bulk modulus not needed) given in Equation (39).
The prediction of the 8-chain model is plotted against the experimental results of Treloar in Figure 4D. Since it is independent of the second invariant, the accuracy of the biaxial loading prediction is low. Nonetheless, the model is observed to be more accurate than some of the phenomenological models such as Mooney-Rivlin.  According to Hossain and Steinmann, 43 the stress expressions for uniaxial tension, equibiaxial tension, and planar shear loading modes considering incompressibility are given in Equation (41).

| Tube model
A comparison of the tube's model predictions to the experimental results by Treloar is plotted in Figure 5A. The model cannot satisfactorily predict the material's behavior. It is mentioned 43 that the reason is that the energy contribution from the chain cross-linking is of NH type and, therefore, the S-shaped curve which is synonymous with rubber-like materials from moderate to high stretch values cannot be realized by the model.

| Extended tube model
The same authors of the tube model, Kaliske and Heinrich, 50 32 The new strain energy function has three parts as shown in Equation (42).  (43).
where λ J λ * = i i −1/3 and the other terms remain as previously defined.
According to Hossain and Steinmann, 43 the stress expressions derived from the strain energy density function for the three loading modes considering incompressibility are given in Equation (44).
The extended tube model's predictions fit the experimental Treloar data to a very high level of accuracy as shown in Figure 5B. The extended tube model is one of the most accurate hyperelastic material models. The material parameters for uniaxial loading mode can be used to predict the behavior of other loading modes.

| Flory-Erman
The Flory and Erman model 51 begins by assuming that the strain energy density of a polymer network exhibiting high elasticity is expressed as a sum of two contributions. The first contribution is from the phantom network (defined as a network whereby the physical effects of the chains between junctions are confined entirely to the forces they exert on the pairs of the junctions to which each is attached) whereas the second contribution is from the constraints emerging from the material properties of chains that are densely distributed in a random network (see Equation 45).
ph ct (45) The subscripts ph and ct represent phantom and constraint, respectively. Furthermore, the phantom network energy is derived from the Gaussian chain statistics hence equivalent to the neo-Hookean strain energy density as shown in Equation (46) The parameter φ in Equation (46) represents the number of chains at a junction whereas μ is the shear modulus. On the other hand, the expression for the constraint contribution of the strain energy is derived from the micromechanics of the chain molecules as is given in Equation represents the strengths of the constraints. Considering incompressibility, the stress-stretch relationship expressions for the three deformation modes are given in Equation (48).
where the subscript ct represents the stress contribution from the con- constitutive model since it is a crucial prerequisite for good numerical predictions. It is, therefore, imperative to present how each of the models reviewed in this study performs in their predictions.
There is no clearly defined method of ranking the models and, therefore, researchers have employed varying methods. Marckmann and Verron 10 ranked 20 hyperelastic models based on their abilities to reproduce two sets of experimental data by Treloar 11 and Kawabata et al. 52 The models that could predict accurately the three modes of loading in the two sets of data were considered as the best whereas those that required a higher number of material parameters were ranked lowly. Badienia, 39 utilized the quality of fit measure whereby the models that reproduced the S-shaped curve of the Treloar's experimental data stress-stretch plot to a higher level of fit were ranked highly. Using Treloar's data, Bergström 32 ranked the models based on the coefficient of determination and also the number of material parameters required. Applying a slightly different approach, Ritto and Nunes 53 utilized two methods in ranking five hyperelastic materials. First, they ranked according to the error between the predicted and the experimental data considering two loading modes; pure and simple shear. Second, they developed a model based on Bayesian statistics which they found to be comprehensive in ranking a given set of hyperelastic models.
In general, two factors are crucial in determining a model's rank. The first and most important factor is the model's capability in replicating the experimental data in all of the three modes of loading namely uniaxial tension, equibiaxial tension, and pure shear. The second is the number of material parameters required by the model. This influences in that a higher number of parameters means more tedious experimental work is required to calibrate the model. This may not be much of an issue these days since codes that accurately approximate the material parameters for hyperelastic models based on uniaxial loading data have been developed such as the MCalibration. 54 In this section, we will present the predictive performances of the hyperelastic models on each of the three loading modes based on the coefficient of determination. Finally, we will propose a new ranking method that is based on the ability of the model to accurately predict the behavior in the three loading modes taking into account the deviation of the R 2 values. Based on our observation, some models are excellent in predicting the behavior in one mode, for example, uniaxial tension but poorly predict the other modes of loading. We pro- pose that a good model should be able to predict all three modes accurately with minimum differences in the R 2 values.

| Uniaxial tension
As shown in Figure 6, the majority of the models discussed pre-

| Equibiaxial tension
The performance of the models in the equibiaxial tension is better than in the uniaxial tension as shown in Figure 7. Even though the

| Pure shear
All the models reviewed in this study predict the planar shear loading behavior pretty accurately. As shown in Figure 8

| Overall predictive capabilities
To deduce the R 2 value for the overall predictive capability for each of the models, a simple average of the R 2 values for the three loading modes was obtained. Without considering the number of material parameters required, the Ogden, Carroll, Extended tube, and Shariff models are the best with over 0.99 R 2 values as shown in Figure 9. with the same predictive performance, the one with fewer parameters will be preferred and ranked higher.
In this section, we propose a ranking scheme based on two important functions of the models and without considering the number of material parameters required. First, the model should not only be able to reproduce the experimental data in the three loading modes to a higher level of accuracy but should also have minimum deviations (found by calculating the SD) between the R 2 values for each of the loading modes.
This means that a model with high predictive capability in one loading mode and low in the other will be ranked low. Second, the R 2 value is taken into consideration by categorizing those that fall between 0.95-1.0 as index 1, 0.90-0.94 as index 2, 0.85-0.89 as index 3, and so on. The ranking coefficient will be the sum of the SD and the index value of the model. The model with the lowest ranking coefficient (meaning both the SD and the index value are the lowest) is taken as the best. The summary of the ranking is presented in Table 1  Note: ET, E-C, HS, 3-C, N-H, M-R, F-E, and V-L represents extended tube, eight-chain, Horgan-Saccomandi, three-chain, Neo-Hookean, Mooney-Rivlin, Flory-Erman, and Valanis-Landel respectively. The UT, EB, and PS represent uniaxial tension, equibiaxial tension, and pure shear, respectively. potential with respect to the strain invariants or the principal stretches.
The material is assumed to be homogenous and isotropic. Incompressibility is mostly assumed to simplify the stress-strain expressions.
For the interest of research and design engineers, this study employed a unique approach to bring about the knowledge of constitutive models for large nonlinear elastic mechanical behavior of elastomeric materials in a succinct manner. The models were categorized based on the approach in which the strain energy density expression is formulated, phenomenological or micromechanical. The former was further grouped into those that are dependent on the strain invariants of the Cauchy-Green deformation tensors or principal stretches. For each model, emphasis was placed on the strain energy density expression, the stress-strain or stress-stretch expression for arbitrary or specific loading conditions, and the material parameters for the classical Treloar's experimental data. Importantly, the predictive performances of the models were ascertained by comparing the predicted stress-stretch or stress-strain data to the experimental data by Treloar where the coefficient of determination R ( ) 2 was employed to determine the closeness of the predicted and the experimental data. The predictive performance was presented in two ways; first for each loading mode (uniaxial tension, equibiaxial extension, and pure shear) and second for the overall or complete behavior in which the average of the coefficient of determination in each loading mode was utilized. Since the performance of the models in different loading modes varies, it is important to establish both the individual and overall behavior. Except for the single-parameter models with a linear dependence on the first invariant, the predictive performance in the uniaxial tension and pure shear was to an acceptable accuracy level. On the other hand, the equibiaxial loading proved to be the most difficult to reproduce by most models. An excellent model that can be utilized to predict the practical behavior of complex 3D components under arbitrary loading conditions should reproduce the experimental data in all the loading modes accurately. Furthermore, a model with a minimum number of material parameters is highly desirable as more will mean increased complexity of the model and the calibration process may lead to instabilities in the parameters. This study proposed a ranking method whereby the model with a high R 2 value (the average of the R 2 values for the three loading modes) and minimum deviations between the R 2 values for each of the loading modes, is considered as the best.
As the computational capabilities of modern computers continue to soar high, the numerical simulations via finite element codes increasingly become indispensable in the component design process and so is the development of the constitutive models. While it is highly encouraged to have models with fewer material parameters for simplicity, stability, and low computational costs, what will matter in the future is the ability of the model to accurately predict the material behavior under arbitrary loading as there are algorithms that extract material parameters from simple uniaxial loading data. On top of achieving models that accurately describe the experimental behavior in all kinds of loading conditions, future research should focus on addressing some of the main drawbacks of the current models that include the number of experimental data required for calibrating phenomenological models and the range of applicability of the models. For stable parameters, the majority of the models require simultaneous fitting to both uniaxial tension and equibiaxial extension data. This is disadvantageous since it increases the number of experiments, the time, and the cost of calibrating models. Furthermore, the equibiaxial loading experimental setup is complicated and prone to errors. On the range of applicability, some models are accurate but only to a limited strain range, for example, moderate strains. Researchers should strive to achieve models that work in the full strain range.

ACKNOWLEDGMENTS
The scholarship to the first author by China Scholarship Council is highly appreciated. This study was supported by the National Natural Science Foundation of China (Grant no. 11772109).