Predicting the behavior of commercially available building foundation isolation materials using a four-parameter fractional Zener model. Colin Bradley 1 Faiz Musafere 2 Josh Havin 3

Matthew Golden 4

Pliteq, Inc. 131 Royal Group Crescent Vaughan, ON, Canada, L4H 1X9

ABSTRACT In this paper we apply a four-parameter fractional Zener model to describe the short- and long-term creep behavior of a viscoelastic recycled-rebonded-crumb-rubber isolating material in noise and vibration mitigation applications. Building from the classical rheological models composed of springs and dashpots suitably combined (Maxwell, Kelvin-Voight, Zener), we include a fractional derivative component in the form of a Scott-Blair element composed of a Caputo derivative of non- integer order. The parameters for each model are identified empirically through measured stress- strain results and hysteresis curves for the material in compression. The models are then compared to the results of a 168h constant load creep test of a commercially available material.

1. INTRODUCTION

As a means of vibration isolation, viscoelastic materials are an attractive choice due to their wide range of mechanical properties and ease of implementation. A detailed understanding of the dynamic behavior of the material is necessary for these parameter-sensitive applications. Viscoelastic materials, by definition, have material properties that are dependent on the applied frequency, static load, dynamic load, and dynamic amplitude, thereby requiring a complex modulus parameter. Here we explore the fundamental behavior of a commercially available recycled-rebonded- crumb-rubber isolation material using Dynamic Mechanical Analysis (DMA) at low to mid frequency ranges in the context of foundation isolation and rail infrastructure. The classic Maxwell, Kelvin- Voigt, and Zener rheological models are considered first and then extended to include a fractional derivative parameter. Using the DMA results, measured based on the procedure and principles of the

1 cbradley@pliteq.com 2 fmusafere@pliteq.com 3 jhavin@pliteq.com 4 mgolden@pliteq.com

standard for measurement of viscoelastic properties, DIN 53513 [1], each model’s parameters are determined and compared to actual dynamic performance and long-term creep. 2. RHEOLOGICAL MODELS FOR VISCOELASTIC MATERIALS

A basic assumption of viscoelasticity is the dependence of stress on strain history in addition to the instantaneous strain of the material. A relaxation function, 𝑒(𝑡) , convoluted over the strain history, 𝜖(𝑡) , of the material captures this assumption and lends it viscoelastic properties, as displayed in Equation 1. 𝜎(𝑡) = ∫𝜖(𝑡−𝜏)𝑒(𝜏)𝑑𝜏

𝑡 0 (1) A Laplace transform of the convolution leads to a complex modulus dependent on frequency, Equation 2. 𝜎(𝑠) = 𝐸 ∗ (𝑠)𝜖(𝑠) = (𝐸 ′ (𝑠) + 𝑖𝐸 ′′ (𝑠))𝜖(𝑠) (2)

Where the complex modulus 𝐸 ∗ is composed of a storage (elastic) modulus 𝐸′ dominated by elastic forces and a loss modulus 𝐸′′ where viscous forces dominate. Using this methodology, viscoelastic problems can be solved as elastic problems. This is known as the elastic/viscoelastic equivalence principle.[2] The component moduli act out of phase with one another, denoted using the imaginary unit ( 𝑖 2 = −1 ). The loss modulus is not to be confused with the loss factor, which is defined as the ratio of the loss modulus to the storage modulus, as it is related to energy dissipation within the material. The value of the dynamic modulus 𝐸 𝑑𝑦𝑛 (𝜎, 𝜔) is simply the magnitude of the complex modulus.

In this study, the material we are testing is a complex composite of recycled-rebonded-crumb- rubber with both microscopic and macroscopic structure, each contributing to the viscoelastic behavior of the bulk material. Chemically, rubber is composed of relatively regular long-carbon chains of monomers such as isoprene (natural rubber), or styrene-butadiene, chloroprene etc. based on petroleum production processes (synthetic rubber). These long-chain polymers are then vulcanized with the addition of sulfur and heat allowing relatively strong connections to establish between chains of polymers, increasing the overall strength, toughness, environmental resistances, and elasticity. [3] The vibration isolation materials tested here are composed of thoroughly sorted rubber particles obtained from grinding tires. The particle size distribution is crucial when recombining the rubber into a bulk material, factors such as the degree of packing, density, particle-particle contact area all affect the overall properties of the rebonded-rubber and can be tuned to the desired values with careful engineering. [4] The macroscopic structure of the vibration isolation product can be adjusted in different ways to allow for an additional parameter depending on the required performance. By introducing geometry to control the contact surface area, for example, the load experienced by the material can be adjusted locally to optimize the materials performance. While advantageous to have access to these parameters when engineering a material, it is a complex system to model in its entirety. Simplifying assumptions regarding the behavior of the system are necessary when considering the viscoelastic behavior of the material. We start by considering the classical Maxwell, Kelvin-Voigt, and Linear Zener models in Figure 1.

Maxwell Kelvin-Voigt Linear Zener

𝜂 1 𝐸 1 𝜎̇ = 𝜂 1 𝜖̇ (3) 𝜎= 𝐸 1 𝜖+ 𝜂 1 𝜖̇ (4) 𝜎+

𝜂 1 ( 𝐸 1 +𝐸 2 )

𝜂 1 𝐸 1 𝜎̇ = 𝐸 2 𝜖+

𝜎+

𝐸 1 𝜖̇ (5)

Figure 1: Schematic and governing equations for several rheological models. Through dynamic mechanical analysis according to DIN 53513, “Determination of Viscoelastic Properties of Elastomers”, we measure the complex modulus and loss factor directly through a procedure involving a frequency sweep from 1 – 120 Hz. [1] In this case, GenieMat  GBV750 is subject to a pre-loading of 150 kPa and a 20% oscillation amplitude, as plotted in Figure 2. The data shows that the Storage Modulus (E’) is just slightly less than the total magnitude of the complex modulus (E*).

Figure 2: Measured dynamic properties of the material. The measured results of the Storage Modulus (E’) have been plotted versus the fitted parameter values for 𝐸 1/2 and 𝜂 1 from each of the three rheological models, Maxwell, Kelvin-Voigt and Zener, as shown in Figure 3. All three models are very close to the measured values with the Zener model showing a slight overprediction. The total Storage Modulus (E’) for the Zener model is calculated from the modeled E 1 and E 2 values, shown as dashed lines in Figure 3.

Figure 3: Comparison of linear rheological model storage moduli vs measured values.

The modeled viscosity ( 𝜂 1 ) used from each of the three rheological models, Maxwell, Kelvin- Voigt, and linear Zener, are shown in Figure 4. The viscosity is an important component of the imaginary portion of the complex modulus. Depending on the model, the remainder of the loss modulus is a function of other parameters. The Zener viscosity lies between Maxwell and Kelvin- Voigt viscosities, which is unsurprising given that it is a composition of these models.

Figure 4: Comparison of linear rheological model viscosities.

3. FRACTIONAL FOUR PARAMETER NON-LINEAR ZENER MODEL

In this section, we begin again by considering the rheological scheme of a fractional Zener and its solution, as expressed in Figure 5. It consists of a spring (with stiffness 𝐸 1 ) in parallel with a spring (with stiffness 𝐸 2 ) in series with a fractional element, colloquially called a spring pot with fractional derivative operator D sp of order α. The fractional derivative operator allows us to take advantage of increased flexibility to describe both retardation and relaxation without the need to maintain a strictly exponential shape. We define the Caputo fractional derivative of order 𝛼 in Equation 6. The

fractional component of the model is implemented as a Scott-Blair element. The procedure for fitting the parameters for a viscoelastic material can be reviewed in references [5,6,7].

𝑑 𝛼

𝑓 ′ ( 𝜏 )

(𝑡−𝜏) 𝛼 𝑑𝜏 𝑥 0 (6)

𝑓(0)

𝑑

𝛼 =

𝐷 𝑠𝑝

𝑑𝑡 𝛼 𝑓(𝑡) =

𝑑𝑡 ∫

Γ(1−𝛼)𝑡 𝛼

Fractional Non-Linear Zener

−1

1

1

𝐸 ∗ = (

(7)

𝐸 1 (𝐷 𝑠𝑝 (𝑖𝜔)) 𝛼 )

𝐸 2 +

Figure 5: Schematic drawing and equation for the complex modulus of a fractional non-linear Zener rheological model.

As the name suggests, the key component is the fractional element which can model the accumulation of energy of an elastic element, and the dissipation of energy of a viscous element. The resulting balance of these two of behaviors is controlled by the alpha parameter. Once again, the four parameters are calibrated by curve-fitting the solution of the Non-Linear Fractional Zener Model expressed in Figure 5 to the complex stiffness of GenieMat GBV750 measured with the constraints described in Section 2. In Figure 6, a comparison of this model to the measurements shows similar variance to that of the linear models.

Figure 6: The measured complex modulus of GenieMat GBV750 (discrete points) compared with a calibrated Non-Linear Fraction Zener Model (continuous or dashed lines)

In conjunction with Figure 3 these values offer a hint as to why this model provides little additional value in comparison to the linear rheological schemes; the complex modulus is primarily controlled by the storage modulus. Although the non-linear model offers a clearer view of the material’s viscosity, this parameter is not dominant to a high enough degree to significantly alter the material’s dynamic stiffness. 4. IMPLICATIONS FOR CREEP BEHAVIOR

As the static and long-term implications are considered, each classical model provides some insight as to where the theoretical models fall relative to the tested materials. Here we will derive the behavior of each model when subjected to a step-change stress and allowed to experience creep. The shapes of these curves are then calibrated and compared to the measured behavior of a high-capacity building isolation material, GenieMat GBV3000.

The theoretical creep response for each linear rheological model is derived in a similar manner. We begin with the linear Zener model by taking its governing equation (Equation 5) to the Laplace domain with the initial condition that the time-derivatives at 𝑡= 0 for both stress and strain are zero and arrive at Equation 8.

𝜂 1 ( 𝐸 1 +𝐸 2 )

𝜂 1 𝐸 1 s σ̅ = 𝐸 2 ϵ̅ +

σ̅ +

𝐸 1 𝑠 𝜖̅ (8)

Where 𝜎̅ and 𝜖 ̅ denote the Laplace transform of stress and strain respectively. Solving for the strain in the Laplace domain and realizing that the Laplace transform of a step-change in stress corresponds to 𝜎(𝑡) = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡→ 𝜎̅ =

𝜎 0

𝑠 , leads us to Equation 9.

𝜎0

𝑠 (1+ 𝜂1𝑠

𝐸2 )

𝜖 (𝑠) ̅̅̅̅̅̅ =

𝐸1 𝑠 (9)

𝐸 2 + 𝜂1(𝐸1+𝐸2)

Reverting to the time domain yields the step-change response outlined in Equation 12, displayed alongside the solutions for the Maxwell and Kelvin-Voigt models.

Maxwell Kelvin-Voigt Linear Zener

1

𝑡

𝜎 0 𝐸 1 (1 −𝑒𝑥𝑝[−𝑡

𝐸 1 𝜂 1 ])

𝜎 0 𝐸 2 (1 −

𝐸 2 𝐸 1 +𝐸 2 𝑒𝑥𝑝[−𝑡

𝐸 1 𝐸 2 (𝐸 1 +𝐸 2 )𝜂 1 ])

𝜖(𝑡) = 𝜎 0 (

𝜖(𝑡) =

𝜖(𝑡) =

𝐸 1 +

𝜂 1 )

(10)

(11)

(12)

Figure 8: Schematic solution for the response to a step-change in stress.

By using measured data to calibrate the parameters within each of these creep functions we can compare the expected response for the theoretical case to the measured strain response. In this section we continue to use a rebonded-recycled rubber, however with a higher bearing capacity. In this case, we are continuously monitoring the deflection of building isolation material, GenieMat GBV3000, as it is subjected to 4.5 MPa of stress for 168 hours. It is important to note that the stress

applied to the GenieMat GBV3000 was increased from ambient to target load over the course of many seconds, as opposed to an instantaneous step change. The measured data is plotted in Figure 9 in comparison to the calibrated models, and the ramp-up period is omitted for the sake of clarity.

Figure 9: Ideal rheological models compared to real viscoelastic creep.

Figure 9 makes it immediately evident that the Maxwell model’s representation of creep has very little correlation with the actual long-term behavior of this material. The Maxwell rheological model is known to be a viscoelastic liquid model, where an applied stress causes the strain to increase linearly with time without bounds. [8] A closer look at the derived step-change response equations for Kelvin-Voigt and linear Zener models (Equations 11,12) shows that as time goes to infinity they approach the deformation of a pure elastic material. The difference between purely elastic behavior and the viscoelastic deflection decays exponentially as time increases. With more parameters to tune, the Zener model can be fit slightly closer to any dataset but due to the same fundamental behavior, the difference is small in this case. The real creep measured is very near 1%/decade while both models underestimate this, likely due to the non-ideal variances already noted above. Due to their parallel structure, linear Zener and Kelvin-Voigt plots offer similar levels of accuracy in long-term creep as time increases to infinity. In this case, the differences between the parameters of the two models have little impact on the overall plot due to stiffness being the dominant parameter of the material, and the quasi-static state of strain. 5. CONCLUSIONS

The response of viscoelastic materials to dynamic stress is complex and requires constant comparison with measured values. The Maxwell, Kelvin-Voight, Zener and Fractional Zener rheological models can all be used to model the Storage Modulus (E’) and Loss Modulus (E’’). For materials with relatively low viscosity, such as the ones tested, all four linear models can approximate the complex dynamic modulus with relative accuracy up to at least 120 Hz. When creep is considered, there are clear differences among the models. The Maxwell model represents a material with unlimited deflection and is best suited for modeling viscoelastic liquids. The Kelvin-Voigt and Zener models have fundamentally similar shapes and can both be used to model long term creep for materials with properties like the ones tested in this paper. Neither model accurately modeled the short-term creep behavior, under 1000 sec.

At this time, none of these models can be used to model long-term creep based on the various static material properties collected so far. Further work will explore this possibility. 6. REFERENCES

1. DIN 53513 - Determination of Viscoelastic Properties of Elastomers.pdf . (n.d.). 2. Etienne, B., & Jean-Michel, L. (2017). Viscoelastic Vibration Toolbox For Use with MATLAB .

1–131. 3. Bhowmick, A. K., & Stephens, H. L. (2001). Handbook of elastomers . New York: M. Dekker.

p.733 4. Bradley, C., Musafere F., Havin, J., (2022). Dynamic mechanical analysis of rubber-based

products in under ballast mat rail applications NOISECON 2022 5. W. Grzesikiewicz, A. Wakulicz, A. Zbiciak, Non-linear problems of fractional calculus in

modeling of mechanical systems, International Journal of Mechanical Sciences 70 (2013) 90– 98, https://doi.org/10.1016/j. ijmecsci.2013.02.007 . 6. Zbiciak, C. Kraskiewicz, A. Al Sabouni-Zawadzka, J. Pełczyn ́ ski, S. Dudziak, A novel

approach to the analysis of under sleeper pads (USP) applied in the ballasted track structures, Materials 13 (11). doi:10.3390/ma13112438.. 7. Feldmann, N., Bause, F., & Henning, B. (2020). B1.4 - Determining fractional Zener model

parameters from low frequency DMA measurements . 172–176. https://doi.org/10.5162/sensor2017/b1.4 8. S. J. Ostrand, J. (2017). Some Deviations From Linear Dynamics Due To More Accurate

Damping Models . Lund University Press www.byggmek.lth.se