## Abstract

Over the past decades, simulation has become vital in vehicle development and virtual testing, especially for autonomous vehicles. High-performance hard real-time simulators, which are crucial for those undertakings, require efficient algorithms to accurately model vehicle behavior within virtual environments. A prime example is tire-ground contact modeling, which is pivotal if we aim to achieve a high level of realism when simulating wheeled vehicles. Contact modeling focuses on an accurate estimation of the parameters needed to compute the forces and torques generated by vehicle-ground interaction. However, the complexity of this task is compounded by the fact that tire-ground contact is a highly nonlinear phenomenon, which is further exacerbated by the need to perform tests to fine-tune state-of-the-art tire-ground contact models. To tackle those challenges, we have developed a novel enveloping model that does not require any fitting of experimental data and is based on the 3D geometry of the intersection between undeformed volumes. In this paper, we provide a detailed description of the algorithm's formulation, the current software implementation (available as an open-source library), as well as the achieved scalability and real-time performance.

## 1 Introduction

From the 1950s onward, vehicle simulation has been crucial for evaluating ride and handling. More recently, as we shift toward autonomous vehicles, there's a growing need for high-quality simulations to develop and test them in various scenarios. Real-world testing is time-consuming and costly, making driving simulators essential for accelerating the transition to autonomous driving, enhancing safety, and reducing expenses. Beyond the hardware, a simulator's essence lies in its software, responsible for shaping the virtual environment based on physics. Aside from computer graphics, this software comprises multiple subprograms, including a Multi-Body System (MBS)—used to integrate the Multi-Body (MB) car system—and several other subprograms that allow the MB system to interact with the environment. In the case of a road vehicle, interactions with the environment primarily occur through tire-ground contact, aerodynamic forces, and, less frequently, through collisions with external objects. To obtain a sufficiently realistic simulation of the vehicle it is therefore of crucial importance to correctly estimate the forces and torques generated by the tire-ground contact [1]. In order to be able to calculate the stress generated by the contact, the first step is to estimate the ground rut pose with respect to the tire reference frame [2]. An incorrect estimate of the local contact plane or surface would inevitably lead to unrealistic or unreliable simulations. Since tire-ground contact evaluation is one of the most time-consuming processes in MB simulations, algorithms designed specifically for this purpose are needed to ensure high efficiency and scalability. Specifically, in the case of Driver/Hardware/Software-In-The-Loop (DIL/HIL/SIL) simulations, this software must guarantee a capability that is at least as fast as real-time. It is also desirable to obtain a faster than real-time capability, in order to speed up costly offline simulations and to allow for the use of less powerful machines. As a matter of fact, the development of faster and at the same time more accurate contact models still remains an open research topic [3,4].

From the scientific literature, it is known that appropriate contact methods must be applied when road unevenness is characterized by wavelengths 2–3 times smaller than the contact patch length. This occurs when riding is done either over cobblestone or Belgian block roads, or simply over road surfaces that present sharp obstacles, e.g., cleats, bumps or potholes [1]. The role of the tire-ground contact model is to represent the excitation of pneumatic tires caused by uneven terrain surfaces. Although the term “uneven ground surface” comprehends any kind of unevenness, attention is mainly paid to short irregularities, as tire deformation is particularly important in those instances. The ability of the tire to deform, i.e., to follow the ground surface shape, is defined as the enveloping property of the tire-ground contact. In specific, the enveloping is the quasi-static component of the tire-ground contact phenomenon, whereas the dynamic component comes from the carcass flexibility and is implemented in the tire force-calculation model [1,5]. The enveloping property depends on the tire geometry and structure and can be represented through physical models such as those presented in Refs. [6–10], or (semi-)empirical models such as Refs. [5] and [11]. In contrast to physical models, the empirical models try to describe the contact through much simpler models designed to mimic the behavior of the tire-ground interaction.

In the early days of vehicle simulation, studies on tire behavior have mostly been limited to the case of flat and asperity-free contact surfaces, focusing on the deformation and dynamics of the tire carcass induced by the contact forces. The importance of accurately describing the local contact area has been highlighted in Refs. [2] and [12], where the influence of obstacles on the overall tire behavior is studied. While the tire is riding over 3D uneven road surfaces the local road plane shows variations in height, and in forward and banking angles. The banking slope angle gives rise to the so-called tire camber thrust, which must be taken into account in order to properly simulate the force generated by the tire. Since then, many enveloping models for arbitrarily uneven surfaces have been developed. Most of the models developed so far are based on heuristic concepts. Without a doubt, the most famous enveloping model is certainly the Swift^{®} model, which allows the Magic Formula (MF) model to be used even on rough road surfaces [5]. Simple and parameterless contact models like that used in TMEasy [11,13] can be applied, but do not always guarantee adequate enveloping properties, which can cause sudden and unrealistic variations of the banking and slope angles in the proximity of sharp cleats or asperities. Another family of more sophisticated enveloping models is based on the physical concept of radial and radial-inter-radial spring tire [6–8]. Models based on the tire radial-spring behavior are a good compromise between computational effort, robustness and physical consistency. Unlike the Swift® model, the radial spring tire models do not show any working range limit in shape, magnitude and orientation of the incoming road obstacles [6]. For this reason, the results are physically consistent even in the case of high cleats. The state-of-the-art approach is referred to as “full-physical tire modeling.” In particular, the commercial software FTire® [14] and CDTire® [3,15] are the leading virtual tire simulation models in this category. They are proven to be multipurpose physics-based tire models that are able to simulate nearly all tire dynamics phenomena. They combine a discrete elements description of both belt and sidewalls elements with a brush type contact, and real-time capabilities [16]. In terms of enveloping behavior, they are also providing physically consistent and extremely realistic results [15,17]. Unfortunately, little information on the actual computational time performance is available for both FTire® and CDTire®.

A very simple physical contact modeling approach we have not yet mentioned before relies on the geometry of intersection between undeformed regions, also known as displaced area models [9]. This kind of model is derived from the physics of the linear radial-spring tire model. As stated also in Ref. [6], a common approach for modeling the 2D vehicle-road contact is to assume that the tire is described by a disk moving on a ground region. Indeed, if we assume that the tire is represented by an infinite number of independent linear radial springs, the force resulting from the contact is linearly proportional to the intersection area and acts along the line passing through the centroid of the intersection region and the wheel center. Despite the simplicity of this approach and the limits that arise from neglecting the tire carcass being a unique cohesive body capable of transmitting shear forces, it provides a good approximation of the contact information. Models based on the geometry of intersection between undeformed regions are already available, but they all lack a full 3D description of the contact. This limits the comparison between other enveloping models to 2D contact scenarios. Moreover, no robust software implementation is freely available and can be used to determine the tire contact point and normal or equivalently a contact plane.

This work aims to present a new software-based algorithm (hereafter called Enve) to model the 3D geometrical contact between tire and ground, represented, respectively, as a generic axial-symmetric surface and a triangular mesh. In particular, the adopted methodology is derived from the physics of the linear radial-spring tire model and is based on the geometry of the intersection between undeformed regions. Indeed, Enve is an enveloping model that is fully parameterizable by the tire's external shape alone and does not require any fitting on experimental data. As will be shown in the following, the proposed approach achieves an efficient and robust modeling of the tire-ground geometrical contact with scalable precision. This model is intended to be coupled with a tire model, i.e., empirical or physical tire model, with or without belt and sidewall dynamics. The tire model calculates tire-ground forces based on geometric contact parameters supplied by the Enve software. In this study, we demonstrate the efficiency and efficacy of the software by coupling Enve with the MF 6.2 tire model.

*Remark 1*. *In this work, the focus will only be on the enveloping properties of the tire-ground contact, i.e., the calculation of the so-called effective road surface. To evaluate dynamic response, one can employ a rigid ring model coupled with a tire model for force calculation, as exemplified by MF-Swift ^{®} [1,5]. In this work, our aim is to determine the application point and direction of contact forces, and as such, we do not delve into the modeling of tire-ground contact forces or the tire carcass dynamics.*

This paper has two main goals. The first is to provide a methodology for the analysis of the geometrical tire-ground contact, which takes into account the geometry of the tire cross section and of the road surface to compute information needed by the tire model to obtain forces. The proposed enveloping model offers advantages over Swift^{®} in terms of algorithm simplicity and scalability. It excels in computing contact on multiple sections, making it suitable for supporting physical-based tire models like the brush model. Additionally, it does not rely on independently identified parameters, further enhancing its effectiveness and immediacy of use. These characteristics make it suitable for hard real-time applications where computational efficiency is the main concern without diminishing the level of accuracy. The second is to evaluate the impact of this new approach coupled with a tire model in an advanced simulation environment with hard real-time scheduling. The software, which is also called Enve, is built on the Acme geometry library [18] and specifically designed to be integrated into vehicle dynamics MB kernel codes as well as into simpler simulation environments [19,20].

## 2 Tire-Ground Enveloping Model

**H**is located in the wheel hub (see Fig. 1). The axes are oriented according to the ISO 8855:2011 standard [21]. The

*x*-axis is directed toward the longitudinal direction of motion, the

*z*-axis points upwards and the

*y*-axis is oriented in such a way that the coordinate system is right-handed. The tire is defined geometrically as an axially symmetric (around the axis $e\u0302y$), convex and closed set $T$

*R*(

*y*) is chosen arbitrarily to properly reconstruct the tire's external morphology. On the other hand, the ground is modeled as a closed half-space $G$

whose interior and boundary are, respectively, denoted with $Go$ and $\u2202G$. The map *g*(*x*, *y*) is a mapping between $\mathbb{R}2$ to $\mathbb{R}$ that defines the height of the track surface. Similarly, the map *f*(*x*, *y*) defines the local friction coefficient scaling factor, which is used to scale the tire performance according to the local ground surface conditions. Notice that the concept behind the friction coefficient scaling factor can be extended to other relevant local properties of the ground surface, e.g., temperature, water film thickness, and asphalt wear.

*A*and the contact volume

*V*are defined, respectively, as

The tire is assumed to be modeled by an infinite number of independent linear radial-springs, whose stiffness *k* is uniform all over the tire section. The inaccuracy introduced by this assumption is recovered by the tire force-calculation model, which employs concentrated parameters to describe the nonuniformity of the carcass stiffness. Before calculating the direction $n\u0302$ and application point **P** of the resultant contact force, we have to make three important observations.

*The force generated by the contact is produced by the deformation of the radial springs. If the terrain is nondeformable the contact force magnitude*$||F||$

*is linearly proportional to the intersection volume V,*i.e.,

Notice that the proportionality constant is the stiffness coefficient *k*, which is a function of the tire's internal structure and pressure. However, these dependencies are considered by the tire force-calculation model through lumped parameters. It should be pointed out that the damping effect is not considered due to the quasi-static nature of the enveloping phenomenon.

Observation 2. *Since the tire is assumed to be modeled by an infinite number of independent radial springs, the contact force*$F$*cannot produce any rolling resistance torque around the wheel hub axis*$e\u0302y$*. In other words, the contact force acts along the direction determined by one of the families of lines passing through the contact point**P**and the hub axis*$e\u0302y$.

It is also important to note that the rolling resistance is not considered in the enveloping models. It is typically considered to be part of the tire contact force model, which must be fed with the information provided by the enveloping model. For instance, the Swift^{®} model is coupled with the MF model to compute the contact forces and torques generated by the tire-ground contact.

**P**is calculated as the weighted average of contact surface points. The weight is proportional to the infinitesimal radial volume interested in the tire-ground contact. This infinitesimal radial volume is calculated as the intersection of the contact volume with the line $\u2113(x,e\u0302y)$ passing through the point $x=[x,\u2009y,\u2009z]T$ and matching the axes of the wheel $e\u0302y$ orthogonally. If $proj(x,\u2202G)$ defines the intersection point of the line $\u2113(x,e\u0302y)$ with $\u2202G$, then

**P**is calculated as

*λ*is calculated as

The solution of the here defined integrals (5)–(7) is not straightforward. This is why, in Sec. 3, we propose a numerical approximation of these integrals.

**P**and axes directions are

The *β _{x}* and

*β*angles, which are used to describe the local ground surface orientation with respect to the wheel hub, are calculated as the Euler angles (

_{y}*zxy*sequence) between the hub reference frame $Hxyz$ and the contact point reference frame $Pxyz$ (see Fig. 3).

## 3 Algorithm Implementation

So far, we have not yet described the actual shapes of the sets $T$ and $G$, which represent the tire and the ground sets, respectively. Since we will deal with complicated obstacle shapes, it is important to restate the presented enveloping model in a formulation that does not limit its numerical efficiency, robustness, and scalability. For this reason, a discretized representation of tire and ground elements is required both to eventually solve the contact problem and to offer a highly scalable framework to work with. More specifically, the tire's external shape roughly corresponds to the external tread pattern profile revolved around the hub axis $e\u0302y$. The tire set is discretized into a series of laterally lumped *disks*, also called *ribs*. On the other hand, the ground boundary is represented by a triangular mesh (see Fig. 4). Both tire and ground descriptions are fully scalable as the density of the ribs and the triangles are adjusted according to the accuracy and execution speed needed for the specific simulation. This representation approach turns out to be useful in the case of hard real-time simulation (where the execution speed represents an important requirement that needs to be satisfied) but also in the case of offline simulations (where the demanded precision is usually higher).

*n*tire ribs

*ϕ*, i.e., $S={\varphi 1,\u2009\varphi 2,\u2009\u2026\u2009\varphi n\u22121,\u2009\varphi n}$. The generic tire rib

*ϕ*is defined geometrically as a disk of radius

*r*, whose center is

**O**, and whose face normal corresponds to $e\u0302y$. In addition, the rib also retains its “virtual” width $w\u2208\mathbb{R}\u22650$, which will be used in the numerical integration processes

*m*triangles

*τ*, i.e., $M={\tau 1,\u2009\tau 2,\u2009\u2026,\u2009\tau m\u22121,\u2009\tau m}$, where

and where $p1,\u2009p2$, and $p3$ are the three vertices of the triangle. Notice that the mesh triangle *τ* also carries a friction coefficient scaling factor $\lambda \u2208\mathbb{R}\u22650$, which is considered constant over the triangle surface.

The intersection between the set of ribs $S$ and the ground mesh triangles $M$ allows us to identify the contact parameters, namely, the contact patch area *A*, the intersection volume *V*, the contact point location **P**, and the contact force direction $n\u0302$. Trivially, if a generic rib *ϕ* is intersected with a generic mesh triangle *τ*, the result will be a set $\sigma =\varphi \u2229\tau $. This intersection can lead to four different cases, in which the set $\sigma $ can be

an empty set, if

*ϕ*does not touch*τ*at all;a point, if

*ϕ*intersects*τ*one of its vertices;a segment, if

*ϕ*intersects*τ*or a portion of it;a convex hull, if

*ϕ*and*τ*are coplanar and intersect.

*A*by restating (3) as the summation

where *w _{i}* is the virtual width of the

*i*th rib.

**O**to a generic point of the intersection segment $\sigma \xaf$. Denote with

*r*the radius of the

_{i}*i*th rib and with $r\xafij(\theta )$ the distance between the

*i*th rib center point and the segment $\sigma \xafij$; where the segment $\sigma \xafij$ is identified as the intersection of the

*i*th rib with the

*j*th triangle. The range of the angle

*θ*for the generic segment $\sigma \xaf$ is denoted with $[\theta a,\theta b]\u2282(\pi ,2\pi )$, which are implicitly defined as

Notice that the intermediate variable $vij(\theta )$ in Eq. (11) represents the length of the generic segment radially connecting the $\sigma \xafij$ to the external tire boundary at the angle *θ* (see Fig. 5). From a physical perspective, $vij(\theta )$ is the deflection of the infinitesimal radial springs at the angle *θ*. It is possible to find a closed-form solution for $r\xafij(\theta )$. However, the analytical expression of the resulting **P** and $n\u0302$ is excessively long. For practical purposes, it is better to use a quadrature formula to numerically approximate the integrals. In this work, Simpson's 1/3 rule is used.

*Remark 2*.

*Simpson's 1/3 rule quadrature formula of*

*f*(

*x*)

*for the interval*$[a,b]$

*reads as*

where $h=b\u2212a,\u2009c=(a+b)/2$, and $\xi \u2208[a,\u2009b]$ [22]. Notice that if $r\u22641,\u2009\theta b\u2212\theta a\u2264\pi /6$, and $1\u2265||O\u2212pa,b||\u22654/5$, the relative error of the numerically computed integrals in Eqs. (12)–(14) is below 1%.

*Remark 3*.

*The midpoint*$r\xaf((\theta a+\theta b)/2)$

*is computed as*

## 4 Software Architecture

The Enve algorithm-based software is written in C++ (2011 standard [23]), which is one of the most widely supported, common and fast among object-oriented general-purpose programing languages. We have only reduced the dependencies on the Eigen3 [24] and Acme [18] libraries: this provides us with a flexible and extensible framework. Eigen3 is a template library for linear algebra. It is very efficient for small vectors and matrices and can exploit LAPACK/BLAS [25] for peak performance when matrices and vectors have a large size. Acme, on the other hand, is an easy-to-use geometry library that aims to be a minimal tool for hard real-time scheduling applications. Furthermore, Matlab^{®} Mex and Simulink^{®} S-Function extensions are provided to allow the integration of Enve in more complex simulation environments. The software is distributed under the BSD-3-Clause License and is freely available online [26].

### 4.1 Data Types.

The software is built on a limited number of classes. We have chosen to keep the library as essential as possible for ease of maintenance and efficiency. For this reason, we use the necessary classes to describe and manipulate virtual ground surfaces and tires. Each geometrical entity representing the tire and ground is organized in a specific class. A brief description of the data types used in the software is now given in this section.

#### 4.1.1 TriangleGround.

The generic ground triangle *τ* is defined by the three points $p1,\u2009p2$ and $p3$. In this way, it geometrically corresponds to the Triangle object described in the Acme library [18]. For this reason, the TriangleGround object publicly inherits from the Acme::Triangle object. Since the TriangleGround object describes the generic triangle representing a local road region, this class also carries a scaling factor for the friction coefficient *λ*. The mathematical representation corresponds to the definition in Eq. (9).

#### 4.1.2 Mesh.

The mesh represents the discretized ground set $M$ on which the tire rolls. It is composed of a multitude of TriangleGround objects. The size of these triangles is determined according to end-use requirements and desired accuracy. As a consequence, the number of triangles can range from a few dozen to tens of thousands. Placing a tire that measures a few tens of centimeters on a road surface of tens of square kilometers poses the delicate challenge of how to perform a computationally efficient intersection between these two elements. When numerous geometric objects are intersected simultaneously, carrying out a sequential analysis requires a lot of time. It is therefore useful to limit intersection operations to those objects that could potentially intersect. For this reason, the approach of the geometry library Acme is based on the construction of a spatial index, also known as a Bounding Volume Hierarchy (BVH), of the type Axis-Aligned Bounding Box (AABB) tree. This data structure makes it possible to partition the space covered by the set of geometric entities hierarchically. The idea behind this type of BVH is to determine simple Bounding Volume Primitives (BVPs) at the different nodes of the hierarchy, roughly covering and fitting the set of geometric primitives contained within the corresponding subtree. The computation of the solution to any collision or intersection problem is accelerated thanks to an a priori Acme's AABB-AABBtree or AABBtree-AABBtree intersection algorithm, which restricts the intersecting candidates only to the candidates contained within the colliding tree leaves [18].

#### 4.1.3 Flat.

*λ*is added to the class description

where **p** is a generic point on the plane and $n\u0302$ is a unit normal vector.

#### 4.1.4 Shape.

The Shape object describes the overall external shape of the tire $T$ and of all the ribs in the set $S$. The external boundary of the tire $\u2202T$ is defined as a surface generated by the revolution of an arbitrary function *R*(*y*). By using the proper surface of revolution, it is possible to represent a wide variety of tire morphologies, ranging from car tires to motorcycle tires (see Fig. 1).

#### 4.1.5 Rib.

As mentioned before, the rib is considered a small section of the overall tire. Geometrically the rib is defined as a disk *ϕ* of radius *r*, whose center **O** is located in $[0,\u2009y,\u20090]T$ in the $Hxyz$ wheel hub reference frame and whose lateral coordinate *y* is constant throughout the whole disk area. The Rib object publicly inherit from the Acme::Disk object. In addition, the Rib class also carries its virtual width $w\u2208\mathbb{R}\u22650$, which is stored as a class member. The mathematical representation of this set corresponds to that in Eq. (8).

#### 4.1.6 Shell.

The shell object encapsulates the set of Rib objects $S$, along with a shared pointer to a Shape object and an object of the type Acme::AABB. This last object represents the bounding box that contains the tire at every instant of time. The Shell object also contains all the attributes and secondary objects that are used to describe the tire and its intersection with the terrain. Whenever an intersection between the tire and the terrain occurs, the affine transformation describing the tire pose as well as the tire AABB is updated. Once these basic operations are completed, the bounding box is intersected with the Acme::AABBtree of the mesh. A list of shared pointers to TriangleGround objects that are internal or that simply touch the tire bounding box is generated. This list is supplied one by one to the Rib objects to find the intersection parameters previously described in Sec. 3. To avoid un-necessary recalculations, the intersection parameters are stored in members that are internal to the Shell class; they are then eventually extracted with methods specifically designed for this purpose.

### 4.2 Algorithm Workflow.

The workflow of the algorithm is highly streamlined. It is divided into three main algorithms that are briefly described below briefly described and that are also reported in the pseudo-codes. Other minor algorithms are used to support the main ones and are not reported here for the sake of brevity.

#### 4.2.1 Algorithm 1.

The initialization of the Shell object is done by passing a Shape object and the number of ribs *n* as input parameters. First, the Shape object is used to reconstruct the tire's external shape; then the Rib objects are instantiated by discretizing the tire into *n* ribs.

Require: Shape: The tire Shape object |

n: The number of ribs |

procedure Shell(Shape$,\u2009\u2009n$) ▷ Shell class constructor |

$S\u2190\u2205$ ▷ Initialize the set of ribs |

$w\u2190\u2009Shape.width()/n$ ▷ The ribs width |

for$i\u21901$ to ndo |

$yi\u2190w/2+(i\u22121)\u2009w$ ▷ The i-th rib lateral coordinate |

$ri\u2190\u2009Shape.radius(yi)$ ▷ The i-th rib radius |

S_{i} ← Rib$(ri,\u2009yi,\u2009w)$ ▷ The i-th rib |

end for |

$F\u2190I4\xd74$ ▷ The default affine transformation |

$B\u2190AABB(S,F)$ ▷ The bounding box of the tire |

store$S,\u2009B,\u2009F$ ▷ Store the data in class members |

end procedure |

Require: Shape: The tire Shape object |

n: The number of ribs |

procedure Shell(Shape$,\u2009\u2009n$) ▷ Shell class constructor |

$S\u2190\u2205$ ▷ Initialize the set of ribs |

$w\u2190\u2009Shape.width()/n$ ▷ The ribs width |

for$i\u21901$ to ndo |

$yi\u2190w/2+(i\u22121)\u2009w$ ▷ The i-th rib lateral coordinate |

$ri\u2190\u2009Shape.radius(yi)$ ▷ The i-th rib radius |

S_{i} ← Rib$(ri,\u2009yi,\u2009w)$ ▷ The i-th rib |

end for |

$F\u2190I4\xd74$ ▷ The default affine transformation |

$B\u2190AABB(S,F)$ ▷ The bounding box of the tire |

store$S,\u2009B,\u2009F$ ▷ Store the data in class members |

end procedure |

#### 4.2.2 Algorithm 2.

The Mesh object is initialized by passing a Road Data Format or Wavefront OBJ extension file containing the triangles' vertices and faces as input. The file is first parsed to instantiate the mesh TriangleGround objects and then used to build the AABBtree, which will later be used to accelerate the intersection process as well as the calculation of contact parameters.

Require:$file$: The mesh file path |

procedure Mesh($file$) ▷ Mesh class constructor |

$V\u2190\xf8$ ▷ Initialize ground mesh vertices |

$F\u2190\xf8$ ▷ Initialize mesh faces |

$\lambda \u2190\xf8$ ▷ Initialize mesh λs |

$V,\u2009F,\u2009\lambda \u2190\u2009Mesh.parse(file)$ ▷ Parse the mesh file |

$M\u2190\xf8$ ▷ Initialize the set of triangles |

$B\u2190\xf8$ ▷ Initialize the set of triangles' AABBs |

$n\u2190F.size()$ ▷ The number of triangles |

for$i\u21901$ to ndo |

$p1\u2190V(F(i,1))$ ▷ The 1st triangle vertex |

$p2\u2190V(F(i,2))$ ▷ The 2nd triangle vertex |

$p3\u2190V(F(i,3))$ ▷ The 3rd triangle vertex |

$Mi\u2190TriangleGround(p1,\u2009p2,\u2009p3,\u2009\lambda (i))$ |

$Bi\u2190AABB(Mi)$ ▷ The i-th triangle AABB |

end for |

$T\u2190AABBtree(M,\u2009B)$ ▷ The mesh's AABBtree |

store$M,\u2009T$ ▷ Store the data in class members |

end procedure |

Require:$file$: The mesh file path |

procedure Mesh($file$) ▷ Mesh class constructor |

$V\u2190\xf8$ ▷ Initialize ground mesh vertices |

$F\u2190\xf8$ ▷ Initialize mesh faces |

$\lambda \u2190\xf8$ ▷ Initialize mesh λs |

$V,\u2009F,\u2009\lambda \u2190\u2009Mesh.parse(file)$ ▷ Parse the mesh file |

$M\u2190\xf8$ ▷ Initialize the set of triangles |

$B\u2190\xf8$ ▷ Initialize the set of triangles' AABBs |

$n\u2190F.size()$ ▷ The number of triangles |

for$i\u21901$ to ndo |

$p1\u2190V(F(i,1))$ ▷ The 1st triangle vertex |

$p2\u2190V(F(i,2))$ ▷ The 2nd triangle vertex |

$p3\u2190V(F(i,3))$ ▷ The 3rd triangle vertex |

$Mi\u2190TriangleGround(p1,\u2009p2,\u2009p3,\u2009\lambda (i))$ |

$Bi\u2190AABB(Mi)$ ▷ The i-th triangle AABB |

end for |

$T\u2190AABBtree(M,\u2009B)$ ▷ The mesh's AABBtree |

store$M,\u2009T$ ▷ Store the data in class members |

end procedure |

#### 4.2.3 Algorithm 3.

The contact parameters are calculated by intersecting the Shell object and the Mesh object. This is done by calling the “setup” method of the Shell object. First, this method intersects the mesh's AABBtree with the tire's AABB and returns a list of triangles that are potentially intersecting the tire. Then, the tire ribs are intersected with the triangles and the contact parameters are calculated as described in Sec. 2.

Require:F: The new tire affine transformation |

$M$: The Mesh object |

procedure Shell.setup($F,\u2009M$) ▷ The intersection method |

B ← Shell$.update_aabb(F)$ ▷ Update the tire AABB |

L ← Mesh$.intersect(B)$ ▷ The intersecting triangles |

$r\u2190S.size()$ ▷ The number of tire ribs |

$n\u2190L.size()$ ▷ The number of intersecting triangles |

$A\u21900$ ▷ Initialize the contact patch area |

$V\u21900$ ▷ Initialize the intersection volume |

$P\u21900$ ▷ Initialize the contact point |

$n\u0302\u21900$ ▷ Initialize the contact normal |

$\lambda \u21900$ ▷ Initialize the friction coefficient scaling factor |

for$i\u21901$ to rdo |

for$j\u21901$ to ndo |

$\sigma \u2190Si\u2229ML(j)$ ▷ The intersection set (segment) |

$A\u2190A+\sigma .area()$ ▷ The contact area (10) |

$V\u2190V+\sigma .volume()$ ▷ The contact volume (12) |

$P\u2190P+\sigma .point()$ ▷ The contact point (13) |

$n\u0302\u2190n\u0302+\sigma .normal()$ ▷ The contact normal (14) |

$\lambda \u2190\lambda +\sigma .lambda()$ ▷ The contact λ(15) |

end for |

end for |

$\lambda \u2190\lambda /V$ ▷ The average λ |

return$A,\u2009V,\u2009P,\u2009n\u0302,\u2009\lambda $ ▷ The contact parameters |

end procedure |

Require:F: The new tire affine transformation |

$M$: The Mesh object |

procedure Shell.setup($F,\u2009M$) ▷ The intersection method |

B ← Shell$.update_aabb(F)$ ▷ Update the tire AABB |

L ← Mesh$.intersect(B)$ ▷ The intersecting triangles |

$r\u2190S.size()$ ▷ The number of tire ribs |

$n\u2190L.size()$ ▷ The number of intersecting triangles |

$A\u21900$ ▷ Initialize the contact patch area |

$V\u21900$ ▷ Initialize the intersection volume |

$P\u21900$ ▷ Initialize the contact point |

$n\u0302\u21900$ ▷ Initialize the contact normal |

$\lambda \u21900$ ▷ Initialize the friction coefficient scaling factor |

for$i\u21901$ to rdo |

for$j\u21901$ to ndo |

$\sigma \u2190Si\u2229ML(j)$ ▷ The intersection set (segment) |

$A\u2190A+\sigma .area()$ ▷ The contact area (10) |

$V\u2190V+\sigma .volume()$ ▷ The contact volume (12) |

$P\u2190P+\sigma .point()$ ▷ The contact point (13) |

$n\u0302\u2190n\u0302+\sigma .normal()$ ▷ The contact normal (14) |

$\lambda \u2190\lambda +\sigma .lambda()$ ▷ The contact λ(15) |

end for |

end for |

$\lambda \u2190\lambda /V$ ▷ The average λ |

return$A,\u2009V,\u2009P,\u2009n\u0302,\u2009\lambda $ ▷ The contact parameters |

end procedure |

## 5 Simulations and Discussion

and where *R _{x}* = 0.313 m,

*R*= 0.11 m,

_{y}*m*= 9,

_{x}*m*= 6, and

_{y}*L*= 0.1025 m.

_{y}### 5.1 Quasi-Static Simulations on Uneven Road Surfaces.

As mentioned before the aim of the first set of simulations is to assess the enveloping capabilities of the presented model at a low speed according to SAE J2731 standard [27]. As motivated in the SAE standard, we want to assess under quasi-static considerations the ability of the model to correctly estimate the contact point height *z*^{2} and the forward and lateral slope angles *β _{x}* and

*β*, respectively. In this regard, we have considered three different road surfaces, characterized by different degrees of unevenness. On each of these surfaces, we have performed a set of simulations with constant vertical loads

_{y}*F*of 2, 4, and 6 kN. It must be pointed out that the TMEasy model output does not change with the vertical load

_{z}*F*and is thus reported with a single line in the plots. Notice that, in these tests, according to the SAE J2731 standard, the wheel hub is only allowed to move along the

_{z}*x*-axis, to advance the tire over the road surface, and along the

*z*-axis, to keep the tire in contact with the road surface.

The first simulation is performed on a series of different obstacles positioned orthogonally to the forward direction of the tire. Specifically, the obstacles' cross section are depicted at the top of Fig. 6. The results of the simulation are reported in the last two plots of Fig. 6. As it can be noticed, the behavior of the presented model is similar to that of the Swift® model, which is considered to be very close to the ground truth (see Ref. [5] for an in-depth discussion on the Swift® model and its accuracy). The performance of Enve depends on the geometry of the obstacles. Indeed, in the presented model the contact parameters are defined only by the geometry of the obstacles, which is not entirely true in the real world. Tire nonlinear radial stiffness and inflation pressure also play a non-negligible role in defining the overall tire-ground contact behavior. However, given the simplicity of the model, the results are quite satisfactory. Conversely, the TMEasy model does not produce realistic results, causing abrupt changes in the contact parameters: *z*, *β _{x}*, and

*β*.

_{y}Although in the previous test TMEasy method is shown not to be particularly accurate, it is certainly an easy-to-implement model, suitable for usage on sufficiently “smooth” road surfaces. To demonstrate its capabilities in estimating the contact point position and the forward and lateral slope angles on simple road surfaces, we perform a second simulation on a road surface defined by the map $z(x,y)=0.025\u2009\u2009sin(7\u2009(x\u2212tan(\pi /18)\u2009y)2)$. The results reported in Fig. 7 show that the TMEasy model does not produce abrupt changes in the contact parameters anymore, but as the road surface wavelength decreases, it fails to estimate the enveloping behavior of the Swift® model. Moreover, the lack of dependency on the vertical load *F _{z}* exacerbates this gap of the TMEasy model behavior. On the other hand, the Enve model results are close to ground truth.

The last simulation is performed on a rough Belgian block road surface depicted in Fig. 1. The results reported in Fig. 8, show that the Enve model behavior is very close to that of the Swift® model in terms of contact point height *z* and forward slope angle *β _{y}*. However, the Enve model is not able to correctly estimate the lateral slope angle

*β*. Even on the rough Belgian block road surface, the TMEasy model is unable to correctly estimate the contact point height

_{x}*z*and the lateral slope angle

*β*. Moreover, it strongly overestimates the forward slope angle

_{x}*β*.

_{y}### 5.2 Full-Vehicle Model Simulation.

To evaluate the behavior of the presented tire-ground enveloping model in a typical use case, it is integrated into a complete real-time vehicle simulation framework, and compared with the Swift® [5] and TMEasy [11,13] enveloping models. The used framework consists of a DIL/SIL/HIL simulator based on a custom high-fidelity vehicle model. Vehicle dynamics is described by a 14 Degree of Freedoms (DOFs) full-vehicle model able to accurately simulate the impact of suspension kinematics and compliance on vehicle behavior within the typical range of frequencies for the ride and handling analysis. Contrary to widespread practice, the vehicle dynamic equations are not linearized. This improves model accuracy in nonlinear regions and vehicle attitude even in abnormal conditions. The described formulation enables the vehicle model for 3D road simulations. It is worth noting that additional DOFs are used to describe internal components' dynamics, such as the powertrain, the steering system or the tire subsystems. The structure of the model is organized in a modular fashion mirroring the mechanical interfaces of vehicle components. This modular structure makes the integration of third-party software or, more generally, of any external system a fairly straightforward process. As a result, it is possible to easily configure the vehicle model for HIL, SIL and/or DIL simulations. The developed enveloping model is thus interfaced, as a subsystem with the vehicle MB system through a predefined model interface. Inputs of the tire-ground contact subsystem are the wheels' reference frames as well as the wheels' geometrical characteristics and the 3D road surface mesh. The outputs of the subsystem are, for each wheel, the equivalent tire contact point, the local contact plane, the contact penetration and the local friction coefficient. Road contact parameters, extracted from the tire-ground enveloping model, are then used by the MF 6.2 tire model to compute contact forces.

The modeled vehicle is a passenger car with sedan-like characteristics. The vehicle features a fully electric powertrain with a front-wheel drive transmission layout. The modeled suspensions are a MacPherson configuration for the front suspensions and a Multi-Link scheme for the rear suspensions. For the sake of brevity, only the most relevant vehicle parameters are listed in Table 1.

Parameter description | Value |
---|---|

Total mass of the vehicle | 1300 kg |

Center-of-mass height | 0.30 m |

Front axle distance from the center-of-mass | 1.25 m |

Rear axle distance from the center-of-mass | 1.45 m |

Wheelbase | 2.70 m |

Yaw inertia | 1400 kg m^{2} |

Track width | 1.50 m |

Steering ratio | 20 |

Maximum torque at wheel | 1200 N m |

Maximum motor power | 150 kW |

Wheel inertia around rotation axis | 1.40 kg m^{2} |

Size specification of the tires | 205/60R15 |

Spring stiffness at wheel | 3530 kN/m |

Damping coefficients at wheel for jounce | 789 N s/m |

Damping coefficients at wheel for rebound | 1578 N s/m |

Parameter description | Value |
---|---|

Total mass of the vehicle | 1300 kg |

Center-of-mass height | 0.30 m |

Front axle distance from the center-of-mass | 1.25 m |

Rear axle distance from the center-of-mass | 1.45 m |

Wheelbase | 2.70 m |

Yaw inertia | 1400 kg m^{2} |

Track width | 1.50 m |

Steering ratio | 20 |

Maximum torque at wheel | 1200 N m |

Maximum motor power | 150 kW |

Wheel inertia around rotation axis | 1.40 kg m^{2} |

Size specification of the tires | 205/60R15 |

Spring stiffness at wheel | 3530 kN/m |

Damping coefficients at wheel for jounce | 789 N s/m |

Damping coefficients at wheel for rebound | 1578 N s/m |

Model accuracy and robustness are widely tested with DIL simulations by driving the vehicle in urban scenarios, characterized by speed bumps, sidewalks or, generally, uneven road surfaces. In this work, we present the results obtained by driving the vehicle at a longitudinal speed of 10 km/h over a 45° oblique and 1 cm high step. The reported results are the chassis accelerations (Fig. 9), contact parameters (Fig. 10), front-left tire forces (Fig. 11), and the tires deflections (Fig. 12). Note that Enve and Swift® enveloping models present a similar behavior even in the case of a simulation of a full-vehicle model. In fact, the generated chassis accelerations, the tires' forces, and deflections are in agreement with each other. The discontinuities in both the contact point height and the banking and forward slope angles of the TMEasy model reflect on the full-vehicle model simulations with abrupt changes in contact forces as well as chassis acceleration not seen in experimental measurements. Notice that according to the results reported in Figs. 7 and 8, the Enve model is able to correctly estimate the contact point height, but it slightly underestimates the banking angle. This causes the correct prediction of the vertical contact force magnitude, but a slight underestimation of the lateral contact force magnitude. This error may be emphasized or reduced depending on the tire's properties.

### 5.3 Real-Time Performance.

To access the computations' performance of the Enve C++ software implementation, we will make use of the Real-Time Factor (RTF) metric. Specifically, the RTF is defined as the ratio between the time needed to process the input and the input duration. In order to consider a system a real-time system, RTF should be $\u2264$1. We performed a test in which the triangles inside the AABB of the tire are increased from 1 up to 10^{4}, while the number of ribs is increased from 1 to 10. The simulation platform is an iHawk™ Concurrent Real-Time provided with 2.5 GHz Intel Xeon Silver 4215 8 Core, 11 MB cache, 32 GB DDR4 RAM, and 64 bit RedHawk Linux RTOS with CentoOS distribution. All tests are performed on a shielded CPU. The results reported in Fig. 13 show that the RTF of Enve is mostly impacted by the number of triangles inside the AABB of the tire rather than the number of ribs. Overall, the RTF of the software is $\u2264$1 for a number of triangles inside the AABB of the tire up to $\u2248$1000, which corresponds to a regular grid discretization of $\u2248$1.5 cm spacing. However, in a typical vehicle simulation, only a part of the time-step is used to perform the tire-ground contact calculations. The available time depends on several factors, such as the complexity of both the vehicle and the tire models. To the authors' knowledge, a good practice is to use roughly 25% of the integration time-step to establish the contact parameters. This means that the RTF of the software should be $\u2264$0.25 to be considered real-time. In this case, the maximum number of triangles inside the tire AABB drops to $\u2248$200, which corresponds to a regular grid discretization of $\u2248$3.5 cm spacing. Overall, these performances are more than sufficient for most real-time applications. Further improvements can be achieved by parallelizing the AABBtree-AABB and the Rib-TriangleGround intersection algorithms, which are currently implemented in a serial fashion. Moreover, different instances of Enve objects can run concurrently in separated CPUs to further improve the overall performance of the full-vehicle model.

## 6 Conclusions

In this paper, we presented a novel approach for tire-ground enveloping modeling. This approach has been implemented in an C++ library called Enve, which has also been embedded in an MBS code to assess the capabilities of the newly introduced algorithm-based software. A thorough comparison between Enve, Swift® and TMEasy has been performed for both quasi-static regime tests and full-vehicle model simulations. The results show that the proposed approach is concurrent with existing state-of-the-art methods. Particularly, Swift® and Enve show similar behaviors, conversely, TMEasy leads to nonphysical behavior of the vehicle model in case of sharp or short wavelength road surface unevennesses. Besides the ability to accurately model the tire enveloping behavior, Enve has also proven to be as fast as real-time or even faster than real-time capable depending on the number of triangles inside the tire AABB. This makes Enve a suitable candidate for real-time applications.

Similarly to TMEasy, Enve is a contact method that can be used for different tire models, e.g., custom physical tire models as well as the MF tire model. Swift®, on the other hand, has been developed to be used in tandem with the MF tire model. Furthermore, Enve is fully parameterizable by the tire's external shape alone. This results in a greater immediacy of use. In addition to the contact point position and the banking and slope angles, it can also provide information about contact patch geometry and local ground properties, such as the effective scaling factor of the friction coefficient. Due to its scalability, the Enve library is also ready to be used in structured MB simulation environments. In the future, thanks to its extensible structure it might be possible to include additional features, which may help with multiple-contact points tire modeling, as well as thermo-mechanical, hydroplaning and road surface noise analyses. This may be carried out by exploiting the contact information on each rib of the discretized tire provided by the Enve library.

## Data Availability Statement

The data and information supporting the findings of this article are freely available online under the open-source BSD 3-Clause License at the link https://github.com/StoccoDavide/enve.

## Footnotes

The contact point height is defined as the *z*-axis component of the contact point **P** in the absolute reference frame. To transform **P** in the absolute reference one might need to exploit the reference frames $Pxyz$ and $Hxyz$, depending on which reference frames **P** is computed.