Empirical Models for the Viscoelastic Complex Modulus with an Application to Rubber Friction

: Up-to-date predictive rubber friction models require viscoelastic modulus information; thus, the accurate representation of storage and loss modulus components is fundamental. This study presents two separate empirical formulations for the complex moduli of viscoelastic materials such as rubber. The majority of complex modulus models found in the literature are based on tabulated dynamic testing data. A wide range of experimentally obtained rubber moduli are used in this study, such as SBR (styrene-butadiene rubber), reinforced SBR with ﬁller particles and typical passenger car tyre rubber. The proposed formulations offer signiﬁcantly faster computation times compared to tabulated/interpolated data and an accurate reconstruction of the viscoelastic frequency response. They also link the model coefﬁcients with critical sections of the data, such as the gradient of the slope in the storage modulus, or the peak values in loss tangent and loss modulus. One of the models is based on piecewise polynomial ﬁtting and offers versatility by increasing the number of polynomial functions used to achieve better ﬁtting, but with additional pre-processing time. The other model uses a pair of logistic-bell functions and provides a robust ﬁtting capability and the fastest identiﬁcation, as it requires a reduced number of parameters. Both models offer good correlations with measured data, and their computational efﬁciency was demonstrated via implementation in Persson’s friction model.


Introduction
Viscoelastic materials exhibit both viscous and elastic characteristics when stressed so that their behaviour is in-between that of a purely viscous liquid and a perfectly elastic solid. When a viscoelastic material is deformed, part of the deformation energy dissipates, and the rest is stored as reversible elastic energy. The properties of viscoelastic materials such as rubber can vary widely according to temperature, frequency and chemical composition [1]. Rubber is a highly deformable material employed in a wide range of products, from everyday necessities, e.g., shoe soles, to industrial applications-notably, vehicle tyres. Such applications demand large deformations, vibration damping [2] and enhanced gripping characteristics (e.g., tyres and conveyor belts). With safety-critical components such as tyres, it is important to be able to predict friction for different types of rubber; however, the non-linear nature of rubber poses significant modelling challenges [3].
The response of a rubber material can be obtained from rheological models. Linear viscoelastic materials can be represented by combinations of springs and dashpots. The most common ones are the Maxwell and Voigt-Kelvin models [4] which use springs and dashpots connected in series or in parallel, respectively (see Figure 1). The standard linear solid (SLS) is an upgraded rheological model in which an extra spring is added in parallel with the Maxwell unit to allow for an accurate representation of a larger group of polymers [5]. Additionally, the Weichert model offers a relaxation spread over a longer period by using a single spring in parallel with several Maxwell units [6,7]. Xu et al. [8] proposed a non-linear rheological model to offer better accuracy in simulating the material's relaxation modulus and viscoelastic responses. A logistic-type function approximation was proposed by [8] to simulate the viscoelastic response of spring-dashpot networks. Another recent study employed a combination of a generalised Maxwell model and a relative fraction derivative model to reproduce viscoelastic material behaviour [9]. Both models proposed herein describe rubber's viscoelastic response in the frequency domain. As shown in [10], rubber friction strongly depends on the rubber's frequency response, and this is accounted for in Persson's friction model [11]. The first model uses a group of polynomials to capture the viscoelastic frequency response with great accuracy, whilst the other method uses a generic empirical formulation with few parameters to expedite the parameter identification process and to offer a compact set of equations. Thereafter, the proposed formulations are applied to Persson's friction model without thermal effects [11] to evaluate their computational efficiency and accuracy. The paper is structured as follows: In Section 2, the measured moduli used as the benchmark are presented. In Section 3, the polynomial equation method and the empirical model are described. In Section 4, the ability of both models to fit measured data and their performances as part of Persson's friction model are evaluated. Section 5 provides a summary of the conclusions.

Measured Viscoelastic Modulus
When a sinusoidal force is applied to a purely elastic material, the stress and deformation occur concurrently, so that both quantities are in phase. On the contrary, when the same sinusoidal force is applied to a purely viscous fluid, the deformation will lag the stress by one-quarter cycle (π/2 radians). This means that when the deformation is maximal, the force is minimal and vice versa. In the case of a viscoelastic material, the phase lag between the force and the deformation lies somewhere between zero and a quarter cycle. Mathematically, a sinusoidal strain can be written using Euler's formula, as shown in Equation (1), where ε 0 is the amplitude of the strain and ω is the frequency: Similarly, the stress can be written as: where σ 0 is the amplitude of the stress and δ is the phase angle. In Figure 2 the stress and strain phasors are shown as they rotate at a frequency ω, with the strain lagging behind the stress by an angle δ. The following relationships can be deduced from Figure 2: With the above definitions, the dynamic or complex modulus will have a real and an imaginary part. The real or storage modulus is defined as the ratio between the real part of the stress and the strain: The imaginary or loss modulus is defined as the ratio between the imaginary part of the stress and the strain: Equations (1)-(6) lead to the broadly used complex or dynamic modulus formula [12]: Finally, the tangent modulus, also known as the loss tangent is defined as: By definition, the modulus of a material is considered as the overall resistance of the material to an applied deformation. Due to its viscoelastic nature, the rubber modulus is split into elastic (storage), E , and viscous (loss), E", components, denoting the ability of the material to store and dissipate energy as heat, respectively. The real and imaginary components of a viscoelastic modulus are collectively referred to as the material's complex modulus. Equally important is the loss tangent or tan δ, which is defined as the ratio of energy loss to energy storage. The area near the peak value denotes high energy loss and is considered as the area where the tyre is designed to be operating most of the time. While spring-dashpot-based analytical models are useful for representing the physics of viscoelastic materials, the experimentally measured complex modulus as a function of excitation frequency is often used instead. In the case of rubber, the complex modulus is shown to be strongly linked to friction [13,14]. Figure 3 shows the typical complex modulus of a rubber used in tyre manufacturing. Such moduli of rubber materials are obtained using oscillatory strain in tension or shear form, while monitoring the resulting stress [15]. When rubber is dynamically stretched and released, the returned energy is less than the energy put into the rubber in the first place. This viscoelastic effect can only be tested dynamically. Dynamic mechanical analysis (DMA) instruments measure the viscoelastic modulus in response to applied oscillating strain. The stress response to the deformation is recorded as a function of time or frequency. The rubber sample is oscillated at different frequencies and the same process is repeated at various temperatures. Results are then shifted according to the time-temperature superposition principle [16], forming a master curve and covering a wide range of frequencies at a chosen reference temperature. The literature was scrutinised and DMA data for different viscoelastic materials were gathered [4,11,13,15,[17][18][19][20][21][22][23][24][25][26][27]. While spring-dashpot networks can be used to reproduce DMA data, the required structure of the network is sometimes difficult to decide a priori, and the demand on parameter identification can be significant. In the following sections it is shown that the DMA responses of different materials can be approximated successfully using our proposed formulations with a limited set of parameters.

Proposed Empirical Models
Two empirical approaches for fitting the viscoelastic complex modulus using piecewise polynomial fitting, or a pair of logistic-bell functions, are presented below.

Approximation Using Polynomials
The first model proposes a piece-wise fitting of both the storage and loss moduli using polynomials. As a first step, the two moduli curves (storage and loss) are divided piecewise into sections that will be described by different polynomials. Typically, both the storage and loss moduli curves include regions that can be described by linear functions, as shown in Figure 4a. These linear functions are estimated first using least squares fitting of all data points between the start and finish of each straight line. Linear regions are typically connected via higher-order polynomials-see Figure 4b. The requirement for any non-linear section is continuity up to-and including-the first derivative with the preceding and succeeding curves. In this context, greatest flexibility is offered by fourthorder polynomials whereby continuity at the start and finish points offers four equations for four parameters. A fifth equation is obtained by requesting that the function passes from a selected midpoint between the start and finish points. A pair of typical first and fourth-order equations used in the process is shown in Equation (9), below. The superscript i indicates the number of the piece-wise curve, with counting starting from the left of the frequency range. This method offers robustness against global optimisation methods and very easy tuning.

The Logistic-Bell Empirical Model
The logistic-bell empirical model (LBEM) is an alternative rheological model that has few parameters that can be easily linked to typifying quantities of the DMA data, but is still fast to compute and can describe the responses of real viscoelastic materials over a wide frequency range. From data gathered in the literature, it can be observed that in the logarithmic scale the storage modulus reassembles the shape of a logistic function and the loss tangent approximates a bell-shaped function. Therefore, a variation of the logistic function is proposed for modelling the storage modulus, as follows: where: • E is the real or storage modulus in (Pa) based on a logistic function. • ry o is a vertical shift. • L is the curve's maximum value. • k the logistic rate of increase or steepness of the curve. • ω is the excitation frequency in rad/s. • rω 0 is the logarithm in base 10 of the frequency value at the sigmoid midpoint.
The loss tangent is modelled using a piecewise bell-shaped function: tan δ(ω) = sy 1 + (p − sy 1 )sech where: • tan δ is the loss tangent. • tω 0 is the logarithm in base 10 of the frequency at the bell's midpoint or peak. • sy 1 is the bell function vertical shift at low frequencies before the peak. • sy 2 is the bell function vertical shift at high frequencies after the peak. • p is the bell function peak value. • ω is the excitation frequency in rad/s. • r 1 is the bell's rate of increase at low frequencies before the peak. • r 2 is the bell's rate of increase at high frequencies after the peak.
Finally, the loss modulus can be derived from Equation (8).

Fitting Accuracy
Several published moduli from the literature were used to assess the quality of fitting of the polynomial method. The storage and loss moduli, and the loss tangent, were fitted separately, and the relevant coefficients are provided in Table 1. For all three materials, three linear parts were combined with two fourth-order polynomials, to obtain the overall curves. In Table 1, any subscript indicates the role of the parameter in Equation (9), and a superscript indicates the sequential number of the respective linear or fourth-order polynomial. Clearly, any two of the three complex modulus curves are sufficient for fully describing any tested material. The correlations between measured and fitted curves for three different rubbers using the polynomial method are shown in Figures 5 and 6. Figure 5 shows that the piecewise polynomial method is capable of following the curvature of the curves faithfully. A quantitative assessment of the quality of fit is provided via the relative percentage error in Figure 6, where it is shown that the maximum relative error for the storage modulus is 0.7%, and for the loss modulus the maximum relative error is 1.49%.  The same dataset was used for assessing the quality of the fit using the LBEM approach. The MATLAB function "lsqnonlin" with the default trust-region-reflective fitting algorithm [28] was used to identify the optimal parameter values of the proposed formulation by minimising the non-linear cost function that effectively consists of the sum of the squared differences between measured and estimated loss tangent and storage modulus. The resulting parameter values for three example datasets have been tabulated in Table 2.  Figure 7 shows a comparison between the measurements and the fitted LBEM models. The results are very satisfactory, as the formula is able to describe all the measured characteristics and is continuous at all times. From Figure 8 it can be appreciated that the relative error for the LBEM is higher than that of the polynomial formulation, as expected. This is because the LBEM approach uses a generic formulation with fewer parameters that can, however, be easily linked to the typifying quantities of the data. For the storage modulus the maximum relative error increased to 1.4%. The loss tangent shows a higher error of up to 27%, but such discrepancies happen at a limited number of discrete points.  The fitting quality of the polynomial and LBEM models was also evaluated by means of the normalised root mean squared error (NRMSE) and the coefficient of determination, denoted R 2 for both the storage modulus and the loss tangent. The results are summarised in Tables 3 and 4.  Table 4. Coefficient of determination for all the compounds using the polynomial and LBEM models.

Results
Results show that the empirical formulations presented herein were able to model all the viscoelastic moduli that were examined with correlation levels similar to those presented in the previous section. More specifically, the two proposed approaches were compared against each other and against other methods from the literature, such as the zero/pole fitting presented in [9] and the fractal derivative generalised Maxwell Model [30]. Specific performance attributes of the models that are addressed include: (a) versatility, (b) number of parameters required, (c) computational efficiency and (d) integration with Persson's friction model [31]. Persson's model was selected as the "carrier" friction model, as it has been heavily relied upon in the literature for the prediction of rubber friction [32][33][34][35][36][37][38].
It was also used because it is a particularly computationally intensive model that will benefit from efficient implementations of the viscoelastic material model. Finally, the effects on model parameters of five SBR (styrene-butadiene rubber) compounds from the study presented in [20] are discussed and the results are linked to the CB (carbon black) content. In line with the findings in [20], the varying levels of CB in the rubber do not affect the peak values of the viscoelastic modulus, but do influence the gradient of the linear sections.

Versatility
Versatility of a model is interpreted as the ability to successfully fit the moduli of different materials. The performances of the models in this respect were demonstrated in the previous section, particularly in Figures 5-8 and Tables 3 and 4. Due to its piecewise nature, the polynomial method is inherently more flexible in terms of fitting any of the three commonly occurring elements of the complex modulus (storage modulus, loss modulus and loss tangent). It is also more likely to fit more accurately a wider range of materials, albeit with increased pre-processing time and an increase in the number of fitting parameters. The LBEM method, on the other hand, requires less pre-processing, does not depend on manual selection of the relevant data segments and relies on fewer parameters than the polynomial method, as described below.

Number of Parameters
Both formulas make use of parameters that can be linked to the critical sections of the data, such as storage modulus transition regions or the peak values of the loss modulus and loss tangent. The LBEM method requires 10 parameters to be defined, whereas the polynomial method needs a minimum of 16 parameters for the storage and loss modulus to be evaluated.
An important quality of the LBEM approach is that the parameters are directly linked to typifying quantities of the data, such as the loss tangent frequency peak location, the curve's growth rate and the loss tangent peak value.

Computational Efficiency
In terms of solution efficiency, both formulas provide a greatly reduced computation time compared with the common data interpolation method, the proposed zero/pole formula of [39] and the FDGM (fractal derivative generalised Maxwell) model [9]. Table 5 compares solution times of Persson's isothermal friction model for a single friction point using a resolution of 200 frequency segments, and several complex modulus implementations. The commonly used data interpolation method was used as a benchmark and had a mean computation time of 0.318 s per friction point. The FDGM model (that utilises a zero/pole formulation) had a computational time of 1.371 s, which is an increase of 332% with respect to the interpolation method. The proposed LBEM model had a reduced computation time of 0.184 s, a reduction of 42% with respect to the interpolation method and 87% with respect to the zero/pole model. On the other hand, the polynomial method resulted in a 93% or 98% more efficient solution compared to interpolation and the zero/pole techniques, respectively, and is also 12% less computationally demanding than the LBEM method. The reduced efficiency of the LBEM approach in comparison to the polynomial method is related to the use of computationally expensive mathematical functions, such as exponentials and the hyperbolic secant.  Figure 9 shows the friction coefficient estimated over a wide range of sliding speeds for the three different materials used previously. Persson's isothermal friction model has been used for this estimation [11] with the model parameters provided in Table 6. Table 6. Commonly used parameters for Persson's isothermal friction model [11], employed in the present study.  All models gave similar friction results, with the LBEM model providing the best agreement with the interpolated data. A discrepancy was produced by the FDGM model, which is most likely a result of extrapolation at frequency extremes. However, the FDGM model may include an increased number of parameters or "zones"; therefore, there is room for improving its fitting accuracy [9].

PHR Content
We used the data published in [20] to study the effect of different concentrations of CB (carbon black) on the DMA data, and by extension, the effects on the identified parameters from both formulations. Five SBR compounds with different parts phr (per hundred rubber) of carbon black, ranging from 0 to 60 phr, are shown in Figure 10. Both models were able to fit the curves with a maximum normalised root mean squared error of 7.75% for the LBEM formulation and 2.12% for the polynomial formulation. From the raw data it can be seen that as the phr of CB increases, so does the storage modulus in the rubbery (low frequency) region. The storage modulus also increases with the phr of CB in the glassy (high frequency) region, albeit at a reduced rate. The peak of the loss tangent is broadly inversely proportional to the concentration of CB. With regard to the polynomial model, Figure 11 demonstrates the linear relationship between the vertical shift of the linear section (see Figure 4a) and the phr concentration of CB in rubber specimens. R-squared values are included in Table 7.  11. Linear regression between the coefficients of the linear polynomial equations (y = ax + b) and the CB concentration. Linear equations reproduce the linear sections of storage and loss modulus, respectively (see Figure 4). Coefficients (a 1 , b 1 ), (a 2 , b 2 ), (a 3 , b 3 ) fit the first, second and third linear sections of the two modulus curves, from left to right.
From Figure 12 it can be observed that the parameters ry 0 , k, sy 2 , r 1 , r 2 and p of the LBEM model have R-squared values equal to or higher than 90%, indicating a high level of correlation with the carbon black concentration. As expected from the raw data analysis, the parameter ry 0 that controls the vertical shift of the storage modulus shows a high correlation. However, the parameter k that controls the logistic rate of increase or steepness of the curve also shows a high correlation, and this was not easy to spot from a visual analysis of the data. In addition, as indicated from the visual analysis of raw data, the parameter p that controls the bell function peak value shows a clear correlation with CB concentration. The loss tangent parameter sy 2 that controls the bell function's vertical shift at high frequencies after the peak, together with r 1 and r 2 that determine the bell's rate of increase at low frequencies and high frequencies, respectively, also have R-squared values above 90%. On the basis of the observed dependency of critical model parameters on CB content, the LBEM model and to a lesser extent the polynomial model, lend themselves to the creation of virtual hypothetical rubber materials to be simulated within friction models, so that the performance of a new tyre with new materials can be predicted.

Conclusions
The proposed empirical models are capable of not only accurately reproducing DMA data but also significantly reducing the computational time. This makes the models ideal for integration into advanced friction models, such as [11,13,14], where the viscoelastic complex modulus is evaluated multiple times and over a wide range of frequencies.
With the polynomial method, good predictions of the storage and loss moduli and the loss tangent can be obtained by selecting the piecewise linear/non-linear segments of the curves. Formulation is versatile in that in can fit a broad range of viscoelastic moduli because of the adjustability of the number of segments. The proposed LBEM formula uses parameters that can be easily linked to the typifying quantities of the data, such us the storage moduli in the rubbery and glassy regions or the frequency where the loss tangent is maximum. The fitting quality of both models was quantitatively evaluated against the data using the NRMSE and R 2 metrics. The R 2 for the polynomial method was above 0.97 for all cases, whereas the performance of the LBEM formula varied more widely, with R 2 varying from 0.99 to 0.75. Comparable results can be obtained with the FDGM model at the expense of computation time. The proposed models are primarily suited to calculations in the frequency domain, but they are not applicable the creep and relaxation tests because the relaxation time cannot be readily deduced from their formulation. With regard to the polynomial model, only a subset of model parameters can be linked to the CB content of the rubber, and by extension, to the chemical composition of the rubber. On the other hand, the LBEM model comes with a smaller parameter set that is linked to material composition; therefore, it can be used for the creation of hypothetical rubber materials to assess/optimise the grip of new tyres ahead of production.