# Numerical modelling of thermomechanic couplings in frictional contact: application to tyre rolling

Ange Kongo-Kondé, Iulian Rosu, Frédéric Lebon, Natacha Cocheteau, Lucie Seguin

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

## To cite this version:

Ange Kongo-Kondé, Iulian Rosu, Frédéric Lebon, Natacha Cocheteau, Lucie Seguin. Numerical modelling of thermomechanic couplings in frictional contact: application to tyre rolling. The Tenth International Conference on Computational Structures Technology, Sep 2010, Valence, Spain. pp.1-18. hal-00483321

HAL Id: hal-00483321

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

Submitted on 24 Mar 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical Modelling of
Thermomechanic Couplings in Frictional Contact:
Application to Tire Rolling

A. Kongo Kondé, I. Rosu, F. Lebon, N. Cocheteau and L. Seguin
LMA-CNRS, Marseille, France

## Abstract

This paper presents the modelling of a body in rolling contact taking into account 3D finite deformation (Mooney-Rivlin or Yeoh), dynamics, heterogeneities of the material, friction, thermomechanical coupling in the contact zone: the material properties do not depend on temperature, but friction law takes thermal effect into account. We consider two numerical approaches: a classical one, using Lagrangian formulation for the mechanical and thermal problems and the steady-state rolling approach (Eulerian/Lagrangian approach) for the mechanical problem associated to a Lagrangian approach for the thermal problem

**Keywords:** finite element, steady state rolling, temperature, thermal-mechanical analysis.

## 1 Introduction

In any vehicle simulation studies, tires are as among the most whimsical components. In vehicle dynamics context, tire models allow to predict tire rolling behavior such as lateral force $\mu_Y$ and self-aligning torque $M_z$ generated as result of the interaction between tires and the road surface. Such models require a combination of contradictory attributes such as computational efficiency, robustness of algorithms, accuracy, and ease of parameter identification. Over the past 50 years, parallel to tire models have been developed for vehicle ride and handling analysis. These models range from clever curve-fitting approaches, such as the “Magic Formula” [10], to relatively detailed models of physics of tire-road interaction. The latter benefitting from the ever-increasing computational power, have tried to integrate a large number of tire behavior attributes including structural deformation, friction, and thermal effects.
Current state-of-art tire models have reached a good level of maturity, but they often struggle to deliver accurate results, especially when operating under conditions substantilally different to those present when the model parameters were identified. One of the possible causes for such discrepancies is the delicate relantionship between tire temperature and road-rubber friction. So far, in the name of computational efficiency, tire models for vehicle dynamic analysis with thermal effects have been based on rather crude representations of heat transfer, involving mostly bulk or one-dimensional heating models. The relationship between the predicted temperatures and the generated friction is usually based on experimental measurements within a limited envelope of operating conditions.

Despite the substantial progress achieved in the area of computational mechanics, finite element analysis of pneumatic tires is still a challenging research topic in the fields of tire engineering. This method, as well as mathematic analytical models like "Magic Formula" [10], mechanical models [2, 1] which define pneumatic as circular ring in which the rigid rim is center at O and connects to a fixed framework by springs and dampers, SEA models based on vibration energies flow transfert, can lead to better understanding of tire behavior under different loading conditions such as inflation [9], static vertical load [3], rolling, cornering [4], thermal simulation [7], without any requirements for producing and testing tire prototypes.

On the one hand the determination of thermomechanical vehicle behavior is mostly high complex and requires hard work due to tire structural complexity (a deep knowledge of the geometry, materials properties, frictional coefficient, dissipative mechanisms, thermal coefficient and other aspects of tire modelling) and on the other hand the requirements of solving a non-linear dynamic thermo-viscoelastic coupling problem with heat sources resulting from internal sources and frictional contact.

Indeed, the ratio $\mu_Y$ (Lateral force over vertical one) debrief the equilibrium which exists between steering lock action and the trajectory of vehicle. During $\mu_Y$ upward phase, there is an equilibrium, the more the driver turn the wheel the more $\beta$ increases. It corresponds to pseudo-sliding in the contact up to $\beta_{max}$ (around $\beta=4$ to $7^\circ$ for automobile pneumatics, around $\beta=6$ to $10^\circ$ for trucks and around $18^\circ$ for aircrafts, see Fig. 1).

Beyond these values, a loss of equilibrium set up: the more the driver try to increase leeway, the more the engin skides. This phase corresponds to pure sliding in the contact. the rubber is warming up, the temperature increases and implies decreasing of the friction coefficient $\mu$. Thus, we observe the decreasing of $\mu_Y$ as shown in Fig.1. The specific values are normalized for analysing the observed phenomena

In skid, the result of friction forces is in the rear of contact area center. This intrinsic property of rubber creates resistant couple which faces drift: this is the
<figure><img src="image_3.png" /><figcaption>Figure 1: $M_Z$ and $\mu_Y$ for aircraft tire.</figcaption></figure>

*self-aligning Torque*, $M_z$. It represents the pneumatic ability to return to his steering direction when the driver releases vehicle controls. It is maximal for slip angle $\beta=6^\circ$ and rapidly vanishes towards zero (around $\beta=18^\circ$) for aircraft tire, see figure 1.

## 2 Theoretical foundations

Finite element analysis of rolling tires is generally carried out using three different numerical algorithms. In the first method, the implicit time integration technique is used and thus the time effect is directly incorporated into the finite element model equations. This technique requires a large amount of computer resources and thus it is generally used in conjunction with special kinematic relations to compute the steady state response of rotating tire. This procedure is known as steady state rolling method. It has been proven that this method gives stable and reliable results and also requires much less CPU time during computation stage. The third approach is the explicit time integration technique for simulation of dynamic behavior of a rolling tire under transient conditions. This technique can model the real time behavior of a running tire under operating conditions. However, in order to obtain non-oscillatory and consistent results, this method requires a very fine mesh as well as small time increments. The computational cost and effort in this method is also much higher than the steady state rolling technique. Consequently, the steady state rolling formulation was employed in the present study to overcome the disadvantages of both transient implicit and explicit methods.

The first and the third methods use Lagrangian formulation while steady state rolling method use a mixed Lagrangian/Eulerian formulation based on the ALE-for-
mulation.

## 2.1 Equilibrium problem including contact with friction in large displacements

Let us consider $\Omega$ a field of $\mathbb{R}^3$ representing an elastic solid of which the boundary $\partial\Omega$ (or $\Gamma$) is shared in $\partial\Omega=\Gamma_F \cup \Gamma_D$ where we impose, respectively, forces and displacements. Forces acting on solid consist of volume forces written down $g$ defined in $\Omega$ and surpeficial forces written down $f$ defined on $\Gamma_F$. We are looking for stable equilibrium positions of this solid.

<figure><img src="image_4.png" /><figcaption>Figure 2: Elasticity problem in large displacement with contact</figcaption></figure>

### 2.1.1 Dynamic problem taking into account inertial effects

Taking inertial effects into account, the dynamic problem formulation is written, at each time, as:

$$ \left\{ \begin{array}{ll} \rho \frac{\partial^2 u}{\partial t^2} = \text{div} T + g & \text{in } \Omega \\ Tn = f & \text{on } \Gamma \end{array} \right. \qquad (1) $$

where $T$ is the stress tensor and $u$ the displacement field. $\rho$ in reference configuration depend on spatial coordinates. This system of evolution is completed by suitable initial conditions.

### 2.1.2 Including contact with friction

For taking the fact that $\Omega$ can come in contact with a rigid base into account, we build a new partition of the boundary $\Gamma = \Gamma_D \cup \Gamma_F \cup \Gamma_C$ by introducing $\Gamma_C$ the potential contact area (see Fig. 2). We suppose that the rigid base is flat, the unit vector normal is $n^0$.

Unilateral contact conditions are:
$$ u_N \le 0 \quad \text{on } \Gamma_C \qquad (2) $$

$$ \lambda_N \le 0 \quad \text{on } \Gamma_C \qquad (3) $$

$$ u_N \lambda_N = 0 \quad \text{on } \Gamma_C \qquad (4) $$

We now consider the dynamic problem in the case where friction forces acting on $\Gamma_C$ following a Coulomb law with a constant friction coefficient c.

This law is applied in the case of dry friction which present a threshold phenomenon: in order that a contact point slips, the tangential stress must reach a fixed value.

Besides, the frictional stress $\lambda_T$ is proportional to normal stress and acts in opposite relative tangential speed direction between both surfaces in contact.

Using the unilateral contact conditions (2), (3) and (4), the conditions related to Coulomb friction are:

$$ \left\{ \begin{array}{ll} u_N < 0 & \text{and } \lambda_T = 0 \\ u_N = 0 & \text{and } \lambda_T > 0 \text{ with } \begin{cases} |\lambda_T| < c|\lambda_N| & \text{and } \dot{u}_T = 0 \\ |\lambda_T| = c|\lambda_N| & \text{and } \exists \theta \ge 0 \text{ such as } \dot{u}_T = -\theta\lambda_T \end{cases} \end{array} \right. \qquad (5) $$

### 2.1.3 Complete Continuum Model

**Problem 1 (Dynamic continuum Problem)** Find displacement u, hydrostatic pressure defined in $\Omega$ and contact forces with friction $\lambda=(\lambda_N, \lambda_T)$ defined on $\Gamma_C$ such as, at each time, we have dynamic equilibrium of structure.

The continuum model is written in its weak form as [8]:

$$ u - \bar{u} \in V, \quad p \in P, \quad \lambda_N \in \Lambda_N, \quad \lambda_T \in \Lambda_T(s) \qquad (6) $$

$$ \int_{\Omega} \rho \ddot{u} \cdot v d\Omega + \int_{\Omega} \frac{\partial W_e}{\partial F}(I + \nabla u) : \nabla v - p \frac{\partial \det}{\partial F}(I + \nabla u) : \nabla v d\Omega - \int_{\Gamma_C} \lambda_T v_T d\Gamma \\ - \int_{\Gamma_C} \lambda_N v_N d\Gamma = \int_{\Omega} g \cdot v d\Omega + \int_{\Gamma_F} f \cdot v d\Omega + \int_{\Gamma_P} \nabla V_P(v) \cdot v d\Gamma \quad \forall v \in V \qquad (7) $$

$$ \int_{\Omega} q(\det(I + \nabla u) - 1) d\Omega = 0 \quad \forall q \in P \qquad (8) $$

$$ \int_{\Gamma_C} (\mu_N - \lambda_N) u_N d\Gamma \ge 0 \quad \forall \mu_N \in \Lambda_N \qquad (9) $$

$$ \int_{\Gamma_C} (\mu_T - \lambda_T) \dot{u}_T d\Gamma \ge 0 \quad \forall \mu_T \in \Lambda_T(s) \qquad (10) $$

$$ s = c|\lambda_N| \quad \text{on } \Gamma_C \qquad (11) $$
Where $\bar{u}$ is a known field of displacement which include imposed displacement conditions $\Gamma_D, \Gamma_P$ is the part of boundary (defining the enclosed volume) where in-flating pressure forces are exerted.

$\bar{u}, f$ and $g$ are time dependent.

This is a continuum formulation of a mechanical problem that we are going to resolve numerically. It is a non-linear dynamical problem which challenges incom-pressible materials in large transformation, unilateral contact with Coulomb’s friction as well non-linear pressure forces. The mathematical analysis of such problem is very difficult and still implies many open questions on the existence and uniqueness of the solution.

## 2.2 ALE-Formulation or Steady State Transport Method

For the computational treatment of rolling contact of deformable bodies a relative kinematic representation is needed, for which an ALE-description has been found to be adequate. As mentioned already in the introduction for different applications alternative ALE-approaches or mixte Lagrangian/Eulerian approaches have been de-veloped [6], all in common with the following characteristics.

* an arbitrary reference configuration is introduced which is neither fixed in space nor connected with the material particles,

* the boundaries of the domain are moving like the material boundaries.

The ALE-description is interpreted as the decomposition of motion into pure rigid body motion and the deformation, see Fig.3. The initial configuration is described by a configuration $\mathfrak{B}(X)$ of the material particles in Euclidian space $\mathfrak{B}^3$. By $\phi(X, t)$ this material configuration id mapped into the current configuration by following the motion of the material particles. Now, this motion is decomposed into a mapping $\chi(X, t)$, which describes the pure rigid body motion of rolling bodies in a Eulerian manner. With respect to this reference configuration the deformation is described by the mapping $\hat{\phi}(\chi, t)$. Thus, the material particles at time $t$ take place at the location $x = \phi(X, t)$. Confining oneself on the observation of the deformation only, $\varphi = \hat{\phi}(\chi, t)$ describes the location where the coincident points of the reference configuration take place at time $t$. Illustratively, the ALE-observer introduced by this mathematically abstract formulation can be interpreted as riding on the axis of the rolling wheel (without rotating).

$X$ are variables related to the material configuration$\mathfrak{B}(x)$, $x$ variables related to the current configuration $\phi(\mathfrak{B})$ and $\hat{\phi}$ which related to the reference configuration.
<figure><img src="image_5.png" /><figcaption>Figure 3: ALE-decomposition of motion</figcaption></figure>

### 2.2.1 Weak Formulation

So far, the problem is stated by the balance of momentum

$$ \dot{\rho}v = \rho b + \mathrm{Div} \hat{\mathbf{P}}, \qquad (12) $$

where $\hat{\mathbf{P}}$ is the Piola-Kirschhoff stress tensor in the reference configuration.

and the boundary conditions : Dirichlet conditions

$$ \psi = \bar{\psi} \quad \text{on} \quad \partial_{\psi} \chi(\mathfrak{B}), \qquad (13) $$

and Neumann conditions

$$ \hat{\mathbf{P}} \cdot \hat{n} = \bar{t} \quad \text{on} \quad \partial_t \chi(\mathfrak{B}), \qquad (14) $$

at disjunct parts of the boundary $\partial\chi(\mathfrak{B}) = \partial_\varphi\chi(\mathfrak{B}) \cup \partial_t\chi(\mathfrak{B})$, $\partial_\varphi\chi(\mathfrak{B}) \cap \partial_t\chi(\mathfrak{B}) = 0$. The contact conditions on $\partial_c\chi(\mathfrak{B}) \in \partial_t\chi(\mathfrak{B})$ which is not known a priori will be introduced later.

Note, that so far the description is independent from the particular material behavior. For solution of real life problems a phenomenological abstract constitutive stress-strain relation has to be introduced in addition, which for a generalized standard material is written formally.

$$ \sigma = \rho \frac{\partial \Psi(F, \alpha, \theta)}{\partial F} \cdot F, \qquad (15) $$

A weak formulation of the problem is obtained by weighting the balance of momentum in addition with the Neumann conditions with a test function $\eta$ and integrating
over the volume of the body.

$$
\begin{equation}
\begin{split}
\int_{\chi(\mathfrak{B})} \hat{\rho} \frac{dv}{dt} \cdot \eta d\hat{v} = & \int_{\chi(\mathfrak{B})} \hat{\rho} \left( \left. \frac{\partial^2 \psi}{\partial t^2} \right|_{\chi} \cdot \eta \right) d\hat{v} + \int_{\chi(\mathfrak{B})} \hat{\rho} \left( \eta \cdot \mathrm{Grad} \left. \frac{\partial \phi}{\partial t} \right|_{\chi} \cdot w - \left. \frac{\partial \phi}{\partial t} \right|_{\chi} \cdot \mathrm{Grad} \eta \cdot w \right) d\hat{v} \\
& - \int_{\chi(\mathfrak{B})} \hat{\rho} \left( (\mathrm{Grad} \psi \cdot w) \cdot (\mathrm{Grad} \psi \cdot w) \right) d\hat{v} \\
& + \int_{\chi(\mathfrak{B})} \hat{\rho} \eta \left( \left. \frac{\partial \phi}{\partial t} \right|_{\chi} + \mathrm{Grad} \psi \cdot w \right) w \cdot \hat{n} d\hat{a}
\end{split}
\tag{16}
\end{equation}
$$

Because of limited computational capabilities even today most of the applications
are focused onto the stationary motion (all relative time derivatives vanish) of rolling
bodies. In this case,

$$
\int_{\chi(\mathfrak{B})} \hat{\rho} \frac{dv}{dt} \cdot \eta d\hat{v} = - \int_{\chi(\mathfrak{B})} \hat{\rho}(\mathrm{Grad}\psi \cdot w) \cdot (\mathrm{Grad}\eta \cdot w) d\hat{v} + \int_{\chi(\mathfrak{B})} \hat{\rho}\eta (\mathrm{Grad}\psi \cdot w) w \cdot \hat{n} d\hat{a} \quad (17)
$$

The contributions from normal and tangential contact are decoupled and the virtual
work of the contact forces is written as

$$
W_c = - \int_{\partial_c \hat{\phi}(\mathfrak{B})} p \delta d_n da - \int_{\partial_c \hat{\phi}(\mathfrak{B})} \tau_\alpha \delta \zeta^\alpha da, \quad (18)
$$

The basic difficulty in solving rolling contact problem is the treatment of the tangential sub-problem. This is because the relative displacement of material particles in contact needed for the evaluation of the constitutive laws is not observed directly within the ALE-description of rolling. The slip distance $\zeta^\alpha$ of two contacting material particles, which is the tangential relative movement that has been appeared since they got into contact at the leading edge, is not computed directly. A direct method to overcome this problem is the computation of the path of motion of the particles and the integration of the history along these lines.

For stationary motions, the tangential contact is written as:

$$
\int_{\partial_c \hat{\phi}(\mathfrak{B})} \tau_\alpha \delta \zeta^\alpha da = \int_{\Gamma_c \phi(\mathfrak{B})} \lambda_\alpha \delta \zeta^\alpha \tilde{c} \cdot n_\Gamma d\Gamma - \int_{\partial_c \hat{\phi}(\mathfrak{B})} \lambda_\alpha \delta \dot{\zeta}_c^\alpha da, \quad (19)
$$

2.3 Thermomechanical coupling

2.3.1 Heat equation

During a contact with friction between two bodies (solids), there is generation of heat
in the contact due to energy dissipation. In a general frame, solve a thermal problem
leads to find solution of the heat equation (20) with boundaries conditions: Dirichlet and Neumann (convectvity, adiabacity and given temperature).

$$ \rho c_p \frac{dT}{dt} = \operatorname{div}(k \cdot \nabla T) + Q \qquad (20) $$

where $\rho$ is the mass density (kg/m³), $c_p$ is the specific heat (Jkg⁻¹K⁻¹), $k$ is the thermal conductivity (Wm⁻¹K⁻¹), $T$ is the temperature (K), $Q$ is the internal heat sources (Wm⁻³) and $t$ the time(s).

### 2.3.2 Internal heat sources

The internal heat sources correspond to internal energy source arise from mechanical contribution. It exists two main sources: plastic strain and the contact with friction. The source heat flux $Q$ can be written as function of a part $\eta$ of plastic strain power and power dissipated in contact.

$$ Q = \eta\sigma : \varepsilon^P - \sigma_c^t \cdot \dot{u}_T \qquad (21) $$

with $\sigma$ the stress tensor, $\varepsilon^P$ the strain rate plastic tensor, $\sigma_c^t$ tangential stress in the contact and $\dot{u}_T$ tangential speed in the contact.

### 2.3.3 Weak formulation

We can now write the variational weak form related to heat equation (20) which is rewritten as:

$$ \int_{\Omega_1 \cup \Omega_2} \rho c_p \dot{T} v d\Omega + \int_{\Omega_1 \cup \Omega_2} k \nabla T \nabla v d\Omega = - \int_{\Gamma_q} k \nabla T \vec{n} v d\Gamma + \int_{\Omega_1 \cup \Omega_2} Q v d\Omega \quad (22) $$

with $v$ a test function and $\vec{n}$ the external normal to the surface $\Gamma_q$.

With Fourier's law (23), we introduce the flux applied on surface by convection.

$$ k \nabla T \vec{n} = q_c = \alpha_c (T - T_{\infty}) \qquad (23) $$

By introducing the mechanical coupling terms, the weak formulation of heat equation becomes

$$ \int_{\Omega_1 \cup \Omega_2} \rho c_p \frac{dT}{dt} v d\Omega + \int_{\Omega_1 \cup \Omega_2} k \nabla T \nabla v d\Omega = - \int_{\Gamma_q} \alpha_c (T - T_{\infty}) v d\Gamma + \int_{\Omega_1 \cup \Omega_2} Q v d\Omega \quad (24) $$

with
$$ \int_{\Omega_1 \cup \Omega_2} Q v d\Omega = \int_{\Omega_1 \cup \Omega_2} \eta(\sigma : \varepsilon^P) v d\Omega - \int_{\Gamma_c} \sigma_c^t \cdot \dot{u}_T v d\Gamma \quad (25) $$

# 3 Modelling and simulation of pneumatic Tires

## 3.1 Constraints to be considered

Modern pneumatic tires consist of a specific combination of rubber compounds, cord and steel belts (see Fig. 4). The main parts of a modern pneumatic tire are its body, sidewalls, beads, and tread. The body is made of rubberized fabric layers, called plies, that give the tire strength and flexibility. The fabric is made of rayon, nylon, or polyester cord. Covering the plies are sidewalls and tread of chemically treated rubber. The sidewalls form the outer walls of the tire. Embedded in the two inner edges of the tire are steel loops, called beads, that hold the tire to the wheel. The rubber components have different characteristics in dependence on their functionality. The tread for example gets in direct contact with the ground and has to be much harder than the sidewalls.

<figure><img src="image_6.png" /><figcaption>Figure 4: Automobile tire components</figcaption></figure>

## 3.2 Finite Element Model

The purpose of the present research is to predict tire rolling behavior of a 175 SR14 tire in contact with a flat rigid surface, subjected to an inflation pressure and a concentrated load on the axle.

In the present study, finite element analyses were performed in four consecutive stages. In the first stage, the tire was modeled under inflation pressure using 2D ax-
isymmetric model. Then a 3D model was developed based on the 2D model for the simulation of the tire under static vertical load (footprint stage). Either the dynamic analysis with acceleration as third stage and cornering rolling as the fourth stage, or the steady state analysis was then carried out in the third stage using the restart option to obtain tire free rolling conditions. In the last stage, slip angle were chosen to simulate the tire under cornering conditions. For taking account of the thermal effects due to contact with friction, a thermomechanical analysis is performed. This stage follows the 3D static analysis.

### 3.2.1 Inflation analysis

Inflation analysis on axisymmetric model (see Fig. 5) is primarily performed because of model integrity checking, finding the deformed shape of the model, finding the stresses in carcass, belts and rubber components and setting up the basis for creation of 3D FE model. The element group CGAX4(H) and CGAX3(H) from the ABAQUS elementary library was selected to represent all rubber parts. These are four and three-nodded bilinear elements with twist and constant pressure to take the incompressibility of rubbery materials and nonzero circumferential component of the displacement tensor into consideration. Reinforcement materials were modelled with rebar in the surface element, SFMGAX1 which is two-nodded linear axisymmetric surface element with twist.

<figure><img src="image_7.png" /><figcaption>Figure 5: 2D axisymmetric model without grooves</figcaption></figure>

### 3.2.2 Analysis of vertically loaded tire

Inflation analysis of axisymmetric model is followed by 3D static analysis of vertically loaded tire. The analysis consists of several steps. In step zero, the 3D model is created, by rotation of axisymmetric model around tire axis, C3D8H and SFM3D4R elements are generated from 2D elements (CGAX4H and SFMGAX1). In step one, tire is inflated, while in step two ground surface gets gradually moved towards tire
axis, until it reaches the goal value. In step three equilibrium is reached, with constant force acting on surface. This step is needed as the basis for subsequent steady state rolling analysis or dynamic rolling analysis. Basic results of vertically loaded tire analysis are shown in Fig.6.

<figure><img src="image_8.png" /><figcaption>Figure 6: 3D axisymmetric model for SS and for RC/TM</figcaption></figure>

**3.3 Methods and Analysis approaches**

Two approaches [4] can be used to solve this kind of problem. We discuss about it above. The first one is called "steady state analysis". This is a mixed Lagrangian/Eulerian approach in which, for an observer in the reference frame, the tire is a fixed set while materials flows through the refined stationary mesh. The second one is called "Classical Dynamic analysis". This Lagrangian approach allows to simulate real rolling. In this approach, for an observer in the reference frame, the tire rolls and, turn after turn, the elements touch and leave the contact area.

**3.3.1 Steady State Analysis "SS"**

* **Steady State rolling analysis** : In the third analysis, the tire was modeled under steady state rolling conditions. This stage is restarted from the results of the second step (footprint analysis). A constant given ground velocity was assumed for the tire. The tire spinning velocity was changed to determine the tire free rolling condition and its dynamic loaded radius. In addition, the friction between tread surface and ground was taken into account using the simple linear Coulomb's law.
* **Rolling analysis with slip angle**: The last analysis step is devoted to the application of the slip angle to the tire for a cornering situation. Slip angle is defined as the angle between the direction of travel and a normal plane to the axis of rotation. Analyses were performed using these values in conjunction with tire free rotational speed found in the previous step.

### 3.3.2 Dynamic analysis or Classical Rolling "RC"

* **Acceleration analysis**: In the third analysis, the tire was modeled under dynamic conditions. This stage is restarted from the results of the second step (footprint analysis). An acceleration is obtained by using a ramp which increase the velocity from zero to the given velocity. This is important to stabilize the model and the velocity before cornering. In addition, the friction between tread surface and ground was taken into account using Coulomb's law.

* **Constant velocity in cornering**: The last analysis was devoted to the application of the constant given velocity at fixed slip angle $\beta$. Slip angle is defined by defining longitudinal and lateral velocities such as the given ground velocity $v_0$ is decomposed in $v_x=v_0cos\beta$ and $v_y = v_0sin\beta$.

### 3.3.3 Thermal-Mechanical analysis

The thermal-Mechanical analysis can be performed in different methods.

* The first method consist in doing a full coupling [12]. In this method, temperature variation during a time increment $\Delta t$ has been taking into account and steps in as additional nodal unknown of the system. Such method is often related to implicit algorithm.

* The second method consist in doing a sequential coupling: at each time increment $\Delta t$, we carry out the resolution of isotherm mechannical system with given temperature at $t - \Delta t$ following by the resolution of thermal system with fixed geometry being deteminated. In this case, we consider that there is no temperature change during the resolution of mechanical system and we suppose that the time increment is small. Such method is used with explicit integration scheme.

# 4 Results and Discussion

## 4.1 Steady State vs Classical Rolling

We have discussed about the advantages and disadvantages of the two methods. Here is a comparison of the two methods. Numerically, both methods are equivalents as
shown in Fig. 8. In term of CPU time, 10 hours are necessary to compute one point in dynamic response while 5 hours are needed to compute total response in steady state. All calculations are done on a HP sever with two Xeon quad-core processors cadenced at 3GHz and using the parallel code version.

<figure><img src="image_9.png" /><figcaption>Figure 7: Lateral Force and Self Aligning Torque in RC</figcaption></figure>

When we use classical rolling method, we just have one point for the chosen loading parameters. As Fig. 7 shows, the tire response in this method is time-dependent with a transient phase at beginning before reaching a stationary one. So, we do not need several tire rotation to reach the stationary phase. $M_Z$ and $\mu_Y$ are evaluated from this stationary phase. We carry out so much simulation as there are loading conditions (points to be drawn). To fulfil all slip angles a lot of simulation and time is necessary while only one simulation is needed in steady state method Fig. 8.

<figure><img src="image_10.png" /><figcaption>Figure 8: Comparison of Steady State and Classical Rolling</figcaption></figure>

## 4.2 Fully coupled thermo-mechanical analysis

### 4.2.1 Pure sliding

The tire is pull on ground along 4 meters. In this case, we observe that the temperature increases significantly when tire is lengthways slided down even in short distance. The
rubber reaches a superficial temperature of 190°C (Fig. 9) close to critical temperature for which the rubber settle on concrete road [11].

<figure><img src="image_11.png" /><figcaption>Figure 9: Result of thermal simulation in pure sliding</figcaption></figure>

### 4.2.2 Cornering rolling

The cornering is simulated by imposing a longitudinal and lateral velocity. The cornering simulation is performed on two loading cases. Following a node on tread, we observe moderate temperature (around $T = 80°C$). In every instance, the rise in temperature is located on tread surface with temperature gradient in the tire thickness (see Fig. 10).

<figure><img src="image_12.png" /><figcaption>Figure 10: Result of thermal simulation in cornering</figcaption></figure>

The reinforced zone seems to constitute a thermal barrier. From now on, ABAQUS allows to carry out a complete simulation including thermal effects in the rolling tire behavior. However, the full coupled analysis remains too expensive. For fixed loading conditions (vertical load, velocity, inflation pressure and slip angle), 3½ days are necessary to cover the given distance. So a real second corresponds to about 604000 s of simulation.

Fig. 11 shows the rise in temperature of tread node for two velocities. The levels represent the node temperature from the exit to the new entry in contact while the slope characterizes the temperature increasing in the contact.
<figure><img src="image_13.png" /><figcaption>Figure 11: Rise in temperature of a tread node at $v=50$ km/h and 100 km/h</figcaption></figure>

In this simulation, we do not consider the convective effects between tire and surrounding air. It should depend on velocity. Naturally, we should have $h = f(V)$. The horizontal levels remind it.

## 4.3 Sequential coupled thermo-mechanical analysis

It is also possible to perform uncoupled analysis of a tire. The heat fluxes can be calculated and applied to a thermal model; then the resulting temperatures can be applied to a stress analysis. However, since the thermal and stress analyses are uncoupled, this approach does not take account of the effect of the thermal deformation on the contact which, in turn, affects the heat generation.

The traditional way of analyzing this kind of problem is to use a Lagrangian approach in which the mesh used to discretize the tire rolls relative on the ground. Since many revolutions are typically required to reach the state of interest to the analyst, this approach is prohibitively expensive and cumbersome. The steady-state transport analysis capability in Abaqus/Standard (Steady-state transport analysis), which uses the Eulerian method in which the finite element mesh of the tire does not roll relative to the ground but the material flows through the mesh, provides a cost-effective alternative approach. The paths that the material points follow through the mesh are referred to as streamlines.

The algorithm for sequential analysis is presented in Fig. 12.

## 5 Conclusion

This present research is a first step of thermomechanical coupling model for vehicle tire. The aim is to develop a cost-effective alternative thermomechanical approach. To develop a frictional law which consider thermal effects, we have to define precisely the temperature in the contact. Indeed, with taking these effects into account in the friction coefficient, we can emphasize the decreasing of $\mu_Y$ beyond $\beta_{max}$ where the
<figure><img src="image_14.png" /><figcaption>Figure 12: Algorithm of sequential analysis</figcaption></figure>

grip potential is reached. In the future, we will work on the proposed approach and
will implement this algorithm in ABAQUS.

References

[1] Dohrmann C.R., "Dynamics of a tire-wheel-suspension assembly", Journal of Sound and Vibration, Vol.210, pp.627-642,1998.

[2] Ferrarese A., Padovese L.R., Costa A.L.A., "Tire Dynamical Models", Computational Methods in Engineering, 1999.

[3] Ghoreish Reza M.H., "Finite Element Analysis of Steel-Belted Radial Tyre with Tread Pattern under Contact Load", Iranian Polymer Journal, Vol.15(8), pp.667-674, 2006.

[4] Kongo Kondé A., Rosu I., Lebon F., Brardo O., Devésa B., "Etude du comportement en roulement d'un pneu d'avion", Colloque National de Calcul de Structures, Giens 2009, Vol.1, pp. 699-704, 2009.

[5] Lahellec N., Mazerolle F., Michel J.C., "Second-order estimate of the macroscopic behaviour of periodic hyperelastic composites: theory and experimental
validation", Journal of the Mechanics and Physics of Solids, Vol.52, pp.27-49, 2004.

[6] Nackenhorst U., "The ALE-formulation of bodies in rolling contact - Theoretical foundations and finite element approach", Comput. Methods Appl. Mech. Engrg, 193, 4299-4322, 2004.

[7] Narasimha Rao K.V., Krishna Kumar R., bohara P.C., Mukhopadhyay R., "A Finite element Algorithm for the Prediction of Steady-State temperatures of Rolling Tires", Tire Science and Technology, TSTCA, 34(3), 195-214, 2006.

[8] Navarro J.P., "Contribution à la Modélisation du Pneumatique de l'Avion", Thèse: Université Toulouse III - Paul Sabatier, 2003.

[9] Olatunbosum O. A., Bolarinwa E. O., "Finite element simulation of the tyre burst test", Proc.Instn Mech.Engrs, 218, 1251-1258, 2006.

[10] Pacejka H.B., Bakker E. and Nyborg L., "Tyre Modelling for Use in Vehicle Dynamics Studies", SAE 870421, Int. Congress and Expo., Detroit, 1987.

[11] Sakai H., Araki K., "Thermal Engineering Analysis of Rubber Vulcanization and Tread Temperatures during Severe Sliding of a Tire", TSTCA, Vol.27, No 1, pp.22-47, January-March 1999.

[12] Van Der Lugt J., Huétink J., "Thermal mechanically coupled finite element analysis in metal-forming processes", Comp. Methods Appl. Mech. Eng., TSTCA, 29 (3), 145-160, 1986.