本篇結合Adams中的TR_rear_pac89.tir文件，介紹一下魔術公式的基本知識。使用魔術公式的輪胎模型主要有Pacejka ’89、Pacejka ’94、MF-Tyre、MF-Swift四種。
1. Pacejka ’89和’94輪胎模型
Pacejka ’89 和’94輪胎模型是以魔術公式主要提出者H. B. Pacejka教授命名的，根據其發布的年限命名。目前有兩種直接被ADAMS引用。
魔術公式是用三角函數的組合公式擬合輪胎試驗數據，用一套形式相同的公式就可以完整地表達輪胎的縱向力Fx、側向力Fy、回正力矩Mz、翻轉力矩Mx、阻力矩My以及縱向力、側向力的聯合作用工況，故稱為“魔術公式”。
魔術公式的一般表達式為：
式中Y(x)可以是側向力，也可以是回正力矩或者縱向力，自變量x可以在不同的情況下分別表示輪胎的側偏角或縱向滑移率，式中的系數B、C、D依次由輪胎的垂直載荷和外傾角來確定。
Pacejka ’89輪胎模型認為輪胎在垂直、側向方向上是線性的、阻尼為常量，這在側向加速度常見范圍≤0.4g，側偏角≤5°的情景下對常規輪胎具有很高的擬合精度。此外，由於魔術公式基於試驗數據，除在試驗范圍的高精度外，甚至在極限值以外一定程度仍可使用，可以對有限工況進行外推且具有較好的置信度。魔術公式正在成為工業標准，即輪胎制造商向整車廠提供魔術公式系數表示的輪胎數據，而不再是表格或圖形。
基於魔術公式的輪胎模型還有較好的健壯性，如果沒有某一輪胎的試驗數據，而使用同類輪胎數據替代仍可取得很好的效果。
圖 基於魔術公式的輪胎模型的輸入和輸出變量
Pacejka ’89輪胎力與力矩的計算
輪胎縱向力計算公式為：
其中X1為縱向力組合自變量：X1=(κ+Sh)，κ為縱向滑移率（負值出現在制動態，-100表示車輪抱死）
C——曲線形狀因子，縱向力計算時取B0值：C = B0
D——巔因子，表示曲線的最大值：
BCD——縱向力零點處的縱向剛度：
B – 剛度因子：B=BCD/(C×D)
Sh——曲線的水平方向漂移：
Sv——曲線的垂直方向漂移：Sv=0
E——曲線曲率因子，表示曲線最大值附近的形狀：
圖 輪胎屬性文件中的縱向力計算系數數據塊
圖 Pacejka ’89輪胎縱向力示例
輪胎側向力計算公式為：
此時的X1為側向力計算組合自變量：X1=(α+Sh)，α為側偏角
C——曲線形狀因子，側向力計算時取A0值：C = A0
D——巔因子，表示曲線的最大值：
BCD——側向力零點處的側向剛度：
B – 剛度因子：B=BCD/(C×D)
Sh——曲線的水平方向漂移：
Sv——曲線的垂直方向漂移：
E——曲線曲率因子，表示曲線最大值附近的形狀：
圖 輪胎屬性文件中的側向力計算系數數據塊
圖 Pacejka ’89輪胎縱向力示例
輪胎回正力矩計算公式為：
此時的X1為回正力矩計算組合自變量：X1=(α+Sh)，α為側偏角
C——曲線形狀因子，回正力矩計算時取C0值：C = C0
D——巔因子，表示曲線的最大值：
BCD——回正力矩零點處的扭轉剛度：
B – 剛度因子：B=BCD/(C×D)
Sh——曲線的水平方向漂移：
Sv——曲線的垂直方向漂移：
E——曲線曲率因子，表示曲線最大值附近的形狀：
圖 輪胎屬性文件中的回正力矩計算系數數據塊
圖 Pacejka ’89輪胎回正力矩示例
側偏剛度（Lateral Stiffness）
側偏剛度在Pacejka ’89和’94輪胎模型中假定是一個常量，在輪胎屬性文件的參數PARAMETER數據段中通過LATERAL_STIFFNESS語句設定。
側向形變De：De＝Fy /LATERAL_STIFFNESS；
翻轉力矩：Mx ＝ -Fz ×De；
縱向力和側偏角聯合作用的回正力矩Mz；
MZ = MZ,MF + Fx×De ，這里MZ,MF為魔術公式計算所得的回正力矩。
滾動阻力（Rolling resistance）
滾動阻力系數RR同樣是在輪胎屬性文件中規定的具體值，滾動阻力矩My：
My = Fz ×Re×RR
這里：Re為輪胎的滾動半徑；RR為滾動阻力系數；Fz垂直載荷（kN）。
平滑過渡（Smoothing）
是否使用平滑過渡也在輪胎屬性文件中規定：
² USE_MODE = 1 或 2：關閉平滑過渡
² USE_MODE = 3 或 4：使用平滑過渡
輪胎屬性文件TR_rear_pac89.tir全文（示例整車模型MDI_Demo_Vehicle.asy使用的）：
$---------------------------------------------------------------------MDI_HEADER
[MDI_HEADER]
FILE_TYPE = 'tir'
FILE_VERSION = 2.0
FILE_FORMAT = 'ASCII'
(COMMENTS)
{comment_string}
'Tire - XXXXXX'
'Pressure - XXXXXX'
'Test Date - XXXXXX'
'Test tire'
'New File Format v2.1'
$--------------------------------------------------------------------------UNITS
[UNITS]
LENGTH = 'mm'
FORCE = 'newton'
ANGLE = 'radians'
MASS = 'kg'
TIME = 'sec'
$--------------------------------------------------------------------------MODEL
[MODEL]
! use mode 1 2 3 4
! -------------------------------------------
! smoothing X X
! combined X X
!
PROPERTY_FILE_FORMAT = 'PAC89' 輪胎模型關鍵詞
FUNCTION_NAME = 'TYR900' 解算器函數
USE_MODE = 4.0 平滑過渡模式
$----------------------------------------------------------------------DIMENSION
[DIMENSION]
UNLOADED_RADIUS = 340.6 輪胎自由半徑
WIDTH = 255.0 輪胎寬度
ASPECT_RATIO = 0.35 高寬比
$----------------------------------------------------------------------PARAMETER
[PARAMETER]
VERTICAL_STIFFNESS = 310.0 縱向剛度系數
VERTICAL_DAMPING = 3.1 縱向阻尼系數
LATERAL_STIFFNESS = 190.0 側偏剛度
ROLLING_RESISTANCE = 0.0 滾動阻力系數
$-----------------------------------------------------------LATERAL_COEFFICIENTS
[LATERAL_COEFFICIENTS]
a0 = 1.65000
a1 = -34.0
a2 = 1250.00
a3 = 3036.00
a4 = 12.80
a5 = 0.00501
a6 = -0.02103
a7 = 0.77394
a8 = 0.0022890
a9 = 0.013442
a10 = 0.003709
a11 = 19.1656
a12 = 1.21356
a13 = 6.26206
$-------------------------------------------------------------------longitudinal
[LONGITUDINAL_COEFFICIENTS]
b0 = 2.37272
b1 = -9.46000
b2 = 1490.00
b3 = 130.000
b4 = 276.000
b5 = 0.08860
b6 = 0.00402
b7 = -0.06150
b8 = 1.20000
b9 = 0.02990
b10 = -0.17600
$----------------------------------------------------------------------aligning
[ALIGNING_COEFFICIENTS]
c0 = 2.34000
c1 = 1.4950
c2 = 6.416654
c3 = -3.57403
c4 = -0.087737
c5 = 0.098410
c6 = 0.0027699
c7 = -0.0001151
c8 = 0.1000
c9 = -1.33329
c10 = 0.025501
c11 = -0.02357
c12 = 0.03027
c13 = -0.0647
c14 = 0.0211329
c15 = 0.89469
c16 = -0.099443
c17 = -3.336941
注意：屬性文件中的單位數據塊[UNITS]不用於魔術公式的系數a,b,c。
2. 代碼
MF5.2的實現代碼：
void RPacejka::CalcMF52() // Pacejka MF5.2 (~2006) { rfloat Fx0,Dx,Cx,Ex,Bx,dfz,Fz0,Fz0_der; // Nominal load, adapted nominal load rfloat kx,k,mux,sign_kx; rfloat Shx,Svx,Kx,gamma,gamma_x; rfloat Fzn; // Fz in Newtons rfloat tmp; // // Global // k=slipPercentage*0.01f; // Slip ratio Fzn=Fz*1000.0f; Fz0=fz0; //4500; Fz0_der=lfz0*Fz0; dfz=(Fzn-Fz0_der)/Fz0_der; gamma=camber/RR_RAD2DEG; // // Fx (pure slip) // Shx=(phx1+phx2*dfz)*lhx; Svx=Fzn*(pvx1+pvx2*dfz)*lvx*lmux; kx=k+Shx; if(kx<0)sign_kx=-1.0f; else sign_kx=1.0f; gamma_x=gamma*lgax; Cx=pcx1*lcx; mux=(pdx1+pdx2*dfz)*(1.0f-pdx3*gamma_x*gamma_x)*lmux; // Different in Pac2006 Ex=(pex1+pex2*dfz+pex3*dfz*dfz)*(1.0f-pex4*sign_kx)*lex; // Limiter on Ex (eq 23) if(Ex>1.0f)Ex=1.0f; Dx=mux*Fzn; Kx=Fzn*(pkx1+pkx2*dfz)*expf(pkx3*dfz)*lkx; // K=BCD (=stiffness) // Low velocity trouble tmp=Cx*Dx; if(fabs(tmp)<0.001f) Bx=Kx/(tmp+0.001f); else Bx=Kx/tmp; tmp=Bx*kx; Fx0=Dx*sinf(Cx*atanf(tmp-Ex*(tmp-atanf(tmp))))+Svx; // // Fy // rfloat Fy0,Dy,Cy,Ey,By; rfloat localMuy; rfloat Shy,Svy,Ky; rfloat say,sign_alpha_y; rfloat gamma_y; rfloat alpha; alpha=sideSlip/RR_RAD2DEG; gamma_y=gamma*lgay; Cy=pcy1*lcy; localMuy=(pdy1+pdy2*dfz)*(1.0f-pdy3*gamma_y*gamma_y)*lmuy; Dy=localMuy*Fzn; Ky=pky1*Fz0_der*sinf(2.0f*atanf(Fzn/(pky2*Fz0_der)))*(1.0f-pky3*fabs(gamma_y))*lfz0*lky; // =BCD=stiffness at slipangle 0 tmp=Cy*Dy; if(fabs(tmp)<0.001f) By=Ky/(tmp+0.001f); else By=Ky/tmp; Shy=(phy1+phy2*dfz)*lhy+phy3*gamma_y; Svy=Fzn*((pvy1+pvy2*dfz)*lvy+(pvy3+pvy4*dfz)*gamma_y)*lmuy; say=alpha+Shy; if(say<0)sign_alpha_y=-1.0f; else sign_alpha_y=1.0f; Ey=(pey1+pey2*dfz)*(1.0f-(pey3+pey4*gamma_y)*sign_alpha_y)*ley; // Limiter on Ey (eq 35) if(Ey>1.0f)Ey=1.0f; // Lateral base force tmp=By*say; Fy0=Dy*sinf(Cy*atanf(tmp-Ey*(tmp-atan(tmp))))+Svy; // // Mz // rfloat Mzr; rfloat t0,Dt,Ct,Bt,alpha_t,Et,gamma_z; rfloat Sht,Shf; rfloat Br,R0; rfloat alpha_r,Dr; const float lr=1.0f; // Not found in paper at lambda section R0=r0; gamma_z=gamma*lgaz; Sht=qhz1+qhz2*dfz+(qhz3+qhz4*dfz)*gamma_z; alpha_t=alpha+Sht; // Avoid division by zero if(fabs(Ky)<0.001f) if(Ky<0)Ky=-0.001f; else Ky=0.001f; Shf=Shy+Svy/Ky; alpha_r=alpha+Shf; Bt=(qbz1+qbz2*dfz+qbz3*dfz*dfz)*(1.0f+qbz4*gamma_z+qbz5*fabs(gamma_y))*lky/lmuy; // Pac2006 adds gamma_y^2 dependence (qbz6) Ct=qcz1; Dt=Fzn*(qdz1+qdz2*dfz)*(1.0f+qdz3*gamma_z+qdz4*gamma_z*gamma_z)*(R0/Fz0_der)*lt; // lt=lamba_t? Et=(qez1+qez2*dfz+qez3*dfz*dfz)*(1.0f+(qez4+qez5*gamma_z)*(2.0f/3.14159265f)*atanf(Bt*Ct*alpha_t)); // <=1 // Clamp Et (eq 51) if(Et>1.0f)Et=1.0f; Br=qbz9*lky/lmuy+qbz10*By*Cy; // Preferred qbz9=0 tmp=Bt*alpha_t; t0=Dt*cosf(Ct*atanf(tmp-Et*(tmp-atan(tmp))))*cosf(alpha); // t_alpha_t in the paper Dr=Fzn*((qdz6+qdz7*dfz)*lr+(qdz8+qdz9*dfz)*gamma_z)*R0*lmuy; // *cosf(alpha_t) for Pac2006 (in MF52 this is still in Mzr=... below) Mzr=Dr*cosf(Ct*atanf(Br*alpha_r))*cos(alpha); // last cos(alpha) is cos(alpha_t) in Pac2006 // No combined aligning moment Mz=-t0*Fy0+Mzr; // // Combined slip // // Combine (page 30+) // Longitudinal float Shxa,Exa,Cxa,Bxa,alpha_s,Gxa0,Gxa; Shxa=rhx1; Exa=rex1+rex2*dfz; Cxa=rcx1; Bxa=rbx1*cosf(atan(rbx2*k))*lxal; // +rbx3*gammaStar*gammaStar) (Pac2006) alpha_s=alpha+Shxa; Gxa0=cos(Cxa*atan(Bxa*Shxa-Exa*(Bxa*Shxa-atan(Bxa*Shxa)))); if(fabs(Gxa0)>D3_EPSILON) Gxa=cos(Cxa*atan(Bxa*alpha_s-Exa*(Bxa*alpha_s-atan(Bxa*alpha_s))))/Gxa0; else Gxa=0; // Or 1 for super grip? Fx=Gxa*Fx0; // Lateral float Dvyk,Svyk,Shyk,Eyk,Cyk,Byk,ks,Gyk0,Gyk; Dvyk=muy*Fz*(rvy1+rvy2*dfz+rvy3*gamma)*cosf(atanf(rvy4*alpha)); Svyk=Dvyk*sin(rvy5*atan(rvy6*k))*lvyka; Shyk=rhy1+rhy2*dfz; Eyk=rey1+rey2*dfz; Cyk=rcy1; Byk=rby1*cosf(atanf(rby2*(alpha-rby3)))*lyka; ks=k+Shyk; Gyk0=cosf(Cyk*atanf(Byk*Shyk-Eyk*(Byk*Shyk-atanf(Byk*Shyk)))); Gyk=cosf(Cyk*atanf(Byk*ks-Eyk*(Byk*ks-atanf(Byk*ks))))/Gyk0; Fy=Gyk*Fy0+Svyk; // Aligning torque float alpha_r_eq,alpha_t_eq,Mzr,Fy_der; float sign_alpha_t,sign_alpha_r; float kk,s,t; if(alpha_t>=0)sign_alpha_t=1.0f; else sign_alpha_t=-1.0f; if(alpha_r>=0)sign_alpha_r=1.0f; else sign_alpha_r=-1.0f; kk=Kx/Ky; kk=(kk*kk*k*k); // kk^2*k^2 alpha_r_eq=sqrtf(alpha_r*alpha_r+kk)*sign_alpha_r; alpha_t_eq=sqrtf(alpha_t*alpha_t+kk)*sign_alpha_t; s=(ssz1+ssz2*(Fy/Fz0_der)+(ssz3+ssz4*dfz)*gamma)*R0*ls; Mzr=Dr*cosf(atanf(Br*alpha_r_eq))*cosf(alpha); Fy_der=Fy-Svyk; // New pneumatic trail tmp=Bt*alpha_t_eq; t=Dt*cosf(Ct*atanf(tmp-Et*(tmp-atanf(tmp))))*cosf(alpha); // Add all aligning forces Mz=-t*Fy_der+Mzr+s*Fx; // Postprocess; negate for Racer?! Fy=-Fy; Mz=-Mz; // Static results latStiffness=-By*Cy*Dy; // There's that '-' again for lateral force (Fy) longStiffness=Bx*Cx*Dx; }
MF6.1實現代碼：
rfloat RPacejka::CalcMz96() // Calculates aligning moment { rfloat Mz; rfloat B,C,D,E,Sh,Sv; rfloat FzSquared; // Calculate derived coefficients FzSquared=Fz*Fz; C=c0; D=c1*FzSquared+c2*Fz; E=(c7*FzSquared+c8*Fz+c9)*(1.0f-c10*fabs(camber)); if((C>-D3_EPSILON&&C<D3_EPSILON)|| (D>-D3_EPSILON&&D<D3_EPSILON)) { B=99999.0f; } else { B=((c3*FzSquared+c4*Fz)*(1-c6*fabs(camber))*expf(-c5*Fz))/(C*D); } Sh=c11*camber+c12*Fz+c13; Sv=(c14*FzSquared+c15*Fz)*camber+c16*Fz+c17; Mz=D*sinf(C*atanf(B*(1.0f-E)*(sideSlip+Sh)+ E*atanf(B*(sideSlip+Sh))))+Sv; return Mz; }
// Fx - longitudinal rfloat RPacejka::CalcFx96() // Pacejka96 model // Calculates longitudinal force (Fx) // From G. Genta's book, page 63 // Note that the units are inconsistent: // Fz is in kN // slipRatio is in percentage (=slipRatio*100=slipPercentage) // camber and slipAngle are in degrees // Resulting forces are better defined: // Fx, Fy are in N // Mz is in Nm { rfloat B,C,D,E; rfloat Fx; rfloat Sh,Sv; rfloat uP; rfloat FzSquared; // Calculate derived coefficients FzSquared=Fz*Fz; C=b0; uP=b1*Fz+b2; D=uP*Fz; // Avoid div by 0 if((C>-D3_EPSILON&&C<D3_EPSILON) || (D>-D3_EPSILON&&D<D3_EPSILON)) { B=99999.0f; } else { B=((b3*FzSquared+b4*Fz)*expf(-b5*Fz))/(C*D); } E=b6*FzSquared+b7*Fz+b8; Sh=b9*Fz+b10; Sv=b11*Fz+b12; // Notice that product BCD is the longitudinal tire stiffness longStiffness=B*C*D; // *RR_RAD2DEG; // RR_RAD2DEG == *360/2PI // Remember the max longitudinal force (where sin(...)==1) maxForce.x=D+Sv; // Calculate result force Fx=D*sinf(C*atanf(B*(1.0f-E)*(slipPercentage+Sh)+E*atanf(B*(slipPercentage+Sh))))+Sv; return Fx; }
// Fy - lateral rfloat RPacejka::CalcFy96() // Pacejka 96 // Calculates lateral force // Note that B*C*D is the cornering stiffness, and // Sh and Sv account for ply steer and conicity forces { rfloat B,C,D,E; rfloat Fy; rfloat Sh,Sv; rfloat uP; // Calculate derived coefficients C=a0; uP=a1*Fz+a2; D=uP*Fz; E=a6*Fz+a7; // Avoid div by 0 if((C>-D3_EPSILON&&C<D3_EPSILON) || (D>-D3_EPSILON&&D<D3_EPSILON)) { B=99999.0f; } else { if(a4>-D3_EPSILON&&a4<D3_EPSILON) { B=99999.0f; } else { // Notice that product BCD is the lateral stiffness (=cornering) latStiffness=a3*sinf(2*atanf(Fz/a4))*(1.0f-a5*fabs(camber)); B=latStiffness/(C*D); } } Sh=a8*camber+a9*Fz+a10; Sv=(a111*Fz+a112)*camber*Fz+a12*Fz+a13; // Remember maximum lateral force maxForce.y=D+Sv; // Calculate result force Fy=D*sinf(C*atanf(B*(1.0f-E)*(sideSlip+Sh)+ E*atanf(B*(sideSlip+Sh))))+Sv; return Fy; }