## 1. Introduction

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, 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 steady-state 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 rate-dependent 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.

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 elasto-slide 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: 25ºC.

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.

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.

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

where,

In Eq. (1), *x*_{0}*=X*_{c,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:

where Ω_{1} and Ω_{2} are correspnding frequencies of 5 and 7.5 Hz, and *X*_{1} and *X*_{2} 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:

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:

where *F*_{c} and *F*_{s} 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:

where, *X*_{c} and *X*_{s} 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:

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.

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, *F*_{s}, and damping force, *F*_{d}, 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:

where

From a given start time, *t*_{0}, when *ϕ(t*_{0}*)*=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, *F*_{s}, in one cycle can be obtained as:

Similarly, from a different start point *ϕ(t*_{0}*’)*=2*kπ*, *k*=0,1,…, the damping force, *F*_{d}, is obtained as follows:

The identified nonlinear stiffness and damping forces are shown in Fig. 5. In Figure 5 a 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, *F*_{s}, and the damping force, *F*_{d}, are denoted as solid lines in the force-displacement 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, *F*_{s}, 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:

where, *K* is the linear stiffness, *δK* is the nonlinear stiffness increment, and *λ*_{k} is used to describe the nonlinear stiffness slope. These parameters can be determined by a nonlinear curve-fitting algorithm, in which an objective function, *J*_{s}, is minimized as follows:

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

As shown in Fig. 5b, the shape of the nonlinear damping force, *F*_{d}, 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:

where, *F*_{y} and *λ* represent the yield force and yield parameter respectively. Similarly, the nonlinear damping parameters also can be determined by minimizing the objective function, *J*_{d}, as follows:

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

The reconstructed force is shown as the solid line in the Fig. 5c, and the analytical force-displacement 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.

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

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 force-displacement time histories under sinusoidal loading.

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, *λ*_{k}, decreases as the frequency increases. Similarly, the identified damping parameters are shown in Fig. 7b. The yield force, *F*_{y}, 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 *F*_{y}. 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, *λ*_{k}, 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.

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.

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

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 force-displacement 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 *f*_{i}^{*}, 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.

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, *k*_{0}, 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 *f*_{e} represents the total resisting force of the elasto-slide elements.

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:

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:

(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, *k*_{0}, the initial loading force is given as

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

Thus, while an elastomer is under a sinusoidal displacement loading, a theoretical force-deflection 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 *x*_{0} such that the slide force at a certain yield region is written as

where *p* is a positive odd integer and *v*_{r} is a characteristic reference velocity. Coupled with the lead spring *kφ*(*f*^{*})d*f*^{*}, the internal displacement *x*_{0} at certain yield region is obtained using the function:

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

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

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 fourth-order 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 rate-dependent 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:

where, M, C, and K are mass matrix, equivalent damping matrix, and stiffness matrix, respectively, *F*_{d} 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 rate-dependent elasto-slide damper model, there are theoretically an infinite number of internal variables, *x*_{0}. 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, *f*_{i}^{*}, the corresponding distribution is equal to the area of that element, *φ*( *f*_{i}^{*})Δ*f*^{*}. Each element has an internal variable denoted as *x*_{0i}, *i=1,…,n*, and each *x*_{0i} satisfies the Eq. 22. Thus, the elastomer force *F*_{d} can be described as:

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

(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, *k*_{0}, 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:

From the initial loading curve Eq. 19, yields

Then

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

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:

Notably, *φ*_{0} is a distribution constant and an index of the post yield force level, and *k* and *k*_{0} 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 force-displacement curve when x→0. 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:

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:

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

For the rate-dependent elasto-slide model, the reference velocity *v*_{r} and the exponent *p* need to be determined. The choice of *v*_{r} 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 *φ*_{0} 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.

### 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 elasto-slide 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.

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.

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:

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 force-displacement 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.

## 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.

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.