Modeling the Full Time-Dependent Phenomenology of Filled Rubber for Use in Anti-Vibration Design.

Component design of rubber-based anti-vibration devices remains a challenge, since there is a lack of predictive models in the typical regimes encountered by anti-vibration devices that are deformed to medium dynamic strains (0.5 to 3.5) at medium strain rates (0.5/s to 10/s). An approach is proposed that demonstrates all non-linear viscoelastic effects such as hysteresis and cyclic stress softening. As it is based on a free-energy, it is fast and easily implementable. The fitting parameters behave meaningfully when changing the filler volume fraction. The model was implemented for use in the commercial finite element software ABAQUS. Examples of how to fit experimental data and simulations for a variety of carbon black filled natural rubber compounds are presented.


Introduction
Suspension is one of the most important systems affecting a vehicle's ride quality, and it is, therefore, a key factor in determining a vehicle's performance. The modelling of bushing force in the automotive industry plays an important role in predicting the dynamic behaviour of the suspension system. The viscous behaviour of elastomers makes the use of rubber components essential to reduce the level of vibrations that are transmitted to the passengers in the vehicle. In particular, natural rubber is the material of choice for these types of engineering applications. For automotive suspension components, carbon black is also typically used as a filler to improve the rubber's mechanical properties such as the fatigue behaviour or the stiffness. For bushing applications, quasi-static deformations up to large strains are of major interest. In these particular conditions, the rubber exhibits non-linear viscoelastic behaviour such as the Mullins Effect, cyclic stress softening, hysteresis and induced anisotropy [1][2][3][4][5]. This characteristic response is related to the molecular microstructure, but it is not yet totally understood [6]. Furthermore, these effects are also present in unfilled strain-crystallizing elastomers as natural rubber [7][8][9][10]. Modelling of filled rubber introduces difficulties related to its nonlinear and incompressible behaviour. Many models try to capture this behavior from a range of different perspectives. The multiscale nature of polymers is reflected by different computational methods for specific length and time scales: quantum (∼10 −10 m, ∼10 −12 s), atomistic (∼10 −9 m, ∼10 −9 -10 −6 s), mesoscopic (∼ 10 −6 m, ∼10 −6 -10 −3 s) or macroscopic scale (∼10 −3 m, ∼1 s) [11][12][13]. At the top level, phenomenological macroscopic models can roughly be classified into three categories [14] damage models, rheological models with serial and parallel combination of elastic and viscous elements and the constitutive equations based on a rubber elasticity model. Damage models [15][16][17] are unable to differentiate between an unloading and a subsequent reloading with one consequence being that they cannot predict the cyclic stress softening phenomenon. Essentially, there is a limitation in fitting the experimental behaviour as a limitation of the background theory. Rheological Framework models are models with elastic and viscous components. The Parallel Rheological Framework ( [18,19]) is capable of reproducing the full range of nonlinear viscoelastic effects that are under examination, but it requires a large number of parameters that are very sensitive to the input test data. The dynamic flocculation model (DFM) assumes the breakdown and reformation of filler aggregates as the main mechanism [20]. Other physically-based works focus on the binding or sliding of polymer chains on the fillers surface [21]. Recent Molecular Dynamics (MD) simulations indicate that both phenomena may play a role [22,23].
Most of these models include softening effects using non-time dependent mathematics, e.g., by softening the material based on maximum strain or stress measures. In contrast, experiments clearly show that softening proceeds by repeatedly stretching the material to a predefined load (cyclic stress softening) or by holding it at constant stretch (stress relaxation) [24]. More specifically, stress decreases logarithmically [25][26][27]. In a recent review article [14] it was found that no predictive model was able to fully reproduce cyclic stress softening to a satisfying quality. In the next sections a new approach based on a recently proposed model [28] is presented. It is able to account for non-linear elasticity and strain history effects.

Free Energy Density
The underlying free energy density is derived on the basis of the extended non-affine tube model of rubber elasticity [29,30] which was shown to be the best compromise between fitting quality and number of parameters [31]. In a simplified form it reads being modified invariants of the left Cauchy Green tensor, expressed by the principal stretches λ i . The straightforward simplifications done to arrive at Equation (1) are outlined in [32]. The parameter G c is called crosslink modulus and scales the contribution of crosslinks to the mechanical response. Accordingly, G e scales proportionally to the number of trapped entanglements of the polymer [33]. The third parameter n measures the number of statistical segments between network nodes, which is a measure of elastically effective chain length. For example, natural rubber has a segment length of about 0.934 nm [34], corresponding to roughly two isoprene units. Fillers are introduced by assuming that they amplify strain heterogeneously within the matrix by local amplification factors X [28]. To avoid problems arising from frame references the amplification factors directly act on the invariants. In the case ofĨ 1 this is reasonable, as it represents the norm of a hypothetical length within the material. The amplified energy density is then calculated as a superposition of differently amplified domains where the distribution of amplification factors is given by The parameter χ gives the width of the distribution and C defines a plateau value around X = 1. For C = 0 Equation (5) is equivalent to its equivalent in [28]. It is normalized to the interval [X min ,X max ], which is motivated from conservation of the number of rubber-filler structures. Carrying out the integral given by Equation (4) generates hypergeometrical functions. Additive splitting of the integrand, as is explained in [32], gives a good approximation containing only elementary functions: From Equation (6) the hyperelastic stresses can be derived, if X min and X max are known. The minimum amplification factor is set to X min = 1, which corresponds to the assumption that there are non-amplified domains (e.g., without filler) within the material.

Stress Softening
In [28] the maximum amplification factor is defined to be a monotonically decreasing function of the all-time maximum of the modified first invariantĨ 1 . This naturally introduces immediate stress softening if the material surpasses the previous maximum strain. As described in the introduction stress softening is not immediate, but proceeds logarithmically or according to a slow powerlaw. Generally, logarithmic relaxation in stress f ∼− log t can be generated by differential equations of the form d f /dt∼− exp | f |. From a physical point of view this can be understood in terms of force-induced hopping over a potential barrier, as formulated by Kramers [35]. It is reasonable to assume that softening is predominantly happening in the most stretched domains. A crude approximation of the stress-scaling in this regime gives where the entanglement part of Equation (1) was neglected, because its influence is small at high strains. The prefactor λ = Ĩ 1 is a frame independent measure of the systems stretch. Heuristically, stress is also directly scaled by the amplification factor, such that d f /dt∼ dX max /dt. Using again the logarithm-generating differential equation d f /dt∼− exp | f | and introducing scaling constants c 1 and c 2 this gives The timescale of zero-load softening given as τ soft = exp(c 2 ) s was put in the exponential for numerical reasons. In practice, it should be sufficiently long that no significant softening happens on the timescale of the simulation, if no load is imposed. This means that forĨ 1 → 0 the result of Equation (8) is usually small and the maximum amplification factor X max stays almost constant. When the model approaches divergence atĨ 1 X max /n → 1 the second summand inside the exponential increases and induces a decrease of X max . It shall be noted that a more careful derivation was carried out in [32] which relates c 1 and c 2 to temperature and allows the modeling of Mullins effect recovery [6] at elevated temperatures.

Viscoelasticity
Hysteresis is modeled by a single Maxwell element. The problem treated in this work requires only a narrow range of strain rates such that this procedure is justified. Viscoelastic stress is evolved according to the differential equation where σ σ σ vis is the viscoelastic stress contribution, τ c is the relaxation time scale, and σ σ σ el is the elastic Cauchy stress derived from the free energy density given by Equation (6).

Implementation
The model was implemented in Matlab for the hypothetical condition of fully incompressible material, considering one material point loaded cyclically in the uniaxial direction and in a compressible full three dimensional form for the finite element analysis. For incompressible materials J = I 3 = 1 and the Cauchy stress is related to the free energy density: The scalar p is the hydrostatic pressure (an indeterminate Lagrange multiplier) that can be determined only from the boundary conditions. B is the left Cauchy-Green deformation tensor. For compressible materials it is useful to split the deformation into distortional and volumetric (dilational) parts via a multiplicative decomposition.
The J 1/3 and J 2/3 terms are related to the volume changes, while the• terms are a modified gradient and strain related to the distortional deformations. These decompositions express the stored energy as W(Ī 1 (I 1 , J),Ī * (I * , J), J), that can be decomposed into volumetric and isochoric (distortional) responses (Equation (14)) and give the Cauchy stress [36] as in Equation (15) The single Prony series element was integrated explicitly and the total stress is defined as: where φ is a fitting constant which is assumed to increase with filler volume fraction. Figure 1 shows the algorithm developed to implement the model as a user-defined material subroutine (VUMAT) coded in Fortran for Abaqus/Explicit. The model requires 10 parameters, plus the bulk modulus in case of compressible material.

Materials and Experiments
Seven compounds supplied by TARRC (Tun Abdul Razak Research Centre, Brickendonbury, Hertford, UK), were examined. Each compound contained a different volume fraction of carbon black (FEF N550) that had been mixed into natural rubber (NR, SMR CV60). The compound formulations in phr (part per hundred of rubber by mass) are given in Table 1. Characterisation of all this materials under cyclic tensile tests is available in [14] and it is briefly recalled in this section for the sake of completeness. Uniaxial tension tests were conducted using an Electropulse Instron Test Machine. Cyclic loading, unloading and reloading tests were performed using four different strain rates of 0.5/s, 1.5/s, 3/s and 6/s. All the specimens had a gauge length of 12 mm and a cross section of 2 mm × 0.5 mm. The nominal strain, ε NOM , was determined as the ratio of the axial displacement to the original length of a specimen (12 mm). This test was difficult to conduct because the large strain amplitude and rate approached the limit of the test machine. To reach the desired strain rates, the specimens had to be short. In addition, the strain rate was too fast to be reliably measured using the optical strain measuring device. The tensile force was measured using a 1 kN load cell.
A subset of the results are shown in Figures 2 and 3 for two of the materials. The results display strong nonlinearity, large hysteresis, complex Mullins behaviour (Figures 2 and 3b), cyclic stress relaxation ( Figure 3a) and permanent set. Full detail of all the experimental results can be found in [14].
In the range of interest the strain rate effect is not a dominant feature, so in the following analysis it was ignored.

Sensitivity Analysis
This section shows how each parameter effects the output stress. The initial parameters shown in Table 2 are used, it is an arbitrary set. The parameter X max,0 = X max (t = 0) is the boundary condition of the differential Equation (8).  Figure 4 shows the effect of three of those parameters G c , G e , X max,0 on the predicted outputs. Figure 5 shows the effect of parameters χ, C and c 1 . Figure 6 shows the effect of the parameters c 2 , n and τ c . Figure 7 shows the effect of the parameter φ. G c and G e , which respectively represent the crosslink modulus and entanglement modulus in the extended tube model, influence the stiffness of the material. The parameter X max,0 determines the value of X max before any deformation in the system, and it influences cyclic stress relaxation phenomena. The parameters χ and C characterise the power law distribution of the amplification factor. They influence the response of the material in the loading paths, modifying the the level of stress and the dissipated energy and defining the level of continuum damage due to cyclic stress softening. n, c 1 ,and c 2 define the relaxation of the maximum amplitude factor to predict the continuum damage, so they operate on the response of the material at a strain higher than 1 and envelop the cyclic stress relaxation process. The constant c 2 scales the zero-load relaxation of the material (which is usually very small), and thus determines, together with c 1 the stress threshold to overcome before softening starts. The parameter φ balances elastic stresses (including softening) and viscous effects. The timescale τ c defines the amount of dissipated energy.  (e) (f) (e) (f) (e) (f)

Model Fit to Experimental Data
Model parameters were determined by minimising the mean square error between the model predictions and experimental data. (c) Figure 8. Comparison of the new model with experimental data loaded with cyclic uniaxial test with uniform maximum strain amplitude [14]. Subfigures (a-c) show respectively the behaviour for NR2, NR40 and NR60.

Effect of Carbon Black Content on the Parameters Used in the New Model
The uniaxial tests, where the compound were stretched cyclically with the step up of maximum strain amplitude reached, represent a perfect test to highlight all the non linear viscoelastic behaviour under examination (Mullins effect, cyclic stress softening and permanent set) (Figure 9). The seven compounds were fitted starting from the material with the smaller amount of filler and using the new set of parameters as the guess set for the following compound. The best-fit parameters for different compounds loaded with this strain history show a robust trend with the true carbon black filler volume. The trend of the parameters with the true Carbon Black (CB) filler volume is presented in Figure 10.
• G c is the modulus resulting from network crosslink and it is relatively constant with the CB volume fraction. It is found to be in a typical range for Natural Rubber vulcanizates. There is a drop in crosslink density at the highest filler loading, whose origin is not clear and may be due to experimental uncertainties and/or parameter correlations. • G e , the modulus resulting from network entanglement, is monotonically increasing with the CB concentration. This is surprising and lacks a clear explanation up to now. Probably this parameter somehow captures the increasing low-strain stiffness (Payne effect), even though it was not designed to fulfil this purpose. • φ, represents the scaling of the elastic and inelastic part, is proportional to the true filler volume fraction.
• n, the distance between the network nodes is in a typical range, too. It is constant for the samples containing less than 30 vol. % of carbon black and then slightly decreases, roughly in accordance with the small drop in crosslink modulus at these filler concentrations. • χ, the exponent of the amplification factor distribution, is approximately constant. • τ c , is decreasing with filler volume fraction. • X max,0 is approximately constant with filler volume fraction. • C is increasing. The parameter was introduced primarily to modify the shape of the virgin loading curve. A value k > 0 generates a rather linear curve, as is observed for many highly filled compounds. • c 2 decreases with volume fraction and appears to asymptotic. From Equation (8) it can be seen that exp c 2 defines the timescale on which X max,0 relaxes without load. The value obtained from fitting here create an optimal model for the timescale of the fit, but it may fail for longer simulation times. An optimal determination of c 2 and c 1 requires a stress relaxation characterisation in the fitting data. • c 1 increases modestly with volume fraction, probably representing that the rubber-filler structures to be broken down during softening become more rigid at higher filler loadings.
Moreover, it can be seen that the parameters G e , C, c 1 , φ and τ c show some kind of discontinuity in the range of 9-13 vol. % carbon black loading. This roughly corresponds to the percolation threshold of the filler network inside the polymer matrix [37].

Benchmark Tests
This section shows the benchmark tests on a single element model with a simple prescribed traction loading. The model is a cube of size 0.5 mm and is meshed with 1 linear brick element with reduced integration (C3D8R). A velocity of 1 mm/s is applied in direction 1 (along the x-axis). The element was loaded for 5 cycles at constant max strain amplitude using 10 steps. The velocity was positive (v = 1 mm/s) in the odd steps to stretch the element and negative (v = −1 mm/s) in the even steps to unload the element. The job was run using double precision. Material parameters are shown in Table 2. Figures 11 and 12 show the typical output for two different cyclic loading path.  Table 2.

Discussion and Conclusions
The new model is based on physical rubber elasticity considerations. The total stress is the weighted sum of two contributions: an amplified, simplified extended tube model, and a hysteretic part. The hysteretic part was implemented using a single Maxwell/Prony element. For compounds with different CB volume fraction the model reproduces the non-linear behaviour under examination: Mullins effect, cyclic stress relaxation, permanent set, and hysteresis. The model requires 10 material parameters. The best-fit sets of parameters show a linear or a polynomial trend with the CB volume fraction. The sensitivity study shows the effect of altering each of the parameters. The FE implementation of the new model reproduces the expected behaviour. Finite Element Analysis was conducted on single-element models and on a realistic component. This preliminary investigation demonstrates that this model is stable when incorporated into a finite element model. The VUMAT works properly with different structures, and it is able to predict the desired non-linear viscoelastic behaviour. The new model was already implemented in the Jaguar Land Rover (JLR) suspension tool. The first job carried out a single cycle at 20 kN on a suspension loaded radially. Preliminary results are reported in Figure 13. Further studies could focus on the modeling of a broader range of deformation rates potentially requiring the use of more Prony elements. Moreover, the idea of a static hysteresis, as implemented using an intrinsic time in [28], deserves further investigation. Author Contributions: F.C. developed the curve fitting algorithm, coding and implementation and performed simulation. J.P. developed the non linear elastic model. F.C. and J.P. prepared the first draft of the manuscript. R.W. sponsored the project. J.B. and M.K. supervised the project throughout and contributed to detailed editing of the manuscript. All authors have read and approved the final version of the manuscript.

Funding:
The project was funded by Jaguar Land Rover, Banbury Road, Gaydon CV35 0RR, UK.