In a previous work, the experimental viscoelastic properties of a renewed rubber obtained from the recycling of end-of-life tires (ELTs) were presented, [1]. The material was tested in the laboratory by means of impact tests that revealed how the behavior of the rubber depends on the strain rate in such a case. Due to the loading conditions of high strains and high strain rates, the nonlinear viscoelastic Bergtröm- Boyce (BB) model, deﬁned by nine parameters, was selected to ﬁt the behavior of the material. To obtain the parameters from the results of the impact tests, diﬀerent optimization methods were used. As many sets of nine parameters were obtained as optimization methods were considered, with notable diﬀerences between them in some cases. The motivation of this work is to know which set of parameters, i.e. which optimization method, is the most appropriate. With this objective, a complete comparative numerical analysis of these sets, forty-eight in total, has been developed by means of simulations performed in LS-DYNA software. As a result of this analysis, the model that best suits the behavior of the material for the load case studied was chosen. This work completes the ﬁrst step to reach the ﬁnal objective of our research, which is to numerically analyze the behavior of systems designed to absorb energy in impact events based on this material. Keywords: End-of-life tires (ELTs), Rubber numerical analysis, Nonlinear viscoelastic models, Numerical rubber model validation, Bergtröm-Boyce (BB) model. 1 Introduction End-of-life tires (ELTs from here on) are those tires that have reached the end of their useful life, becom- ing, since this moment, a big environmental problem. The reason for that environmental problem is its main component: rubber designed to withstand the most extreme conditions. Consequently, this rubber is reduced, after its useful life, to a waste diﬃcult to eliminate. Currently, the two main ways to manage this problem are: (i) using this waste as fuel due to its high caloriﬁc value or (ii) to reuse this material either for a diﬀerent purpose, as raw material for other products or in any other way. 1 The numerical analysis of this type of materials is of great interest since, in this way, it is possible to study their applicability to diﬀerent uses. For this reason, there are numerous recent studies on the matter. Some of them are listed below. Pan et al. [2] developed a mesoscale model of rubberised concrete with diﬀerent crumb rubber contents up to 30% to evaluate the dynamic characteristics of the rubberised concrete. The existing experimental data proved the accuracy of the numerical model. Gao et al. [3] proposed a new form of double layer rubber dam to improve the water-retaining ability of conventional rubber dam. Numerical studies using FLAC2Dsoftware were carried out to analyze the behavior of the proposed method. The numerical model was veriﬁed by the results from the laboratory model tests and analytical studies in literature. Dhir et al. [4] have published a study about a novel mesoscale approach for simulating the response of reinforced concrete frames equipped with masonry inﬁll walls and rubber joints under horizontal loads such as those induced by earthquakes. Yang et al. [5] developed an aircraft tire and grooved runway models according to ﬁnite element (FE) simulations. The interaction model for the aircraft tire and grooved runway was developed using the Power Spectrum Density (PSD) and the viscoelastic properties of tread rubber. The FE model was validated for skid resistance by conducting runway ﬁeld experiments. Leng et al. [6] presented an experimental and numerical method on stress relaxation and unrecoverable damage characteristics of rubber materials. Liu et al. [7] carried out a study focused on the heterogeneous deformation and heat build-up of the crack tip region in a rubber fatigue crack growth rate (FCGR) test. Zhuo et al.[8] made numerical simulations and model tests conducted to study the inﬂuence of rubber fraction by volume (RF) on the macro/micro mechanical characteristic of lightweight ﬁll material (LFM) and its engineering applications. The inﬂuence of LFM buried depth on the deformation of cut-and-cover tunnel (CCT), earth pressure around CCT, and vertical displacement of backﬁll sand above CCT were also discussed. Šulc et al. [9] presented a new proposal for modeling the damping of hard rubbers at large deformations and quasi-harmonic excitation in the low frequency range. A torsion stand was designed and manufactured for rubber testing. According to the results of the experiments, the hyperelastic proportional damping (HPD) model based on the generalized Rivlin hyperelastic material model is suﬃciently accurate and, in addition, is easy to implement in a FE code. Nguyen et al. [10] developed analytical formulations and non-linear numerical models to study the be- haviour of natural rubber bearings incorporating U-shaped damper devices under cyclic loadings. Sheikhi et al [11] evaluated a natural rubber bearing system (NRB) equipped with a damper under lateral loading. U-shaped dampers with shape memory and steel alloy materials were used. The eﬀect of pa- rameters such as shear strain amplitude, the optimal thickness ratio of shaped memory alloy (SMA) to structural steel, NRB shape factor, thickness of SMA dampers, and the eﬀect of lateral loading directions on the performance characteristics of bearing system were investigated. Mooney-Rivlin and Prony models were selected for modeling the rubber material in ANSYS ®FE software. Modheg et al. [12] modeled a visco-hyperelastic damper (VHD) according to steel and rubber test results using high axial damping rubber in VHD under harmonic loading. With the broader aim of working towards the development of a class of comprehensive hyperelastic mod- els, Anssari et al. [13] utilised a three-parameter strain energy function for application to various types of isotropic rubber-like materials. In a recent work [1], the authors of the present one analyzed experimentally a new recycled rubber from ELTs produced by the Danish company genan ®. The recycled material, thanks to a re-vulcanization process, becomes a renewed rubber. Unfortunately, due to the manufacturing procedure its shape is restricted to small pellets, as it is shown in Figure 1. This limiting aspect causes its use needs the help of a resin acting as binder. In [1] the material was calibrated, considering the nonlinear viscoelastic Bergström-Boyce (BB) model, by mean of compression impact tests under diﬀerent impact velocities. Due to its agglomerated typical use, 2 Figure 1: Appearance of longest pellets manufactured from ELTs by genan ®. the tests were done over pellets in its original presentation and coated with a polyurethane resin. Through- out these tests, the energy dissipation capacity of this material was also evaluated. The nine material properties deﬁning the BB model were computed via several optimization techniques, applied to the ex- perimental records, using MCalibration ®software. As it is expectable, each optimization technique gives a diﬀerent set of material constants capable, in theory, to capture the viscoelastic behavior of the mate- rial. The present work focuses on the numerical analysis of the problem investigating how accurate are really those sets of constants by mean of their implementation in a FE model. With this aim the impact tests are simulated using the LS-DYNA ®software according two versions: •The native approach of LS-DYNA ®. •A modiﬁed approach (running under LS-DYNA ®) called PolyUMod ®speciﬁcally designed, among oth- ers, for the considered nonlinear viscoelastic model, developed by PolymerFEM ®. The following section, Section 2, presents the two BB model formulations considered in the present work. Section 3 shows the diﬀerent sets of constants obtained for both BB formulations using MCalibration ® software. As the work developed in [1] yielded to sets of material parameters according to only one of the formulations examined here, a new calibration procedure to obtain new parameter sets, coherent with the new formulation, has been done. This section also summarizes a comparative analysis of the diﬀerent sets of parameters obtained for both approaches. Section 4 is devoted to the numerical analysis. It describes the used model and the obtained results. Three cases of impact velocities are investigated, covering a strain rates range from 133.33s-1to400s-1. Finally, the manuscript ends with concluding remarks in Section 5. 2 Bergström-Boyce model formulations. Since its publication in 1998, [14], the non linear viscoelastic behavior deﬁned by the Bergstöm-Boyce model has been slightly modiﬁed to include more aspects of polymer’s behavior, making it more and more gen- eralist, see [15]. This model is based on the eight chains (EC) hyperelastic model published by Arruda and Boyce in 1993 to describe the elastic behavior of rubbers, [16]. In the present work two slightly diﬀerent approaches are considered for the model. Each one brings us to a diﬀerent solver running, both, under LS-DYNA ®software. One solver, in fact the native LS-DYNA ®solver, is based on the work published in 2009 by Dal and Kaliske, [17]. The other solver, called PolyUmod ®, has been developed by PolymerFEM ®. In the following sections both approaches are brieﬂy depicted. 2.1 PolyUMod ®formulation Conceptually, BB model is represented by two networks (or branches) acting in parallel (Fig. 2, left). One branch (let’s say branch a) contain only hyperelastic information while the other (let’s say branch b) cap- ture the non-linear viscoelastic behavior in series with hyperelastic information as well. The hyperelastic 3 behavior of each of these branches must be understood under the eight-chain macromolecular networks concept by Arruda and Boyce, [16]. Due to the parallel arrangement of branches aandb, the deformation gradient acting on the whole model is the same for both branches (Fig. 2, right). For the same reason, the total stress of the model is the sum of elastic stresses of each branch. Mathematically: F=Fa=Fb (1) σ=σa+σb (2) In eq. (1), Fstands for deformation gradient tensor and, in eq. (2), σdenotes the Cauchy stress tensor. Note: as usual, boldface has been used for tensors and normal text for scalars. Eight-chain model (hyperelastic) Viscoelastic flowEight-chain model (hyperelastic)Branch 𝑎: Equilibrium response Branch 𝑏: Time dependent responseReference configuration Relaxed configurationCurrent configuration𝐅=𝐅𝑎=𝐅𝑏 𝐅𝑏𝝂 𝐅𝑏𝒆 Figure 2: Conceptual representation of Bergström-Boyce model (left) and consequent multiplicative de- composition of deformation (right). In accordance with the constitutive equation of the eight-chain model, [16], the Cauchy stress tensors for the networks aandbhave the following form: σ=µ Jλ∗L−1 λ∗/λlock L−1(1/λlock)dev[b∗] +k(J−1)I (3) In eq. (3): •J=det(F), where det(·)stands for tensor’s determinant. •b∗is the distortional (constant volume) part of left Cuchy-Green tensor. •µis the initial shear modulus of the material. •kis the bulk modulus. •λ∗= [tr(b∗)/3]1/2, being tr(·)is the trace of the tensor. •Lis the Lanvegin function: L(x) =coth(x)−1/x. •Istands, as usual, for the identity matrix. Of course, the time dependence of the model must be in the viscoelastic branch b. The BB model decom- poses its deformation gradient tensor into the product of its elastic ( e) and viscoelastic ( v) components, as shown in eq. (4). Fb=Fe bFv b (4) In order to complete the model, see [15] for further details, it is necessary to consider the spatial velocity gradient of branch b(Lb), related with deformation gradient tensor as is indicated in eq. (5). Lb=˙FbF−1 b(5) Following [15], this tensor is decomposed in elastic and viscous components, eq. (6), Lb=Lbe+Dbv (6) 4 and the viscous component Dbvis expressed as in eq. (7). Dbv= ˙γ0 λbv−1 +ξC Rτ fvτbase−ˆτcutmdev(σB) τ(7) Where •R(x) = (x+|x|)/2, is the ramp function. •fvis anetwork interaction factor to capture the observed behavior for most elastomers having dif- ferent strain-rate dependence during loading than unloading. •τis the eﬀective stress driving the viscous ﬂow deﬁned as τ=||dev[σB]||, being || ||the Frobenius norm here and hereinafter. •˙γ0is a constant introduced for dimensional consistency called reference strain rate. • λbv−1 +ξCis a term introduced by reptational dynamics. This term captures a stretch depen- dence of the eﬀective viscosity, see [18]. •ξ,C,τbase,ˆτcutandm; are material constants. In particular, ξis small positive constant introduced to eliminate a removable singularity in the ﬂow rate in the undeformed state [15]. In summary, the constants set necessary to fully deﬁne a BB material following this approach, i.e., to solve the model using the PolyUMod ®solver, are the ones shown in table 1. Table 1: Material constants for PolyUMod ®solver, [19]. Mat. param. Name/Meaning µ[MPa ] Shear modulus of network a k[MPa ] Bulk modulus λL[−]∗Locking stretch s[−] Relative stiﬀness of network b ξ[−] Strain adjustment factor C[−] Strain exponential τbase[MPa ]Flow resistance m[−] Stress exponential ˆτcut[−] Normalized cut-oﬀ stress for ﬂow ∗Note: [−]stands for non-dimensional parameters. 2.2 LS-DYNA ®formulation The native formulation of LS-DYNA ®for BB model, based on [17] (of course, founded on the theory of paral- lel branches), follows the idea that the volumetric part is elastic and governed by the bulk modulus ( k) and the deviatoric Kirchhoﬀ stress tensor (τ)is the sum of an elastic (τe)and a viscoelastic part (τv). These deviatoric terms, following again the EC model in [16], are, in this case, written as following: τe=G 33−λ2 r 1−λ2r b−tr b 3I! , (8) τv=Gv 33−λ2 v 1−λ2v be−tr be 3I! , (9) being •b=J−2/3FFTwithJ=det(F) •be=FeFT ewithJ−1/3F=FeFi, where subindices e and i stands for elastic and inelastic parts, respectively. •GandGvare the elastic and viscoelastic shear modulus, respectively. •λ2 r=tr b /(3N)andλ2 v=tr be /(3Nv), where NandNvare the numbers of elastic and viscoelastic segments, respectively. 5 The time evolution is introduced in this approach via the inelastic rate-of-deformation tensor (Di)ex- pressed as in eq. (10). Di= ˙γ0(λi−0.999)c||τv|| ˆτ√ 2mτv ||τv||(10) Where •˙γ0is the reference strain rate. •ˆτis the reference Kirchhoﬀ stress. •λi=tr FT iFi /3. •candmare, respectively, the inelastic strain and stress exponents. As it has been stated above, the volumetric part of the model is elastic and governed by the bulk modulus (k). According to that, the pressure ( p) is given by: p= J−1−1 k. (11) In consequence, the constants set necessary to solve the BB model according to this formulation are sum- marized in table 2. Table 2: Material constants for LS-DYNA ®native solver, [20]. Mat. param. Name/Meaning k[MPa ] Elastic bulk modulus G[MPa ] Elastic shear modulus Gv[MPa ] Viscoelastic shear modulus N[−]∗Elastic segment number Nv[−] Viscoelastic segment number C[−] Inelastic strain exponent m[−] Inelastic stress exponent ˙γ0[1/s] Reference strain rate ˆτ[MPa ] Reference Kirchhoﬀ stress ∗Note: [−]stands for non-dimensional parameters. 3 Material calibrations Hereinafter, when necessary, the diﬀerent sets of constants will be denoted by LS-set and PUM-set for LS-DYNA ®and PolyUmod ®solvers, respectively. In the previous work, [1], the authors presented the experimental procedure to obtain the PUM-sets via compression impact test over original pellets, i.e., pellets as obtained from the recycled procedure, and over pellet with coating polyurethane resin 8wt%. In the present work, for comparison purposes, the same procedure has been used to obtain LS-sets. As stated above, software MCalibration ®, by PolymerFEM ®, has been used to perform the model calibra- tion. The utilized version (6.2.0) oﬀers twelve diﬀerent ways of calibration by combining three objective functions and four calibration algorithms, see [1] for a complete description of the calibration process. Tak- ing into account that we have two kind of specimens, coated (with polyurethane) and uncoated (in raw form) and two sets of parameters, LS and PUM sets, we have a total of forty-eight combinations. For the sake of completeness, all sets of experimental parameters have been summarized in the Appendix: tables 3, 4, 5 and 6. Looking those tables, some diﬀerences between properties, obtained from diﬀerent optimization proce- dures, appear to be very signiﬁcant. Looking for any conclusion about it, in the following, the results for the used calibration methods will be represented with respect to their averaged values. Trying to draw any trend regarding the optimization process, all sets of properties have been plotted in normalized way. With this aim, each parameter has been normalized by mean of the averaged values obtained using the four algorithms, for each objective function. Of course, this normalization procedure 6 has not physical meaning, but is able to show, at a glance, how diﬀerent are the results given for diﬀerent algorithms. The normalized parameters, denoted by and asterisk, are sumarized in Figs. 3 - 6. These plots have been grouped by objective function ( R2,MSD andNMAD ) and represented vs. the calibration algorithm ( GSM , EAM 1,EAM 2andQAM ). See [1] for more information about calibration algorithms and objective func- tions. 00,71,42,12,83,5 GSM EAM1 EAM2 QAM00,30,60,91,21,5 GSM EAM1 EAM2 QAM00,30,60,91,21,5 GSM EAM1 EAM2 QAM 00,30,60,91,21,5 GSM EAM1 EAM2 QAM00,61,21,82,43 GSM EAM1 EAM2 QAM00,40,81,21,62 GSM EAM1 EAM2 QAM 00,91,82,73,64,5 GSM EAM1 EAM2 QAM00,40,81,21,62 GSM EAM1 EAM2 QAM00,71,42,12,83,5 GSM EAM1 EAM2 QAM𝜇∗ 𝜆∗ 𝜅∗s∗ 𝜉∗ C∗ 𝜏𝑏𝑎𝑠𝑒∗ m∗ 𝜏𝑐𝑢𝑡∗ GSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMGSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMGSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMR2MSD NMAD Figure 3: Normalized PUM-set constants obtained for all the considered optimization procedures. Un- coated rubber. 00,61,21,82,43 GSM EAM1 EAM2 QAM00,511,522,5 GSM EAM1 EAM2 QAM00,30,60,91,21,5 GSM EAM1 EAM2 QAM 00,40,81,21,62 GSM EAM1 EAM2 QAM00,71,42,12,83,5 GSM EAM1 EAM2 QAM00,30,60,91,21,5 GSM EAM1 EAM2 QAM 00,91,82,73,64,5 GSM EAM1 EAM2 QAM00,30,60,91,21,5 GSM EAM1 EAM2 QAM00,81,62,43,24 GSM EAM1 EAM2 QAM𝜇∗ 𝜆∗ 𝜅∗s∗ 𝜉∗ C∗ 𝜏𝑏𝑎𝑠𝑒∗ m∗ 𝜏𝑐𝑢𝑡∗ GSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMGSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMGSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMR2MSD NMAD Figure 4: Normalized PUM-set constants obtained for all the considered optimization procedures. Coated rubber. 7 00,30,60,91,21,5 GSM EAM1 EAM2 QAM00,511,522,5 GSM EAM1 EAM2 QAM00,91,82,73,64,5 GSM EAM1 EAM2 QAM 00,511,522,5 GSM EAM1 EAM2 QAM00,511,522,5 GSM EAM1 EAM2 QAM00,40,81,21,62 GSM EAM1 EAM2 QAM 00,91,82,73,64,5 GSM EAM1 EAM2 QAM00,61,21,82,43 GSM EAM1 EAM2 QAM00,30,60,91,21,5 GSM EAM1 EAM2 QAM𝜅∗ 𝐺∗ 𝐺𝑣∗N∗ 𝑁𝑣∗ C∗ሶγ0∗ GSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMGSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMGSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMR2MSD NMAD m∗ Ƹ𝜏∗Figure 5: Normalized LS-set constants obtained for all the considered optimization procedures. Uncoated rubber. 00,30,60,91,21,5 GSM EAM1 EAM2 QAM00,61,21,82,43 GSM EAM1 EAM2 QAM00,61,21,82,43 GSM EAM1 EAM2 QAM 00,40,81,21,62 GSM EAM1 EAM2 QAM00,61,21,82,43 GSM EAM1 EAM2 QAM00,40,81,21,62 GSM EAM1 EAM2 QAM 00,511,522,5 GSM EAM1 EAM2 QAM00,40,81,21,62 GSM EAM1 EAM2 QAM00,511,522,5 GSM EAM1 EAM2 QAM𝜅∗ 𝐺∗ 𝐺𝑣∗N∗ 𝑁𝑣∗ C∗ሶγ0∗ GSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMGSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMGSMEAM1EAM2QAM GSMEAM1EAM2QAM GSMEAM1EAM2QAMR2MSD NMAD m∗ Ƹ𝜏∗ Figure 6: Normalized LS-set constants obtained for all the considered optimization procedures. Coated rubber. As it can be observed, exceptionally, the value of any of the materials parameters remain quite stable, i.e., its normalized value remain near to one, for example µ∗orm∗in Fig. 3. However, most of the values show big (and random) variations that can reach values 300% or even more: ξ∗in Figs. 3 or 4, κ∗and N∗ νin Figs. 5 and 6, respectively, among many others. Consequently, it is clear that no pattern can be observed at all. Taking this into account and on the basis of the numerical results that will be shown latter on, the conclusion is that attempting to analyze, individually, each parameter is a mistake. Every group of parameters has to be understood as a whole. So that, two groups of very diﬀerent parameters can give rise to a very similar viscoelastic behavior, as will be seen in the next section. 8 4 Numerical study In this section the numerical results obtained for all material properties summarized in the Appendix are presented. The compression impact tests in [1] were applied to rubber cylinders by an aluminum-made cylindrical hammer, see Figure 7 left. To record the reaction during the impact, the specimens were placed over a load cell, see Figure 7 right. Figure 7: Specimen over load cell (left) and hammer over the specimen (right). Obviously, its numerical implementation does not oﬀer any geometric diﬃculty. As it has been stated, for impact test simulation LS-DYNA ®software has been selected. The calculations are made by the explicit method. For symmetry reasons, a quarter of the geometry has been modeled as it can be see in Figure 8. Although the problem posses axial symmetry, no 2D modelization was possible because it is not imple- mented for BB material model. Figure 8: One-quarter pellet-hammer-rigid plane model. Figure 8 shows, the cylinder (BB material) stays over a rigid surface playing the role of the load cell, i.e., the reaction over this rigid surface can be compared with the results obtained from the load cell during the tests. The hammer is deﬁned by another cylindrical body with aluminum elastic properties and a ﬁctitious value of density in order to reproduce the total falling mass. For all contacts, sliding hypothesis have been considered. The computation starts just in the contact mo- ment, therefore the nominal impact velocities considered, i.e. 2,4and6m/s (minimum, medium and maximum test velocities in [1]) are established as initial conditions. The nominal dimensions of the rubber cylinder are 4.2mm of diameter and 15mm long. 9 4.1 Convergence analysis The problem is strongly non linear because contact, non-linear material, high impact velocities and large strains are under consideration. The necessary convergence analysis must take into account spacial ( ∆L) and temporal ( ∆t) discretizations. At any event, it must be considered that LS-DYNA ®assumes ( ∆t) as initial time step and changes it as needed. For this purpose, we have compared the results obtained with three diﬀerent models. The considered models were deﬁned as: •Model A: ∆L = 0 .5mm,∆t= 0.01ms. •Model B: ∆L = 0 .2mm,∆t= 0.01ms. •Model C: ∆L = 0 .2mm,∆t= 0.001ms. For this analysis the same material constants set has been considered to avoid the inﬂuence of these pa- rameters over the convergence conclusions. Figure 9 shows that the engineering stress vs. engineering strain results are, essentially, the same for all the analyzed models. In such a situation, we must conclude that the model is stable enough regarding spacial and temporal discretization. Although computational time is not a strong limitation, since the computation time is kept in the order of minutes, model A is about 50 times faster than B or C. For all the indicates considerations, the model A has been ﬁxed for the rest of the computations of the work. Figure 9: Convergence analysis regarding time and spacial discretizations. 4.2 Numerical results Figures 10 to 13 show the numerical results obtained (together the experimental results, in grey color) for uncoated and coated pellets, for the two BB formulations investigated. To show the inﬂuence of the calibration method, the plots have been grouped according to the calibration algorithm. As it can be seen, an excellent ﬁt has been obtained for all LS-sets and for all strain rates (Figures 10 to 12). These results show that all the calibration procedures have led to parameter sets that faithfully reproduce the behavior of the material, regardless the objective function and/or optimization algorithm used. As indicated above, this suggests that the separate consideration of the behavior parameters, for comparison purposes, for instance, makes not sense. They must be considered as part of a whole. This "whole" is a closed group of parameters that bring us to a speciﬁc viscoelastic material behavior. On the contrary, the PUM-sets (Figure 13) have shown, for a strain rate of 400s-1, much more rigid results than the experimental ones, except for the GSM algorithm whose results are good enough. In such a scenario, we have decide do not use PUM-sets for the rest of velocities considering the dispersion of the results for that strain rate. 10 Figure 10: Numerical results, via LS-DYNA ®native solver, for coated and uncoated LS-sets, considering a strain rate of 133.33s−1. 11 Figure 11: Numerical results, via LS-DYNA ®native solver, for coated and uncoated LS-sets, considering a strain rate of 266.66s−1. 12 Figure 12: Numerical results, via LS-DYNA ®native solver, for coated and uncoated LS-sets, considering a strain rate of 400s−1. 13 Figure 13: Numerical results, via PolyUmod ®solver, for coated and uncoated PUM-sets, considering a strain rate of 400s−1. 14 5 Conclusions The sets (forty-eight in total) of material parameters, deﬁning a viscoelastic Bergström-Boyce behavior model for a renewed rubber from end of life tires, obtained in [1] by mean of an experimental procedure, have been compared from a numerical point of view. To this aim, the software LS-DYNA ®has been used to simulate compression impact tests applying impact velocities of 2,4and6m/s; that means strain rates of 133.33,266.66, and 400s−1, respectively. Two diﬀerent solvers have been tested with this software using slightly diﬀerent approaches of the Bergström- Boyce model: (i) the native LS-DYNA ®solver based on the work published in 2009 by Dal and Kaliske [17] and the solver developed by PolymerFEM ®called PolyUmod ®, [19], based on Bergström formulation [21]. Our results indicate that, using the native solver of LS-DYNA ®, all the material parameters sets give us material behaviors ﬁtting perfectly with experimental results for all considered strain rates. On the contrary, more rigid material behaviors are found using the PolyUmod ®solver by PolymerFEM ®, except for GSM optimization algorithm that give, as well, a material behavior ﬁtting well enough whit our experimental ﬁndings. In summary, in a frame where a strong limiting shape of specimens does not allow to make more than one kind of test, we have demonstrated that the compression impact test, followed an optimization procedure (for instance, the ones include in MCalibration ®software), can give us very good results to obtain the nine parameters necessary to deﬁne the BB viscoelastic model of a renewed recycled. That is, a useful tool has been developed to facilitate the implementation of this kind of materials in a ﬁnite element code for further numerical computations. 6 Acknowledgments The authors want to show their gratefully to: •Consejería de Economía y Conocimiento of Junta de Andalucía and the European Regional Develop- ment Fund through the project UMA20-FEDERJA-028 for its ﬁnancial support. •Spanish Ministry of Financial and Competitiveness for the funds granted by mean of the projects with refs. UNMA13-1E-1901 and UNMA15-CE-3571. These grants allowed to acquire the experimental equipment used in this work.