Constitutive framework of a new hyperelastic model for isotropic rubber-like materials for finite element implementation

The strain energy function of hyperelastic material models must fulfill some mathematical conditions to satisfy requirements such as numerical stability and physically plausible mechanical behavior. In the framework of computer simulations with Newton-type methods, numerical stability is assured by the positive definiteness of the tangent operator. The Baker-Ericksen inequalities, on the other hand, are sufficient and necessary conditions in order to guarantee that the material behaves in a physically plausible way, although they are rarely taken into account during parameter identification. The present work proposes a modification in the strain energy function of a previously developed model for isotropic rubber-like materials. The new expression for W allows the satisfaction of both of the aforementioned requirements. The complete constitutive framework for its implementation in a Finite Element code is provided. Representative examples are analyzed and to show the superior performance when compared to well-known models found in the specialized literature both for homogeneous and non-homogeneous cases of deformation.


INTRODUCTION
The increasing application of rubber-like materials in engineering components and structures throughout the last decades has led engineers to better understand the inherently non-linear mechanical behavior of this class of materials. Since the early works such as that of Mooney (1940), in which he proposed a strain energy function to relate stress and strain for elastomers, several authors have proposed hyperelastic models based either on phenomenological observations or micromechanical considerations.
It is possible to characterize the mechanical behavior of rubber-like materials via hyperelastic models in a wide variety of engineering applications, ranging from conveyor belts, tires, gaskets, hoses, etc. to biological applications such as tendons, muscles, skin, etc. Not only the nonlinearity between stress and strain take place, but usually these materials undergo large deformation in service, which makes it necessary the use of a constitutive framework within the theory of finite elasticity.
The strain energy function of metallic materials under linear elastic regime is well known and assumes a quadratic form (Hooke's law), which leads to a linear relationship between stress and strain. When dealing with hyperelastic materials, however, due to the large variety of different materials, together with the fact that both their chemical composition and manufacturing process strongly affect their mechanical behavior in service, it is virtually impossible to define one single strain energy function that would work to any hyperelastic material. In this sense, the ultimate goal of researchers in the field is to propose a strain energy function that models accurately the largest possible majority of different rubber-like materials, while calibrating them from a set of experimental data and predicting with acceptable accuracy their mechanical behavior when subjected to different deformation modes.
The main objectives of the present communication are (i) to propose a novel constitutive model suitable for isotropic rubber-like materials with the complete necessary framework for Finite Element (FE) implementation together with (ii) the set of restrictions that should be satisfied by its parameters in order to fulfill the Baker-Ericksen (B-E) inequalities (Baker and Ericksen, 1954) and ensure positive definiteness of its consistent tangent operator, as well as to (iii) assess the new model by quantitative comparison with some well-known models, both for homogeneous and nonhomogeneous cases of deformation.
In Section 3 we present the constitutive framework with all the necessary equations for a correct implementation in a code based on the finite element method (FEM). Section 4 brings the mathematical restrictions that must be satisfied to fulfill B-E inequalities and positive definiteness of the tangent operator. In Section 5 we propose a mathematical modification of a strain energy function aiming to comply with the aforementioned conditions together with the necessary equations for a FE implementation. Results are provided for different sets of representative examples. Finally, Section 6 contains a summary of the results and concludes with an assessment of the advantages of the new model.

LITERATURE REVIEW
Rubber elasticity has been governed mostly by statistical (thermodynamic) and phenomenological strain energy function proposals (Treloar, 1975).
Statistical models do not capture well the stiffening effects characteristic of rubber-like materials due to the assumption of Gaussian chain statistics, not valid under large strains (Yeoh and Fleming, 1997). Non-Gaussian statistical (Langevin chain) models led to more realistic proposals for adding finite-extensibility of the polymer chains, especially the well-known four-chain (Wang and Guth, 1952) and eight-chain (Arruda and Boyce, 1993) models. Examples of recent advances in modeling hyperelasticity using this approach can be found in Arora et al. (2020), Morch et al. (2019), Li et al. (2019). Since non-Gaussian models become mathematically very complex when applied to realistic atomic distribution, phenomenological models are still the best options to describe particular behaviors observed in experiments. Yeoh's model (Yeoh, 1993) captures the cubic behavior of shear modulus at small strains by adding an exponential term of the first strain invariant to the classic Rivlin's polynomial series (Rivlin and Saunders, 1951), which adds two constants to the Rivlin's model. Gent's model incorporated the limiting stretch through a logarithmic term of the first strain invariant (Gent, 1996) so that stress at large strains can be better represented, at the cost of a single constant if the limiting stretch is known. A combination of these two models led to the Yeoh-Fleming's model (Yeoh and Fleming, 1997), making it a good choice when the transition from the small strain-softening regime to the large strain-stiffening regime must be captured. The use of mathematical terms empirically added to strain energy proposals has been explored during the last two decades. In particular, exponential terms have been used in Fung's type strain energy function Fung, 1967), Kuhn-Grün type energy function (Beatty, 2003;Boyce and Arruda, 2000), molecular models (Kilian, 1981;Horgan and Saccomandi, 2006;Fried, 2002), in alternative energy functions with asserted convexity (Mansouri and Darijani, 2014;Ghiba et al., 2015) and several known hyperelastic models, to name a few.
However, for elastomers showing significant monotonic stretch between these two characteristic behaviors the Yeoh-Fleming's may fail to fit well experimental results for uniaxial and biaxial tensile tests. One potential solution is the inclusion of power-law type terms in the strain energy function. Power-law functions (Knowles, 1977) have been explored in the past, and are known for reducing to Gent-type strain energy models, besides performing well in non-homogeneous deformations (Horgan and Saccomandi, 2006).
The authors' experience (Hoss et al., 2011) shows that many of the celebrated models fail to provide good theoretical predictions to deformation modes other than the one(s) used to calibrate the model or, and when they succeed, they do not perform accurately in the whole deformation range. This opens room for a hyperelastic model that conjugates exponential, logarithmic, and power-law functions. The original Hoss-Marczak model was proposed for incompressible materials (Hoss et al., 2011) and provided good results for a variety of elastomeric materials, and performs well in all deformation ranges. This model was also extended to include non-conservative effects (Wrubleski and Marczak, 2018).
However, Hoss-Marczak's model has not yet been tested in a finite element environment, having only being used for analytical calculations of cases of homogeneous deformations for parameter identification. In this work, we show that a small (but essential) modification to its strain energy function must be performed in order to simultaneously (i) fulfill the B-E inequalities and (ii) derive a positive-definite tangent operator.
The B-E inequalities aim to guarantee that for an isotropic elastic body the maximum principal stress always occurs in the same direction of the maximum principal stretch. They are a necessary condition for the material to behave in a physically plausible way. According to Neff et al. (2015), "the Baker-Ericksen inequalities are arguably an absolutely necessary requirement for reasonable material behavior". It must be pointed out, however, that cases in which the material undergoes phase transition (for natural rubber, for example, when cold crystallization or stretch-induced crystallization take place), the B-E inequalities could be violated (Marsden and Hughes, 2012). In this paper we consider that these phenomena do not occur during the process of deformation. In Section 4.1 we derive the necessary mathematical conditions for the satisfaction of the B-E inequalities.
In the structure of numerical simulations via Newton-type methods, the mathematical treatment of basic boundary value problems is based on variational methods, and in this context the constitutive equations, represented by the strain energy function must not only reflect the adequate mechanical behavior of the material but also fulfill some general requirements in order to guarantee numerical stability (Balzani et al., 2006). In these cases, one crucial condition to achieve the expected quadratic convergence and numerical stability is the positive-definiteness of the Hessian (here the consistent tangent operator).
In Section 4.3 we show a case in which a well-known hyperelastic model with a positive-definite tangent operator, with accurate calibration to experimental data, which does not satisfy the B-E inequalities. Numerical results show a mechanical behavior that is not physically plausible.
In Section 5 we show that the hyperelastic model proposed by Hoss et al. (2011) is not capable of fulfilling the B-E inequalities and present a positive-definite tangent operator simultaneously. Dealing strictly with curve fitting and parameter identification (analytical solutions), does not impose any problem, but for the purpose of implementing it numerically (as in a finite element code, for instance) this inconsistency must be overcome. This issue is also addressed in Section 5, in which we propose a modification to its strain energy function.

CONSTITUTIVE FRAMEWORK
The framework detailed here is standard hyperelasticity theory from Nonlinear Continuum Theory and further details can be found in specialized textbooks (Holzapfel, 2000;Ogden, 1984). For a body 3    subjected to a finite deformation 3 :     one defines the deformation gradient , in which X is a material point in the reference, undeformed configuration. The deformation gradient F can be multiplicatively decomposed into (Flory, 1961) in which F  refers to the volumetric part and F to the isochoric part, respectively defined as Latin American Journal of Solids and Structures, 2021, 18(1), e346 4/17 (2) (3) From Equation 1, the isochoric right Cauchy-Green tensor is defined as Now, the strain energy function W of a compressible, isotropic hyperelastic model can be additively decoupled into isochoric W and volumetric Ŵ parts depending on the strain invariants of C as It is important to point out that it is also possible to write W in terms of the invariants I B and II B of the left Cauchy-Green deformation tensor B (defined by T  B FF ) and, although C and B are different tensors (and, consequently, C and B too), their invariants are identical, which leads to the fact that the equation for W written in terms of both C or B have the same format.
The second Piola-Kirchhoff stress tensor S is written in terms of the Cauchy stress tensor σ as in which Following the same decomposition criterion of the strain energy function, the second Piola-Kirchhoff stress tensor can be additively decoupled as (Holzapfel, 2000) 2 where S is the isochoric portion of S and S  is the volumetric portion, computed, respectively, by in which  is a fourth-order projection tensor and Š is a fictitious stress tensor respectively, in which The fourth-order tangent operator is given as the Hessian of W . It is proportional to the derivative of S with respect to C and can also be decoupled additively into isochoric and volumetric parts (Holzapfel, 2000): in which the volumetric portion   is computed as with the notation The computation of the isochoric part  of Equation 17 is as follows: with the definitions

The Baker-Ericksen inequalities
The application of the B-E inequalities to the strain energy function of an isotropic material can be derived from Equation 8 which, according to Atkin and Fox (1980), can be written as in which the coefficients i a are scalar functions of the invariants of B . Making use of the Cayley-Hamilton theorem and a few mathematical rearrangements, it is possible to get the most general form of a constitutive relation for an elastic isotropic material subjected to finite deformation (Truesdell and Noll, 1965): in which, assuming incompressibility of the material, 0 the hydrostatic pressure). Equation 28 can be written as (Atkin and Fox, 1980;Bonnet and Wood, 1997) Recalling that the B-E inequalities guarantee that for an isotropic elastic body the maximum principal stress i  always occurs in the same direction of the maximum principal stretch i  , it means mathematically that which can be written in terms of the coefficients i b of Equation 28 for isotropic incompressible materials as in which i and j and k are permutations of 1 , 2 and 3 . Taking into account that 1 2 Truesdell and Noll (1965) conclude by inspection that in order to satisfy the B-E inequalities

Positive-definiteness of the tangent operator
One important requirement for the numerical treatment of the finite element method using Newton-type methods to present stability and quadratic convergence is the positive definiteness of the Hessian, i.e. the tangent operator  . Writing  in terms of the strain invariants I C , II C and III C its Hessian is computed according to 2 2 2 2 2 2 I II I whose positive defineteness is assured if its determinant is positive, or

Violation of the B-E inequalities
In this section we show a case of a simple hyperelastic structure, modeled with the 5-terms Mooney-Rivlin model (Rivlin and Saunders, 1951) I  3  II  3  I  3 II  3  I  3  II 3 we identify the parameters to fit the experimental data presented in Meier et al. (2003) for biaxial extension of a silicone elastomer (Figure 1). The parameters obtained for Equation 36   A one-mesh element was subjected to an unsymmetrical biaxial stretch according to the boundary conditions depicted in Figure 3 (left). In the horizontal axis, the element was stretched up to a 150% of strain, while in the vertical direction the maximum strain was 75% . Results in terms of horizontal and vertical stresses throughout the simulation (right-hand side of Figure 3) show that the vertical stress presents an impossible physical behavior, becoming negative at approximately 50% of vertical strain.

MODIFICATION OF A HYPERELASTIC MODEL
With Equations 33, 34 and 35 at hand, we investigate the hyperelastic model proposed by Hoss et al. (2011) concerning the fulfillment of the mathematical restrictions analyzed in this work. In their original work, they proposed the following equation for the isochoric portion of the strain energy function of incompressible hyperelastic model based on the first and second invariants of C : in which 1 C , 2 C , 3 C , 4 C , 5 C , and 6 C are parameters to be identified according to a set of experimental data.
In order to comply with the restrictions mentioned in An easy way to determine which values the parameters 1 2 6 { ; ... } C C C can assume in order to fulfill Equations 38, 39 and 40 is to impose conditions 1 0 C  , 2 0 C  , 3 5 0 Now, for the new, modified strain energy function of Equation 41, the following conditions must hold:  I  3  1 1  I  2  I II  I II We can now present the following set of consistent conditions to be satisfied by the material parameters of the model of Equation 41 in order to fulfill the Baker-Ericksen inequalities and to present a positive-definite tangent operator: Section 6.1 brings all necessary quantities for the numerical implementation of the constitutive model of Equation 41 in the framework of the Finite Element Method.

Quantities for numerical implementation
Following the notation and concepts introduced in Section 3, the consistent implementation of the model of Equation 41 in a FEM code requires the computation of the second Piola-Kirchhoff stress tensor S and the fourth-order tangent operator  .    (Yeoh, 1990) and MR5 (Equation 36) had their parameters identified to fit the Treloar's experimental data (Treloar, 1973) by means of an optimization process through which the sum of the error for uniaxial, biaxial and pure shear are minimized simultaneously. This way, we are able to define one single family of constitutive parameters that fits these three distinct deformation cases at the same time solving the problem where t C is the vector containing the parameters i t C ; T n , P n and B n are the number of experimental points for uniaxial tension, pure shear and biaxial tension, respectively; T t and E t are the numerical and experimental stresses, respectively; and T , P and B refer to the uniaxial tension, pure shear and biaxial tension tests, respectively.
The material parameters t C must satisfy the constraints of Equations 45 to 49, which makes of Equation 57 a constrained optimization problem. We chose for its solution the method of Powell (1994) named COBYLA (the source code is obtainable via the internet, see also Powell (1998)).
For the solution of the optimization procedure, we developed a FORTRAN program that calls an external finite element-based software (in our case, the general-purpose finite element software FEAP version 8.2 (Taylor, 1999)) that runs the cases of uniaxial tension, pure shear and biaxial extension consecutively and returns data of stress versus strain as outputs. These outputs are used in Equation 57 via the terms ( ) For the numerical implementation, we used the following well-known volumetric part   where  is the material's bulk modulus. Equations 11 and 18 become, respectively, Table 1 presents the parameters obtained for each model via the constrained optimization problem applied to Treloar's data together with the value of the optimized objective funcion. The bulk modulus  was chosen large enough so incompressibility could be considered. Figures 4 and 5 plot inequalities of Equations 33, 34 and 35 applied to the proposed model of Equation 41 showing that both the Baker-Ericksen inequalities and the positive-definiteness of  are satisfied. For models Y3 and MR5, the positiveness of the parameters was enforced during the optimization procedure of Equation 57, therefore assuring the satisfaction of conditions oimposed by Equations 33, 34 and 35. Figure 6 shows the numerical results of each case of homogeneous deformation compared to the experimental data. The proposed model performed better than both Y3 and MR5, which can be concluded by the analysis of the value of the minimized objective function in Table 1.    The same three models were used to numerically replicate the experimental data for an unfilled silicone rubber of Meunier et al. (2008). The optimization procedure of Equation 57 led to the material parameters presented in Table 2. Again, the proposed model results in the smaller value of the optimized objective function, meaning a more accurate fit to the experimental data (see Figure 7). All parameters of Table 2

Non-homogeneous deformation
From the results of Section 5.2.1 it can be seen that the proposed model is able to reproduce homogeneous cases of deformation with acceptable accuracy, performing better that Mooney-Rivlin's and Yeoh's model (for the experimental data considered). The present section has the goal of assessing also the capacity of the proposed model in simulating cases of non-homogeneous deformation.

Rubber bushing
Here a second type of numerical case is performed, in which we consider a situation of non-homogeneous deformation. The problem was first proposed by Zdunek and Bercovier (1986) and consists of a rubber bushing 44mm long with an outer radius of 30mm and an inner radius of 15mm (Figure 8 (left)). In the original paper, the authors used data of pure shear experiment and a 2-parameter Mooney-Rivlin model, having identified the following constitutive parameters: . It is worth noticing that these parameters fulfill Equations 33, 34 and 35.
We used the experimental data of Zdunek and Bercovier (1986) to identify the parameters of the model of Equation 41 using again the COBYLA optimization procedure, but with only the set of experimental data of pure shear. Figure 8 (right) shows the comparison between experimental and numerical results for the parameter identification and Table 3 presents the coefficients i C obtained. The bulk modulus used was identical to that of the original paper:  The inner cylinder is displaced vertically 4.8mm while the outer cylinder is fixed and results of vertical reaction force (sum of the reaction force on the inner cylinder at different levels of deformation) versus vertical displacement are obtained. The finite element mesh is shown in Figure 9 (left), in which it can be seen that 12 elements were used in the circumferential direction and 6 elements in the radial direction. Due to symmetry, only one half of the structure was modeled. Figure 9 (right) shows the comparison between the original results and those obtained in the present paper, in which a good agreement is observed, even though they are simulated with distinct hyperelastic models. It is important to point out that in this particular case we do not have experimental data of the rubber bushing in service (Figure 9 (right)) in order to check the accuracy of the proposed model versus 2-parameter Mooney-Rivlin's. However, it can be seen from the results of Figure  9 that non-homogeneous cases of deformation can also be simulated with the proposed model.

The Poynting effect
For the second and last case of non-homogeneous deformation, we capture the occurrence of the Poynting effect on a rubber band. The Poynting effect is a well-known case of multiaxial deformation consisting of combined simple shear and triaxial stretch. For elastomers, a positive Poynting effect occurs if, during a simple shear experiment, the sheared faces tend to spread apart, so a compressive normal force is needed in order to keep them parallel to each other. The model of Equation 41 with parameters of Table 1 was used to replicate the experiment of Gent et al. (2007), in which a block with an aspect ratio of 10 is subjected to a unitary amount of simple shear deformation. Among several other analyses, in the original paper, the authors perform an assessment on the sensitivity of the normal stresses for different Poisson's ratios. The goal here is not to reproduce their results, but to test if the proposed model is able to reproduce the Poyting effect qualitatively. Figure 10 depicts the expected uniform and compressive normal stress (second Piola-Kirchhoff stress) in a large portion of the center of the block. As it can be seen, a positive Poyting effect is captured by the proposed model for the material and boundary conditions applied, indicating that it can be considered for use in complex finite element simulations.

CONCLUSIONS
A hyperelastic model providing simultaneously numerical stability and consistent physical behavior is proposed for rubber-like materials in the context of FEM. This was accomplished by a mathematical modification in the strain energy function of the model proposed by Hoss et al. (2011), as it was detected that their original model was not able to fulfill both requirements at the same time. The results presented in Section 4.3 are new in the literature, and show that the lack of satisfaction of B-E inequalities may lead to results that are not physically possible.
It is important to keep in mind, however, that depending on the parameters obtained during the calibration of the model, one or both of the restrictions might be violated, leading to either inconsistent results or lack of quadratic convergence during FE calculations, or both. For that reason, it is suggested for the analyst to make use of constrained optimization procedures during parameter identification, taking into account the inequalities given by Equations (45) to (49).
The complete constitutive framework for the numerical implementation of the proposed model presented in Equation 41 in FEM code is derived and presented in Equations 9 to 13, 50, 51, and 59 for the second Piola-Kirchhoff stress tensor and in Equations 17, 19 to 22, 52 to 55 and 60 for the fourth-order tangent operator. In the present work, the implementation was performed in the general-purpose FEM code FEAP and the simulations depicted in Section 5.2 demonstrate that the model is capable of accurately modeling the mechanical behavior of isotropic rubber-like materials both in homogeneous and non-homogeneous cases of deformation, even showing better performance than the wellknown models of Rivlin and Saunders (1951) and Yeoh (1990) for the homogeneous cases of deformation analyzed.
As suggestions for future work, we mention the possibility to additively couple the proposed model to an anisotropic function, allowing the simulation of fiber reinforced rubber-like materials. Moreover, it is possible to study the inclusion of degradation, damage or account for inelastic strains in order to account for energy dissipation phenomena, such as viscoelasticity. Combining these two features (anisotropy and viscoelasticity) opens room for the simulation of biological soft tissues.