# An adapted two-steps approach to simulate nonlinear vibrations of solid undergoing large deformation in contact with rigid plane - Application to a grooved cylinder

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-steps approach to simulate nonlinear vibrations of solid undergoing large deformation in contact with rigid plane - Application to a grooved cylinder. Applied Sciences, 2022, 12 (3), pp.1447. 10.3390/app12031447. hal-03604354

HAL Id: hal-03604354

https://hal.science/hal-03604354

Submitted on 10 Mar 2022

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>

Article

# An Adapted Two-Steps Approach to Simulate Nonlinear Vibrations of Solid Undergoing Large Deformation in Contact with Rigid Plane—Application to a Grooved Cylinder

Zakaria Knar<sup>1,2,†</sup>, Jean-Jacques Sinou<sup>1,3,*,†</sup>, Sébastien Besset<sup>1,†</sup>, Vivien Clauzon<sup>2,†</sup>

**Abstract:** Simulating the full dynamic response of a rolling sculpted tire requires not only taking into account various non-linearities but also considering the multi-scale nature of the dynamic response itself. On one hand, there is the macroscopic rolling dynamic behavior that operates around the rotating frequency with relatively high amplitudes. On the other hand, the vibratory response operates in a larger frequency window with relatively low amplitudes. In contrast to a straightforward strategy that consists of using an energy-conserving stable time integrator to predict the multi-scale dynamic response, the proposed strategy is based on a two-steps approach to separate the dynamics operating at different scales. This methodology is applied to simulate the nonlinear vibrations of a hyperelastic solid undergoing large deformations in contact with a rigid plane. In order to illustrate the potential of the proposed numerical method, the nonlinear vibrations response of a grooved cylinder rolling on a rigid plane is investigated.

**Keywords:** non-linear dynamics; contact mechanics; multi-scale modeling; grooved cylinder; tire

## 1. Introduction

Tire/road noise is becoming a major source of vehicle noise as the transition towards electrical vehicles is underway. This has resulted in tire manufacturers refining their tire design to help reduce traffic noise. The basic physical principles of tire/pavement noise are well understood and documented [1]. For the interested reader, a literature review of design considerations on parameters influencing tire/pavement interaction noise can be found in [2] and a complete review of models on tire-pavement interaction noise is also provided in [3]. The tire/pavement noise generation mechanisms can be divided into two categories. The first one is an air-flow mechanism that generates noise at high frequencies (above 1000 Hz) and is commonly referred to as air-pumping. The second one is the tire vibrations mechanism that generates noise at low and medium frequencies (below 1000 Hz). These noise generation mechanisms are well described in [4]. This paper focuses on the latter mechanism.

The tire solid vibrations are generated by two excitation sources. The first one is the road roughness for which the resulting vibrations could be well predicted for axisymmetric tires using an ALE-formulation [5,6]. This approach consists of using a fixed mesh through which the tire material particles flow. This approach gives a stationary state of the loaded rotating tire. Afterward, the small vibrations due to the road excitation could be computed using a modal superposition technique [7,8]. Most industrial models found in the literature are based on this approach. The second excitation source is the tire's groove resonance. The small scale tread deformation has an impact on the noise emitted by the tire [9]. The main goal of the grooves is to ensure sufficient grip in wet conditions in order

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

### check for updates

**Citation:** Knar, Z.; Sinou, J.-J.; Besset, S.; Clauzon, V. 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.* **2022**, *12*, 1447. https://doi.org/10.3390/app12031447

Academic Editor: Jörg Wallaschek

Received: 15 December 2021

Accepted: 25 January 2022

Published: 28 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, École Centrale de Lyon, 36 Avenue Guy de Collongue, CEDEX, 69134 Écully, France; zakaria.knar@ec-lyon.fr (Z.K.); sebastien.besset@ec-lyon.fr (S.B.)</footnote>

<footnote n="2">² Michelin Campus RDI, Rue Bleue, 63118 Cébazat, France; vivien.clauzon@michelin.com</footnote>

<footnote n="3">³ Institut Universitaire de France, 75231 Paris, France</footnote>

<footnote>* Correspondence: jean-jacques.sinou@ec-lyon.fr</footnote>

<footnote>† The authors contributed equally to this work.</footnote>
to avoid hydroplaning; however, the tire groove geometry parameters (length, width, angle, depth...) and the material proprieties change the vibratory behavior of the tire [10,11]. The prediction of the vibrations induced by this second mechanism can be more complex to achieve and requires the implementation of a specific methodology. Indeed, obtaining a stationary state is no longer possible and full Lagrangian vision should be adopted to treat the tire's large deformations. Hence, a time integration process seems unavoidable. The vibrations small scale requires using energy conserving schemes [12,13] but remains costly from a CPU time perspective.

The main contribution of the proposed study is to promote a numerical methodology based on a two-step approach for solving such problem [14]. Firstly, a non-linear quasi-static simulation is performed to avoid using a high-order stable numerical scheme capable of handling the considered non-linearities. Afterward, the non-stationary equations of motion are linearized using the quasi-static response as a reference point. The method was originally proposed by Valyaev [15] for linear elastic materials. The effectiveness of the proposed approach is illustrated by predicting the non-linear vibrations of a grooved cylinder rolling on a rigid plane with frictionless contact in the first case and with frictional contact in the second case. Even if special attention is given to the mechanical contact, since the tire vibration noise is mainly caused by the impact between the tire and the pavement [4], the main objective is not to propose a reality of the physical phenomena and the interaction between the vibratory response of the tire and the road, but only to demonstrate the capacity of the numerical methodology put forward to correctly predict the vibratory response of solid undergoing large deformation with frictional contact.

This paper is organized as follows. Section 2 presents the mathematical framework used to describe the mechanical behavior of a hyperelastic solid undergoing large deformations in contact with a rigid plane. Afterward, the linearization process is detailed and several hypotheses are made and justified. The efficiency of the proposed strategy is demonstrated by applying it to simulate the non-linear vibrations of a grooved cylinder. The prototype structure and the numerical results are presented in Section 3. Finally, the concluding remarks are given in Section 4.

## 2. Problem Setting

This section sets the mathematical framework used to develop the proposed methodology. Firstly, the mechanical quantities needed to use the large deformation theory are introduced. Afterward, the mathematical formulation of the frictional contact is presented. Then, the equations of motion are given with the corresponding weak formulations. Finally, the perturbation of the quasi-static simulation that gives the vibration response is detailed.

### 2.1. Notations

An open set $\Omega \subset \mathbb{R}^3$ represents the reference state of a continuum body with material points $\mathbf{X} \in \Omega$, $\mathbf{x} = \varphi(\mathbf{X}, t) \in \mathbb{R}^3$ is the deformation field of the body relative to its reference state at time $t \in [0, +\infty[$ and $\mathbf{u}(\mathbf{X}, t) = \mathbf{x} - \mathbf{X}$ is the displacement field. The boundary of the material domain $\partial\Omega$ is composed of $\Gamma_u$ where the displacement is imposed, $\Gamma_f$ where the external force is imposed and $\Gamma_c$ the potential contact zone.

The gradient of a quantity in the reference configuration is noted $\nabla$. To describe the deformation, several mechanical quantities are introduced. The identity matrix is denoted by $\mathbf{I}$ and the deformation gradient is represented by $\mathbf{F} = \mathbf{I} + \nabla\mathbf{u}$. The Jacobian of $\varphi$ is denoted by $J = \det(\mathbf{F})$. $\mathbf{C} = \mathbf{F}^T\mathbf{F}$ represents the Cauchy-Green tensor and $\mathbf{E} = \frac{1}{2}(\mathbf{C} - \mathbf{I})$ is the Green-Lagrange strain tensor. $\sigma$ is the Cauchy stress tensor, $\mathbf{P} = J\sigma\mathbf{F}^{-T}$ is the first Piola-Kirchhoff stress tensor and $\mathbf{S} = J\mathbf{F}^{-1}\sigma\mathbf{F}^{-T}$ 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 through $\mathbf{E}$ [16]. The second Piola-Kirchhoff stress tensor is given by:

$$ \mathbf{S} = \frac{\partial W}{\partial \mathbf{E}}(\mathbf{E}) \qquad (1) $$
with corresponding fourth-order elasticity tensor

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

The fourth-order material elasticity tensor $C$ is necessary within typical finite element procedures and is used in the linearization procedure later on.

## 2.2. Contact Mechanics

A unilateral contact with rigid plane is considered, as illustrated in Figure 1. $\Gamma_c$ is a restricted part of $\partial\Omega$ and represents the surface where contact/friction phenomena may occur. $N$ (resp. $\mathbf{n}$) is the unit outward normal in the reference (resp. current) configuration. The motion of the body is restricted to the upper half-space $z \ge 0$ from the rigid plane $z = 0$.

<figure><img src="image_5.png" /><figcaption>Figure 1. The reference and the deformed configuration of a deformable body in contact with a rigid plane.</figcaption></figure>

To treat the geometrical non linearity, the equilibrium equations are written in the reference configuration. Thus, the second definition of the force will be privileged. The traction vector is decomposed in the deformed configuration as follows:

$$\mathbf{PN} = \Pi_n \mathbf{n} + \Pi_t \qquad (3)$$

where $\Pi_n \mathbf{n}$ is the normal stress vector and $\Pi_t$ is the tangential stress vector at the interface. $\Pi_t$ is calculated as follows:

$$\Pi_t = (\mathbf{I} - \mathbf{n} \otimes \mathbf{n}) \cdot \mathbf{PN} \qquad (4)$$

The contact constraint's mathematical formulation requires also the displacement field decomposition in the contact patch

$$\mathbf{u} = u_n \mathbf{n} + \mathbf{u}_t \qquad (5)$$

where $u_n$ (resp. $\mathbf{u}_t$) is the normal (resp. tangential) displacement.

### 2.2.1. Signorini Conditions

The non-penetration condition in the contact surface is mathematically formulated using the normal gap function $g_n$ defined as:

$$g_n = -\mathbf{x} \cdot \mathbf{n} \ge 0 \qquad (6)$$
The gap function could be written using the displacement **u**

$$ 
\begin{aligned}
g_n(\mathbf{u}) &= -(\mathbf{X} + \mathbf{u}) \cdot \mathbf{n} \\
&= g_0 - u_n
\end{aligned} 
\qquad (7) $$

where $g_0 > 0$ is the initial gap. To prevent the body from penetrating the plane, a non-positive contact pressure arises in the contact zone.

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

Finally, if a point is not in contact, then $\Pi_n = 0$, meanwhile $\Pi_n < 0$ if not. This leads to the non-penetration non-adhesion condition:

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

The force acting on an element of area in the deformed configuration is defined via

$$ df = \sigma \mathbf{n} d\Gamma^t = \mathbf{P} \mathbf{N} d\Gamma \qquad (10) $$

The set of the presented conditions are called Hertz-Signorini-Morea or Karush-Kuhn-Tucker conditions [17]:

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

These conditions provide the basis to treat frictionless contact problems in the context of constraint optimization.

### 2.2.2. Coulomb Friction Law

Friction on the contact interface can be introduced in the general form [18]:

$$ \Pi_t = \Pi_t(\Pi_n, \dot{\Pi}_n, \dot{\mathbf{u}}_t, t, \dots) \qquad (12) $$

The simplest friction law states that the tangential resistance of the contact interface depends on the contact pressure and that the direction of the tangential stress vector is given by the sliding direction **s**.

$$ \Pi_t = \Pi_t(\Pi_n, \mathbf{s}) \quad \text{where} \quad \mathbf{s} = \begin{cases} \frac{\dot{\mathbf{u}}_t}{\|\dot{\mathbf{u}}_t\|} & \text{if } \|\dot{\mathbf{u}}_t\| > 0 \\ \mathbf{0} & \text{if } \|\dot{\mathbf{u}}_t\| = 0 \end{cases} \qquad (13) $$

where $\dot{\mathbf{u}}_t$ is the sliding velocity. The relative motion of the point along the rigid plane is confined by the frictional shear force in the way that the point remains at the initial contact location if the shear stress vector is smaller than the critical frictional stress $\Pi_t^c$, which in Coulomb's friction law is proportional to the contact pressure, so the point does not move if

$$ \|\Pi_t\| \le \Pi_t^c = \mu |\Pi_n| \qquad (14) $$

where the coefficient of proportionality $\mu$ is the coefficient of friction. The maximal permitted tangential stress in equilibrium straight motion is limited by $\mu|\Pi_n|$. The contact condition in the slip regime is

$$ \Pi_t - \mu |\Pi_n| \mathbf{s} = \mathbf{0} \qquad (15) $$

The former condition could be generalized to the stick regime by multiplying it by the norm of the sliding velocity $\|\dot{\mathbf{u}}_t\|$

$$ \|\dot{\mathbf{u}}_t\| \|\Pi_t - \mu |\Pi_n| \dot{\mathbf{u}}_t\| = 0 \qquad (16) $$

A complementary condition is added to distinguish stick and slip states

$$ \|\mathbf{s}\| \|\Pi_t - \mu |\Pi_n| \mathbf{s}\| = 0 \qquad (17) $$
Union of all fore mentioned conditions forms the set of conditions for Coulomb's friction law

$$ ||\mathbf{\Pi}_t|| \leq \mu|\Pi_n|, \quad ||\dot{\mathbf{u}}_t||\mathbf{\Pi}_t - \mu|\Pi_n|\dot{\mathbf{u}}_t = 0, \quad ||\mathbf{s}||\mathbf{\Pi}_t - \mu|\Pi_n|\mathbf{s} = 0 \quad \text{on } \Gamma_c \quad (18) $$

The Coulomb's friction law is sufficient to describe the macroscopic dry friction treated in this paper. For the interested reader, more complex friction laws based on tribological considerations at contact interfaces could be found in [19].

### 2.3. Equations of Motion

The contact conditions together with the relevant boundary conditions complement the local balance of momentum and angular momentum ($\sigma = \sigma^T$). The inertia effects are neglected and a quasi-static process will be considered to deal with Coulomb friction. Hence, the equations of motion are written as follows:

$$ \left\{ \begin{array}{ll} \nabla \cdot \mathbf{P} + \mathbf{f} = \mathbf{0} & \text{in } \Omega \\ \mathbf{P} \mathbf{N} = \mathbf{F}_0 & \text{on } \Gamma_f \\ \mathbf{u} = \mathbf{u}_d & \text{on } \Gamma_u \\ g_n \ge 0, \quad \Pi_n \le 0, \quad g_n \Pi_n = 0 & \text{on } \Gamma_c \\ ||\mathbf{\Pi}_t|| \le \mu|\Pi_n|, \quad ||\dot{\mathbf{u}}_t||\mathbf{\Pi}_t - \mu|\Pi_n|\dot{\mathbf{u}}_t = 0, \quad ||\mathbf{s}||\mathbf{\Pi}_t - \mu|\Pi_n|\mathbf{s} = 0 & \text{on } \Gamma_c \end{array} \right. \quad (19) $$

$\mathbf{f}$ defines the volume external forces, $\mathbf{F}_0$ is a prescribed traction (Neumann boundary conditions) and $\mathbf{u}_d$ is a prescribed displacement (Dirichlet boundary conditions).

In order to obtain the weak formulation, which is the starting point of the finite element method, the first equation of (19) is multiplied by a test function $\mathbf{v}$ and integration by parts is then performed. The weak formulation is stated as follows:

Find $\mathbf{u} \in \mathcal{U} = \{\mathbf{u} \in H^1(\Omega) | \mathbf{u} = \mathbf{u}_d \text{ on } \Gamma_u\}$

$$ \int_{\Omega} \mathbf{P} : \nabla \mathbf{v} \, d\Omega - \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \, d\Omega - \int_{\partial\Omega} \mathbf{P}\mathbf{N} \cdot \mathbf{v} \, d\Gamma = 0 \quad \forall \mathbf{v} \in \mathcal{V} \quad (20) $$

where $\mathcal{V} = \{\mathbf{u} \in H^1(\Omega) | \mathbf{u} = \mathbf{0} \text{ on } \Gamma_u\}$. The integral defined on the boundary are separated as follows:

$$ \int_{\partial\Omega} \mathbf{P}\mathbf{N} \cdot \mathbf{v} \, d\Gamma = \int_{\Gamma_f} \mathbf{F}_0 \cdot \mathbf{v} \, d\Gamma + \int_{\Gamma_c} \mathbf{P}\mathbf{N} \cdot \mathbf{v} \, d\Gamma \quad (21) $$

The first term takes into account the Neumann boundary condition. The second one is related to the contact constraint and could be written as follows

$$ \begin{aligned} \int_{\Gamma_c} \mathbf{P}\mathbf{N} \cdot \mathbf{v} \, d\Gamma &= \int_{\Gamma_c} (\Pi_n \mathbf{n} + \mathbf{\Pi}_t) \cdot (\mathbf{v}_n \mathbf{n} + \mathbf{v}_t) \, d\Gamma \\ &= \int_{\Gamma_c} (\Pi_n \mathbf{v}_n + \mathbf{\Pi}_t \cdot \mathbf{v}_t) \, d\Gamma \end{aligned} \quad (22) $$

The weak formulation could now be written as follows:

$$ \left\{ \begin{array}{ll} \text{Find } \mathbf{u} \in \mathcal{U} & \\ \int_{\Omega} \mathbf{P} : \nabla \mathbf{v} \, d\Omega - \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \, d\Omega - \int_{\Gamma_f} \mathbf{F}_0 \cdot \mathbf{v} \, d\Gamma = \int_{\Gamma_c} \Pi_n \mathbf{v}_n + \mathbf{\Pi}_t \cdot \mathbf{v}_t \, d\Gamma & \forall \mathbf{v} \in \mathcal{V} \quad (23) \\ g_n \ge 0, \quad \Pi_n \le 0, \quad g_n \Pi_n = 0 & \text{on } \Gamma_c \\ ||\mathbf{\Pi}_t|| \le \mu|\Pi_n| \quad ||\dot{\mathbf{u}}_t||\mathbf{\Pi}_t - \mu|\Pi_n|\dot{\mathbf{u}}_t = 0 \quad ||\mathbf{s}||\mathbf{\Pi}_t - \mu|\Pi_n|\mathbf{s} = 0 & \text{on } \Gamma_c \end{array} \right. $$

Several formulations can be applied to incorporate the contact constraints into the variational formulation. Most standard finite element codes that are able to handle contact problems use either the penalty or the Lagrange multiplier method. Each of the methods has its own advantages and disadvantages. The methods are designed to fulfill the constraint equations in the normal direction in the contact interface. For the tangential part, a constitutive relation is needed when stick/slip motion occurs [19].
The penalty method is used to treat the normal contact condition. An active set strategy is employed to obtain the active contact zone. The normal contact condition contribution to the weak formulation is approximated as follows [20].

$$ \int_{\Gamma_c} \Pi_n v_n \, d\Gamma = \int_{\Gamma_c} \varepsilon_n g_n^{-}(\mathbf{u}) v_n \, d\Gamma \qquad (24) $$

where $\varepsilon_n$ is the normal penalty parameter and $g_n^{-} = \min(g_n, 0)$ is the penetration function. As $\varepsilon_n \to \infty$, the gap function $g_n \to 0$ on $\Gamma_c$, the constraint is increasingly well-satisfied. However, large values of $\varepsilon_n$ will lead to an ill-conditioned numerical problem.

The Coulomb friction law is regularized using a linear approximation in the stick regime. The tangential contact contribution to the weak formulation is written as in [21]:

$$ \int_{\Gamma_c} \Pi_t \cdot \mathbf{v}_t \, d\Gamma = \int_{\Gamma_c} \mu |\Pi_n| \begin{cases} \frac{1}{\varepsilon_r} (2 - \frac{\|\dot{\mathbf{u}}_t\|}{\varepsilon_r}) \dot{\mathbf{u}}_t \cdot \mathbf{v}_t \, d\Gamma & \text{if } \|\dot{\mathbf{u}}_t\| < \varepsilon_r \\ \frac{1}{\|\dot{\mathbf{u}}_t\|} \dot{\mathbf{u}}_t \cdot \mathbf{v}_t \, d\Gamma & \text{if } \|\dot{\mathbf{u}}_t\| \ge \varepsilon_r \end{cases} \qquad (25) $$

where $\varepsilon_r$ is the regularization coefficient. It determines the boundary between the stick and the slip regime. If $\varepsilon_r \to 0$, we obtain the classical Coulomb law. Since the problem is written in a quasi-static framework, the sliding velocity $\dot{\mathbf{u}}_t$ is approximated using a finite difference method. The analytical contribution of each integral term to the tangent matrix and to the nonlinear residual could be found in [21].

## 2.4. Dynamic Response

The weak formulation is solved using finite element method. Afterward, the quasi-static response is used to linearize the dynamic response. Several hypotheses are made and explained and the linearization process is detailed.

### 2.4.1. Main Hypothesis

The quasi-static problem, as formulated in the former section, is given by:

$$ \int_{\Omega} \mathbf{P}_1 : \nabla \mathbf{v} \, d\Omega - \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \, d\Omega - \int_{\Gamma_f} \mathbf{F}_0 \cdot \mathbf{v} \, d\Gamma = \mathbf{R}_c(\mathbf{u}_1) \quad \forall \mathbf{v} \in \mathcal{V} \qquad (26) $$

where $\mathbf{R}_c$ is the virtual work of the contact forces. The virtual work due to the inertial forces is added to the weak formulation. The dynamic problem is given by:

$$ \int_{\Omega} \rho \ddot{\mathbf{u}}_2 \cdot \mathbf{v} \, d\Omega + \int_{\Omega} \mathbf{P}_2 : \nabla \mathbf{v} \, d\Omega - \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \, d\Omega - \int_{\Gamma_f} \mathbf{F}_0 \cdot \mathbf{v} \, d\Gamma = \mathbf{R}_c(\mathbf{u}_2, \dot{\mathbf{u}}_2) \quad \forall \mathbf{v} \in \mathcal{V} \qquad (27) $$

where $\rho$ is the density. The dynamic response $\mathbf{u}_2$ is supposedly the sum of the quasi-static response $\mathbf{u}_1$ and the vibrations field $\mathbf{w}$. In other words, the inertia effects are only considered for the vibrations problem.

$$ \mathbf{u}_2(\mathbf{X}, t) = \mathbf{u}_1(\mathbf{X}, t) + \mathbf{w}(\mathbf{X}, t) + o(\mathbf{w}) \qquad (28) $$

An homogeneous Dirichlet boundary condition is applied on $\Gamma_u$ since the prescribed displacement is the same for both problems.

$$ \mathbf{w} = \mathbf{0} \quad \text{on } \Gamma_u \qquad (29) $$

The active contact zone $\Gamma_c^t$ and the virtual work of the contact forces are supposedly the same for the dynamic and quasi-static problem.

$$ \mathbf{R}_c(\mathbf{u}_2, \dot{\mathbf{u}}_2) \approx \mathbf{R}_c(\mathbf{u}_1) \qquad (30) $$

The active contact zone is not exactly the same since the inertia terms will change the contact stiffness and thus alter the contact patch [22]. Another difference between the
virtual works is the friction law formulation. The sliding velocity is approximated in the quasi-steady regime whereas it is equal to the tangent velocity in the dynamic regime.

Given the hypothesis that the inertia effect does not change the contact patch, the normal displacement in the quasi-steady regime does not change in the dynamic regime since it is determined by the rigid plane position. Therefore, the normal vibrations field is null on the contact active zone.

$$ \mathbf{w} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_c^t \qquad (31) $$

This boundary condition is an interpretation of the contact constraint as a Dirichlet-like boundary condition. The active contact zone $\Gamma_c^t$ is not related to a material domain, which is problematic from a mathematical point view; however, it is easily implemented after a finite element discretization.

### 2.4.2. Linearization

The linearization process starts with expressing the mechanical quantities of the dynamic problem as function of the mechanical quantities of the quasi-static problem with neglecting high order terms

$$ \begin{aligned} \mathbf{F}_2 &= \mathbf{F}_1 + \nabla \mathbf{w} \\ \mathbf{E}_2 &= \mathbf{E}_1 + \frac{1}{2}(\mathbf{F}_1^\top \nabla \mathbf{w} + \nabla^\top \mathbf{w} \mathbf{F}_1 + \underbrace{\nabla^\top \mathbf{w} \nabla \mathbf{w}}_{\approx 0}) \\ &\approx \mathbf{E}_1 + \frac{1}{2}(\nabla \mathbf{w} + \nabla^\top \mathbf{w} + \nabla^\top \mathbf{u}_1 \nabla \mathbf{w} + \nabla^\top \mathbf{w} \nabla \mathbf{u}_1) \\ &\approx \mathbf{E}_1 + \mathbf{E}_{\mathbf{w}, \mathbf{u}_1} \end{aligned} \qquad (32) $$

The constitutive material law is linearized as follows

$$ \begin{aligned} \mathbf{S}_2 &= \frac{\partial W}{\partial \mathbf{E}}(\mathbf{E}_2) = \frac{\partial W}{\partial \mathbf{E}}(\mathbf{E}_1 + \mathbf{E}_{\mathbf{w},\mathbf{u}_1}) \\ &\approx \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{w},\mathbf{u}_1} \\ &\approx \mathbf{S}_1 + C_1 : \mathbf{E}_{\mathbf{w},\mathbf{u}_1} \end{aligned} \qquad (33) $$

The deformation integral term is linearized as follows

$$ \begin{aligned} \mathbf{P}_2 : \nabla \mathbf{v} &= \mathbf{S}_2 : (\mathbf{F}_2^\top \nabla \mathbf{v}) \\ &= (\mathbf{S}_1 + C_1 : \mathbf{E}_{\mathbf{w},\mathbf{u}_1}) : (\mathbf{F}_1^\top \nabla \mathbf{v} + \nabla^\top \mathbf{w} \nabla \mathbf{v}) \\ &= \mathbf{P}_1 : \nabla \mathbf{v} + \mathbf{S}_1 : \nabla^\top \mathbf{w} \nabla \mathbf{v} + (C_1 : \mathbf{E}_{\mathbf{w},\mathbf{u}_1}) : (\mathbf{F}_1^\top \nabla \mathbf{v}) \end{aligned} \qquad (34) $$

The weak formulation governing the vibrations response is given by

$$ \int_{\Omega} \rho \ddot{\mathbf{w}} \cdot \mathbf{v} \, d\Omega + \int_{\Omega} \mathbf{S}_1 : \nabla^\top \mathbf{w} \nabla \mathbf{v} + (C_1 : \mathbf{E}_{\mathbf{w},\mathbf{u}_1}) : (\mathbf{F}_1^\top \nabla \mathbf{v}) \, d\Omega = - \int_{\Omega} \rho \ddot{\mathbf{u}}_1 \cdot \mathbf{v} \, d\Omega \quad \forall \mathbf{v} \in \mathcal{V} \qquad (35) $$

The deformation integral term depends on the gradient of the quasi-static displacement. The quasi-steady acceleration $\ddot{\mathbf{u}}_1$ appears as the excitation force. It is approximated using a finite difference scheme. The second derivative numerical approximation generates high frequency noise. It could be handled by using low-pass filter or by choosing an adequate discretization. Equation (35) is linear with respect to $\mathbf{w}$ and could be solved using a time integration scheme.

## 3. Numerical Results

In this section, a numerical example intended to demonstrate the utility and performance of the proposed method is presented. The quasi-static problem is solved using a finite element code called MEF++, which is developed by the GIREF of Laval University (Groupe Interdisciplinaire de Recherche en Éléments Finis). The studied prototype is presented and the simulation numerical parameters are given. Afterward, the quasi-static
simulation results are analyzed. Finally, the vibrations field is extracted using the presented method. The coefficient of friction impact on the quasi-static problem and on the vibrations is investigated.

### 3.1. Prototype

The method presented in the former section is applied to predict the dynamic response of the grooved cylinder shown in Figure 2. The cylinder is composed of 10 blocks. The chosen mesh is a result of a convergence study with respect to the quasi-static displacement and the contact forces. The mesh is composed of 35,640 nodes and 29,440 linear elements.

<figure><img src="image_6.png" /><figcaption>Figure 2. The prototype mesh.</figcaption></figure>

The material domain is composed of the potential contact surface (in red) $\Gamma_c$ and the rim zone $\Gamma_u$ (in blue) where the following Dirichlet boundary condition is applied:

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

The structure is loaded kinematically by applying a vertical displacement $d$ in the rim zone. The load could be thus calculated by summing the vertical reactions in the contact patch. $V$ is the translation velocity, and $\mathbf{R}$ is the rotation matrix associated with the angular velocity $\omega$. The translation velocity depends on the angular velocity via the effective radius (also known as the rolling radius) $r_\omega = V/\omega$. A first order approximation of $r_\omega$ is given as a function of the geometrical radius $R_e$ (see [23] for more details) as follows

$$ r_{\omega} \approx R_{e} - \frac{d}{3} \quad (37) $$

The tire rubber is commonly modeled by using an hyperelastic material constitutive law. The simplest hyperelastic material model is the Saint Venant–Kirchhoff model, which is just an extension of the geometrically linear elastic material model to the geometrically nonlinear regime.

$$ \mathbf{S} = \lambda \mathrm{tr}(\mathbf{E})\mathbf{I} + 2\mu\mathbf{E} \quad (38) $$

where $\lambda$ and $\mu$ are the Lamé constants. The numerical results presented in the next section are obtained using the geometrical and the material constitutive law properties given in Table 1. The rolling conditions and the simulation numerical parameters are also given in Table 2. The chosen parameters are the result of a preliminary convergence study that is not presented. In fact, the resulting quasi-static displacement is practically invariant with respect to higher values of $\epsilon_n$ and lower values of $\epsilon_r$ using the chosen mesh and the mechanical properties given in Table 1.
<table><caption>**Table 1.** Prototype mechanical properties.</caption><tr><td>External radius</td><td>R<sub>e</sub> = 0.1 m</td></tr><tr><td>Internal radius</td><td>R<sub>i</sub> = 0.02 m</td></tr><tr><td>Width</td><td>w = 0.1 m</td></tr><tr><td>Young Modulus</td><td>E = 2 × 10<sup>8</sup> Pa</td></tr><tr><td>Poisson coefficient</td><td>ν = 0.3</td></tr><tr><td>Density</td><td>ρ = 900 kg/m<sup>3</sup></td></tr></table>

<table><caption>**Table 2.** Simulation parameters.</caption><tr><td>Deflection</td><td>d = 0.01 m</td></tr><tr><td>Rotation speed</td><td>ω = 150 rad/s</td></tr><tr><td>Normal penalization</td><td>ε<sub>n</sub> = 10<sup>12</sup></td></tr><tr><td>Regularization coefficient</td><td>ε<sub>r</sub> = 10<sup>-6</sup></td></tr><tr><td>Time step</td><td>Δt = 10<sup>-4</sup> s</td></tr></table>

## 3.2. Quasi-Static Rolling

The quasi-static rolling results are now presented. In order to study the deformation of the structure, the displacement field is extracted for different points as shown in Figure 3.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Nodes selection: (left) $r = R_e (\bullet)$, $r = (R_e + R_i)/2 (\bullet)$, $r = R_i (\bullet)$. (right) $y = width/2 (\bullet)$, $y = width/4 (\bullet)$, $y = 0 (\bullet)$.</figcaption></figure>

The first selection provides insight into the displacement variation in the radial direction while the second one provides an insight into the displacement variation in the lateral direction.

### 3.2.1. Preample-Frictionless Contact

The analysis of the quasi-static rolling is initiated by considering a frictionless contact. Although this first case studied, based on frictionless contact, does not make sense for the reality of tire-road interaction, as it does not correspond to the real phenomena of tire contact, this preliminary study is carried out only to illustrate the feasibility and validity of the proposed two-steps approach on the grooved cylinder under study. In addition, this first application with frictionless contact will serve as a reference for a qualitative comparison with the following case, which considers the impact of friction. So the main objective of this section is to study the kinematics and the contact pressure of the grooved cylinder in order to validate both the mechanical modeling and the associated numerical methodology.

Figure 4 shows the time evolution of the quasi-static displacement components. The displacement field is, as expected, periodic. The longitudinal displacement $U_x$ slows down at the entry of the contact patch ($t \approx 0.05$ s) and speeds up at the exit of the contact patch ($t \approx 0.055$ s) to reach the rim longitudinal position. The longitudinal deformation decreases as we move away from the contact zone. The same kinematical behavior is
observed for tires in [23]. In fact, the radial distance between the tread pattern and the rotation axis is minimal at the contact patch and reach its maximum at the opposite of the contact patch. The vertical displacement $U_z$ shows that the non-penetration condition is well respected ($0.05 \text{ s} \le t \le 0.055 \text{ s}$) on the contact patch. The lateral displacement amplitude $U_y$ decreased as we move away from the contact patch due to the Poisson effect. The lateral deformation is obviously greater in the contact active zone and is null at the opposite of the contact patch ($t \approx 0.07 \text{ s}$). Furthermore, the lateral displacement has a symmetric distribution; the deformation at the entry of the contact zone is the same at the exit of the contact patch and the maximum is reached at the contact patch center ($t \approx 0.052 \text{ s}$). The displacement of the node occupying the rim zone shows that the Dirichlet boundary conditions are well respected.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Time evolution of the quasi-static displacement for $r = R_e$ (—), $r = (R_e + R_i)/2$ (—) and $r = R_i$ (—) with $y = \text{width}/2$.</figcaption></figure>

Figure 5 shows that the lateral position does not impact the longitudinal and vertical displacement. The lateral displacement amplitude decreases toward 0 as we move to the contact patch center. This result was expected since the plane ($y = 0$) is a symmetry plane of the structure.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Time evolution of the quasi-static displacement for $y = \text{width}/2$ (—), $y = \text{width}/4$ (—) and $y = 0$ (—) with $r = R_e$.</figcaption></figure>

Equation (35) shows that the quasi-steady acceleration appears as the excitation force of the vibrations response. The acceleration evolution in the circumferential direction and in the frequency domain is displayed in order to have an idea on how the vibrations will be excited.

Figure 6 shows the time evolution of the acceleration in the rotating frame. The angle $\theta = 0$ deg denotes the contact patch center and $\theta < 0$ deg designates the contact patch exit. From a time perspective, the acceleration's maximum is reached when the block edge enters ($t \approx 0.052 \text{ s}$) and leaves ($t \approx 0.0563 \text{ s}$) the contact patch. From a circumferential point of view, the maximum of the vertical and lateral acceleration is located at the contact patch boundaries and have a similar distribution at the entry ($\theta \approx 20$ deg) and at the exit ($\theta \approx -17$ deg) of the contact active zone. The vertical acceleration is null in the
contact zone since the vertical displacement is constant in the contact patch. The lateral acceleration decreases towards 0 as we move away from the contact patch. The longitudinal acceleration has a shear distribution in the contact patch. An acceleration at the entry of the contact patch and a deceleration at the exit patch are observed. The balance between the acceleration and the deceleration varies with respect to the contact patch form.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Time evolution of the quasi-steady acceleration in the circumferential direction ($y = width/2$ and $r = R_e$).</figcaption></figure>

Figure 7 shows the associated quasi-steady acceleration spectrum plotted with respect to the order. The order is defined as the frequency over the rotation speed and it is used to analyze the spectrum distribution independently of the rotation frequency. The fundamental frequency is equal to the rotation frequency (Order = 1). Two families of harmonics appear since various non-linearities are taken into account in the model. The first family of harmonics operates around the multiples of the rotation frequency $n\omega$, $n \in N$. The second family of harmonics operates around $10\omega$, $k \in N$, with a modulation effect. It concerns the contact forces frequency since the structure is composed of 10 blocks. This modulation phenomena is due to the rotation effect [24,25]. In fact, the kinematical quantities are calculated with respect to the fixed frame while the contact forces are connected to the rotating frame. This phenomena is not observed in the lateral direction since it represents the rotation axis direction. The vertical and longitudinal accelerations have the same order of magnitude, which is greater than the one of the lateral acceleration. Finally, the acceleration amplitude decreases as we move away from the contact patch.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Frequency spectrum of the quasi-steady acceleration for $r = R_e$ (-), $r = (R_e + R_i)/2$ (-) and $r = R_i$ (-) with $y = width/2$.</figcaption></figure>

Figure 8 shows the contact pressure distribution with respect to the position of the groove (leading and trailing edge). When a block enters the contact zone, the contact
pressure is concentrated at the entry. As the block advances, the contact pressure moves to the block center and becomes similar to the distribution of the cylindrical contact pressure described in Hertz theory [26]. As the block exits the contact patch, the pressure becomes concentrated in the trailing edge with the same distribution as the one observed for the leading edge. A small variation of the total contact force is observed. Since the deflection *d* is constant and the contact patch form is changing, the load supported by the structure varies. The last figure shows the spectrum of the total force variation *δF*<sub>z</sub>, which represents the difference between the total contact force and its mean value. It shows that the variation spectrum is composed of harmonics located around the multiples of the blocks frequency. In the following, the Coulomb friction law is now considered and its impact on the quasi-static rolling is investigated.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Time evolution of the normal contact pressure (top), the total contact force (middle) and the frequency spectrum of the contact force variation (bottom). Node selection for the contact pressure analysis: block entry θ<sub>e</sub>(—), block center θ<sub>c</sub>(—), block exit θ<sub>s</sub>(—).</figcaption></figure>

3.2.2. Extension by Considering the Impact of the Friction Coefficient

The Coulomb friction is now taken into account using the regularization method pre-
sented in Section 2. The friction coefficient impact on the quasi-static rolling is investigated.

Figure 9 shows the displacement evolution in the contact zone for various values of μ. We observe that increasing μ favors the stick regime. The longitudinal displacement remains constant for μ ≥ 0.3 and the lateral displacement variation decreases as the coefficient of friction increases. Indeed, the tangential pressure needed to switch to the slip regime increases knowing that the structure is subjected to the same deflection. The lateral displacement maximum has moved towards the exit of the contact patch when the friction is added. Furthermore, the lateral displacement at the entry of the contact patch is no longer similar to the one at the exit of the contact patch. The vertical displacement is more or less the same since the friction is mainly acting on the tangential components.

Figure 10 shows the normal contact pressure distribution at different stages of the contact. Similar to the frictionless contact, the contact pressure maximum is located at the entrance (resp. exit) for the nodes located at the entry (resp. exit) of the block. Unlike the frictionless case, the distribution at the exit is different than the distribution at the entry since the friction law is characterized by a time-dependent path and the mechanical system is no longer conservative. Increasing μ leads also to the saturation of the normal contact pressure.
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Time evolution of the quasi-static displacement in the contact patch ($y$ = width/2 and $r = R_e$) for frictionless contact (—), $\mu = 0.1$ (—), $\mu = 0.3$ (—), $\mu = 0.5$ (—) and $\mu = 1.0$ (—).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Time evolution of the normal contact pressure in the contact patch ($y$ = width/2, $r = R_e$) for frictionless contact (—), $\mu = 0.1$ (—), $\mu = 0.3$ (—), $\mu = 0.5$ (—) and $\mu = 1.0$ (—) ($\theta = \theta_e$ (left); $\theta = \theta_c$ (middle); $\theta = \theta_s$ (right)).</figcaption></figure>

Figure 11 shows the tangential contact pressure is increasing with respect to $\mu$ unlike the normal contact pressure where a saturation effect is noticed. The tangential contact pressure has the same evolution as the normal pressure for the nodes located at the middle of the block $\theta = \theta_c$ and at the exit of the block $\theta = \theta_s$; however, an increase in the norm of the tangential pressure of the leading edge at the exit of the contact patch is observed. The differences observed between the tangential pressure for various values of $\mu$ is explained by the change of the slope in the stick regime and the change of the tangential pressure value in the slip regime.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Time evolution of the tangential contact pressure norm in the contact patch ($y$ = width/2 and $r = R_e$) for frictionless contact (—), $\mu = 0.1$ (—), $\mu = 0.3$ (—), $\mu = 0.5$ (—) and $\mu = 1.0$ (—) ($\theta = \theta_e$ (left); $\theta = \theta_c$ (middle); $\theta = \theta_s$ (right)).</figcaption></figure>

Figure 12 shows the total contact force variation. As for the frictionless case, the total lateral force $F_y$ is null since the plane ($y = 0$) is a symmetry plane of the structure. The variation of the vertical load is a little bit greater as the coefficient of friction $\mu$ increases until reaching an invariant evolution. In fact, $F_z$ remains the same for $\mu = (0.3 – 0.5 – 1.0)$.
This is coherent with the saturation phenomena observed for the normal contact pressure. The total longitudinal contact force amplitude decreases as we increase $\mu$. High frequency oscillations appear for larger values of $\mu$. These oscillations represent a consequence of the tangent pressure variation described in Figure 11.

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Time evolution of the total contact force ($y$ = width/2 and $r = R_e$) for frictionless contact (—), $\mu = 0.1$ (—), $\mu = 0.3$ (—), $\mu = 0.5$ (—) and $\mu = 1.0$ (—).</figcaption></figure>

### 3.3. Vibration Response

The linearization method is used to calculate the vibrations field. The vibrations problem is solved using a time integration technique via an implicit Euler scheme. Since the prototype has a relatively simple geometry and constitutive material law, the dynamic problem could also be solved using the same integration scheme. The validity of the linearization process is then verified. Afterward, the coefficient of friction influence on the vibratory response is investigated.

#### 3.3.1. Preample-Frictionless Contact

As for the quasi-static rolling, a frictionless contact is considered first. As a reminder, this first preliminary study based on the prediction of the nonlinear vibration response with frictionless contact only aims to illustrate the feasibility of the proposed two-steps approach. A comparison between the linearization method solution and the reference solution, which represents the difference between the dynamic and the quasi-static response is presented. Then, the vibrations velocity spectrum is analyzed. Finally, the vibrations spatial evolution is visualized to assess the global vibration behavior of the grooved cylinder.

Figure 13 illustrates the comparison between the linearized solution and the reference solution. Both solutions are, as expected, periodic. Both solutions have the same vibration frequency with a small difference with respect to the phase. The vertical vibrations are null in the contact patch for both methods; therefore, treating the normal contact constraint as a Dirichlet-like boundary condition seems to be the right choice. The longitudinal vibrations maximum amplitude is located in the contact patch and have the same order of magnitude as the vertical vibrations. The lateral vibrations are located in the contact zone and decrease towards 0 as we move away from the contact zone. The magnitude of the lateral vibrations are much lower than the vertical and longitudinal components.
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Time evolution of the vibrations displacement for ($y = \text{width}/2, r = R_e$) of the linearized solution (—) and the reference solution (—) (full simulation (top) and zoom around one period (bottom)).</figcaption></figure>

Figure 14 shows the vibrations velocity spectrum. Both solutions give more or less the same vertical and longitudinal velocities. The main difference between the methods are observed in the lateral direction. A set harmonics around the multiples of the blocks frequency with the same modulation effect described in the quasi-static rolling results is observed for both solutions. The lateral vibrations velocity spectrum has a frequency comb distribution. The harmonics around the multiples of the rotation frequency have vanished for the longitudinal and vertical vibrations. Hence, the linearization method allows us to separate the dynamic behaviors operating at different scales. The vertical and longitudinal vibration velocity magnitude is more or less the same and greater than the magnitude of the lateral velocity. The same observation was found for the quasi-steady acceleration.

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Frequency spectrum of the vibration velocity for ($y = \text{width}/2, r = R_e$) of the linearized solution (—) and the reference solution (—).</figcaption></figure>

Figure 15 illustrates the vibrations variation in the radial and lateral directions. The vibrations amplitude decreases as we move away from the tread pattern in the radial direction. The lateral vibrations amplitude decreases as we move towards the contact center in the lateral direction whereas the vertical and longitudinal vibrations remain constant. The same observation was found for the quasi-static displacement in Figures 4 and 5.
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Time evolution of the vibration displacement of the linearized solution for $r = R_e$ (—), $r = (R_e + R_i)/2$ (—) and $r = R_i$ (—) $-y = \text{width}/2$ (top) and for $y = \text{width}/2$ (—), $y = \text{width}/4$ (—) and $y = 0$ (—) $-r = R_e$ (bottom).</figcaption></figure>

Figure 16 shows the vibrations velocity variation in the circumferential direction. Unlike the quasi-steady acceleration, the mechanical behavior at the entry of the contact patch is no longer similar to the behavior at the exit of the contact patch. The vertical vibrations are, as expected, null in the contact patch and at the opposite of the contact patch ($\theta = 180$ deg). The lateral vibrations maximum is situated at the contact patch entry and decreases towards 0 as we move away from the contact active zone. The longitudinal vibrations velocity has a more or less symmetric distribution with respect to the plane ($z = R_i$).

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Time evolution of the vibration velocity in the circumferential direction ($r = R_e$, $y = \text{width}/2$) using the linearization method.</figcaption></figure>

### 3.3.2. Coefficient of Friction Impact

The aim of this section is to study the coefficient of friction impact on the linearization process and on the vibrations response. A comparison between the displacement components of the linearized solution with the reference solution for various values of $\mu$ is displayed. A new boundary condition is proposed to improve the linearization method. Finally, the coefficient of friction impact on the vibrations amplitude is investigated.

Unlike the frictionless case, a greater difference between the vibrations is observed in Figure 17 for all values of $\mu$. The differences are especially observed in the contact active zone for the longitudinal and lateral vibrations displacement. Increasing the value of $\mu$ makes the vibrations constant in the contact patch for the reference solution. Indeed, the quasi-static rolling results show that increasing $\mu$ favors the stick regime. Hence, the longitudinal and lateral vibrations displacements remain constant in the contact patch. The constant value is not equal to 0 since the inertia changes the contact stiffness [22]. The boundary condition (31) used to treat the contact constraint is no longer valid for
friction. A simple change of this condition is to impose an homogeneous displacement in the contact active zone as follows

$$w = 0 \quad \text{on} \quad \Gamma_c^t \qquad (39)$$

The results of the boundary conditions change are given in Figure 18 that illustrates the time evolution of the vibrations field for various values of $\mu$. A clear improvement of the linearized solution with respect to the reference solution is observed whether inside or outside the contact patch and especially for high values of $\mu$. The biggest difference is observed for $\mu = 0.1$. Indeed, for this value of $\mu$, the nodes in the contact patch operate in the slip regime, which is clearly incompatible with the new boundary condition but still better than the previous one especially outside the contact zone. A more precise boundary condition to treat the frictional contact should depend on the contact regime. Now, we compare between the vibrations for different values of $\mu$ for both solutions.

Figure 19 shows the vibrations sensitivity with respect to the coefficient of friction. Both solutions give more or less the same observations expect for $\mu = 0.1$ due to the incompatibility between the contact boundary condition and the friction slip regime. Larger vibrations are observed for the frictionless case for both solutions. Increasing $\mu$ results in damping the vibrations amplitude until reaching an invariant vibration response for $\mu > 0.3$. This translates the saturation effect observed for the quasi-static rolling. The same saturation phenomenon is observed in [27] for advancing pin-loaded contacts with friction where the equations of motion are solved using the Green's functions.

<figure><img src="image_21.png" /><figcaption>Figure 17. Time evolution of the vibrations displacement of the linearized solution (—) and the reference solution (—) for various values of $\mu$ ($y = width/2, r = R_e$).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 18.** Time evolution of the vibrations displacement of the linearized solution (—) with the new boundary condition (—) and the reference solution (—) for various values of $\mu$ ($y = \text{width}/2, r = R_e$).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 19.** Time evolution of the vibrations displacement of the linearized solution with the new boundary condition (**right**) (—) and the reference solution (**left**) for frictionless contact (—), $\mu = 0.1$ (—), $\mu = 0.3$ (—), $\mu = 0.5$ (—) and $\mu = 1.0$ ($y = \text{width}/2, r = R_e$).</figcaption></figure>

4. Conclusions

In this paper, an adapted two-step approach to simulate the nonlinear vibrations of a hyperelastic solid undergoing large deformations in contact with rigid plane is presented. The methodology allows not only the separation between the dynamics operating at
different scales but also shows numerical robustness and consistency. The Coulomb friction law is used in the contact interface. This numerical approach is applied to a finite element model corresponding to a grooved cylinder. The linearization process seems to work very well for frictionless contact; however, significant differences can be observed for the case with frictional contact, especially for higher values of $\mu$. The change of the boundary condition in the contact active zone seems to improve the linearized solution. A much more precise condition in the contact active zone would probably depend on the friction (slip or stick) regime.

As previously mentioned, the example proposed and discussed in this paper corresponds to a grooved cylinder undergoing large deformation in contact with rigid plane. One of the most crucial next steps is to be able to provide an illustration of the feasibility of the proposed approach on a tire structure. Numerous future research can be carried out to increase the understanding of non-linear phenomena on tire structures. Although it is not possible to give an exhaustive list of topics of interest, some studies that seem to have a certain priority and interest may be considered:

* First of all, it can be noted that damping has been neglected in modeling; however, experiments show that vibration energy is not transported to the opposite of contact patch due to strong damping of tire structure for such relatively high frequencies. The vibration energy is thus concentrated close to the entrance and the exit of contact patch. Moreover the vibrations evolution in the circumferential direction shows the opposite observation; therefore, the addition of a suitable damping model seems important and is one of the priority issues to be explored. The inclusion of the Kelvin-Voigt model, often used to treat the viscosity of tire rubber, seems to be an interesting option [28].

* In a similar way, the modeling of friction at the contact interface should be completed by the consideration of the coefficient of rolling friction [29] in order to reproduce more realistic physical phenomena involved in tire/pavement interaction vibration and noise. Indeed, the consideration of only the coefficient of sliding friction force as proposed in this study may be not sufficient to accurately reproduce real phenomena of tire vibration response.

* It should be interesting to apply the proposed two-step approach on a real industrial tire structure and to compare the results with full experimental tests. One of the most important challenges may be to identify the different sources of nonlinear vibrations in order to identify and quantify the first-order physical phenomena involved. In this context, more realistic modeling including both tire rubber viscosity and Hertzian contact for polymorph materials or the model of contact 3D multi-layered structures with orthotropic material should be developed to reproduce more specifically the design of a real wheel tire.

* Extension of the prediction of non-linear vibrations in a real environment and in real conditions would also be interesting with the consideration of a more realistic pavement with imperfections. The road roughness excitation mechanism could be added to the proposed methodology by integrating the road texture in the normal gap function calculation. In this context, one of the main challenges for engineers working on the developments of optimization procedures for tire structures in a real environment should be the inclusions of modeling errors and uncertainties.

* The developments and combination of optimization procedures, genetic algorithms or artificial neural in relation with the proposed methodology should also be of prior interest to conduct parametric studies leading to increase the reliability and safety of tire structures subjected to complex and various external excitations.

**Author Contributions:** Conceptualization, Z.K., J.-J.S., S.B. and V.C.; Investigation, Z.K., J.-J.S., S.B. and V.C.; Methodology, Z.K., J.-J.S., S.B. and V.C.; Validation, J.-J.S., S.B. and V.C.; Visualization, Z.K.; Writing—original draft, Z.K.; Writing—review and editing, J.-J.S., S.B. and V.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.
**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** J.-J. Sinou acknowledges the support of the Institut Universitaire de France.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Nomenclature**

<table><tr><td>Ω</td><td>Reference configuration</td></tr><tr><td>Ω<sup>t</sup></td><td>Current configuration</td></tr><tr><td>X</td><td>Coordinates of reference configuration</td></tr><tr><td>x</td><td>Coordinates of current configuration</td></tr><tr><td>Γ<sub>u</sub></td><td>Dirichlet boundary condition zone in the reference configuration</td></tr><tr><td>Γ<sub>u</sub><sup>t</sup></td><td>Dirichlet boundary condition zone in the current configuration</td></tr><tr><td>Γ<sub>f</sub></td><td>Neumann boundary condition zone in the reference configuration</td></tr><tr><td>Γ<sub>f</sub><sup>t</sup></td><td>Neumann boundary condition zone in the current configuration</td></tr><tr><td>Γ<sub>c</sub></td><td>Potential contact zone</td></tr><tr><td>Γ<sub>c</sub><sup>t</sup></td><td>Contact active zone</td></tr><tr><td>φ</td><td>Mapping of reference to current configuration</td></tr><tr><td>u</td><td>Displacement field</td></tr><tr><td>∇</td><td>Gradient with respect to X</td></tr><tr><td>I</td><td>Identity second order tensor</td></tr><tr><td>F</td><td>Deformation gradient</td></tr><tr><td>C</td><td>Cauchy-Green tensor</td></tr><tr><td>E</td><td>Green-Lagrange strain tensor</td></tr><tr><td>σ</td><td>Cauchy stress tensor</td></tr><tr><td>P</td><td>First Piola-Kirchhoff stress tensor</td></tr><tr><td>S</td><td>Second Piola-Kirchhoff stress tensor</td></tr><tr><td>W</td><td>Strain energy function</td></tr><tr><td>C</td><td>Fourth order elasticity tensor</td></tr><tr><td>N</td><td>Unit outward normal in the reference configuration</td></tr><tr><td>n</td><td>Unit outward normal in the current configuration</td></tr><tr><td>Π<sub>n</sub></td><td>Normal contact pressure</td></tr><tr><td>Π<sub>t</sub></td><td>Tangential stress vector</td></tr><tr><td>u<sub>n</sub></td><td>Normal displacement</td></tr><tr><td>u<sub>t</sub></td><td>Tangential displacement</td></tr><tr><td>g<sub>n</sub></td><td>Normal gap function</td></tr><tr><td>s</td><td>Sliding direction</td></tr><tr><td>μ</td><td>Coulomb friction coefficient</td></tr><tr><td>u̇<sub>t</sub></td><td>Sliding velocity</td></tr><tr><td>V</td><td>Translation velocity</td></tr><tr><td>ω</td><td>Rotation speed</td></tr><tr><td>r<sub>ω</sub></td><td>Rolling radius</td></tr><tr><td>d</td><td>Deflection</td></tr><tr><td>ε<sub>n</sub></td><td>Normal penalization</td></tr><tr><td>ε<sub>r</sub></td><td>Regularization coefficient</td></tr><tr><td>Δt</td><td>Time step</td></tr><tr><td>R<sub>e</sub></td><td>External radius of the prototype</td></tr><tr><td>R<sub>i</sub></td><td>Internal radius of the prototype</td></tr><tr><td>w</td><td>Width of the prototype</td></tr><tr><td>(λ, μ)</td><td>Lamé constants</td></tr><tr><td>E</td><td>Young Modulus</td></tr><tr><td>ν</td><td>Poisson coefficient</td></tr><tr><td>ρ</td><td>Density</td></tr><tr><td>H<sup>1</sup>(Ω)</td><td>First order Sobolev space</td></tr></table>

**References**

1. Sandberg, U.; Ejsmont, J.A. *Tyre/Road Noise Reference Book*; Informex: Kisa, Sweden, 2002.

2. Li, T. Influencing parameters on tire-pavement interaction noise: Review, experiments and design considerations. *Designs* **2018**, *2*, 38. [CrossRef]
3. Li, T.; Burdisso, R.; Sandu, C. Literature review of models on tire-pavement interaction noise. *J. Sound Vib.* **2018**, *2*, 357–445.
[CrossRef]

4. Ling, S.; Yu, F.; Sun, D.; Sun, G.; Xu, L. A comprehensive review of tire-pavement noise: Generation mechanism, measurement methods, and quiet asphalt pavement. *J. Clean. Prod.* **2021**, *287*, 125065. [CrossRef]

5. Nackenhorst, U. The ALE-formulation of bodies in rolling contact: Theoretical foundations and finite element approach. *Comput. Methods Appl. Mech. Eng.* **2004**, *193*, 4299–4322. [CrossRef]

6. Brinkmeier, M.; Nackenhorst, U.; Volk, H. A Finite Element Approach to the Transient Dynamics of Rolling Tires with Emphasis on Rolling Noise Simulation. *Tire Sci. Technol.* **2007**, *35*, 165–182. [CrossRef]

7. De Gregoriis, D.; Naets, F.; Kindt, P.; Desmet, W. Development and validation of a fully predictive high-fidelity simulation approach for predicting coarse road dynamic tire/road rolling contact forces. *J. Sound Vib.* **2019**, *452*, 147–168. [CrossRef]

8. Gonzalez Diaz, C.; Kindt, P.; Middelberg, J.; Vercammen, S.; Thiry, C.; Close, R.; Leyssens, J. Dynamic behaviour of a rolling tyre: Experimental and numerical analyses. *J. Sound Vib.* **2016**, *364*, 147–164. [CrossRef]

9. Pinay, J.; Saito, Y.; Mignot, C.; Gauterin, F. Understanding the contribution of groove resonance to tire-road noise on different surfaces under various operating conditions. *Acta Acust. 2020*, *4*, 1–12. [CrossRef]

10. Ejsmont, J.A.; Sandberg, U.; Taryma, S. Influence of tread pattern on tire/road noise. *SAE Transactions* **1984**, *93*, 632–640.
[CrossRef]

11. Li, T. Influence of tread pattern on tire/road noise. *Automot. Tire Noise Vib.* **2020**, 27–41. [CrossRef]

12. Haaret, P. Méthodes Numériques Pour la Dynamique des Structures non Linéaires Incompressibles à Deux Échelles. Ph.D Thesis,
Ecole Polytechnique, Palaiseau, France, 2004.

13. Gonzalez, O. Exact energy and momentum conserving algorithms for general models in nonlinear elasticity. *Comput. Methods Appl. Mech. Eng.* **2000**, *190*, 1763–1783. [CrossRef]

14. Biot, M.A. *Mechanics of Incremental Deformations*; John Wiley & Sons, Inc.: New York, NY, USA; London, UK; Sydney, Australia, 1965. [CrossRef]

15. Valyaev, V. On approaches to simulation of rolling tire vibrations. *Michelin Intern. Rep.* **2016**.

16. Ogden, R.W. *Non-Linear Elastic Deformations*; Dover Publications Inc.: Mineola, NY, USA, 1997.

17. Signorini, A. Questioni di elasticita non linearizzata e semilinearizzata. *Rend. Mat. Delle Sue Appl.* **1959**, *18*, 95–139.

18. Yastrebov, V. Computational Contact Mechanics: Geometry, Detection and Numerical Techniques. Ph.D Thesis, École Nationale Supérieure des Mines de Paris, Paris, France, 2011.

19. Wriggers, P. Constitutive Equations for Contact Interfaces. In *Computational Contact Mechanics*; Springer: Berlin/Heidelberg, Germany, 2006; pp. 69–108. [CrossRef]

20. Simo, J.C.; Laursen, T.A. An augmented lagrangian treatment of contact problems involving friction. *Comput. Struct.* **1992**, *42*, 97–116. [CrossRef]

21. Fortin, A.; Garon, A. Les Éléments Finis: De la Théorie à la Pratique. GIREF 2020. Available online: https://giref.ulaval.ca/ afortin/elements_finis.pdf (accessed on 14 December 2021).

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

23. Jazar, R.N. Tire Dynamics. *Veh. Dyn.* **2014**, 95–163. [CrossRef]

24. Chipato, E.T.; Shaw, A.D.; Friswell, M.I. Nonlinear rotordynamics of a MDOF rotor-stator contact system subjected to frictional and gravitational effects. *Mech. Syst. Signal Process.* **2021**, *159*, 107776. [CrossRef]

25. Zhou, P.; Du, M.; Chen, S.; He, Q.; Peng, Z.; Zhang, W. Study on intra-wave frequency modulation phenomenon in detection of rub-impact fault. *Mech. Syst. Signal Process.* **2019**, *122*, 342–363. [CrossRef]

26. Mlika, R.; Renard, Y.; Chouly, F. An unbiased Nitsche's formulation of large deformation frictional contact and self-contact. *Comput. Methods Appl. Mech. Eng.* **2017**, *325*, 265–288. [CrossRef]

27. Sundaram, N.; Farris, T. Mechanics of advancing pin-loaded contacts with friction. *J. Mech. Phys. Solids* **2010**, *58*, 1819–1833.
[CrossRef]

28. Scaraggi, M.; Persson, B. Friction and universal contact area law for randomly rough viscoelastic contacts. *J. Phys. Condens. Matter* **2015**, *27*, 105102. [CrossRef] [PubMed]

29. Alaci, S.; Muscà, I.; Pentiuc, S.G. Study of the Rolling Friction Coefficient between Dissimilar Materials through the Motion of a Conical Pendulum. *Materials* **2020**, *13*, 5032. [CrossRef] [PubMed]