车辆运动控制(3)轮胎模型
1. 简介

 在 《车辆运动控制(2)车辆横摆动力学建模》 中根据单轨模型分析得到2自由度的车辆横摆动力学微分方程:
 
     
      
       
        
         
         
          
           
            {
           
           
            
             
              
               
                
                 
                  
                   v
                  
                  
                   ˙
                  
                 
                 
                  y
                 
                
                
                 =
                
                
                 −
                
                
                 
                  v
                 
                 
                  x
                 
                
                
                 
                  φ
                 
                 
                  ˙
                 
                
                
                 +
                
                
                 
                  1
                 
                 
                  m
                 
                
                
                 [
                
                
                 
                  F
                 
                 
                  
                   y
                  
                  
                   f
                  
                 
                
                
                 cos
                
                
                 
                
                
                 (
                
                
                 
                  δ
                 
                 
                  f
                 
                
                
                 )
                
                
                 +
                
                
                 
                  F
                 
                 
                  
                   y
                  
                  
                   r
                  
                 
                
                
                 ]
                
               
              
             
            
            
             
              
               
              
             
            
            
             
              
               
                
                 
                  φ
                 
                 
                  ¨
                 
                
                
                 =
                
                
                 
                  1
                 
                 
                  
                   I
                  
                  
                   z
                  
                 
                
                
                 (
                
                
                 
                  l
                 
                 
                  f
                 
                
                
                 
                  F
                 
                 
                  
                   y
                  
                  
                   f
                  
                 
                
                
                 cos
                
                
                 
                
                
                 (
                
                
                 
                  δ
                 
                 
                  f
                 
                
                
                 )
                
                
                 −
                
                
                 
                  l
                 
                 
                  r
                 
                
                
                 
                  F
                 
                 
                  
                   y
                  
                  
                   r
                  
                 
                
                
                 )
                
               
              
             
            
           
          
         
         
         
          
           (16)
          
         
        
       
       
         \begin{cases} \dot{v}_y=-v_x\dot{\varphi} +\frac{1}{m}[F_{yf}\cos(\delta_f)+F_{yr}]\\\\ \ddot{\varphi}=\frac{1}{I_z}(l_fF_{yf}\cos(\delta_f)-l_rF_{yr}) \end{cases} \tag{16} 
       
      
     ⎩⎪⎨⎪⎧v˙y=−vxφ˙+m1[Fyfcos(δf)+Fyr]φ¨=Iz1(lfFyfcos(δf)−lrFyr)(16)
其中, F y f , F y r F_{yf},F_{yr} Fyf,Fyr 车辆前、后轴上轮胎侧向力的合力还需要进一步通计算而得
在车辆的运动过程中
 轮胎所受的纵向力、侧向力、垂直力及回正力矩对汽车的操纵稳定性和安全性起着重要作用
 由于轮胎结构复杂,动力学性能呈非线性
 选择符合实际又便于使用的轮胎模型是建立车辆动力学模型的关键
2. 轮胎模型
目前,主要的轮胎模型可以分为:
- 理论轮胎模型
 - 经验轮胎模型
 - 物理轮胎模型
 
虽然轮胎模型各有不同
 但是这些模型都力求更加准确地 表达轮胎和地面的函数关系
 并且都是 描述轮胎的输人和输出之间的关系
3. Pacejka 轮胎模型
其中,最常见的是 Pacejka 提出的以 魔术公式(Magic Formula, MF) 为基础的半经验轮胎模型
 此模型运用三角函数的组合公式拟合轮胎试验数据
 描述轮胎的纵向力
    
     
      
       
        
         F
        
        
         x
        
       
      
      
       F_x
      
     
    Fx、侧向力
    
     
      
       
        
         F
        
        
         y
        
       
      
      
       F_y
      
     
    Fy、回正力矩
    
     
      
       
        
         M
        
        
         z
        
       
      
      
       M_z
      
     
    Mz,、翻转力矩
    
     
      
       
        
         M
        
        
         x
        
       
      
      
       M_x
      
     
    Mx,、阻力矩
    
     
      
       
        
         M
        
        
         y
        
       
      
      
       M_y
      
     
    My,与侧偏角
    
     
      
       
        α
       
      
      
       \alpha
      
     
    α、滑移率之间的定量关系以及纵向力、侧向力的联合作用工况
 能够表达不同驱动情况时的轮胎特性
魔术公式的一般表达式为:
 
     
      
       
        
         
         
          
           
            Y
           
           
            (
           
           
            x
           
           
            )
           
           
            =
           
           
            D
           
           
            sin
           
           
            
           
           
            {
           
           
            C
           
           
            arctan
           
           
            
           
           
            [
           
           
            B
           
           
            x
           
           
            −
           
           
            E
           
           
            (
           
           
            B
           
           
            x
           
           
            −
           
           
            arctan
           
           
            
           
           
            (
           
           
            B
           
           
            x
           
           
            )
           
           
            )
           
           
            ]
           
           
            }
           
          
         
         
         
          
           (17)
          
         
        
       
       
        Y(x) = D\sin\{ C\arctan[Bx - E(Bx - \arctan ( Bx))]\} \tag{17}
       
      
     Y(x)=Dsin{Carctan[Bx−E(Bx−arctan(Bx))]}(17)
| 符号 | 含义 | 符号 | 含义 | 
|---|---|---|---|
| Y Y Y | 输出变量 轮胎的纵向力 F x F_x Fx 或 侧向力 F y F_y Fy 或 回正力矩 M z M_z Mz  | x x x | 输入变量 轮胎的侧偏角 α \alpha α 或 纵向滑移率  | 
| D D D | 峰值因子 | C C C | 形状因子 | 
| B B B | 刚度因子 | E E E | 曲率因子 | 
以上四个因子由轮胎的 垂向载荷 F z F_z Fz和 外倾角 γ \gamma γ 确定
Pacejka 轮胎模型 的输人量与输出量之间的关系如图所示:

4. 轮胎侧向力
4.1. 根据魔术公式计算
在横摆动力学建模中,主要关注轮胎的侧向力 F y F_y Fy
利用魔术公式计算轮胎侧向力的公式如下:
 
     
      
       
        
         
         
          
           
            
             F
            
            
             y
            
           
           
            =
           
           
            D
           
           
            sin
           
           
            
           
           
            {
           
           
            C
           
           
            arctan
           
           
            
           
           
            [
           
           
            B
           
           
            x
           
           
            −
           
           
            E
           
           
            (
           
           
            B
           
           
            x
           
           
            −
           
           
            arctan
           
           
            
           
           
            (
           
           
            B
           
           
            x
           
           
            )
           
           
            )
           
           
            ]
           
           
            }
           
           
            +
           
           
            
             S
            
            
             v
            
           
          
         
         
         
          
           (18)
          
         
        
       
       
        F_y = D\sin\{ C\arctan[Bx - E(Bx - \arctan ( Bx))]\} + S_v \tag{18}
       
      
     Fy=Dsin{Carctan[Bx−E(Bx−arctan(Bx))]}+Sv(18)
 其中,
    
     
      
       
        x
       
       
        =
       
       
        a
       
       
        +
       
       
        
         S
        
        
         h
        
       
      
      
       x=a+S_h
      
     
    x=a+Sh
| 符号 | 含义 | 符号 | 含义 | 符号 | 含义 | 
|---|---|---|---|---|---|
| α \alpha α | 轮胎侧偏角 | S h S_h Sh, S v S_v Sv | 曲线的水平、垂直方向漂移 | E E E | 曲率因子 | 
| D D D | 峰值因子 | C C C | 形状因子 | B B B | 刚度因子 | 
| F z F_z Fz | 轮胎垂向载荷 | γ \gamma γ | 轮胎外倾角 | 
魔术公式的参数
    
     
      
       
        A
       
      
      
       A
      
     
    A可以根据轮胎实验数据拟合得到
 当轮胎与道路之间的附着系数为1.0时,取值如下表:

水平方向漂移:
    
     
      
       
        
         S
        
        
         h
        
       
       
        =
       
       
        
         A
        
        
         9
        
       
       
        
         F
        
        
         z
        
       
       
        +
       
       
        
         A
        
        
         10
        
       
       
        +
       
       
        
         A
        
        
         8
        
       
       
        γ
       
      
      
       S_h=A_9F_z+A_{10}+A_{8}\gamma
      
     
    Sh=A9Fz+A10+A8γ
 垂直方向漂移:
    
     
      
       
        
         S
        
        
         v
        
       
       
        =
       
       
        
         A
        
        
         11
        
       
       
        
         F
        
        
         z
        
       
       
        γ
       
       
        +
       
       
        
         A
        
        
         12
        
       
       
        
         F
        
        
         z
        
       
       
        +
       
       
        
         A
        
        
         13
        
       
      
      
       S_v=A_{11}F_z\gamma+A_{12}F_z+A_{13}
      
     
    Sv=A11Fzγ+A12Fz+A13
 刚度因子:
    
     
      
       
        B
       
       
        =
       
       
        
         A
        
        
         3
        
       
       
        sin
       
       
        
       
       
        [
       
       
        2
       
       
        arctan
       
       
        
       
       
        (
       
       
        
         F
        
        
         z
        
       
       
        /
       
       
        
         A
        
        
         4
        
       
       
        )
       
       
        ]
       
       
        ∗
       
       
        (
       
       
        1
       
       
        −
       
       
        
         A
        
        
         5
        
       
       
        ∣
       
       
        γ
       
       
        ∣
       
       
        )
       
       
        /
       
       
        (
       
       
        C
       
       
        ∗
       
       
        D
       
       
        )
       
      
      
       B = A_3\sin[2\arctan(F_z/A_4)] * (1 - A_5 | \gamma |)/(C* D)
      
     
    B=A3sin[2arctan(Fz/A4)]∗(1−A5∣γ∣)/(C∗D)
 形状因子:
    
     
      
       
        C
       
       
        =
       
       
        
         A
        
        
         0
        
       
      
      
       C=A_0
      
     
    C=A0
 峰值因子:
    
     
      
       
        D
       
       
        =
       
       
        
         A
        
        
         1
        
       
       
        
         F
        
        
         z
        
        
         2
        
       
       
        +
       
       
        
         A
        
        
         2
        
       
       
        
         F
        
        
         z
        
       
      
      
       D=A_1F_z^2+A_2F_z
      
     
    D=A1Fz2+A2Fz
 曲率因子:
    
     
      
       
        E
       
       
        =
       
       
        
         A
        
        
         6
        
       
       
        
         F
        
        
         z
        
       
       
        +
       
       
        
         A
        
        
         7
        
       
      
      
       E=A_6F_z +A_7
      
     
    E=A6Fz+A7
在 《车辆运动控制(2)车辆横摆动力学建模》 的 理想化假设3 只考虑 纯侧偏角 特性下
 根据公式18计算不同轮胎垂向载荷 
    
     
      
       
        
         F
        
        
         z
        
       
      
      
       F_z
      
     
    Fz 所对应的轮胎侧向力
    
     
      
       
        
         F
        
        
         y
        
       
      
      
       F_y
      
     
    Fy:
import matplotlib.pyplot as plt
import math
# 魔术公式参数
A = [1.65, -34, 1250, 3036, 12.8, 0.00501, -0.02103,  0.77394,
     0.0022890, 0.013442, 0.003709, 19.1656, 1.21356, 6.26206]
# 轮胎垂向载荷 kN
Fz_list = [2.5, 5, 8.5, 14]
# 轮胎侧偏角 °
alpha_list = range(-8, 9)
def getFy(Fz, alpha, gamma=0):
    '''
    计算轮胎侧向力
    Fz:轮胎垂向载荷
    alpha:轮胎侧偏角
    gamma:轮胎外倾角,假设为0
    '''
    # 水平方向漂移
    Sh = A[9] * Fz + A[10] + A[8] * gamma
    # 垂直方向漂移
    Sv = A[11] * Fz * gamma + A[12] * Fz + A[13]
    # 形状因子
    C = A[0]
    # 峰值因子
    D = A[1] * Fz**2 + A[2] * Fz
    # 侧向力零点处的侧向刚度
    Ca = A[3] * math.sin(2 * math.atan(Fz / A[4])) * (1 - A[5] * abs(gamma))
    # 刚度因子
    B = Ca / (C * D)
    # 曲率因子
    E = A[6] * Fz**2 + A[7]
    # 输入
    x = alpha + Sh
    # 轮胎侧向力
    Fy = D * math.sin(C * math.atan(B * x - E * (B * x - math.atan(B * x)))) + Sv
    return Fy / 1000
FyOutput = [[0 for c in range(len(alpha_list))]for r in range(len(Fz_list))]
for x in range(len(Fz_list)):
    for y in range(len(alpha_list)):
        FyOutput[x][y] = getFy(Fz_list[x], alpha_list[y])
# 设置正常显示中文
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 设置曲线数值
for x in range(len(Fz_list)):
    plt.plot(alpha_list, FyOutput[x], label='Fz={}kN'.format(Fz_list[x]))
plt.xticks(alpha_list)  # 设置X轴坐标数值标识
plt.xlabel('轮胎侧偏角/(°)')  # 设置X轴的名字
plt.ylabel('轮胎侧向力/(kN)')  # 设置Y轴的名字
plt.title("根据公式18计算的不同垂向载荷下轮胎侧向力图")   # 设置标题
plt.legend()  # 设置图例
plt.show()  # 显示图表
 

4.2. 根据刷子轮胎模型计算
Pacejka 提出的 刷子轮胎模型 也经常被用来计算轮胎侧向力
    
     
      
       
        
         F
        
        
         y
        
       
      
      
       F_y
      
     
    Fy,其计算公式为:
 
     
      
       
        
         
         
          
           
            
             F
            
            
             y
            
           
           
            =
           
           
            
             {
            
            
             
              
               
                
                 
                  
                   C
                  
                  
                   a
                  
                 
                 
                  tan
                 
                 
                  
                 
                 
                  (
                 
                 
                  α
                 
                 
                  )
                 
                 
                  −
                 
                 
                  
                   
                    C
                   
                   
                    a
                   
                   
                    2
                   
                  
                  
                   
                    3
                   
                   
                    μ
                   
                   
                    
                     F
                    
                    
                     z
                    
                   
                  
                 
                 
                  ∣
                 
                 
                  tan
                 
                 
                  
                 
                 
                  (
                 
                 
                  α
                 
                 
                  )
                 
                 
                  ∣
                 
                 
                  tan
                 
                 
                  
                 
                 
                  (
                 
                 
                  α
                 
                 
                  )
                 
                 
                  +
                 
                 
                  
                   
                    C
                   
                   
                    a
                   
                   
                    3
                   
                  
                  
                   
                    27
                   
                   
                    
                     μ
                    
                    
                     2
                    
                   
                   
                    
                     F
                    
                    
                     z
                    
                    
                     2
                    
                   
                  
                 
                 
                  
                   
                    tan
                   
                   
                    
                   
                  
                  
                   3
                  
                 
                 
                  (
                 
                 
                  α
                 
                 
                  )
                 
                 
                  ,
                 
                 
                  ∣
                 
                 
                  α
                 
                 
                  ∣
                 
                 
                  <
                 
                 
                  arctan
                 
                 
                  
                 
                 
                  (
                 
                 
                  
                   
                    3
                   
                   
                    μ
                   
                   
                    
                     F
                    
                    
                     z
                    
                   
                  
                  
                   
                    C
                   
                   
                    a
                   
                  
                 
                 
                  )
                 
                
               
              
             
             
              
               
                
               
              
             
             
              
               
                
                 
                  μ
                 
                 
                  
                   F
                  
                  
                   z
                  
                 
                 
                  s
                 
                 
                  g
                 
                 
                  n
                 
                 
                  (
                 
                 
                  α
                 
                 
                  )
                 
                 
                  ,
                 
                 
                  ∣
                 
                 
                  α
                 
                 
                  ∣
                 
                 
                  ≥
                 
                 
                  arctan
                 
                 
                  
                 
                 
                  (
                 
                 
                  
                   
                    3
                   
                   
                    μ
                   
                   
                    
                     F
                    
                    
                     z
                    
                   
                  
                  
                   
                    C
                   
                   
                    a
                   
                  
                 
                 
                  )
                 
                
               
              
             
            
           
          
         
         
         
          
           (19)
          
         
        
       
       
         F_y= \begin{cases} C_a\tan(\alpha)-\frac{C_a^2}{3\mu F_z} | \tan(\alpha) | \tan(\alpha)+\frac{C_a^3}{27\mu^2 F_z^2}\tan^3(\alpha),| \alpha | < \arctan(\frac{3\mu F_z}{C_a})\\\\ \mu F_z sgn(\alpha),| \alpha | ≥ \arctan(\frac{3\mu F_z}{C_a}) \end{cases} \tag{19} 
       
      
     Fy=⎩⎪⎨⎪⎧Catan(α)−3μFzCa2∣tan(α)∣tan(α)+27μ2Fz2Ca3tan3(α),∣α∣<arctan(Ca3μFz)μFzsgn(α),∣α∣≥arctan(Ca3μFz)(19)
| 符号 | 含义 | 符号 | 含义 | 符号 | 含义 | 
|---|---|---|---|---|---|
| α \alpha α | 轮胎侧偏角 | C a C_a Ca | 轮胎的侧偏刚度 | μ \mu μ | 轮胎与道路之间的附着系数 | 
根据公式19,不同 垂直载荷 F z F_z Fz 下 轮胎侧向力 F y F_y Fy 与 轮胎侧偏角 α \alpha α 的关系如图:
import matplotlib.pyplot as plt
import math
# 魔术公式参数
A = [1.65, -34, 1250, 3036, 12.8, 0.00501, -0.02103,  0.77394,
     0.0022890, 0.013442, 0.003709, 19.1656, 1.21356, 6.26206]
# 轮胎垂向载荷 kN
Fz_list = [2800, 8000, 14000, 19000]
# 轮胎侧偏角 °
alpha_list = range(-12, 13)
def sgn(a):
    if a > 0:
        return 1
    elif a < 0:
        return -1
    else:
        return 0
def getFy(Fz, alpha, mu=1, gamma=0):
    '''
    计算轮胎侧向力
    Fz:轮胎垂向载荷
    alpha:轮胎侧偏角
    mu:轮胎与道路之间的附着系数
    gamma:轮胎外倾角,假设为0
    '''
    # 角度转换
    alpha = alpha / 180 * math.pi
    # 轮胎的侧偏刚度
    Ca =  A[3] * math.sin(2 * math.atan(Fz / 1000 / A[4])) * (1 - A[5] * abs(gamma)) * 180 / math.pi
    if abs(alpha) < math.atan((3 * mu * Fz) / Ca):
        Fy = Ca * math.atan(alpha) - Ca**2 /(3 * mu * Fz) * abs(math.atan(alpha)) * math.atan(alpha) + Ca**3 / (27 * mu**2 * Fz**2) * math.atan(alpha)**3
    else:
        Fy = mu * Fz * sgn(alpha)
    return Fy
FyOutput = [[0 for c in range(len(alpha_list))]for r in range(len(Fz_list))]
for x in range(len(Fz_list)):
    for y in range(len(alpha_list)):
        FyOutput[x][y] = getFy(Fz_list[x], alpha_list[y])
# 设置正常显示中文
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 设置曲线数值
for x in range(len(Fz_list)):
    plt.plot(alpha_list, FyOutput[x], label='Fz={}N'.format(Fz_list[x]))
plt.xticks(alpha_list)  # 设置X轴坐标数值标识
plt.xlabel('轮胎侧偏角/(°)')  # 设置X轴的名字
plt.ylabel('轮胎侧向力/(N)')  # 设置Y轴的名字
plt.title("根据公式19计算的不同垂向载荷下轮胎侧向力图")   # 设置标题
plt.legend()  # 设置图例
plt.show()  # 显示图表
 

5. 线性化
根据上两个图所示的不同载荷下轮胎侧向力与轮胎侧偏角的关系曲线可以看出
 当轮胎侧偏角较小时,轮胎侧向力
    
     
      
       
        
         F
        
        
         y
        
       
      
      
       F_y
      
     
    Fy 可以近似表示为 轮胎侧偏角
    
     
      
       
        α
       
      
      
       \alpha
      
     
    α 的线性函数:
 
     
      
       
        
         
         
          
           
            
             F
            
            
             y
            
           
           
            =
           
           
            
             
              C
             
             
              ˉ
             
            
            
             a
            
           
           
            α
           
          
         
         
         
          
           (20)
          
         
        
       
       
        F_y = \bar{C}_a \alpha\tag{20}
       
      
     Fy=Cˉaα(20)
其中, C ˉ a \bar{C}_a Cˉa 为轮胎的线性侧偏刚度
此轮胎线性化模型在 侧向加速度 
    
     
      
       
        
         α
        
        
         y
        
       
       
        ≤
       
       
        0.4
       
       
        g
       
      
      
       \alpha_y ≤ 0.4g
      
     
    αy≤0.4g 、轮胎侧偏角 
    
     
      
       
        α
       
       
        <
       
       
        5
       
       
        °
       
      
      
       \alpha<5°
      
     
    α<5° 的情景下
 对常规轮胎具有较高的拟合精度
 因此应用公式20时需对轮胎侧偏角的范围进行约束
谢谢










