**Anelastic Behavior in Filled Elastomers Under Harmonic Loading Using Distributed Rate-Dependent Elasto-Slide Elements**

Wei Hu and Norman M. Wereley

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/2784

### **1. Introduction**

16 Will-be-set-by-IN-TECH

We have obtained the circularly polarized reflectances and transmittances as function of the chiral order parameter of a cholesteric elastomer immersed in a racemic solvent. We have found considerably changes in the bandwidth of the reflectance for left- and right-copolarized light under the presence of the solvent which are susceptible to be detected experimentally. We have obtained a *RRR* bandwidth larger than that for an undistorted elastomer (*α* = 0) for *α* = 0.4 and obliquely incidence for angles larger than 60◦. Also in this value we have observed a thin band reflection for both polarizations *RRR* and *RLL* for oblique incidence.

The modifications of the transmittance and reflectance spectra during a preferentially absorbing process suggest to utilize the optical spectra as an indirect method to determine

*Depto de Fisica, Universidad de Sonora, Apdo Post. 1626, Hermosillo, Son. Mexico 83 000, México*

*Depto. de Investigacion en Fisica, Universidad de Sonora, Hermosillo, Son., Mexico, Apdo P. 5-088,*

[1] A. Jakli, *et al*. *One and two dimensional fluids: Properties of smetics, lamellar and columnar liquid crystals (Condensed Matter Physics)*. Taylor and Francis: Florida, EU, (2006).

[4] M. Warner, *et al. Liquid Crystal Elastomers*. Oxford Science Publications: Oxford, EU,

[9] A. Lakhtakia,W.S. Weiglhofer. *Proc. R. Soc. Lond. A.,* 453, 93, (1997); correction: 454, 3275

the concentration of preferentially absorbed molecules during a segregation process.

**6. Conclusions**

**Author details**

Jesus Manzanares-Martinez

[2] B. Singh. *J. Phys. D: Appl. Phys*., 40, 584 (2007). [3] Y. Mao, M. Warner. *Phys. Rev. Lett*., 86, 5309, (2001).

[5] M. Rivera, J. A. Reyes*. Appl. Phys. Lett*., 90, 023513, (2007).

[8] Y. Mao,M. Warner. *Phys. Rev. Lett*., 84, 5335, (2000).

[10] N. Marcuvitz,J. Schwinger. *J. Appl. Phys*., 22, 806, (1951).

[6] P. Cicuta, A. R. Tajbakhsh, E.M Terentjev. *Phys. Rev. E.*, 70, 011703, (2004). [7] S. Courty, A. R.Tajbakhsh, E. M. Terentjev. *Phys. Rev. Lett.*, 91, 085503, (2003).

[11] P. Cicuta, A. R. Tajbakhsh. E. M. Terentjev. *Phys. Rev. E.*, 65, 051704, (2002).

P. Castro-Garay

*México 83 190.*

**7. References**

(2007).

(1998).

Elastomeric materials are broadly used for stiffness and damping augmentation in various applications due to their simple design, low weight, and high reliability. However, filled elastomeric materials exhibit significant nonlinear behavior [1]. Such nonlinear characteristics include a) non-elliptical shear strain vs. stress hysteresis diagram under sinusoidal excitation, b) stiffness and damping dependent on amplitude, frequency, temperature and even preload, c) Mullins effect with reduction of stiffness at small strain following cyclic deformation at large strains, and d) low stress relaxation and creep rates. Specifically, the large reduction of damping with increasing amplitude of harmonic displacement excitation leads to excessive size and weight of dampers in order to accommodate all operating conditions. It was also found that highly damped elastomeric dampers demonstrated low loss factors at low amplitudes resulting in unacceptable limit cycle oscillations [2]. Therefore, a precise analytical model is necessary to describe the nonlinear behavior of an elastomer and to determine its dynamic characteristics.

Some prior research introduced nonlinear terms into conventional Kelvin or Zener models. The complex modulus is usually used to characterize viscoelastic materials under harmonic excitation: the storage (or in-phase) modulus is a measure of the energy stored over a cycle or period of oscillation, and the loss (or quadrature) modulus is measure of the energy dissipated over a cycle. This model can be represented as a spring and a dashpot in parallel (Kelvin chain). However, the complex modulus is a linearization method in the frequency domain that represents the nonlinear hysteresis cycle as an equivalent ellipse, and is only applicable to steady harmonic forced response analysis. Some researchers have extended the basic Kelvin chain to complicated mechanism-based modeling approaches in order to explain the nonlinear behavior of elastomers. To display basic behavioral characteristics,

© 2012 Hu and Wereley, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

such as creep and relaxation, a viscoelastic solid can be represented as a spring in series with Kelvin elements (if one Kelvin element is used, then a Zener model results [3]). Gandhi and Chopra [4] developed a nonlinear viscoelastic solid model in which a nonlinear lead spring was used in series with a single linear Kelvin chain. Using this model, the variation of complex moduli with oscillation amplitude closely matched experimental data. Felker *et al..* [5] further developed a nonlinear complex modulus model based on a single Kelvin chain, in which the spring force was a nonlinear function of the displacement, and the damping force was a nonlinear function of displacement and velocity. This model was used to describe the amplitude dependent moduli and to study dual frequency damper motions. The parameters in all of these models were identified using amplitude-dependent complex modulus data. As a result, the nonlinear hysteresis behavior of could not be captured using these models.

Elastomeric materials typically demonstrate tribo-elastic behavior [1]. Physically, the amplitude dependent behavior exhibited by the elastomer is results from the interaction of the fillers with the rubber or elastomeric matrix materials [6]. Before large deformation of a filled elastomeric damper, an intact filler structure displays a large stiffness and low loss factor for small amplitudes. As the amplitude increases, the filler structure breaks resulting in a stiffness reduction. However, the breaking of filler structures, which is similar to friction, increases the loss factor. As amplitude increases further and the frictional effect is fully manifested, both stiffness and loss factor are further reduced, which are then maintained relatively constant by the remaining polymer chains. Motivated by the above physical mechanisms, a number of researchers have combined springs and frictional slides to represent the filler and rubber compound in filled rubbers or elastomers.

In the model developed by Tarzanin *et al.* [7], the elastomeric behavior was represented by a nonlinear spring and a nonlinear Coulomb friction damper. This model was based on single frequency elastomer data and matched the value of energy dissipation per cycle. Panda *et al.* [2] replaced the Coulomb friction damping element with a variable friction damping element whose force was calculated based on the peak displacement of excitation when the velocity was zero. This model correlated well with the experimental hysteresis data, but the effectiveness of this model over a range of amplitudes and frequencies has not been demonstrated in the literature. As early as 1930, Timoshenko [8] suggested that general hysteretic systems consist of a large number of ideal elasto-plastic elements with different yield levels. Iwan [9, 10] further developed a distributed-element model to study the steadystate dynamic response of a hysteretic system. Instead of specifying a distribution function numerically to agree with experimental data, a constant band-limited statistical function was used to define yield properties of the slide elements in this model. This model proved successful in predicting steady-state frequency response of a hysteretic system using a method of linearization. However, the excitation amplitude must be known as *a priori* while this model is applied in the analysis, which is only applicable in steady-state response prediction. The theory of triboelasticity [11] also stated that the behavior of a filled elastomer can be represented by a large or infinite number of alternate springs and frictional slides in series, and each slide has a constant yield force and each spring has a constant stiffness. Coveney *et al.* [1] developed a three-parameters standard triboelastic solid (STS) model based on the theory of triboelasticity and further developed a four-parameters ratedependent triboelastic (RT) model. These models gave a satisfactory representation of the material behavior. However, because the yield force was fixed along different slides, these models showed less flexibility in representing the amplitude dependent behavior for different filled-level materials.

308 Advanced Elastomers – Technology, Properties and Applications

these models.

such as creep and relaxation, a viscoelastic solid can be represented as a spring in series with Kelvin elements (if one Kelvin element is used, then a Zener model results [3]). Gandhi and Chopra [4] developed a nonlinear viscoelastic solid model in which a nonlinear lead spring was used in series with a single linear Kelvin chain. Using this model, the variation of complex moduli with oscillation amplitude closely matched experimental data. Felker *et al..* [5] further developed a nonlinear complex modulus model based on a single Kelvin chain, in which the spring force was a nonlinear function of the displacement, and the damping force was a nonlinear function of displacement and velocity. This model was used to describe the amplitude dependent moduli and to study dual frequency damper motions. The parameters in all of these models were identified using amplitude-dependent complex modulus data. As a result, the nonlinear hysteresis behavior of could not be captured using

Elastomeric materials typically demonstrate tribo-elastic behavior [1]. Physically, the amplitude dependent behavior exhibited by the elastomer is results from the interaction of the fillers with the rubber or elastomeric matrix materials [6]. Before large deformation of a filled elastomeric damper, an intact filler structure displays a large stiffness and low loss factor for small amplitudes. As the amplitude increases, the filler structure breaks resulting in a stiffness reduction. However, the breaking of filler structures, which is similar to friction, increases the loss factor. As amplitude increases further and the frictional effect is fully manifested, both stiffness and loss factor are further reduced, which are then maintained relatively constant by the remaining polymer chains. Motivated by the above physical mechanisms, a number of researchers have combined springs and frictional slides

In the model developed by Tarzanin *et al.* [7], the elastomeric behavior was represented by a nonlinear spring and a nonlinear Coulomb friction damper. This model was based on single frequency elastomer data and matched the value of energy dissipation per cycle. Panda *et al.* [2] replaced the Coulomb friction damping element with a variable friction damping element whose force was calculated based on the peak displacement of excitation when the velocity was zero. This model correlated well with the experimental hysteresis data, but the effectiveness of this model over a range of amplitudes and frequencies has not been demonstrated in the literature. As early as 1930, Timoshenko [8] suggested that general hysteretic systems consist of a large number of ideal elasto-plastic elements with different yield levels. Iwan [9, 10] further developed a distributed-element model to study the steadystate dynamic response of a hysteretic system. Instead of specifying a distribution function numerically to agree with experimental data, a constant band-limited statistical function was used to define yield properties of the slide elements in this model. This model proved successful in predicting steady-state frequency response of a hysteretic system using a method of linearization. However, the excitation amplitude must be known as *a priori* while this model is applied in the analysis, which is only applicable in steady-state response prediction. The theory of triboelasticity [11] also stated that the behavior of a filled elastomer can be represented by a large or infinite number of alternate springs and frictional slides in series, and each slide has a constant yield force and each spring has a constant stiffness.

to represent the filler and rubber compound in filled rubbers or elastomers.

Alternative elastomeric models were developed using internal variable or nonlinear integral equations. Strganac [12] used a stress shift function to formulate a nonlinear time domain model for elastomers, but the nonlinear integral formulation in the model was difficult to implement in traditional aeromechanical analysis. Lesieutre and Bianchini [13] developed the anelastic displacement field (ADF) method to describe the frequency-dependent behavior of viscoelastic materials. It was based on the notion of scalar internal variables or augmenting thermodynamic fields (ATF) [14] that described the interaction of the displacement field with irreversible processes occurring at the material level. In the ADF approach, the effects of the thermodynamic processes were focused on the displacement field, which consists of both elastic and anelastic fields. The anelastic part may be further subdivided to consider the effects of multiple relaxation processes. Although there is no explicit physical interpretation when multi-anelastic elements are involved, one single ADF model is mechanically analogous to the Zener model. In order to capture the characteristic nonlinear hysteretic behavior of elastomeric materials, Govindswamy *et al.* [15] developed a nonlinear ADF model, in which the linear ADF parameters were replaced with nonlinear terms. The model captured the variations of the complex modulus with amplitude, and performed as well in matching the strain vs. stress hysteresis cycle. Furthermore, other functional forms for the ADF parameters were introduced in order to improve hysteresis loop predictions. Brackbill *et al.* [16] improved the nonlinear ADF model by adding rate independent nonlinearity, in which friction-damping and linear-spring elements in parallel with the baseline nonlinear ADF model were used to provide additional amplitude dependent relaxation behavior. As many as sixteen parameters were used to construct the model. Although the complex moduli were fitted well in certain amplitude and frequency ranges, the performance of the model in predicting nonlinear hysteresis behavior could still be improved. Moreover, the process to determine model parameters was complicated by the fact that some parameters were chosen by empirical observation. In a recent study, Ramrakhyani *et al.* [17] developed an ADF based model containing nonlinear fractional derivatives and frictional elements. This model used eight parameters instead of sixteen parameters to capture the amplitude-dependent and mild frequency-dependent modulus. However, the prediction of hysteresis loops was not noticeably improved, and the determination of model parameters remains complicated.

In this chapter, a hysteresis model is introduced to study the anelastic behavior of elastomers under harmonic excitation. Using this model, the behavior of the elastomer is analogous to the behavior of a nonlinear Kelvin element, in which the stiffness is a nonlinear monotonic function of displacement and the damping is a monotonic rate dependent hyperbolic tangent function. Although, the hysteresis model can capture the hysteresis

behavior of the elastomer, the use of simple nonlinear spring and damping elements are not sufficient because model parameters, at the very least, are amplitude dependent. Thus, a computationally efficient and precise model for the elastomer is developed from a profound understanding of the damping mechanism within the damping material. As mentioned before, the anelastic behavior demonstrated by the elastomer is mostly based on the interaction between fillers and rubber compound inside the filled elastomeric materials. Based on this physical mechanism, a non-uniform distribution of rate-dependent elastoslide elements is used to emulate filler structure behavior and a parallel linear spring (and linear viscous damping) is used to represent the remaining polymer stiffness (and damping). Extensive testing including single frequency and dual frequency testing is conducted for material characterization, identification of model parameters and validation of the model. It is shown that this material model can be used for damping element design and can be integrated into numerical analysis for dynamic systems. It is also shown that the non-uniform distributed rate-dependent elastomer model is applicable in complex loading conditions without any *priori* information, and the flexibility in determining model parameters provides a potential to improve the model and to apply the model for elastomers with different filler structures.

### **2. Characterization of elastomers**

To characterize elastomeric materials under harmonic excitation, dynamic tests were conducted for two different elastomer configurations. The first elastomeric configuration was a double lap shear specimen as shown in Fig. 1a, or flat linear bearing, hereinafter called elastomeric specimen 1. The second elastomeric configuration was a linear concentric tubular bearing as shown in Fig. 1b, hereafter denoted the elastomer specimen 2. Both specimens were characterized under pure shear deformation. Testing was carried out with varying excitation amplitudes and frequencies, and all tests were conducted at room temperature: 25C.

As shown in Fig. 1a, elastomer specimen 1 was a double lap shear speciment that is comprised of three parallel brass plates between which the elastomeric material is sandwiched symmetrically. To study the effect of preload on the behavior of the elastomer, elastomeric specimen testing was conducted with and wihtout preload. Preload was applied to the elastomeric specimen by compressing the double lap shear specimen 10% of the width of the specimen using a simple vise. The testing setup for the elastomer specimen 1 is shown in Fig. 2a. A 24.466 kN servo-hydraulic MTS test machine was used to characterize the specimen. Fixtures and grips were designed and machined to hold the specimen in place. A hydraulic power supply (HPS) unit supplied the servo fluid to the testing machine for power, and the specimen was loaded and tested on the load frame. An actuator provided the sinusoidal loading, and an LVDT sensor measured displacement while a load cell measured force. Single and dual frequency tests were conducted using this load frame. In the dual frequency test, an HP 8904A multi-function synthesizer was used to generate and sum the sinusoidal signals for both frequencies. The single frequency test was conducted with displacement control for excitation amplitude ranging from 0.25 mm to 5 mm, i.e. 2.5% to 50% shear strain, in increments of 0.25 mm. The frequencies were chosen as 2.5 Hz, 5 Hz and 7.5 Hz. The dual frequency testing was carried out at a combination of 5 Hz and 7.5 Hz, and the amplitudes for 7.5 Hz were 0.5, 1.5, 2.5, 3.5, and 4.5 mm, respectively, while the amplitude for 5 Hz was maintained the same as for the single frequency tests.

(a) Specimen 1: Double Lap Shear

(b) Specimen 2: Cylindrical Tubular Shear

**Figure 1.** Elastomer Specimen (Unit: mm)

310 Advanced Elastomers – Technology, Properties and Applications

elastomers with different filler structures.

**2. Characterization of elastomers** 

temperature: 25C.

behavior of the elastomer, the use of simple nonlinear spring and damping elements are not sufficient because model parameters, at the very least, are amplitude dependent. Thus, a computationally efficient and precise model for the elastomer is developed from a profound understanding of the damping mechanism within the damping material. As mentioned before, the anelastic behavior demonstrated by the elastomer is mostly based on the interaction between fillers and rubber compound inside the filled elastomeric materials. Based on this physical mechanism, a non-uniform distribution of rate-dependent elastoslide elements is used to emulate filler structure behavior and a parallel linear spring (and linear viscous damping) is used to represent the remaining polymer stiffness (and damping). Extensive testing including single frequency and dual frequency testing is conducted for material characterization, identification of model parameters and validation of the model. It is shown that this material model can be used for damping element design and can be integrated into numerical analysis for dynamic systems. It is also shown that the non-uniform distributed rate-dependent elastomer model is applicable in complex loading conditions without any *priori* information, and the flexibility in determining model parameters provides a potential to improve the model and to apply the model for

To characterize elastomeric materials under harmonic excitation, dynamic tests were conducted for two different elastomer configurations. The first elastomeric configuration was a double lap shear specimen as shown in Fig. 1a, or flat linear bearing, hereinafter called elastomeric specimen 1. The second elastomeric configuration was a linear concentric tubular bearing as shown in Fig. 1b, hereafter denoted the elastomer specimen 2. Both specimens were characterized under pure shear deformation. Testing was carried out with varying excitation amplitudes and frequencies, and all tests were conducted at room

As shown in Fig. 1a, elastomer specimen 1 was a double lap shear speciment that is comprised of three parallel brass plates between which the elastomeric material is sandwiched symmetrically. To study the effect of preload on the behavior of the elastomer, elastomeric specimen testing was conducted with and wihtout preload. Preload was applied to the elastomeric specimen by compressing the double lap shear specimen 10% of the width of the specimen using a simple vise. The testing setup for the elastomer specimen 1 is shown in Fig. 2a. A 24.466 kN servo-hydraulic MTS test machine was used to characterize the specimen. Fixtures and grips were designed and machined to hold the specimen in place. A hydraulic power supply (HPS) unit supplied the servo fluid to the testing machine for power, and the specimen was loaded and tested on the load frame. An actuator provided the sinusoidal loading, and an LVDT sensor measured displacement while a load cell measured force. Single and dual frequency tests were conducted using this load frame. In the dual frequency test, an HP 8904A multi-function synthesizer was used to generate and sum the sinusoidal signals for both frequencies. The single frequency test was conducted with displacement control for excitation amplitude ranging from 0.25 mm to 5 mm, i.e. 2.5%

An important effect of filler materials in the filled elastomeric specimen is stress-softening. If an elastomeric sample is stretched for the first time to 100% followed by a release in the strain and then stretched again to 200%, there is a softening in the strain of up to 100% after which it continues in a manner of following the first cycle. This stress softening effect was first discovered by Mullins and is called the 'Mullins Effect' [18]. To account for this phenomenon, the test samples were first cycled and loosened before the actual tests by exciting them at 1 Hz frequency and 5 mm for 300 cycles since 5 mm is the maximum amplitude during tests. Stress relaxation is also shown in the case of dynamic loading. As the material is subjected to cycling loading, energy dissipation in the material heats up the material and results in elevated temperature softening. Usually, material self-heating and

other unsteady effects require about 250 seconds to stabilize and reach a steady state. Hence, in order to ensure temperature stabilization and consistency of data, during a normal test, the elastomeric sample was typically excited at the test frequency and amplitude for 300 seconds before collecting data. For simplification, the stress softening and relaxation effects were not considered in the modeling process such that the model parameters were independent of the loading level and temperature.

(a) Testing Setup for Specimen 1

(b) Testing Setup for Specimen 2

#### **Figure 2.** Testing Setup

As shown in Fig. 1b, elastomer specimen 2 was fabricated from two concentric cylindrical metal tubes, with an elastomeric layer sandwiched between the outer and inner tubes. The volume enclosed by the inner tube forms a cylindrical inner chamber, and a threaded trapezoidal column is attached to one end of the inner tube. As shown in Fig. 2b, the specimen was installed in the MTS testing machine, the outer tube was attached to the load cell on the fixed end of the MTS machine, and the inner tube was connected to the actuator through an adapter. Thus, the axial translation of the actuator induced a relative translation between the inner tube and the outer tube which in turn led to a shear deformation of the elastomer along the tube length. The specimen was excited in displacement control by a sinusoidal signal, and the displacement and force were measured by the LVDT sensor and load cell of the MTS machine. The excitation amplitude ranged from 0.25 mm to 1 mm in increments of 0.25 mm (approximately 5% to 20% shear) at three different frequencies of 2.5, 5.0, and 7.5 Hz, respectively.

All test data were collected using a high sampling frequency (2048 Hz) such that most higher harmonic components in the measured nonlinear force were included. To reduce the noise of the sinusoidal displacement signal, a Fourier series was used to reconstruct the input displacement. The reconstructed displacement signal was then differentiated to obtain the velocity signal. The Fourier series expansion of the input displacement, *x(t)*, is

$$\mathbf{x}\left(t\right) = \frac{\mathbf{x}\_0}{2} + \sum\_{k=1}^{\Leftrightarrow} \left[ X\_{c,k} \cos\left(k\alpha t\right) + X\_{s,k} \sin\left(k\alpha t\right) \right] \tag{1}$$

where,

312 Advanced Elastomers – Technology, Properties and Applications

independent of the loading level and temperature.

**Figure 2.** Testing Setup

other unsteady effects require about 250 seconds to stabilize and reach a steady state. Hence, in order to ensure temperature stabilization and consistency of data, during a normal test, the elastomeric sample was typically excited at the test frequency and amplitude for 300 seconds before collecting data. For simplification, the stress softening and relaxation effects were not considered in the modeling process such that the model parameters were

As shown in Fig. 1b, elastomer specimen 2 was fabricated from two concentric cylindrical metal tubes, with an elastomeric layer sandwiched between the outer and inner tubes. The

(b) Testing Setup for Specimen 2

(a) Testing Setup for Specimen 1

$$\begin{aligned} X\_{c,k} &= \frac{\alpha}{\pi} \int\_0^{2\pi} x\left(t\right) \cos\left(k\alpha t\right) dt\\ X\_{s,k} &= \frac{\alpha}{\pi} \int\_0^{2\pi} x\left(t\right) \sin\left(k\alpha t\right) dt \end{aligned} \tag{2}$$

In Eq. (1), *x0=Xc,0*. For single frequency data processing, any bias and higher harmonics were filtered, so that only the frequency of interest, , remained, that is 2.5, 5.0 and 7.5 Hz. For dual frequency testing, the general equation for the input dual displacement signal is written as:

$$\mathbf{x}\begin{pmatrix} t \\ \end{pmatrix} = X\_1 \sin\left(\Omega\_1 t\right) + X\_2 \cos\left(\Omega\_2 t\right) \tag{3}$$

where 1 and 2 are correspnding frequencies of 5 and 7.5 Hz, and *X1* and *X2* are the amplitudes at each frequency. The signal is periodic with a frequency corresponding to the highest common factor of both harmonics, i.e., 2.5 Hz. The displacement signal was filtered using 2.5 Hz as the base frequency. The first three harmonics were needed to reconstruct the dual frequency displacement signal in order to capture 1 and Ω2. Due to the nonlinearity of elastomers, the higher harmonics of the measured force were not filtered.

A typical approach used for characterizing elastomer behavior is the complex stiffness. The linearized complex stiffness, *K\** , is composed of an in-phase or storage stiffness, *K*, and a quadrature or loss stiffness, *K*, as follows:

$$K^\* = K' + jK'' \tag{4}$$

Therefore, the elastomer force can be written as the summation of an in-phase spring force and a quadrature damping force, and the elastomer force can be approximated by the first Fourier sine and cosine components at the charterized frequencies, i.e. 2.5, 5.0 and 7.5 Hz:

$$\begin{aligned} F\left(t\right) &= F\_c \cos\left(\alpha t\right) + F\_s \sin\left(\alpha t\right) \\ &= K' \mathbf{x}\left(t\right) + \frac{K''}{\alpha} \dot{\mathbf{x}}\left(t\right) \end{aligned} \tag{5}$$

where *Fc* and *Fs* are the first harmonic Fourier coefficients of the measured force. The storage stiffness, *K*, and the loss stiffness, *K*, are determined using the following equations:

$$\begin{aligned} K'(\alpha) &= \frac{F\_c X\_c + F\_s X\_s}{X\_c^2 + X\_s^2} \\ K''(\alpha) &= \frac{F\_c X\_s - F\_s X\_c}{X\_c^2 + X\_s^2} \end{aligned} \tag{6}$$

where, *Xc* and *Xs* are the first harmonic Fourier coefficients of *x(t)*. The loss factor, , is also used to measure the relative levels of the loss stiffness to the storage stiffness. The ratio is written as:

$$
\eta = \frac{K''}{K'} \tag{7}
$$

Typical linear characterization results for elastomeric specimen 1 in the absence of preload are shown in Fig. 3. These plots indicate that the linearized storage and loss stiffness of the specimen are highly amplitude dependent at low amplitudes. For smaller amplitudes, the rate of change of the storage stiffness and the loss stiffness is much greater than one for larger amplitudes. However, the complex stiffness does not change substantially over the narrow frequency range tested. The loss factor, , is also strongly dependent on amplitude. The maximum value of the loss factor is as high as 1.025 for this filled elastomer. This elastomeric material performs most effectively as a damping material within the amplitude range of 0.76 - 1.27 mm 7.5% to 12.5% shear strain), in which the loss factor, *>1*. Similar results were exhibited by the elastomer under preload. In general, the linearized behavior of the elastomeric specimen is highly amplitude dependent and weakly frequency dependent in this frequency range.

The measured complex modulus and loss factor of the elastomer specimen 2 are shown in Fig. 4. Both in-phase (storage) and quadrature (loss) stiffness demonstrate moderate amplitude dependence and weak frequency dependence. In contrast to the characteristics of the elastomeric specimen 1, the in-phase stiffness of specimen 2 is much greater than its quadrature stiffness, and both in-phase and quadrature stiffness vary with the displacement amplitude at a similar rate. Thus, the loss factor of the specimen is quite low (around 0.25 to 0.3) and almost constant over the range of amplitude and frequency tested.

**Figure 3.** Linear Characterization of Specimen 1

, and the loss stiffness, *K*

narrow frequency range tested. The loss factor,

in this frequency range.

stiffness, *K*

written as:

\* *K K jK* (4)

(5)

(6)

(7)

, is also strongly dependent on amplitude.

*>1*. Similar

, is also

, are determined using the following equations:

Therefore, the elastomer force can be written as the summation of an in-phase spring force and a quadrature damping force, and the elastomer force can be approximated by the first Fourier sine and cosine components at the charterized frequencies, i.e. 2.5, 5.0 and 7.5 Hz:

cos sin *c s Ft F t F t <sup>K</sup> Kx t x t*

2 2

*cc ss c s cs sc c s*

*FX FX*

*X X FX FX*

2 2

*X X*

where *Fc* and *Fs* are the first harmonic Fourier coefficients of the measured force. The storage

used to measure the relative levels of the loss stiffness to the storage stiffness. The ratio is

*K K* 

Typical linear characterization results for elastomeric specimen 1 in the absence of preload are shown in Fig. 3. These plots indicate that the linearized storage and loss stiffness of the specimen are highly amplitude dependent at low amplitudes. For smaller amplitudes, the rate of change of the storage stiffness and the loss stiffness is much greater than one for larger amplitudes. However, the complex stiffness does not change substantially over the

The maximum value of the loss factor is as high as 1.025 for this filled elastomer. This elastomeric material performs most effectively as a damping material within the amplitude

results were exhibited by the elastomer under preload. In general, the linearized behavior of the elastomeric specimen is highly amplitude dependent and weakly frequency dependent

The measured complex modulus and loss factor of the elastomer specimen 2 are shown in Fig. 4. Both in-phase (storage) and quadrature (loss) stiffness demonstrate moderate amplitude dependence and weak frequency dependence. In contrast to the characteristics of the elastomeric specimen 1, the in-phase stiffness of specimen 2 is much greater than its quadrature stiffness, and both in-phase and quadrature stiffness vary with the displacement amplitude at a similar rate. Thus, the loss factor of the specimen is quite low (around 0.25 to

range of 0.76 - 1.27 mm 7.5% to 12.5% shear strain), in which the loss factor,

0.3) and almost constant over the range of amplitude and frequency tested.

*K*

*K*

where, *Xc* and *Xs* are the first harmonic Fourier coefficients of *x(t)*. The loss factor,

**Figure 4.** Linear Characterization of Specimen 2

The single frequency linear characterization can capture the general trends of the in-phase and quadrature stiffness for a filled elastomer. However, this linear analysis cannot be used to accurately reconstruct the nonlinear hysteresis behavior exhibited by the elastomer [19]. Therefore, nonlinear modeling methods are requried to accurately describe the material behavior of elastomers, and two modeling methods are decribed in the following sections

#### **3. Nonlinear hysteresis model**

#### **3.1. Modeling approach**

The mechanical properties of a linear viscoelastic material are represented by a Kelvin model, which consists of a spring and dashpot in parallel. The spring and damping coefficients are constants, and the Kelvin model is equivalent to a complex modulus approach. Although the Kelvin model cannot describe the relaxation process after a constant strain is applied to a material specimen, it can successfully characterize stiffness and damping during steady-state harmonic excitation. Therefore, it is the simplest approach to describe the nonlinear behavior of elastomeric materials based on a Kelvin model.

As mentioned before, when a linear viscous material is subjected to sinusoidal loading, the stiffness force is an in-phase force with a constant stiffness, and the damping force with a constant damping is a quadrature force. For an elastomeric specimen, both stiffness and damping are not constant. To extract the stiffness force, *Fs*, and damping force, *Fd*, from the experimental force data, *F*, the displacement phase angle, *(t)*, at arbitrary time *t* should be known firstly from the reconstructed displacement signal. A single frequency sinusoidal displacement can be written as:

$$\mathbf{x}\left(t\right) = \sqrt{\mathbf{X}\_c^2 + \mathbf{X}\_s^2} \sin\left(\phi\left(t\right)\right) \tag{8}$$

where

$$\phi\left(t\right) = \alpha t + \arctan\left(\frac{X\_c}{X\_s}\right) \tag{9}$$

From a given start time, *t0*, when *(t0)*=2*k*+/2, *k*=0,1,.., one cycle of the force data is used to determine the stiffness force. Since the stiffness force has the same phase angle as *(t)* and the damping force is /2 ahead of *(t)*, the stiffness force, *Fs*, in one cycle can be obtained as:

$$\begin{aligned} F\_s \left[ \phi \left( t \right) \right] &= \frac{F \left[ \phi \left( t \right) \right] + F \left[ 2\pi - \phi \left( t \right) \right]}{2} \\ F\_s \left[ 2\pi - \phi \left( t \right) \right] &= F\_s \left[ \phi \left( t \right) \right] \\ \phi \left( t \right) &= \phi \left( t\_0 \right), ..., \phi \left( t\_0 \right) + \pi \end{aligned} \tag{10}$$

Similarly, from a different start point *(t0')*=2*k*, *k*=0,1,…, the damping force, *Fd*, is obtained as follows:

$$\begin{aligned} F\_d\left[\phi\left(t\right)\right] &= \frac{F\left[\phi\left(t\right)\right] + F\left[2\pi - \phi\left(t\right)\right]}{2} \\ F\_d\left[2\pi - \phi\left(t\right)\right] &= F\_d\left[\phi\left(t\right)\right] \\ \phi\left(t\right) &= \phi\left(t\_0\right), \dots, \phi\left(t\_0\right) + \pi \end{aligned} \tag{11}$$

The identified nonlinear stiffness and damping forces are shown in Fig. 5. In Figure 5a and b, the experimental force-displacement and force-velocity hysteresis cycles are plotted using dotted lines for a sinusoidal displacement excitation at 5Hz and 1.5mm amplitude. The stiffness force, *Fs*, and the damping force, *Fd*, are denoted as solid lines in the forcedisplacement plane and the force-velocity plane, respectively. Clearly, the stiffness force is a nonlinear monotonic function of displacement, and the damping force is a nonlinear monotonic function of velocity. Then, in order to establish a nonlinear Kelvin model, effective functional forms must be found to represent both the nonlinear stiffness and damping forces, respectively.

In Fig. 5a, it is shown that the stiffness force, *Fs*, is almost linear but increases with increasing displacement. By observing different stiffness lines at different displacement amplitudes, it is noted that the nonlinear behavior of stiffness always occurs while the velocity is close to zero. It implies that the nonlinearity of the stiffness is due to velocity but not displacement. Therefore, the nonlinear stiffness model is given by:

$$\hat{F}\_s = \left[K + \delta K \exp\left(\lambda\_k \left|\dot{\mathbf{x}}\right|^{1.5}\right)\right] \mathbf{x} \tag{12}$$

where, *K* is the linear stiffness, *K* is the nonlinear stiffness increment, and *<sup>k</sup>* is used to describe the nonlinear stiffness slope. These parameters can be determined by a nonlinear curve-fitting algorithm, in which an objective function, *Js*, is minimized as follows:

$$\min\_{\mathcal{K}, \delta \mathcal{K}, \lambda\_K} J\_s \left( K, \delta \mathcal{K}, \lambda\_K \right) = \sum\_{i=1}^n \left[ \hat{F}\_s \left( t\_i \right) - F\_s \left( t\_i \right) \right]^2 \tag{13}$$

where, *n* is the number of data points.

316 Advanced Elastomers – Technology, Properties and Applications

displacement can be written as:

From a given start time, *t0*, when

Similarly, from a different start point

damping forces, respectively.

/2 ahead of

the damping force is

as follows:

where

experimental force data, *F*, the displacement phase angle,

As mentioned before, when a linear viscous material is subjected to sinusoidal loading, the stiffness force is an in-phase force with a constant stiffness, and the damping force with a constant damping is a quadrature force. For an elastomeric specimen, both stiffness and damping are not constant. To extract the stiffness force, *Fs*, and damping force, *Fd*, from the

known firstly from the reconstructed displacement signal. A single frequency sinusoidal

arctan *<sup>c</sup>*

*Ft F t*

The identified nonlinear stiffness and damping forces are shown in Fig. 5. In Figure 5a and b, the experimental force-displacement and force-velocity hysteresis cycles are plotted using dotted lines for a sinusoidal displacement excitation at 5Hz and 1.5mm amplitude. The stiffness force, *Fs*, and the damping force, *Fd*, are denoted as solid lines in the forcedisplacement plane and the force-velocity plane, respectively. Clearly, the stiffness force is a nonlinear monotonic function of displacement, and the damping force is a nonlinear monotonic function of velocity. Then, in order to establish a nonlinear Kelvin model, effective functional forms must be found to represent both the nonlinear stiffness and

*Ft F t*

0 0

*(t0')*=2*k*

 ' ' 0 0

*tt t*

*F tFt*

*d d*

,...,

 

*F t Ft tt t*

*s s*

,...,

 

*t t*

determine the stiffness force. Since the stiffness force has the same phase angle as

*(t0)*=2*k*+

2

2

*F t*

*d*

*F t*

*s*

*s X*

*X* 

2 2

2 2

  2 2 sin *c s xt X X t* (8)

*(t)*, the stiffness force, *Fs*, in one cycle can be obtained as:

*(t)*, at arbitrary time *t* should be

(9)

(10)

(11)

*(t)* and

/2, *k*=0,1,.., one cycle of the force data is used to

, *k*=0,1,…, the damping force, *Fd*, is obtained

As shown in Fig. 5b, the shape of the nonlinear damping force, *Fd*, is reminiscent of a friction damping. This implies that filled elastomeric materials exhibit anelasticity instead of viscoelasticity. In some references, an inverse hyperbolic tangent [20] or exponential [21] functional approximation has been suggested to describe the tribo-elastic behavior in filled elastomers. In our model, a hyperbolic tangent function is used due to its simplicity and computational efficiency, so that the analytical friction damping force is assumed to behave as:

$$
\hat{F}\_d = F\_y \tanh\left(\lambda \dot{x}\right) \tag{14}
$$

where, *Fy* and represent the yield force and yield parameter respectively. Similarly, the nonlinear damping parameters also can be determined by minimizing the objective function, *Jd*, as follows:

$$\min\_{F\_{y'},\lambda} J\_d\left(F\_{y'}\lambda\right) = \sum\_{i=1}^n \left[\hat{F}\_d\left(t\_i\right) - F\_d\left(t\_i\right)\right]^2\tag{15}$$

The total predicted force due to the displacement excitation is the summation of the stiffness model and damping model as follows:

$$
\hat{F} = \hat{F}\_d + \hat{F}\_s \tag{16}
$$

The reconstructed force is shown as the solid line in the Fig. 5c, and the analytical forcedisplacement hysteresis matches the experimental data very well. Moreover, using the identified friction function, a force-displacement hysteresis due to friction damping is reconstructed. It is known that energy dissipation due to the damping is proportional to the enclosed area inside the damping force-displacement hysteresis loop. The area enclosed by

the total force-displacement hysteresis loop is equal to the area enclosed by the damping hysteresis. This proves that the nonlinear damping function precisely describes the energy dissipation while the elastomer specimen is under sinusoidal loading.

**Figure 5.** Hysteresis Used for Decoupling of Stiffness and Damping

Briefly, this nonlinear Kelvin model is based on a hysteresis modeling approach developed from damper modeling efforts. Since the nonlinear hysteresis loops of an elastomeric damper are described by a nonlinear monotonic stiffness and a nonlinear monotonic damping function respectively, the model parameters can be determined separately and efficiently. The modeling results can precisely capture the force-displacement time history data of the elastomer. Meanwhile, the introduction of the friction damping physically emphasizes the anelasticity of elastomeric materials. Since the hysteresis cycles for different amplitudes and frequencies are different, it is helpful to study the model parameter variations and get insight into the nonlinear behavior of elastomers by characterizing the different hysteresis cycles.

#### **3.2. Modeling results**

318 Advanced Elastomers – Technology, Properties and Applications

the total force-displacement hysteresis loop is equal to the area enclosed by the damping hysteresis. This proves that the nonlinear damping function precisely describes the energy

(a) (b)

dissipation while the elastomer specimen is under sinusoidal loading.

**Figure 5.** Hysteresis Used for Decoupling of Stiffness and Damping

different hysteresis cycles.

Briefly, this nonlinear Kelvin model is based on a hysteresis modeling approach developed from damper modeling efforts. Since the nonlinear hysteresis loops of an elastomeric damper are described by a nonlinear monotonic stiffness and a nonlinear monotonic damping function respectively, the model parameters can be determined separately and efficiently. The modeling results can precisely capture the force-displacement time history data of the elastomer. Meanwhile, the introduction of the friction damping physically emphasizes the anelasticity of elastomeric materials. Since the hysteresis cycles for different amplitudes and frequencies are different, it is helpful to study the model parameter variations and get insight into the nonlinear behavior of elastomers by characterizing the

(c)

The model characterization was obtained based on three sets of force-displacement time history data of the elastomer specimen 1 at three different frequencies, i.e. 2.5Hz, 5Hz and 7.5Hz. For every frequency, there were twenty displacement amplitudes ranging from 0.25mm to 5mm. After the model parameters were determined, the output force was predicted using known displacement input data, and then the force-displacement hysteresis was reconstructed. In Fig. 6, the reconstructed hysteresis curves and test data are shown for different amplitudes (0.5mm, 1.5mm, 2.5mm, and 3.5mm) at 5Hz frequency with zero preloading. It can be seen that at low amplitudes, the experimental hysteresis plots are nearly elliptical in shape, while for higher amplitudes there is a deviation from this elliptical behavior. However, compared to the elastomer hysteresis models developed by Krishnan [22] and Snyder [23], the current nonlinear model more accurately captures nonlinear forcedisplacement time histories under sinusoidal loading.

**Figure 6.** Modeling Results of Hysteresis Model

Furthermore, model parameter variation as a function of the amplitude at different frequencies was studied. As shown in in Fig. 7a, all three nonlinear stiffness parameters are amplitude dependent. However, the linear stiffness, *K*, and the nonlinear stiffness increment, *K*, maintain nearly the same value at different frequencies, but the nonlinear stiffness slope, *<sup>k</sup>*, decreases as the frequency increases. Similarly, the identified damping parameters are shown in Fig. 7b. The yield force, *Fy*, for amplitudes above 1.5mm is nearly constant with both amplitude and frequency. At low amplitudes (<1.5mm), the optimized yield force shows large deviations since the damping is almost linear at small amplitude excitations and the model is insensitive to *Fy*. The yield parameter, , shows strong amplitude and frequency dependent characteristics.

Some interesting characteristics are noted in Fig. 8a, in which the stiffness is shown as a function of velocity amplitude. The nonlinear stiffness slope parameters of all three frequencies follow exactly the same curve, which is inversely proportional to the velocity amplitude. This implies that the nonlinear stiffness slope, *<sup>k</sup>*, only depends on the velocity amplitude. Similarly, as shown in Fig. 8b, the nonlinear damping yield parameter, , is also dependent on the acceleration amplitude. Both characteristics make this nonlinear model nearly independent of frequency and thus useful to predict dual frequency response.

**Figure 7.** Hysteresis Model Parameters as a Function of Amplitude and Frequency

(a) Stiffness Slope VS Velocity Amplitude (b) Yield Parameter VS Acceleration Amplitude

**Figure 8.** Hysteresis Model Parameters as a Function of Velocity and Acceleration Amplitude

To capture the behavior of the elastomer under dual frequency excitation, the elastomer model parameters were determined without frequency information because the excitation amplitude was known a *priori*. Specifically, both the linear stiffness and stiffness increment are inversely proportional to the displacement amplitude, and the nonlinear slope can be interpolated using the known velocity amplitude. Meanwhile, for nonlinear damping, the yield force is almost a constant and the yield parameter can be determined using the acceleration amplitude. As a result, the nonlinear dual frequency behavior of the elastomer could be predicted. In Fig. 9, the dual frequency displacement excitation was a summation of two single frequency inputs, i.e., 5 Hz and 7.5 Hz, and the reconstructed hysteresis cycles using the hysteresis model are compared to experimental data. It is shown that the modeling results can accurately capture the hysteresis behavior exhibited by the elastomer.

**Figure 9.** Dual Frequency Modeling Results Using Hysteresis Model

As shown in the nonlinear hysteresis modeling effort, the nonlinear forced response of elastomers under harmonic excitation consists of uncoupled nonlinear stiffness force and nonlinear damping force. Thus, this model is mechanically analogous to a nonlinear Kelvin model where the stiffness is a nonlinear monotonic function of displacement and the damping is a monotonic rate dependent friction function. Since it accurately describes the characteristics of the hysteresis loops, this model can predict steady state or harmonic forced response very well. However, since the model parameters are still amplitude dependent, this model cannot be easily used to describe the transient or stress relaxation behavior of the elastomer.

### **4. Distributed rate-dependent elasto-slide model**

### **4.1. Model development**

320 Advanced Elastomers – Technology, Properties and Applications

dependent on the acceleration amplitude. Both characteristics make this nonlinear model

(a) Stiffness Parameters (b) Damping Parameters

nearly independent of frequency and thus useful to predict dual frequency response.

**Figure 7.** Hysteresis Model Parameters as a Function of Amplitude and Frequency

**Figure 8.** Hysteresis Model Parameters as a Function of Velocity and Acceleration Amplitude

To capture the behavior of the elastomer under dual frequency excitation, the elastomer model parameters were determined without frequency information because the excitation amplitude was known a *priori*. Specifically, both the linear stiffness and stiffness increment are inversely proportional to the displacement amplitude, and the nonlinear slope can be interpolated using the known velocity amplitude. Meanwhile, for nonlinear damping, the yield force is almost a constant and the yield parameter can be determined using the acceleration amplitude. As a result, the nonlinear dual frequency behavior of the elastomer could be predicted. In Fig. 9, the dual frequency displacement excitation was a summation of two single frequency inputs, i.e., 5 Hz and 7.5 Hz, and the reconstructed hysteresis cycles

(a) Stiffness Slope VS Velocity Amplitude (b) Yield Parameter VS Acceleration Amplitude

The distributed rate-dependent elasto-slide model is shown in Fig. 10, in which a series of elasto-slide elements is combined in parallel with a constant linear spring. The model can be applied either in force-displacement relations or in stress-strain relations, but only the forcedisplacement formulation will be used in this study. Each elasto-slide element consists of a leading spring with stiffness *k/N* in series with a slide which has a yield force *fi \** , where *N* is the total number of elements. The yield force for each element is different, and the stiffness for each leading spring is assumed to be a constant. The ideal model assumes that the yield force is a Coulomb force, which is a constant as the slide moves at any speed. However, the rate-dependent elasto-slide model represents real friction behavior by representing the yield force of the slide as a function of the slide velocity. This function will be discussed later. For simplicity, the description of the model starts from using an ideal slide (which has an ideal

Coulomb force), and then the behavior of the model will be studied when the ideal slide is replaced by a rate-dependent slide.

**Figure 10.** Distributed Elasto-Slide Model

First, we will apply a displacement, *x*, to the ideal elasto-slide element, which assumes a constant Coulomb force. At the beginning, the displacement is small such that the consequent leading spring force is smaller than the yield force, so that only the spring is deformed. After the spring force reaches the yield force, the slide yields and a motion is induced, and the resisting force of the element remains constant with the same value as the yield force. Since each elasto-slide element in the model is assumed to have a different yield force level, this model presents gradual stiffness reduction as amplitude increases until such a condition as all elements have yielded. At that time, only the parallel spring, *k0*, remains to represent the polymer stiffness of the elastomer. Since the elastomer is a continuum, the total number of elements, *N*, approaches infinity, and, in the limit, the discrete yield force for a single slide is replaced with a distributed density within a certain yield force range. Alternatively, a mathematically equivalent model is shown in Fig. 11, in which *fe* represents the total resisting force of the elasto-slide elements.

**Figure 11.** Mathematically Equivalent Elasto-Slide Model

**Figure 12.** Effectiveness of Elasto-Slide Model

replaced by a rate-dependent slide.

**Figure 10.** Distributed Elasto-Slide Model

the total resisting force of the elasto-slide elements.

**Figure 11.** Mathematically Equivalent Elasto-Slide Model

Coulomb force), and then the behavior of the model will be studied when the ideal slide is

First, we will apply a displacement, *x*, to the ideal elasto-slide element, which assumes a constant Coulomb force. At the beginning, the displacement is small such that the consequent leading spring force is smaller than the yield force, so that only the spring is deformed. After the spring force reaches the yield force, the slide yields and a motion is induced, and the resisting force of the element remains constant with the same value as the yield force. Since each elasto-slide element in the model is assumed to have a different yield force level, this model presents gradual stiffness reduction as amplitude increases until such a condition as all elements have yielded. At that time, only the parallel spring, *k0*, remains to represent the polymer stiffness of the elastomer. Since the elastomer is a continuum, the total number of elements, *N*, approaches infinity, and, in the limit, the discrete yield force for a single slide is replaced with a distributed density within a certain yield force range. Alternatively, a mathematically equivalent model is shown in Fig. 11, in which *fe* represents

The following simulation using the model will show that the existence of the filler structures inside the elastomer can lead to hysteretic behaviors when the damper is cycled between fixed deflection limits. Analytically, a distribution function of the yield force is denoted as (*f \** ) such that the density of the slides with yield force *f \** is expressed as (*f \** )d*f \** . Using the ideal slide assumption, the slide behaves as a Coulomb friction element. Thus, upon initial loading, the leading spring at a certain yield element stretches with the displacement *x* until the spring force reaches the maximum slide yield force, that is:

$$\begin{aligned} \text{d}f &= k \text{x} \text{q} \left( f^\* \right) \text{d}f^\*, \dot{\text{x}} > 0, \text{} 0 \le \text{x} \le \frac{f^\*}{k} \\ \text{d}f &= f^\* \text{q} \left( f^\* \right) \text{d}f^\*, \dot{\text{x}} > 0, \text{} \ge \frac{f^\*}{k} \end{aligned} \tag{17}$$

where, *k* is the stiffness of the leading spring. If the direction of loading is reversed, the force-deflection relation is more complicated. Including the yielded and non-yielded elements, the force-deflection relation becomes:

$$\begin{aligned} \text{d}f &= \left[ \text{f}^\* - k \left( A - \text{x} \right) \right] \text{\text{\textquotedblleft}} \left( \text{f}^\* \right) \text{df}^\*, \text{\textquotedblright} < \text{0}, A - 2 \frac{\text{f}^\*}{k} \le \text{x} \le A, \text{\textquotedblleft} f^\* \le kA \\ \text{d}f &= -\text{f}^\* \text{\textquotedblleft} \left( \text{f}^\* \right) \text{df}^\*, \text{\textquotedblright} < \text{0}, \text{x} \le A - 2 \frac{\text{f}^\*}{k}, \text{f}^\* \le kA \\ \text{d}f &= k \text{x} \text{q} \left( \text{f}^\* \right) \text{df}^\*, \text{\textquotedbl} < \text{0}, \text{f}^\* > kA \end{aligned} \tag{18}$$

where *A* is the maximum deflection of the elastomer. An expression similar to Eq. 18 is obtained when the loading reaches the minimum deflection and is reversed again. This process continues until the loading is terminated. Clearly, at each time, only some of the elasto-slide elements have yielded. The total resisting force can be obtained by integrating all of the elasto-slide forces along the yield distribution region and by adding the spring force due to the residual polymer stiffness. Using Eq. 17 and adding the effect of the polymer stiffness, *k0*, the initial loading force is given as

$$f = \int\_0^{kx} f^\* \, \mathfrak{g}\left(f^\*\right) \mathrm{d}f^\* + k \mathrm{x} \int\_{kx}^{x} \mathfrak{g}\left(f^\*\right) \mathrm{d}f^\* + k\_0 \mathrm{x}, \dot{\mathrm{x}} > 0 \tag{19}$$

Similarly, the resisting force due to the reversed loading is obtained by integrating Eq. 18:

$$\begin{split} f &= \int\_{0}^{\underline{k(A-x)}} -f^\* \bullet \Big(f^\*\Big) \mathrm{d}f^\* + \int\_{\underline{k(A-x)}}^{\underline{k}A} \Big[ kx - \left(kA - f^\*\right) \Big] \mathrm{q}\Big(f^\*\Big) \mathrm{d}f^\* + kx \Big\|\_{\underline{k}A}^{\circ\prime} \mathrm{q}\Big(f^\*\Big) \mathrm{d}f^\* + k\_0 x\_{\prime} \\ \dot{x} &< 0, x \le A \end{split} \tag{20}$$

Thus, while an elastomer is under a sinusoidal displacement loading, a theoretical forcedeflection hysteresis cycle is shown as dashed line in Fig. 12, where *A* is 2.5 mm and frequency is 2.5 Hz. Compared to the experimental hysteresis shown as dotted line, the model prediction gives a good match to the experimental result.

The distributed elasto-slide model resembles the physical mechanism of an elastomer, so it can account for the nonlinear characteristics of the behavior demonstrated by the elastomer under a cyclic loading either in single frequency or multi frequencies. However, using the ideal slide with a Coulomb force, the maximum and minimum displacement of the excitation must be known for response calculation, which makes it impossible for the model to describe elastomer behavior under complex loading conditions. The ideal elasto-slide is also incapable of modeling frequency dependent properties and non-hysteretic behavior in the time domain such as stress relaxation or creep.

Actually, the Coulomb slide is only an idealized friction model. The practical friction behavior includes a preyield slip and a postyield steady resistance leading to a rate dependent damping effect [21, 24]. Thus, a rate-dependent elasto-slide model is introduced to improve the modeling performance. In the rate-dependent elasto-slide model, the Coulomb slide is replaced with a non-Coulombic friction function and the coupling between a slide and a leading spring is described by an internal displacement denoted as *x0* such that the slide force at a certain yield region is written as

$$\mathbf{d}f = f^\* \mathbf{q} \left( f^\* \right) \mathbf{d}f^\* \left( \frac{\dot{\mathbf{x}}\_0}{\upsilon\_r} \right)^{\frac{1}{p}} \tag{21}$$

where *p* is a positive odd integer and *vr* is a characteristic reference velocity. Coupled with the lead spring *k*(*f\** )d*f \** , the internal displacement *x0* at certain yield region is obtained using the function:

$$\frac{\dot{\mathbf{x}}\_0}{\upsilon\_r} = \left[ \frac{k}{f^\*} (\mathbf{x} - \mathbf{x}\_0) \right]^p \tag{22}$$

By integration over the whole yield force region, the total force due to the elasto-slide element is obtained as:

Anelastic Behavior in Filled Elastomers Under Harmonic Loading Using a Distributed Rate-Dependent Elasto-Slide Model 325

$$f\_e = \int\_0^\alpha k\left(x - x\_0\right) \wp\left(f^\*\right) \mathrm{d}f^\* \tag{23}$$

Adding the spring force due to the polymer stiffness, the damper force due to any deflection loading, *x*, is determined as:

324 Advanced Elastomers – Technology, Properties and Applications

0,

the lead spring *k*

element is obtained as:

the function:

(*f\** )d*f \**

*x xA*

*kA x*

 \* \*\* \*\* <sup>0</sup> <sup>0</sup> d d ,0 *kx kx f f f f kx f f k x x*

Similarly, the resisting force due to the reversed loading is obtained by integrating Eq. 18:

 \* \*\* \* \*\* \*\* <sup>2</sup> 0 0

Thus, while an elastomer is under a sinusoidal displacement loading, a theoretical forcedeflection hysteresis cycle is shown as dashed line in Fig. 12, where *A* is 2.5 mm and frequency is 2.5 Hz. Compared to the experimental hysteresis shown as dotted line, the

The distributed elasto-slide model resembles the physical mechanism of an elastomer, so it can account for the nonlinear characteristics of the behavior demonstrated by the elastomer under a cyclic loading either in single frequency or multi frequencies. However, using the ideal slide with a Coulomb force, the maximum and minimum displacement of the excitation must be known for response calculation, which makes it impossible for the model to describe elastomer behavior under complex loading conditions. The ideal elasto-slide is also incapable of modeling frequency dependent properties and non-hysteretic behavior in

Actually, the Coulomb slide is only an idealized friction model. The practical friction behavior includes a preyield slip and a postyield steady resistance leading to a rate dependent damping effect [21, 24]. Thus, a rate-dependent elasto-slide model is introduced to improve the modeling performance. In the rate-dependent elasto-slide model, the Coulomb slide is replaced with a non-Coulombic friction function and the coupling between a slide and a leading spring is described by an internal displacement denoted as *x0* such that

*<sup>x</sup> ff f f*

\* \*\* <sup>0</sup> d d *<sup>p</sup>*

where *p* is a positive odd integer and *vr* is a characteristic reference velocity. Coupled with

 <sup>0</sup> \* 0

By integration over the whole yield force region, the total force due to the elasto-slide

*x x*

*r*

*x k*

*v f*

1

, the internal displacement *x0* at certain yield region is obtained using

(22)

(21)

*r*

*p*

*v* 

 

*kA x kA <sup>f</sup> f f f kx kA f f f kx f f k x*

2

*kA*

model prediction gives a good match to the experimental result.

the time domain such as stress relaxation or creep.

the slide force at a certain yield region is written as

(19)

(20)

d d d,

$$f = f\_e + k\_0 \mathbf{x} = \int\_0^\infty k \left(\mathbf{x} - \mathbf{x}\_0\right) \mathbf{q} \left(f^\*\right) \mathbf{d}f^\* + k\_0 \mathbf{x} \tag{24}$$

This relation can be shown in Fig. 11. Eq. 22 is a typical well-posed initial-value problem, and numerical solution for this differential equation can be obtained given an initial condition. The simplest way to guarantee a stable solution for such a stiff initial-value problems is to adopt a predictor-corrector approach with the corrector iterated to convergence (PECE) [25]. In this approach, the numerical algorithm is based on the Adams-Bashforth four-step method as the predictor step and one iteration of the Adams-Moulton three-step method as the corrector step, with the starting values obtained from a fourthorder Runge-Kutta method. In accordance with the ratio of the yield force and the stiffness, *k/f\** , and an appropriate choice of *p*, the Adams-Bashforth-Moulton method gives relatively stable and fast convergence of a solution within a limited number time steps.

For an elastomer under a sinusoidal displacement excitation, the steady-state response predicted by the rate-dependent elasto-slide model is shown as the solid line in Fig. 12. The predicted hysteresis cycle correlates much better with the experimental data, especially at turning points of the loading deflection. It should be noted that there is no requirement for excitation amplitude information in this modeling process. Thus, the distributed ratedependent elasto-slide model can predict time domain forced response of an elastomer under a sinusoidal displacement excitation.

In order to apply the elastomer model to a dynamic system, a numerical method using MATLAB ODE algorithm was also evaluated. For a dynamic system with a governing equation:

$$\mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} + F\_d\left(\mathbf{x}, \dot{\mathbf{x}}\right) = F\left(t\right) \tag{25}$$

where, M, C, and K are mass matrix, equivalent damping matrix, and stiffness matrix, respectively, *Fd* is used to describe the force due to an elastomer, and *F* is a loading vector applied to the system. The size of the matrix depends on the degree of freedom of the system. For simplicity, only a 1-DOF system is considered now. In the distributed ratedependent elasto-slide damper model, there are theoretically an infinite number of internal variables, *x0*. However, the distributed yield stress usually falls within a limited range. Therefore, according to the form of a distribution function, the continuous yield force distribution area can be uniformly decomposed into *n* discrete elements from minimum to maximum yield force and each element has a yield force range, *f \** . At each yield stress range, *fi \** , the corresponding distribution is equal to the area of that element, ( *fi \** )*f\** . Each element has an internal variable denoted as *x0i*, *i=1,…,n*, and each *x0i* satisfies the Eq. 22. Thus, the elastomer force *Fd* can be described as:

$$F\_d = k \sum\_{i=1}^n \Phi\left(f\_i^\*\right) \Delta f^\*\left(\mathbf{x} - \mathbf{x}\_{0i}\right) + k\_0 \mathbf{x} \tag{26}$$

Rewriting Eq. 25 using a first order form and combining Eq. 22 and 26, the state equation of the system is expressed as

$$\begin{aligned} \dot{\bar{x}}\_{1} &= \dot{\bar{x}}\_{1} \\ \dot{\bar{x}}\_{2} &= \frac{F\left(t\right)}{M} - \frac{1}{M} \left[ \left(K + k\_{0}\right)\mathbf{x} - k \sum\_{i=1}^{n} \mathbf{q}\left(f\_{i}^{\*}\right) \Delta f^{\*}\left(\mathbf{x} - \mathbf{x}\_{0i}\right) \right] - \frac{\mathbf{C}}{M} \dot{\bar{x}}\_{1} \\ \dot{\bar{x}}\_{01} &= \left[\frac{k}{f\_{1}^{\*}} \left(\mathbf{x} - \mathbf{x}\_{01}\right)\right]^{p} v\_{r} \\ \dot{\bar{x}}\_{1} &\vdots \vdots \\ \dot{\bar{x}}\_{0n} &= \left[\frac{k}{f\_{n}^{\*}} \left(\mathbf{x} - \mathbf{x}\_{0n}\right)\right]^{p} v\_{r} \end{aligned} \tag{27}$$

This is a *(n+2)*th order state function. Using the ODE23 algorithm in MATLAB, the forced response or the transient response due to initial displacement and velocity can be solved numerically. For the system with more degrees of freedom, the state function can easily accommodate these additions by using additional number of states.

#### **4.2. Model parameters determination**

As seen in the construction of the model, the major parameters to be determined are the leading spring, *k*, and yield force distribution function, (*f\** ), for the distributed elasto-slide element, and the parallel spring, *k0*, for the remaining polymer stiffness. In the absence of the knowledge of the elastomer structure, the selection of these parameters is only based on experimental data in this stage. One possible selection of methods would be to make use of an experimentally determined initial loading curve.

The definition of the distribution function implies that (*f\** ) has to obey the following three constraints:

$$\begin{aligned} \int\_0^\infty \mathfrak{g}\left(f^\*\right) \mathrm{d}f^\* &= 1\\ \mathfrak{g}\left(f^\*\right) &\ge 0\\ 0 \le f^\* &\le \infty \end{aligned} \tag{28}$$

From the initial loading curve Eq. 19, yields

$$\begin{cases} \frac{\text{d}f}{\text{d}x} = k \int\_{kx}^{\infty} \text{\c}\big(f^{\*}\big) \text{d}f^{\*}\\ f^{\*} = kx \end{cases} \tag{29}$$

Then

326 Advanced Elastomers – Technology, Properties and Applications

*x x*

**4.2. Model parameters determination** 

constraints:

the system is expressed as

\* \*

*F k f f x x kx*

Rewriting Eq. 25 using a first order form and combining Eq. 22 and 26, the state equation of

*i*

This is a *(n+2)*th order state function. Using the ODE23 algorithm in MATLAB, the forced response or the transient response due to initial displacement and velocity can be solved numerically. For the system with more degrees of freedom, the state function can easily

As seen in the construction of the model, the major parameters to be determined are the

element, and the parallel spring, *k0*, for the remaining polymer stiffness. In the absence of the knowledge of the elastomer structure, the selection of these parameters is only based on experimental data in this stage. One possible selection of methods would be to make use of

\* \*

*f f*

0

d 1

\* \*

*f*

\*

<sup>d</sup> <sup>d</sup> d *kx*

 

*<sup>f</sup> k ff*

*f*

0 \*

0

\*

*x f kx* (*f\**

(*f\**

*F t C <sup>x</sup> Kkxk f f xx x <sup>M</sup> M M*

*di i*

1

1 *<sup>n</sup>*

*r*

*p*

*p*

*i*

 

01 \* 01 1

*k x xx v f*

0 0 \*

leading spring, *k*, and yield force distribution function,

an experimentally determined initial loading curve.

From the initial loading curve Eq. 19, yields

The definition of the distribution function implies that

*k x xx v f*

accommodate these additions by using additional number of states.

*n n r n*

 

*n*

0 0

\* \* 0 0 1

*i i*

(26)

(27)

), for the distributed elasto-slide

) has to obey the following three

(28)

(29)

$$\frac{\mathbf{d}^2 f}{\mathbf{d}x^2} = -k^2 \boldsymbol{\varrho} \left( f^\* \right) \tag{30}$$

Thus, the distribution function would be related to the curvature of the initial loading curve by the following formula

$$\exp\left(f^\*\right) = -\frac{1}{k^2} \frac{\mathrm{d}^2 f}{\mathrm{d}x^2} \tag{31}$$

Determination of the distribution function relies on the identification of the initial loading curve from the experimental data.

In the view of the distributed elasto-slide model using an ideal Coulomb slide, the initial loading curve is independent of loading rate such that the maximum force in the initial loading curve responds to the maximum displacement in cyclic loading as seen in Eq. 19. As a result, an initial loading curve can be obtained using a series of experimental hysteresis loops at different amplitudes. An example of the initial loading curve is shown in Fig. 13, the initial loading curve at three different frequencies are obtained from hysteresis cycles of the elastomeric specimen by identifying the force at corresponding maximum displacements. The analytical initial loading curve is determined by considering the influence of the rate-dependent slide. This curve appears elasto-plastic behavior, which is described as:

$$f = \frac{1}{\mathfrak{sp}} \left( 1 - e^{-k\mathfrak{q}\_0 \cdot x} \right) + k\_0 \mathfrak{x} \tag{32}$$

Notably, *<sup>0</sup>* is a distribution constant and an index of the post yield force level, and *k* and *k0* are the stiffness of the leading spring and the remaining polymer stiffness respectively. Summation of the leading spring and the polymer stiffness is just the slope of the forcedisplacement curve when x0. This is conceivable since there only exists the influence of the springs while all slide elements are not yielded. Substituting Eq. 32 into Eq. 31, yields a very simple distribution function as:

$$\text{op}\begin{pmatrix} f^\* \end{pmatrix} = \text{op}\_0 e^{-\text{q}\_0 f^\*} \tag{33}$$

The distribution area for the elastomeric specimen is shown as the shaded area in Fig. 14. It is easily shown that the distribution function satisfies all the properties of Eq. 28.

As the distribution function is determined, for the distributed elasto-slide model using the Coulomb slide, the steady-state forced response of the elastomer under cyclic loading will be predicted as follows:

$$\begin{split} f &= \frac{1}{\mathfrak{op}\_0} \left( 1 + e^{-k\mathfrak{q}\_0 \mathbf{x}\_{\text{max}}} - 2e^{-k\mathfrak{q}\_0 \frac{\mathbf{x} - \mathbf{x}\_{\text{min}}}{2}} \right) + k\_0 \mathbf{x}\_\star \dot{\mathbf{x}} > 0 \\ f &= -\frac{1}{\mathfrak{op}\_0} \left( 1 + e^{-k\mathfrak{q}\_0 \mathbf{x}\_{\text{max}}} - 2e^{-k\mathfrak{q}\_0 \frac{\mathbf{x}\_{\text{max}} - \mathbf{x}}{2}} \right) + k\_0 \mathbf{x}\_\star \dot{\mathbf{x}} < 0 \end{split} \tag{34}$$

where *x*max and *x*min are the maximum and minimum amplitude of the cyclic loading, respectively.

**Figure 13.** Initial Loading Curve

**Figure 14.** Yield Distribution

For the rate-dependent elasto-slide model, the reference velocity *vr* and the exponent *p* need to be determined. The choice of *vr* was based on steady state force-velocity curves in which the boundary between the pre-yield and post-yield region was approximated. Analytically, *p* should be as large as possible such that the post-yield force rapidly transitions to a constant value, which is similar to friction behavior. However, large values of *p* result in a stiffer system. Thus, *p* was chosen by a tradeoff between both factors. For the elastomer specimen 1, the determined model parameters are shown in Table 1, in which the elastomeric specimen has two preload conditions. Notably, the distribution constant *<sup>0</sup>*at 10% preload is lower than that without preload. This implies that the yield force level can be increased with the normal force in the preload condition. Similarly, the preload force also can increase the stiffness. As a result, the addition of a preload perpendicular to the loading axis tends to increase the equivalent stiffness and damping over the entire amplitude range. This effect is due to the compressive preload increasing the friction response of the filler in the elastomer, and is not reflected in the distributed elasto-slide model. Thus, the model parameters are different for different preload conditions. For the elastomer specimen 2, the determined model parameters are shown in Table 2. Notably, the elastomer specimen 2 is much stiffer than the elastomer specimen 1 since the stiffness of the leading spring and remaining spring is much higher. Thus, the loss factor of the elastomer specimen 2 appears much lower than the elastomer specimen 1 though the yield force level of the specimen 2 is higher than the specimen 1.


**Table 1.** Model Parameters for Elastomer Specimen 1

328 Advanced Elastomers – Technology, Properties and Applications

respectively.

**Figure 13.** Initial Loading Curve

100

200

300

Force(N)

400

500

600

**Figure 14.** Yield Distribution

0

0

2.5Hz 5.0Hz 7.5Hz Determined min

0

(34)

0

2

max

2

<sup>0</sup> 0 max

 

*fee kxx*

*x x <sup>k</sup> k x*

<sup>1</sup> 1 2 , 0

*x x <sup>k</sup> k x*

where *x*max and *x*min are the maximum and minimum amplitude of the cyclic loading,

<sup>0</sup> 0.5 <sup>1</sup> 1.5 <sup>2</sup> 2.5 <sup>3</sup> 3.5 <sup>4</sup> 4.5 <sup>5</sup> <sup>0</sup>

Displacement Amplitude(mm)

<sup>1</sup> 1 2 , 0

<sup>0</sup> 0 max

*f ee kxx*


**Table 2.** Model Parameters for Elastomer Specimen 2

### **4.3. Modeling results and validation**

As stated before, the distributed rate-dependent elasto-slide model is reminiscent of the behavior of filler structures in the elastomer such that it can predict the forced harmonic response of an elastomer in the time domain. In this section, single frequency and dual frequency steady-state hysteresis data are used to validate the model. To assess the model's capability in describing elastomer behavior under complex loading conditions, the response under dual frequency loading with slowly varying amplitude is also correlated with model predictions.

For the elastomer specimen 1, three sets of single frequency hysteresis cycle data were used to assess model fidelity. Each set of data was obtained by measuring the forced response while the elastomeric specimen was under sinusoidal displacement excitation at 2.5 Hz, 5.0 Hz and 7.5 Hz, respectively. At each frequency, the displacement amplitude was chosen as 1 mm, 2 mm, 3 mm and 4 mm. In Fig. 15, the experimental data at three frequencies are shown compared to the modeling results. Generally, the modeling results correlate quite well with the experimental results while the displacement amplitude is in the moderate amplitude range, i.e. 2<x<5 mm. In the small amplitude range, i.e. x<2 mm, the analytical model under-predicts the area enclosed by the hysteresis cycle. The reason for that is partly because the lower yield region for the elasto-slide element was replaced with a non-zero constant yield force for numerical consideration and the influence of this approximation was amplified at small deflection loading.

The complex modulus determined by the analytical model is also compared to the experimental result. As shown in solid lines in Fig. 16, the predicted storage and loss stiffnesses using the model have the same amplitude dependent trend as the experimental result. The experimental moduli are well matched with the analytical moduli at moderate amplitude range except that the moduli over small amplitude range are under-predicted especially for loss stiffness. Model predictions are also compared to the experimental data for the loss factor. The predicted loss factor represent common features of the elastomeric response. Clearly, at small amplitude, most of the filler structures, or corresponding elastoslide elements, have not yielded, so that the loss factor is small. As the amplitude increases, breaking filler structures or yielding of slides leads to a rise in the loss factor. After all of the slide elements have yielded, the loss factor decreases again. In Fig. 16, it also shows that both experimental and predicted moduli are weakly dependent on frequency. This phenomenon is consistent with the tribo-elastic mechanism of elastomeric materials [11].

Similar single frequency modeling results for the elastomer specimen 2 are shown in Fig. 17 as force-displacement diagrams for different amplitudes at 2.5 Hz (Fig. 17a), 5 Hz (Fig. 17b) and 7.5 Hz (Fig. 17c), respectively. Clearly, the analytical model captures the amplitude and frequency dependent behavior of the elastomer specimen 2.

In some applications, the elastomer would experience multi-frequency excitation. Under such a circumstance, the potential loss of damping at the lower frequency due to limitation of stroke is well known [5], so it is important to predict the response of the elastomer under dual frequency excitation. Experimental dual frequency force-displacement data of the elastomer specimen 1 were used to evaluate the adaptability of the model under complex loading conditions.

330 Advanced Elastomers – Technology, Properties and Applications

As stated before, the distributed rate-dependent elasto-slide model is reminiscent of the behavior of filler structures in the elastomer such that it can predict the forced harmonic response of an elastomer in the time domain. In this section, single frequency and dual frequency steady-state hysteresis data are used to validate the model. To assess the model's capability in describing elastomer behavior under complex loading conditions, the response under dual frequency loading with slowly varying amplitude is also correlated with model

For the elastomer specimen 1, three sets of single frequency hysteresis cycle data were used to assess model fidelity. Each set of data was obtained by measuring the forced response while the elastomeric specimen was under sinusoidal displacement excitation at 2.5 Hz, 5.0 Hz and 7.5 Hz, respectively. At each frequency, the displacement amplitude was chosen as 1 mm, 2 mm, 3 mm and 4 mm. In Fig. 15, the experimental data at three frequencies are shown compared to the modeling results. Generally, the modeling results correlate quite well with the experimental results while the displacement amplitude is in the moderate amplitude range, i.e. 2<x<5 mm. In the small amplitude range, i.e. x<2 mm, the analytical model under-predicts the area enclosed by the hysteresis cycle. The reason for that is partly because the lower yield region for the elasto-slide element was replaced with a non-zero constant yield force for numerical consideration and the influence of this approximation was

The complex modulus determined by the analytical model is also compared to the experimental result. As shown in solid lines in Fig. 16, the predicted storage and loss stiffnesses using the model have the same amplitude dependent trend as the experimental result. The experimental moduli are well matched with the analytical moduli at moderate amplitude range except that the moduli over small amplitude range are under-predicted especially for loss stiffness. Model predictions are also compared to the experimental data for the loss factor. The predicted loss factor represent common features of the elastomeric response. Clearly, at small amplitude, most of the filler structures, or corresponding elastoslide elements, have not yielded, so that the loss factor is small. As the amplitude increases, breaking filler structures or yielding of slides leads to a rise in the loss factor. After all of the slide elements have yielded, the loss factor decreases again. In Fig. 16, it also shows that both experimental and predicted moduli are weakly dependent on frequency. This phenomenon is consistent with the tribo-elastic mechanism of elastomeric materials [11].

Similar single frequency modeling results for the elastomer specimen 2 are shown in Fig. 17 as force-displacement diagrams for different amplitudes at 2.5 Hz (Fig. 17a), 5 Hz (Fig. 17b) and 7.5 Hz (Fig. 17c), respectively. Clearly, the analytical model captures the amplitude and

In some applications, the elastomer would experience multi-frequency excitation. Under such a circumstance, the potential loss of damping at the lower frequency due to limitation of stroke is well known [5], so it is important to predict the response of the elastomer under

frequency dependent behavior of the elastomer specimen 2.

**4.3. Modeling results and validation** 

amplified at small deflection loading.

predictions.

**Figure 15.** Single Frequency Modeling Results for Specimen 1

**Figure 16.** Complex Modulus for Specimen 1

**Figure 17.** Single Frequency Modeling Results for Specimen 2

(a) In-phase Stiffness

(b) Quadrature Stiffness

(c) Loss Factor

**Figure 16.** Complex Modulus for Specimen 1

The dual frequency test data were obtained while the amplitudes at both 5 and 7.5 Hz frequencies were held constant. For each test condition, the amplitude for 5 Hz and 7.5 Hz frequencies ranged from 0.25 mm to 5 mm, and the sum of both amplitudes must not exceed 5 mm, which corresponds to the maximum allowable strain of 50%. The modeling result at each dual frequency loading condition was correlated with the corresponding experimental result. Some of these dual frequency modeling results are presented in Fig. 18 and 19. The figures are grouped according to the 7.5 frequency amplitude. Fig. 18 shows the modeling results for 2.5 mm amplitude at 7.5 Hz and at four different amplitudes at 5 Hz. As the displacement amplitude at 7.5 Hz is 2.5mm, the experimental dual frequency behavior can be matched quite well with the modeling results. Comparatively, Fig. 19 shows the modeling results for 0.5 mm amplitude at 7.5 Hz and at four different amplitudes at 5 Hz. Notably, the model under-predicts the forced response as the total amplitude at 5 Hz and 7.5 Hz is below 2.5 mm since the high yield force region in Fig. 14 is not well described by the numerical algorithm of the model. In general, the distributed rate-dependent elasto-slide model performs well in the moderate amplitude range except it over-predicts the inner loop in some cases. The model also should be improved to predict the response over the small amplitude range.

**Figure 18.** Dual Frequency Modeling Results for Specimen 1

**Figure 19.** Dual Frequency Modeling Results for Specimen 1

**Figure 18.** Dual Frequency Modeling Results for Specimen 1

amplitude range.

The dual frequency test data were obtained while the amplitudes at both 5 and 7.5 Hz frequencies were held constant. For each test condition, the amplitude for 5 Hz and 7.5 Hz frequencies ranged from 0.25 mm to 5 mm, and the sum of both amplitudes must not exceed 5 mm, which corresponds to the maximum allowable strain of 50%. The modeling result at each dual frequency loading condition was correlated with the corresponding experimental result. Some of these dual frequency modeling results are presented in Fig. 18 and 19. The figures are grouped according to the 7.5 frequency amplitude. Fig. 18 shows the modeling results for 2.5 mm amplitude at 7.5 Hz and at four different amplitudes at 5 Hz. As the displacement amplitude at 7.5 Hz is 2.5mm, the experimental dual frequency behavior can be matched quite well with the modeling results. Comparatively, Fig. 19 shows the modeling results for 0.5 mm amplitude at 7.5 Hz and at four different amplitudes at 5 Hz. Notably, the model under-predicts the forced response as the total amplitude at 5 Hz and 7.5 Hz is below 2.5 mm since the high yield force region in Fig. 14 is not well described by the numerical algorithm of the model. In general, the distributed rate-dependent elasto-slide model performs well in the moderate amplitude range except it over-predicts the inner loop in some cases. The model also should be improved to predict the response over the small

> The behavior of the elastomer under a dual frequency excitation with a slowly-varying amplitude modulated periodic loading can also be predicted using the analytical model. For simplicity, the analytical and experimental simulation results are only shown for one scenario, in which the amplitude for 7.5 Hz is 1.5 mm and the amplitude for 5 Hz is assumed to be as below:

$$A\_{5Hz} = 1.5 \left[ 1 + 0.2 \sin \left( 0.2 \pi t \right) \right] \left( \text{mm} \right) \tag{35}$$

The predicted forced response data shown in Fig. 20a and 20b compared well to the experimental results for two different time scales. Similarly, the modeling forcedisplacement hysteresis cycle is also matched well with the experimental data as shown in Fig. 20c. The predicted damper response due to the slowly varying displacement excitation exhibits the same varying trend as the observed elastomer behavior, and also the force value

is tracked quite well. Clearly, the proposed elastomeric model performs fairly well in predicting dual frequency response, and especially the distributed elasto-slide model can predict the behavior of the elastomer under slowly-varying amplitude modulated periodic loadings.

**Figure 20.** Modeling Results for Dual Frequency with Slowly Varying Amplitude

#### **5. Summary**

Modeling methods for describing elastomeric material behavior were investigated. Most prior models introduced nonlinear terms into the conventional Kelvin model or Zener model. Because filled elastomers are anelastic materials, a friction mechanism damping element proves useful to model rate-independent damping. Nonlinearity in tested elastomeric materials manifested in two ways. First, the forced response of an elastomer subjected to harmonic displacements was nonlinear (non-elliptical), which meant that the response could not be predicted by linear differential or integral equation. Second, the stiffness and damping of elastomers varied as a function of amplitude and frequency. While some models capture the amplitude dependent complex moduli very well using constant parameters, such models cannot predict stress-strain or force-displacement hysteresis accurately. On the other hand, most hysteresis models can predict non-elliptical hysteresis quite well, but their parameters are usually amplitude and frequency dependent. The methods require amplitude and frequency as prior information when these models are implemented.

336 Advanced Elastomers – Technology, Properties and Applications

loadings.

**5. Summary** 

is tracked quite well. Clearly, the proposed elastomeric model performs fairly well in predicting dual frequency response, and especially the distributed elasto-slide model can predict the behavior of the elastomer under slowly-varying amplitude modulated periodic

(a) Force Response (0-9 s) (b) Force Response (7-9 s)

**Figure 20.** Modeling Results for Dual Frequency with Slowly Varying Amplitude

Modeling methods for describing elastomeric material behavior were investigated. Most prior models introduced nonlinear terms into the conventional Kelvin model or Zener model. Because filled elastomers are anelastic materials, a friction mechanism damping element proves useful to model rate-independent damping. Nonlinearity in tested

(c) Force-Displacement Hysteresis

A nonlinear hysteresis model was developed to characterize the nonlinear behavior of the elastomer under a cyclic loading. This model is mechanically analogous to a nonlinear Kelvin model where the stiffness is a nonlinear monotonic function of displacement and the damping is a monotonic rate dependent friction function. Since it accurately describes the characteristics of the hysteresis loops, this model can predict steady state or harmonic forced response very well. However, the model parameters are still amplitude dependent. A challenge still remained to describe transient or stress relaxation behavior using this type of mechanisms-based model.

Therefore, a distributed rate-dependent elasto-slide elastomeric model was used to describe the amplitude dependent characteristics of an elastomer. This physically motivated damper model resembles the behavior of filler structures in the elastomer under cyclic loading. A method to determine the model parameters was presented. It was found that a unique exponential function could be used to describe the yield force distribution for elastomers. Numerical algorithms were developed for model applications. Dynamic test were conducted on a double lap shear elastomeric specimen and a linear concentric tubular elastomeric specimen, respectively, and the measured data were used to evaluate the modeling method. The fidelity of the model was verified by the good correlation between predicted single and dual frequency force-displacement hysteresis and the experimental results except that the damping at lower amplitude range cannot be fully predicted by the model. Since the proposed model is a time domain model, the adaptability of the model in predicting damper response under a slowly varying displacement excitation was evaluated. The predicted force response of the elastomeric specimen under this slowly varying displacement excitation correlated quite well with the corresponding experimental data.

In conclusion, the distributed rate-dependent elasto-slide elastomeric damper model is a time-domain modeling approach to capture nonlinear behavior of the elastomer. The damper model, formulated as a state space model has the advantage that it could easily be implemented into dynamic system models. Because the model is physically motivated, the flexibility in determining the distribution function provides means to improve the model performance especially over the low amplitude range. Although only a one-dimensional

elastomeric model is described in this paper, the distributed elasto-slide model can also be extended into a three-dimensional form such that it can be implemented easily into a finite element analysis for a complex elastomeric damper configuration.

### **Author details**

Wei Hu and Norman M. Wereley\*

*Department of Aerospace Engineering, University of Maryland, College Park, MD, USA* 

### **6. References**


<sup>\*</sup> Corresponding Author

[12] Strganac, T.W., "An Experiment and Analytical Methodology to Characterize Nonlinear Elastomeric Lag Dampers," Proceeding of the 38th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Material Conference, Orlando, Florida, April 1997.

338 Advanced Elastomers – Technology, Properties and Applications

**Author details** 

**6. References** 

516.

1948.

53.

617.

Corresponding Author

 \*

Wei Hu and Norman M. Wereley\*

No. 4, 1995, pp. 660-670.

N.Y., 1933, pp. 679-680.

element analysis for a complex elastomeric damper configuration.

elastomeric model is described in this paper, the distributed elasto-slide model can also be extended into a three-dimensional form such that it can be implemented easily into a finite

[1] Coveney, V.A., Johnson D.E., and Turner D.M., "A Triboelastic Model for The Cyclic Mechanical Behavior of Filled Vulcanizates," *Rubber Chemistry and Technology*, Vol. 68,

[2] Panda, B., Mychalowycz, E. and Tarzanin, F.J., "Application of Passive Dampers to Modern Helicopters," *Smart Materials and Structures*, Vol. 5, No. 5, 1996, pp. 509-

[3] Zener, C.M., *Elasticity and Anelasticity of Metals*, University of Chicago Press, Chicago, IL.,

[4] Gandhi, F., and Chopra, I., "A Time-domain Non-linear Viscoelastic Damper Model,"

[5] Felker, F., Lau, B., McLaughlin, S. and Johnson, W., "Nonlinear Behavior of an Elastomeric Lag Damper undergoing Dual-Frequency Motion and its Effect on Rotor Dynamics," *Journal of the American Helicopter Society*, Vol. 32, No. 4, 1987, pp. 45-

[6] Payne, A.R. and Whittaker, R.E., "Low Strain Dynamic Properties of Filled Rubbers,"

[7] Tarzanin F.J., and Panda, B., "Development and Application of Nonlinear Elastomeric and Hydraulic Lag Damper Models," Proceedings of the 36th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials

[8] Timoshenko, S.P., *Strength of Materials*, Part 2, D. Van Nostrand Company, New York,

[9] Iwan, W.D., "A Distributed-Element Model for Hysteresis and Its Steady-State Dynamic

[10] Iwan, W.D., "On a Class of Models for the Yielding Behavior of Continuous and Composite Systems," *Journal of Applied Mechanics*, Vol. 34, September 1967, pp. 612-

[11] Turner, D.M., "A Triboelastic Model for the Mechanical Behavior of Rubber," *Plastics* 

Response," *Journal of Applied Mechanics*, Vol. 33, No. 4, 1966, pp. 893-900.

*and Rubber, Processing and Applications*, Vol.9, No. 4, 1988, pp. 197-201.

*Smart Materials and Structures*, Vol. 5, No. 5, 1996, pp. 517-528.

*Rubber Chemistry and Technology*, Vol. 44, No. 2, 1971, pp. 440-478.

Conference, New Orleans, Louisiana, April 1995.

*Department of Aerospace Engineering, University of Maryland, College Park, MD, USA* 


	- [25] Faires, J.D. and Burden, R., *Numerical Methods*, Brooks/Cole Publishing Company, Pacific Grove, California, 1998.

**Chapter 14** 

## **Modelling Friction and Abrasive Wear of Elastomers**

Nándor Békési

340 Advanced Elastomers – Technology, Properties and Applications

Pacific Grove, California, 1998.

[25] Faires, J.D. and Burden, R., *Numerical Methods*, Brooks/Cole Publishing Company,

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50498

### **1. Introduction**

Tribological modelling of elastomer parts has high importance in engineering practice due to their widespread industrial use. Generally, the experimental investigation of the wear behaviour is time consuming and expensive, which led to the development of numerical techniques.

In order to properly model the friction and wear of elastomer materials, it is essential to understand some specific mechanical and tribological features. Mechanical characteristics of elastomer materials are different from those of metals and other relatively hard materials. Elastomers can bear large deformations, the stress-strain curves are non-linear, time and temperature dependent. When the elastomer is deformed, a part of the strain energy is dissipated. This feature makes the elastomer parts fulfil damping functions. The contact mechanical aspects of elastomers are also important for modelling. When an elastomer part is pressed against a counter-surface, the contact area would be larger than that estimated by the Hertzian theory, because of the effect of adhesion [1]. For the same reason, when the elastomer is contacting a rough surface, the real area of contact would be bigger than it is common for harder materials, as the elastomer tends to fill the valleys of the rough surface.

Friction of elastomers is usually considered as a joint effect of the friction on the surface by adhesion (and fluid film shearing in lubricated friction) and the internal friction generated by the hysteresis of the deformation by the counterpart [2]. Since the hysteretic contribution is caused by the visco-elastic properties of the material, it depends on the temperature and on the exciting frequency. When the elastomer is sliding on a hard, rough counter-surface, the asperities of the rough surface repeatedly deform the elastomer and so cause hysteresis loss. Thus the hysteresis contribution is dependent on the sliding speed. The hysteresis caused by the surface roughness of the counterpart is usually referred as micro-hysteresis. Hysteresis can occur on the macro level as well, where it can be a significant part of the friction.

© 2012 Békési, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

When designing elastomer parts for tribological applications, modelling friction is inevitable in order to determine the contact conditions, which usually have an important effect on the functionality of the part. Wear, on the other hand, has less importance in short term use (an average passenger car tire can run up to 50 000 km before the need to change). However, in most cases, the reason for the elastomer parts to be replaced is the wear. In long term or when the lubrication is not ideal wear can dramatically change the geometry of the part and thus the nature of the contact. Therefore wear modelling is needed to estimate the product lifespan and to schedule maintenance operations. Beside experimental tests, which are cost and time consuming, numerical wear simulation techniques were developed in the last decades and became more and more commonly used. Although there are wear simulation methods incorporating various numerical methods, such as the boundary element method [3] or the discrete element method [4], the most successful and popular approaches deal with the finite element method, since it is a general method for mechanical stress analysis.

The most popular FE based wear simulation method was developed by Pödra and Andersson [5]. In this iterative method, firstly the contact pressure distribution is determined by FEM and the nodal wear increments are calculated from the wear equation of Archard [6]. Then the contacting nodes are moved with respect to the nodal wear values. Finally the FE contact calculation is carried out again with the modified mesh and the cycle is repeated according to the simulated wear process. This method was improved by several researchers (e.g. [7,8]) for various applications, different geometries and materials, such as metals, ceramics, polymers and composites. Kónya and Váradi in [9] improved the method in order to consider heat generation and time dependent material properties in the wear simulation.

Even though this chapter shows methods for modelling and simulation of elastomers, some experimental work is inevitable for the calculations. The parameters of the material models are usually determined based on experimental data. In my simulations I used the material model described by Pálfi et al. in [10]. It is a Mooney-Rivlin model with an attached generalized Maxwell model. This material model was constructed based on stress relaxation measurements. The tribological properties, like the *coefficient of friction* and the *specific wear rate* also require experimental tests. The basic mechanical and tribological properties of the investigated materials are summarised in [11] and in [12]. The coefficient of friction depends on many factors, yet the lubrication has one of the most important effects. Figure 1 shows the result of a *starvation* test, which was carried out at the Institute for Composite Materials, Kaiserslautern, Germany on a shaft-on-plate test rig (Dr. Tillwich GmbH, Horb-Ahldorf, Germany). In this configuration the rubber sheet was pressed against the rotating steel shaft. The normal and the tangential force were measured online. In this experiment lubricant (Hydraulan 407-1, BASF, Ludwigshafen, Germany) was applied only at the beginning of the test. Later, as the amount of lubricant decreased between the rubber and the steel shaft, the nature of the lubrication changed. In the beginning, a mostly constant coefficient of friction with a value of µ = 0.06 lets us assume *soft EHD* lubrication regime (the green section in the figure) [13]. In the next two phases, the coefficient of friction is increasing, which means more and more solid body contacts as the lubricant decreases (blue and orange). In the last phase, the coefficient of friction has an approximately constant value that is higher than one, which implies technically dry friction (red).

The presented finite element (FE) models were generated with MSC.Marc Mentat software package; the calculations were done by the MSC.Marc Solver. The developed wear simulation algorithms operate the MSC.Marc model and result files, however they are independent of the FE environment and can easily be adapted to other FE systems.

**Figure 1.** The coefficient of friction with starving lubrication (v = 250 mm/s Fn = 28 N)

### **2. Modelling internal friction of elastomers**

During deformation of elastomer materials, because of the visco-elastic material behaviour, a part of the applied strain energy is transformed to heat as a result of internal friction, hysteresis [14]. The amount of this energy loss can be characterized by the loss factor (tan(δ)), which is the ratio of the loss and the storage moduli of the material. When repeated loads are present, the hysteresis contribution to the friction is more significant. In case of sliding friction between elastomer and a rough rigid counter-surface, the elastomer is subjected to repeated, cyclic deformation by the asperities of the rough surface. This kind of internal heat generation can lead to fatigue wear [15].

### **2.1. Modelling micro-hysteresis**

342 Advanced Elastomers – Technology, Properties and Applications

general method for mechanical stress analysis.

simulation.

When designing elastomer parts for tribological applications, modelling friction is inevitable in order to determine the contact conditions, which usually have an important effect on the functionality of the part. Wear, on the other hand, has less importance in short term use (an average passenger car tire can run up to 50 000 km before the need to change). However, in most cases, the reason for the elastomer parts to be replaced is the wear. In long term or when the lubrication is not ideal wear can dramatically change the geometry of the part and thus the nature of the contact. Therefore wear modelling is needed to estimate the product lifespan and to schedule maintenance operations. Beside experimental tests, which are cost and time consuming, numerical wear simulation techniques were developed in the last decades and became more and more commonly used. Although there are wear simulation methods incorporating various numerical methods, such as the boundary element method [3] or the discrete element method [4], the most successful and popular approaches deal with the finite element method, since it is a

The most popular FE based wear simulation method was developed by Pödra and Andersson [5]. In this iterative method, firstly the contact pressure distribution is determined by FEM and the nodal wear increments are calculated from the wear equation of Archard [6]. Then the contacting nodes are moved with respect to the nodal wear values. Finally the FE contact calculation is carried out again with the modified mesh and the cycle is repeated according to the simulated wear process. This method was improved by several researchers (e.g. [7,8]) for various applications, different geometries and materials, such as metals, ceramics, polymers and composites. Kónya and Váradi in [9] improved the method in order to consider heat generation and time dependent material properties in the wear

Even though this chapter shows methods for modelling and simulation of elastomers, some experimental work is inevitable for the calculations. The parameters of the material models are usually determined based on experimental data. In my simulations I used the material model described by Pálfi et al. in [10]. It is a Mooney-Rivlin model with an attached generalized Maxwell model. This material model was constructed based on stress relaxation measurements. The tribological properties, like the *coefficient of friction* and the *specific wear rate* also require experimental tests. The basic mechanical and tribological properties of the investigated materials are summarised in [11] and in [12]. The coefficient of friction depends on many factors, yet the lubrication has one of the most important effects. Figure 1 shows the result of a *starvation* test, which was carried out at the Institute for Composite Materials, Kaiserslautern, Germany on a shaft-on-plate test rig (Dr. Tillwich GmbH, Horb-Ahldorf, Germany). In this configuration the rubber sheet was pressed against the rotating steel shaft. The normal and the tangential force were measured online. In this experiment lubricant (Hydraulan 407-1, BASF, Ludwigshafen, Germany) was applied only at the beginning of the test. Later, as the amount of lubricant decreased between the rubber and the steel shaft, the nature of the lubrication changed. In the beginning, a mostly constant coefficient of friction with a value of µ = 0.06 lets us assume *soft EHD* lubrication regime (the green section in the figure) [13]. In the next two phases, the coefficient of friction is increasing, which means

In a series of numerical calculations the micro-hysteresis was investigated. First the hysteretic friction caused by a single asperity was studied. The asperity was modelled by a rigid cylinder with a radius of 250 µm (Figure 2). The asperity was pressed and rubbed against a rubber block with varying sliding speed. In this model the surface friction was neglected, so only the friction caused by the hysteresis was considered.

**Figure 2.** The model of one cylindrical asperity sliding on rubber

Figure 3 shows the results of the numerical simulations. The coefficient of friction was calculated as the ratio of the tangential and the normal reaction forces. It can be seen that the velocity dependence of the hysteretic coefficient of friction is in agreement with the measured loss factor curve.

**Figure 3.** The hysteretic coefficient of friction vs. the exciting frequency and the corresponding loss factor curve

#### **2.2. Modelling macro-hysteresis**

A 3D finite element model was created in order to model a Pin-on-Plate test described in reference [12]. A 4 mm thick rubber sheet was modelled by 7410 3D brick elements (Figure 4). The hemispherical steel slider (r = 5 mm) was modelled as a rigid body. The rubber sheet was secured on a steel substrate, which was modelled as a rigid body too. The assumption of rigid bodies is reasonable if one considers the fact that the steel has four orders of magnitude higher elastic modulus than the rubber specimen.

**Figure 2.** The model of one cylindrical asperity sliding on rubber

measured loss factor curve.

factor curve

**2.2. Modelling macro-hysteresis** 

0.00

0.05

0.10

Coefficient of Friction

0.15

0.20

magnitude higher elastic modulus than the rubber specimen.

Figure 3 shows the results of the numerical simulations. The coefficient of friction was calculated as the ratio of the tangential and the normal reaction forces. It can be seen that the velocity dependence of the hysteretic coefficient of friction is in agreement with the

CoF tan(δ)

**Figure 3.** The hysteretic coefficient of friction vs. the exciting frequency and the corresponding loss

1.E-06 1.E-03 1.E+00 1.E+03 1.E+06 1.E+09 1.E+12 1.E+15

0

0.2

0.4

Loss Factor (tan(δ))

0.6

0.8

Exciting Frequency [Hz]

A 3D finite element model was created in order to model a Pin-on-Plate test described in reference [12]. A 4 mm thick rubber sheet was modelled by 7410 3D brick elements (Figure 4). The hemispherical steel slider (r = 5 mm) was modelled as a rigid body. The rubber sheet was secured on a steel substrate, which was modelled as a rigid body too. The assumption of rigid bodies is reasonable if one considers the fact that the steel has four orders of

**Figure 4.** 3D finite element model of a Pin-on-Plate configuration (the elements are blue; the surfaces considered to be rigid are magenta). Adaptive local remeshing was applied in the contact area.

The length of the elements varied between 0.5 and 2 mm. The hemispherical counter surface was moving along a circular track of 33 mm diameter at various sliding velocities. Fn = 100 N normal load acted on the substrate, which pressed the rubber sheet against the slider.

No displacement was possible between the substrate and the specimen; the coefficient of friction specified between the specimen and the counter surface was µ = 0.06 in order to take into consideration the friction developing on the lubricated rubbing surfaces (See Figure 1).

A series of calculations were run on this FE model with different sliding velocities (in a range between 3 and 1036 mm/s).

After the calculations, the friction corresponding to each velocity value was determined from the finite element result files. Changes in the coefficient of friction are shown in Figure 5. The figure also indicates the corresponding measured results.

A slight decrease in coefficient of friction can be observed both in the calculation and in the test results. Note that these calculations did not consider the changing viscosity of the lubricant. Any changing of the friction can be attributed to the changing of the internal friction, since the surface friction was modelled by a constant prescribed coefficient of friction. The decreasing friction can be explained by the decreasing penetration depth of the slider and the decreasing contact area. The higher the sliding velocity, the lower the penetration depth, because of the visco-elastic behaviour of the material (the rubber seems harder at higher strain rates). This causes decreasing in the excited volume, and the hysteresis part of the friction depends highly on the excited volume. Of course the hysteresis – as it was already mentioned – depends on the exciting frequency as well, however in the investigated sliding velocity range the resultant frequency (approximately from 0.3 Hz to 103.6 Hz) does not change the tan(δ) and the coefficient of friction significantly.

Figure 6 shows the contact area and the penetration depth of the counter surface for the sliding velocities examined. It can be established that the counter surface moving at a higher velocity can be pressed into the rubber surface to a lesser extent; thereby the contact area will also be smaller. There is a 0.25 mm difference between the penetration values corresponding to the lowest and the highest velocity.

**Figure 5.** Coefficient of friction as a function of sliding velocity

**Figure 6.** Contact area and penetration depth of the counter surface in function of sliding velocity

### **3. Modelling wear of elastomers**

Three different wear simulation approaches are presented in the followings. The first one is the method of moving the contact nodes according to the nodal wear increments as proposed by Pödra and Andersson in [5]. It will be shown, that this method is suitable for modelling wear, which is much smaller than the elements in the FE model. The second and the third methods can model relatively large wear by global remeshing and by deactivation of certain elements, respectively.

#### **3.1. Wear simulation by moving the nodes**

The wear simulation method of Pödra and Andersson was adapted to model the wear process of an elastomer part considering its non-linear material properties. The simulation process is shown in Figure 7. In the modelled configuration a rubber sheet (thickness: 2 mm) was pressed against a rotating steel cylinder (diameter: 10 mm) with a normal force of 16 N. The rubber part was modelled in 2D assuming plain strain conditions (Figure 8). The rotating steel cylinder was modelled as a rigid body. A coefficient of friction with a value of 1 was prescribed between the rubber and the cylinder, according to the dry wear test results.

After the FE contact calculation, the wear increments were determined based on the stress distribution. The *Δhi* nodal wear of the *i*-th node in the contact area by Archard [16] is defined as:

$$
\Delta \mathbf{t}\_i = \mathbf{V}\_s p\_i \upsilon \Delta t \tag{1}
$$

where

346 Advanced Elastomers – Technology, Properties and Applications

corresponding to the lowest and the highest velocity.

**Figure 5.** Coefficient of friction as a function of sliding velocity

0

0

5

10

**Contact Area [mm^2]**

15

20

25

0.05

0.1

**Coefficient of Friction**

0.15

0.2

**Figure 6.** Contact area and penetration depth of the counter surface in function of sliding velocity

0 0.2 0.4 0.6 0.8 1 **Sliding velocity [m/s]**

0

0.2

0.4

**Penetration depth [mm]**

0.6

0.8

1

Figure 6 shows the contact area and the penetration depth of the counter surface for the sliding velocities examined. It can be established that the counter surface moving at a higher velocity can be pressed into the rubber surface to a lesser extent; thereby the contact area will also be smaller. There is a 0.25 mm difference between the penetration values

Measured Calculated Prescribed CoF

0 0.2 0.4 0.6 0.8 1 **Sliding velocity [m/s]**

Contact Area Penetration

*Ws* is the specific wear rate,

*pi* is the nodal contact pressure at node *i*,

*v* is the sliding velocity,

*Δt* is the time increment.

In the next step of the simulation the contacting nodes were moved with respect to the nodal wear increments. Then the cycle was restarted with a new contact calculation using the modified geometry.

Figure 9 shows the result of the wear simulation after one minute of sliding. The profile of the worn specimen is also plotted in the figure. It can be seen that the wear is not symmetric. The simulation and the experimental results are in good agreement both by the means of the amount of wear and the worn shape. However, in longer term simulations and with finer FE mesh, this method failed to accurately estimate the experimental results.

**Figure 7.** The flowchart of the wear simulation

**Figure 8.** The finite element mesh of the shaft-on-plate configuration

**Figure 9.** The calculated and measured worn profile of the rubber sheet (t = 60 s, FN = 16 N, ω = 50 rad/s)

Despite the widespread use of the Pödra-Andersson method, the wear simulation is highly limited; only the top layer of elements can be worn. This limitation does not affect the usability of the method in case of relatively hard materials such as metals or ceramics, since minimal wear can have a great effect on the contact area and pressure distribution and thus on the performance of these parts (e.g. bearings). In case of rubber parts, wear in the magnitude of some µm does not change the contact conditions significantly. The wear needed for a rubber part to malfunction is much greater than those of metal parts; therefore the wear of the top layer is insufficient in rubber applications. To increase the volume of the wear to be modelled, the size of the elements can be increased; however it is not recommended, since it decreases the accuracy of the calculation. In the followings two methods are proposed to simulate wear regardless of the element size. It will be shown that these methods can model the material loss due to wear even in a scale comparable with the size of the part.

#### **3.2. Wear simulation by global remeshing**

348 Advanced Elastomers – Technology, Properties and Applications

Initial geometry

FE contact calculation

Wear calculation

Contact Pressure Distribution

Specific wear rate

Nodal wear depth

New geometry

t < tmax ? End Yes No

Moving the nodes

**Figure 7.** The flowchart of the wear simulation




**Wear depth [mm]**


0

**Figure 8.** The finite element mesh of the shaft-on-plate configuration

**Figure 9.** The calculated and measured worn profile of the rubber sheet (t = 60 s, FN = 16 N, ω = 50 rad/s)


Measurement Simulation

In order to model wear that is larger than the elements of the FE mesh, a wear simulation procedure was developed using global remeshing. To demonstrate this method a wear simulation study of a reciprocating sliding seal (Figure 10) is presented. In the investigated application a rubber seal is coupled with an aluminium rod (diameter: ø22.2 mm). Since this method integrates the FE mesh generation, the geometric model is the starting point of the wear simulation. The seal was modelled assuming axisymmetry.

The counter-surface was modelled as rigid body, since its elastic modulus is some orders of magnitude larger than the elastomer material of the seal. The seal was mounted axially by fitting it in the housing, which was considered ideally rigid. The rod was moving with alternating motion at a speed of 20 mm/s with 9 mm amplitude. The contour of the seal section was modelled by line segments (Figure 11). A fine approximation in the vicinity of the lip as well as the ridges was required; in the region of the ridges and the lip the average length of the line segments was 10 µm.

The seal section was meshed by the built-in automatic meshing procedure using three-node axisymmetric triangular finite elements (Figure 12). In the FE calculations the interference fit between the seal and the rod was also taken into consideration: the unloaded inner diameter of the seal was 1.8 mm smaller than the diameter of the rod. This caused the seal to be stretched in the first load step. In the next step the alternating motion was modelled considering the friction during the wear simulation. According to the working conditions, in the inward strokes of the rod a pressure load of p = 2 MPa was applied on the right side of the seal, modelling the sealed pressure. There was no pressure applied in the outward strokes. Between the seal and the rod a prescribed coefficient of friction with a value of µ = 0.1 was defined in order to model the lubricated friction.

In the wear simulation process, the FE mesh was created based on the line segments describing the geometry. The automatic mesher was set to create elements of different sizes, so the elements in the vicinity of the contact zone were smaller than those deep inside the material. In

the initial mesh approximately 4000 elements were distributed in a way that the regions, that are important for the wear calculations, were modelled more accurately with the small elements. Even so the total number of the elements remained in a reasonably low range, so the CPU time for one calculation was short enough to handle the iterative simulation process.

**Figure 10.** Section view of the investigated reciprocating seal

**Figure 11.** The geometry of the seal consisting of line segments and points

**Figure 10.** Section view of the investigated reciprocating seal

**Figure 11.** The geometry of the seal consisting of line segments and points

the initial mesh approximately 4000 elements were distributed in a way that the regions, that are important for the wear calculations, were modelled more accurately with the small elements. Even so the total number of the elements remained in a reasonably low range, so the CPU time for one calculation was short enough to handle the iterative simulation process.

**Figure 12.** Deformed shape of the seal in the first outward stroke (without pressure). The curves represent the rigid bodies (the rod and the housing). The element size varies from 10 to 200 µm

The flowchart of the simulation can be seen in Figure 13. In the first step the FE mesh was created based on the initial geometry. After adjusting the simulation parameters (coefficient of friction, sliding velocity, material properties, applied pressure, time increment, etc.) the contact calculation was run. In the first stroke the rod was moving rightward. The wear calculation was performed using Equation 1. The moving of the contacting nodes according to the nodal wear increments was the next step, however in a slightly different way. The surface nodes were attached to the contour points of the seal, so if a point was moved the attached node would also move with respect to the nodal wear value of the attached node. After moving the points, the whole FE mesh was deleted. Based on the new geometry the automatic mesher created the new elements and the cycle ran over again. The direction of the rod motion was changed in each cycle, so one simulation cycle represented one outward (left) or one inward (right) stroke. The applied pressure was also changed stroke by stroke, as for the inward strokes the working pressure of 2 MPa was applied, while in the outward strokes the working pressure was turned off in the simulation process. 200 simulation cycles were calculated in the frame of this study.

Figure 14 shows the deformed shape of the seal in the first outward stroke, when the seal was not pressurized. The resulting contact pressure comes from the stretching of the seal on the rod. It can be seen that the contact area is reduced: only five ridges on the left side, one on the right side and the lip are in contact with the rod. The maximum of the contact pressure can be observed at the lip. The values of the contact pressure are generally low, since there is no working pressure applied to press the seal surface to the rod.

**Figure 13.** The flowchart of the wear simulation process

The calculated worn profiles of the seal are shown in Figure 15 after 0, 100 and 200 cycles of the simulation. It can be seen that the left side ridges of the seal and the lip edge wear first. Figure 15 also shows the contact pressure distributions in the investigated wear phases. Note that the FE mesh is different in each state, since the model was remeshed in every cycle; however the element sizes are distributed identically. It should be mentioned that the total simulated wear depth reaches 30 µm in some regions (e.g. at the lip), which is three times higher than the size of the elements in the contact area. This wear depth is so large that the Pödra-Andersson method can not be applied here.

Basic simulation parameters

**Figure 13.** The flowchart of the wear simulation process

that the Pödra-Andersson method can not be applied here.

The calculated worn profiles of the seal are shown in Figure 15 after 0, 100 and 200 cycles of the simulation. It can be seen that the left side ridges of the seal and the lip edge wear first. Figure 15 also shows the contact pressure distributions in the investigated wear phases. Note that the FE mesh is different in each state, since the model was remeshed in every cycle; however the element sizes are distributed identically. It should be mentioned that the total simulated wear depth reaches 30 µm in some regions (e.g. at the lip), which is three times higher than the size of the elements in the contact area. This wear depth is so large

End

Yes

t ≥ t max ?

No

FE frictional contact calculation

FE meshing

Start

Initial geometry

Wear calculation

Nodal wear depths

Stress distribution

Moving the points

New geometry

**Figure 14.** Deformed shape of the seal in the first outward stroke (without applied pressure) and its contact pressure distribution.

**Figure 15.** The calculated worn profiles after 0, 100 and 200 cycles of the simulation in inward strokes (p = 2 MPa) and the corresponding contact pressure distributions along the axis

Based on the change of contact pressure distribution the followings can be established. The area of contact increased, while the maximum values of the contact pressure decreased. This phenomenon can also be seen in Figure 16. for the contact pressure distribution of the lip over the wear process. One can see that the seal at the beginning has about three times higher contact pressure than later. The pressure reduction is caused by the wear process, which also increases the contact area. It is remarkable that at the left side of the lip the contact pressure is increasing at first, and then slightly decreasing. It can be explained by the wear of the lip. As the lip has the highest contact pressure in the beginning of the simulation, it will suffer the most severe wear in the first cycles. As the lip wears, the contact pressure distribution becomes flatter, which causes the slight increase in the lower pressure region. Later this mostly uniform contact pressure remains, only the average pressure decreases.

**Figure 16.** The change of the contact area and the contact pressure distribution at the lip over the wear process (after 0, 50 and 100 simulation cycles)

#### **3.3. Wear simulation by deactivation of elements**

In hydraulic systems the seals are well lubricated in most cases, so those parts show minimal wear and they can function for years without any sign of damage. However, in some special cases the rubber seals can work under mixed or boundary lubrication conditions. It is suspected that sometimes (typically after the machine is stopped for a longer period) the seal can squeeze the lubricant out of the space between the surfaces and some parts come to direct contact. If it happens, adhesion can occur, which significantly increases the frictional force and the wear mechanism would change completely. The different wear mechanisms of rubbers are examined by several studies e.g. in [17, 18]. In general it can be said that dry sliding wear of rubber can be traced back to the initiation and propagation of cracks in the rubber material. The modelling of this type of wear requires advanced simulation techniques.

A popular approach for modelling material detachment is the deactivation or deleting elements from the FE mesh (sometimes also referred as element death) based on the characteristic damage criterion of the given material. It is commonly used for simulation of manufacturing processes e.g. in [19, 20]. The principle of wear – from the simulation point of view – is similar to the manufacturing processes. Authors in [21] modelled the particle detachment process based on the strain history of the elements in dry sliding.

The wear process of an experimental test of the sliding seal described above was modelled. The seal was cut and a section of it was straightened and fixed in the clamping. The aluminium rod was then pressed and rubbed against the surface of the seal. The detailed description of the test can be found in [22]. In the series of wear tests, different lubrication conditions were studied, and it was found that when the amount of lubricant was decreased, the seal surface worn rather intensively. The surface of the seal was inspected with a MicroProf white light profilometry device (WLP, Fries Research & Technology GmbH., Bergisch Gladbach, Germany) before and after the test, so the wear of the seal can be evaluated and compared with the simulation results (See Figure 17). One can see that the most wear occurs in the symmetry plane of the seal, where the contact pressure was maximal; this region shows a rough surface texture, which is probably caused by the wear particles that were not removed from the contact region, since there are peaks that are higher than the original surface. In order to avoid the false results caused by the wear debris, a plane with slightly less but still significant wear was investigated, which was clear of wear debris (indicated in Figure 17).

**Figure 17.** Worn surface of the seal

354 Advanced Elastomers – Technology, Properties and Applications

decreases.

Based on the change of contact pressure distribution the followings can be established. The area of contact increased, while the maximum values of the contact pressure decreased. This phenomenon can also be seen in Figure 16. for the contact pressure distribution of the lip over the wear process. One can see that the seal at the beginning has about three times higher contact pressure than later. The pressure reduction is caused by the wear process, which also increases the contact area. It is remarkable that at the left side of the lip the contact pressure is increasing at first, and then slightly decreasing. It can be explained by the wear of the lip. As the lip has the highest contact pressure in the beginning of the simulation, it will suffer the most severe wear in the first cycles. As the lip wears, the contact pressure distribution becomes flatter, which causes the slight increase in the lower pressure region. Later this mostly uniform contact pressure remains, only the average pressure

0 50 100

**Figure 16.** The change of the contact area and the contact pressure distribution at the lip over the wear

2.5 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 x direction [mm]

In hydraulic systems the seals are well lubricated in most cases, so those parts show minimal wear and they can function for years without any sign of damage. However, in some special cases the rubber seals can work under mixed or boundary lubrication conditions. It is suspected that sometimes (typically after the machine is stopped for a longer period) the seal can squeeze the lubricant out of the space between the surfaces and some parts come to direct contact. If it happens, adhesion can occur, which significantly increases the frictional force and

process (after 0, 50 and 100 simulation cycles)

0

1

2

3

4

Contact pressure [MPa]

5

6

7

**3.3. Wear simulation by deactivation of elements** 

The previous wear simulation methods were not suitable to model the compound wear mechanism, which was caused by boundary lubrication. Beside the wear simulation method of moving the nodes, the deactivation of elements was integrated in the simulation procedure. The flowchart of the simulation can be seen in Figure 18.

**Figure 18.** Flowchart of the simulation [22]

The test can be divided into repeated cycles of the inward (right) and outward (left) strokes. Since the stress distributions in the inward and outward strokes are different, the wear calculations were performed four times in every cycle (see Figure 19).

Starting with the initial geometry a frictional contact calculation is carried out assuming inward stroke. Based on the results, the normal stress of each element in the model was compared to the predefined critical stress limit, the tensile strength, which was 19.0 MPa for the material of the seal. If the tensile stress of an element was greater than the critical value, the element became deactivated. The damage check was done in every element in the following way. If *σi* > *σcrit* (where *σi* is the normal stress in the *i*-th element, and *σcrit* is the ultimate tensile strength of the material), the element *i* became deactivated.

356 Advanced Elastomers – Technology, Properties and Applications

Specific wear rate

Initial geometry

**Figure 18.** Flowchart of the simulation [22]

procedure. The flowchart of the simulation can be seen in Figure 18.

The previous wear simulation methods were not suitable to model the compound wear mechanism, which was caused by boundary lubrication. Beside the wear simulation method of moving the nodes, the deactivation of elements was integrated in the simulation

FE contact calculation

Stress distribution

Damage check *σ<sup>i</sup>* > *σcrit ?*

No Yes

Start

Wear calculation

Nodal wear depths

Deactivation of elements

Moving the nodes

New geometry

The test can be divided into repeated cycles of the inward (right) and outward (left) strokes. Since the stress distributions in the inward and outward strokes are different, the wear

End

Yes

t ≥ t max ?

No

calculations were performed four times in every cycle (see Figure 19).

**Figure 19.** The motion of the rod (continuous red) and the points of the wear calculations in a cycle (blue dots)

Wear calculation was applied to determine the nodal wear of the nodes in the contact area using Equation (1). Once the nodal wear values were known, the nodes, which were in contact, were moved to their new position as in the Pödra-Andersson technique.

The seal section was modelled by four-node 2D finite elements assuming plain strain condition. The housing, the clamping, which fixed the seal to the housing, and the counterpart, was considered ideally rigid. In the first load step of the FE calculation, the seal was mounted in the holder, and the normal load was applied to the counterpart, according to the wear test. After this, the rod started its motion. There was a prescribed coefficient of friction between the seal and the rod with a value of µlub = 0.3, representing the boundary lubricated friction; and another between the seal and the holder parts with a value of µdry = 1, since the friction condition was dry in this case. Values of coefficient of friction in case of lubricated and dry friction were determined based on the results in [12] and [11], respectively. Note that the lubrication was taken into consideration only by the lower value of coefficient of friction. Other effects – like the lubricant separating the surfaces reducing real contact and wear – were neglected. The FE mesh with the ideally rigid parts and with the applied prescribed coefficients of friction is shown in Figure 20.

The calculated worn profiles are shown in Table 1. It can be seen that the ridges of the seal and the lip edge wear first, which is in accordance with the experimental findings. One can see that during the wear simulation the wear occurs not only in the top layer, but in several

layers of elements. The maximal wear depth is about ten times bigger than the typical element edge length in the contact area.

**Figure 20.** FE model of the seal section before (left) and after mounting in the holder (right)

The height of the elements in contact gets reduced due to wear and over time the contacting elements become more and more distorted and can even turn inside out, which can make the calculations unable to run. To avoid this and to maintain the stability of the simulation process, the following procedure was implemented. A filter was applied in the simulation algorithm, which checked the distortion of the elements before each FE calculation. Distortion in the FE software is determined at each node of the element as the cosine of the angle between the element edges at the nodes. The distortion of the element is defined as the maximum of these cosines. If the normal vector to the element changes sign from one node to the next (inside-out element), the distortion is one. If the distortion value of an element exceeded a critical value (0.99) before the load was applied, the element got deactivated. With this method the inside-out elements and the elements that were too flattened due to the wear, were became deactivated.

The wear process of one ridge can be tracked in the successive images of Figure 21. One can see the effect of the wear algorithm on the height of the elements. The dark shaded elements were deactivated by the damage check module of the simulation.

The measured and the calculated worn profiles of the seal are compared in Figure 22. The calculated worn profiles are mostly identical with the measured worn surfaces. However, in the vicinity of the lip, the simulation over predicts the wear. It is suspected that the deviation is caused by the assumption of linear wear theory in the wear calculations, while in reality wear is a non-linear function of contact pressure. This simplification nevertheless appears to be suitable for the regions with less contact pressure.

element edge length in the contact area.

the wear, were became deactivated.

layers of elements. The maximal wear depth is about ten times bigger than the typical

**Figure 20.** FE model of the seal section before (left) and after mounting in the holder (right)

The height of the elements in contact gets reduced due to wear and over time the contacting elements become more and more distorted and can even turn inside out, which can make the calculations unable to run. To avoid this and to maintain the stability of the simulation process, the following procedure was implemented. A filter was applied in the simulation algorithm, which checked the distortion of the elements before each FE calculation. Distortion in the FE software is determined at each node of the element as the cosine of the angle between the element edges at the nodes. The distortion of the element is defined as the maximum of these cosines. If the normal vector to the element changes sign from one node to the next (inside-out element), the distortion is one. If the distortion value of an element exceeded a critical value (0.99) before the load was applied, the element got deactivated. With this method the inside-out elements and the elements that were too flattened due to

The wear process of one ridge can be tracked in the successive images of Figure 21. One can see the effect of the wear algorithm on the height of the elements. The dark shaded elements

The measured and the calculated worn profiles of the seal are compared in Figure 22. The calculated worn profiles are mostly identical with the measured worn surfaces. However, in the vicinity of the lip, the simulation over predicts the wear. It is suspected that the deviation is caused by the assumption of linear wear theory in the wear calculations, while in reality wear is a non-linear function of contact pressure. This simplification nevertheless

were deactivated by the damage check module of the simulation.

appears to be suitable for the regions with less contact pressure.

**Table 1.** Evolution of the calculated worn profile (The deactivated elements are not plotted.)

**Figure 21.** The wear process of one ridge of the seal. The dark elements are deactivated by the damage check algorithm. The images show the left ridge of the seal model after 0, 15, 30, 45 and 60 minutes of sliding, respectively.

**Figure 22.** The measured and calculated worn profile of the reciprocating seal

The developed wear simulation algorithm was able to model the wear process even in the case, when the amount of wear was much larger than the size of the elements used in the FE calculations.

### **4. Conclusion**

Friction of elastomers has a hysteretic part, which is caused by the energy dissipation in the visco-elastic material. The hysteretic friction depends mainly on the sliding velocity and the geometry. It was shown that the effect of sliding velocity on the coefficient of friction can be modelled numerically both in macro and in micro level. Furthermore, the penetration depth of the counter surface and the size of the contact area can be modelled as functions of sliding velocity.

Three different wear simulation techniques were presented for different purposes. The first method can be useful for modelling relatively small amount of wear. The second can model wear even bigger than the elements in the model by global remeshing. The third method deactivates the damaged elements, thus the effect of the surface rupture of the elastomer part can be modelled.

### **Author details**

360 Advanced Elastomers – Technology, Properties and Applications

sliding, respectively.

calculations.

velocity.

**4. Conclusion** 

0

0.1

0.2

0.3

Normal size [mm]

0.4

0.5

0.6

part can be modelled.

**Figure 21.** The wear process of one ridge of the seal. The dark elements are deactivated by the damage check algorithm. The images show the left ridge of the seal model after 0, 15, 30, 45 and 60 minutes of

Initial Measured worn Calculated worn

**Figure 22.** The measured and calculated worn profile of the reciprocating seal

The developed wear simulation algorithm was able to model the wear process even in the case, when the amount of wear was much larger than the size of the elements used in the FE

0 0.5 1 1.5 2 2.5 3 Axial size [mm]

Friction of elastomers has a hysteretic part, which is caused by the energy dissipation in the visco-elastic material. The hysteretic friction depends mainly on the sliding velocity and the geometry. It was shown that the effect of sliding velocity on the coefficient of friction can be modelled numerically both in macro and in micro level. Furthermore, the penetration depth of the counter surface and the size of the contact area can be modelled as functions of sliding

Three different wear simulation techniques were presented for different purposes. The first method can be useful for modelling relatively small amount of wear. The second can model wear even bigger than the elements in the model by global remeshing. The third method deactivates the damaged elements, thus the effect of the surface rupture of the elastomer Nándor Békési *Department of Machine and Product Design, Budapest University of Technology and Economics, Budapest, Hungary* 

### **5. References**

	- [16] Archard J.F. (1953) Contact and rubbing of flat surface, J Appl. Phys., 24: 981–988
	- [17] Schallamach, A. (1971) How does rubber slide? Wear 17: 301-312
	- [18] Fukahori, Y., Yamakazi, H. (1994) Mechanism of rubber abrasion. Part I: Abrasion pattern formation in natural rubber vulcanizate, Wear 171: 195-202
	- [19] Yan, J., and Strenkowski, J. S., (2006), A Finite Element Analysis of Orthogonal Rubber Cutting, Journal of Materials Processing Technology, 174: 102-108
	- [20] Ko D-C., Kim B-M., Choi J-C. (1997) Finite-element simulation of the shear process using the element-kill method, Journal of Materials Processing Technology, 72:1: 129- 140, doi: http://dx.doi.org/10.1016/j.bbr.2011.03.031
	- [21] Eleőd A., Devecz J., Balogh T. (2000) Numerical modelling of the mechanical process of particle detachment by finite element method, Periodica Polytechnica Transportation Engineering 28:1-2: 77-90 (available:http://www.pp.bme.hu/tr/2000\_1/tr2000\_1\_07.html)
	- [22] Békési N., Váradi K., Felhős D. (2011) Wear simulation of a reciprocating seal, Journal of Tribology 133:3: 031601-1-031601-6, doi: http://dx.doi.org/10.1115/1.4004301
