A Study of the Thermal Runaway Phase Generated During Cyclic Loading of Viscoelastic Materials Accounting for the Prestress

THE GOOD DAMPING PERFORMANCE AND INHERENT STABILITY OF VISCOELASTIC MATERIALS IN RELATIVELY BROAD FREQUENCY BANDS, BESIDES COST EFFECTIVE-NESS, OFFERS MANY POSSIBILITIES FOR PRACTICAL ENGINEERING APPLICATIONS. HOWEVER, FOR VISCOELASTIC DAMPERS SUBJECTED TO DYNAMIC LOADINGS SUPERIMPOSED ON STATIC PRELOADS, ESPECIALLY WHEN GOOD ISOLATION CHA-RACTERISTICS ARE REQUIRED AT HIGH FREQUENCIES, TRADITIONAL DESIGN GUIDELI-NES CAN LEAD TO POOR DESIGNS DUE TO THE RAPIDLY INCREASING RATE OF TEM-PERATURE CHANGE INSIDE THE MATERIAL. THIS PAPER IS DEVOTED TO THE NU-MERICAL AND EXPERIMENTAL INVESTIGATION IN THE DEGRADATION OF THE STIFFNESS AND CAPACITY OF A VISCOELASTIC MATERIAL INDUCED BY THE THERMAL RUNAWAY PHASE, WHEN IT IS SUBJECTED TO DYNAMIC AND STATIC LOADS SIMUL-TANEOUSLY. AFTER THE THEORETICAL BACKGROUND, THE OBTAINED RESULTS IN TERMS OF THE TEMPERATURE EVOLUTIONS AT DIFFERENT POINTS WITHIN THE VOLUME OF THE MATERIAL AND THE HYSTERESIS LOOPS FOR VARIOUS STATIC PRE-LOADS ARE COMPARED AND THE MAIN FEATURES OF THE PROPOSED STUDY ARE HIGHLIGHTED.


INTRODUCTION
As compared to other control strategies, the viscoelastic materials present some important advantages such as inherent stability, effectiveness in broad frequency bands and moderate development and maintenance costs (Nashif et al., 1985;Rao, 2001;Samali and Kwok, 1995). As a result, they have been incorporated in many engineering applications such as tall buildings and aircraft and spacecraft structures in various types of configurations (Rao, 2001;Samali and Kwok, 1995;Drake and Soovere, 1984;Merlette, 2005). Although better damping performance is obtained by using traditional design procedures of viscoelastic dampers by assuming uniform and timeindependent temperature distribution within the material, whose value is generally chosen to be coincident with the ambient temperature under operation (Cazenove et al., 2012;Schapery, 1964;Palmeri and Ricciardelli, 2006;Lesieutre and Govindswamy, 1996), such strategies do not take into account the thermal runaway phase.
Depending on the applied loading conditions, the vibration energy of the viscoelastic material is converted to heat at a rate faster than the heat is conducted away, leading to a rapidly increasing rate of local temperature change known as thermal runaway phase. As a result, it affects significalntly the stiffness and damping proprerties of the viscoelastic materials (Lesieutre and Govindswamy, 1996), especially in applications in which the viscoelastic materials are subjected to dynamic loadings superimposed on static preloads, such as engine mounts and tall buildings (Palmeri and Ricciardelli, 2006).
As discussed by Koyalenko and Karnaukhov (1972), in the context of the heat transfer analysis of viscoelastic polymers, the interrelation between thermal and mechanical quantities leads to a highly nonlinear coupled problem for which exact solutions cannot be easily obtained. Thus, approximate numerical solutions must be searched (Schapery, 1964). The authors in Cazenove et al. (2012) have proposed a first hybrid numerical-experimental strategy to investigate the influence of frequency and amplitude of the excitation on the self-heating phenomenon in a simple two-dimensional viscoelastic damper subjected to a pure shear deformation. They observed that, as the frequency and the amplitude of the cyclic loading increase, the temperature inside the viscoelastic media increases significantly, leading to a great reduction in the damping capacity of the damper. However, this early work was limited to two-dimensional stress states and the effects of the static strains on the self-heating nor the runaway phase have not been addressed. Gopalakrishna and Lai (1998) have concluded that the local temperature values inside the viscoelastic materials normally used in tall buildings as passive control strategy can increase up to 10°C in a few seconds during a storm. However, this work was also limited to two-dimensional stress states and the prestressed has not been investigated. Rittel (1999) has demonstrated experimentally that the conversation ratio used to measure the mechanical energy transformed into heat for polycarbonates deformed in a range of strain rates is strongly dependent on strain and strain rate. Also, as these later are very difficult to estimate, their values have been assumed rather arbitrarily in the range [0.1; 1.0], such as in the study reported by Rittel and Rubin (2000) in which it is performed thermal analyses of the cyclic compression of Polymethylmethacrylate (PMMA) and polycarbonate cylindrical specimens.
More recently, de Lima et al. (2015) have suggested a three-dimensional numerical approach for describing the self-heating in viscoelastic dampers accounting for the complex state of stress, but nothing was reported about the thermal runaway phase nor the combined effects of the dynamic Latin American Journal of Solids and Structures 13 (2016) 2834-2851 and static strains on it. Lesieutre and Govindswamy (1996) have used the Anelastic Displacement Field (ADF) model incorporating the temperature dependence for thermorheologically-simple viscoelastic materials in order to capture the essential features of the mechanical and thermal responses for a simple one-dimensional viscoelastic device in simple shear. Rodas et al. (2014) have used a three-dimensional large strain thermo-viscoelastic constitutive model to describe the self-heating of rubber materials during low-cycle fatigue response, using a Zener rheological representation for the mechanical part. However, the combined effects of static and dynamic loadings and the thermal runaway phase have not been addressed by the authors. Mortazavian et al. (2015) have performed an interesting experimental study to evaluate the effects of frequency and self-heating on fatigue life of reinforced and unreinforced thermo-plastic polymers. Braeck and Podladchikov (2007) have performed a theoretical investigation of the thermal runaway failure mechanism of viscoelastic materials without considering the prestressed. In a subsequent work (Braeck and Podladchikov, 2009), the authors have suggested a model based on a continuum formulation of viscoelastic materials with Arrhenius dependence of viscosity on temperature to investigate the thermal runaway. This is thought to be a limiting assumption as it has not been considered the complete three-dimensional strain state to compute the heat source nor the effects of the prestressed on the thermal runaway phenomenon. Jue et al. (2006) have used the Split Hopkinson pressure bar thecnique to characterize the stress uniformity for viscoelastic materials. However, they do not investigate the influence of the temperature on the stress distribution along the thickness of the material. Neto and Rosa (2008) have suggested an interesting strategy to incorpórate metrological aspects to characteize the input parameters uncertainty of a viscoelastic model through measument processes. Again, nothing was reported about the self-heating phenomenon.
Although heat transfer analyses of elastomers has been already conducted by some other authors (Katunin, 2012;Katunin and Gnatowski, 2012;Kirichok, 2013;Lesieutre and Mingori, 1990), the contribution intended for the present study is to investigate the degradation in stiffness and damping capacity of viscoelastic materials due to the thermal runaway phase induced by cyclic loadings superimposed on static preloads. Also, by comparing the obtained results with those appearing in the open literature, the following features are to be highlighted: (i) the thermal-viscoelasticmechanical modeling as presented here can be applied to more complex viscoelastic systems due to the formulation of a complete strain state; (ii) the effects of the static preloads on the thermal runaway phase and its effects on the stiffness degradation and damping capacity of the viscoelastic material is also considered as contribution; (iii) in view of practical applications by using the available commercial FE codes from the engineering perspective, it was performed a harmonic prestressed analyses to take into account the static preloads superimposed on the dynamic excitations; (iv) an optimization problem has been formulated in order to identify the unknown thermal parameters describing the experimental scenarios.

Modeling of Viscoelastic Systems Using Commercial FE Codes
It is well-kwon that the mechanical behavior of viscoelastic materials depends strongly on certain environmental and operation factors, among which the most relevant are the excitation frequency and the temperature (Nashif et al., 1985;Christensen, 1982). In this context, some mathematical models have been developed to represent this behavior and have been shown to be particularly suitable to be used in combination with FE discretization. Among those models, it should be mentioned the so-called Fractional Derivative Model (FDM) (Bagley and Torvik, 1979;Bagley and Torvik, 1983), the Golla-Hughes-McTavish Model (GHM) (Golla dn Hughes, 1985;MacTavish and Hughes, 1993), the Anelastic Displacement Fields Model (ADF) suggested by Lesieutre and Bianchini (1995) and the use of Prony series in order to characterize the viscoelastic relaxation modulus (Pacheco et al., 2014). All these models can represent the viscoelastic behavior both in time and frequency domain. However, when the viscoelastic materials are subjected to dynamic loads superimposed on static preloads, it is necessary to evaluate the Complex Modulus to describe this combined effect. Thus, based on the fact that the static and dynamic properties of linear viscoelastic materials can be measured independently of each other (Nashif et al., 1985), and taking into account the Frequency-Temperature Superposition Principle (FTSP), the Complex Modulus can be expressed as (de Lima et al., 2015): is the reduced frequency,  is the actual excitation frequency, T indicates the temperature, 0 T is a reference value of temperature,  is the static strains, and   are the storage modulus, loss modulus and loss factor, respectively.  (Christensen, 1982), which establishes an equivalence between the excitation frequency,  , and temperature, T , on the viscoelastic material properties. This implies that the viscoelastic behavior at different temperatures can be related to each other by shifts in the actual values of the frequency. Figure 1 illustrates the FTSP, showing that having the modulus and loss factor of an arbitrary viscoelastic material for different temperature values, 1  T , 0 T , 1 T , if horizontal shifts along the frequency axis are applied to each of these curves, all of them can be combined into a single one, named master curve (Nashif et al., 1985). The horizontal shift is given by and depends on the temperatura. Drake and Soovere (1984) suggest the analytical expression (2) for the complex modulus for the 3M™ ISD112 viscoelastic material defined in the following temperature and frequency intervals K T 360 210   and Hz From the reduced temperature nomogram generated by the computation of Eq. (2), it is possible to obtain the complex modulus and the loss factor at any given temperature into a frequency band of interest. However, the derivation of    F is more complicated and the difficulty in predicting general accurate analytical expressions comes from the fact that it could take several forms depending on the relation between the static preload,  , and engineering strain,  , involved in the stress-strain relationship for a viscoelastic device. Moreover, a supplementary difficulty arises from the large number of combinations of dynamic and static loading parameters involved in the experimental characterization of viscoelastic materials (Mooney, 1940;Rivlin, 1947). Hence, in the present study, the dynamic viscoelastic material data as a function of frequency and temperature are provided by Eq. (2), and the static strains are directly applied in the viscoelastic damper by performing a prestressed harmonic analysis using the ANSYS FE sofware (2015) in order to compute the dynamic responses of the prestressed viscoelastic damper. However, to make possible the modeling of the viscoelastic behavior using this FE code, some modifications in the equations of motion of the viscoelastic system must be done before.
From the equations of motion for a viscoelastic system in the time-domain (de Lima et al., is the stiffness matrix related to the purely elastic part of the structure, is the frequency-and temperature-dependent viscoelastic stiffness matrix, and is the mass matrix. N is the number of degrees-of-freedom (DOF's). Thus, taking into account the real and imaginary parts of the complex modulus function (2), the dynamic stiffness matrix can be rewritten as: Equation (3) shows that the viscoelastic behavior can be easily introduced into commercial FE codes based on the assumption of a viscous damping matrix proportional to the stiffness matrix, with a frequency-and temperature-dependent coefficient given by . This strategy is considered here to perform the prestressed harmonic analyses of the viscoelastic damper.

Heat Generation Rate for Viscoelastic Materials
According to Rittel and Rubin (2000), for a viscoelastic material, the basic transient heat balance equation, expressing the relationship between a heat generation term and the spatial and temporal variations of the temperature field is given by: where k is the thermal conductivity,  is the density, p c is the specific heat and T is the temperature.
If the material is subjected to cyclic stresses superimposed on static preloads,  , the heat generation term assumes the form (de Lima et al., 2015), the thermal conversion ratio defined as a fraction of the mechanical power dissipated by the material (Rittel, 1999) and C is the frequency-and temperature-independent material properties matrix for which the complex modulus,   T G ,  , has been factored-out. By integrating g q over a cycle of vibration for a sinusoidal strain response it is possible to demonstrate that the contribution of the real part of the complex modulus to the dissipated mechanical energy (associated to the elastic power stored in the material) vanishes. However, the imaginary part associated to the viscous power leads to the following average value of the heat generation (per volume): where   and 0  represent, respectively, the amplitudes of the applied static and dynamic strains.
Hence, in the context of the FE modeling, the heat generation term for a viscoelastic element can be obtained by integrating Eq. (5) over the volume of it: is the real part of the viscoelastic stiffness matrix, B is the matrix formed by the differential operators appearing in the strain-displacement relations ( lastic material due to the dynamic loadings accounting for the prestress. However, as demonstrated by Koyalenko and Karnaukhov (1972), the exact solutions of the resulting nonlinear thermoviscoelastic problem cannot be obtained easily. As a result, an iterative resolution scheme has been implemented here by the proper inclusion of the strain components into Eq. (6) in the ANSYS FE code.

ITERATIVE PROCEDURE
Since direct coupling between thermal and structural fields would result in prohibitive computational costs as discussed in de Lima et al. (2015), the problem was solved by assuming weak coupling between thermal and mechanical fields by using a sequential iterative scheme implemented in ANSYS FE software. The main steps of the procedure are summarized as follows: 1. firstly, the viscoelastic stiffness, the ambient temperature, respectively; 6. finally, the viscoelastic material properties are updated and if the convergence criterion based on temperature variations between two consecutive iterations is not satisfied within a specified tolerance, a new iteration is initiated and the structural and thermal analysis are performed based the updated viscoelastic stiffness and equivalent viscous damping matrices.

NUMERICAL RESULTS
To verify the thermal runaway phase and its strongly influence on the stiffness degradation and damping capacity of a viscoelastic material subjected simultaneously to dynamic and static loadings, numerical simulations were performed using the viscoelastic damper illustrated in Fig. 2. It consists of two viscoelastic layers made of 3M  ISD112 (2015) viscoelastic material inserted between three steel blocks. The associated FE mesh shown only for one-half of the viscoelastic damper is depicted in the same figure. Also, it was assumed the same thickness for both the viscoelastic cores and the steel blocks in z-direction of 30mm. In the discretization procedure, it was used the following three-dimensional coupled mechanical and thermal filed elements SOLID45 and SOLID70, respectively. For the thermal boundary condition, it is assumed a heat transfer by natural convection between the surfaces of the outer-steel plates and the surrounding air with, . Also, for the thermal runaway phase, it is not possible to identify a progressive stabilization of the temperatures in the loading phase and although not clearly shown in Fig. 3, the temperature does approach a higher equilibrium temperature value. In this application, the cyclic displacement was superimposed on different values of static preloads according to the following scenarios: (1) without preload; . By analyzing the temperature distributions for points A, B and C for all test scenarios shown in Fig. 4, it is possible to note immediately that the temperature inside the viscoelastic core is not constant and changes from one point to the other, owing to the different heat transfer conditions to which each point is subjected (conduction and natural convection). Also, as the static preload increases, the thermal runaway characteristics become more pronounced, as shown in Fig. 5, which enables to compare the temperature evolutions at point A for different values of static preloads. As a result, it is observed a significantly increasing in the temperature of the viscoelastic material. Thus, it is expected that such levels of temperature variations will influence significantly the stiffness and damping capacity of the viscoelastic material or even its complete failure in practical applications.   (1) and (3). It can be clearly perceived that the temperature gradient is oriented towards the center of the viscoelastic core, but due to the fact that the thermal boundary conditions are not symmetric and the low thermal conductivity of the viscoelastic material, the zone of the higher temperatures is concentrated exclusively in the viscoelastic core. Again, it is possible to observe that the temperature distribution changes significantly with the applied static preloads in contrast with the assumption of uniform temperature within the volume of the material, which is frequently found in the analysis and design of viscoelastic dampers. scenario (1) scenario (3) . 1207 s t  Figure 7 shows the dissipated, W , and storage, p W , energies of the viscoelastic material as functions of the number of cycles and Fig. 8 represents the storage modulus and loss factor as functions of time for all the test scenarios. As the static preload incresases, the thermal runaway characteristics is more pronounced, resulting in a great influence on the viscoelastic material properties, since it is observed that the temperature rise leads to maximum decreases of the loss factor of 22% for the thermal equilibrium and of 83.90%, 84.40% and 88.58% for the three thermal runaway situations, respectively. Clearly, variations of such magnitudes can lead to significant decrease of the viscoelastic damping capacity in practical applications.  By comparing the stress-strain trajectory (hysteresis loops) shown in Fig. 9 for all the test scenarios, it can be noted that a significant reduction in the amplitudes of the tip displacements of the viscoelastic damper shown in Fig. 2 is achieved when the thermal runaway phase becomes more pronounced. This is due to the stiffness degration and the reduction in the damping capacity of the viscoelastic material for high temperature values, resulting in an augmentation of the transverse displacement of it. Moreover, the changes in slopes and shapes of these trajectories indicate the occurrence of initial instabilities due to transition from stiff, loss behavior to softer (Lesieutre and Govindswamy, 1996), demonstrating the great influence of the thermal runaway phase on the behavior of the viscoelastic damper. However, although initial instabilities appear in both test scenarios investigated herein, the trajectories converge to a single closed curve, which indicates that the response of the system is stable.

EXPERIMENTAL INVESTIGATION
This section describes the experimental setup used to evaluate experimentally the influence of the cyclic loading combined with static preloads on the thermal runaway phase in viscoelastic materials. The viscoelastic damper was constructed according to the geometrical dimensions depicted in Fig. 2, composed by two 5mm-thick layers of 3M  ISD112 rubber-like material between three rigid steel blocks attached to a rigid frame mounted on a universal test machine MTS-800, as depicted in Fig.  10. The temperatures at points A and C of the specimen were measured by using thermocouples and its signal was acquired and processed by a signal analyzer Agilent  34970. Also, in all tests, a vertical displacement,     To compare the simulations obtained previously with the acquired experimental results, the static displacement of the viscoelastic material was computed based on the tangent stiffness matrix concept (de Lima et al., 2015a). Also, it is assumed that only the viscoelastic cores are deformed during the application of the static displacements by the screws. Thus, the static displacements to be applied on the specimen is, is the cross-sectional area of the viscoelastic core, MPa G 430700 0  designates the static value of the modulus at a low frequency range, is the thickness of viscoelastic layer, and  denotes the value of the static preloads used in the simulations. It enables to define the following test scenarios to be investigated: (a) At this time, it is important to emphasize that the effects of pre-cycling the viscoelastic damper prior to experiments have not been performed. However, the experimental setup has been reset from one test to another by alleviating the screws used to apply the static displacements. Also, an ideal period of time has been respected accounting for the creep phenomenon normally exhibited by the viscoelastic materials in order to guarantee that the computed changes in the temperature evolutions due to thermal runaway phase are mainly due to the static strain variations on the viscoelastic material. Figure 11 presents the temperature evolutions obtained experimentally for the scenario (1), demonstrating a first phase in which there is a rapid temperature increase; and the subsequent behavior characterized by a progressive stabilization of the temperature field which indicates a tendency of convergence for the temperature of the material. The experimental observations are in agreement with those observed in simulations.
The results appearing in Fig. 12 confirm experimentally the influence of the static preload on the thermal runaway phase and the difference of the temperature distribution within the viscoelastic material from one point to the other, as verified in the simulations. Again, the qualitative comparison between the experimental and numerical observations is satisfactory, but the direct comparison of it is not possible, as shown in Fig. 13, since for the simulations the values of the thermal conversion ratio,  , and the natural convection coefficient, h , have been assumed arbitrarily.
Cleary, it should be reminded that these parameters vary significantly with the geometry of the Latin American Journal of Solids and Structures 13 (2016) 2834-2851 damper and the loading conditions (Rittel and Rubin, 2000). Thus, it is expected that the values of  and h used in the simulations do not represent the experimental conditions, and it was found to be interesting to develop a curve-fitting procedure to identify these parameters for each experiment.   The experiments were performed in similar environmental conditions and it results to the same natural convection coefficient for all the tests. However, since the thermal conversion ratio is strongly dependent on the loading conditions, its valued vary according to the test scenario. Thus, an optimization problem was formulated here in order to fit the thermophysical parameters: (i) first, at point A, for test scenario (b),  and h are identified simultaneously by using the Non-dominated Sorting Genetic Algorithm (Srinivas and Deb, 1993), by minimizing the function, The total number of generations of the NSGA was limited to 10, which means that the maximum number of evaluations of the cost function was 300. Figure 14 compares the temperature evolutions predicted by the adjusted model using the identified values and the experiments for point A, considering tests scenarios (a) and (b). For all the test scenarios, correlations between numerical and experimental measured temperatures are very satisfactory. Moreover, it can be seen that as the static preload increases, the thermal conversion ratio decreases, since a complementary part of the dissipated power is stored in the viscoelastic material through microstructural changes. As a result, as the static strains increase, the microstructural modifications increase accordingly, leading to a lower fraction of the dissipated energy that is converted into heat. These observations are qualitatively very consistent with those reported by Rittel and Rabin (2000).

CONCLUDING REMARKS
The effects of the thermal runaway phase in viscoelastic materials have been investigated by means of numerical simulations and experimental tests. The FE-based modelling strategy has been implemented in ANSYS code by using an iterative scheme to solve the system of equations for the thermal and structural fields. One of its advantages is the ability to incorporate the viscoelastic behaviour in commercial FE codes based on the concept of equivalent viscous damping matrix and the complex modulus function, whose values are provided by the material manufacturers, accounting for the frequency-and temperature-dependent behaviour of viscoelastic materials. Moreover, due to the difficulty in estimating some of the thermal properties, a curve fitting procedure based on the formulation of an optimization problem has been implemented.
The numerical and experimental results obtained confirm the strongly influence of the static preloads on the temperature distribution within the viscoelastic materials. Moreover, the results have demonstrated that the rapidly increasing in the temperature field in the viscoelastic media has the characteristics of the thermal runaway phase although it is observed that the temperature does approach a higher equilibrium value. Consequently, it is expected that this temperature increases can affect significantly the mechanical properties of the viscoelastic materials and jeopardize their damping capacity.
The relevance of the obtained results enables to conclude that the use of traditional design procedures of viscoelastic dampers for vibration mitigation, in which constant and uniform temperature distribution is assumed, can be inadequate in certain loading situations. As a result, it leads to an unexpected damping performance or even complete failure of the viscoelastic damper due to the stiffness degradation of the viscoelastic material induced by the thermal runaway phase.