Modal Analysis of a Discrete Tire Model with a Contact Patch and Rolling Conditions Using the Finite Difference Method

Abstract

Obtaining the modal parameters of a tire with ground contact and rolling conditions represents a challenge due to the complex vibration characteristic behaviors that cause the distortion of the tire’s symmetry and the bifurcation phenomena of the natural frequencies. An in-plane rigid–elastic-coupled tire model was used to examine the 200 DOF finite difference method (FDM) modal analysis accuracy under non-ground contact and non-rotating conditions. The discrete in-plane rigid–elastic-coupled tire model was modified to include the contact patch restriction, centrifugal force and Coriolis effect, covering a range from 0 to 300 Hz. As a result, the influence of the contact patch and the rotating tire conditions on the natural frequencies and modes were obtained through modal analysis.

1. Introduction

To describe in-plane tire mode shapes, wave propagation through the treadband/belt must first be explored. The treadband exposes three types of waves: bending, longitudinal, and rotational. The deformations of these waves are shown in Figure 1. The radial displacement of the belt is affected by the bending wave, whereas the tangential displacement is a consequence of the longitudinal wave; however, the rotational wave appears at a much higher propagation speed of approximately 2500 Hz. Thus, the rotational wave can be ignored if the analysis is in the frequency range of 0–500 Hz [1].

The rubber viscoelastic behavior limits the modal analysis to 500 Hz because the rubber loss factor increases as the frequency increases, so at a high-frequency range, the propagation waves generated from the excitation decay too quickly (see Figure 2). In an out-of-plane tire, the belt/treadband cross-section can vibrate; hence, at each cross-sectional wavelength number, there is an infinite number of circumferential modes, as shown in Figure 3.

Wheeler [3] used a notation for the tire structure modes labeled in two indices, $[c,a]$. The circumferential index, c, represents the number of circumferential bending wavelengths of the belt or the circumferential mode number. The cross-section index, a, represents the number of half-wavelengths in the axial direction of the treadband at a circumferential location where the shape is at extreme radial displacement.

When $n=0$, the breathing mode is the tangential zero mode, which can be categorized as a constant radial deformation of the belt. The breathing mode is a function of the belt membrane’s and radial sidewall’s stiffness. In the first circumferential mode, the rigid mode, the belt moves as a rigid structure [4].

Kindt [6] et al. extended the mode description to apply to the loaded tire modes. As mentioned previously, the loaded tire breaks the symmetry, causing the split of doubled poles that occur on the unloaded tire into two single poles. Two additional labels were introduced, 0 and $extr$. The label 0 is added to the mode with a radial displacement of zero in the middle of the tire/road contact area. The label $extr$ is added for the mode with an extreme radial displacement in the middle of the tire/road contact area. An illustration is presented in Figure 4.

For the non-contact tire, there are two radial symmetric mode shapes rotated by $\frac{180}{2n}$ at each eigenvalue represented by a harmonic function, where n is the mode wavenumber along the circumference. This symmetry of mode shapes will be destroyed for a deformed tire and will generate a unique radial mode shape at each natural frequency (see Figure 5 and Figure 6). Soedel and Prasad [7] extracted the tire mode shapes that are in ground contact from the non-contact tire’s mode shapes and natural frequencies. This approach is based on a small contact patch length (very low normal load) to ignore the shifting of the natural frequencies caused by the loaded tire. In Kozhevnikov’s study, using the wheel with a reinforced tire consisted of a rigid disk attached to sidewalls and tread without slipping. The natural frequency for the loaded tire is a function of the contact patch length. As the contact patch length increases, the tire mass attached to that region will be excluded from the tire vibration mass; thus, the natural frequencies will rise up. Additionally, the two identical mode shapes will be subdivided into two unique mode shapes, symmetric and antisymmetric, under the deformed tire condition [8]. Matsubara et al. performed an experimental modal analysis and applied a cylindrical shell model to realize the effect of the contact patch restriction, which showed that the modal shape and the natural frequency variation are due to contact tire boundary conditions. The natural frequencies are greater in contact tires than in non-contact tires; this variation grows as the radial mode number increases, and the half mode radial number ($n=1.5,2.5,\cdots$) can be used in the cylindrical shell natural frequency as well as the harmonic eigenfunction to predict the modal characteristics of a contact tire [9].

The boundary condition of the one-direction movement of the rim in the vertical direction compared with a fixed rim shows an impact on the first radial natural frequencies and has no effect on tangential natural frequencies. The first radial natural frequencies (rigid mode) are split into two distinct frequencies under the rim movement boundary condition. The split of the rigid mode observed by the elastic ring model is supported by the elastic foundation and 12 DOF doubly curved, geometrically non-linear, thin shell finite element model [10]. Tire mode shapes and natural frequencies depend on the patch and spindle constraints as well as the tire’s construction properties. For the fixed patch, Richards et al. found no distinction in the radial or tangential modes between fixed and pinned wheel spindle constraints, but for a freewheel spindle boundary condition compared with the fixed or pinned constraints, the first radial mode is different, illustrating the effect of the mass wheel movement. The mass wheel movement becomes irrelevant for higher radial modes. In torsional modes, a difference occurs among the pinned, fixed, and free wheel, caused by inertia [11]. Rim movement was found to influence the breathing mode ($n=0$) and the first radial mode ($n=1$).

A simpler rigid ring model, presented by Zegelaar, can be used to obtain the vibration transmission between the contact point and the rim via the transfer function matrix. The transfer function matrix relates the displacement response to the input forces. A vertical residual stiffness at the contact patch is attached to the rigid ring in order to involve the contribution of the higher-order radial modes in the total deformation of the tire. Additionally, horizontal residual stiffness was introduced to obverse the slip contact model (brush model). This method used a rigid ring transfer function matrix to relate the contact force to the rim movement, with an assumption that the contribution of the higher order mode would be residual stiffness [12].

The influence of the tire parameters on the natural frequency was explored using a thin cylindrical shell theory with the Rayleigh method. The contribution ratio of the tire substructure parameters against the natural frequency showed that the first tire modes were affected by the sidewall stiffness only because these modes are the rigid modes. For higher modes, the contribution of the tension term becomes dominant. The torsion and bending tread element does not influence tire vibration [13].

Shokouhfar and El-Gindy generated a 3D finite element model of a radial-ply truck tire in LS-Dyna to predict the natural frequencies and corresponding mode shapes. The nonlinear large deformation finite element theory includes the contribution of the stress from the vertical load to tire stiffness matrix. The inflation pressure has an almost proportional relationship with the natural frequencies, yielding greater impact on higher tire modes. Likewise, the natural frequencies tend to increase with the greater vertical load. Additionally, the rolling conditions show an effect on the dynamic tire stiffness [14].

The Coriolis acceleration and centrifugal force need to be explored in the local frame of reference. For a non-rotating flexible ring model, there exists only one frequency for each mode. However, under rotation conditions, there are two frequencies for every mode. This phenomenon is called “the bifurcation of the natural frequencies”, and is created by the Coriolis effect if observed from the rotating coordinate [15] Huang and Soedel [16] showed that the bifurcation is independent of the elastic foundation stiffness. The Coriolis effect shows a higher impact on the natural frequencies compared with the centrifugal force. The centrifugal force affects the tire stiffness but without any bifurcation of the natural frequencies. It has two contributions on a rotating tire: first, it increases the radial load on the rubber ring so that the pre-tension force on the tire is increased [17]. Second, it created a spin-softening effect, acting as a negative contribution on the generalized stiffness matrix [18].

The FEA model by Lardner [19] used a virtually tested tire on a 2.5 diameter circular drum with a cleat acting as an excitation force to perform the vibration analysis under various load and rolling conditions. The cylindrical shell theory model by Matsubara [20] obtains the vibration characteristics of the contact and rolling conditions by applying the receptance method, as well as the frequency equation and mode shape function.

The rotating elastic beam under tension on an elastic foundation model [21] is used to study the vibration characteristics under rotating conditions and to extract the mode shape and natural frequencies from the eigenfunction and natural frequency equation. Tsotras [4] introduced two discretized models, namely a truss and a beam-based model for a tire ring model. The truss model is insufficient for radial analysis due to the lack of bending stiffness, except at very low frequencies. The beam model performed better than the truss model, but the inflation pressure effect as a pre-tension force cannot be implemented.

In this study, a 200 DOF discrete tire model was generated from an in-plane rigid–elastic-coupled tire model using the FDM. The central difference theory discretized the tire model equation of motion to generate stiffness, mass, and damping matrices. The model comprises the tread as an Euler–Bernoulli beam, the sidewall as lumped masses with two-sectional radial stiffness, and the inflation pressure simulated as the axis force. The coupling between the sidewall lumped masses with the rim and the tread extends the frequency band to 300 Hz; this makes the model more precise for use on roads with relatively short and sharp unevenness.

The main issue with full modal models (FEA, for example) is the inability to integrate the tire modal model with simulation software such as ADAMS and CarSim. Since the majority of vehicle dynamics tuning (including suspension and steering system settings) is achieved through simulations, having an accurate reduced order discrete modal model allows the design of the experiment to tune the vehicle frequencies for the best ride performance.

This study aimed to develop the discrete compact model of a pneumatic tire that is able to be integrated with a complete vehicle simulation that covers the range 0–300 Hz, predicting the change in the tire’s modal parameters caused by tire deformation and the rolling condition. It has the ability to implement the contact patch length restriction, as well as the Coriolis and centrifugal forces. The proposed model improved the in-plane rigid–elastic-coupled tire model by including the ground contact restriction and rolling condition. The extra radial stiffness was introduced and attached to the contact patch elements to resolve the tread-ground contact difficulty. The Coriolis force generates a skew-symmetric damping matrix to be added to the system, which gives a better understanding of the concept of the bifurcation of the natural frequencies of tires. However, the centrifugal force affects the generalized stiffness matrix by changing the pre-tension force and acting as a spin softener. The mass, stiffness, and damping matrices produced from the FDM were applied to build a tire model with 200 DOF under ground contact, rolling, and fixed rim boundary conditions. In order to estimate the natural frequencies and mode shapes, a modal analysis simulation was performed.

2. The Undeformed Tire

Liu et al. [22] presented the in-plane rigid–elastic-coupled tire model, which consists of three main structure treads simulated as an Euler–Bernoulli beam, a sidewall simulated as the two-sectional radial stiffness, and the inertia force and inflation pressure simulated as the axis force of the Euler–Bernoulli beam, as shown in Figure 7. Coupling the belt with the rim using the sidewall stiffness will give a precise transfer function prediction within the 180 Hz range. However, adding the sidewall as a lumped mass with two-sectional sidewall radial stiffness will expand the transfer function’s prediction to 300 Hz.

2.1. The Equation of Motion for the Euler–Bernoulli Beam Coupled with the Sidewall Lumped Mass

The equation of motion of the Euler–Bernoulli beam is subject to crucial assumptions that the cross-section is perpendicular to the bending line, there is no rotation of the cross-section due to shear deformation, and the stress in terms of thickness is ignored. The dynamics of the tread based on force and moment balance (see Figure 8) are expressed in the following equation where P is the axial (tension) load, M is the bending moment, f is the external force per unit length, V is the shear force, and $\rho$A$dx$ is the mass of the differential element:

$$V-(V+\frac{\partial V}{\partial x}dx)+fdx+(P+\frac{\partial P}{\partial x}dx)sin(\Theta +d\Theta )-Psin\Theta =\rho Adx\frac{{\partial }^2{u}_r}{\partial t^2}$$

(1)

For small angles and constant inflation pressure P, the force balance becomes:

$$-\frac{\partial V}{\partial x}+f+P\frac{{\partial }^2{u}_r}{\partial x^2}=\rho A\frac{{\partial }^2{u}_r}{\partial t^2}$$

(2)

Applying the moment balance equation:

$$(M+dM)-M-(V+dV)\left(dx\right)-fdx\left(\frac{dx}{2}\right)=0$$

(3)

$$\frac{\partial V}{\partial x}=\frac{{\partial }^{2}M}{\partial x^2}$$

(4)

From the elementary theory of bending at beams where I is a constant moment of inertia and E is a beam elastic modulus:

$$\frac{{\partial }^{2}M}{\partial x^2}=EI\frac{{\partial }^4{u}_r}{\partial x^4}(x,t)$$

(5)

The Euler–Bernoulli beam equation of motion under tension with no internal damping:

$$\rho A\frac{{\partial }^2{u}_r(x,t)}{\partial t^2}+EI\frac{{\partial }^4{u}_r(x,t)}{\partial x^4}-P\frac{{\partial }^2{u}_r(x,t)}{\partial x^2}=f(x,t)$$

(6)

Let the length of the beam $L=2\pi R$, where R is the inflated tire radius, so $x=R\theta$:

$$\rho A\frac{{\partial }^2{u}_r(\theta ,t)}{\partial t^2}+\frac{EI}{R^4}\frac{{\partial }^4{u}_r(\theta ,t)}{\partial {\theta }^4}-\frac{P}{R^2}\frac{{\partial }^2{u}_r(\theta ,t)}{\partial {\theta }^2}=q(\theta ,t)$$

(7)

where (7) is the equation of motion expressed in polar coordinates. Now, the Euler–Bernoulli beam is attached to the lumped sidewall masses $m_s$ with two different values of sectional sidewall stiffness $k_{r1}$ and $k_{r2}$. As shown in Figure 7, $q(\theta ,t)$ represents the forces from the elastic foundation of the beam, where $u_{r,j}$ is the belt radial displacement and $u_{sr,j}$ is the sidewall lumped mass radial displacement:

$$q(\theta ,t)=k_{r1}(u_{sr,j}-u_{r,j})$$

(8)

The axial pretension force of the flexible tread/belt due to inflation pressure derived from Figure 9 is shown as follows:

$$P=F_N=\frac{1}{2}{\int }_0^{\pi }\left(pb\right)Rsin\theta d\theta =p_{o}Rb$$

(9)

where p is the inflation pressure and $p_o$ the total pressure applied to the tire width b. The in-plane rigid–elastic-coupled tire model equation:

$$\rho A\frac{{\partial }^2{u}_r(\theta ,t)}{\partial t^2}+\frac{EI}{R^4}\frac{{\partial }^4{u}_r(\theta ,t)}{\partial {\theta }^4}-\frac{p_{o}b}{R}\frac{{\partial }^2{u}_r(\theta ,t)}{\partial {\theta }^2}=k_{r1}(u_{sr,j}-u_{r,j})$$

(10)

Applying Newton’s second law to the sidewall lumped masses will produce a second equation coupled with the Euler–Bernoulli governing equation:

$$m_s\frac{{\partial }^2{u}_{sr,j}}{\partial t^2}=k_{r1}(u_{r,j}-u_{sr,j})-k_{r2}{u}_{sr,j}$$

(11)

2.2. Analytical Solution for the Mode Shape and Natural Frequencies

For proportional damping, Equations (10) and (11) were used to obtain the natural frequencies and mode shapes. Let $u_r=R_r\left(\theta \right)e^{i{\omega }_{n}t}$, where $R_r\left(\theta \right)=Acos\left(n(\theta -\varphi )\right)$ is the tread harmonic eigenfunction, and $u_{sr}=R_{sr}\left(\theta \right)e^{i{\omega }_{n}t}$, where $R_{sr}\left(\theta \right)=Bcos\left(n(\theta -\varphi )\right)$ is the sidewall harmonic eigenfunction. ${\omega }_n$ is the natural frequency, $A_{amp}$ and $B_{amp}$ are the mode shape amplitudes, n is the modal number and $\varphi$ is the arbitrary angle, which depends on the excitation orientation. Thus, the equations are re-written in matrix form as follows:

$$\left[\begin{array}{cc}\frac{EI}{R^4}{n}^4+\frac{p_{o}b}{R}{n}^2+k_{r1}-\rho A{\omega }_n^2& -k_{r1}\\ -k_{r1}& -m_s{\omega }_n^2+k_{r1}+k_{r2}\end{array}\right]\times \left[\begin{array}{c}{A}_{amp}\\ B_{amp}\end{array}\right]=\left[\begin{array}{c}0\hfill \\ 0\hfill \end{array}\right]$$

(12)

The characteristic equation is:

$$\begin{array}{c}\hfill A\rho m_s{\omega }_n^4-\left(\frac{EI{n}^4}{R^4}+k_{r1}+\frac{p_{o}b{n}^2}{R}\right)m_S+\rho A\left(k_{r1}+k_{r2}\right){\omega }_n^2+\\ \hfill \left(\frac{EI{n}^4}{R^4}+\frac{p_{o}b{n}^2}{R}+k_{r1}\right)\left(k_{r1}+k_{r2}\right)=0\end{array}$$

(13)

The solution of Equation (13) for every mode number $n\ge 0$ offers two frequency values and two mode shapes (in-phase and out-phase). Substituting the calculated natural frequency functions from (13) in Equation (12) yields the mode shape amplitude for the tread and sidewall stiffness. The positive and negative ratio between these two amplitudes decides which natural frequency functions represent the in-phase or out-phase frequencies:

$$\rho A{m}_S{\omega }^4-B{\omega }^2+C=0$$

(14)

$$C=\left(k_{\mathrm{r}1}+k_{\mathrm{r}2}\right)\left(\frac{EI{n}^4}{R^4}+\frac{p_{o}b{n}^2}{R}+k_{\mathrm{r}1}\right)-m_S\left(\frac{EI{n}^4}{R^4}+\frac{p_{o}b{n}^2}{R}+k_{\mathrm{r}1}\right)$$

(15)

$$B=-\rho A\left(k_{\mathrm{r}1}+k_{\mathrm{r}2}\right)$$

(16)

$${\omega }_{n,in-phase}=\frac{\sqrt{2}\sqrt{\rho A{m}_S\left(B-\sqrt{-4\rho A{m}_{S}C+B^2}\right)}}{2\rho A{m}_S}$$

(17)

$${\omega }_{n,out-phase}=\frac{\sqrt{2}\sqrt{\rho A{m}_S\left(B+\sqrt{-4\rho A{m}_{S}C+B^2}\right)}}{2\rho A{m}_S}$$

(18)

2.3. Approximate Solution Using the Finite Difference Method

In mathematics, FDMs are numerical methods for solving differential equations by approximating them with difference equations. Finite difference equations approximate the derivatives based on the Tylor series expansion. Thus, FDMs are discretization methods. FDMs convert a linear (non-linear) ODE/PDE into a system of linear (non-linear) equations, which can then be solved by matrix algebra techniques. There are three main techniques for FDM: backward, forward, and central difference methods [25]. Using the central FDM:

$$\frac{\partial u_r}{\partial \theta }=\frac{u_{r,j+1}-u_{r,j-1}}{\delta \theta }$$

(19)

$$\frac{{\partial }^2{u}_r}{\partial {\theta }^2}=\frac{u_{r,j+1}+u_{r,j-1}-2u_{r,j}}{\delta {\theta }^2}$$

(20)

$$\frac{{\partial }^4{u}_r}{\partial {\theta }^4}=\frac{u_{r,j-2}-4u_{r,j-1}+6u_{r,j}-4u_{r,j+1}+u_{r,j+2}}{\delta {\theta }^4}$$

(21)

where $\delta \theta =\frac{2\pi R}{N}$ and N is the number of elements around the tire’s circumference. Substituting Equations (19) through (21) in Equation (10) yields the discrete equation of the beam in the form of ordinary differential equations:

$${\alpha }_1{u}_{r,j+2}+{\alpha }_2{u}_{r,j+1}+{\alpha }_3{u}_{r,j}+{\alpha }_2{u}_{r,j-1}+{\alpha }_1{u}_{r,j-2}-u_{sr,j}{k}_{r1}+\rho A\frac{{\partial }^2{u}_{r,j}}{\partial t^2}=0$$

(22)

where ${\alpha }_1=\frac{EI}{\delta {\theta }^4{R}^4}$, ${\alpha }_2=-(\frac{p_{o}b}{R\delta {\theta }^2}+\frac{4EI}{R^4\delta {\theta }^4})$ and ${\alpha }_3=\frac{6EI}{R^4\delta {\theta }^4}+\frac{2p_{o}b}{R\delta {\theta }^2}+k_{r1}$.

Rewriting Equations (22) and (11) in a matrix form yields the following:

$$\left[\begin{array}{cc}{M}_t\hfill & 0\hfill \\ 0\hfill & M_s\hfill \end{array}\right]\left[\begin{array}{c}{\ddot{u}}_{r,j}\\ {\ddot{u}}_{sr,j}\end{array}\right]+\left[\begin{array}{cc}{K}_{ur}\hfill & K_{urs}\hfill \\ K_{urs}\hfill & K_{us}\hfill \end{array}\right]\left[\begin{array}{c}{u}_{r,j}\hfill \\ u_{sr,j}\hfill \end{array}\right]=0,$$

$$M_t=\left[\begin{array}{ccccc}\rho A\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & \rho A\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & \rho A\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & \rho A\hfill & 0\hfill \\ .\hfill & .\hfill & .\hfill & .\hfill & .\hfill \end{array}\right]$$

$$M_s=\left[\begin{array}{ccccc}{m}_s\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & m_s\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & m_s\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & m_s\hfill & 0\hfill \\ .\hfill & .\hfill & .\hfill & .\hfill & .\hfill \end{array}\right],$$

$$K_{urs}=\left[\begin{array}{ccccc}-k_{r1}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & -k_{r1}\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & -k_{r1}\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & -k_{r1}\hfill & 0\hfill \\ .\hfill & .\hfill & .\hfill & .\hfill & .\hfill \end{array}\right]$$

$$K_{ur}=\left[\begin{array}{ccccccc}{\alpha }_3\hfill & {\alpha }_2\hfill & {\alpha }_1\hfill & 0\hfill & 0\hfill & {\alpha }_1\hfill & {\alpha }_2\hfill \\ {\alpha }_2\hfill & {\alpha }_3\hfill & {\alpha }_2\hfill & {\alpha }_1\hfill & 0\hfill & 0\hfill & {\alpha }_1\hfill \\ {\alpha }_1\hfill & {\alpha }_2\hfill & {\alpha }_3\hfill & {\alpha }_2\hfill & {\alpha }_1\hfill & 0\hfill & 0\hfill \\ .\hfill & .\hfill & .\hfill & .\hfill & .\hfill & .\hfill & .\hfill \\ .\hfill & .\hfill & .\hfill & .\hfill & .\hfill & .\hfill & .\hfill \\ {\alpha }_1\hfill & 0\hfill & 0\hfill & {\alpha }_1\hfill & {\alpha }_2\hfill & {\alpha }_3\hfill & {\alpha }_2\hfill \\ {\alpha }_2\hfill & {\alpha }_1\hfill & 0\hfill & 0\hfill & {\alpha }_1\hfill & {\alpha }_2\hfill & {\alpha }_3\hfill \end{array}\right],$$

$$K_{us}=\left[\begin{array}{ccccc}{k}_{r2}+k_{r1}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & k_{r2}+k_{r1}\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & k_{r2}+k_{r1}\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & k_{r2}+k_{r1}\hfill & 0\hfill \\ .\hfill & .\hfill & .\hfill & .\hfill & .\hfill \end{array}\right]$$

Let $M_G=\left[\begin{array}{cc}{M}_t\hfill & 0\hfill \\ 0\hfill & M_s\hfill \end{array}\right]$,$K_G=\left[\begin{array}{cc}{K}_{ur}\hfill & K_{urs}\hfill \\ K_{urs}\hfill & K_{us}\hfill \end{array}\right]$ and $C_G$ is the proportional damping matrix.

Let $\ddot{u}=\left[\begin{array}{c}{\ddot{u}}_{r,j}\hfill \\ {\ddot{u}}_{sr,j}\hfill \end{array}\right]$, $\dot{u}=\left[\begin{array}{c}{\dot{u}}_{r,j}\hfill \\ {\dot{u}}_{sr,j}\hfill \end{array}\right]$ and $u=\left[\begin{array}{c}{u}_{r,j}\hfill \\ u_{sr,j}\hfill \end{array}\right]$.

The general equation of motion in matrix form becomes:

$$\left[M_G\right]\ddot{u}+\left[C_G\right]\dot{u}+\left[K_G\right]u=0$$

(23)

Substituting $q={\left[M_G\right]}^{1/2}u\left(t\right)$ in Equation (23) and multiplying with ${\left[M_G\right]}^{-1/2}$ gives:

$$I\ddot{q}+{\left[M_G\right]}^{-1/2}\left[C_G\right]{\left[M_G\right]}^{-1/2}\dot{q}+{\left[M_G\right]}^{-1/2}\left[K_G\right]{\left[M_G\right]}^{-1/2}q=0$$

(24)

$$I\ddot{q}+{\tilde{C}}_G\dot{q}+{\tilde{K}}_{G}q=0$$

(25)

The proportional damping system can be represented with real modes where ${\tilde{K}}_G={\left[M_G\right]}^{-1/2}\left[K_G\right]{\left[M_G\right]}^{-1/2}$ is the mass normalized stiffness matrix in k/m. ${\tilde{K}}_G$ is a symmetric, positive definite, and generates real and distinct eigenvalues. The eigenvalues of ${\tilde{K}}_G$ are the squares of the system’s natural frequencies. However, in this case, the eigenvectors of ${\tilde{K}}_G$ are related to the mode shapes by ${\varphi }_i=M^{1/2}{p}_i$ where ${\varphi }_i$ are the mode shapes and $p_i$ are the eigenvectors of ${\tilde{K}}_G$ [26].

3. The Deformed Tire

In practice, the contact forces in the contact patch are unknown because they are determined by the treadband displacement. To resolve the treadband contact difficulty, a minor modification to the physical model is required. A small change to tire tread rubber in the contact patch area, and an extra component (radial stiffness, $k_{cw}$) is added to the model in the radial direction around the outside surface of the tread circumference. The height of these added components is assumed to be zero. The extra radial stiffness component that linked the tread with the road can be considered extremely high because the road surface has sufficiently higher stiffness than the tire (see Figure 10).

The discrete equations introduced in (22) and (11) were used to implement ${\alpha }_3=\frac{6EI}{R^4\delta {\theta }^4}+\frac{2p_{o}b}{R\delta {\theta }^2}+k_{r1}$ for non-contact elements and ${{\alpha }_3}^{\prime }=\frac{6EI}{R^4\delta {\theta }^4}+\frac{2p_{o}b}{R\delta {\theta }^2}+k_{r1}+k_{cw}$ for contact patch elements in the stiffness matrix $K_{ur}$. The matrix was resolved for eigenvalues and eigenvectors to obtain the natural frequencies and mode shapes under the deformed tire. It was assumed that $k_{cw}=5.93\times {10}^8$ N/m was large enough to constrain a zero deflection of the tire mode shapes in the contact patch area.

4. The Rotating Tire in Body Frame

In this section, we consider a rotating tire with a constant rotational velocity $\mathrm{\Omega }$. Additional complexity is added to the system due to fictitious forces generated in a non-inertial frame (i.e., a rotating body frame attached to the tire origin). The three fictitious forces in the non-inertial reference frame are the centrifugal force, the Coriolis force, and the Euler force. In tire research, a constant velocity of the tire is always assumed, so the Euler force is not considered.

The equation of motion in the rotating reference frame can be derived by transforming the stationary reference frame equation of motion into the moving coordinate frame (see Figure 11). This transformation relates the acceleration in the rotating reference frame to the acceleration in the stationary reference frame. Let $\left\{u\right\}={\left[\begin{array}{c}{u}_{r,j}{u}_{sr,j}\end{array}\right]}^T$ be the position vector in a fixed frame of reference for a non-rotating system; $\left\{q\right\}={\left[\begin{array}{c}{q}_{r,j}{q}_{sr,j}\end{array}\right]}^T$ is the position vector in a rotating body frame. The transformation from a stationary frame to a rotating frame is given by

$$\left\{\ddot{q}\right\}=\left\{\ddot{u}\right\}-2\left\{\mathrm{\Omega }\right\}\times \left\{\dot{q}\right\}-\left\{\mathrm{\Omega }\right\}\times (\left\{\mathrm{\Omega }\right\}\times \left\{q\right\})-\left\{\dot{\mathrm{\Omega }}\right\}\times \left\{q\right\}$$

(26)

The Euler force is zero for a constant spin speed:

$$F_E=\left\{\dot{\mathrm{\Omega }}\right\}\times \left\{q\right\}=0$$

(27)

The Coriolis force is:

$$F_c=2\left\{\mathrm{\Omega }\right\}\times \left\{\dot{q}\right\}=2\left[G\right]\left\{\dot{q}\right\}$$

(28)

where G = blockdiag ([$\mathrm{\Omega }$]) with size $2n\times 2n$ and $\left[\mathrm{\Omega }\right]=\left[\begin{array}{cc}0& -{\mathrm{\Omega }}_3\\ {\mathrm{\Omega }}_3& 0\end{array}\right]$

For centrifugal force:

$$F_g=\left\{\mathrm{\Omega }\right\}\times (\left\{\mathrm{\Omega }\right\}\times \left\{q\right\})=\left[D\right]\left\{q\right\}$$

(29)

where D = blockdiag $\left(\left[\widehat{\mathrm{\Omega }}\right]\right)$ with size $2n\times 2n$, $\left[\widehat{\mathrm{\Omega }}\right]=\left[\begin{array}{cc}-{\mathrm{\Omega }}_3^2& 0\\ 0& -{\mathrm{\Omega }}_3^2\end{array}\right]$

The discrete equation in the rotating frame:

$$\left[M_G\right]\left\{\ddot{q}\right\}+(\left[C_G\right]+2\left[M_G\right]\left[G\right])\left\{\dot{q}\right\}+(\left[K_G^{\ast }\right]+\left[M_G\right]\left[D\right])\left\{q\right\}=\left\{F\right\}$$

(30)

$$\left[K_G^{\ast }\right]=\left[\begin{array}{cc}{K}_{ur}^{\ast }\hfill & K_{urs}\hfill \\ K_{urs}\hfill & K_{us}\hfill \end{array}\right]$$

(31)

where ${\alpha }_3^{\ast }=\frac{6EI}{R^4\delta {\theta }^4}+\frac{2P^{\ast }}{R^2\delta {\theta }^2}+k_{r1}$ and $P^{\ast }=p_{o}bR+\rho A{R}^2{\mathrm{\Omega }}^2$ is the pre-tension force on the tread due to the inflation pressure and the centrifugal force:

$$K_{ur}^{\ast }=\left[\begin{array}{ccccccc}{\alpha }_3^{\ast }\hfill & {\alpha }_2\hfill & {\alpha }_1\hfill & 0\hfill & 0\hfill & {\alpha }_1\hfill & {\alpha }_2\hfill \\ {\alpha }_2\hfill & {\alpha }_3^{\ast }\hfill & {\alpha }_2\hfill & {\alpha }_1\hfill & 0\hfill & 0\hfill & {\alpha }_1\hfill \\ {\alpha }_1\hfill & {\alpha }_2\hfill & {\alpha }_3^{\ast }\hfill & {\alpha }_2\hfill & {\alpha }_1\hfill & 0\hfill & 0\hfill \\ .\hfill & .\hfill & .\hfill & .\hfill & .\hfill & .\hfill & .\hfill \\ .\hfill & .\hfill & .\hfill & .\hfill & .\hfill & .\hfill & .\hfill \\ {\alpha }_1\hfill & 0\hfill & 0\hfill & {\alpha }_1\hfill & {\alpha }_2\hfill & {\alpha }_3^{\ast }\hfill & {\alpha }_2\hfill \\ {\alpha }_2\hfill & {\alpha }_1\hfill & 0\hfill & 0\hfill & {\alpha }_1\hfill & {\alpha }_2\hfill & {\alpha }_3^{\ast }\hfill \end{array}\right]$$

$$\left[M_G\right]\left\{\ddot{q}\right\}+\left(\left[C_G\right]+2\left[M_G\right]\left[G\right]\right)\left\{\dot{q}\right\}+\left(\left[K_G^{\ast }\right]+\left[M_G\right]\left[D\right]\right)\left\{q\right\}=0$$

(32)

The state space form:

$$\left\{\dot{q}\right\}=\left[\begin{array}{c}{\dot{q}}_1\\ {\dot{q}}_2\end{array}\right]=\left[\begin{array}{cc}0& I\\ -M_G^{-1}\left(\left[K_G^{\ast }\right]+\left[M_G\right]\left[D\right]\right)& -M_G^{-1}\left(2\left[M_G\right]\left[G\right]+\left[C_G\right]\right)\end{array}\right]\left[\begin{array}{c}{q}_1\\ q_2\end{array}\right]=Aq$$

(33)

Equation (33) can be solved to find the natural frequencies and mode shapes, where $\left[M_G\right]$,$\left[C_G\right]$ and $\left[K_G\right]$ are the symmetric matrix orthogonal to the mode shape $\mathrm{\Phi }$:

$$\left[C_G\right]=\alpha \left[M_G\right]+\beta \left[K_G\right]$$

(34)

5. Results

The results shown in this section are based on the parameters presented in

Table 1. Known geometrical and structural parameters of a GL073A tire [27].

Parameters	Symbol	Unit	Value
Tread width	b	$\mathrm{m}$	0.35
Inflation pressure	$p_0$	$\mathrm{N}/{\mathrm{m}}^2$	$7.9\times {10}^5$
Tire radius	R	$\mathrm{m}$	0.65
Density per rad of sidewall	$m_s$	$\mathrm{kg}/\mathrm{rad}$	10
Density per line of tread	$\rho \mathrm{A}$	$\mathrm{kg}/\mathrm{m}$	19.64

and

Table 2. Under-identified geometrical and structural parameters of the GL073A tire [27].

Parameters	Symbol	Value	Unit
Radial stiffness connecting the sidewall and tread	$k_{r1}$	$6.686\times {10}^6$	$\mathrm{N}/\mathrm{m}$
Radial stiffness connecting the sidewall and rim	$k_{r2}$	$4.431\times {10}^6$	$\mathrm{N}/\mathrm{m}$
Bending stiffness of tread	E·I	$25.697$	$\mathrm{N}/\mathrm{m}$
Proportional coefficient of mass matrix	$\alpha$	$8.2\times {10}^{-5}$	$\mathrm{N}/(\mathrm{m}/\mathrm{s})$
Proportional coefficient of stiffness matrix	$\beta$	$1.1\times {10}^{-4}$	$\mathrm{N}/(\mathrm{m}/\mathrm{s})$

.

5.1. Undeformed Tire

The finite difference method’s (FDM) accuracy was tested with the analytical solution for an undeformed non-rotating tire. The FDM model was built with 200 DOF (100 elements for the tread and 100 elements for the sidewall). Figure 12 compares the mode shapes between FDM and analytical harmonic eigenfunctions, representing the in-phase mode 3 and mode 4; however, all in-phase modes are located between 0 and 180 Hz. The outer line expresses the deformation of the tread, whereas the inner line is the sidewall deformation. For lower modes, using a 200 DOF FDM model estimates a good harmonic shape for both the sidewall and tread. On the other hand, Figure 13 shows the comparison between modes for the out-of-phase state, and all out-of-phase modes are placed between 180 and 300 Hz. As shown in Figure 14, the number of elements plays a vital role in the higher natural frequency values. For 40 tread elements, the natural frequencies below mode 5 for the in-phase and out-phase show a good approximation with the exact solution; however, above mode 5, the number of elements must be increased to approximately 200 segments to provide a good approximation for higher modes. Furthermore, the mode shapes and the natural frequencies were experimentally obtained by Liu et al. [22] and agreed with the FDM and analytical results obtained here.

5.2. Deformed Tire

In Section 3, the mode shapes and natural frequencies for the deformed tire are illustrated. The addition of an extra element (radial stiffness, $k_{cw}$) in the contact patch length modified some diagonal elements of the global stiffness matrix, as discussed in Section 2. The contact patch length specifies the number of diagonal elements that must be adjusted. As shown in Figure 15 and Figure 16, the two identical mode shapes for each eigenvalue can be subdivided into two unique mode shapes: symmetric and anti-symmetric. This phenomenon can be observed in both in-phase and out-of-phase modes. Figure 17 identifies the increase in natural frequencies due to the deformed condition, calling attention to the fact that the increase in natural frequencies depends on the contact patch length; thus, the vertical load. The symmetric mode shows a change of approximately $2\%$ compared with the undeformed condition and approximately 10–$12\%$ in the anti-symmetric mode. These results are validated by the study of Lopes et al. [28]. Another factor that controls the number of elements, in addition to the modes of interest, is the contact patch length. The precise contact patch length requires more elements to accurately predict the shifting natural frequencies.

5.3. Rotating Tire

Figure 18 shows the dispersion relation obtained by solving the system under rotation. Note that the x axis is the radial wave number, $k_{\varphi }$, related to the radial mode number, n, by $k_{\varphi }=n/R$. The natural frequencies in the y axis are observed from the local body frame when the rotation speed value = $20\pi$ rad/s. The frequency of the deformed tire is higher than that of the undeformed condition at the same wavenumber, $k_{\varphi }$, which means that the tire’s deformation increases its stiffness. The phase velocity of the flexural wave moves faster in stiffer tires; this aspect is observed by comparing the slope of the dispersion curve between the deformed and undeformed tires. However, adding the rotation effect to the system increases the pre-tension force due to the centrifugal force; thus, the tire stiffness increases. Although there is a spin-softening effect, the positive centrifugal impact is higher on the pre-tension. The rotation causes the dispersion curve to be asymmetric; this is called the “bifurcation” effect. This means that if the rotating tire is excited at any given point, two opposite wave propagations throughout the tread will be created with various phase speeds (positive wave and negative wave). The change in the wave propagation speed breaks its physical symmetry, although the rotating tire’s geometric symmetry is still intact; as a result, standing waves cannot be produced from a single frequency under rotation conditions, as shown in Figure 19.

6. Conclusions

• A compact, discrete, in-plane rigid–elastic-coupled tire model was developed and modified to include the contact patch length restriction and the tire rolling conditions.
• The results compare favorably with those found in the literature for undeformed, deformed, and rotating tires.
• The proposed model is able to be integrated with vehicle models and covers a 0–300 Hz frequency range without ignoring the change in the tire’s modal parameters caused by tire deformation and rolling conditions. It has the ability to implement the contact patch length restriction as well as the Coriolis and centrifugal forces.
• The proposed tire model has several advantages: (1) it can give a prediction of the changes in natural frequencies under rolling and ground contact conditions; and (2) it is easy to change the tire’s boundary and operating conditions and to allow the user to define the number of discrete elements to be used.
• Predicting the change in the tire’s modal parameters under various conditions gives a better estimation of the transfer function between the road and vehicle when attached to a complete vehicle model. In future studies, the model can be used to build a modal model for the tire/vehicle to investigate the complete dynamic behavior of the vehicle under various road and operating conditions.