Hyperelasticity with rate-independent microsphere hysteresis model for rubberlike materials

The mechanical behavior of elastomers strongly differs from one to another. Among these differences, hysteresis upon cyclic load can take place, and can be either rate-dependent or rate-independent. In the present paper, a microsphere model taking into account rate-independent hysteresis is proposed and applied to model ﬁlled silicone rubbers behavior. The hysteresis model is based on a combination of monodimensional constitutive equations distributed in space. The behavior of each direction is described by a collection of parallel spring slider elements. The sliders are Coulomb dampers with non-zero break-free force in tension. This model is tested on a ﬁlled silicone rubber by the way of uniaxial tensile and pure shear tests. The mechanical response of the material is well predicted for such tests. Finally, the constitutive equations are implemented in the ﬁnite element software ABAQUS. Calculation results highlight good performances of the proposed model.


Introduction
The study and modeling of the mechanical behavior of rubberlike materials have been widely studied in last decades due to the increasing number of industrial applications, such as vibration isolators, tires or shock absorbers, non exhaustively.Classically, rubbers are filled by mineral fillers in order to improve their physical properties.The addition of fillers typically implies an increase of stiffness and a reinforcement of crack growth resistance [1,2].However, it also induces numerous additional effects [1,2].Among them, one can cite the stress softening, which mainly occurs between the first and second loads, and often called the Mullins effect [3] (which can rarely be observed in unfilled rubbers too [4]), the stress relaxation and the mechanical hysteresis (unload different from load).
The load and unload responses of filled rubber differ during cyclic tests.Even if this evolution is mainly due to the Mullins effect during the first cycle, a difference between load and unload responses is still observed once the material is softened, i.e. after the first cycle.This phenomenon, so-called hysteresis, can depend on the strain rate [5], the crosslinks density [6] or the temperature [7].
Several micromechanism-inspired approaches of the hysteresis phenomenon were proposed in the literature, based on different physical considerations.To the opinion of Bergström and Boyce [6], hysteresis is induced by the viscous reptation of elastically inactive macromolecules.This type of micromechanism has also motivated a micro-sphere approach [8], where hysteresis is considered to be due to the fillers entanglement.Some authors considered that a part of the broken cross-links can be re-formed upon unload [9].More recently, the idea of a successive breakdown of filler clusters taking place during a load, and of a complete re-aggregation of the filler particles taking place during an unload was proposed [10].This implies filler-induced hysteresis.
According to either the physical or the phenomenological approach used, different viscoelastic models with or without damage were proposed to take into account the hysteresis (for example [11][12][13]).A classical hypothesis of the proposed constitutive models is the multiplicative split of the deformation gradient into elastic and inelastic parts [6,[14][15][16][17].Other approaches used a history dependent function to account for the hysteresis [18].
Even though most of the papers focused on the viscoelastic behavior of the material, some of them studied the rate-independent hysteresis [9,11,18,19].To distinguish stress-softening and hysteresis, most of the authors removed the Mullins effect from the material behavior by carrying out beforehand several loading cycles [6,15].This paper focuses on the quasi-static modeling (i.e. the rate independent behavior) of hysteresis effect after stress-softening.The aim of this work is to propose a directional method able to represent the difference between load and unload, easily extendable to anisotropy and well-adapted to the implementation in a finite element code.
The different mechanical tests carried out on a filled silicone rubber and the results obtained are presented in Section 2. The theoritical background to model hysteresis is introduced in Section 3. In Section 4, experimental data are compared with simulations and results are discussed.Section 5 presents model implementation in a finite element code with calculation results obtained in case of heterogeneous test and capacities of the model.Finally, concluding remarks close the paper.

Materials
The material considered here is a filled silicone rubber (Bluestar RTV 3428).Its mechanical properties were previously investigated [7,20,21] and it was shown that this material exhibits stress-softening, occurring mainly during the first load, hysteresis and temperature dependent behavior.

Loading conditions
Two types of classical tests are carried out, uniaxial tensile tests and pure shear tests.These tests are performed using a Gabo Eplexor 500 N.For the uniaxial tensile tests, the length, width and thickness of the specimens are 12, 2 and 2 mm, respectively and a 25 N load cell is used.For the pure shear tests, the length, width and thickness of the specimens are 2, 40 and 2 mm, respectively and a 500 N load cell is used.Several tests were carried out at different strain rates on the considered material, and it was observed that no significant difference can be observed for strain rates lower than _ k ¼ 1:67 Â 10 À3 s À1 .As a consequence, to assume quasi-static conditions, the tests of this paper are performed at a low strain rate, i.e. _ k ¼ 1:67 Â 10 À3 s À1 .In order to distinguish mechanical hysteresis from Mullins effect, a preconditioning test is carried out on each specimen to remove the stress softening from the mechanical behavior.This procedure is classically applied in the literature to ensure test repeatability [6,15,22].In the present study, preconditioning consists in performing 5 cycles of load unload to a stretch of k ¼ 3. Results of uniaxial tensile tests are presented in Fig. 1(a), in terms of the nominal stress versus the stretch k.Fig. 1(b) shows that there is no more stress softening between the 4th and 5th cycles proving that 5 cycles are sufficient to remove the Mullins effect from the mechanical behavior of the silicone rubber.The results are quite similar for planar tensile tests and are not reported here.

Experimental results
Uniaxial tensile test is first carried out with increasing load cycles.These cycles are performed to k ¼ 1:5; 2 and 2:5 (see histogram in Fig. 2).It can be observed in Fig. 2(a) that there is a difference between loads and unloads, and that this phenomenon increases with increasing stretch.It is worth noting that there is still a very few residual strain.
A similar planar tensile test is performed, with three cycles to k ¼ 1:5; 2 and 2:5.The results of this test are presented in Fig. 2(b).The curve is quite similar as for the uniaxial tensile test, the only difference is the stress level that is higher for the planar tensile test.
Two additional non-classical uniaxial tensile and planar tensile tests are performed to study the hysteresis phenomenon.The first test (called 'nc1' in the following) is a load-unload major loop with two hysteresis loops during the load and the unload.The second test ('nc2') is also a load-unload tensile test with a hysteresis sub-loop inside a first loop during the load.Fig. 3 presents the histograms and results of these tests.Once again, it can be seen that the stress level is higher for planar tensile tests.In the 'nc2' test, the size of the hysteresis sub-loop is smaller than the size of the hysteresis loop, whereas for the 'nc1' test the four hysteresis loops seem to have the same size.

Motivation
As mentioned above, hysteresis is a phenomenon whose physical mechanisms are still not clearly understood.Similarly to the other elastomers, silica filled silicone rubbers exhibit a hysteresis loop in terms of the stress-strain relationship [7].The hysteresis loop is not observed when the silicone rubber is unfilled [7], showing that hysteresis is not due to macromolecular network deformation but to fillers, more especially to interactions between fillers and rubber matrix [23].Similar results are observed in other elastomers, for which hysteresis can also be wrongly attributed to viscosity.Indeed, phenomena that are rate-independent in a given range of strain rates and whose effects are different between load and unload induce a hysteresis loop.This is for instance the case in natural rubber, for which hysteresis loop is due to strain-induced crystallization and not to viscosity [24,25].

Decomposition of the material behavior
According to the previous motivations, the model for filled silicone rubber should take into account the macromolecular network elasticity and the hysteresis generated by the difference in the fillers/rubber matrix interactions between load and unload.Classically, models reported in the literature are based on the decomposition of the material behavior into two or three parts.A first hyperelastic part is completed with a viscous part, a damage part or both of them.Even though such a representation is usually used for elastomers, it was also used for other materials, for example to study metallic alloys and magnetization phenomena [26,27].
In the present paper, this approach is adapted to represent time independent hysteresis in elastomers.For this purpose, the model is only composed of a hyperelastic part, represented by a nonlinear spring, and a hysteretic part, represented by infinity of spring-friction slider associations.
The hyperelastic Cauchy stress r hyper and the hysteresis Cauchy stress r hyst are added to obtain the mechanical response of the material.By assuming that the material is incompressible, the total stress writes as follows: where p is the hydrostatic pressure, introduced for taking account the incompressibility, and I is the identity tensor.

Hyperelastic part
The hysteresis part described in the next section can be used with any hyperelastic constitutive equations.As an example, in the present paper, a classical model with a Biderman strain energy [28] is used.This constitutive equation is given by: where c 10 ; c 20 ; c 30 and c 01 are the material parameters, I 1 and I 2 are respectively the first and second strain invariants of the right Cauchy-Green tensor C ¼ F T F, and F is the deformation gradient.These invariants are defined as I 1 ¼ trðCÞ and I 2 ¼ 1 2 ðtrðCÞ 2 À trðC 2 ÞÞ.

General form
For an initially isotropic material, the chains constituting the silicone rubber are assumed to be equidistributed in space.In this case, an infinity of directions in the space is used to model the material behavior [29].However, an infinite number of directions leads to integration problems, and consequently a discrete number of directions is a better solution.This is the reason why the 42 directions proposed by Baz ˘ant and Oh [30] are used in the present study.The directional repartition is classically used to model rubber-like materials, especially to describe Mullins effect [31][32][33][34].Nevertheless, any other spatial repartition can be chosen.By the way, initial directions unit vectors a The main advantage of using a microsphere repartition is that only monodimensional constitutive equations are needed to represent a tridimensional behavior.No tridimensional tensorial generalization of constitutive equations is required as in [19].
The time independent model used for the hysteretic part is shown in Fig. 4. In this representation, the effect of interactions between fillers and the rubber matrix is considered to be different during load and unload.As a consequence, a non-classical friction slider is introduced and its behavior is presented in Fig. 5. First, when loaded in tension, the stretch and stress increase until the stress reaches the limit value of the slider, and stays then equal to this value.When it is relaxed, the stretch and stress decrease until the stress reaches zero and stays constant to this value.So the system comes back to its initial configuration by creating a hysteresis loop as illustrated in Fig. 5(a).Second, when the network is initially compressed, the stretch decreases but the stress stays equal to zero.When the load in compression is relaxed, the network starts to resist and the stress increases until the sliding limit is reached (cf.Fig. 5(b)).This represents a friction slider that is only used in one loading direction.
Each direction is represented by the scheme presented in Fig. 4. The sum of an infinity of linear springs with different sliders can be represented by a global equation.Here, the repartition proposed by Favier and Guélin [26] is used, leading to a global formulation given by: where r ðiÞ is the stress in the considered direction i, E is equivalent to an initial slope, r 0 the maximum reachable stress by the hysteresis part of each direction, and I The stress tensor is finally obtained by summing the contributions of each direction: where x i is the weight corresponding of each direction [30].

Inversion point
To create a difference between load and unload, inversion points are introduced.A similar method was already used by [19,[35][36][37][38][39]].An inversion point is defined as a change in the load evolution (load to unload or conversely).The creation of inversion points is detected by means of an evaluation of the sign of the expression ðr ðiÞ À r ðiÞ inv ÞdI ðiÞ 4 , where dI ðiÞ 4 is the fourth invariant increment.An algorithm introducing the conditions of add and erase of inversion points is presented in Fig. 6 [36].This algorithm is applied for each direction of the microsphere hysteresis model.It is worth noting that the first point of a test is considered as  inversion point, so there is always at least one inversion point.
Once the inversion points introduced, the monodimensional constitutive equation presented in Eq. ( 4) becomes: Two examples are presented to illustrate the inversion points.Fig. 7(a) shows a load-unload curve with creation of an inversion point.At the end of the first load, an inversion point is created (point A).Thus, the unload stress-strain curve AO is the same as the load curve OA, but in a frame with point A as origin and the axes oriented so that r ðiÞ À r ðiÞ inv and ln ffiffiffiffiffi ffi q are positive on AO.Thus, the unload curve is the same as the load curve by a similarity of ratio À1.By this way, the scheme presented by the friction sliders is totally respected and allows us to extend easily the model to more complex cyclic paths.Let us consider a subloop in a first loop, as represented in Fig. 7(b).The first load (OA) is defined as previously by Eq. ( 6).There is next a first unload (AB), defined by Eq. ( 6), with the point A as new inversion point.Then, a second load occurs (BC), where the point B is a new inversion point, followed by a second unload (CD, with C a new inversion point).These second load and unload are also governed by Eq. ( 6), with for inversion point the points B and C respectively.Finally, a third and last load is performed (DE).This reload is defined by Eq. ( 6), considering the point D as inversion point until joining the point C. At this point (C), the subloop CDC is closed and points C and D are no more considered as inversion points.The third load takes now the point B for inversion point in Eq. ( 6) and is the extension of BC until joining point A. At point A, the subloop ABA is closed and there is only O as inversion point.As a consequence, the end of this reload AE is the extension of OA and is governed by Eq. ( 6).The load DE is so made by three different arcs.Table 1 summarizes the equations and inversion points governing the different segments of this curve.This kind of method was already used to model hysteresis in rubberlike materials [19,40], where similar inversion points were implemented in a tridimensional model, with a tridimensional criterion, but it was validated on uniaxial tensile and compression tests only.The originality of the present approach is that the inversion points are only controlled by a monodimensional criterion.
The model is based on the use of friction pads, with inversion points.This idea was already used for other kind of materials, and the thermodynamical dissipation of the scheme was already studied by Guélin [36].The principle stays the same for the present study.

Results and discussion
The constitutive equations previously described are used to simulate a filled silicone behavior with the Matlab software.The parameters are identified from experimental data obtained during the uniaxial tensile test presented in Section 2. They are reported in Table 2.These parameters are then used to simulate the other experimental tests.
The simulation of the first uniaxial tensile test with three increasing cycles is presented in Fig. 8(a).This figure shows that the silicone behavior is well predicted by the model for this experiment as the curves are nearly confounded.The prediction of the hysteresis loops size (calculated as the difference between load and unload) presented in Fig. 8(b) is satisfactory for this test.Fig. 9(a) shows the simulation of the planar tensile test with three increasing cycles.It can be seen that the stress level reached is overestimated in this case, which is due to the hyperelastic Table 1 Equations and inversion points in Fig. 7(b).

Segment Equation List of inversion points
Inversion point used in Eq. ( 6) Values of the material parameters fitted on experimental data.constitutive equation, overestimating the stress for the planar tension for elongation between k ¼ 2 and 4.This was previously explained by Marckmann and Verron [41].However, the hysteresis loop size prediction is acceptable despite the overestimation of the stress, as presented in Fig. 9(b).A more complete constitutive equation could be used to obtain best predictions of the hyperelastic behavior of the material, but the aim of this paper is to show that the hysteresis part can be adapted for all the existing hyperelastic constitutive equations, so the Biderman model is kept as it gives reasonable predictions.The results of the simulations of 'nc1' for uniaxial and planar tensile tests are presented in Fig. 10(a) and (b) respectively.The sizes of the hysteresis loops are very well predicted, as shown in Fig. 10(c) and (d).Fig. 10(e) and (f) show respectively the results of the simulation of the test 'nc2' for uniaxial and planar tensile tests.Even if the size of the greater hysteresis loop is not perfectly predicted in the two cases, the size of the little one is the same for the model as for the experiments and the prediction of the material behavior is acceptable.

Finite element implementation
The model proposed in this paper is developed in the framework of isotropic and anisotropic invariants.To take into account the material incompressibility, a quasi-incompressible formulation of the strain energy function is needed: where 4 are the incompressible invariants and J 2 ¼ I 3 ¼ detðCÞ the square of the volume variation.The model is implemented into the finite element software ABAQUS via a UMAT.Details about the numerical implementation of I 1 ; I 2 ; I 4 constitutive equations can be found in [42][43][44][45].The finite element implementation needs the calculation of the second Piola-Kirchhoff stress tensor and the Lagrangian deviatoric tangent modulus.They are calculated by means of the strain energy derivation.The second Piola-Kirchhoff stress tensor is defined as: The volume variation function considered here is the classical function UðJÞ ¼ KðJ À 1Þ 2 , used in ABAQUS, where K is a material parameter.The general form of the Lagrangian tangent modulus is presented in [46].The isochoric part of the Lagrangian deviatoric tangent modulus is calculated by:  In the case of a strain energy function of the form where W ;i is the strain energy derivation with respect to the ith incompressible invariant of the right Cauchy-Green tensor.This tangent modulus is expressed as a function of the derivatives of the strain energy function with respect to the first, second and fourth invariants.The results of the derivations with respect to the first and second invariants are immediate, they are the derivatives of W hyper of Eq. ( 2).For the fourth invariant the first and second derivatives are given without inversion points respectively by: These derivations change when an inversion point occurs, and in this case the first and second derivatives of the strain energy function with respect to the fourth invariant are given respectively by: It is to note that the tangent modulus has to be expressed in the Eulerian configuration.For this purpose, push-forward operation is carried out [46].
In order to check the model performance and convergence, a non-classical test was carried out on a diabolo structure, meshed with approximately 6000 elements (8-node linear bricks with hybrid formulation, i.e. linear in displacement and constant in pressure) and using the same material parameters as previously (cf.Table 2).It is subjected to a torsion load of 300°, next partially unloaded and reloaded, and finally totally unloaded (see the histogram in Fig. 11).The moment function of the angle is plotted in Fig. 11.
In each element, each direction a ðiÞ 0 of the hysteresis model experiences inversion points during the simulated test.In order to illustrate that, the direction a ðiÞ 0 ¼ z is considered in the initial configuration.Pictures 1-5 represent the number of inversion points in that direction in each element at different times (see the histogram in Fig. 11).It is worth noting that before the end of the first load (picture 1), some elements have already one or two inversion points meaning that the considered direction was loaded and unloaded several times during the torsion load, due to the important rotation of the material directions and to the volume conservation of the considered elements.An increase of the inversion points number is observed during the first three steps, and it decreases after the closure of the hysteresis loop.It should be noted that a maximum principal elongation value of 2.5 is obtained (the parameters were identified for this elongation during the uniaxial tensile test in the previous section), proving the good performance of the proposed model.

Higher hysteresis simulation
The material used in this paper to validate the proposed model presents a weak rate-independent hysteresis.In order to show the possibilities of fitting of higher hysteresis, a simulation with greater parameters of hysteresis is presented in this section.The parameters for the hyperelastic constitutive equation are kept (Table 2), and the hysteretic parameters E and r 0 are taken equal to 5 MPa and 1 MPa respectively.Three increasing cycles of load-unload tensile are simulated, from a stretch equal to 1 to a stretch equal to 1.5, 2 and 2.5 successively.Results of this simulation is presented in Fig. 12, and it can be seen that the difference between load and unload is more important than previously.

Conclusions
In this paper, a microsphere model was proposed to take into account rate-independent hysteresis.This method is based on the decomposition of rubber-like mechanical behavior into a first usual hyperelastic model and a second hysteresis model.This method can be adapted to any existing hyperelastic constitutive equations.This was illustrated by associating the hysteresis part with a classical Biderman constitutive equation.To describe the hysteretic behavior, a non classical friction slider was introduced, representing a different behavior in tension and compression, and it must be initially loaded in tension to bring a contribution.
The predictions of the model are compared to experimental data in case of uniaxial tensile and planar tensile tests.The hysteresis is well predicted for these tests, even though the hyperelastic part cannot perfectly predict the elastomer mechanical behavior for planar tensile tests.
A finite element implementation was carried out on the software ABAQUS.A test on a non-classical geometry with a complex loading and an important number of elements shows the good efficiency of the model.A simulation with high hysteresis loops shows the possibilities of adaptation of the model for different materials.Further investigations about the association of hyperelastic mechanical response, Mullins effect and hysteresis are currently carried out by the authors.

Fig. 1 .
Fig. 1.Preconditioning test to remove the stress softening during an uniaxial tensile test, with (a) the five first cycles and (b) the fourth and fifth cycle, at a stretch of k ¼ 3.

ðiÞ0A
are defined, and initial direction tensors A ðiÞ 0 , current direction vectors a ðiÞ and tensors A ðiÞ and normalized current direction tensor A ðiÞ n are respectively calculated by: ðiÞ ¼ a ðiÞ a ðiÞ A ðiÞ n ¼ a ðiÞ a ðiÞ a ðiÞ k k a ðiÞ k k

Fig. 2 .
Fig. 2. Evolution of the hysteresis during a tensile test for (a) uniaxial tensile test and (b) planar tensile test.

ðiÞ 4 is
the fourth invariant of the right Cauchy-Green tensor C. It is defined by I ðiÞ 4 ¼ tr CA ðiÞ 0 , and it represents the square of the elongation k ðiÞ in direction a ðiÞ 0 .It is worth noting that by affecting a different value of E and r 0 in each direction it is possible to induce anisotropy.

Fig. 3 .
Fig. 3. Results of (a) uniaxial load-unload tensile test and (b) planar tensile load-unload test with four hysteresis loops and of (c) uniaxial load-unload tensile test and (d) planar tensile load-unload test with a hysteresis loop inside another one.
the stress and the square of the elongation of the inversion point.

Fig. 6 .Fig. 7 .
Fig. 6.Management of inversion points algorithm, where r is the current stress, r ðiÞ I is the stress of the inversion point, I ðiÞ is the number of inversion points and dI ðiÞ 4 is the difference between the values of I ðiÞ 4 at the current and the previous step time.

Fig. 8 .
Fig. 8. Uniaxial tensile test with increasing cycles.Experiments (full lines) and simulation (dashed lines) (a) and comparison of the experimental and simulated hysteresis loops size (b).

Fig. 9 .
Fig. 9. Planar tensile test with increasing cycles.Experiments (full lines) and simulation (dashed lines) (a) and comparison of the experimental and simulated hysteresis loops size (b).

Fig. 10 .
Fig. 10.Uniaxial tensile (a) and planar tensile (b) tests with four hysteresis loops and uniaxial tensile (e) and planar tensile (f) tests with a hysteresis loop inside another one.Comparison of the experimental and simulated hysteresis loops size (c), (d), (g) and (h).

Fig. 11 .
Fig. 11.Results of a numerical torsion test on the diabolo.Evolution of the number of inversion points along the cyclic torsion test.

Fig. 12 .
Fig. 12. Simulation of a cyclic tensile test with a great hysteresis.