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

Abstract

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.

Keywords:
Hyperelasticity; Baker-Ericksen inequalities; Constitutive framework; Finite Element Method

Graphical Abstract

1 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 non-homogeneous 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.

2 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 (Holzapfel et al., 2000; 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).

As examples of recent works proposing new hyperelastic models, we can cite Bien-aimé et al. (2020), Külcü (2020), Yaya and Bechir (2018), and Korba and Barkey (2017). Significant reviews and comparison between different models can be found in Mihai and Goriely (2017), Steinmann et al. (2012), Bien-aimé et al. (2019) and the references therein.

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.

3 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 $\Omega \subset {ℝ}^3$ subjected to a finite deformation $\phi :\Omega \to {ℝ}^3$ one defines the deformation gradient $F=\partial \phi (X)/\partial X$ with the Jacobian $J=\det F>0$, in which $X$ is a material point in the reference, undeformed configuration. The deformation gradient $F$ can be multiplicatively decomposed into (Flory, 1961)

$F=\widehat{F}\overline{F},$ (1)

in which $\widehat{F}$ refers to the volumetric part and $\overline{F}$ to the isochoric part, respectively defined as

$\widehat{F}=J^{1/3}\mathrm{1,}with\det \widehat{F}=\det F=Jand$ (2)

$\overline{F}=J^{-1/3}F,with\det \overline{F}=1.$ (3)

From Equation 1, the isochoric right Cauchy-Green tensor is defined as

$\overline{C}={\overline{F}}^T\overline{F}={\overline{U}}^2=J^{-2/3}C$ (4)

with $\det \overline{C}=1$, $J^2=\det C$ and $C=F^{T}F$. The main invariants of $\overline{C}$ are

${\text{I}}_{\overline{C}}=\text{tr}\overline{C},{\text{II}}_{\overline{C}}=\frac{1}{2}\left({(\text{tr}\overline{C})}^2-\text{tr}{\overline{C}}^2\right)and{\text{III}}_{\overline{C}}=\det \overline{C}=1.$ (5)

Now, the strain energy function $W$ of a compressible, isotropic hyperelastic model can be additively decoupled into isochoric $\overline{W}$ and volumetric $\widehat{W}$ parts depending on the strain invariants of $\overline{C}$ as

$W=\overline{W}(I_{\overline{C}},I{I}_{\overline{C}})+\widehat{W}(J).$ (6)

It is important to point out that it is also possible to write $W$ in terms of the invariants $I_{\overline{B}}$ and $I{I}_{\overline{B}}$ of the left Cauchy-Green deformation tensor $\overline{B}$ (defined by $\overline{B}=\overline{F}{\overline{F}}^T$) and, although $C$ and $B$ are different tensors (and, consequently, $\overline{C}$ and $\overline{B}$ too), their invariants are identical, which leads to the fact that the equation for $W$ written in terms of both $\overline{C}$ or $\overline{B}$ have the same format.

The second Piola-Kirchhoff stress tensor $S$ is written in terms of the Cauchy stress tensor $\sigma$ as

$S=\left(\det F\right)F^{-1}\sigma F^{-T},$ (7)

in which

$\sigma =2J^{-1}F\frac{\partial W(C)}{\partial C}{F}^T=2J^{-1}F\frac{\partial W(B)}{\partial B}{F}^T.$ (8)

Following the same decomposition criterion of the strain energy function, the second Piola-Kirchhoff stress tensor can be additively decoupled as (Holzapfel, 2000)

$S=2\frac{\partial W}{\partial C}=\overline{S}+\widehat{S},$ (9)

where $\overline{S}$ is the isochoric portion of $S$ and $\widehat{S}$ is the volumetric portion, computed, respectively, by

$\overline{S}=J^{-2/3}\mathbb{P}:\check{S},$ (10)

$\widehat{S}=J\frac{\partial \widehat{W}}{\partial J}{C}^{-1},$ (11)

in which $\mathbb{P}$ is a fourth-order projection tensor and $\check{S}$ is a fictitious stress tensor $\check{S}=\partial \overline{W}(\overline{C})/\partial \overline{C}$, computed as

$\mathbb{P}=\mathbb{I}-\frac{1}{3}{C}^{-1}\otimes C,$ (12)

$\check{S}=\frac{\partial \overline{W}}{\partial \overline{C}}=\overline{{\gamma }_1}1+\overline{{\gamma }_2}\overline{C},$ (13)

respectively, in which

$\overline{{\gamma }_1}=2\left(\frac{\partial \overline{W}(I_{\overline{C}},I{I}_{\overline{C}})}{\partial I_{\overline{C}}}+I_{\overline{C}}\frac{\partial \overline{W}(I_{\overline{C}},I{I}_{\overline{C}})}{\partial I{I}_{\overline{C}}}\right),$ (14)

$\overline{{\gamma }_2}=-2\frac{\partial \overline{W}(I_{\overline{C}},I{I}_{\overline{C}})}{\partial I{I}_{\overline{C}}},$ (15)

${\mathbb{I}}_{ikjl}={\delta }_{ik}{\delta }_{jl}.$ (16)

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):

$ℂ=2\frac{\partial S}{\partial C}=\overline{ℂ}+\widehat{ℂ},$ (17)

in which the volumetric portion $\widehat{ℂ}$ is computed as

$\widehat{ℂ}=J\left(\frac{\partial \widehat{W}}{\partial J}+J\frac{{\partial }^2\widehat{W}}{\partial J^2}\right)\left(C^{-1}\otimes C^{-1}\right)-2J\frac{\partial \widehat{W}}{\partial J}(C^{-1}\odot C^{-1}),$ (18)

with the notation ${(C^{-1}\odot C^{-1})}_{ijkl}=\frac{1}{2}(C_{ik}^{-1}{C}_{jl}^{-1}+C_{il}^{-1}{C}_{jk}^{-1})$.

The computation of the isochoric part $\overline{ℂ}$ of Equation 17 is as follows:

$\overline{ℂ}=\mathbb{P}:\check{ℂ}:{\mathbb{P}}^T+\frac{2}{3}Tr(J^{-2/3}\check{S})\tilde{\mathbb{P}}-\frac{2}{3}(C^{-1}\otimes \overline{S}+\overline{S}\otimes C^{-1}),$ (19)

with the definitions

$\tilde{\mathbb{P}}=C^{-1}\odot C^{-1}-\frac{1}{3}{C}^{-1}\otimes C^{-1},$ (20)

$Tr(*)=(*):C,$ (21)

$\check{ℂ}=J^{-4/3}({\overline{\delta }}_{1}1\otimes 1+{\overline{\delta }}_2(1\otimes \overline{C}+\overline{C}\otimes 1)+{\overline{\delta }}_3\overline{C}\otimes \overline{C}+{\overline{\delta }}_4\mathbb{I}),$ (22)

in which ${\overline{\delta }}_1$, ${\overline{\delta }}_2$, ${\overline{\delta }}_3$ and ${\overline{\delta }}_4$ are defined as

${\overline{\delta }}_1=4\left(\frac{{\partial }^{2}W}{\partial I_{\overline{C}}^2}+2I_{\overline{C}}\frac{{\partial }^{2}W}{\partial I_{\overline{C}}\partial I{I}_{\overline{C}}}+\frac{\partial W}{\partial I{I}_{\overline{C}}}+I_{\overline{C}}^2\frac{{\partial }^{2}W}{\partial I{I}_{\overline{C}}^2}\right),$ (23)

${\overline{\delta }}_2=-4\left(\frac{{\partial }^{2}W}{\partial I_{\overline{C}}\partial I{I}_{\overline{C}}}+I_{\overline{C}}\frac{{\partial }^{2}W}{\partial I{I}_{\overline{C}}^2}\right),$ (24)

${\overline{\delta }}_3=4\frac{{\partial }^{2}W}{\partial I{I}_{\overline{C}}^2},$ (25)

${\overline{\delta }}_4=-4\frac{\partial W}{\partial I{I}_{\overline{C}}}.$ (26)

4 MATHEMATICAL RESTRICTIONS TO THE STRAIN ENERGY FUNCTIONS

4.1 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

$\sigma =a_{0}I+a_{1}B+a_2{B}^2,$ (27)

in which the coefficients $a_i$ 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):

$\sigma =b_{0}I+b_{1}B+b_2{B}^{-1},$ (28)

in which, assuming incompressibility of the material, $b_0=-p$, $b_1=2\partial W/\partial I_{\overline{C}}$ and $b_2=-2\partial W/\partial I{I}_{\overline{C}}$ (with $p$ being the hydrostatic pressure). Equation 28 can be written as (Atkin and Fox, 1980; Bonnet and Wood, 1997)

$\sigma =-pI+2\frac{\partial W}{\partial I_{\overline{C}}}B-2\frac{\partial W}{\partial I{I}_{\overline{C}}}{B}^{-1}.$ (29)

Recalling that the B-E inequalities guarantee that for an isotropic elastic body the maximum principal stress ${\sigma }_i$ always occurs in the same direction of the maximum principal stretch ${\lambda }_i$, it means mathematically that

$\left({\sigma }_i-{\sigma }_j\right)\left({\lambda }_i-{\lambda }_j\right)>0for{\lambda }_i\ne {\lambda }_j$ (30)

which can be written in terms of the coefficients $b_i$ of Equation 28 for isotropic incompressible materials as

$b_1>{\lambda }_j{}^2{b}_{2}for{\lambda }_i\ne {\lambda }_k,$ (31)

$b_1\ge {\lambda }_j{}^2{b}_{2}for{\lambda }_i={\lambda }_k,$ (32)

in which $i$ and $j$ and $k$ are permutations of $1$, $2$ and $3$. Taking into account that $b_1=2\partial W/\partial I_{\overline{C}}$ and $b_2=-2\partial W/\partial I{I}_{\overline{C}}$, Truesdell and Noll (1965) conclude by inspection that in order to satisfy the B-E inequalities (Equation 30), the following conditions must hold:

$b_1>0or\frac{\partial W}{\partial I_{\overline{C}}}>\mathrm{0,}$ (33)

$b_2\le 0or\frac{\partial W}{\partial I{I}_{\overline{C}}}\ge \mathrm{0,}$ (34)

4.2 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_{\overline{C}}$, $I{I}_{\overline{C}}$ and $II{I}_{\overline{C}}$ its Hessian is computed according to

$\left[\begin{array}{cc}\frac{{\partial }^{2}W}{\partial I_{\overline{C}}^2}& \frac{{\partial }^{2}W}{\partial I_{\overline{C}}\partial I{I}_{\overline{C}}}\\ \frac{{\partial }^{2}W}{\partial I_{\overline{C}}\partial I{I}_{\overline{C}}}& \frac{{\partial }^{2}W}{\partial I{I}_{\overline{C}}^2}\end{array}\right]$

whose positive defineteness is assured if its determinant is positive, or

$\left(\frac{{\partial }^{2}W}{\partial I_{\overline{C}}^2}\right)\left(\frac{{\partial }^{2}W}{\partial I{I}_{\overline{C}}^2}\right)-{\left(\frac{{\partial }^{2}W}{\partial I_{\overline{C}}\partial I{I}_{\overline{C}}}\right)}^2>0.$ (35)

4.3 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) (MR5), whose constitutive parameters don't satisfy the B-E inequality of Equation 34, in order to show what kind of problem can occur in the results of FE simulations if these restrictions are not satisfied.

Making use of the five-terms truncation of a Rivlin series

$W=C_1\left(I_C-3\right)+C_2\left(I{I}_C-3\right)+C_3\left(I_C-3\right)\left(I{I}_C-3\right)+C_4{\left(I_C-3\right)}^2+C_5{\left(I{I}_C-3\right)}^2,$ (36)

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 were $C_1=1.787$, $C_2=-1.013$, $C_3=-0.0047$, $C_4=0.286$, and $C_5=0.000237$. In the diagram of Figure 2, in which Equation 34 is plotted, it can be seen that $\partial W/\partial I{I}_{\overline{C}}$ assumes only negative values in the ranges of $I_{\overline{C}}$ and $I{I}_{\overline{C}}$ considered.

Figure 1
Theoretical fitting of model MR5 to experimental data of biaxial extension presented in Meier et al. (2003).

Figure 2
Violation of the B-E inequality of Equation 34 by the MR5 model.

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.

Figure 3
Violation of the B-E inequalities: (left) mesh and boundary conditions and (right) results in terms of stresses through time.

5 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 $\overline{C}$:

$\overline{W}=\frac{C_1}{C_2}\left(1-e^{-C_2(I_{\overline{C}}-3)}\right)+\frac{C_5}{2C_3}\left[{\left(1+\frac{C_3(I_{\overline{C}}-3)}{C_4}\right)}^{C_4}-1\right]+C_6\ln \frac{I{I}_{\overline{C}}}{3},$ (37)

in which $C_1$, $C_2$, $C_3$, $C_4$, $C_5$, and $C_6$ are parameters to be identified according to a set of experimental data.

In order to comply with the restrictions mentioned in Sections 4.1 and 4.2 we apply Equations 33, 34 and 35 to Equation 37 and arrive, respectively, at

$\frac{\partial \overline{W}}{\partial I_{\overline{C}}}=C_1{e}^{-C_2\left(I_{\overline{C}}-3\right)}+\frac{C_5}{2}{\left(1+\frac{C_3\left(I_{\overline{C}}-3\right)}{C_4}\right)}^{C_4-1}>0,$ (38)

$\frac{\partial \overline{W}}{\partial I{I}_{\overline{C}}}=\frac{C_6}{I{I}_{\overline{C}}}\ge 0,$ (39)

$\left(\frac{{\partial }^{2}W}{\partial I_{\overline{C}}^2}\right)\left(\frac{{\partial }^{2}W}{\partial I{I}_{\overline{C}}^2}\right)-{\left(\frac{{\partial }^{2}W}{\partial I_{\overline{C}}\partial I{I}_{\overline{C}}}\right)}^2=\frac{-C_6}{I{I}_{\overline{C}}^2}\left(-C_1{C}_2{e}^{-C_2(I_{\overline{C}}-3)}+\frac{C_5{C}_3}{2C_4}\left(C_4-1\right){\left(1+\frac{C_3\left(I_{\overline{C}}-3\right)}{C_4}\right)}^{C_4-2}\right)>0.$ (40)

An easy way to determine which values the parameters $\{C_1;C_2\mathrm{...}{C}_6\}$ can assume in order to fulfill Equations 38, 39 and 40 is to impose conditions $C_1>0$, $C_2<0$, $C_3{C}_5>0$ and $C_4>2$. Regarding$C_6$, however, we face an inconsistency. Having in mind that both $I_{\overline{C}}\ge 3$ and $I{I}_{\overline{C}}\ge 3$, a negative value of $C_6$ violates Equation 39 (plausible mechanical behavior could not be guaranteed), but a positive value of $C_6$ violates Equation 40 (numerical stability could not be assured).

With that in mind, and by inspection, in this work we propose a modification to the hyperelastic model of Equation 37 in order to overcome this problem by multiplying its last term by $I{I}_{\overline{C}}$, arriving at

$\overline{W}=\frac{C_1}{C_2}\left(1-e^{-C_2(I_{\overline{C}}-3)}\right)+\frac{C_5}{2C_3}\left[{\left(1+\frac{C_3(I_{\overline{C}}-3)}{C_4}\right)}^{C_4}-1\right]+C_{6}I{I}_{\overline{C}}\ln \frac{I{I}_{\overline{C}}}{3}$ (41)

Now, for the new, modified strain energy function of Equation 41, the following conditions must hold:

$\frac{\partial \overline{W}}{\partial I_{\overline{C}}}=C_1{e}^{-C_2\left(I_{\overline{C}}-3\right)}+\frac{C_5}{2}{\left(1+\frac{C_3\left(I_{\overline{C}}-3\right)}{C_4}\right)}^{C_4-1}>0,$ (42)

$\frac{\partial \overline{W}}{\partial I{I}_{\overline{C}}}=C_6\ln \left(\frac{I{I}_{\overline{C}}}{3}\right)+C_6\ge 0,$ (43)

$\left(\frac{{\partial }^{2}W}{\partial I_{\overline{C}}^2}\right)\left(\frac{{\partial }^{2}W}{\partial I{I}_{\overline{C}}^2}\right)-{\left(\frac{{\partial }^{2}W}{\partial I_{\overline{C}}\partial I{I}_{\overline{C}}}\right)}^2=\frac{C_6}{I{I}_{\overline{C}}}\left(-C_1{C}_2{e}^{-C_2(I_{\overline{C}}-3)}+\frac{C_5{C}_3}{2C_4}\left(C_4-1\right){\left(1+\frac{C_3\left(I_{\overline{C}}-3\right)}{C_4}\right)}^{C_4-2}\right)>0$ (44)

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:

$C_1>0,$ (45)

$C_2<0,$ (46)

$C_3{C}_5>0,$ (47)

$C_4>2,$ (48)

$C_6>0.$ (49)

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.

5.1 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 $ℂ$.

The isochoric part $\overline{S}$ of the second Piola-Kirchhoff stress tensor is computed according to Equations 10 to 15, in which $\overline{{\gamma }_1}$ and $\overline{{\gamma }_2}$, applied to Equation 41 are computed, respectively, as

$\overline{{\gamma }_1}=2\left(C_1{e}^{-C_2(I_{\overline{C}}-3)}+\frac{C_5}{2}{\left(1+\frac{C_3}{C_4}(I_{\overline{C}}-3)\right)}^{C_4-1}+I_{\overline{C}}\left(C_6+C_6\ln \frac{I{I}_{\overline{C}}}{3}\right)\right),$ (50)

$\overline{{\gamma }_2}=-2\left(C_6+C_6\ln \frac{I{I}_{\overline{C}}}{3}\right).$ (51)

Equations 19 to 26 are used to write the isochoric part $\overline{ℂ}$ of the tangent operator, in which the quantities $\overline{{\delta }_1}$, $\overline{{\delta }_2}$, $\overline{{\delta }_3}$, and $\overline{{\delta }_4}$ are computed, respectively, as

$\overline{{\delta }_1}=4\left(-C_1{C}_2{e}^{-C_2\left(I_{\overline{C}}-3\right)}+\frac{C_5{C}_3\left(C_4-1\right)}{2C_4}{\left(1+\frac{C_3\left(I_{\overline{C}}-3\right)}{C_4}\right)}^{C_4-2}+C_6+C_6\ln \frac{I{I}_{\overline{C}}}{3}+\frac{I_{\overline{C}}^2{C}_6}{I{I}_{\overline{C}}}\right),$ (52)

$\overline{{\delta }_2}=-4I_{\overline{C}}\frac{C_6}{I{I}_{\overline{C}}},$ (53)

$\overline{{\delta }_3}=4\frac{C_6}{I{I}_{\overline{C}}},$ (54)

$\overline{{\delta }_4}=-4\left(C_6+C_6\ln \frac{I{I}_{\overline{C}}}{3}\right).$ (55)

5.2 Representative examples

5.2.1 Homogeneous deformation

5.2.1.1 Experimental data of Treloar (1973)

The first set of numerical examples reconstructed three different cases of homogeneous deformation: uniaxial tension, pure shear and biaxial extension. The proposed model of Equation 41 together with Yeoh's model (Y3) (Yeoh, 1990)

$W=C_1\left(I_C-3\right)+C_2{\left(I_C-3\right)}^2+C_3{\left(I_C-3\right)}^3$ (56)

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

$\underset{C_t}{\min }\left[{\displaystyle \sum _{i=1}^{n_T}{(t_{T_i}-t_{E_i})}_T^2}+{\displaystyle \sum _{j=1}^{n_P}{(t_{T_j}-t_{E_j})}_P^2}+{\displaystyle \sum _{k=1}^{n_B}{(t_{T_k}-t_{E_k})}_B^2}\right],$ (57)

where $C_t$ is the vector containing the parameters $C_{t_i}$; $n_T$, $n_P$ and $n_B$ are the number of experimental points for uniaxial tension, pure shear and biaxial tension, respectively; $t_T$ and $t_E$ 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 $C_t$ 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 ${(t_{T_i})}_T$, ${(t_{T_i})}_P$ and ${(t_{T_i})}_B$ together with the experimental data ${(t_{E_i})}_T$, ${(t_{E_i})}_P$ and ${(t_{E_i})}_B$ to construct the objective function. The program then calls COBYLA, which computes a new set of material parameters $C_t$ that are used to re-enter FEAP and re-calculate the new theoretical stresses and strain. This procedure is repeated until the difference of the newly computed $C_t$ and the previous set reaches a prescribed tolerance.

For the numerical implementation, we used the following well-known volumetric part $\widehat{W}\left(J\right)$

$\widehat{W}\left(J\right)=\frac{1}{2}\kappa {\left(J-1\right)}^2$ (58)

where $\kappa$ is the material's bulk modulus. Equations 11 and 18 become, respectively,

$\widehat{S}=J\kappa \left(J-1\right)C^{-1},$ (59)

$\widehat{ℂ}=J\left(2\kappa J-\kappa \right)\left(C^{-1}\otimes C^{-1}\right)-2\kappa J\left(J-1\right)\left(C^{-1}\odot C^{-1}\right).$ (60)

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 $\kappa$ 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.

Table 1
Hyperelastic parameters fitted to Treloar's experimental data.

	Proposed model	Y3	MR5
$C_1$	1.20E-01MPa	1.59E-01 MPa	1.42E-01 MPa
$C_2$	-6.80E-06	1.00E-06 MPa	2.10E-03 MPa
$C_3$	1.30E-01	2.51E-05 MPa	1.00E-06 MPa
$C_4$	3.00	-	1.00E-06 MPa
$C_5$	4.50E-02 MPa	-	4.72E-05 MPa
$C_6$	1.65E-04	-	-
$\kappa$	1300 MPa	1300 MPa	1300 MPa
Objective function	1.04	1.42	40.85

Figure 4
Satisfaction of the proposed model to (left) Equation 33 and (right) Equation 34.

Figure 5
Satisfaction of the proposed model to Equation 35.

Figure 6
FE simulations of the homogeneous deformation cases using the three models to replicate Treloar's experimental data for (t) uniaxial tension, (b) biaxial extension and (p) pure shear.

5.2.1.2 Experimental data of Meunier et al. (2008)

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 fulfill Equations 33, 34 and 35, since they were enforced during the constrained optimization algorithm.

Table 2
Hyperelastic parameters fitted to the experimental data of Meunier et al. (2008)

	Proposed model	Y3	MR5
$C_1$	2.00E-02MPa	1.55E-01 MPa	1.64E-01 MPa
$C_2$	-3.80E-01	5.60E-03 MPa	2.47E-04 MPa
$C_3$	3.60E-03	1.10E-03 MPa	1.00E-06 MPa
$C_4$	2.60	-	1.00E-06 MPa
$C_5$	2.85E-02 MPa	-	1.69E-03 MPa
$C_6$	1.00E-03	-	-
$\kappa$	1300 MPa	1300 MPa	1300 MPa
Objective function	2.67E-02	5.70E-02	2.80E-02

Figure 7
FE simulations of the homogeneous deformation cases using the three models to replicate the experimental data ofMeunier et al. (2008) for (t) uniaxial tension, (b) biaxial extension and (p) pure shear.

5.2.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.

5.2.2.1 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: $C_{10}=0.177$, $C_{01}=0.045$ and bulk modulus $\kappa =666N/m{m}^2$. 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 $C_i$ obtained. The bulk modulus used was identical to that of the original paper: $\kappa =666N/m{m}^2$ .

Figure 8
Rubber bushing: (left) cross-section of the structure and (right) results of the parameter identification procedure.

Table 3
Hyperelastic parameters of the proposed model calibrated to pure shear experimental data.

$C_1$	$C_2$	$C_3$	$C_4$	$C_5$	$C_6$	$\kappa$
1.66E-01 MPa	-1.77E-03	4.70E-03 MPa	2.00	8.49E-02 MPa	1.00E-03	666 MPa

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.

Figure 9
Rubber bushing: (left) finite element mesh and (right) results.

5.2.2.2 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.

Figure 10
Compressive normal stress (second Piola-Kirchhoff stress) denoting a positive Poynting effect at a block subjected to simple shear.

6 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 well-known 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.

Acknowledgements

F. T. Stumpf acknowledges CNPq/Brazil through project 407375/2018-0 for the financial support during the development of this research. R. J. Marczak acknowledges CNPq/Brazil through project 310649/2017-0 for the financial support during the development of this research.

References

• Arora N, Li J, Slesarenko V, Rudykh S. (2020) Microscopic and long-wave instabilities in 3D fiber composites with non-Gaussian hyperelastic phases. International Journal of Engineering Science; 157, 103408
• Arruda E, Boyce M. (1993) A Three-dimensional Constitutive Model for the Large Stretch Behavior of Rubber Elastic Materials. J of the Mech Phys Solids; 41: 389-412.
• Atkin R, Fox N. (1980) An Introduction to the Theory of Elasticity. New York: Longman.
• Baker M, Ericksen J. (1954) Inequalities Restricting the Form of the Stress-Deformation Relation for Isotropic Elastic Solids and Reiner-Rivlin Fluids. Journal of the Washington Academy of Science; 44: 33-35.
• Balzani D, Neff P, Schroeder J, Holzapfel G. (2006) A Polyconvex Framework for Soft Biological Tissues. Adjustment to Experimental Data. Int J Solids and Structures; 43: 6059-6070.
• Beatty M. (2003) An Average-Stretch Full-Network Model for Rubber Elasticity. Journal of Elasticity; 70: 65-86.
• Bien-aimé LKM, Blaise BB, Beda T. (2019) Comparison of Continuum Constitutive Hyperelastic Models based on Exponential Forms. International Journal of Innovative Science and Research Technology; 4(7): 1360-1367
• Bien-aimé LKM, Blaise BB, Beda T. (2020) Characterization of hyperelastic deformation behavior of rubberlike materials. SN Applied Sciences; 2, p.648
• Bonnet J, Wood R. (1997) Nonlinear Continuum Mechanics for Finite Element Analysis. New York: Cambridge University Press.
• Boyce M, Arruda E. (2000) Constitutive Models of Rubber Elasticity: A Review. Rubber Chemistry and Technology; 73(3):504-523.
• Flory P. (1961) Thermodynamic relations for high elastic materials. Trans Faraday Soc; 57: 829-838.
• Fried E. (2002) An Elementary Molecular-Statistical Basis for the Mooney and Rivlin-Saunders Theories of Rubber Elasticity. J Mechanics and Physics of Solids; 50: 571-582.
• Fung Y. (1967) Elasticity of soft tissues in sample elongation. American Journal of Physiology; 213: 1532-1540.
• Gent A, Suh J, Kelly III S. (2007) Mechanics of rubber shear springs. International Journal of Nonlinear Mechanics; 42:241-249.
• Gent A. (1996) A New Constitutive Relation for Rubber. Rubber Chemistry and Technology; 69(1): 59-61.
• Ghiba I, Neff P, Silhavy M. (2015) The Exponentiated Hencky-Logarithmic Strain Energy. Improvement of Planar Polyconvexity. Int J Non-Linear Mechanics; 71: 48-51.
• Holzapfel G, Gasser T, Ogden R. (2000) A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models. Journal of Elasticity; 61: 1-48.
• Holzapfel G. (2000) Nonlinear Solid Mechanics: a Continuum Approach for Engineering. John Wiley and Sons.
• Horgan C, Saccomandi G. (2006) Phenomenological Hyperelastic Strain-Stiffening Constitutive Models for Rubber. Rubber Chemistry and Technology; 79: 152-169.
• Hoss L, Stumpf F, Bortolli D, Marczak R. (2011) Constitutive Models for Rubber-Like Materials: Comparative Analysis, Goodness of Fitting Measures, and Proposal of a New Model. In: International Rubber Conference; Sao Paulo.
• Kilian H. (1981) Equation of state of real networks. Polymer; 22: 209-217.
• Knowles J. (1977) The Finite Anti-Plane Shear Field Near the Tip of a Crack for a Class of Incompressible Elastic Solids. Int J Fract. 13: 611-639.
• Korba AG, Barkey ME. (2017) New Model for Hyper-Elastic Materials Behavior With an Application on Natural Rubber. ASME 2017 12th International Manufacturing Science and Engineering Conference collocated with the JSME/ASME 2017 6th International Conference on Materials and Processing
• Külcü ID. (2020) A hyperelastic constitutive model for rubber-like materials. Archive of Applied Mechanics; 90: 615-622
• Li T, Yang Q-S, Ma L-H, Zhang YX. (2019) A non-Gaussian model for the chemo-mechanical coupling behavior of largely deformed hydrogels. Soft Materials; 17(2): 126-136
• Mansouri M, Darijani H. (2014) Constitutive Modeling of Isotropic Hyperelastic Materials in an Exponential Framework Using a Self-Contained Approach. Int J Solids and Structures; 51: 4316-4326.
• Marsden J, Hughes T. (2012) Mathematical Foundations of Elasticity. Prentice-Hall Personal Computing SeriesPrentice-Hall.
• Meier P, Khader S, Preuss R, Dietrich J, Voges D. (2003) Uniaxial and equi-biaxial tension tests of silicone elastomer. In: European Conference on Constitutive Models for Rubber.; September 15-17; London.
• Meunier L, Chagnon G, Favier D, Orgéas L, Vacher P. (2008) Mechanical experimental characterisation and numerical modeling of an unfilled silicone rubber. Polymer Testing; 27: 765-777.
• Mihai L, Goriely A. (2017) How to characterize a nonlinear elastic material: A review on nonlinear constitutive parameters in isotropic finite elasticity. Proceedings of the Royal Society A; 473(20170607).
• Mooney M. (1940) A theory of large elastic deformation. Journal of Applied Physics; 11: 582-592.
• Morch A, Astruc L, Witz J-F, Lesaffre F, Lecomte-Grosbras P, Soulat D, Brieu M. (2019) Modeling of anisotropic hyperelastic heterogeneous knitted fabric reinforced composites. Journal of the Mechanics and Physics of Solids; 127: 47-61
• Neff P, Ghiba ID, Lankeit J. (2015) The exponentiated Hencky-logarithmic strain energy. Part i: constitutive issues and rank-one convexity. Journal of Elasticity; 121(2): 143-234.
• Ogden RW. (1984) Nonlinear Elastic Deformations. John Wiley & Sons, New York.
• Powell M. (1994) A direct search optimization method that models the objective and constraint functions by linear interpolations. In: Gomez S, Hennart J., eds. Advances in Optimization and Numerical AnalysisDordrecht: Kluwer Academic. (pp.51-67).
• Powell M. (1998) Direct search algorithms for optimization calculations. Acta Numerica; 7: 287-336.
• Rivlin R, Saunders D. (1951) Large Elastic Deformations of Isotropic Materials. VII. Experiments on the Deformation of Rubber. Philos Trans R Soc Lon; 243(A): 251-288.
• Steinmann P, Hossain M, Possart G. (2012) Hyperelastic models for rubber-like materials: consistent tangent operators and suitability for Treloar’s data. Archive of Applied Mechanics; 82: 1183-1217.
• Taylor R. (1999) FEAP - A Finite Element Analysis Program, Version 7.1. University of California at Berkeley.
• Treloar L. (1973) The Elasticity and Related Properties of Rubbers. Reports on Progress in Physics; 36: 755-829.
• Treloar L. (1975) The Physics of Rubber Elasticity. Oxford: Clarendon Press .
• Truesdell C, Noll W. (1965) The Nonlinear Field Theories of Mechanics. Berlin/Heidelberg: Springer.
• Wang M, Guth E. (1952) Statistical Theory of Networks of Non - Gaussian Flexible Chains. J Chem Phys; 20: 1144-1157.
• Wrubleski E, Marczak R. (2018) A new pseudo-energy function to simulate the Mullins effect. Journal of Elastomers and Plastics; 50(6): 554-575.
• Yaya K, Bechir H. (2018) A new hyper-elastic model for predicting multi-axial behaviour of rubberlike materials: formulation and computational aspects. Mechanics of Time-Dependent Materials; 22: 167-186
• Yeoh O, Fleming P. (1997) A New Attempt to Reconcile the Statistical Phenomenological Theories of Rubber Elasticity. Journal of Polymer Science - Part B: Polymer Physics; 35: 1919-1931.14
• Yeoh O. (1990) Characterization of Elastic Properties of Carbon Black Filled Rubber Vulcanizates. Rubber Chemistry and Technology; 63: 792-805.
• Yeoh O. (1993) Some Forms of the Strain Energy Function for Rubbers. Rubber Chemistry and Technology; 66(5): 754-771.
• Zdunek A, Bercovier M. (1986) Numerical Evaluation of Finite Element Methods for Rubber Parts. SAE Transactions; 95(4):257-269.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  15 Mar 2021
• Date of issue
  2021

History

• Received
  02 Dec 2020
• Reviewed
  17 Dec 2020
• Accepted
  22 Dec 2020
• Published
  08 Jan 2021