Large strain rate-dependent response of elastomers at different strain rates: convolution integral vs. internal variable formulations

Two different viscoelastic frameworks adapted to large strain rate-dependent response of elastomers are compared; for each approach, a simple model is derived. Within the Finite Linear Viscoelasticity theory, a time convolution integral model based on an extension to solid of the K-BKZ model is proposed. Considering the multiplicative split of the deformation gradient into elastic and inelastic parts, an internal variable model based on a large strain version of the Standard Linear Solid model is considered. In both cases, the strain energy functions involved are chosen neo-Hookean, and then each model possesses three material parameters: two stiffnesses and a viscosity parameter. These parameters are set to ensure the equivalence of the model responses for uniaxial large strain quasi-static and infinitely fast loading conditions, and for uniaxial rate-dependent small strain loading conditions. Considering their responses for different Eulerian strain rates, their differences are investigated with respect to the strain rate; more specifically, both stiffness and dissipative properties are studied. The comparison reveals that these two models differ significantly for intermediate strain rates, and a closing discussion highlights some issues about their foundations and numerical considerations.


Introduction
Elastomers are often used for damping parts in industrial applications because of their remarkable dissipative properties. Because they can undergo severe mechanical loading conditions, i.e. simultaneous large strain and large strain rates, they are particularly welladapted to shock absorption applications. Nevertheless, in the case of dynamic loading conditions that correspond to moderate or high strain rates, the mechanical response of elastomers is not well-known: it can vary from purely rubber-like to glassy (Yi et al. 2006;Sarva et al. 2007), but also exhibits wave propagation (Niemczura and Ravi-Chandar 2011a, 2011b, 2011c. Modelling of elastomers has been a subject of numerous studies. On the one hand, large strain hyperelastic constitutive equations are well-suited for static, but are also considered for very large strain rate responses as shown by Hoo Fatt and Bekar (2004), for example. In both cases, classical hyperelastic models are employed (neo-Hookean, Mooney-Rivlin, Ogden, etc.; see Holzapfel 2000;Marckmann and Verron 2006). On the other hand, the rate-dependent small strain response of elastomers is usually predicted in the general framework of linear viscoelasticity. In this case, integral models issued from the Boltzmann superposition principle and the rheological Maxwell, Kelvin, Zener, Poynting-Thomson models are equivalent (Ferry 1980;Christensen 1982;Tschoegl 1989;Wineman and Rajagopal 2000). Modelling the large strain rate-dependent response of rubber-like materials is highly more complex, and there is not a unique framework for large strain viscoelasticity; then numerous models have been derived. Figure 1 summarizes those different loading conditions. Phenomenological viscoelastic rubber-like constitutive equations are derived following two different frameworks. Here, only the basics of these approaches are recalled; for a recent and comprehensive state of the art, the reader is referred to the remarkable introduction proposed by Amin et al. (2006). The first approach is founded on the Green and Rivlin (1957) expansion theory and the Finite Linear Viscoelasticity theory which extends the Boltzmann superposition principle to finite strain (Coleman and Noll 1961). The total stress is the sum of the long-term low strain rate response and the viscous overstress which is expressed in terms of time convolution (or hereditary) integrals of the material history. One of the most recognized model founded on this approach is the K-BKZ model which has been widely used in rheology (Kaye 1962;Bernstein et al. 1963;Zapas and Craft 1965;Tanner 1988). Solid models have also been proposed by Christensen (1980); Chang et al. (1976); Morman (1988); Sullivan (1987);Holzapfel (1996); Haupt and Lion (2002). It is to note that a comparison of hereditary integral models for viscoelasticity has been recently proposed by Ciambella et al. (2010). The second approach is based on the semi- Fig. 1 Applicability ranges of models nal paper of Green and Tobolsky (1946) who derived the viscoelastic counterpart of the neo-Hookean model. Their idea was rationalized by introducing internal variables: following the work of Lee (1969) for plasticity, authors introduced the internal variables through the split of the deformation into an elastic and a viscous parts (Sidoroff 1974;Lubliner 1985). Thanks to the development of relevant algorithms for finite element analysis (Simo and Hughes 1998), this approach has been extensively used in the bibliography because of its thermodynamical consistency (Lion 1997;Reese and Govindjee 1998;Huber and Tsakmakis 2000;Bonet 2001).
The aim of the present paper is the comparison of the two above-mentioned theories for intermediate strain rates as depicted in Fig. 1. In this way, we derive two simple models for large strain viscoelasticity: the first one referred to as a Convolution Integral Model (CIM) is a solid extension of the K-BKZ model, and the second one, referred to as an Internal Variable Model (IVM), is the large strain counterpart of the Zener model as proposed by Huber and Tsakmakis (2000). Then, for equivalent sets of material parameters, their uniaxial responses, i.e. their stiffness and associated dissipation, are compared and discussed with respect to the strain rate.

Two simple models for large strain viscoelasticity of rubber-like materials
We consider an elastomeric material; it is assumed homogeneous, isotropic and incompressible. We focus on the viscoelastic nature of its mechanical response and we consider two simple models: the former based on the convolution integral approach, and the latter based on internal variables theory.

General derivation
The model considered here is the sum of a large strain hyperelastic model for the elastic part and a fluid K-BKZ model for the viscous part. It is to note that this model is different from the solid version of the K-BKZ formulation that was proposed in the paper of Bernstein et al. (1963).
The key mechanical quantity of the convolution integral approach is the relative deformation gradient F t (τ ), i.e. the deformation gradient at the current time τ in the current configuration (C τ ) with respect to the final deformed configuration (C) at time t , as described in Fig. 2: In a similar way, we define F τ (t), the deformation gradient in (C) with respect to the deformation gradient at the current time τ : Introducing the right Cauchy-Green strain tensor C τ (t), and the Green-Lagrange strain tensor E τ (t): where I is the 3 × 3 identity tensor and (·) T stands for the transposition, the K-BKZ model leads to the following expression of the Cauchy stress tensor (Bernstein et al. 1963) in which p stands for the Lagrange multiplier due to the incompressibility constraint, referred to as the hydrostatic pressure, and U 1 is a strain energy density per unit of volume in the current configuration (C τ ). Considering objectivity and incompressibility, this strain energy is a function of the two first strain invariants of the right Cauchy-Green strain tensor C τ (t), I Cτ (t) and II Cτ (t) . Thus where the two first invariants of a tensor X are I X = trX and II X = 1 2 (trX) 2 − trX 2 .
Moreover, as the material is isotropic, the strain energy U 1 can be written as a function of the left Cauchy-Green strain tensor B τ (t) that is equal to C −1 t (τ ), i.e. in terms of their invariants I C −1 t and II C −1 t . Finally, recalling that (see, for example, Eq. (5.93) in Holzapfel 2000) and after some simple algebraic manipulations, in which p (t) is the modified hydrostatic pressure that will be simply denoted p(t) in the following. Equation (8) represents the classical fluid K-BKZ constitutive equation. In order to consider solids, we simply add a hyperelastic contribution defined by a strain energy function written in terms of the two first strain invariants of the left Cauchy-Green strain tensor of the whole deformation, B(t) = F(t)F T (t), U 0 (I B , II B ); finally, the corresponding Cauchy stress tensor is:

A simple model
In order to define the simplest model based on the previous derivation, we consider that -The elastic part obeys the neo-Hookean constitutive equation: -For the viscous part, the effects of time and strain are separable. Thus, the strain energy U 1 can be chosen as the product of a decreasing exponential time function with only one relaxation time τ R and of a neo-Hookean strain energy density in terms of C −1 t (τ ) Thus, our CIM involves only three material parameters: two stiffnesses g 0 and g 1 (in MPa), and one relaxation time τ R (in seconds); and the stress-strain relationship reduces to This equation can be integrated by parts to explicitly introduce the influence of the strain rate ∂C −1 (τ )/∂τ :

General derivation
The model is based on the multiplicative split of the whole deformation gradient F between elastic, F e , and inelastic, F i , parts This split leads to the definition of an intermediate configuration (C i ) between the reference undeformed configuration (C 0 ) and the final deformed configuration (C), as shown in Fig. 3(a). The intermediate configuration is assumed stress-free and represents the state of the material once loading is infinitely fast removed. The simplest internal variable constitutive equations for the viscoelastic elastomers were derived by Huber and Tsakmakis (2000). The authors studied the two large strain counterparts to the classical Standard Linear Solids: model A stands for the Maxwell form of this model, i.e. a Maxwell model and a dashpot in parallel, and model B for the Zener form, i.e. a Kelvin-Voigt model and a spring in series. In the present paper, we consider the former one, because the model B does not reduce to a hyperelastic model when strain rate tends As rubber is assumed incompressible, det F = 1; moreover, we also consider that each deformation is incompressible: det F e = det F i = 1. We now introduce the total strain energy density per unit of undeformed volume Ψ ; it depends on both F and F e . With the incompressibility assumption, the volume in (C 0 ) and (C i ) remains the same and then Ψ can be written as where Ψ 1 is the strain energy involved in the deformation between (C 0 ) and (C), and Ψ 2 is the strain energy involved in the deformation between (C i ) and (C); they correspond to the two springs in Fig. 3(b). With the objectivity principle and the isotropic nature of the material, these strain energy densities can be written as and where B e is the left Cauchy-Green strain tensor associated with the elastic deformation, i.e. F e F T e . Application of the Second Principle of Thermodynamics leads to where D int is the internal dissipation and D is the rate of deformation tensor. With the specific form of the strain energy density, it leads to After some algebraic manipulations, one establisheṡ where L is the velocity gradient defined asḞF −1 , anḋ where D i is the inelastic strain rate tensor defined as With these expressions and basic algebraic manipulations, the second principle Eq. (19) becomes Next we apply the Coleman and Noll (1963) procedure that consists in satisfying this inequality for all possible deformation processes, i.e. for all possible F and F i which satisfy the incompressibility constraints. It leads to the stress-strain relationship and the dissipation Similarly as in Huber and Tsakmakis (2000), we consider the simplest sufficient condition to satisfy Eq. (25): where (·) D stands for the deviatoric operator, i.e. (·) − tr(·)I/3. This equation gives the following evolution equation for the elastic straiṅ Considering now the isotropic nature of the material, each strain energy density depends on the two first invariants of the corresponding strain tensors and the model reduces to

A simple model
In order to define the simplest model based on the previous derivation, we consider that both strain energy functions are of the neo-Hookean type and Thus, our IVM involves only three material parameters: two stiffnesses C 1 and C 2 (in MPa) and a viscosity η (in MPa.s); the stress-strain relationship reduces to and the evolution equation of B e iṡ 2.3 Comparison of the models: choice of the material parameters In order to compare the two above-mentioned simple models, it is necessary to consider material parameters that render them comparable. In this way, we investigate three limit responses of the material: 1. The quasi-static limit response, i.e. very low strain rates, 2. The infinitely fast limit response, i.e. very high strain rates, 3. The linear viscoelastic limit response, i.e. small strain.
These three limit responses will lead to three equations to relate the two sets of material parameters, i.e. (g 0 , g 1 , τ R ) and (C 1 , C 2 , η).

Quasi-static response
We consider infinitely slow deformation processes such that the strain-stress path corresponds to the static path of equilibrium states. For the CIM, we consider the method proposed by Christensen (1980) that consists in simulating accelerated or decelerated strain histories with a change in time scale of the strain: C(αt) with α > 1 corresponds to accelerated strain histories and α < 1 corresponds to decelerated strain histories. In this case, the Cauchy stress tensor Eq. (13) is then which can also be written as Then, the quasi-static response is obtained by considering the infinitely decelerated strain history, i.e. α → 0: For the IVM, the quasi-static response is defined by the nullity of the strain rate tensor: As this equation must be satisfied for all possible deformation processes, i.e. for all possible F e and F i , all strain rate tensorsḞ e ,Ḟ i , L,Ḃ e ,. . . must be null. Thus, Eq. (33) reduces to which means that B e is a spherical tensor q being a scalar quantity. Thus, the quasi-static limit response of the IVM is Finally, comparing Eqs. (36) and (40) gives the following relationship between the material parameters
For the IVM, when the strain rate D tends to infinity,Ḟ tends also to infinity because F is finite. As viscous effects do not take place, F i tends to I andḞ i tends to zero, thus F e tends to F, and F e tends to infinity. Then, as B e tends to B, we have σ | D→∞ = −pI + 2(C 1 + C 2 )B.
Finally, comparing Eqs. (42) and (43), and recalling Eq. (41) leads to the following relationship between the material parameters

Linear viscoelastic response
In order to compare the small strain responses of the two models, we consider the uniaxial extension defined by the deformation gradient where λ is the stretch ratio, e 1 is the extension direction, and e 2 and e 3 are the transverse directions. With the small strain assumption, the kinematical (linearized) quantities are where the subscript (·) l stands for the linearized counterpart of the given tensor and ε is the true strain in the extension direction. Similar results apply to F e . Once the linearized quantities defined, we now consider dynamic processes and introduce the classical complex notation for both strain and stress in the extension direction ε = ε 0 exp(iωt) and σ = σ 0 exp i(ωt + δ) where ε 0 and σ 0 are the strain and stress amplitude, ω is the loading frequency and δ is the phase angle.
Considering that the transverse stresses σ 22 and σ 33 are null, the CIM stress-strain relationship Eq. (13) reduces to or with the complex quantities, to Integrating the last right-hand side term we easily obtain the complex modulus E * Thus, the storage and loss moduli, i.e. the real and imaginary parts of E * , respectively, are and With the aforementioned linearized quantities, the IVM stress-strain relationship Eqs. (32) and (33) simplify into the following set of equations Introducing the complex representation of the elastic deformation ε e = ε e0 exp i ωt + δ (59) in which ε e0 is the elastic strain amplitude and δ its phase angle with ε, Eq. (58) reduces to and permits expressing the elastic strain in terms of the total strain: Reporting this equation into Eq. (57) leads to the expression of the complex modulus and then to the storage modulus and to the loss modulus: Finally, with the help of the previously established relationships between the material parameters Eqs. (41) and (44), equalizing the storage moduli Eqs. (55) and (63), and of the loss moduli Eqs. (56) and (64) gives

Results and discussion
In order to compare the two models, we study here their stress-strain responses in uniaxial extension. First, we define the sets of material parameters that will permit this comparison. By examining the two models and their response at low and high strain rates, it appears that the parameters g 0 in the CIM and C 1 in the IVM represent the stiffness for quasi-static loading conditions, and that the parameters g 1 in the CIM and C 2 in the IVM represent the additional stiffness reached at high strain rates. The third parameters, respectively τ R for the CIM and η for the IVM, account for the dependence of the behaviour on the strain rate: they drive the increase in stiffness at high strain rates, and both size and shape of the hysteresis loop during loading-unloading response. In order to respect the previously  Table 1. These values are not founded on experimental measurements, and the large values for g 1 and C 2 are chosen in order to emphasize the phenomena (hysteresis, strengthening) and to highlight the differences between the two approaches.
As mentioned above, we consider uniaxial cyclic tension. In order to investigate the influence of strain rate on the models, we consider cycles at constant true strain rate, i.e. in terms of the strain rate tensor D. In fact, we consider the strain rate that is actually borne by the material in the actual deformed configuration. Noting λ(t) the stretch ratio in the extension direction e 1 and recalling the incompressibility assumption, the deformation gradient is and the strain rate tensor defined as reduces to Considering that this tensor is constant, we introduce the strain rate α that is equal tȯ λ(t)/λ(t); and then the prescribed stretch ratio is λ(t) = e αt during loading, λ max e −α(t−tmax) during unloading, where λ max is the maximum stretch during the cycle and t max is the corresponding time given by log λ max /α. Loading curves for λ max = 4 and different values of α are presented in Fig. 4. The stretch history being defined, we now determine the stress for both models. In uniaxial extension, the unique non-zero stress is σ 11 and the equation σ 22 = σ 33 = 0 is used to determine the value of the hydrostatic pressure p. Thus, recalling Eq. (13), the CIM stressstrain response is and recalling Eqs. (32) and (33), the IVM stress-strain response is The (Cauchy) stress-strain curves obtained for different strain rates α are shown in Fig. 5 for both the CIM and the IVM; the limit curves of quasi-static and infinitely fast loading conditions are also shown in this figure.
At first sight, the two models exhibit the same response: strengthening and a change in the hysteresis loop as the strain rate increases. Nevertheless, quantitative differences exist. Firstly, strengthening differs: as an example the stress at λ = 4 is larger for the CIM than for IVM for all strain rates but this is not the case at λ = 1.5. Secondly, the hysteresis loop is larger for IVM than for CIM. These differences will be discussed in the two following sections that focus on the stiffness of the models and on their dissipative properties, respectively.

Comparison of the loading curves
In this section, only the loading part of the stress-strain curves are discussed. In order to study the difference in the responses, we consider the ratio of stress of both models, i.e. σ 11IVM /σ 11CIM , for different strain rates as shown in Fig. 6. First, it is to note that the quasistatic and the infinitely fast ratios are equal to 1 for all stretch values because the models coincide in these two cases; they are not shown in the figure. For intermediate strain rates, the three curves differ with no similarity in their evolution. We can only mention that the three curves tend to 1 as the stretch tends to 1, i.e. for small strain. This result is obvious because we choose the material parameters such as the two models coincide for small strain.
As the two models cannot be easily compared, we choose to investigate the influence of the strain rate on their responses, individually. Figure 7 presents the evolution of the reduced Cauchy stress with respect to the stretch for CIM (Fig. 7(a)) and IVM ( Fig. 7(b)); in order to render comparable the results, the stress is normalized with respect to the quasi-static (α = 0) stress of the same model. These curves lead to the following comments: -For both models, all curves tend to 10 as the stretch tends to 1; it is easily explained by recalling that the initial stiffness of the models is the one of the elastic infinitely fast response which is ten times larger than the stiffness of the quasi-static response, i.e. (g 0 + g 1 )/g 0 = 10 for CIM and (C 1 + C 2 )/C 1 = 10 for IVM. For low strain rates, the stress ratios of both models quickly decrease to attain the horizontal asymptote σ (α)/σ (0) = 1, indeed the stress-strain response tends to the quasi-static limit response. As the strain rate increases, the decrease in the response becomes less important, and CIM and IVM exhibit different results that are emphasized in the following. -For the CIM, all curves monotonically decrease without inflection point. This behavior can be explained by the decreasing exponential function of time in the integral: for high strain rates, the time scale is revealed very small as compared with the relaxation time, and then the stress does not noticeably decrease in the range of stretch. -For the IVM, all curves also decrease; nevertheless, we observe an inflection point on the curve α = 100 s −1 at about λ = 1.5. This can be explained by the behavior of the dashpot during stretching. At the beginning of loading, the dashpot is not stretched because a certain duration is necessary before its activation, and the spring in series with the dashpot is highly stretched (see Fig. 3(b)). Thus, the corresponding stress-strain curve tends to the infinitely fast response. Afterwards, as the dashpot becomes active, the spring is not stretched more, the stress in it remains constant and the corresponding elastic stretch ratio, i.e. λ e drops.

Dissipative properties
Dissipation of the models are quantified by the shape and size of the hysteresis loop between loading and unloading stress-strain curves. As emphasized above, both shape and size change with the strain rate. Dissipation is null as the strain rate tends to zero and to infinity because the two models reduces to their hyperelastic counterparts. In order to discuss the dissipation, we calculate the energy dissipated over one cycle as follow where the time at the beginning of the cycle is set to 0 and t e is the time of the end of the cycle. Similarly, we define the elastic energy stored in the material during loading, i.e. the area under the loading curve as Figure 8 presents the ratio D/E as a function of the strain rate for both models. First, we notice that the maximum value of the ratio is about 1, which would mean that all the stored energy is dissipated. Obviously, it is not representative of real elastomers; it is due to the large value of the second stiffness parameters, i.e. g 1 and C 2 , and also because the energy stored during unloading has been neglected. Nevertheless, the ratio represents a relevant measure of dissipation. As expected, the energy dissipated tends to zero for very low and very high strain rates, and increases for intermediate strain rates. For strain rates smaller than 5 s −1 , the responses of the two models are nearly superimposed. For higher strain rates, differences take place: the IVM exhibits more dissipation and over a wider range of strain rate than the CIM. We also investigate the influence of the maximum stretch λ max on the dissipative properties. As shown in Fig. 9, changing the maximum stretch of a cycles leads to different ratios of dissipated energy for both the CIM ( Fig. 9(a)) and the IVM ( Fig. 9(b)). For a given maximum stretch, the ratio is identical for both models; and for a given model, it differs as the maximum stretch changes. The maximum dissipation decreases as the maximum stretch increases for both models. Nevertheless, there is an important difference between the two approaches. For a given dissipation ratio, the CIM admits a limit value of strain rate, e.g. D/E < 0.4 when α > 80 s −1 . Oppositely, the curves of dissipation ratio for the IVM admit a decreasing maximum value and are only shifted to the right as λ max increases.

Closing remarks
As a summary, we compare two large strain viscoelastic frameworks, i.e. convolution integral and internal variable models, by studying the response of two simple constitutive equations that are equivalent for very small and large strain rates, and small strain. It has been shown that the responses of these models differ for intermediate strain rates. Nevertheless, it is quite difficult to determine which model must be preferred for dynamic applications. To help the choice, we can invoke two aspects that were not discussed previously: the foundations of the approaches and numerical considerations.
The convolution integral approach is based on the extension to large strain of the wellestablished Boltzmann superposition principle: the time dependence is accounted for by a decreasing time function in a convolution integral (Christensen 1982). This approach is the most natural extension of linear viscoelasticity to large strain. A large relaxation spectrum can be easily handled by considering multiple relaxation times and a Prony series for the relaxation function. One of the major criticism is that it is not founded on the calculation of the dissipation; thus the material functions must be carefully chosen and the second principle of thermodynamics must be verified a posteriori. The internal variable approach is founded on a completely different point of view: the split of the deformation gradient was proposed by Sidoroff (1974) as a transposition to viscoelasticity of the plastic split due to Lee (1969). Obviously, there is no problem with the dissipation: the constitutive equations are directly derived from the application of the second principle of thermodynamics. Nevertheless, in the case of viscoelasticity the split of the deformation gradient is questionable: the intermediate reference configuration defined by F i is not an equilibrium configuration, and thus it tends to relax to the reference configuration. Indeed, it is considered as an equilibrium configuration if the time scale of the phenomena is highly smaller than the relaxation time. Indeed, if the time scale of the phenomena evolves, the inelastic part of the deformation gradient may become elastic! Solutions have been proposed to consider large relaxation spectra for quasistatic problems (Knauss and Emri 1981;Govindjee and Simo 1993;Holzapfel 2000).
From a numerical point of view, on the one hand, it is accepted that IVMs are well-suited for implicit (static in most of the cases) finite element analysis considering similar methods as those employed for plasticity (see, for example, Simo and Hughes 1998). On the other hand, the CIMs approach is well-adapted to dynamic explicit finite element analysis, it is quite easy to implement explicit recurrence formula for the computation of the convolution integrals (Feng 1992;Shrivastava and Tang 1993;Verron et al. 2001;Marckmann et al. 2001).
Thus, with these remarks, the choice of a constitutive equation framework for large strain viscoelasticity of rubber-like materials must be discussed in regards with the range of strain rates that is involved in the problem; for dynamical applications, we believe that the CIM framework provides the best compromise between theoretical requirements and numerical efficiency.