Analysis of the compressible, isotropic, neo‑Hookean hyperelastic model

The most widely-used representation of the compressible, isotropic, neo-Hookean hyperelastic model is considered in this paper. The version under investigation is that which is implemented in the commercial finite element software ABAQUS, ANSYS and COMSOL. Transverse stretch solutions are obtained for the following homogeneous deformations: uniaxial loading, equibiaxial loading in plane stress, and uniaxial loading in plane strain. The ground-state Poisson’s ratio is used to parameterize the constitutive model, and stress solutions are computed numerically for the physically permitted range of its values. Despite its broad application to a number of engineering problems, the physical limitations of the model, particularly in the small to moderate stretch regimes, are not explored. In this work, we describe and analyze results and make some critical observations, underlining the model’s advantages and limitations. For example, a snap-back feature of the transverse stretch is identified in uniaxial compression, a physically undesirable behavior unless validated by experimental data. The domain of this non-unique solution is determined in terms of the ground-state Poisson’s ratio and the state of stretch and stress. The analyses we perform are essential to enable the understanding of the characteristics of the standard, compressible, isotropic, neo-Hookean model used in ABAQUS, ANSYS and COMSOL. In addition, our results provide a framework for the parameter-fitting procedure needed to characterize this standard, compressible, isotropic neo-Hookean model in terms of experimental data

incompressible approximation for those materials is often adopted for simplicity. We note that there exist numerous isotropic, incompressible, hyperelastic models, and researchers continue to propose new models to represent elastic deformations in the finite strain regime more accurately. Comparisons of existing models can be found in various papers [5][6][7][8]. It is difficult to choose a particular model as the best one to represent the nonlinear elastic behavior of different rubber-like materials. There are simpler models with fewer material parameters, but highly complex phenomenological versions also exist involving numerous material parameters. Some models are proposed for moderate but finite strain, whereas others can be used for very high strains at the cost of additional material parameters. The incompressible assumption can be used to obtain closed-form expressions of the transverse stretches for most models during standard homogeneous loadings (uniaxial loading, equibiaxial loading in plane stress, uniaxial loading in plane strain, pure shear loading, simple shear loading). Consequently, the stretch-stress relations for these incompressible, hyperelastic models can be expressed in closed-form analytical formulae, which can then be used in a parameter-fitting process to characterize a given rubber-like material in terms of experimental data.
The simplest isotropic, incompressible, hyperelastic model is undoubtedly the neo-Hookean formulation that was proposed by Rivlin [9]. It contains only one material parameter, and is utilized mainly in the moderate strain regime. It should be emphasized that it is widely used due to its simplicity. Although most researchers use Rivlin's original representation, there exists a new, two-parameter generalized version [10,11].
In cases where material compressibility cannot be neglected, the incompressible, isotropic, neo-Hookean hyperelastic model can be extended by addition of a term in the strain energy accounting for volumetric deformation. The compressible version of the resulting neo-Hookean model is not unique, and there exist several forms in the literature. The particular version implemented in the commercial finite element (FE) software ABAQUS [12], ANSYS [13], and COMSOL [14] is the most widely-used and widelyaccepted. We reference this model in the following as the standard, compressible, isotropic, neo-Hookean model. However, it should be borne in mind that other models have been suggested [15][16][17][18][19].
Compressible, isotropic, neo-Hookean hyperelastic models (CNH) are used for various applications. Brock and Hanson [20] examine crack growth in a CNH material, deriving analytical formulae for wave speeds in the crack plane of their problem. It is important to emphasize that these authors use the CNH model for a highly compressible material having ground-state Poisson's ratio of 0.25. Horgan and Saccomandi [16] investigate the variation of the apparent Poisson's ratio, the Poisson function, in terms of the applied stretch for uniaxial tension in a nonlinear CNH model. They find that the Poisson function decreases with increasing stretch; the largest groundstate Poisson's ratio they use is 0.45. The internal stability of a CNH model in homogeneous deformation is studied by Clayton and Bliss [17], using the formulation proposed by Simo and Pister [21]. This stability analysis is performed for both positive and negative ground-state Poisson's ratios. Begley, Creton and McMeeking [22] obtain the displacement solution around a Mode I crack tip in plane strain for a CNH model, while Budday, Raybaud and Kuhl [23] model human brain tissue in the form of a stretch-driven growing inner core and a morphogenetically growing outer surface. The cortex of the brain is modeled using a CNH formulation with ground-state Poisson's ratio of 0.458, and the authors infer that mechanical deformation plays a crucial role in brain development. Their material parameters are also used by Morin et al. [24]. Zhang [25] derives a lattice model to represent the mechanical behavior of a CNH material, demonstrating the lattice model's performance and accuracy for ground-state Poisson's ratios ranging from 0.417 to 0.495. Behera, Roy and Madenci [18] use a CNH model with a ground-state Poisson's ration equal to 0.452 to study large deformation and rupture in rubber-like materials, employing the bondassociated weak form of peridynamics with nonuniform horizon in their analysis. We note that CNH models can be extended to characterize transversely isotropic porous solids. For example, Guo and Caner [26] use a novel anisotropic constitutive model developed for materials having aligned cylindrical pores, and Chen and Ravi-Chandar [27] develop an innovative CNH constitutive model to capture the nonlinear poroviscoelastic behavior of a gelatin-based hydrogel.
The complete understanding of the standard CNH model is essential not only for mechanical engineers but also for graphics and computer vison engineers working in computer graphics. In this application, a fundamental goal is to simulate natural phenomena realistically. This task includes the simulation of the movement and deformation of characters composed of soft solids (e.g., humans, animals, imaginary creatures). Bender et al. [28] utilize a CNH model to simulate the deformation of deformable solids and clothing, paying particular attention to the simulation of cloth, where achieving realistic wrinkling is crucial. The CNH model is used by these authors for a very low ground-state Poisson's ratio, namely 0.1, demonstrating that CNH models are used for materials with significant compressibility. Smith, De Goes and Kim [19] conclude that to capture the fleshy appearance of virtual characters it is essential to use hyperelastic material models. These authors propose a new version of the compressible neo-Hookean model, and demonstrate that the new model behaves well for ground-state Poisson's ratios in the range from 0 to 0.5. We note that Smith et al. [19] emphasize that biological tissues are nearly incompressible, but numerical simulations within this regime are extremely challenging. Sheen, Larionov and Pai [29] underline the popularity of the CNH model in computer graphics for enforcing incompressibility using a ground-state Poisson's ratio close to 0.5. These authors comment that accurate simulation of the human body in contact with its environment is crucial in visual effects and apparel design. They propose a novel approach to alleviate numerical instabilities and locking in coarsely discretized bodies, and demonstrate the performance of their CNH material model in numerical examples. Macklin and Müller [30] emphasize that the neo-Hookean model has recently become popular in computer graphics, and employ the standard CNH formulation as implemented in ABAQUS, ANSYS and COMSOL. They present numerical examples with a range of ground-state Poisson's ratios from 0.4 to 0.4995.
Although the expression for the strain energy associated with the standard, compressible, isotropic, neo-Hookean model is relatively simple, closed-form stress solutions cannot be obtained for homogeneous deformations, not even for simple uniaxial tension. Consequently, the determination of material parameters via a parameter-fitting procedure using experimental data is more complicated than for the incompressible models. One must use a numerical solution strategy to calculate the transverse stretch and the longitudinal stress in uniaxial loading. Furthermore, stress/stretch solutions for the standard CNH model in homogeneous loading, e.g., uniaxial loading, equibiaxial loading in plane stress, uniaxial loading in plane strain, are not yet reported in the literature. There is no paper presenting the effect of compressibility on stress/stretch solutions and on the uniqueness properties of the model. The theory guides and manuals for the commercial FE software ABAQUS, ANSYS, and COMSOL do not provide detailed information about the behavior of the compressible neo-Hookean model implemented in these codes. In particular, the effect of the value of the initial ground-state Poisson's ratio on stress/stretch solutions is unclear, and the domains of uniqueness for stress and stretch are not identified. This is a serious challenge for the application of such models to understand experimental data obtained using compressible material systems, since, as will be demonstrated below, the standard CNH material model can yield multiple stress-stretch solutions in certain ranges of the applied stretch in uniaxial compression, and exhibits some other behavior that is unphysical unless validated by comparison with experimental data. In addition, the uniqueness, or its lack, and the stability of hyperelastic models are important, and must be understood completely in order to avoid physically unrealistic responses. Peng and Li [31] identified stable domains of volumetric deformation for a few compressible, hyperelastic material models, including an extended CNH model differing from the standard version. However, analysis of the standard CNH model for the standard loading cases is not included in their investigation. Pence and Gou [15], provide a detailed analysis of the behavior of three different CNH models, but exclude the standard CNH model. The authors also address the stability analysis of the non-standard CNH models only in uniaxial loading.
This paper provides a detailed analysis of the characteristics of the standard CNH model. Three commonly used, simple, homogeneous deformations (uniaxial loading, equibiaxial loading in plane stress, uniaxial loading in plane strain) are considered. It is shown that the transverse stretches parallel to zero stress directions cannot be expressed in closed-form, and are obtained using numerical calculations. The ground-state Poisson's ratio, 0 , is used to quantify the ratio of the two material parameters involved in the constitutive model, and to parameterize the constitutive model and the results of the analysis. Computational results are presented for the entire permissible range, −1 < 0 < 0.5 . The variation of the transverse stretches in terms of the specified stretches is analyzed for each loading case. It is found that the model can yield multiple stress solutions for the same longitudinal stretch in uniaxial compression. The domain for unique solutions is determined, and stress solutions are calculated for the homogeneous deformations under investigation. These solutions are compared to those for the incompressible material to quantify numerically the differences resulting from the presence of compressibility. These new results give insight into the nature of the standard CNH model. Furthermore, these new solutions can be used in parameterfitting procedures using experimental data, where compressibility of the material is not negligible.

Constitutive equations for the standard, compressible, isotropic, neo-Hookean model
The strain energy for isotropic, hyperelastic solids is expressed with use of the principal invariants I 1 , I 2 , I 3 of the left Cauchy-Green deformation tensor b = FF T [2]: where F is the deformation gradient, and i (with i = 1, 2, 3) are the principal stretches, whereas the determinant of the deformation gradient defines the volume ratio as J = det(F) . The derivatives of the principal scalar invariants with respect to b read The Cauchy stress tensor is defined by the relation U ,i (with i = 1, 2, 3) denote the partial derivative of U with respect to the ith principal scalar invariant. Using the relation U ,3 = U ,J ∕(2J) and Eq. (5), we find that the Cauchy stress tensor can be formulated as The modified principal invariants are given by the relations The decoupled representation of the strain energy function is assumed to be written as where the term U dev represents the strain energy resulting from distortional (volume preserving) deformation and U vol denotes the energy associated with volumetric deformation. If the distortional part of U is independent of I 2 then the Cauchy stress solution reduces to The most common form of the strain energy for the compressible, isotropic, neo-Hookean (CNH) model (implemented in ABAQUS, ANSYS and COMSOL) is written as where G and K are two material parameters, equal to the ground-state shear and bulk moduli, respectively, as discussed in "Appendix 2". The partial derivatives required for this model are Consequently, the general form of the Cauchy stress tensor for the CNH model is given by where as the first Piola-Kirchhoff stress tensor can be obtained from the Cauchy stress tensor as If a deformation is imposed that is volume preserving,(J = 1) , the result in Eq. (13) . This is a valid prediction of the stress arising from volume preserving loading such as occurs in the idealized simple shear deformation of the double lap shear test. However, due to end effects, homogeneous deformation does not occur in practice in the double lap shear test. Furthermore, the single lap shear test is complicated by an applied moment as well as end effects, and thus this test lacks homogeneous deformation. Consequently, the analysis of the CNH model for lap shear tests is excluded from the paper, and we restrict our attention to tests that can have long gauge lengths or large gauge areas in which homogeneous deformations occur. These cases are uniaxial loading, equibiaxial loading in plane stress, and uniaxial loading in plane strain, where, in each case, volume change is involved during deformation. We note that the analysis of uniaxial loading in plane strain is a valid treatment of the deformation arising in the test of what is known as the pure shear specimen, as long as its height (i.e., in the direction subject to applied load) is large compared to its unconstrained thickness, and the length of the specimen is large compared to its height.

Uniaxial loading
Uniaxial loading is depicted in Fig. 1a. The deformation gradient and the left Cauchy-Green deformation tensors are expressed with the help of the principal stretches as The longitudinal stretch is and the transverse stretch is denoted T . The volume ratio and the first principal invariant are J = 2 T and I 1 = 2 + 2 2 T . Inserting these quantities into Eq. (13), (14) we can express the nominal stress components as: The zero transverse stress constraint, P 22 = P 33 = 0 , cannot be solved analytically to obtain a closed-form expression for the stretch relation T = T ( ) . However, numerical schemes can be used to compute the transverse stretch. The solutions depend on the ratio K∕G , i.e. on the groundstate Poisson's ratio 0 . The definition of the groundstate Poisson's ratio is explained in "Appendix 2".  One of the most critical observations in the results is the occurrence of "snap-back" of the transverse stretch in the compression regime for values of 0 > 0.296 . It can be seen that in compression for < 0.4 there are multiple solutions for the transverse stretch for a given longitudinal stretch . This domain is therefore a regime of non-unique and possibly unstable behavior. The results show that if a material with such a non-unique characteristic is subjected to a uniaxial stress in compression with the longitudinal stretch being controlled monotonically, the component will at first thicken in the transverse direction. However, at a critical longitudinal stretch, the material will snap from being transversely thicker than its original width to being thinner than it. However, the exact behavior will depend on whether the longitudinal stretch or the axial load is being controlled, with different behavior occurring when the axial load is controlled monotonically, as will be described below. Note that for all values of ground-state Poisson's ratio other than 0.5 the material eventually experiences transverse thinning for small longitudinal stretch ratios. Therefore, a positive value of ground-state Poisson ratio, other than 0.5, does not guarantee that transverse thickening will always occur in axial compression. However, for ground-state Poisson's ratio values less than 0.296 the transition from transverse thickening to transverse thinning occurs smoothly and with continuity, and there is no snap-through. Note also that for longitudinal stretch approaching zero, the transverse stretch also approaches zero other than for 0 = 0.5 . Therefore, the volume of the material approaches zero in this circumstance, in violation of quantum exclusion.
If 0 < 0.296 the phenomenon of multiple solutions is absent. The domain for multiple solutions is determined in the − 0 plane and is shown in Fig. 3.
The coordinates at the bottom-right corner of the domain in Fig. 3 are (0.3967, 0.296) . Consequently, if the ground-state Poisson's ratio is smaller than 0.296 ( K∕G ≈ 2.118 ) then the model provides a unique solution in the entire regime of the possible stretches.
In uniaxial tension, the model for a material with a positive ground-state Poisson ratio behaves as  Fig. 2c where the transverse stretch is plotted versus the longitudinal stretch out to very large axial strains. The largest value for the ground-state Poisson's ratio used in the plot is 0.4999 (red curve), whereas the smallest value is − 0.99 (purple curve). The black dashed line represents T = 1 . For negative ground-state Poisson's ratios, Fig. 2c shows that the transverse stretch has a maximum when the longitudinal stretch is increased, and thereafter decreases monotonically. Eventually, the plots for negative ground-state Poisson's ratios cross the axis for T = 1 and the material begins to thin as the longitudinal stretch is increased. The longitudinal stretch values at which T = 1 in uniaxial tension are shown in Fig. 4 for different values of 0 .
Once the transverse stretch is determined, the nominal stress component P 11 can be obtained using Eq. (16). The result for P 11 ∕G is shown in Fig. 5a for different values of 0 . The largest value for the ground-state Poisson's value used in the plot is 0.4999 (red curve), whereas the smallest value is − 0.99 (purple curve). It can be observed there are domains on the compression side of Fig. 5a where multiple stress solutions exist for the same longitudinal stretch value. This plot reveals that if the longitudinal stretch is controlled monotonically, there will be a "snap-back" feature in compression in some cases as is reduced. It is also observed that the limit values for the nominal stress at = 0 is zero other than for the incompressible case, where infinite stress occurs. Due to quantum exclusion, a zero value for the nominal stress at = 0 is physically unrealistic as it would predict that zero force is needed to compress the material completely. Note that when the axial load is controlled monotonically during uniaxial stressing, other than for the incompressible material, the behavior is such that a  Fig. 5a goes through a minimum. For magnitudes of the nominal stress in compression greater than the stress magnitude at the minima shown in Fig. 5a, the solution for stress ceases to exist. Such a behavior implies, though does not model, material failure.
We introduce the dimensionless parameter measuring the ratio of the compressible P comp 11 and incompressible P inc 11 nominal stress results. The incompressible stress solutions are summarized in "Appendix 1". With the help of this parameter, we can illustrate the effect of compressibility on the stress. The ratio is shown in Fig. 5b as a function of for different values of 0 . The largest value for the ground-state Poisson's ratio used in the plot is 0.4999 (red curve), whereas the smallest value is -0.99 (purple curve). It is noted that is undetermined at = 1 , but its limit can be calculated. One can conclude that is always smaller than 1.
Results for the Cauchy stress are provided in "Appendix 3".

Equibiaxial loading in plane stress
The schematic for equibiaxial loading in plane stress is shown in Fig. 1b. The deformation gradient and the left Cauchy-Green deformation tensors are expressed as The volume ratio and the first principal invariant are J = 2 T and I 1 = 2 2 + 2 T . Inserting the kinematic quantities in Eq. (13), (14) we can express the nominal stress components as: The zero transverse stress constraint, P 33 = 0 , cannot be solved in closed-form for the transverse stretch, T , as a function of the equibiaxial stretch, . Instead, the transverse stretch is obtained using a numerical method, and the results are presented in Fig. 6 It is important to observe that the solution for transverse stretch is unique for each value of in the entire domain of 0 , as shown in Fig. 6. However, there is a local maximum for the transverse stretch as is reduced in the compression regime if 0 is positive. As is reduced in compression, the transverse stretch eventually transitions so that the transverse thickness of the material has reduced below T = 1 instead of T being greater than unity, i.e., thinning instead of thickening. In addition, the transverse stretch tends to zero as → 0 . This is a physically unrealistic phenomenon unless validated by experiment, as in the uniaxial stress case. A local maximum in T also exists in tension in Fig. 6 for negative ground-state Poisson's ratios.
The equibiaxial nominal stress is shown in Fig. 7a as a function of the equibiaxial stretch. The largest value for the ground-state Poisson's ratio used in the plot is 0.49999 (red curve), whereas the smallest value is − 0.9999 (purple curve). In this case the equibiaxial nominal stress is unique for a given equibiaxial stretch. However, the equibiaxial stretch is not unique for a given equibiaxial nominal stress, i.e., the stress-stretch relationship is not monotonic. It can be seen that there is a minimum value for the equibiaxial nominal stress in compression, similar to the uniaxial loading case. Below this minimum value, the solution for equibiaxial nominal stress ceases to exist, with implications similar to those holding in the uniaxial loading case. It is notable that the limit equibiaxial nominal stress at = 0 is zero, as in the uniaxial stress case.
The ratio of the compressible stress solution and the incompressible stress solution for equibiaxial stretch in plane stress is denoted by the parameter as in the uniaxial stress case. The variation of as a function of biaxial stretch is shown in Fig. 7b. The largest value for the ground-state Poisson's ratio used in the plot is 0.49999 (red curve), whereas the smallest value is − 0.9999 (purple curve). The results in Fig. 7b reveal that the stress solutions in the compressible case are always smaller in magnitude than the incompressible stress solutions.
Results for the Cauchy stress are provided in "Appendix 3".

Uniaxial loading in plane strain
Uniaxial loading in plane strain is depicted in Fig. 1c It is axial stressing subject to plane strain constraint in direction 2. The deformation gradient and the left Cauchy-Green deformation tensors are expressed as The volume ratio and the first principal invariant are J = T and I 1 = 1 + 2 + 2 T . Substituting the kinematic variables into Eq. (13), (14), we can formulate the nominal stress components as: The boundary conditions imply that P 33 = 0 . Due to its highly nonlinear nature, this equation cannot be solved in closed-form for the transverse stretch. We use a numerical scheme to determine T for a given and the results are plotted in Fig. 8. For 0 > 0 , we see the same phenomenon in compression as in the uniaxial and equibiaxial cases: the transverse stretch increases at first as is reduced, goes through a maximum and then decreases. It is concluded that the transverse stretch tends to a finite value for = 0 . This limiting value is T = 1∕ √ 2 , in contrast to the zero limit observed in uniaxial loading and in equibiaxial loading in plane stress.
The axial nominal stress, P 11 , is shown in Fig. 9a as a function of the axial stretch. The largest value for the ground-state Poisson's ratio used in the plot is 0.49999 (red curve), whereas the smallest value is -0.9999 (purple curve). These stress versus stretch (23) This ratio for uniaxial loading in plane strain is shown in Fig. 9b. The largest value for the ground-state Poisson's ratio used in the plot is 0.49999 (red curve), whereas the smallest value is − 0.9999 (purple curve). The stresses in the compressible case are smaller in magnitude than in the incompressible case in the entire domain.
The numerical results for the transverse stress, P 22 , are plotted in Fig. 10a. The largest value for the ground-state Poisson's ratio used in the plot is 0.49999 (red curve), whereas the smallest value is -0.9999 (purple curve). In axial compression for positive values of the ground-state Poisson's ratios, P 22 at first decreases as is reduced, goes through a minimum and then increases. For positive values of the ground-state Poisson ratio, the limit value for P 22 is +∞ at = 0 . This result is physically unrealistic as a negative transverse stress is expected in axial compression when the ground-state Poisson's ratio is positive.
The transverse nominal stress in axial tension in plane strain also exhibits unrealistic characteristics, but at relatively large axial stretch. Figure 10b presents the results for P 22 for higher axial stretches. The largest value for the ground-state Poisson's ratio used in the plot is 0.499999 (red curve), whereas the smallest value is − 0.999999 (purple curve). In the incompressible case, P 22 approaches G as → ∞ , consistent with the analytical solution in the limit of Eq. (30). If material compressibility is allowed, then the limiting value of P 22 is zero as → ∞ . As a result, as is increased from unity, P 22 at first increases, goes through a maximum and then decreases to zero as → ∞. Results for the Cauchy stress are provided in "Appendix 3".

Discussion
These results have important implications for the application of compressible, isotropic, neo-Hookean hyperelastic material models to the analysis of experimental data. The ability to robustly model the deformations of compressible materials in the moderate strain regime is required for broad classes of engineering materials, including most biological and biomimetic materials, which typically have a Poisson ratio in the range of 0.2-0.5 and exhibit nonlinear elastic responses [33][34][35][36][37][38]. The ability to accurately relate observed deformations to applied stress fields is particularly important for the development of traction force microscopy methods, which allow determination of cell-generated forces to understand the fundamental mechanobiology, as well as improved disease modeling and design of soft tissue replacements [39][40][41][42][43][44][45].
In this study, the behavior of the standard version of the compressible, isotropic, neo-Hookean hyperelastic material model for modeling such complex materials is analyzed. In uniaxial loading conditions, a stable regime producing unique solutions is found for 0 < 0.296 . For larger values of 0 , multiple solutions are observed in a particular domain. This suggests the possibility of an unstable, but physically realizable, regime that may provide opportunities to leverage snap-through instabilities, which have been demonstrated to be useful in actuation with applications to soft robotics [46][47][48]. In the case of equibiaxial loading in plane stress, the nominal stress is unique for a given stretch; however, the stress-stretch relationship is not monotonic, just as it was not monotonic in the uniaxial loading case. In the case of uniaxial loading in plane strain the transverse stretch tends to a finite value T = 1∕ √ 2 and the axial nominal stress tends toward −∞ as approaches 0. This contrasts with the zero limit observed in uniaxial loading and in equibiaxial loading in plane stress. For uniaxial loading in plane strain, the numerical results for the transverse stress are not meaningful across all conditions, however. In particular, the limit value of the transverse stress at zero stretch is physically unrealistic for positive values of the ground-state Poisson's ratio, and the transverse nominal stress in axial tension exhibits unrealistic characteristics at sufficiently large axial stretch values. Taken together, these observations place important bounds on the conditions under which such models produce physically meaningful results. The limit values of the transverse stretch, the nominal stress and the Cauchy stress are summarized in Fig. 11 for each loading case. The monotonicity of each result is also indicated in the table. Note that the Cauchy and nominal stress components, 22 and P 22 , are equal in uniaxial loading with plane strain. This stress component tends to +∞ as approaches 0, whereas the limit value is 0 for → ∞ as discussed in the previous section.

Conclusions
The standard version of the compressible, isotropic, neo-Hookean hyperelastic material model is considered in this paper. It is demonstrated that analytical, closed-form solutions for the model for transverse stretch and axial stress cannot be obtained for simple, homogenous deformation modes due to the Funding Open access funding provided by Budapest University of Technology and Economics.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that influenced or could have appeared to influence the work reported in this paper.

Appendix 1: Results for the incompressible case
The strain energy in the incompressible case is defined as The closed-form Cauchy stress solutions, including the added hydrostatic stress, are where U , B , P are the stresses in uniaxial loading, equibiaxial loading in plane stress and uniaxial loading in plane strain. The stress P is the axial stress while PC denotes the transverse stress in the direction of plane strain constraint. The limit value for PC is G at = ∞.
The closed-form nominal stress solutions are given by formulae (29) and (30) and shown in Fig. 12 Appendix 2: Linearized solutions in uniaxial loading Inserting Eq. (15) in Eq. (13), we can express the Cauchy stress components as Linearization of these stress components are formulated as (29) where the engineering strains e = − 1 and e T = T − 1 are introduced for simplicity. The zero transverse stress constraint ̃2 2 = 0 is solved for the transverse engineering strain, giving: The parameter 0 denotes the ground-state Poisson's ratio, which is related to the ratio of parameters K and G as Insertion of the transverse strain Eq. (37) into Eq. (35) gives the stress along the loading direction: where E 0 represents the ground-state Young's modulus. Combining Eq. (37) and Eq. (35), one finds that the material parameters in the CNH model are the ground-state shear and bulk moduli: G = G 0 and K = K 0 . The dependence of the ratio K∕G on 0 is visualized in Fig. 13. The limit value for ratio K∕G is 0 at 0 = −1 , whereas it is ∞ at 0 = −0.5 (incompressible case).
Some numerical values are summarized in Table  1. The value 0 = 0.495 (where K∕G is ~ 100) has a particular importance because ABAQUS/Standard recommends the use of solid continuum hybrid elements to avoid difficulties with convergence for materials with ground-state (or initial) Poisson's ratio greater than this value. Furthermore, if the Poisson's ratio is set to 0.5 (incompressible case) for an elastic material in ABAQUS/Standard, then the software changes it to 0.4999999 (where K∕G is ~ 5 × 10 6 ).

Appendix 3: Cauchy stress solutions
Once the nominal stress components are obtained for the compressible cases, the Cauchy stress components are also calculated according to the formula in Eq. (14). The Cauchy stress 11 in uniaxial loading is shown in Fig. 14a for different values of 0 . The largest value for the ground-state Poisson's ratio used in the plot is 0.4999 (red curve), whereas the smallest value is − 0.99 (purple curve). It can be observed there are domains on the compression side of Fig. 14a, where multiple Cauchy stress solutions exist for the same applied longitudinal stretch value. This plot reveals that if the longitudinal stretch is controlled monotonically, there will be a jump in Cauchy stress in compression in some cases as is reduced. However, if the axial Cauchy stress is controlled, the behavior will be smooth and continuous in compression as far as is concerned, although in some cases the material will go through a transient stage of incremental stretching even as the compressive stress magnitude is increased. It is also observed that the limit values for the Cauchy stress at = 0 is finite other than for the incompressible case. Due to quantum exclusion, a finite value for the Cauchy stress at = 0 is physically unrealistic. Numerical analysis reveals that the limiting value for 11 at = 0 in each case is exactly −3K.   Cauchy stress results for equibiaxial loading in plane stress are plotted in Fig. 14b. The largest value for the ground-state Poisson's ratio used in the plot is 0.49999 (red curve), whereas the smallest value is − 0.9999 (purple curve). It can be concluded that the equibiaxial Cauchy stress is unique for a given equibiaxial stretch. Similarly, the equibiaxial stretch is also unique for a given equibiaxial stress, i.e., the stress-stretch relationship is monotonic for the Cauchy stress. It is interesting to note that the limit Cauchy stress at = 0 is finite, as in the uniaxial loading case (see Fig. 14a). Numerical analysis reveals that the limiting value for 11 = 22 at = 0 is exactly −3K∕2.
The axial Cauchy stress, 11 , is shown in Fig. 13c as a function of the axial stretch in uniaxial loading in plane strain. The largest value for the groundstate Poisson's ratio used in the plot is 0.49999 (red curve), whereas the smallest value is − 0.9999 (purple curve). These Cauchy stress versus stretch plots are monotonic and the stress tends toward −∞ as → 0 . This limiting behavior is physically realistic. The transverse Cauchy stress, 22 , is identical to the transverse nominal stress, P 22 , in uniaxial loading in plane strain. Therefore, the discussions and conclusions for the transverse nominal stress presented in Sect. 3.3 are also valid for the transverse Cauchy stress.