# An adapted two-step approach to simulate nonlinear vibrations of patterned tires rolling on a smooth surface

Zakaria Knar, Jean-Jacques Sinou, Sébastien Besset, Vivien Clauzon

► To cite this version:

Zakaria Knar, Jean-Jacques Sinou, Sébastien Besset, Vivien Clauzon. An adapted two-step approach to simulate nonlinear vibrations of patterned tires rolling on a smooth surface. Applied Mathematical Modelling, 2024, 135, pp.477-503. 10.1016/j.apm.2024.06.032. hal-04649187

HAL Id: hal-04649187

https://hal.science/hal-04649187v1

Submitted on 16 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# An adapted two-step approach to simulate nonlinear vibrations of patterned tires rolling on a smooth surface

Zakaria Knar<sup>a,b</sup>, Jean-Jacques Sinou<sup>a,*</sup>, Sébastien Besset<sup>a</sup>, Vivien Clauzon<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Rolling noise

Nonlinear dynamics

Contact mechanics

Multiscale modeling

## ABSTRACT

The tire/road noise is one of the major problems facing the tire industry due to the nuisance felt by the vehicle's driver and passengers. Moreover, it is expected in the coming years that this sound nuisance will be one of the main sources of vehicle noise due to the transition from combustion engine driven vehicles to electric vehicles. Tire manufacturers have therefore refined the design of their tire structures to find technological solutions to reduce traffic noise. Tire/road noise is also generated by various mechanisms which depend on different parameters such as the properties of the tire, road texture and driving conditions. This noise is partly caused by the acoustic radiation induced by the tire vibrations due to contacts with the road. The simulation and analysis of the tire's vibrations remains a challenge for the tire industry. Indeed, the simulation of the full dynamic response of a rolling patterned tire requires not only taking onto account various nonlinearities but also the multi-scale nature of the generated dynamic response. Contrary to a straightforward strategy that consists in using a time integrator to predict the multi-scale dynamic response, the strategy proposed in this study is based on a two-step approach to separate the dynamics occurring at different scales. The mathematical formulation of the proposed method is detailed as are the different modeling choices to simulate tire rolling under imposed load. The sensitivity of the tire's vibrations response is analyzed under different rolling conditions and with respect to certain key tire tread pattern parameters.

## 1. Introduction

Urban noise is a factor of stress whose effects on public health demand attention. Road traffic is one of the main sources of this nuisance. Progress made by automobile constructors has led to a considerable reduction of engine noise and the advent of the electric vehicle is bringing about a genuine breakthrough. In parallel, the improved aerodynamics of vehicles and the trend towards lower speed limits in outer urban areas contribute towards reducing aerodynamic noise. Increasing attention is being given to the rolling noise produced by tires whose relative contribution has increased to the point of becoming predominant [1]. Therefore, the standards imposed by public authorities and the specifications of automobile constructors focused on the noise emitted by tires are more and more severe [2].

<footnote><sup>a</sup> École Centrale de Lyon, LTDS UMR CNRS 5513, 36 avenue Guy de Collongue, 69134 Ecully Cedex, France</footnote>

<footnote><sup>b</sup> Michelin Campus RDI, Rue Bleue, 63118 Cébazat, France</footnote>

<footnote>* Corresponding author at: École Centrale de Lyon, LTDS UMR CNRS 5513, 36 avenue Guy de Collongue, 69134 Ecully Cedex, France.
E-mail address: jean-jacques.sinou@ec-lyon.fr (J.-J. Sinou).</footnote>
Mechanical impact is the main structural mechanism involved in the generation of rolling noise. Indeed, the vibrations generated by the impact between the tire tread and the road are propagated in the tire whose surface radiates acoustic waves [3]. On the other hand, air pumping is considered to be the main aerodynamic mechanism generating rolling noise. This mechanism results from fluctuations of air entering and exiting the contact patch. Indeed, the air captured at the entry of the contact patch is compressed in the cavities of the tire tread pattern and the texture of the road. These small volumes of air are then suddenly released at the exit of the contact patch [4,5]. The vibrations generated by these two mechanisms are magnified by resonance phenomena due to the geometric effects of rolling (pipe resonance, horn effect) [6–10]. The rolling conditions, the tires’ properties and the road’s texture have a significant impact on rolling noise [6,11].

From the industrial standpoint, tire design is complex since its performances (grip, wear, rolling resistance, noise, etc.) are often antagonistic and require making choices upstream to achieve the best compromise possible with respect to the targeted performances. This phase of predicting tire's performance relies mainly on numerical modeling and simulation. Indeed, simulating the dynamics of tire rolling permits, for example, predicting the grip of the tire tread and obtaining the distribution of stresses in the tire, information that cannot be obtained through experimental measurements. In addition, numerical simulation permits reducing the number of certification tests intended to validate the reliability and performances of tire structures and thus testing a considerable number of designs to finally choose a tire that satisfies the performances desired as well as possible.

However, modeling a tire system and the associated simulation of the dynamic response present numerous challenges. For example, taking onto account large deformations requires processing the geometric nonlinearity. The simulation of the dynamic response can no longer be done with the modal superposition method as in the framework of linear elastodynamics. Moreover, the mixture of rubbers present in the tire is modelized in practice by a nonlinear visco-hyperelastic constitutive law [12,13]. Moreover, the stiffness of the metal and textile reinforcements embedded in the layers of the tire architecture is far higher than that of rubber. This can lead to additional complexity when modeling such systems. Furthermore, the excitation force imposed on the tire results from its contact with the road surface. Contact and friction phenomena are characterized by nonlinear laws and therefore require specific numerical treatment [14,15]. Moreover, the two-scale nature of the tire's geometry is an additional source of difficulty during the numerical treatment of this type of problem. Indeed, the characteristic dimension of the tire corresponds to its diameter (generally about 50 centimeters) and controls the macroscopic rolling whose frequency coincides with the rotation frequency. Furthermore, the vibratory response of tires is excited partially by the shape of the tread pattern whose size is in the region of several millimeters. This vibratory response occurs at medium and high frequencies with lower amplitudes.

Therefore, these different complexities involved in modeling and resolving the mechanical problem to be treated require a con-
trolled and efficient numerical strategy in order to simulate and predict the vibratory behavior of a tire, and of the rolling noise
radiated by the structural vibrations of the tire. A classical approach for dealing with this type of problem consists in solving the
dynamic problem with a time integration strategy using a conservative and stable numerical scheme, capable of providing the entire
dynamic response of the tire system for a predefined excitation [16,17]. However, the considerable cost of computation inherent
to this type of numerical strategy means that it is ill-adapted at present for industrial applications. The second approach consists in
developing better adapted numerical techniques based on understanding the physics of the problem and defining a configuration of
the state of the tire through modeling adapted to the excitation representing the physical phenomena of noise generation. This is the
case for example of the arbitrarily Lagrangian-Eulerian formulation (ALE) [18,19] used to simulate the internal noise representative
of passenger comfort and generated mainly by the roughness of the road [20–22]. However, this approach is only adapted for ax-
isymmetric tire systems. Thus, it cannot be used in the simulation of external noise given the importance of the tire tread pattern
geometry. Also, one of the originalities of the proposed study is to make a scientific contribution to developing an adapted numerical
strategy to model the tire’s vibratory response induced by the impact of the pattern of the tire tread. One of the benefits of the
proposed numerical approach is to allow an efficient and robust estimation of the multi-scale dynamic response by decoupling the
dynamics occurring at different scales. This step is an essential component for estimating and predicting at a later stage rolling noise
which can be calculated by solving the problem of the associated acoustic radiation whose boundary conditions correspond to the
vibrations velocity [22,23]. It should be noted that the calculation and analysis of the acoustic problem will not be dealt with in this
study, which focuses in particular on the problem of the vibratory behavior of a tire system. More specifically, the main challenges
and objectives of the study are to demonstrate the feasibility of an efficient numerical strategy that not only provides a reliable
estimation of the evolution of the vibratory behavior of the tire but also allows determining the impact of the different characteristic
parameters of a tire structure, given that from the industrial standpoint rolling noise is often studied at a more advanced stage when
many parameters of the full mechanical system are already determined and the optimization of this performance is generally carried
out by modifying only the geometry of the tire tread pattern.

This paper is structured as follows. Section 2 is dedicated to the mathematical development of the modeling strategy employed to simulate the vibratory response of tires. Then, the influence of the rolling conditions and design parameters of the tread pattern of a tire on the vibratory response are analyzed using the numerical method developed previously in Section 3. The last part focuses on the extension of the model of a tire system by considering the case of rolling under an imposed load. The mathematical formulation and the numerical problem are discussed and the tire’s vibratory response is analyzed and compared to the previous results in order to determine the potential contributions linked to this type of model of a tire rolling under an imposed load.

## 2. Modeling

The objective of this section is to define the mechanical problem under consideration and to outline the proposed methodology for solving the associated vibrations problem. Initially, the mechanical quantities required to use the theory of large deformations
<figure><img src="image_4.png" /><figcaption>Fig. 1. Deformation of solid Ω in unilateral contact with the rigid foundation Γ₀.</figcaption></figure>

are introduced. Then, the mathematical formulation of the unilateral frictionless contact is presented after which the equations of motion are given with the corresponding weak formulation. Lastly, the numerical strategy used to simulate the vibratory response of the tire system is detailed.

## 2.1. Notations

This section is dedicated to the mathematical description of the deformation of a solid in a time interval [0, T] in which its reference configuration is denoted Ω (see Fig. 1). The adherence Ω coincides with the domain occupied by the structure at time t = 0. A Lagrangian description is used to treat the geometrical nonlinearity. The kinematic quantities are described with respect to the material points of the reference configuration. The position of a material point X ∈ Ω is an independent variable and the resolution of the equations of motion consists in finding the position φ(X, t) = x of the material point in the deformed configuration or the displacement field u(X, t) = x − X.

The gradient of a quantity in the reference configuration is denoted ∇. To describe the deformation of Ω, several mechanical quantities are introduced. The second order identity tensor is denoted I and the deformation gradient tensor is represented by F = I + ∇u. The Jacobian of φ is denoted J = det(F). C = F†F represents the Cauchy-Green tensor and E = 1/2(C − I) is the Green-Lagrange strain tensor. σ is the Cauchy stress tensor, P = JσF† is the first Piola-Kirchhoff stress tensor and S = JF⁻¹σF⁻¹ is the second Piola-Kirchhoff stress tensor. To treat rubber-like materials, a general hyperelastic constitutive law is considered. This law is derived from the strain energy function W that depends on the deformation E [24]. The second Piola-Kirchhoff stress tensor is given by:

$$S = \frac{\partial W}{\partial E}(E)$$

with the corresponding fourth-order elasticity tensor:

$$C = \frac{\partial S}{\partial E} = \frac{\partial^2 W}{\partial E \partial E}$$

This tensor is necessary in the resolution of the nonlinear variational problem resulting from the local equilibrium equations with Newton's method and will also be used in setting up the strategy to resolve the dynamic problem.

## 2.2. Contact mechanics

A unilateral frictionless contact with a rigid foundation Γ₀, whose unit outward normal is n_c, is considered as shown in Fig. 1. The significant quantitative difference between the stiffness of the road and the tire's rubber justifies the consideration of a rigid foundation. Γ_c is a restricted part of the boundary ∂Ω and represents the surface where contact phenomena may occur.

To define the condition of non-interpenetration between the deformable solid and the surface Γ₀, the Euclidian projection of point x = φ(X, t) on the surface Γ₀ is defined as follows:

$$\tilde{Y}(X,t) = \arg \min_{Y \in \Gamma_o} \|\varphi(X,t) - Y\| \quad (1)$$

The surface Γ₀ is assumed to be sufficiently regular to ensure the unicity of the projection Y. It is worthwhile recalling that other methods exist for measuring the distance between the points of the potential contact zone and the rigid foundation such as the ray-tracing method. The latter is more suitable than the projection method in the case of non-regular surfaces [25]. The gap normal function measuring the normal distance between the position x and its projection Y is defined by:

$$g_n(X, t) = n_c \cdot (\varphi(X, t) - \tilde{Y}(X, t)) \quad (2)$$
The condition of non-interpenetration is formulated mathematically as follows:

$$g_n(\mathbf{X}, t) \geq 0 \qquad (3)$$

The contact mechanics phenomena impose a relation between the kinematic quantities and the stresses on the contact zone. Therefore, the stress vector **T** is expressed using the first Piola-Kirchhoff tensor since the equations of motion are written in the reference configuration in the framework of large deformations:

$$\mathbf{T}(\mathbf{X}, t) = \mathbf{P}(\mathbf{X}, t)\mathbf{N}(\mathbf{X}) \qquad (4)$$

The full formulation of the contact conditions requires the decomposition of the stress vector in the basis defined by the unit outward normal **n** in the contact zone:

$$\mathbf{T}(\mathbf{X}, t) = \Pi_n \mathbf{n} + \Pi_t \qquad (5)$$

It should be noted that the tangential pressure $\Pi_t$ is null since the contact considered is frictionless. To prevent the penetration of the deformable solid in the rigid foundation $\Gamma_o$, a non-positive contact pressure arises in the contact zone so that

$$\Pi_n \le 0 \qquad (6)$$

If a point is in contact with the rigid foundation, then $g_n = 0$. Otherwise, if the point is not in contact, then $\Pi_n = 0$. This allows defining the following non-penetration and non-adhesion condition:

$$g_n \Pi_n = 0 \qquad (7)$$

The set of conditions defined by Equations (3), (6) and (7) correspond to the normal contact constraints. These conditions correspond to Hertz-Signorini-Morea (HSM) conditions or those of Karush-Kuhn-Tucker (KKT) [26].

$$g_n \ge 0, \quad \Pi_n \le 0, \quad g_n \Pi_n = 0 \quad \text{in} \quad \Gamma_c \qquad (8)$$

The Signorini conditions not only define a nonlinear contact law but also a multivalued equation for $g_n = 0$. The form of the conditions is identical to the KKT first-order optimality conditions present in the constrained optimization problems [27].

## 2.3. Equations of motion

The equations of motion of a deformable solid $\Omega$ in unilateral contact with a rigid foundation $\Gamma_o$ are obtained by adding the Signorini conditions described previously to the conservation equations of the linear and the angular momentums. The boundary of the material domain $\partial\Omega$ is composed of $\Gamma_u$ where the displacement is imposed, $\Gamma_\sigma$ where the external forces are imposed, and $\Gamma_p$ where a follower pressure ($\sigma \cdot \mathbf{n} = \mathbf{p}\mathbf{n}$) is imposed. The latter follows the direction of the normal of the actual configuration (contrary to a classical Neumann condition) and will be used in what follows to model the tire inflation pressure. The associated mathematical problem can be written in the following form:

$$\left\{ \begin{array}{lll} \nabla \cdot \mathbf{P} + \mathbf{g}_0 = \rho \ddot{\mathbf{u}} & \text{in} & \Omega \times [0, T] \\ \mathbf{u} = \mathbf{u}_d & \text{on} & \Gamma_u \times [0, T] \\ \mathbf{P}\mathbf{N} = \mathbf{f}_0 & \text{on} & \Gamma_\sigma \times [0, T] \\ \mathbf{P}\mathbf{N} = \rho J (\mathbf{F}^{-1}\mathbf{N}) & \text{on} & \Gamma_p \times [0, T] \\ g_n \ge 0, \quad \Pi_n \le 0, \quad g_n \Pi_n = 0 & \text{on} & \Gamma_c \times [0, T] \end{array} \right. \qquad (9)$$

where $\rho$ is the density and $g_0$ represents the external volume forces vector. The weak formulation of the problem is obtained by multiplying the equation of the local equilibrium by a test function $\mathbf{w}$ and by performing an integration by parts as follows:

$$\int_{\Omega} \nabla \cdot \mathbf{P} \cdot \mathbf{w} d\Omega = \int_{\partial\Omega} \mathbf{P}\mathbf{N} \cdot \mathbf{w} d\Gamma - \int_{\Omega} \mathbf{P} : \nabla \mathbf{w} d\Omega = \int_{\Gamma_{\sigma}} \mathbf{f}_0 \cdot \mathbf{w} d\Gamma + \int_{\Gamma_{p}} pJ(\mathbf{F}^{-1}\mathbf{N}) \cdot \mathbf{w} d\Gamma + \int_{\Gamma_{c}} \mathbf{P}\mathbf{N} \cdot \mathbf{w} d\Gamma - \int_{\Omega} \mathbf{P} : \nabla \mathbf{w} d\Omega \qquad (10)$$

The virtual work of the contact forces (on the potential contact zone $\Gamma_c$) can be written in the following form:

$$\int_{\Gamma_c} \mathbf{P}\mathbf{N} \cdot \mathbf{w} d\Gamma = \int_{\Gamma_c} (\Pi_n \omega_n + \mathbf{\Pi}_t \cdot \mathbf{w}_t) d\Gamma = \int_{\Gamma_c} \Pi_n \omega_n d\Gamma \qquad (11)$$

The gradient of the test function appears after the development of the integration by parts. The functional spaces, to which the solution $\mathbf{u}$ and the test functions $\mathbf{w}$ belong, are introduced so that the integrals are well-defined.

$$\mathcal{V} = \{ \mathbf{u} \in H^1(\Omega) | \mathbf{u} = \mathbf{u}_d \text{ on } \Gamma_u \} \\ \mathcal{V} = \{ \mathbf{u} \in H^1(\Omega) | \mathbf{u} = \mathbf{0} \text{ on } \Gamma_u \} \qquad (12)$$
where $H^1(\Omega)$ is the first-order Sobolev space of the functions whose first weak derivative belongs to the space of the finite energy functions ($L^2(\Omega)$). $\mathcal{U}$ designates the space of the kinematically admissible displacements and depends on time since the imposed displacement $\mathbf{u}_d$ can evolve in time. The vector space $\mathcal{V}$ is the tangent space of $\mathcal{U}$. The weak formulation of the mechanical problem is given by:

$$ egin{cases} 
\text{Find } \mathbf{u} \in \mathcal{U} \\ 
\begin{aligned} 
\int_{\Omega} \rho \ddot{\mathbf{u}} \cdot \mathbf{w} d\Omega + \int_{\Omega} \mathbf{P} : \nabla \mathbf{w} d\Omega &= \int_{\Omega} \mathbf{g}_0 \cdot \mathbf{w} d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f}_0 \cdot \mathbf{w} d\Gamma + \int_{\Gamma_{\rho}} p J(\mathbf{F}^{-\top} \mathbf{N}) \cdot \mathbf{w} d\Gamma + \int_{\Gamma_c} \Pi_n w_n d\Gamma \quad \forall \mathbf{w} \in \mathcal{V} \\ 
g_n \ge 0, \quad \Pi_n \le 0, \quad g_n \Pi_n = 0 \quad \text{on} \quad \Gamma_c \times [0, T] 
\end{aligned} 
\end{cases} 
\quad (13) $$

The penalization method is chosen to treat the normal contact conditions. Accordingly, the virtual work of the contact forces is approximated as follows [28]:

$$ \int_{\Gamma_c} \Pi_n w_n \, d\Gamma = \int_{\Gamma_c} \varepsilon_n g_n^{-}(\mathbf{u}) w_n \, d\Gamma \quad (14) $$

with $\varepsilon_n$ being the penalization parameter and $g_n^{-} = \min(g_n, 0)$ the penetration function. The definition of the penetration function permits penalizing only the points belonging to the active zone of the contact and therefore avoid using the active set strategy [15]. Furthermore, the solution obtained with a penalization $\varepsilon_n = \infty$ is equal to the exact solution obtained with Lagrange multipliers [29]. However, large values of $\varepsilon_n$ will lead to an ill-conditioned tangent matrix. The analytical contribution of each term of the weak formulation to the tangent matrix and to the nonlinear residual can be found in [13].

## 2.4. Vibratory response

The main idea of the methodology proposed is to linearize the dynamic problem around the quasi-static configuration. The mathematical formulation of the method is described in this section. The hypotheses applied are explained and justified. All the mathematical quantities presented are indexed by 1 (and 2, respectively) for the quasi-static (and dynamic one respectively) problem. A unilateral frictionless contact is considered between a deformable solid, whose reference configuration is noted $\Omega$, and a rigid foundation whose outward unit normal is $\mathbf{n}_c$. The notations introduced in the previous parts are used in the developments that follow.

The weak formulation of the quasi-static problem is given by:

$$ \left\{ \begin{array}{ll} \text{Find } \mathbf{u}_1 \in \mathcal{U} & \\ \displaystyle \int_{\Omega} \mathbf{P}_1 : \nabla \mathbf{w} \, d\Omega = \mathbf{R}_e + \mathbf{R}_1^c & \forall \mathbf{w} \in \mathcal{V} \end{array} \right. \quad (15) $$

$\mathbf{R}_e$ represents the variational term of the external forces and $\mathbf{R}_1^c$ is the variational term of the contact forces. The weak formulation of the dynamic problem is given by:

$$ \left\{ \begin{array}{ll} \text{Find } \mathbf{u}_2 \in \mathcal{U} & \\ \displaystyle \int_{\Omega} \rho \ddot{\mathbf{u}}_2 \cdot \mathbf{w} \, d\Omega + \int_{\Omega} \mathbf{P}_2 : \nabla \mathbf{w} \, d\Omega = \mathbf{R}_e + \mathbf{R}_2^c & \forall \mathbf{w} \in \mathcal{V} \end{array} \right. \quad (16) $$

$\mathbf{R}_e$ represents the variational term of the external forces which is the same as that of the quasi-static problem. $\mathbf{R}_2^c$ is the variational term of the contact forces. The mathematical expression of $\mathbf{R}_1^c$ and $\mathbf{R}_2^c$ depends on the type of contact and the numerical method used to treat the contact constraints (penalization, augmented Lagrangian...) [14].

The main hypothesis of the methodology proposed in this study is that the dynamic response $\mathbf{u}_2$ results from the quasi-static response $\mathbf{u}_1$ and from the vibratory response $\mathbf{u}_v$. The terms of the higher order of $\mathbf{u}_v$ are neglected ($\|\mathbf{u}_v\| \ll 1$):

$$ \mathbf{u}_2(\mathbf{X}, t) = \mathbf{u}_1(\mathbf{X}, t) + \mathbf{u}_v(\mathbf{X}, t) + o(\mathbf{u}_v) \quad (17) $$

The inertia effect is taken into account only in the calculation of the vibrations and a quasi-static process is assumed to be sufficient to describe the macroscopic rolling of the structure. Indeed, the characteristic time of the inertia effects is assumed to be lower than that of the excitation of the vibratory response and higher than that of the macroscopic rolling (i.e. $t_{\text{vibration}} < t_{\text{inertia}} = \sqrt{\frac{m_e}{k_e}} < t_{\text{rolling}}$ with $m_e$, and $k_e$ respectively, a quantity characterizing the global mass and the global stiffness of the structure as the norm of the mass matrix and stiffness matrix, respectively, for example).

The Dirichlet boundary conditions are the same for both problems. Thus, a null vibratory displacement is imposed on $\Gamma_u$:

$$ \mathbf{u}_v = \mathbf{0} \quad \text{on} \quad \Gamma_u \quad (18) $$

The active contact zone and the variational terms of the contact forces are assumed to be invariant between the two problems (i.e. $\Gamma_{c,1}^t \approx \Gamma_{c,2}^t \approx \Gamma_c^t$ and $\mathbf{R}_2^c \approx \mathbf{R}_1^c$). It is noteworthy that the active contact zone is not exactly the same between the two problems
since the inertia changes the stiffness of the contact and thus modifies its shape [30]. However, the change induced by the variation of the stiffness is assumed to be negligible.

The normal component of displacement in the active contact zone is determined by the coordinates of the rigid foundation for the two problems. Thus, the normal component of the vibrations response, on the active contact zone $\Gamma_c^t$ given by the solution of the quasi-static problem, is null:

$$ \mathbf{u}_v \cdot \mathbf{n}_c = 0 \quad \text{on } \Gamma_c^t \qquad (19) $$

$\mathbf{n}_c$ is the outward unit normal of the rigid foundation. It should be noted that this boundary condition is an interpretation of the contact conditions as a Dirichlet boundary condition. This is problematic from a mathematical point of view given that the active contact zone $\Gamma_c^t$ is not attached to a fixed material domain. Indeed, the contact zone is spatially included in the potential contact zone's deformed configuration $\varphi(\Gamma_c, t)$. However, its numerical implementation does not present a particular problem following a finite element discretization.

The mechanical quantities of the dynamic problem can be expressed as a function of the quasi-static mechanical quantities while neglecting the higher order terms. The linearization of the deformation tensors is done as follows:

$$ \mathbf{F}_2 = \mathbf{F}_1 + \nabla \mathbf{u}_v \qquad (20) $$

$$ \mathbf{E}_2 = \mathbf{E}_1 + \frac{1}{2}(\mathbf{F}_1^\mathrm{T}\nabla\mathbf{u}_v + \nabla^\mathrm{T}\mathbf{u}_v\mathbf{F}_1 + \underbrace{\nabla^\mathrm{T}\mathbf{u}_v\nabla\mathbf{u}_v}_{\approx 0}) = \mathbf{E}_1 + \frac{1}{2}(\nabla\mathbf{u}_v + \nabla^\mathrm{T}\mathbf{u}_v + \nabla^\mathrm{T}\mathbf{u}_1\nabla\mathbf{u}_v + \nabla^\mathrm{T}\mathbf{u}_v\nabla\mathbf{u}_1) = \mathbf{E}_1 + \mathbf{E}_{\mathbf{u}_v, \mathbf{u}_1} \qquad (21) $$

The constitutive law is linearized as follows:

$$ \mathbf{S}_2 = \frac{\partial W}{\partial \mathbf{E}}(\mathbf{E}_2) = \frac{\partial W}{\partial \mathbf{E}}(\mathbf{E}_1 + \mathbf{E}_{\mathbf{u}_v, \mathbf{u}_1}) = \frac{\partial W}{\partial \mathbf{E}}(\mathbf{E}_1) + \frac{\partial^2 W}{\partial \mathbf{E} \partial \mathbf{E}}(\mathbf{E}_1) : \mathbf{E}_{\mathbf{u}_v, \mathbf{u}_1} = \mathbf{S}_1 + \mathbf{C}_1 : \mathbf{E}_{\mathbf{u}_v, \mathbf{u}_1} \qquad (22) $$

With the expressions obtained in Equations (20), (21) and (22), the linearization of the variational deformation term is performed as follows:

$$ \mathbf{P}_2 : \nabla \mathbf{w} = \mathbf{S}_2 : (\mathbf{F}_2^\mathrm{T} \nabla \mathbf{w}) \qquad (23) $$

which can be expressed as

$$ \mathbf{P}_2 : \nabla \mathbf{w} = \mathbf{P}_1 : \nabla \mathbf{w} + \mathbf{S}_1 : \nabla^\mathrm{T} \mathbf{u}_v \nabla \mathbf{w} + (C : \mathbf{E}_{\mathbf{u}_v, \mathbf{u}_1}) : (\mathbf{F}_1^\mathrm{T} \nabla \mathbf{w}) \qquad (24) $$

By replacing Equation (17) in the dynamic variational problem (16) and using the expression obtained with the linearization of the variational deformation term, the vibrations response is obtained by solving the following variational problem:

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{u}_v \in \mathcal{V} \\ \displaystyle \int_{\Omega} \rho \ddot{\mathbf{u}}_v \cdot \mathbf{w} \, d\Omega + \int_{\Omega} \mathbf{S}_1 : \nabla^\mathrm{T} \mathbf{u}_v \nabla \mathbf{w} + (C_1 : \mathbf{E}_{\mathbf{u}_v, \mathbf{u}_1}) : (\mathbf{F}_1^\mathrm{T} \nabla \mathbf{w}) \, d\Omega = - \int_{\Omega} \rho \ddot{\mathbf{u}}_1 \cdot \mathbf{w} \, d\Omega \quad \forall \mathbf{w} \in \mathcal{V} \end{array} \right. \qquad (25) $$

with

$$ \mathbf{E}_{\mathbf{u}_v, \mathbf{u}_1} = \frac{1}{2}(\nabla \mathbf{u}_v + \nabla^\mathrm{T} \mathbf{u}_v + \nabla^\mathrm{T} \mathbf{u}_1 \nabla \mathbf{u}_v + \nabla^\mathrm{T} \mathbf{u}_v \nabla \mathbf{u}_1) $$

This problem is therefore linear with respect to $\mathbf{u}_v$. The deformation variational term naturally represents the tangent stiffness matrix of the quasi-static problem [13]. The excitation force of the vibratory problem depends on the quasi-static acceleration $\ddot{\mathbf{u}}_1$. A finite differences scheme is used to approximate this acceleration. The error of approximation induced by the numerical derivation is the result of a rounding error determined by the number of significant digits and a truncation error determined by the discretization step. In addition, the numerical derivation amplifies the numerical noise at high frequencies. Several solutions can be used to treat these different problems. A low-pass filter can be used to filter this noise. Its cutoff frequency is chosen according to the limit of the studied frequency domain. Another strategy consists in interpolating the displacement with derivable functions (polynomials, cubic splines) before carrying out the numerical derivation [31].

The initial conditions of the vibratory problem are extracted from the quasi-static response and the initial conditions of the dynamic problem:

$$ \begin{aligned} \mathbf{u}_v(\mathbf{X}, 0) &= \mathbf{u}_2(\mathbf{X}, 0) - \mathbf{u}_1(\mathbf{X}, 0) \\ \dot{\mathbf{u}}_v(\mathbf{X}, 0) &= \dot{\mathbf{u}}_2(\mathbf{X}, 0) - \dot{\mathbf{u}}_1(\mathbf{X}, 0) \end{aligned} \qquad (26) $$

The initial quasi-static velocity is also approximated with a finite difference scheme. The linear vibratory problem can be solved, after finite element discretization, with a time integration scheme. It depends on the results of the quasi-static problem. The time discretization of the vibratory problem is therefore the same as that used in the quasi-static problem.
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Finite elements model of the tire (a) Three-dimensional finite element model with reference frame (b) Meridian section with measurement points on the inner liner (IL) - Center (•) and Shoulder (•).</figcaption></figure>

<table><caption>**Table 1**

Mechanical properties of the tire.</caption><tr><td>Geometric radius</td><td>R<sub>e</sub> = 268 mm</td></tr><tr><td>Rim radius</td><td>R<sub>i</sub> = 180 mm</td></tr><tr><td>Width</td><td>w = 77 mm</td></tr><tr><td>Young's modulus</td><td>E = 2.10<sup>8</sup> Pa</td></tr><tr><td>Poisson's ratio</td><td>ν = 0.3</td></tr><tr><td>Density</td><td>ρ = 900 kg/m<sup>3</sup></td></tr></table>

## 3. Numerical results

The method developed in Section 2.4 is applied to simulate the vibratory response of a patterned tire rolling on a rigid plane. The validity of this numerical technique has already been discussed in a previous study of the case of a grooved cylinder rolling on a rigid plane [32]. In the following, the mechanical and geometrical properties of the tire structure are given as well as the simulation's numerical parameters. Several types of sensitivity studies will be carried out. Firstly, the rolling conditions influence (i.e. the rolling speed and the deflection) on the tire vibrations response will be studied. This first analysis will give indications on the importance of the operating conditions on the tire's non-linear vibrations. Afterwards, the vibrations' sensitivity with respect to the design parameters of the tire tread pattern will be examined. The influence of the grooves thickness, the inclination angle of the lateral grooves with respect to the lateral direction of the tire structure, and the rib shifting will be more specifically discussed. This second analysis will make it possible to highlight the impact of design modifications on sculpted tires and thus to provide indications on the different possibilities to allow tire manufacturers to better control and attenuate the levels of non-linear vibration responses of tire structures and thus the generation of rolling noise.

The validity and applicability of the strategy proposed will be highlighted by illustrating the possibility of solving the problem by decoupling the amplitudes contributions of the macroscopic rolling and the vibratory response. In the following, the numerical results will be supported by not only sensitivity analysis under different rolling conditions and with respect to certain key tire tread pattern parameters, but also by comparison with experimental results available in the literature.

### 3.1. Tire structure

Figs. 2(a) and (b) illustrate the three-dimensional finite element model studied and a meridian section of the tire, respectively. The geometric and mechanical properties of the tire are given in Table 1. The tire tread (shown in red) is composed of 36 lateral grooves with a thickness of 10 mm. This zone represents the potential contact zone $Γ_c$. The tire tread is composed of a central rib and two shoulder ribs. For the interested reader, the model proposed and the associated physical parameters correspond to data for one conventional tire used on passenger cars.

The finite element model is composed of 15660 nodes, 10944 linear hexahedral elements and 216 linear prismatic elements. The nodes, used for the vibrations response's evaluation, are located on the inner liner (IL), as shown on the tire's meridian section. The first point (shown in red) is located below the middle of the central rib and the second point (shown in blue) is located below the
<table><caption>**Table 2**
Numerical parameters of the simulation.</caption><tr><td>Inflation pressure</td><td>p = 1.5 bar</td></tr><tr><td>Normal penalization</td><td>ε<sub>n</sub> = 10<sup>6</sup></td></tr><tr><td>Time step</td><td>Δt = 3.10<sup>-5</sup> s</td></tr></table>

shoulder rib. These two points are taken at the level of the inner liner since the mesh of the latter remains identical for tires having different tread patterns. The vibrations' evaluation at these two points also allows the study of its evolution with respect to the lateral position. The yellow area represents the inflation zone where a static follower pressure *p* is applied to take into account the tire's inflation. The structure is loaded kinematically by applying a deflection *d* in the vertical direction on the rim zone (shown in blue). The load, carried by the tire, can be calculated by summing the contact pressure over the contact patch (CP). The rim zone, denoted Γ<sub>u</sub>, is also used to drive kinematically the tire rolling with the following boundary condition:

$$ \mathbf{u}(\mathbf{X}, t) = [Vt \ 0 \ -d]^T + \mathbf{R}(\omega t)\mathbf{X} - \mathbf{X} \quad \text{on } \Gamma_u \qquad (27) $$

*V* is the translation velocity and **R** is the rotation matrix associated with the angular velocity *ω*.

$$ \mathbf{R}(\omega t) = \begin{bmatrix} \cos \omega t & 0 & \sin \omega t \\ 0 & 1 & 0 \\ -\sin \omega t & 0 & \cos \omega t \end{bmatrix} \qquad (28) $$

The translation velocity depends on the angular velocity via the rolling radius (also known as the effective radius) *r*<sub>*ω*</sub> = *V*/*ω*. The rolling radius is defined as the radius of the rigid structure having the same rolling speed as the deformable structure. A first order approximation of the rolling radius is given as a function of the geometric radius *R*<sub>e</sub> and the deflection *d* as follows [33]:

$$ r_{\omega} \approx R_{e} - \frac{d}{3} \qquad (29) $$

The tire's material is homogeneous, isotropic and its constitutive law is determined by the Saint Venant-Kirchhoff model.

$$ S = \lambda \operatorname{tr}(\mathbf{E})I + 2\mu\mathbf{E} \qquad (30) $$

*λ* and *μ* denote the Lamé coefficients defined by:

$$ \lambda = \frac{E\nu}{(1 + \nu)(1 - 2\nu)} \qquad \mu = \frac{E}{2(1 + \nu)} \qquad (31) $$

*E* is the Young's modulus and *ν* represents Poisson's ratio. A frictionless contact with a rigid plane is considered and treated using the penalty method. Thus, the impact of the tread pattern is the only mechanism generating vibrations that is taken into account. The air-pumping mechanism and the visco-hyperelastic proprieties of the tire's architecture are not integrated in the model. Therefore, the objective of the following studies is to examine the tread pattern geometry influence on the vibratory response generated by the impact mechanism using the methodology described in Section 2.4.

Regarding the resolution of the problem, the simulation of the quasi-static rolling is performed with the finite elements solver MEF++ [13]. This generic solver was developed jointly by the Groupe Interdisciplinaire de Recherche en Éléments Finis (GIREF) of Université Laval and Michelin. The quasi-static simulation is performed in three successive steps. The first step consists in simulating the inflation by applying a static follower pressure on the inflation zone. The inflated tire configuration is then used to initialize the loading simulation which consists in imposing a vertical displacement on the rim zone. The third step consists in simulating quasi-static rolling by imposing the kinematic condition (27) on the rim zone. The latter step is initialized with the loaded tire configuration. A time integration method with an implicit Euler scheme is then used to solve the vibrations problem (25) given its unconditional stability. A classical convergence study, not presented in this study for the sake of brevity, was first performed in order to calibrate the numerical parameters, such as the mesh size, the normal penalization and the time step, used in the previously described solvers for the resolution of the problem. The selected values of the normal penalization and time step are given in Table 2 (i.e. the simulation results will be invariant with respect to finer mesh, higher values ε<sub>n</sub> and lower values of Δ*t*).

In order to assess the relevance and credibility of both the mathematical modeling proposed for the sculpted tire system and the numerical strategy, the trends obtained with the numerical results will be qualitatively compared with experimental results obtained by Michelin [34] or from the literature. The experimental set-up of the Michelin tire prototype is schematically shown in Fig. 3 with the position of the three measurement microphones at the front and rear of the tire (i.e. microphones located at the inlet and outlet of contact zone defined by points 1 and 2, respectively) and on the lateral side of the tire (point 3). Fig. 4 illustrates the associated acoustic pressure measurements. These experimental results make it possible to observe the evolution of rolling noise as a function of different parameters such as the rolling speed (see all Fig. 4), the grooves thickness (see Fig. 4(a)), the inclination angle of the lateral grooves with respect to the lateral direction of the tire structure (see Fig. 4(b)), and the rib shifting (see Fig. 4(c)). In the following sections, a combined discussion based on numerical findings and comparison with experiments will be proposed to offer a more cohesive understanding of the study's implications and findings. Moreover It should be noted that we are only interested in comparing the evolution of trends as a function of the various parameters, given that the experimental tests were not initially carried
<figure><img src="image_6.png" /><figcaption>Fig. 3. Experimental set-up of the Michelin tire prototype (a) tire used in the measurements and (b) experimental setup.</figcaption></figure>

out in connection with the proposed simulations, and that the tires studied are not necessarily exactly the same in terms of geometry and tread pattern.

## 3.2. Influence of rolling speed

The influence of rolling speed on the vibratory response is performed by considering three cases of quasi-static rolling simulations with three angular velocities $\omega = \{50, 100, 150\}$ rad/s for a deflection of $d = 10$ mm. These three angular velocities correspond to the translation velocities $V = \{48, 96, 144\}$ km/h calculated using the approximation of the rolling radius (29). Figs. 5 and 6 illustrate the frequency spectrum of the vibrations velocity components calculated with the three angular velocities. The spectra are plotted as a function of the order which is equal to the frequency divided by the rotation frequency. It is noteworthy that using the order in the spectral analysis permits normalizing the results with respect to the rotation speed. For the three tested configurations, the frequency spectra are mainly composed of harmonics located around the number of lateral grooves multiplied by the rotation frequency (which corresponds to Order = 36 in Figs. 5 and 6). It can be observed that the amplitudes of speed according to the three directions differ according to the spatial directions and the measurement point on the inner liner of the tire. Indeed, the major contributions correspond to vibratory responses according to directions x and z whereas the amplitudes of speed according to direction y are very clearly lower. It also notable that the lateral vibration velocities $v_y$ calculated at the center of the Inner Liner (IL) are almost null whatever the rolling velocity, which is a result that was expected due to the lateral symmetry of the tire structure. However, the levels of the lateral vibrations $v_y$ calculated on the shoulder of the IL are present, non-null but of low amplitudes in comparison to the other directions x and z (to be compared with the levels of lateral vibrations $v_x$ and $v_z$). They increase with the augmentation of the rotation speed. Thus it appears that the proposed numerical methodology efficiently separates the two scales of the dynamic response of the tire and highlights the specificity of the tire's vibratory behavior resulting from the impact between the pattern of the tire tread and the rigid plane. On examining more specifically the influence of the rolling speed, it seems that the levels of vibrations are very low for $\omega = 50$ rad/s compared to the results for the rolling simulations at $\omega = \{100, 150\}$ rad/s. The levels of the contributions of the different harmonics for the two configurations at $\omega = \{100, 150\}$ rad/s are generally similar for the two directions x and z. It can also be seen that the harmonics appearing around the grooves impact frequency are more numerous for $\omega = 150$ rad/s in comparison to $\omega = \{50, 100\}$ rad/s.

For interested readers, it has been shown experimentally that rolling speed is one of the most important parameters having an influence on the external rolling noise. As the vehicle speed and the number of tire revolutions increase, the speed with which the tire comes in contact with the road surface increases resulting in higher tire/road impact and also an increased pumping in and out of air from the contact patch edge [35]. Fig. 4 illustrates the evolution of rolling noise as a function of the rolling speed. The rolling noise level increases noticeably with increasing rolling speed between 20 km/h and 80 km/h, whatever the tire tread pattern parameters. A saturation of the rolling noise is then observed for speeds around 90 km/h. These experimental results are in agreement with the numerical results. Thus increasing rotation speed shifts vibration levels towards higher frequencies, as the frequency of impact between the tire tread and the ground depends on rotation speed. Simulation results and experimental measurements show that rolling speed has a significant impact on vibration levels.

## 3.3. Influence of deflection

A second parameter of interest during tire rolling is the load carried by the tire. Indeed, the load carried mainly determines the shape of the contact patch and the mechanisms generating external noise are located around the contact patch.
<figure><img src="image_7.png" /><figcaption>Fig. 4. Evolution of the acoustic pressure measured experimentally [34] as a function of velocity for different (a) groove thickness values (b) inclination angles $\psi_2$ of the lateral grooves with respect to the lateral direction of the tire and (c) rib shiftings.</figcaption></figure>

The simulation of the vibratory response of the tire is performed for three different deflections $d = \{5, 10, 15\}$ mm by considering a rolling speed $V = 96$ km/h. It should be recalled that the tire's loading is simulated kinematically by imposing a vertical deflection on the rim zone (the issue of modeling quasi-static rolling with an imposed load will be dealt with in Section 4).

The frequency spectra of the vibration velocities of the node located on the shoulder of the IL are shown in Fig. 7. In the vicinity of the harmonic contributions corresponding to a number of lateral grooves multiplied by the rotation frequency (i.e. order 36), it appears that the contributions are lower for the case corresponding to the largest deflection (i.e. $d = 15$ mm) whereas the contributions for the two other cases considered (deflection of $d = 10$ mm and $d = 5$ mm) are quite similar. In addition, the amplitude of the harmonics around twice the impact frequency (i.e. order = 72 in Fig. 7) becomes comparable to that of the harmonics appearing around order 36 for the case of a deflection $d = 5$ mm, a phenomenon not seen for the two other configurations. To better understand this tendency and the physical phenomena at the level of the contact patch, the shape of the CP obtained for each deflection is examined. Fig. 8 shows the distribution of the normal contact pressure calculated for each of the three configurations. In the case of $d = 5$ mm, the grooved central rib is the only element of the tire in contact with the rigid plane. On the contrary, the contact zone estimated with a deflection $d = 10$ mm shows that the three ribs (centered rib and the shoulder ribs on each end) carry the load, naturally with a contact symmetry in relation to the center of the tire tread. Regarding the configuration with a deflection $d = 15$ mm, it can be noted that the vertical load is mainly carried by the shoulder ribs. For those readers interested, the frequency spectra also show that the lateral grooves of the center rib generate most of the vibrations (results not shown for the sake of concision). Thus, the tendencies
<figure><img src="image_8.png" /><figcaption>Fig. 5. Frequency spectrum of the vibration velocity of the central node of the IL for $V = 48$ km/h (■), $V = 96$ km/h (■) and $V = 144$ km/h (■).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Frequency spectrum of the vibration velocity of the shoulder node of the IL for $V = 48$ km/h (■), $V = 96$ km/h (■) and $V = 144$ km/h (■).</figcaption></figure>

given by the model can be explained by the fact that the support added by the non-grooved shoulder ribs at each end damps the vibration levels generated by the lateral grooves.

Experimental tests have shown that, generally speaking, a larger contact area increases noise levels [11]. Similarly, the work of Iwao and Yamazaki [36] illustrated that acoustic sound pressure generated from tire rolling reaches a peak for a given load value when the speed of rotation increases. It may also be noted that in some studies, it is analysed that this increase remains marginal in comparison with other factors [37]. On this issue, although the numerical studies clearly show that the evolution of normal contact pressure in regard to different load carried by the tire can be complex (see Fig. 8), the trends in nonlinear response levels are not in complete agreement with the above experimental observations. While the contributions for the two configurations tested numerically (i.e. $d = 10$ mm and $d = 5$ mm) are fairly similar (in line with Tong's experimental observations [37]), the contributions of nonlinear vibrations are lower for the case corresponding to the greatest deflection (i.e. $d = 15$ mm with the smallest contact surface), which is
<figure><img src="image_10.png" /><figcaption>Fig. 7. Frequency spectrum of the vibration velocity of the shoulder node of the IL for $d = 5$ mm (■), $d = 10$ mm (■) and $d = 15$ mm (■).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Normal contact pressure calculated with a deflection $d = \{5, 10, 15\}$ mm.</figcaption></figure>

not in agreement with the experimental observations of Li et al. [11]. It should be noted that the issue of appropriate mathematical modeling of quasi-static rolling with an imposed load will be addressed in Section 4.

### 3.4. Sensitivity with respect to the grooves thickness

This part is devoted to studying the influence of the grooves thickness on the tires vibrations response. The studied tread patterns are shown in Fig. 9. A deflection $d = 10$ mm and a rolling speed of $V = 96$ km/h are taken as rolling conditions. The vibrations response of the nodes located on Inner liner (IL) will be compared and analyzed since this part of the tire structure remains invariable in the tires meshes used for this sensitivity study (contrary to the nodes located on the tire tread, indicated previously in red in Fig. 2).

Figs. 10 and 11 show the vibrations velocities frequency spectra of the center and shoulder nodes, respectively. The spectra of the longitudinal and vertical vibration velocities $v_x$ and $v_z$ are almost the same and composed of harmonics around the number of grooves of the wheel with a modulation effect. In addition, a reduction of the harmonic contributions can be seen (in the vicinity of
<figure><img src="image_12.png" /><figcaption>Fig. 9. Tread patterns used in the study of the sensitivity of the vibratory response in relation to lateral groove thickness.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Frequency spectrum of the vibration velocity of the center node of the IL for $e = 5$ mm (■), $e = 10$ mm (■), $e = 15$ mm (■) and $e = 20$ mm (■).</figcaption></figure>

the harmonics corresponding to a number of lateral grooves multiplied by the rotation frequency, i.e. order = 36) with an increase of groove thickness. This reduction of the levels of amplitude of each harmonic contribution can be explained physically by the reduction of the impact surface on the rigid plane due to the increase in the thickness of the grooves. It is also notable that the contributions linked to the lateral vibrations $v_y$ at the center (i.e. Fig. 10) are null, which is physically consistent given the lateral symmetry of tire structures. On the contrary, the harmonic contributions of the lateral vibrations $v_y$ calculated at the end of the tire (i.e. shoulder rib associated with Fig. 11) are non-null, which can be explained physically by the fact that the lateral vibrations follow the same tendencies as the longitudinal $v_x$ and vertical $v_z$ vibrations calculated at the center due to the Poisson effect.

The experimental results, illustrated in Fig. 4(a), show trends that are not necessarily consistent with those obtained with the numerical model. This can be explained by the air pumping mechanism, which tends to increase with larger grooves, since the quantity of air captured in the contact area increases with the size of the grooves. It should be noted on this point that the experimental measurements found in the literature do not give a clear trend on the evolution of external noise in relation to the thickness of the lateral grooves [6,11] because the impact and air pumping mechanisms are in competition. Similarly, it seems difficult to understand these trends since the variation in the load carried by the tyre induced by the variation in thickness is not negligible. Moreover, the mathematical modeling of quasi-static rolling with an imposed load will (as proposed in the following in Section 4) would be more appropriate.

### 3.5. Sensitivity with respect to the angle $\psi_2$

A parameter used classically in tire tread design is the inclination angle of the lateral grooves with respect to the lateral direction of the tire structure. This angle is illustrated in Fig. 12 by angle $\psi_2$. In what follows, the influence of this parameter on the tire
<figure><img src="image_14.png" /><figcaption>Fig. 11. Frequency spectrum of the vibration velocity of the shoulder node of the IL for $e = 5$ mm (—), $e = 10$ mm (—), $e = 15$ mm (—) and $e = 20$ mm (—).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Tire treads used to study the sensitivity of the vibratory response in relation to the inclination angle of the lateral grooves.</figcaption></figure>

vibratory response will be examined by considering the four tire treads configurations shown in Fig. 12. The rolling conditions remain unchanged compared to the previous part (i.e. $V = 96$ km/h and $d = 10$ mm). The vibrations responses are compared at the level of the points of the inner liner.

Figs. 13 and 14 display the evolution of harmonic contributions in the frequency domain of vibration velocities at the center node and shoulder node, respectively, for the four configurations of the lateral groove inclination angle.

A global attenuation of harmonic amplitudes is observed by increasing the inclination angle values $\psi_2$ whatever the direction (with however amplitudes for $v_y$ weaker than those for the longitudinal $v_x$ and vertical $v_z$) velocities. This attenuation of the harmonics amplitudes with respect to the inclination angle $\psi_2$ can be physically explained by the progressive entries and exits of the grooves in the contact patch contrary to the lateral grooves ($\psi_2 = 0$) where the entry and exit occur instantaneously. When focusing more specifically on the configuration where $\psi_2 = 0$ and the lateral velocity $v_y$, it appears that the harmonic contributions are very weak: more precisely we find a null contribution for the evaluation point located at the center of the IL and a weak contribution for evaluation point located at the shoulder of the IL. These results are in perfect agreement with the previous observations and analysis proposed in the studies of the previous sections. In addition, the presence of the contribution of non-negligible harmonics of the lateral velocity $v_y$ for configurations where $\psi_2 > 0$ can be explained physically by the breakup of the tire lateral symmetry.

The experimental measurements shown in Fig. 4(b) and those found in the literature [6,11] show that the acoustic noise follows the same trend as the numerical results with respect to the inclination angle $\psi_2$ of the lateral grooves with respect to the lateral direction of the tire. In addition, these results can be explained physically by the fact that the air pumping mechanism at the tire has the same sensitivity to angle $\psi_2$ as the impact mechanism, since the compression (respectively the relaxation) of the air at the inlet (the outlet, respectively) of contact zone occurs more gradually as the value of angle $\psi_2$ increases.
<figure><img src="image_16.png" /><figcaption>Fig. 13. Frequency spectrum of the velocity f the vibrations of the center node of the IL for $\psi_2 = 0$ deg (—), $\psi_2 = 22$ deg (—), $\psi_2 = 35$ deg (—) and $\psi_2 = 45$ deg (—).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Frequency spectrum of the velocity of vibrations of the shoulder node of the IL for $\psi_2 = 0$ deg (—), $\psi_2 = 22$ deg (—), $\psi_2 = 35$ deg (—) and $\psi_2 = 45$ deg (—).</figcaption></figure>

### 3.6. Sensitivity with respect to rib shifting

The last parametric analysis proposed concerns the notion of “rib shifting” which considers the case where the shoulder ribs are also grooved. The problem therefore consists in determining the impact of the alignment of the grooves of the shoulder ribs with respect to the grooves of the center rib. In what follows, the shoulder ribs are composed of 36 lateral grooves and three configurations are studied, as shown in Fig. 15: the grooves of the three ribs (center and two ends) are aligned in the first configuration; a shift of 25% is applied to one of the shoulder ribs in the second configuration (as shown in the image at the center of Fig. 15 for the rib of the upper end); the third configuration consists in keeping the shift of 25% for the rib of the upper end while applying a shift of 50% on the other shoulder rib (as shown in the right hand image of Fig. 15 for the rib of the lower end). The simulation of the vibratory
<figure><img src="image_18.png" /><figcaption>Fig. 15. Treads used in the study of the sensitivity of the vibratory response with respect to rib shifting.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. Frequency spectrum of the vibration velocity of the center node of the IL of config. 1 (—), of config. 2 (—) and of config. 3 (—).</figcaption></figure>

response is done using the same rolling conditions as for the previous studies (i.e. $V = 96$ km/h, $d = 10$ mm). The vibratory behavior is analyzed by observing the vibration velocity at the measurement points located on the inner liner.

Figs. 16 and 17 illustrate the evolution of harmonic contributions in the frequency domain of vibration velocities at the center node and shoulder node, respectively, across the three ‘rib shifting’ configurations. It can be seen that the frequency spectra are very similar for the three configurations, even though a slight attenuation of the levels of harmonic contributions for configurations 2 and 3 can be detected for the longitudinal $v_x$ and vertical $v_z$ vibrations. This attenuation can be explained by the shift created by the rib shifting at the entry and exit of the grooves in the contact zone. Indeed, the aim of the rib shifting is to break the synchrony of the grooves impact of the three ribs on the road. Furthermore, it seems that the introduction of “rib shifting” generates lateral velocities $v_y$ at the center of the wheel (contrary to all the studies performed in the previous sections for which the harmonic contributions of the lateral velocities $v_y$ are null). This can be easily explained physically by the fact that the presence of “rib shifting” breaks the lateral symmetry of the tire. It can be noted that the contributions of the lateral velocity $v_y$ estimated on the shoulder rib for the two shifted rib configurations (i.e. configurations 2 and 3) are higher than those calculated with the tread with aligned lateral grooves (i.e. configuration 1).

The experimental measurements shown in Fig. 4(c) indicate that the rib shifting attenuates rolling noise levels in the case of 50% shift while it has a negligible impact in the case of 25% shift. The numerical results show that the rib shifting induce a slight decrease of the vibrations. As noted in some studies [6,11], the degree of attenuation depends also on driving and operating conditions and other tire tread parameters. As already explained, this can be physically interpreted by the fact that the attenuation of rolling noise may be due to the fact that the offset of the ribs results in a desynchronization of the impact of the ribs at the entry and exit of the grooves in the contact zone.
<figure><img src="image_20.png" /><figcaption>Fig. 17. Frequency spectrum of the vibration velocity of the shoulder node of the IL of configuration 1 (—), of configuration 2 (—) and of configuration 3 (—).</figcaption></figure>

## 4. Rolling under imposed load

The tire loading simulation was performed in what preceded by imposing a vertical displacement on the rim zone. This control of the quasi-static rolling via the deflection can be considered as being an approximation of the real operating conditions of a tire given that tire systems roll while bearing the load of the vehicle. However, it should be remembered that the fact of imposing a vertical force on the rim zone makes the mechanical problem presented previously ill-posed.

Therefore, if one wants to get closer to more realistic operating conditions that takes into account the load of a vehicle and integrate its impact on the tire's vibratory behavior, it is necessary to introduce an additional equation in the mathematical formulation proposed previously to model the system, so as to calculate the deflection, ensuring that the load carried by the tire is equal to the target load at each moment. In what follows, we first propose discussing the ill-posed nature of the mechanical problem due to taking into account a vertical force applied on the tire rim zone. Secondly, the additional equation permitting the integration of the load in the problem is explained by emphasizing in particular the modifications made in the formulation of the problem to control rolling with an imposed load. Then, the method is applied to simulate quasi-static rolling of the tire while taking account an imposed load. To this end, the vibratory response is calculated by considering this new mathematical formulation using the strategy developed previously in 2.4.

### 4.1. Ill-posed problem

The notion of a well-posed problem in the meaning of its mathematical formulation was introduced by Hadamard [38] as early as 1923. It is a problem for which the solution exists (existence of the solution of the problem), the solution is unique (uniqueness of the solution of the problem) and it depends continually on data (stability of the solution of the problem). In what follows, we analyze the formulation associated with taking into account an imposed load and determining the associated mathematical formulation to be applied to treat this type of problem. We consider a vertical force $F_z$ applied on the tire rim zone. It is noteworthy that the introduction of this new operating condition under an imposed load naturally substitutes taking into account an imposed deflection chosen initially (as proposed in the previous sections of this work. The longitudinal $u_x$ and lateral $u_y$ displacements are assumed to be null on the rim zone $\Gamma_u$ as shown in Fig. 18 (rim zone identified by the red circle of the tire). The solution of the problem $\mathbf{u}$ is therefore sought in the functional space $\mathcal{V}_1$ defined by:

$$ \mathcal{V}_1 = \{ \mathbf{u} \in H^1(\Omega) \mid u_x = u_y = 0 \ \text{sur} \ \Gamma_u \} \qquad (32) $$

The contact condition can be interpreted as a Neumann condition where the contact forces are applied on the contact zone [39]. Thus, the variational formulation of the equations of motion modeling the static loading of the tire on a rigid plane under the effect of force $F_z$ is given by:

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{u} \in \mathcal{V}_1 \\ \displaystyle \int_{\Omega} \mathbf{P} : \nabla \mathbf{w} \, d\Omega = \int_{\Gamma_p} pJ(F^{-1}N) \cdot \mathbf{w} d\Gamma + \int_{\Gamma_u} F_z w_z \, d\Gamma + \int_{\Gamma_c} \Pi_n w_n \, d\Gamma \quad \forall \mathbf{w} \in \mathcal{V}_1 \end{array} \right. \qquad (33) $$
<figure><img src="image_21.png" /><figcaption>Fig. 18. Mapping of crushing with load $F_z$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. Parameterization of the ground position.</figcaption></figure>

If $\mathbf{u}$ is a solution of the variational problem, then the function $\mathbf{u}_e = \mathbf{u} + c\mathbf{e}_z$ is also a solution of the problem (with $c$ being a non-null constant and $\mathbf{e}_z$ the vertical direction unit vector):

* The kinematic boundary conditions imposed on $\Gamma_u$ are verified by $\mathbf{u}_e$ since its longitudinal and lateral components are the same as those of $\mathbf{u}$.

* The boundary conditions on the follower pressure $p$ and on the contact forces are verified by $\mathbf{u}_e$ since its gradient is the same as that of $\mathbf{u}$.

Thus, the solution of the problem (33) is not unique and therefore it is ill posed. Remember that this type of analysis of non-uniqueness of a solution is a classical conclusion in solid mechanics. Indeed, it corresponds to the fact that the solution of a mechanical problem is characterized by a rigid body motion when a boundary condition in displacement is not imposed on all the directions of the considered space. This difficulty can be bypassed by keeping control of the tire loading with the imposed deflection and by adding an additional equation to the mechanical problem. This technique will be detailed in the following section.

## 4.2. Load equation

The vertical load borne by the structure is calculated by integrating the normal pressure in the contact patch. To express the gap function gap $g_n$ analytically as a function of the deflection $d$, the nodal displacement on the rim zone is null and the position of the ground $z_p$ is variable and equal to the value of the deflection, as shown in Fig. 19. By considering the penalization method chosen to manage the non-interpenetration condition, load $F_z$ can be calculated as follows:

$$F_z(\mathbf{u}, d) = \int_{\Gamma_c^t(d)} p_n d\Gamma = \int_{\Gamma_c^t(d)} \varepsilon_n g_n d\Gamma = \int_{\Gamma_c^t(d)} \varepsilon_n (g_0 + u_n - d) d\Gamma \quad (34)$$
<table><caption>**Table 3**

Algorithm simulating quasi-static rolling with an imposed load.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1 - Initialization</td><td>Initialization of displacement and deflection using the loading calculation (u, d)<br>Choice of convergence thresholds (ε₁, ε₂) and the maximum number of iterations N</td></tr><tr><td>2 - Time loop at time t<sub>k</sub></td><td>Iterative loop i<br>Calculate u<sub>k</sub> solution of the quasi-static problem at time t<sub>k</sub> with deflection d<sub>i</sub><br>Calculate load F<sub>z</sub>(d<sub>i</sub>)</td></tr><tr><td></td><td>Convergence test<br>If (|F<sub>z</sub>(d<sub>i</sub>) - C<sub>z</sub>| < ε₁) or (|Δd<sub>i</sub>| < ε₂) or (i = N)<br>Exit loop<br>Otherwise</td></tr><tr><td></td><td>d<sub>i+1</sub> = d<sub>i</sub> + Δd<sub>i</sub><br>Repeat (i = i + 1).</td></tr><tr><td></td><td>Continue to time t<sub>k+1</sub></td></tr></tbody></table>

with $g_0$ being the initial normal gap and $u_n$ the normal displacement component. Rolling is controlled with the boundary conditions defined previously in equation (27) but by considering in addition a variable deflection dependent on time. This corresponds to the zero of the following additional equation Eq(d):

$$Eq(d) = (F_z(d) - C_z)^2 \qquad (35)$$

where $C_z$ represents the imposed load. Newton's method is used to find the zero of equation (35). Thus, after having chosen an initial value of the deflection $d_0$, it is updated with the following iterative scheme:

$$d_{i+1} = d_i + \Delta d_i = d_i - \frac{Eq'(d_i)}{Eq(d_i)} \qquad (36)$$

The calculation of the derivative of the additional equation (i.e. Eq'(d)) introduced in the previous equation is therefore necessary when implementing the iterative scheme given by the Newton method. An approximation of this derivative can be obtained using the analytical expression of the gap function $g_n$. It becomes:

$$Eq'(d) = 2(F_z(d) - C_z) \frac{\partial F_z}{\partial d} \approx -2\epsilon_n |\Gamma_c^t(d)| (F_z(d) - C_z) \qquad (37)$$

where $|\Gamma_c^t(d)|$ represents the surface of the contact patch. Note that the calculation of the derivative of the expression of the gap in relation to the deflection is straightforward. However, the derivation of the integration domain $\Gamma_c(d)$ in relation to the deflection is more complicated to perform (this part of the derivative is not taken into account by assuming that it is negligible in the vicinity of the solution of the problem being dealt with). The solution of the load equation is performed with a quasi-Newton method. Its convergence, with a lower rate than that of the Newton method, remains influenced by the initialization of the deflection. The algorithm simulating the quasi-static rolling with an imposed load is given in Table 3.

By analyzing the structure of the algorithm proposed in Table 3, it can be seen clearly that the computation time of the quasi-static rolling simulation with an imposed load is longer than that with an imposed deflection since several resolutions of the quasi-static problem are performed at each time step because the Newton loop is embedded in the quasi-Newton loop of the deflection. A weak coupling between displacement u and deflection d is therefore considered. It can be seen that the convergence of the iterative loop can be accelerated by multiplying the correction $\Delta d_i$ by a coefficient $\alpha > 1$.

### 4.3. Numerical application

The calculation method and strategy presented in the previous section are applied to simulate the quasi-static rolling under imposed load of the tire shown in Fig. 2. The mechanical properties and the numerical parameters of the simulation are the same as those used for the numerical simulations performed previously. A load of 500 daN and a rolling speed of $V = 96$ km/h are imposed as operating conditions on the tire system.

Fig. 20 shows the temporal evolution of the load borne by the tire, the temporal evolution of the deflection and the frequency spectrum of the variation of the deflection (see the blue curves for the present case, referred to as “initial rolling condition”). It appears that the load borne by the tire is equal to the imposed load (as illustrated in Fig. 20(a)). It can also be mentioned that variable $d$ evolved in time (see Fig. 20(b)) which characterizes the fact that deflection $d$ varies through time to conserve a vertical load equal to the target load during the resolution of the problem. In addition, this variation can be explained physically by the variation of the shape of the contact patch due to the lateral grooves composing the tread pattern. The frequency spectrum of the variation of deflection $\delta d$ defined as being the difference between deflection $d$ and its average value $\bar{d}$ is composed of harmonics around the multiples of the impact frequency of the grooves on the contact patch (order = 36k, $k \in N$) with a strong predominance of the principal order ($k = 1$), as shown in Fig. 20(c)).

Figs. 21 and 22 show the evolution of the vibratory response of the tire (at the center and shoulder nodes, respectively) whose quasi-static rolling is performed with an imposed load (see the blue curves for the present case, referred to as “initial rolling condition”).
<figure><img src="image_23.png" /><figcaption>Fig. 20. Time evolution (a) of the load borne by the tire (b) of the deflection. (c) Frequency spectrum of the deflection variation (—) with the initial rolling condition (—) with the new rolling condition (to be noted that the blue and red curves are identical).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 21. Time evolution of the vibratory displacement of the center node of the IL with the initial rolling condition (—) and with the new rolling condition (—).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 22. Time evolution of the vibratory displacement of the shoulder node of the IL with the initial rolling condition (—) and with the new rolling condition (—).</figcaption></figure>

A numerical amplification (nonphysical) in time is observed at the level of the longitudinal vibrations $u_x$ and the lateral vibrations $u_y$, directly linked to the model of the quasi-static rolling. To investigate the source of this amplification, the quasi-static acceleration representing the excitation signal of the vibratory problem is examined at the level of the two evaluation points chosen previously (i.e. center node and shoulder node of IL), as shown in Figs. 23 and 24. An amplification of the longitudinal quasi-static acceleration $a_x$ can be seen clearly at the two evaluation points. This amplification effect was not however observed on the lateral accelerations $a_y$. Thus, it can be concluded that the amplification of the lateral vibrations $u_y$ (blue curves given in Figs. 21 and 22) are the consequence
<figure><img src="image_26.png" /><figcaption>Fig. 23. Time evolution of the quasi-static acceleration of the center node of the IL with an imposed load of $C_z = 500$ daN.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 24. Time evolution of the quasi-static acceleration of the shoulder node of the IL with an imposed load of $C_z = 500$ daN.</figcaption></figure>

of the Poisson effect. This amplification can be explained moreover by the condition used to control the advance of the tire in time with an imposed load. To demonstrate this, we propose in what follows to explain and discuss the expressions of the longitudinal displacement $u_x(x, y, z, t)$ (imposed on the rim zone) and its derivatives $\dot{u}_x(x, y, z, t)$ and $\ddot{u}_x(x, y, z, t)$. The expressions $u_x, \dot{u}_x$ and $\ddot{u}_x$ are given by:

$$ 
\begin{aligned}
u_x(x, y, z, t) &= (R_e - \frac{d(t)}{3})\omega t + x \cos(\omega t) + z \sin(\omega t) - x \\
\dot{u}_x(x, y, z, t) &= (R_e - \frac{d(t)}{3})\omega - \frac{\dot{d}(t)}{3}\omega t - \omega x \sin(\omega t) + z \omega \cos(\omega t) \\
\ddot{u}_x(x, y, z, t) &= -\frac{2\dot{d}(t)}{3}\omega - \frac{\ddot{d}(t)}{3}\omega t - \omega^2 x \cos(\omega t) - z \omega^2 \sin(\omega t)
\end{aligned}
\quad (38) $$

where x, y and z represent the coordinates of the position vector X of a point of the reference configuration. Equation (38) clearly shows that the variation of the deflection amplifies the quasi-static longitudinal acceleration imposed in the control of the quasi-static rolling. Indeed, this amplification is characterized by the terms $\frac{d(t)}{3}\omega t$ and $\frac{\dot{d}(t)}{3}\omega t$ for the contribution in velocity $\dot{u}_x(x, y, z, t)$ and acceleration $\ddot{u}_x(x, y, z, t)$, respectively. This amplification on the longitudinal component of the velocity and acceleration therefore impacts the quasi-static acceleration of the entire structure. Also, to correct this non-physical amplification, a simple strategy to be implemented consists in replacing the deflection $d(t)$ by its average value $\bar{d}$, which permits getting rid of the time dependency of the deflection variable in relation to time $t$ (and thus canceling the two terms $d(t)$ and $\dot{d}(t)$ and avoiding the numerical increase in time of the velocity and the longitudinal acceleration). The new longitudinal displacement $u_x(x, y, z, t)$ imposed on the rim zone is therefore given by:

$$ u_x(x, y, z, t) = (R_e - \frac{\bar{d}}{3})\omega t + x \cos(\omega t) + z \sin(\omega t) - x \quad (39) $$

The results of the correction made to the kinematic condition imposed to control the quasi-static rolling are given in Fig. 20 (see the red curves for this new kinematic condition, referred to as “new rolling condition”). The load carried out by the tire is equal to the imposed load. The evolution of deflection $d$ and the frequency spectrum of the variation of deflection $\delta d$ are the same as those obtained with the initial rolling condition. Consequently, the correction permits keeping the same contact patch while eliminating the numerical amplification of the quasi-static acceleration induced by the variation of the deflection. This therefore validates the mathematical formulation and the proposed numerical strategy.

The vibratory response of the tire under imposed load is now calculated with the resolution strategy proposed and using the quasi-static response obtained with the new control condition of rolling. Results are provided in Figs. 21 and 22 for displacements according to three directions for the node located at the center and at one end of the tire of the IL, respectively. These figures illustrate
<figure><img src="image_28.png" /><figcaption>Fig. 25. Frequency spectrum of the vibration velocity of the center node of the IL with an imposed load $C_z = 500$ daN.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 26. Frequency spectrum of the vibration velocity of the shoulder node of the IL with an imposed load $C_z = 500$ daN.</figcaption></figure>

the comparison between the vibratory responses by considering the initial formulation of the rolling condition (blue curves) and the modified proposal (red curves), for displacements according to three directions for the node located at the center and at one end of the tire of the IL, respectively.

First, it can be clearly observed that the new rolling condition proposed via the modified longitudinal displacement permits eliminating the numerical amplification observed at the level of the longitudinal vibrations and at the level of the lateral vibrations calculated at the center and at the shoulder of the IL. Moreover, it can be noted that it leads to an attenuation of the amplitude of the lateral $u_y$ and longitudinal $u_x$ vibrations while the vertical vibrations $u_z$ are slightly affected. It can also be seen that the vibration frequency does not vary due to the correction made to the rolling condition. Thus, the correction of the rolling condition makes it possible to obtain a periodic vibratory behavior consistent with the physical phenomena integrated in the model of the tire system (i.e. loaded structure rolling on a rigid plane).

The frequency spectra of the vibration velocities calculated with an imposed load (and by considering the new rolling control condition) are given in Figs. 25 and 26 for the center and shoulder nodes of the IL, respectively. These results should be compared
with those proposed in the previous sections (results with the imposed deflection condition). The main difference in comparison to the simulations with the imposed deflection condition is observed for the vertical vibration frequencies $v_z$. Indeed, the main contribution can be observed around the grooves impact frequency (order = 36) with less marked secondary contributions identified as harmonics located around the multiples of the grooves impact frequency (order = $36k, k \in \mathbb{N}^*$). It can also be seen that the maximum contribution is identified for the vertical velocity with amplitudes of maximum contributions approximately 3 times higher than those detected for the transversal velocity $v_x$, whereas for the imposed deflection configuration the intensities of harmonic contributions were quite similar for vertical $v_z$ and longitudinal $v_x$ velocities. This change of tire behavior compared with the imposed deflection configuration can be explained by the variation of the vertical quasi-static displacement brought about by the variation of the deflection. In addition, the amplitudes ratio between the vertical velocity $v_z$ and the longitudinal $v_x$ velocity corresponds to the impact of the deflection on the rolling radius defined in Equation (29).

Regarding the longitudinal velocities of $v_x$ it can also be noted that they do not vary greatly, for the two measurement points considered, in comparison to the calculations with an imposed deflection, and that the global composition of the associated spectra is mainly composed of harmonics located around frequency of groove passages in the contact patch. To conclude, it can be observed that the lateral velocity $v_y$ is canceled at the center of the IL which can be explained as stated previously by the tire's lateral symmetry. However, it is interesting to note the presence of harmonics around the multiples of contact efforts in the case of this new configuration under imposed load with an effect of modulation in the spectrum of lateral vibration velocities $v_y$.

In conclusion, these numerical results, which allow considering the case of a tire through controlled quasi-static rolling with an imposed load, show the non-negligible effect of the choice of model of the problem of a tire rolling on a rigid plane. It appears without any possible doubt that a considerable change of tire vibratory behavior occurs regarding in particular vertical vibrations $v_z$ both in terms of frequency content and the intensity of different preponderant harmonic contributions. This evolution of vibratory behavior in comparison to the imposed deflection configuration is consistent with the variation of the shape of the contact patch between the two configurations tested (i.e. in deflection and imposed load, respectively). We recall that the latter configuration tested under imposed load normally corresponds to a more realistic case of tire loading and, due to the differences of behavior in comparison to the first configuration tested, it should be privileged for the tire's rolling simulation.

### 4.4. CPU time

The aim of this last part is to analyze the computation time of the different simulations presented in this paper. This point is fundamental for industrial applications given that the size of the finite element models treated is considerable (degrees of freedom $\approx 10^6$). The aspects of the computation sequence used to calculate the rolling noise developed in this paper are shown in Fig. 27. The results examined in this article concern the part framed in red in Fig. 27.

The architecture of the MEF++ solver makes it possible to perform high performance calculations (HPC). Indeed, the assembly of matrices and the resolution of the linear systems (PETSc library) required in the computation of quasi-static rolling are parallelized. This permits the efficient use of computer resources and therefore reducing computation time. The assembly of the stiffness matrix of the vibratory problem, performed at each time step because it depends on quasi-static displacement, is also paralleled in the script developed on Matlab to solve the vibratory problem. The computation time, examined in what follows, concerns the times of the simulations performed with the mesh shown in Fig. 2. The final simulation time corresponds to 6000 time steps and permits the structure to make three wheel rotations at a rolling speed of $V = 96$ km/h. A line search algorithm [40] is added to the quasi-static Newton calculation loop to increase the robustness of convergence.

Table 4 shows the different CPU times obtained for the two steps corresponding to the quasi-static rolling (MEF++) and to the evaluation of the vibrations (Matlab) in the case of the imposed deflection configuration. It can be noted that the scientific developments proposed in terms of mathematical modeling and computer methods, as well as the calculation times announced to solve complex engineering problems are compatible with use in the design phase for manufacturers working on sculpted tire design problems. Thus the proposed methodology can serve as an aid to a better understanding of complex and realistic problems in the multi-scale nature of the dynamic response of an industrial rolling sculpted tire.

It can be seen very clearly that increasing the number of computation cores used leads to a considerable reduction in the calculation time of the quasi-static rolling and vibratory response simulations. It is almost proportional to the number of cores used in the case of quasi-static rolling. This proves the scalability of the MEF++ solver. The decrease in the computation time of the vibratory response is less substantial than that of the quasi-static rolling but nonetheless remains non-negligible. This all demonstrates without ambiguity the potential of the strategy proposed in this study for utilization in an industrial context and on large-scale finite element models.

Also, Table 5 shows the evolution of the computation time of quasi-static rolling under imposed load as a function of the number of computation cores. Increasing the number of cores also leads to considerable reduction of the computation time. Computation times with rolling under an imposed load are far larger than those of rolling under an imposed deflection. This can be expected since the Newton loop of the quasi-static calculation is embedded with the quasi-Newton loop that updates the deflection. The latter ensures that the load borne by the tire is equal to the target load $C_z$ at each time step.

In conclusion, the computation times announced previously for the two rolling configurations in imposed deflection and imposed load show without the least ambiguity the potential of the strategy proposed in this study for utilization in an industrial context and on large scale finite element models. The evolution of the calculation time as a function of the computer resources used in the simulation shows the importance of parallelization. For interested readers, it is clear that the computation times could be optimized still further by integrating the calculation of the vibratory behavior in the MEF++ solver since the stiffness matrix of the vibratory problem is the tangent stiffness matrix of the quasi-static problem. All the conclusions regarding the CPU times are encouraging in
<figure><img src="image_30.png" /><figcaption>Fig. 27. Rolling noise computing sequence.</figcaption></figure>

<table><caption>**Table 4**
Evolution of computation time as a function of the number of cores - Rolling under imposed deflection.</caption><thead><tr><th></th><th>1 core</th><th>6 cores</th><th>12 cores</th></tr></thead><tbody><tr><td>Quasi-static rolling (MEF++)</td><td>72 h</td><td>14 h</td><td>7 h 30 min</td></tr><tr><td>Vibrations (Matlab)</td><td>27 h</td><td>7 h 30 min</td><td>5 h 18 min</td></tr></tbody></table>

<table><caption>**Table 5**
Evolution of computation time as a function of the number of cores – Rolling under imposed load.</caption><thead><tr><th></th><th>6 cores</th><th>12 cores</th><th>20 cores</th></tr></thead><tbody><tr><td>Quasi-static rolling (MEF++)</td><td>83 h</td><td>45 h</td><td>33 h</td></tr></tbody></table>

view to the industrialization of the numerical strategy developed to overcome the problem of tire design and determining the optimal design by proposing intensive computations regarding the possible evolution of certain physical parameters of the tire structure and different operational conditions.

## 5. Conclusions and perspectives

As a reminder, one of the difficulties in simulating the complete dynamic response of a rolling sculpted tire is the multi-scale nature of the dynamic response, which decomposes into a macroscopic rolling dynamic behavior around the rotational frequency, and the vibratory response in a wider frequency window. A key point is that the amplitudes of the macroscopic dynamic behavior are relatively low, while the amplitudes of the vibratory response are relatively high. From this standpoint the first contribution of the proposed study is to describe a strategy to solve the vibratory problem of a patterned tire rolling on a rigid plane. The latter was based on the decomposition of the problem in two main steps with the definition of the configuration of the state of the tire to be perturbed by an excitation representing the physical phenomena of noise generation, linked to successive making and losing contact at the tire/rigid plane contact patch. The major advantage of the proposed strategy is therefore that it enables us to solve the problem in two stages to separate the dynamics operating at different scales, unlike a simple strategy which would involve using a stable energy-conserving time integrator to predict the dynamic response at several scales, and would therefore be less efficient in terms of computation speed, data storage and potentially robust numerical results.

One second major contribution is to promote increased insights into real-world problems and more particularly in the field of engineering systems and tire structures through mathematical and numerical modeling of an industrial rolling sculpted tire and the prediction of the dynamic response at several scales versus different geometric parameters of the sculpted tire. To illustrate this, the formulation and the modeling of a structure with a geometry representing an industrial tire rolling on a rigid plane were described. The numerical simulations highlighted the influence of the rolling conditions on the vibration levels, and the sensitivity of vibrations in relation to several design parameters of the tire tread pattern. All of these results demonstrate not only the relevance of the proposed approach to efficiently and quickly predict the response dynamic at several scales on a real industrial sculpted tire but can also constitute a first step with the aim of optimizing the characteristics of the tire according to different operating conditions and so advanced mathematical modeling and computational tools for engineering leading to future innovations and novel technologies of sculpted tire. Moreover it appears that the proposed numerical results globally follow the trends observed experimentally, even
if it is sometimes difficult to make a direct comparison due to the complexity of the physical mechanisms involved in rolling noise generation. This illustrates that the proposed mathematical modeling and the strategy based on a two-step approach to separate the dynamics occurring at different scales are effective for use as design aid tools when specifying tire tread pattern dimensions and the problems of attenuating non-linear vibrations generated by pavement-road contact and the associated rolling noise.

The second part of this study developed in Section 4 was devoted to the quasi-static rolling simulation with an imposed load. The need to propose an extension of the first mathematical formulation proposed was discussed and validated. Then, the vibratory response was simulated and the differences of behavior compared to the imposed deflection configuration were analyzed. The main contribution of this last part of the study is to propose new scientific developments for mathematical and computational modeling adapted to a complex realistic problem. The new two-step approach proposed in this study then proves to be entirely relevant in the solution of practical engineering problems (i.e. the nonlinear vibrations of patterned tires rolling dynamics subjected to an imposed load).

A large number of perspectives can be considered following this study. Several paths for future developments and research works linked to modeling, computation strategy and numerical resolution are provided in the following:

* First, regarding the modeling of the contact phenomena, we emphasize that a frictionless contact was considered between the tire and the rigid foundation in the framework of this study as the methodology had already been validated for this type of contact [32]. The integration of the frictional contact would not only allow studying the sensitivity of the vibratory response with respect to the grip conditions but also integrating the stick-slip and stick-snap phenomena in the model and thus covering a wider range of structural mechanisms generating rolling noise. Although the efficiency and pertinence of the strategy proposed in this study had already been validated on an academic structure (rolling and frictional contact), a major issue in the framework of an industrial tire system is to efficiently model the vibratory response in the contact patch as a function of the friction regime.

* More specifically, the quasi-static rolling configuration under imposed load charge was employed by assuming a weak coupling between the deflection and the displacement. This formulation is easy to implement but penalizes the rate of convergence. Therefore, a formulation with strong coupling, where the derivative of the additional equation with respect to the displacement is integrated in the tangent matrix, would allow simultaneously calculating the displacement and the deflection in a single Newton loop. Strong coupling would make it possible to obtain the same magnitude of computation time as for the configuration with an imposed deflection. The sensitivity of the vibratory response with respect to the rolling conditions and the design parameters of the tire tread could then be studied with an optimized computation time.

* Most tire vibrations are located close to the entry and exit of the contact patch given the viscous damping properties of tire rubber. The integration of a damping model in the model is crucial to obtain a realistic circumferential distribution of the vibratory response. Thus, one of the priority paths in terms of additional modeling in comparison to the existing model would be the inclusion of a viscous dissipation equation from rheological models (e.g., Kelvin-Voigt or Zener type) [41].

* Tire rolling was controlled by simultaneously imposing the rotation and the translation velocities. They were linked by a first order approximation of the rolling radius. This purely kinematic control is stable and easy to implement. The rolling could also be controlled with imposed torque with an additional equation linking the imposed rotation speed on the rim zone and the torque applied to the wheel. Furthermore, realistic reproduction of the rolling control would consist in imposing only rotation (translation, respectively) by allowing the contact tangential forces to control the advance (rotation, respectively) of the tire. This control could be implemented by imposing the kinematic on a fictional point defined geometrically as the center of the wheel and linked rigidly with the rim zone.

* The impact of the road roughness mechanism could be taken into account in the model by including the profile of the texture in the calculation of the normal gap function. This would make it possible to obtain a more realistic response from the tire structure due to the introduction of an excitation more consistent with the reality of a rolling system in operation. Regarding the problem of the strategy of solving such a problem, the question raised is that of knowing whether the texture should be integrated in the quasi-static rolling or in the vibration problem. The first option would require a fine mesh capable of reproducing the micro-textures of the road so as to capture the entire excitation frequency signature of the vibration problem. The second option would require refining the mesh only for the linear vibration problem whose solution is less costly than that of the nonlinear quasi-static problem. In this case, the problem of the road would be imposed kinematically on the contact patch given by the quasi-static rolling which would be performed on a smooth ground.

## Nomenclature

* $\Omega$ Reference configuration

* $\Omega'$ Current configuration

* $\mathbf{X}$ Coordinates of reference configuration

* $\mathbf{x}$ Coordinates of current configuration

* $\Gamma_u$ Dirichlet boundary condition zone in the reference configuration

* $\Gamma_u'$ Dirichlet boundary condition zone in the current configuration

* $\Gamma_\sigma$ Neumann boundary condition zone in the reference configuration

* $\Gamma_\sigma'$ Neumann boundary condition zone in the current configuration

* $\Gamma_p$ Follower pressure boundary condition zone in the reference configuration

* $\Gamma_p'$ Follower pressure boundary condition zone in the current configuration
$\Gamma_c$ Potential contact zone

$\Gamma_c^t$ Contact active zone

$\varphi$ Mapping of reference to current configuration

$\mathbf{u}$ Displacement field

$\nabla$ Gradient with respect to $\mathbf{X}$

$\mathbf{I}$ Identity second order tensor

$\mathbf{F}$ Deformation gradient

$\mathbf{C}$ Cauchy-Green tensor

$\mathbf{E}$ Green-Lagrange strain tensor

$\boldsymbol{\sigma}$ Cauchy stress tensor

$\mathbf{P}$ First Piola-Kirchhoff stress tensor

$\mathbf{S}$ Second Piola-Kirchhoff stress tensor

$\mathbf{W}$ Strain energy function

$\mathbf{C}$ Fourth order elasticity tensor

$\mathbf{N}$ Outward unit normal in the reference configuration

$\mathbf{n}$ Outward unit normal in the current configuration

$\Pi_n$ Normal contact pressure

$\Pi_t$ Tangential stress vector

$u_n$ Normal displacement

$g_n$ Normal gap function

$V$ Traveling speed

$\omega$ Rotation speed

$r_\omega$ Rolling radius

$d$ Deflection

$\epsilon_n$ Normal penalization

$\rho$ Inflation pressure

$\Delta t$ Time step

$R_g$ Tire geometrical radius

$R_j$ Rim radius

$\omega$ Tire width

($\lambda, \mu$) Lamé constants

$E$ Young Modulus

$\nu$ Poisson coefficient

$\rho$ Density

$H^1(\Omega)$ First order Sobolev space

$u_x$ Longitudinal displacement

$u_y$ Lateral displacement

$u_z$ Vertical displacement

$v_x$ Longitudinal velocity

$v_y$ Lateral velocity

$v_z$ Vertical velocity

$a_x$ Longitudinal acceleration

$a_y$ Lateral acceleration

$a_z$ Vertical acceleration

$e$ Groove width

$\psi_2$ Groove inclination angle

$C_z$ Imposed load

$\bar{d}$ Deflection mean value

$\delta d$ Deflection variation

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The authors do not have permission to share data.
## References

[1] P.R. Donavan, R. Schumacher, Exterior noise of vehicles – traffic noise prediction and control, in: Handbook of Noise and Vibration Control, John Wiley & Sons, Ltd, 2007, pp. 1427–1437, section: 120.

[2] EU, Directive 2001/43/EC of the European Parliament and of the Council of 27 June 2001 amending Council Directive 92/23/EEC relating to tyres for motor vehicles and their trailers and to their fitting, Jun. 2001.

[3] Y. Azizi, Chapter 6 – Generation mechanisms of tire/road noise, in: X. Wang (Ed.), Automotive Tire Noise and Vibrations, Butterworth-Heinemann, 2020, pp. 91–114.

[4] R.E. Hayden, Roadside noise from the interaction of the rolling tire with the road surface, J. Acoust. Soc. Am. 50 (1971) 113.

[5] F. Conte, Modélisation CFD du phénomène acoustique de pompage d'air dans un contact pneumatique/chaussée, Ph.D. thesis, INSA Lyon, Décembre 2008.

[6] U. Sandberg, J. Ejsmont, Tyre/Road Noise Reference Book, 1st edition, INFORMEX Ejsmont and Sandberg Handelsbolag Harg, 2002.

[7] T. Li, A review on physical mechanisms of tire-pavement interaction noise, SAE Int. J. Veh. Dyn. Stab. NVH 3 (1) (2019) 1–26.

[8] W. Kropp, K. Larsson, S. Barrelet, The influence of belt and tread band stiffness on the tire noise generation mechanisms, J. Acoust. Soc. Am. 103 (5) (1998) 2919.

[9] M. Jennewein, M. Bergmann, Investigations concerning tyre/road noise sources and possibilities of noise reduction, Proc. Inst. Mech. Eng., Part D: Transp. Eng. 199 (3) (1985) 199–206.

[10] R.A.G. Graf, C.Y. Kuo, A.P. Dowling, W.R. Graham, On the horn effect of a tyre/road interface, part I: experiment and computation, J. Sound Vib. 256 (3) (2002) 417–431.

[11] T. Li, Influencing parameters on tire-pavement interaction noise: review, experiments, and design considerations, Designs 2 (4) (2018) 38.

[12] P. Le Tallec, Numerical methods for nonlinear three-dimensional elasticity, in: Handbook of Numerical Analysis, vol. 3, Elsevier, 1994, pp. 465–622.

[13] A. Fortin, A. Garon, Les éléments finis: de la théorie à la pratique, GIREF.

[14] P. Wriggers, Computational contact mechanics, Comput. Mech. 3 (2003) 141.

[15] T.A. Laursen, Computational Contact and Impact Mechanics, Springer Berlin Heidelberg, Berlin, Heidelberg, 2003.

[16] O. Gonzalez, Exact energy and momentum conserving algorithms for general models in nonlinear elasticity, Comput. Methods Appl. Mech. Eng. (2000) 1763–1783.

[17] P. Hauret, Méthodes numériques pour la dynamique des structures non-linéaires incompressibles à deux échelles, PhD Thesis. Ecole polytechnique.

[18] U. Nackenhorst, The ale-formulation of bodies in rolling contact: theoretical foundations and finite element approach, Comput. Methods Appl. Mech. Eng. 193 (39) (2004) 4299–4322.

[19] C. Rahier, Analyse numerique du roulage stationnaire de structures viscoelastiques non lineaires en grandes deformations, These de doctorat, Paris 6, Jan. 1991.

[20] C.G. Diaz, et al., Dynamic behaviour of a rolling tyre: experimental and numerical analyses, J. Sound Vib. 364 (2016) 147–164.

[21] D.D. Gregoriisab, F. Naetsbc, P. Kindta, W. Desmetbc, Development and validation of a fully predictive high-fidelity simulation approach for predicting coarse road dynamic tire/road rolling contact forces, J. Sound Vib. 452 (2019) 147–168.

[22] M. Brinkmeier, U. Nackenhorst, H. Volk, A finite element approach to the transient dynamics of rolling tires with emphasis on rolling noise simulation, Tire Sci. Technol. 35 (3) (2007) 165–182.

[23] J.-C. Autrique, F. Magoulès, Studies of an infinite element method for acoustical radiation, Appl. Math. Model. 30 (7) (2006) 641–655.

[24] R. Ogden, Non-linear Elastic Deformations, Ellis Horwood and John Wiley, 1985.

[25] K. Poulos, Y. Renard, An unconstrained integral approximation of large sliding frictional contact between deformable solids, Comput. Struct. 153 (2015) 75–90.

[26] A. Signorini, Questioni di elasticita non linearizzata e semilinearizzata, Rend. Mat. Appl. (1959) 95–139.

[27] J. Nocedal, S.J. Wright, Numerical Optimization, Springer Series in Operations Research, Springer, New York, 1999.

[28] J.C. Simo, T.A. Laursen, An augmented Lagrangian treatment of contact problems involving friction, Comput. Struct. 42 (1) (1992) 97–116.

[29] N. Kikuchi, Y. Song, Penalty/finite-element approximation of a class of unilateral problems in linear elasticity, Q. Appl. Math. 39 (1981) 1–22.

[30] J. Winroth, P. Andersson, W. Kropp, Importance of tread inertia and damping on the tyre/road contact stiffness, J. Sound Vib. 333 (2014) 5378–5385.

[31] I. Knowles, R.J. Renka, Methods for numerical differentiation of noisy data, Electron. J. Differ. Equ. (2012) 12.

[32] Z. Knar, J.-J. Sinou, S. Besset, V. Clauzon, An adapted two-steps approach to simulate nonlinear vibrations of solid undergoing large deformation in contact with rigid plane—application to a grooved cylinder, Appl. Sci. 12 (2022) 1447.

[33] R. Jazar, Vehicle Dynamics: Theory and Application, Springer Science+Business Media, New York, 2017.

[34] Y. Saito, Noise3d: Analytical study of tread pattern design for tire exterior noise, Michelin internal report.

[35] T. Li, J. Feng, R. Burdisso, C. Sandu, Effects of speed on tire-pavement interaction noise (tread-pattern-related noise and non-tread-pattern-related noise), Tire Sci. Technol. 46 (2) (2018) 54–77.

[36] K. Iwao, I. Yamazaki, A study on the mechanism of tire/road noise, JSAE Rev. 17 (2) (1996) 139–144.

[37] G. Tong, Q. Wang, K. Yang, X.C. Wang, An experiment investigation to the radial tire noise, Adv. Mater. Res. 694–697 (2013) 361–365.

[38] J. Hadamard, Lectures on Cauchy's Problem in Linear Partial Differential Equations, Yale University Press, 1953.

[39] V. Yastrebov, Computational contact mechanics: geometry, detection and numerical techniques, PhD Thesis, École Nationale Supérieure des Mines de Paris.

[40] J.E. Dennis, R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Classics in Applied Mathematics, vol. 16, Society for Industrial and Applied Mathematics, Philadelphia, 1996.

[41] P. Le Tallec, C. Rahier, A. Kaiss, Three-dimensional incompressible viscoelasticity in large strains: formulation and numerical approximation, Comput. Methods Appl. Mech. Eng. 109 (3) (1993) 233–258.