# A Brush-Type Tire Model with Nonsmooth Representation

Authors: Ryo Kikuuwe
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper proposes a brush-type tire model with a new mathematical representation. The presented model can be seen as a generic model that describes the distributed viscoelastic force and Coulomb-like friction force, which are balancing each other at each point, in the contact patch. The model is described as a partial differential algebraic inclusion (PDAI), which involves the set-valuedness to represent the static friction. A numerical integration algorithm for this PDAI is derived through the implicit Euler discretization along both space and time. Some numerical comparisons with Magic Formula and a LuGre-based tire model are presented. The results show that, with appropriate choice of parameters, the proposed model is capable of producing steady-state characteristics similar to those of Magic Formula. It is also shown that the proposed model realizes a proper static friction state, which is not realized with a LuGre-based tire model.


---

## Body

## 1. Introduction

Tire-road interaction plays a crucial role in the dynamic behaviors of wheeled machinery, such as automobiles and wheeled mobile robots. It is a complicated physical process that includes the viscoelastic deformation and the mixture of static friction and kinetic friction in a single contact patch. Better models of such a physical interaction will contribute to many applications, such as computer-aided manufacturing of automobiles, physics-based animation and gaming, and control of wheeled mobile robots and autopilot vehicles.One main factor of tire-road interaction is the distributed Coulomb-like friction, which inherently includes the nonsmooth (discontinuous) relation between the friction force and the relative velocity. In the kinetic friction state, the force is determined by the relative velocity, while in the static friction state, the relative velocity maintains zero as long as the friction force is below the maximum static friction force. That is, the causality between the relative velocity and the friction force is reversed between the two states. Moreover, these two states coexist in the single contact patch; even while a portion of the contact patch is in the static friction state, the other portions may not be. Such features pose difficulty in the construction of a physics-based tire model.One typical approach to the modeling of tires is called a brush-type tire model [1–3]. In this model, the rigid wheel is in contact with the road through brush-like elastic bristles, which represent the deformable tire surface. The Coulomb-like friction takes place between the bristle tips and the road. The most impactful assumption here is that the bristles do not interact with one another. This assumption allows for the simplicity of the model but contradicts the fact that the tire surface is a continuum, which produces compressive and tensile stress in the tangential direction. Moreover, especially in the pneumatic tire of automobiles, there exists an elastic layer called carcass, between the rigid wheel and the contact surface, which is neglected in the brush-type model. Despite this flaw, the brush-type model has been recognized as a candidate for a tire model [1–4] and also as a basic framework to which modifications can be made [3, 4]. One can also say that the brush-type model may be valid for nonpneumatic tires, such as those used for mobile robots.Some mathematical representations of brush-type tire models have been developed based on the friction model called the LuGre model. The LuGre model [5] is originally a one-dimensional, single-state dynamic friction model that has been developed for system control purposes. The main advantage of the LuGre friction model is that it is described as a smooth ordinary differential equation (ODE) involving no discontinuities. Canudas-de-Wit et al. [6] extended the LuGre friction model for the modeling of the longitudinal tire friction force, and Velenis et al. [7] further extended it into a two-dimensional model. Deur et al. [1, 8] also constructed two-dimensional LuGre-based tire models, which differ from Velenis et al.’s [7] model in the treatment of anisotropy and in the way of two-dimensional extension. One drawback of the original, one-dimensional LuGre friction model is that it does not reproduce the exact static friction, as pointed out by many previous researchers [9, 10]. This flaw is rooted in the smoothness of the model, neglecting the nonsmooth nature of real friction, and is inherited by the LuGre-based tire models.This paper presents a new brush-type tire model described as a partial differential algebraic inclusion (PDAI). The new model is an extension of the nonsmooth friction model presented by the author and his colleagues [10], which is described as a differential algebraic inclusion (DAI) that combines a set-valued representation of the Coulomb friction and the linear viscoelasticity. This paper extends this DAI-based friction model into a new PDAI-based tire model, which involves distributed viscoelasticity and Coulomb-like friction. It should be noted that this paper does not argue the superiority of the brush-type tire model over different types of sophisticated tire models [11–13]. Instead, this paper presents a new generic representation of a rolling object with a viscoelastic surface. Although its accuracy may not be high with respect to pneumatic tires, it may be usable as a model of nonpneumatic tires, such as those of relatively stiff rubber, and also as a basic model upon which some modifications can be made. One strength of the proposed approach is that it realizes the static friction state, which is not realized by LuGre-based tire models. Another notable point may be that it does not involve singularity at the zero velocity [14], which should be properly handled in the application of Pacejka’s Magic Formula [15, 16].The remainder of this paper is organized as follows. Section2 discusses the previous work on general friction models and tire friction models. Section 3 presents the new mathematical representation of the brush-type tire model and its spatio-temporal integration scheme based on the implicit Euler discretization. Section 4 presents some illustrative numerical examples, including comparison with Magic Formula [15, 16] and a LuGre-based tire model. Section 5 provides concluding remarks.


## 2. Preliminaries

Throughout this paper,ℝ denotes the set of all real numbers and ℝ+ denotes the set of strictly positive real numbers. The statements ϕ:X⟶Y and ϕ:X⇉Y mean that ϕ is a singled-valued function and a set-valued function, respectively, from the set X to the set Y. The symbol ∂xfx denotes the generalized gradient [17] of the function f with respect to x. We also let 0 be the zero or the zero vector of an appropriate dimension.We let the set-valued functionssgn:ℝ⇉ℝ and SgnX:ℝ2⇉ℝ2 be defined by(1)sgnx=Δxx,ifx≠0,−1,1,ifx=0,SgnXx=ΔX2xXx,ifx≠0,Xe∈ℝ2e≤1,ifx=0,respectively, where X is a 2×2 positive definite matrix. Note that (2)∂xXx=SgnXx.We also write Sgnx=ΔSgnIx, where I is the 2×2 identity matrix.


### 2.1. One-Dimensional Friction Models

Many friction models have ever been proposed for control and simulation purposes [18, 19]. Among them, the LuGre friction model [5] is a model that has been extended into some tire models. The original LuGre friction model can be described as follows:(3)e˙=v−KveΨv,f=Ke+Be˙.Here, v∈ℝ is the input velocity, f∈ℝ is the output friction force, e∈ℝ is a state variable that is interpreted as the average deflection of contact bristles, and K and B are positive constants representing the stiffness and damping coefficients, respectively, of the contact bristles. The function Ψ:ℝ⇉ℝ is a function representing a velocity-dependent friction law, which is typically given as follows:(4)Ψv=ΔsgnvFs−Fs−Fks1vVs,γ,where(5)s1u,γ=Δ1−exp−uγ,and Fs, Fk, Vs, and γ are positive constants. The constants Fs and Fk can be interpreted as the maximum static friction force and the kinetic friction force, respectively, which usually satisfy Fs≥Fk. With this function Ψ, the magnitude of the friction force decreases from Fs to Fk as the velocity increases along a curve characterized by Vs and γ. The function Ψ defined above is set-valued but the right-hand side of the first equation of (3) is always single-valued and continuous because it is zero when v=0. That is, the LuGre friction model (3) can be seen as an ODE, which is smooth and mathematically convenient.The form (3) of the model is inherited from the Dahl model [20], which has been chosen so that(6)f≈Ψv≈Keis achieved when the input velocity v maintains a constant nonzero value. A physical interpretation of the model (3) is often represented as in Figure 1. In this model, there is an elastic deflection e in the frictional interface and the friction force f is obtained as the viscoelastic force determined by the second equation of (3).Figure 1
A physical model upon which the LuGre friction model [5] and Kikuuwe et al.’s [10] friction model are built.Because the LuGre friction model neglects the nonsmooth nature of friction, it does not make the distinction of two different states: the static friction state and the kinetic friction state. In this model, the state variablee can be understood as the elastic deflection of the bristle, but its rate-of-change is not equal to the input velocity v. In other words, the velocity of the bristle tip cannot be the exact zero, and thus the static friction state cannot be realized. As a consequence, the object’s displacement is not purely elastic even in the presliding regime and unbounded drift can happen even under external forces smaller than the maximum static friction force. This point has been formally pointed out by Theorem 1 of Dupont et al. [9]. Some friction models to circumvent this problem have been proposed so far [21–25], with special attention on reproducing hysteresis behavior in the elastic displacement in the static friction state.In contrast to the smooth LuGre model, Kikuuwe et al.’s [10] friction model is a nonsmooth model. Their friction model is described as the following DAI:(7)f∈Ψv−e˙,f=Ke+Be˙.This can be viewed as a pair of simultaneous differential equations that determines f and e˙ so that both equations are satisfied (Mathematical expressions like (7), which involves the symbol “∈” and derivatives, are referred to as differential inclusions (DIs) [26] instead of differential equations. Moreover, because the derivative e˙ in (7) cannot be isolated in the left-hand side, expression (7) can also be referred to as a differential algebraic inclusion (DAI)). This model can also be interpreted as in Figure 1 and can be said to be more faithful to Figure 1 than the LuGre friction model because it takes into account the exact balance between the friction force, f in the first equation of (7), and the viscoelastic force, f in the second equation of (7). Kikuuwe et al. [10] presented an approach to implement this model with the backward Euler discretization, and Xiong et al. [25, 27] extended it into the continuous-time domain. By setting B=0, Kikuuwe et al.’s [10] friction model reduces to Hayward and Armstrong’s [28] friction model, which is defined only in the discrete-time domain.


### 2.2. LuGre-Based Brush-Type Tire Models

As overviewed in Section1, many of the brush-type tire models are based on the LuGre friction model (3). In their approach, the input velocity v is the velocity of the rigid wheel (or of the deformed carcass centerline, as will be discussed below) relative to the road and the output f corresponds to the distributed tangential force per length (hereafter, tangential traction) acting on the contact patch. The total longitudinal and lateral forces on the contact patch are obtained by the spatial integration of the traction f.Because the LuGre friction model is inherently a smooth model, LuGre-based tire models either do not realize the exact static friction and their presliding displacement cannot be purely elastic. The practical disadvantage of this particular feature has not been pointed out in the tire-modeling literature. It however could cause some unrealistic simulation results, such as those in which a parked vehicle slowly drifts away when it is subjected to small external forces, even if it is below the maximum static friction force. Such behavior will be illustrated through numerical examples in Section4.2.Another concern regarding the original LuGre friction model (3) is that the function Ψ, which determines the velocity dependency of the friction force, resides in the denominator. Due to this structure, the multidimensional extension including the anisotropy is not straightforward. With Deur et al.’s [1] two-dimensional LuGre model, the friction force is determined in a component-wise manner, in which the friction force in each direction is determined only by the velocity in the same direction. This approach is not consistent with the principle of the maximum dissipation [29]. In contrast, Velenis et al. [7] take the principle of the maximum dissipation into account in two extreme cases: the cases of the zero velocity and the sufficiently high velocity. It uses the same Stribeck curve to represent the velocity dependency of the friction force in both x and y directions, which may not accord with reality.The central assumption of the brush-type tire model is that each bristle is attached to the rigid wheel and interacts with the road individually, without interacting with neighboring bristles. In pneumatic tires, however, the rigid wheel is covered with a compliant carcass, and its deformation is much larger than the microscopic deflection of the surface bristles. It has been pointed out that the carcass compliance strongly influences the cornering stiffness [30] and lateral forces [8]. Precise modeling of the carcass deformation is not easy because of its continuum nature. Some researchers proposed models involving the carcass compliance [8, 30–32], in which the brush-type models, some of which are LuGre-based tire models, are modified so that the bristles protrude from the deformed carcass centerline, instead of the rigid wheel. Previous studies have been successful with some strong assumptions for the deformed carcass shape, such as piece-wise straight lines [8] and quadratic curve [31, 32], which depend on the total lateral force. It should be emphasized that brush-type models have been useful as the basis for such extensions.


### 2.3. Effects of Distribution of the Vertical Pressure

In the application of conventional friction models to tire modeling, one needs to take into account the fact that the vertical pressure varies according to the location within the contact patch. Because the LuGre friction model does not explicitly consider the variable vertical force, some modifications need to be made.One reasonable assumption, which has been adopted by most of the previous work, is that, at each point, the tangential traction is proportional to the vertical pressure via the kinetic friction coefficient. Besides this, one can also assume that, at each point and in the steady state, the tangential traction is proportional to the tangential deflection of the surface. These assumptions can be described as follows:(8)f=μfz≈Ke,where fz denotes the vertical pressure, f denotes the tangential traction, μ denotes the kinetic friction coefficient, e denotes the tangential deflection of the contact patch, and K denotes the tangential stiffness of the contact patch. Note that, in the SI units, f and K are measured in N/m and N/m2, respectively.In the original one-dimensional LuGre friction model, where the effect of the normal pressure is not considered, the output friction force is obtained asf≈σ0e in the steady state, where e is a state variable and σ0>0 is a stiffness factor. With the effect of the normal pressure taken into account, the majority of the previous LuGre-based tire models [1, 6–8] are constructed so that the following is satisfied:(9)μ≈σ0e.It results in the following:(10)f=μfz≈fzσ0e.This relation can be seen in many previous works, e.g., [7] (equation (28)), [1] (equation (18)), and [8] (equation (11)). Comparing (10) to (8), one can see that(11)K=fzσ0.This means that, if σ0 is assumed to be constant as has been assumed in the previous works [1, 7, 8], the tangential stiffness coefficient K is proportional to the vertical pressure fz.Another possible approach is to keepK in (8) constant, instead of σ0. Comparison of these two schemes has been discussed in the context of one-dimensional friction models [33, 34] and tire models [1]. Marques et al. [34] compared such two versions of LuGre friction models and concluded that, with abrupt changes in the normal force, the constant-K approach results in unrealistic friction force. Many of the previous LuGre-based tire models [1, 6–8, 35] use the constant-σ0 approach, which may be justifiable [1] by the fact that the tire material actually becomes stiffer as the vertical pressure increases. Some researchers [36, 37] pointed out that this approach results in unrealistic distribution of the deflection e near the trailing edge, presumably because the stiffness becomes close to zero there. Liang et al. [37] have provided a method to circumvent this flaw, in which the input velocity is multiplied by a normalized vertical pressure. Its physical meaning, however, is not fully clear.


## 2.1. One-Dimensional Friction Models

Many friction models have ever been proposed for control and simulation purposes [18, 19]. Among them, the LuGre friction model [5] is a model that has been extended into some tire models. The original LuGre friction model can be described as follows:(3)e˙=v−KveΨv,f=Ke+Be˙.Here, v∈ℝ is the input velocity, f∈ℝ is the output friction force, e∈ℝ is a state variable that is interpreted as the average deflection of contact bristles, and K and B are positive constants representing the stiffness and damping coefficients, respectively, of the contact bristles. The function Ψ:ℝ⇉ℝ is a function representing a velocity-dependent friction law, which is typically given as follows:(4)Ψv=ΔsgnvFs−Fs−Fks1vVs,γ,where(5)s1u,γ=Δ1−exp−uγ,and Fs, Fk, Vs, and γ are positive constants. The constants Fs and Fk can be interpreted as the maximum static friction force and the kinetic friction force, respectively, which usually satisfy Fs≥Fk. With this function Ψ, the magnitude of the friction force decreases from Fs to Fk as the velocity increases along a curve characterized by Vs and γ. The function Ψ defined above is set-valued but the right-hand side of the first equation of (3) is always single-valued and continuous because it is zero when v=0. That is, the LuGre friction model (3) can be seen as an ODE, which is smooth and mathematically convenient.The form (3) of the model is inherited from the Dahl model [20], which has been chosen so that(6)f≈Ψv≈Keis achieved when the input velocity v maintains a constant nonzero value. A physical interpretation of the model (3) is often represented as in Figure 1. In this model, there is an elastic deflection e in the frictional interface and the friction force f is obtained as the viscoelastic force determined by the second equation of (3).Figure 1
A physical model upon which the LuGre friction model [5] and Kikuuwe et al.’s [10] friction model are built.Because the LuGre friction model neglects the nonsmooth nature of friction, it does not make the distinction of two different states: the static friction state and the kinetic friction state. In this model, the state variablee can be understood as the elastic deflection of the bristle, but its rate-of-change is not equal to the input velocity v. In other words, the velocity of the bristle tip cannot be the exact zero, and thus the static friction state cannot be realized. As a consequence, the object’s displacement is not purely elastic even in the presliding regime and unbounded drift can happen even under external forces smaller than the maximum static friction force. This point has been formally pointed out by Theorem 1 of Dupont et al. [9]. Some friction models to circumvent this problem have been proposed so far [21–25], with special attention on reproducing hysteresis behavior in the elastic displacement in the static friction state.In contrast to the smooth LuGre model, Kikuuwe et al.’s [10] friction model is a nonsmooth model. Their friction model is described as the following DAI:(7)f∈Ψv−e˙,f=Ke+Be˙.This can be viewed as a pair of simultaneous differential equations that determines f and e˙ so that both equations are satisfied (Mathematical expressions like (7), which involves the symbol “∈” and derivatives, are referred to as differential inclusions (DIs) [26] instead of differential equations. Moreover, because the derivative e˙ in (7) cannot be isolated in the left-hand side, expression (7) can also be referred to as a differential algebraic inclusion (DAI)). This model can also be interpreted as in Figure 1 and can be said to be more faithful to Figure 1 than the LuGre friction model because it takes into account the exact balance between the friction force, f in the first equation of (7), and the viscoelastic force, f in the second equation of (7). Kikuuwe et al. [10] presented an approach to implement this model with the backward Euler discretization, and Xiong et al. [25, 27] extended it into the continuous-time domain. By setting B=0, Kikuuwe et al.’s [10] friction model reduces to Hayward and Armstrong’s [28] friction model, which is defined only in the discrete-time domain.


## 2.2. LuGre-Based Brush-Type Tire Models

As overviewed in Section1, many of the brush-type tire models are based on the LuGre friction model (3). In their approach, the input velocity v is the velocity of the rigid wheel (or of the deformed carcass centerline, as will be discussed below) relative to the road and the output f corresponds to the distributed tangential force per length (hereafter, tangential traction) acting on the contact patch. The total longitudinal and lateral forces on the contact patch are obtained by the spatial integration of the traction f.Because the LuGre friction model is inherently a smooth model, LuGre-based tire models either do not realize the exact static friction and their presliding displacement cannot be purely elastic. The practical disadvantage of this particular feature has not been pointed out in the tire-modeling literature. It however could cause some unrealistic simulation results, such as those in which a parked vehicle slowly drifts away when it is subjected to small external forces, even if it is below the maximum static friction force. Such behavior will be illustrated through numerical examples in Section4.2.Another concern regarding the original LuGre friction model (3) is that the function Ψ, which determines the velocity dependency of the friction force, resides in the denominator. Due to this structure, the multidimensional extension including the anisotropy is not straightforward. With Deur et al.’s [1] two-dimensional LuGre model, the friction force is determined in a component-wise manner, in which the friction force in each direction is determined only by the velocity in the same direction. This approach is not consistent with the principle of the maximum dissipation [29]. In contrast, Velenis et al. [7] take the principle of the maximum dissipation into account in two extreme cases: the cases of the zero velocity and the sufficiently high velocity. It uses the same Stribeck curve to represent the velocity dependency of the friction force in both x and y directions, which may not accord with reality.The central assumption of the brush-type tire model is that each bristle is attached to the rigid wheel and interacts with the road individually, without interacting with neighboring bristles. In pneumatic tires, however, the rigid wheel is covered with a compliant carcass, and its deformation is much larger than the microscopic deflection of the surface bristles. It has been pointed out that the carcass compliance strongly influences the cornering stiffness [30] and lateral forces [8]. Precise modeling of the carcass deformation is not easy because of its continuum nature. Some researchers proposed models involving the carcass compliance [8, 30–32], in which the brush-type models, some of which are LuGre-based tire models, are modified so that the bristles protrude from the deformed carcass centerline, instead of the rigid wheel. Previous studies have been successful with some strong assumptions for the deformed carcass shape, such as piece-wise straight lines [8] and quadratic curve [31, 32], which depend on the total lateral force. It should be emphasized that brush-type models have been useful as the basis for such extensions.


## 2.3. Effects of Distribution of the Vertical Pressure

In the application of conventional friction models to tire modeling, one needs to take into account the fact that the vertical pressure varies according to the location within the contact patch. Because the LuGre friction model does not explicitly consider the variable vertical force, some modifications need to be made.One reasonable assumption, which has been adopted by most of the previous work, is that, at each point, the tangential traction is proportional to the vertical pressure via the kinetic friction coefficient. Besides this, one can also assume that, at each point and in the steady state, the tangential traction is proportional to the tangential deflection of the surface. These assumptions can be described as follows:(8)f=μfz≈Ke,where fz denotes the vertical pressure, f denotes the tangential traction, μ denotes the kinetic friction coefficient, e denotes the tangential deflection of the contact patch, and K denotes the tangential stiffness of the contact patch. Note that, in the SI units, f and K are measured in N/m and N/m2, respectively.In the original one-dimensional LuGre friction model, where the effect of the normal pressure is not considered, the output friction force is obtained asf≈σ0e in the steady state, where e is a state variable and σ0>0 is a stiffness factor. With the effect of the normal pressure taken into account, the majority of the previous LuGre-based tire models [1, 6–8] are constructed so that the following is satisfied:(9)μ≈σ0e.It results in the following:(10)f=μfz≈fzσ0e.This relation can be seen in many previous works, e.g., [7] (equation (28)), [1] (equation (18)), and [8] (equation (11)). Comparing (10) to (8), one can see that(11)K=fzσ0.This means that, if σ0 is assumed to be constant as has been assumed in the previous works [1, 7, 8], the tangential stiffness coefficient K is proportional to the vertical pressure fz.Another possible approach is to keepK in (8) constant, instead of σ0. Comparison of these two schemes has been discussed in the context of one-dimensional friction models [33, 34] and tire models [1]. Marques et al. [34] compared such two versions of LuGre friction models and concluded that, with abrupt changes in the normal force, the constant-K approach results in unrealistic friction force. Many of the previous LuGre-based tire models [1, 6–8, 35] use the constant-σ0 approach, which may be justifiable [1] by the fact that the tire material actually becomes stiffer as the vertical pressure increases. Some researchers [36, 37] pointed out that this approach results in unrealistic distribution of the deflection e near the trailing edge, presumably because the stiffness becomes close to zero there. Liang et al. [37] have provided a method to circumvent this flaw, in which the input velocity is multiplied by a normalized vertical pressure. Its physical meaning, however, is not fully clear.


## 3. New Nonsmooth Representation

### 3.1. PDAI Representation

Here, we present a new mathematical representation of the brush-type tire model, taking the nonsmooth nature of friction into account. In this model, the wheel is assumed to be a rigid circular ring. The surface of the tire is viscoelastically deformable in both the longitudinal and lateral directions. The existence of the carcass, which intervenes between the tire surface and the wheel in a pneumatic tire, is totally neglected. The tire is rolling in the longitudinal direction and is also slipping both in the longitudinal and lateral directions. The friction force on the tire surface produces the deformation of the surface. The contact patch is modeled as a one-dimensional line segment as shown in Figure2. The distribution of the pressure along the contact patch is not uniform. The surface viscoelasticity and the friction are allowed to be anisotropic.Figure 2
Contact patch and the coordinate system. (a) Tire and the coordinate system. (b) Contact patch, top view. The yaw rateωz is measured positive counterclockwise from the top view.
(a)
(b)LetR∈ℝ+ be the radius of the tire and L∈ℝ+ be the length of the contact patch. As shown in Figure 2, the origin of the coordinate system Σt is set at the leading edge of the contact patch, and the x direction, measured by ξ∈0,L, is set along the contact patch. The z direction is in the vertical direction, positive upward, and the y direction is set parallel to the axle, i.e., perpendicular to both the x and z directions.Let the relative velocity of the ground with respect to the tire be denoted byVt∈ℝ2, which is measured in Σt. Let ωzt∈ℝ be the yaw rate of the tire and ωyt∈ℝ be the rolling angular velocity of the tire around its axle. Then, the velocity of the particles (i.e., bristles) at the location ξ on the contact patch can be described as the following function vξ,t of the location ξ and the time t:(12)vξ,t=Vt−Rωytξ−Ξcωzt,where Ξc∈0,L is the location of the center of the yaw rotation of the tire, which is typically set as Ξc=L/2.The vertical force applied to the contact patch is denoted byFz∈ℝ+. The pressure (the vertical force per length) on the contact patch can be written as Fzf^zξ/L/L where f^z:0,1⟶ℝ+ is a function that satisfies(13)∫01f^zξ^dξ^=1.For simplicity, we here assume that the distribution does not vary along time, although it in reality may depend on the vehicle velocity as suggested in [38].Now let us discuss the traction vectorfξ,t∈ℝ2, which is a function of the location ξ∈0,L and the time t. Let eξ,t∈ℝ2 be the horizontal deflection of the bristle at the location of ξ∈0,L and time t∈ℝ+. Then, the traction vector fξ,t is determined by the Coulomb-like friction law and also by the viscoelasticity of the surface bristle, which are described as the following PDAI:(14)fξ,t∈Fzf^zξ/LLΨvξ,t−DDteξ,t,fξ,t=Keξ,t+BDDteξ,t,where D/Dt denotes the material derivative, which can be defined as follows:(15)DDteξ,t=Rωyt∂∂ξeξ,t+∂∂teξ,t.Here, Ψ:ℝ2⇉ℝ2 is a friction coefficient vector function, which is set-valued at the origin. The simplest example of the function Ψ is Ψx=μSgnx, which represents the pure isotropic Coulomb friction. In (14), K and B are 2×2 positive definite matrices representing stiffness and viscosity coefficients, respectively, which are distributed along the contact patch. The values of the functions fξ,t and eξ,t are determined by PDAI (14) based on the inputs vξ,t and ωyt. The boundary condition is e0,t=0 when ωyt>0 and eL,t=0 when ωyt<0.Based on the traction vectorfξ,t determined by (14), the total force and the total moment acting on the tire can be obtained as follows:(16)FxtFyt=∫0Lfξ,tdξ,Mzt=∫0Lξ−Ξcfyξ,tdξ,where fyξ,t is the y component of fξ,t, Fx and Fy are the total longitudinal and lateral forces, respectively, and Mz is the total moment acting on the tire, which is often referred to as the self-aligning torque.Remark 1.
One important feature of the model (14) is that, in the steady states, the bristle deflection eξ,t is proportional to the tangential traction fξ,t via a fixed stiffness matrix K. That is, the model (14) adopts the constant-K approach mentioned in Section 2.3. It is in contrast to the constant-σ0 approach (e.g., [1, 6–8]), which erroneously results in zero stiffness at the trailing edge. On the other hand, Marques et al. [34] pointed out that the constant-K approach applied to the LuGre friction model is not feasible to deal with fast variation in the vertical load. One way to avoid such problems in the model (14) may be setting K and B as functions of the vertical pressure fzξ, which would be consistent with the fact that the tire material stiffens when it is compressed. The model (14), based on the constant-K approach, reserves the possibility of such an extension.


### 3.2. Choice of FunctionΨ: Two-Dimensional Stribeck Effect

Here we discuss the choice of the nonlinear functionΨ. Let μs∗ and μk∗ denote the static and kinetic friction coefficients, respectively, in the ∗∗∈x,y direction and Vsx, Vsy, and γ denote positive constants characterizing Stribeck curves. One requirement for the function Ψ:ℝ2⇉ℝ2 is that it should satisfy the following:(17)Ψv,0T=sgnvμsx−s1vVsx,γμsx−μkx0,Ψ0,vT=0sgnvμsy−s1vVsy,γμsy−μky,where s1 is the one defined in (5). That is, the function Ψ needs to produce Stribeck curves both in x and y directions. In addition, it is preferable that the function Ψ should possess a dissipation potential [29], which is a scalar function U:ℝ2⟶ℝ+ that satisfies the following:(18)∂vUv=Ψv.The existence of such a function Uv is not physically a requirement but is preferred for the convenience of computation, as shown in Section 3.3.Now, let us define the following functions and matrices:(19)s2u,γ=Δ1−ΓL1/γ,uγuγ,s3u,γ=Δ∂s2u,γ∂u=s1u−s2uu,Vs=ΔdiagVsx,Vsy∈ℝ2×2,Mk=Δdiagμkx,μky∈ℝ2×2,Ms=Δdiagμsx,μsy∈ℝ2×2,where ΓL:ℝ×ℝ⟶ℝ is the lower incomplete Gamma function defined as follows: (20)ΓLa,z=Δ∫0zξa−1exp−ξdξ.With these functions and matrices, this paper proposes the following dissipation potential to realize a two-dimensional Stribeck effect:(21)Uv=ΔMsv−s2Vs−1v,γMsv−Mkv.A straightforward derivation shows that this scalar function Uv satisfies the following:(22)Ψv=Δ∂vUv=SgnMsv−s2Vs−1v,γSgnMsv−SgnMkv−s3Vs−1v,γMsv−MkvSgnVs−1v.It can be seen that this function Ψv satisfies (17). With Mk=Ms=μI, the function Ψv in (22) reduces to the pure isotropic Coulomb friction law, which is Uv=μv and Ψv=μSgnv. It should be noted that, with γ=1, the function ΓL reduces to the following:(23)ΓL1γ,uγ=∫0uexp−ξdξ=1−exp−u.Therefore, the implementation is computationally inexpensive in the special case of γ=1 than the general case.Remark 2.
As discussed in Section2.2, there have been some different extensions of the LuGre friction model into the two-dimensional case. As far as the author is aware, there has been no multidimensional friction force law that possesses a dissipation potential and produces the Stribeck effect.


### 3.3. Simulation Algorithm

Let us discuss the discretization of PDAI (14) for temporal and spatial integration. Let T∈ℝ+ and Z∈ℝ+ be the step sizes for the time t and the space ξ, respectively. The two terms of the right-hand side of (15) can be discretized as follows:(24)ωyt∂∂ξeξ,t≈ωyteξ,t−eξ−sgnωytZ,tZ,∂∂teξ,t≈eξ,t−eξ,t−TT.Thus, (15) can be approximated as follows:(25)DDteξ,t≈eξ,t−e∗ξ,tT∗t,where(26)e∗ξ,t=ΔZeξ,t−T+TRωyteξ−sgnωytZ,tZ+TRωyt,T∗t=ΔTZZ+TRωyt.Here, the set-valuedness of the part sgnωyt does not take effect because it is multiplied by ωyt in the second term of the numerator.By using (25), PDAI (14) can be discretized as follows:(27)fξ,t∈Fzf^zξ/LLΨvξ,t−eξ,t−e∗ξ,tT∗t,fξ,t=B^teξ,t−Be∗ξ,tT∗t,where(28)B^t=ΔKT∗t+B.Eliminating eξ,t from (27) yields the following:(29)fξ,t∈Fzf^zξ/LLΨvξ,t+B^t−1Ke∗ξ,t−fξ,t,which includes fξ,t in both sides. In order to obtain fξ,t, one needs to solve the algebraic inclusion (generalized equation) (29) with respect to fξ,t.Here, let us define another functionΨ^:ℝ2×ℝ2×2⟶ℝ2 as follows:(30)Ψ^v,X=Δys.t.y=Ψv−X−1y,ifXv∉Ψ0,Xv,ifXv∈Ψ0.This definition of the function Ψ^ is not of the closed form, but even if it is not obtained analytically, it can be obtained numerically or approximately. With this function, the solution fξ,t of (25) can be obtained as follows:(31)fξ,t=Fzf^zξ/LLΨ^v∗ξ,t,LB^tFzf^zξ/L,where(32)v∗ξ,t=Δvξ,t+B^t−1Ke∗ξ,t.Based on this, eξ,t can also be computed as follows:(33)eξ,t=B^t−1T∗tfξ,t+Be∗ξ,t.If definition (22) of Ψ is adopted, definition (30) of Ψ^ is specialized into the following:(34)Ψ^v,X=Δys.t.y=Ψv−X−1y,ifMs−1Xv>1,Xv,ifMs−1Xv≤1.In the first line, we need to solve the nonlinear equation y=Ψv−X−1y with respect to y. By setting x=v−X−1y and considering the dissipation potential Uv, this algebraic equation can be rewritten as(35)∂xUx−Xv−x=0.It can be seen that the solution of (35) coincides with the solution of the following minimization problem:(36)minimizex∈ℝ2Ux−vTXx+xTXx2.This problem can be solved with some numerical schemes such as the steepest gradient method.In conclusion, the spatio-temporal integration of PDAI (14) can be performed by Algorithm 1. The output ekk∈0,N+1 at a timestep is used as one of the inputs in the next timestep.Algorithm 1:Proposed algorithm.
RequireT,N,R,Ξc,L,K,B,f^z,Ψ^
(1)
function ALG1 V,ωy,ωz,Fz,ekk∈0,N+1
(2)
Fx≔0;Fy≔0;Mz≔0;
(3)
Z≔L/N
(4)
T∗≔TZ/Z+TRωy
(5)
B^≔KT∗+B
(6)
if ωy≥0 then kd≔1else kd≔−1 endif
(7)
k≔N+11−kd/2
(8)
ek≔0
(9)
loop
(10)
k+=kd
(11)
if k<1 or N<kthen break endif
(12)
ξ≔k−1/2Z
(13)
e∗≔Zek+TRωyek−kd/Z+TRωy
(14)
v∗≔V−Rωyξ−Ξcωz+B^−1Ke∗
(15)
f≔Fzf^zξ/L/LΨ^v∗,LB^/Fzf^zξ/L
(16)
ek≔B^−1T∗f+Be∗
(17)
FxFy+=fZ
(18)
Mz+=ξ−ΞcfyZ
(19)
end loop
(20)
return Fx,Fy,Mz,ekk∈0,N+1
(21)
end function


## 3.1. PDAI Representation

Here, we present a new mathematical representation of the brush-type tire model, taking the nonsmooth nature of friction into account. In this model, the wheel is assumed to be a rigid circular ring. The surface of the tire is viscoelastically deformable in both the longitudinal and lateral directions. The existence of the carcass, which intervenes between the tire surface and the wheel in a pneumatic tire, is totally neglected. The tire is rolling in the longitudinal direction and is also slipping both in the longitudinal and lateral directions. The friction force on the tire surface produces the deformation of the surface. The contact patch is modeled as a one-dimensional line segment as shown in Figure2. The distribution of the pressure along the contact patch is not uniform. The surface viscoelasticity and the friction are allowed to be anisotropic.Figure 2
Contact patch and the coordinate system. (a) Tire and the coordinate system. (b) Contact patch, top view. The yaw rateωz is measured positive counterclockwise from the top view.
(a)
(b)LetR∈ℝ+ be the radius of the tire and L∈ℝ+ be the length of the contact patch. As shown in Figure 2, the origin of the coordinate system Σt is set at the leading edge of the contact patch, and the x direction, measured by ξ∈0,L, is set along the contact patch. The z direction is in the vertical direction, positive upward, and the y direction is set parallel to the axle, i.e., perpendicular to both the x and z directions.Let the relative velocity of the ground with respect to the tire be denoted byVt∈ℝ2, which is measured in Σt. Let ωzt∈ℝ be the yaw rate of the tire and ωyt∈ℝ be the rolling angular velocity of the tire around its axle. Then, the velocity of the particles (i.e., bristles) at the location ξ on the contact patch can be described as the following function vξ,t of the location ξ and the time t:(12)vξ,t=Vt−Rωytξ−Ξcωzt,where Ξc∈0,L is the location of the center of the yaw rotation of the tire, which is typically set as Ξc=L/2.The vertical force applied to the contact patch is denoted byFz∈ℝ+. The pressure (the vertical force per length) on the contact patch can be written as Fzf^zξ/L/L where f^z:0,1⟶ℝ+ is a function that satisfies(13)∫01f^zξ^dξ^=1.For simplicity, we here assume that the distribution does not vary along time, although it in reality may depend on the vehicle velocity as suggested in [38].Now let us discuss the traction vectorfξ,t∈ℝ2, which is a function of the location ξ∈0,L and the time t. Let eξ,t∈ℝ2 be the horizontal deflection of the bristle at the location of ξ∈0,L and time t∈ℝ+. Then, the traction vector fξ,t is determined by the Coulomb-like friction law and also by the viscoelasticity of the surface bristle, which are described as the following PDAI:(14)fξ,t∈Fzf^zξ/LLΨvξ,t−DDteξ,t,fξ,t=Keξ,t+BDDteξ,t,where D/Dt denotes the material derivative, which can be defined as follows:(15)DDteξ,t=Rωyt∂∂ξeξ,t+∂∂teξ,t.Here, Ψ:ℝ2⇉ℝ2 is a friction coefficient vector function, which is set-valued at the origin. The simplest example of the function Ψ is Ψx=μSgnx, which represents the pure isotropic Coulomb friction. In (14), K and B are 2×2 positive definite matrices representing stiffness and viscosity coefficients, respectively, which are distributed along the contact patch. The values of the functions fξ,t and eξ,t are determined by PDAI (14) based on the inputs vξ,t and ωyt. The boundary condition is e0,t=0 when ωyt>0 and eL,t=0 when ωyt<0.Based on the traction vectorfξ,t determined by (14), the total force and the total moment acting on the tire can be obtained as follows:(16)FxtFyt=∫0Lfξ,tdξ,Mzt=∫0Lξ−Ξcfyξ,tdξ,where fyξ,t is the y component of fξ,t, Fx and Fy are the total longitudinal and lateral forces, respectively, and Mz is the total moment acting on the tire, which is often referred to as the self-aligning torque.Remark 1.
One important feature of the model (14) is that, in the steady states, the bristle deflection eξ,t is proportional to the tangential traction fξ,t via a fixed stiffness matrix K. That is, the model (14) adopts the constant-K approach mentioned in Section 2.3. It is in contrast to the constant-σ0 approach (e.g., [1, 6–8]), which erroneously results in zero stiffness at the trailing edge. On the other hand, Marques et al. [34] pointed out that the constant-K approach applied to the LuGre friction model is not feasible to deal with fast variation in the vertical load. One way to avoid such problems in the model (14) may be setting K and B as functions of the vertical pressure fzξ, which would be consistent with the fact that the tire material stiffens when it is compressed. The model (14), based on the constant-K approach, reserves the possibility of such an extension.


## 3.2. Choice of FunctionΨ: Two-Dimensional Stribeck Effect

Here we discuss the choice of the nonlinear functionΨ. Let μs∗ and μk∗ denote the static and kinetic friction coefficients, respectively, in the ∗∗∈x,y direction and Vsx, Vsy, and γ denote positive constants characterizing Stribeck curves. One requirement for the function Ψ:ℝ2⇉ℝ2 is that it should satisfy the following:(17)Ψv,0T=sgnvμsx−s1vVsx,γμsx−μkx0,Ψ0,vT=0sgnvμsy−s1vVsy,γμsy−μky,where s1 is the one defined in (5). That is, the function Ψ needs to produce Stribeck curves both in x and y directions. In addition, it is preferable that the function Ψ should possess a dissipation potential [29], which is a scalar function U:ℝ2⟶ℝ+ that satisfies the following:(18)∂vUv=Ψv.The existence of such a function Uv is not physically a requirement but is preferred for the convenience of computation, as shown in Section 3.3.Now, let us define the following functions and matrices:(19)s2u,γ=Δ1−ΓL1/γ,uγuγ,s3u,γ=Δ∂s2u,γ∂u=s1u−s2uu,Vs=ΔdiagVsx,Vsy∈ℝ2×2,Mk=Δdiagμkx,μky∈ℝ2×2,Ms=Δdiagμsx,μsy∈ℝ2×2,where ΓL:ℝ×ℝ⟶ℝ is the lower incomplete Gamma function defined as follows: (20)ΓLa,z=Δ∫0zξa−1exp−ξdξ.With these functions and matrices, this paper proposes the following dissipation potential to realize a two-dimensional Stribeck effect:(21)Uv=ΔMsv−s2Vs−1v,γMsv−Mkv.A straightforward derivation shows that this scalar function Uv satisfies the following:(22)Ψv=Δ∂vUv=SgnMsv−s2Vs−1v,γSgnMsv−SgnMkv−s3Vs−1v,γMsv−MkvSgnVs−1v.It can be seen that this function Ψv satisfies (17). With Mk=Ms=μI, the function Ψv in (22) reduces to the pure isotropic Coulomb friction law, which is Uv=μv and Ψv=μSgnv. It should be noted that, with γ=1, the function ΓL reduces to the following:(23)ΓL1γ,uγ=∫0uexp−ξdξ=1−exp−u.Therefore, the implementation is computationally inexpensive in the special case of γ=1 than the general case.Remark 2.
As discussed in Section2.2, there have been some different extensions of the LuGre friction model into the two-dimensional case. As far as the author is aware, there has been no multidimensional friction force law that possesses a dissipation potential and produces the Stribeck effect.


## 3.3. Simulation Algorithm

Let us discuss the discretization of PDAI (14) for temporal and spatial integration. Let T∈ℝ+ and Z∈ℝ+ be the step sizes for the time t and the space ξ, respectively. The two terms of the right-hand side of (15) can be discretized as follows:(24)ωyt∂∂ξeξ,t≈ωyteξ,t−eξ−sgnωytZ,tZ,∂∂teξ,t≈eξ,t−eξ,t−TT.Thus, (15) can be approximated as follows:(25)DDteξ,t≈eξ,t−e∗ξ,tT∗t,where(26)e∗ξ,t=ΔZeξ,t−T+TRωyteξ−sgnωytZ,tZ+TRωyt,T∗t=ΔTZZ+TRωyt.Here, the set-valuedness of the part sgnωyt does not take effect because it is multiplied by ωyt in the second term of the numerator.By using (25), PDAI (14) can be discretized as follows:(27)fξ,t∈Fzf^zξ/LLΨvξ,t−eξ,t−e∗ξ,tT∗t,fξ,t=B^teξ,t−Be∗ξ,tT∗t,where(28)B^t=ΔKT∗t+B.Eliminating eξ,t from (27) yields the following:(29)fξ,t∈Fzf^zξ/LLΨvξ,t+B^t−1Ke∗ξ,t−fξ,t,which includes fξ,t in both sides. In order to obtain fξ,t, one needs to solve the algebraic inclusion (generalized equation) (29) with respect to fξ,t.Here, let us define another functionΨ^:ℝ2×ℝ2×2⟶ℝ2 as follows:(30)Ψ^v,X=Δys.t.y=Ψv−X−1y,ifXv∉Ψ0,Xv,ifXv∈Ψ0.This definition of the function Ψ^ is not of the closed form, but even if it is not obtained analytically, it can be obtained numerically or approximately. With this function, the solution fξ,t of (25) can be obtained as follows:(31)fξ,t=Fzf^zξ/LLΨ^v∗ξ,t,LB^tFzf^zξ/L,where(32)v∗ξ,t=Δvξ,t+B^t−1Ke∗ξ,t.Based on this, eξ,t can also be computed as follows:(33)eξ,t=B^t−1T∗tfξ,t+Be∗ξ,t.If definition (22) of Ψ is adopted, definition (30) of Ψ^ is specialized into the following:(34)Ψ^v,X=Δys.t.y=Ψv−X−1y,ifMs−1Xv>1,Xv,ifMs−1Xv≤1.In the first line, we need to solve the nonlinear equation y=Ψv−X−1y with respect to y. By setting x=v−X−1y and considering the dissipation potential Uv, this algebraic equation can be rewritten as(35)∂xUx−Xv−x=0.It can be seen that the solution of (35) coincides with the solution of the following minimization problem:(36)minimizex∈ℝ2Ux−vTXx+xTXx2.This problem can be solved with some numerical schemes such as the steepest gradient method.In conclusion, the spatio-temporal integration of PDAI (14) can be performed by Algorithm 1. The output ekk∈0,N+1 at a timestep is used as one of the inputs in the next timestep.Algorithm 1:Proposed algorithm.
RequireT,N,R,Ξc,L,K,B,f^z,Ψ^
(1)
function ALG1 V,ωy,ωz,Fz,ekk∈0,N+1
(2)
Fx≔0;Fy≔0;Mz≔0;
(3)
Z≔L/N
(4)
T∗≔TZ/Z+TRωy
(5)
B^≔KT∗+B
(6)
if ωy≥0 then kd≔1else kd≔−1 endif
(7)
k≔N+11−kd/2
(8)
ek≔0
(9)
loop
(10)
k+=kd
(11)
if k<1 or N<kthen break endif
(12)
ξ≔k−1/2Z
(13)
e∗≔Zek+TRωyek−kd/Z+TRωy
(14)
v∗≔V−Rωyξ−Ξcωz+B^−1Ke∗
(15)
f≔Fzf^zξ/L/LΨ^v∗,LB^/Fzf^zξ/L
(16)
ek≔B^−1T∗f+Be∗
(17)
FxFy+=fZ
(18)
Mz+=ξ−ΞcfyZ
(19)
end loop
(20)
return Fx,Fy,Mz,ekk∈0,N+1
(21)
end function


## 4. Properties of the Model: Numerical Results

This section presents some simulation results illustrating the behaviors of the proposed model composed of (14) and (16) and Algorithm 1.


### 4.1.Fx, Fy, and Mz Curves

The proposed model composed of (14) and (16), implemented as Algorithm 1, is now fit to the earliest version of Magic Formula [15], which is a conventional tire model that is empirically known to be accurate for pneumatic tires. The parameter values of the Magic Formula were obtained from the original Bakker et al.’s report [15]. The results of Magic formula are often represented by three curves: the κ-Fx, α-Fy, and α-Mz curves, where κ is the slip ratio in the percentage, α is the slip angle in degrees, and Fx, Fy, and Mz are the total force and moments acting on the tire. The parameter values for the proposed model were chosen as in Table 1, which were obtained in a procedure similar to the one presented in the author and his colleague’s previous conference publication [39] and some handcrafting. This paper does not attempt to seek further optimality.Table 1
Parameters used for the proposed model with four different vertical loadsFz. The parameters qz and nz are for the function f^z, appearing in (14), according to the definition described as equation (50) of [39], which is originally proposed in [38]. The discretization was performed with T=0.001s and N=10.
Symbol
Unit
F
z
=
2000 N
F
z
=
4000 N
F
z
=
6000 N
F
z
=
8000 N
L
m
0.083
0.176
0.250
0.345
R
m
0.294
Ξ
c
m
L
/
2
γ
—
1.0
V
s
m/s
9.0
I
q
z
—
0.321
0.273
0.184
0.147
n
z
—
30.
30.
30.
29.3
K
x
N
/
m
2
1.75
×
10
7
7.90
×
10
6
6.85
×
10
6
4.92
×
10
6
K
y
N
/
m
2
1.17
×
10
7
3.84
×
10
6
2.01
×
10
6
9.77
×
10
5
B
Ns
/
m
2
0.0003
K
μ
s
x
—
1.18
1.17
1.11
1.06
μ
s
y
—
1.08
0.92
0.91
0.88
μ
k
x
—
0.75
0.60
0.60
0.49
μ
k
y
—
0.63
0.89
0.79
0.75Simulations with the new algorithm were performed in two scenarios, Scenarios A and B. Scenario A was a straight, pure-braking motion with a slowly varying slip ratioκ, and Scenario B was a pure-cornering motion with slowly varying slip angle α. The duration of each scenario was TL=Δ100 s. In Scenario A, the inputs to the algorithm were as follows:(37)Vt=V00T,ωyt=V0100−κt100R,ωzt=0,κt=100tTL,where V0=16.67m/s (≈60km/h). In this scenario, the data of κt and Fxt were recorded. In Scenario B, the inputs to the algorithm were as follows:(38)Vt=V0cosπαt180V0sinπαt180T,ωyt=V0Rcosπαt180,ωzt=πα˙t180,αt=20sin2πtTL,where V0=19.44m/s (≈70km/h). In this scenario, the data of Fyt, Mzt, and αt were recorded. The values of V0 in these scenarios were taken from the original Bakker et al.’s report [15]. The vertical load Fz was chosen as 2000, 4000, 6000, and 8000 N because Bakker et al.’s report [15] shows the Magic Formula for these values of the vertical load Fz.The results are shown in Figure3. The κ-Fx and α-Fy curves of the proposed model well match the Magic Formula curves, reproducing the smooth transition between the sticking (low κ or low α) to the slipping (high κ or high α). In contrast, the α-Mz curves of the proposed model do not match well with Magic Formula. Such mismatch can be attributed to the fact that the brush-type model inherently does not take the elastic carcass into account. It however still possesses quantitatively similar features to the Magic Formula α-MZ curves, such as the smooth transition from the sticking to the slipping and the negative aligning torque in the large slip angle.Figure 3
Steady-state results of the proposed Algorithm1 (red) and Magic Formula (gray thick) with the parameter values in Table 1 of [15] (the κ-Fx graphs were obtained by Scenario A and the other graphs were obtained by Scenario B).It should be emphasized that the purpose of this comparative study is to clarify how closely the new algorithm is capable of approximating pneumatic tires. This paper does not attempt to argue the superiority of the model over existing models in terms of the accuracy with respect to the real automobile behaviors because the brush-type tire model itself does not capture important mechanisms, such as carcass deformation, of the real pneumatic tires.


### 4.2. Bicycle Model: Behavior under External Forces

To investigate the applicability of the proposed model to vehicle simulation, a planar bicycle model illustrated in Figure4, which is a typical simplified model of a symmetric four-wheeled vehicle, was constructed. The parameters for both front and rear tires were chosen as those of Fz=4000N in Table 1.Figure 4
Planar bicycle model used in the simulations of Scenarios C to H. The dimensions of the vehicle were set asLf=Lr=1.215m and the mass was set as 800 kg.The most important feature of the proposed model is its set-valuedness for representing the static friction, in contrast to the continuity of the LuGre-based tire models [1, 7]. The previous LuGre-based tire models have quite different structures from the proposed model, which make the direct comparison difficult. Thus, we here constructed another LuGre-based tire model that is comparable with the proposed model. Let us recall that the original DAI-based friction model is (7) and the original one-dimensional LuGre model is (3). Considering that the DAI-based model (7) is extended into the tire model (14), one can extend the LuGre friction model (3) into the following tire model:(39)fΨξ,t=Fzf^zξ/LLΨvξ,t,Deξ,tDt=vξ,t−Kxvxξ,texξ,tfΨ,xξ,tKyvyξ,teyξ,tfΨ,yξ,tT,fξ,t=Keξ,t+BDDteξ,t.The set-valuedness of the function Ψ is not a concern here because the second term of the left-hand side of the second equation of (39) is zero if vξ,t=0. This model, specifically the second equation of (39), is the axis-wise application of the one-dimensional LuGre model (3) involving the material derivative, in the same way as in equation (17) of [1] and equation (28) of [7]. Implicit Euler discretization of the above equations leads to Algorithm 1 with the lines 14 to 16 being replaced by the following ones:(40)v≔V−Rωyξ−Ξcωc,fΨ≔Fzf^zξ/LLΨv,H≔diagKxvxξ,tfΨ,x,Kyvyξ,tfΨ,y,ei≔IT∗+H−1v+e∗T∗,f≔Kei+Bei−e∗T∗.We performed simulations in which the bicycle model was parked being subjected to an oscillatory external force acting on the chassis. The wheel axles were locked, and the steering angle was fixed asθs=0. Simulations of two scenarios, Scenarios C and D, were performed. An oscillatory external force was applied in the longitudinal direction in Scenario C and in the lateral direction in Scenario D. The results are shown in Figure 5. In both scenarios, the external force acted from t=1s to 8 s, changing its magnitude between 936 N and 3120 N sinusoidally with the frequency of 2 Hz. It is shown that the drift is much smaller with the proposed model than the LuGre-based model. In Figure 5(a), one can see that the proposed model does produce a small drift under the longitudinal external force; by seeing that the displacement after the external force vanishes is around 0.001 m, not the exact zero. This small drift is likely attributed to the error in discretization (25) of the material derivative, which should be addressed in the future study. In Figure 5(b), under the lateral external force, one can see that the vehicle does not drift with the proposed model. This feature can be counted as an advantage of the proposed algorithm over LuGre-based tire models.Figure 5
Results of Scenarios C and D: response of a parked vehicle against oscillatory external forces, with the fixed steering angleθs=0 and the initial conditions x=y=0 and θ=0. (a) Scenario C, longitudinal displacement in the presence of longitudinal external force and (b) Scenario D, lateral displacement in the presence of lateral external force.
(a)
(b)In another set of simulations, Scenario E, a driving torque of 1005 N·m, was constantly applied to the front wheel axle of the vehicle and an oscillatory external force was applied to the chassis in the positivey direction, which is leftward from the chassis, mimicking the situation of driving in a side wind. The magnitude of the external force was changed in the same way as in Scenarios C and D. The results are shown in Figure 6. One can see that the vehicle deviates from the line y=0 and that the proposed model produces a smaller deviation. This kind of behaviors would happen in reality because, as the tires rotate, the lateral elastic deformations in the tires are cumulatively converted into the lateral drift of the vehicle. Considering the fact the LuGre-based tire model drifts even in the parked situations, the behavior of the proposed model is inferred to be closer to reality.Figure 6
Scenario E: Trajectory of the vehicle with the fixed steering angleθs=0 and the initial conditions x=y=0 and θ=0 subjected to a lateral (positive y-directional) oscillatory external force.


### 4.3. Bicycle Model: Locked Cornering

Some simulations were performed to test the cornering behaviors of the proposed tire model. Three scenarios, Scenarios F, G, and H, were performed with the same vehicle model as the one in Section4.2. In all scenarios, the front axle was driven with the driving torque of 1005 N·m and the steering angle θs was changed from zero to θs=−5 deg at t=5s. In Scenario F, the vehicle continued running after t=5s. In Scenario G, both axles were locked at t=7s. In Scenario H, only the rear axle was locked at t=7s.The results are shown in Figure7. Figure 7(a) shows that, in Scenario F, the vehicle properly performs a cornering motion. This cornering trajectory is also included in Figures 7(b) and 7(c) in gray as reference. Figure 7(b) shows that, in Scenario G, the tire locking in both axles results in understeer, deviating from the normal cornering trajectory indicated in gray, as is the case with reality. Figure 7(c) shows that the rear-tire locking, while the front tire is not being locked, causes a sharp turn of the vehicle pivoting around the front tire, which is known as a “handbrake turn.” It should be noted the load transfer between tires was not considered here; i.e., the load to each tire was constant even during the braking. These results suggest the potential usefulness of the proposed model for the applications to computer graphics and simple simulation.Figure 7
Results of Scenarios F (a), G (b), and H (c). Black curves indicate the trajectories of the vehicle. The gray curves indicate the trajectory of Scenario F without braking, for reference. The rectangles indicate the locations and the orientations of the vehicle at the intervals of 0.5 (s) in blue when the either or both axles were locked and in red otherwise.
(a)
(b)
(c)


## 4.1.Fx, Fy, and Mz Curves

The proposed model composed of (14) and (16), implemented as Algorithm 1, is now fit to the earliest version of Magic Formula [15], which is a conventional tire model that is empirically known to be accurate for pneumatic tires. The parameter values of the Magic Formula were obtained from the original Bakker et al.’s report [15]. The results of Magic formula are often represented by three curves: the κ-Fx, α-Fy, and α-Mz curves, where κ is the slip ratio in the percentage, α is the slip angle in degrees, and Fx, Fy, and Mz are the total force and moments acting on the tire. The parameter values for the proposed model were chosen as in Table 1, which were obtained in a procedure similar to the one presented in the author and his colleague’s previous conference publication [39] and some handcrafting. This paper does not attempt to seek further optimality.Table 1
Parameters used for the proposed model with four different vertical loadsFz. The parameters qz and nz are for the function f^z, appearing in (14), according to the definition described as equation (50) of [39], which is originally proposed in [38]. The discretization was performed with T=0.001s and N=10.
Symbol
Unit
F
z
=
2000 N
F
z
=
4000 N
F
z
=
6000 N
F
z
=
8000 N
L
m
0.083
0.176
0.250
0.345
R
m
0.294
Ξ
c
m
L
/
2
γ
—
1.0
V
s
m/s
9.0
I
q
z
—
0.321
0.273
0.184
0.147
n
z
—
30.
30.
30.
29.3
K
x
N
/
m
2
1.75
×
10
7
7.90
×
10
6
6.85
×
10
6
4.92
×
10
6
K
y
N
/
m
2
1.17
×
10
7
3.84
×
10
6
2.01
×
10
6
9.77
×
10
5
B
Ns
/
m
2
0.0003
K
μ
s
x
—
1.18
1.17
1.11
1.06
μ
s
y
—
1.08
0.92
0.91
0.88
μ
k
x
—
0.75
0.60
0.60
0.49
μ
k
y
—
0.63
0.89
0.79
0.75Simulations with the new algorithm were performed in two scenarios, Scenarios A and B. Scenario A was a straight, pure-braking motion with a slowly varying slip ratioκ, and Scenario B was a pure-cornering motion with slowly varying slip angle α. The duration of each scenario was TL=Δ100 s. In Scenario A, the inputs to the algorithm were as follows:(37)Vt=V00T,ωyt=V0100−κt100R,ωzt=0,κt=100tTL,where V0=16.67m/s (≈60km/h). In this scenario, the data of κt and Fxt were recorded. In Scenario B, the inputs to the algorithm were as follows:(38)Vt=V0cosπαt180V0sinπαt180T,ωyt=V0Rcosπαt180,ωzt=πα˙t180,αt=20sin2πtTL,where V0=19.44m/s (≈70km/h). In this scenario, the data of Fyt, Mzt, and αt were recorded. The values of V0 in these scenarios were taken from the original Bakker et al.’s report [15]. The vertical load Fz was chosen as 2000, 4000, 6000, and 8000 N because Bakker et al.’s report [15] shows the Magic Formula for these values of the vertical load Fz.The results are shown in Figure3. The κ-Fx and α-Fy curves of the proposed model well match the Magic Formula curves, reproducing the smooth transition between the sticking (low κ or low α) to the slipping (high κ or high α). In contrast, the α-Mz curves of the proposed model do not match well with Magic Formula. Such mismatch can be attributed to the fact that the brush-type model inherently does not take the elastic carcass into account. It however still possesses quantitatively similar features to the Magic Formula α-MZ curves, such as the smooth transition from the sticking to the slipping and the negative aligning torque in the large slip angle.Figure 3
Steady-state results of the proposed Algorithm1 (red) and Magic Formula (gray thick) with the parameter values in Table 1 of [15] (the κ-Fx graphs were obtained by Scenario A and the other graphs were obtained by Scenario B).It should be emphasized that the purpose of this comparative study is to clarify how closely the new algorithm is capable of approximating pneumatic tires. This paper does not attempt to argue the superiority of the model over existing models in terms of the accuracy with respect to the real automobile behaviors because the brush-type tire model itself does not capture important mechanisms, such as carcass deformation, of the real pneumatic tires.


## 4.2. Bicycle Model: Behavior under External Forces

To investigate the applicability of the proposed model to vehicle simulation, a planar bicycle model illustrated in Figure4, which is a typical simplified model of a symmetric four-wheeled vehicle, was constructed. The parameters for both front and rear tires were chosen as those of Fz=4000N in Table 1.Figure 4
Planar bicycle model used in the simulations of Scenarios C to H. The dimensions of the vehicle were set asLf=Lr=1.215m and the mass was set as 800 kg.The most important feature of the proposed model is its set-valuedness for representing the static friction, in contrast to the continuity of the LuGre-based tire models [1, 7]. The previous LuGre-based tire models have quite different structures from the proposed model, which make the direct comparison difficult. Thus, we here constructed another LuGre-based tire model that is comparable with the proposed model. Let us recall that the original DAI-based friction model is (7) and the original one-dimensional LuGre model is (3). Considering that the DAI-based model (7) is extended into the tire model (14), one can extend the LuGre friction model (3) into the following tire model:(39)fΨξ,t=Fzf^zξ/LLΨvξ,t,Deξ,tDt=vξ,t−Kxvxξ,texξ,tfΨ,xξ,tKyvyξ,teyξ,tfΨ,yξ,tT,fξ,t=Keξ,t+BDDteξ,t.The set-valuedness of the function Ψ is not a concern here because the second term of the left-hand side of the second equation of (39) is zero if vξ,t=0. This model, specifically the second equation of (39), is the axis-wise application of the one-dimensional LuGre model (3) involving the material derivative, in the same way as in equation (17) of [1] and equation (28) of [7]. Implicit Euler discretization of the above equations leads to Algorithm 1 with the lines 14 to 16 being replaced by the following ones:(40)v≔V−Rωyξ−Ξcωc,fΨ≔Fzf^zξ/LLΨv,H≔diagKxvxξ,tfΨ,x,Kyvyξ,tfΨ,y,ei≔IT∗+H−1v+e∗T∗,f≔Kei+Bei−e∗T∗.We performed simulations in which the bicycle model was parked being subjected to an oscillatory external force acting on the chassis. The wheel axles were locked, and the steering angle was fixed asθs=0. Simulations of two scenarios, Scenarios C and D, were performed. An oscillatory external force was applied in the longitudinal direction in Scenario C and in the lateral direction in Scenario D. The results are shown in Figure 5. In both scenarios, the external force acted from t=1s to 8 s, changing its magnitude between 936 N and 3120 N sinusoidally with the frequency of 2 Hz. It is shown that the drift is much smaller with the proposed model than the LuGre-based model. In Figure 5(a), one can see that the proposed model does produce a small drift under the longitudinal external force; by seeing that the displacement after the external force vanishes is around 0.001 m, not the exact zero. This small drift is likely attributed to the error in discretization (25) of the material derivative, which should be addressed in the future study. In Figure 5(b), under the lateral external force, one can see that the vehicle does not drift with the proposed model. This feature can be counted as an advantage of the proposed algorithm over LuGre-based tire models.Figure 5
Results of Scenarios C and D: response of a parked vehicle against oscillatory external forces, with the fixed steering angleθs=0 and the initial conditions x=y=0 and θ=0. (a) Scenario C, longitudinal displacement in the presence of longitudinal external force and (b) Scenario D, lateral displacement in the presence of lateral external force.
(a)
(b)In another set of simulations, Scenario E, a driving torque of 1005 N·m, was constantly applied to the front wheel axle of the vehicle and an oscillatory external force was applied to the chassis in the positivey direction, which is leftward from the chassis, mimicking the situation of driving in a side wind. The magnitude of the external force was changed in the same way as in Scenarios C and D. The results are shown in Figure 6. One can see that the vehicle deviates from the line y=0 and that the proposed model produces a smaller deviation. This kind of behaviors would happen in reality because, as the tires rotate, the lateral elastic deformations in the tires are cumulatively converted into the lateral drift of the vehicle. Considering the fact the LuGre-based tire model drifts even in the parked situations, the behavior of the proposed model is inferred to be closer to reality.Figure 6
Scenario E: Trajectory of the vehicle with the fixed steering angleθs=0 and the initial conditions x=y=0 and θ=0 subjected to a lateral (positive y-directional) oscillatory external force.


## 4.3. Bicycle Model: Locked Cornering

Some simulations were performed to test the cornering behaviors of the proposed tire model. Three scenarios, Scenarios F, G, and H, were performed with the same vehicle model as the one in Section4.2. In all scenarios, the front axle was driven with the driving torque of 1005 N·m and the steering angle θs was changed from zero to θs=−5 deg at t=5s. In Scenario F, the vehicle continued running after t=5s. In Scenario G, both axles were locked at t=7s. In Scenario H, only the rear axle was locked at t=7s.The results are shown in Figure7. Figure 7(a) shows that, in Scenario F, the vehicle properly performs a cornering motion. This cornering trajectory is also included in Figures 7(b) and 7(c) in gray as reference. Figure 7(b) shows that, in Scenario G, the tire locking in both axles results in understeer, deviating from the normal cornering trajectory indicated in gray, as is the case with reality. Figure 7(c) shows that the rear-tire locking, while the front tire is not being locked, causes a sharp turn of the vehicle pivoting around the front tire, which is known as a “handbrake turn.” It should be noted the load transfer between tires was not considered here; i.e., the load to each tire was constant even during the braking. These results suggest the potential usefulness of the proposed model for the applications to computer graphics and simple simulation.Figure 7
Results of Scenarios F (a), G (b), and H (c). Black curves indicate the trajectories of the vehicle. The gray curves indicate the trajectory of Scenario F without braking, for reference. The rectangles indicate the locations and the orientations of the vehicle at the intervals of 0.5 (s) in blue when the either or both axles were locked and in red otherwise.
(a)
(b)
(c)


## 5. Concluding Remarks

This paper has proposed a new tire friction model based on nonsmooth, set-valued Coulomb-like friction law. In the new model, the contact patch is modeled as a line segment on which the viscoelastic displacement and Coulomb-like friction force are distributed. The model is described as a PDAI that represents the balance between the linear viscoelastic force and the Coulomb-like friction force. A numerical integration algorithm has been derived based on the implicit Euler discretization. Some numerical results have been presented.The main focus of the paper has been placed on constructing a generic representation of the nonsmooth nature of the conventional brush-type model structure. As an intrinsic limitation of the brush-type structure, the presented model may not be accurate for pneumatic tires, which involve the carcass deformation. The presented model would be useful for applications such as gaming and computer graphics, which do not require very high accuracy, and also for simulations of vehicles with nonpneumatic tires such as mobile robots. To be able to be applied to precise engineering simulation of automobiles, the presented model would need extensions to include other factors, such as carcass dynamics, nonzero camber angle, and ply steer. Such extensions will be achieved through an elaborate but relatively straightforward application of existing approaches as they do not involve mathematical difficulties such as nonsmoothness. Guidelines for parameter tuning should also be sought in the future study.



---

*Source: 9747605-2019-12-17.xml*