魔术公式轮胎建模

同济大学汽车学院1133090使用魔术公式的轮胎模型一、魔术公式简介魔术公式是用三角函数的组合公式拟合轮胎试验数据，用一套形式相同的公式就可以完整地表达轮胎的纵向力Fx、侧向力Fy、回正力矩Mz、翻转力矩Mx、阻力矩My以及纵向力、侧向力的联合作用工况，故称为“魔术公式”。魔术公式的一般表达式为：BxBxEBxCDxYarctanarctansin式中Y(x)可以是侧向力，也可以是回正力矩或者纵向力，自变量x可以在不同的情况下分别表示轮胎的侧偏角或纵向滑移率，式中的系数B、C、D依次由轮胎的垂直载荷和外倾角来确定。在侧向加速度常见范围≤0.4g，侧偏角≤5°的情景下对常规轮胎具有很高的拟合精度。此外，由于魔术公式基于试验数据，除在试验范围的高精度外，甚至在极限值以外一定程度仍可使用，可以对有限工况进行外推且具有较好的置信度。基于魔术公式的轮胎模型的输入和输出变量二、建立模型本程序根据魔术公式的基本结构，在给定参数的基础上，在合理范围内拟定了一部分参数，得到了仿真结果，从总体上来看复合魔术公式的结果和实际情况。模型程序见附件。拟定参数：****0.91110;11xExHxyCxCy1111111KyaKyrtxykvyks[10,10][0.5,0.5]k输入变量纵向滑移率κ侧偏角α[rad]外倾角γ[rad]车轮垂直载荷FZ[kN]输出变量纵向力FX[N]侧向力FY[N]翻转力矩MX[Nm]滚动阻力矩MY[Nm]回正力矩MZ[Nm]魔术公式同济大学汽车学院1133090三、仿真结果分析根据所建模型进行仿真，得仿真结果如下图所示。-1-0.8-0.6-0.4-0.200.20.40.60.81-5000-4000-3000-2000-1000010002000300040005000纵向滑移率纵向力/（N）纵向力(纯纵滑)-10-9-8-7-6-5-4-3-2-1012345678910-4000-3000-2000-100001000200030004000侧偏角/（度）侧偏力/（N）侧偏力(纯侧滑)-10-9-8-7-6-5-4-3-2-1012345678910-80-70-60-50-40-30-20-1001020304050607080侧偏角/（度）回正力矩/（N*m）回正力矩(纯侧滑)-10-9-8-7-6-5-4-3-2-1012345678910-5000-4000-3000-2000-1000010002000300040005000侧偏角/（度）纵向力/（N）纵向力(组合滑移)-1-0.8-0.6-0.4-0.200.20.40.60.81-3000-2000-10000100020003000纵向滑移率侧向力/（N）侧偏力(组合滑移)-10-9-8-7-6-5-4-3-2-1012345678910-60-50-40-30-20-100102030405060侧偏角/（度）回正力矩/（N*m)回正力矩(组合滑移)1、纯侧滑工况从纵向力-纵向滑移率曲线可以看到，力与滑移率之间具有很强的非线性关系；但可以看出当轮胎的滑移率低于3%时，轮胎所传递的纵向力和滑移率之间近似成线性关系。而在车辆正常行驶中，轮胎的滑移率很少超过2%，因此在这一范围内可以定义纵向刚度()xvRFCv其中F为纵向力，xC为纵向刚度。当纵向滑移率超过3%时，纵向力Fx随滑移率的增大呈非线性关系增大，当纵向滑移率增大到15%左右时，纵向力达到最大值；此后滑移率继续增大，纵向力随之减小，最终趋于稳定值，约为0.7Fz。从侧偏力-侧偏角曲线可以看到，由于正的侧偏力产生负的侧偏角，因此曲线和魔术公式中的曲线y值刚好相反，假如在画图时画成-Fy则结果和魔术公式中的一致。在侧偏角接近0度的范围内，侧偏力和侧偏角同样存在近似的线性关系，同理可以求得侧偏刚度。1时的力称为测偏功率pC。侧偏角继续增大时，侧偏力-侧偏角曲线的斜率逐渐减小，这时轮胎在接地处已经发生部分侧滑。最后侧偏力达到附着极限时，整个轮胎侧滑。同济大学汽车学院1133090从回正力矩-侧偏角曲线可以看到,回正力矩在起始阶段逐步增大，在侧偏角达到4度左右时达到最大值，而后减小。回正力矩的最大值随轮胎与路面间的摩擦系数的变化而变化。在侧偏角0时的斜率称为回正刚度，它是关系到转向是否沉重、前轮摆动等的重要因素。2、组合滑移工况从纵向力-侧偏角曲线可以看到，在侧偏角小于2度的范围内，纵向力随着侧偏角的增大迅速增大，一旦侧偏角超过2度，纵向力便开始减小，并最终趋于稳定值。从纵向滑移率-侧偏力曲线上可以看出，在纵向滑移率接近0的小区间内，滑移率的值变化不会引起测偏力的变化，测偏力始终等于0；当滑移率不断增大时，测偏力开始跟随滑移率呈非线性增长。从侧偏角-回正力矩曲线图上可以看出，在侧偏角小于1度的范围内，回正力矩迅速增长，且在侧偏角约为1度的时候达到最大值，随后随着侧偏角的增大回正力矩减小。3、总结从纯滑移和组合滑移的仿真结果对比可以看出，两者的曲线趋势基本上是一致的，区别只在于横坐标的初始阶段两者稍有区别。纵坐标为测偏力时，测偏力随横坐标的变化稍有滞后；纵坐标为纵向力或者回正力矩时，纵坐标上升的速度更快。附件：模型文件%%MagicFormulamodel%ParametervaluesofMagicFormulainAppendix3Fz0=4000;%nominal(rated)load(0）,NR0=0.313;%unloadedtyreradius(=ro)，mFz=4000;%verticalload,Nr=0.1*pi/180;%camberanglePhx1=-0.002;Phx2=0.002;Pdx1=1.21;Pdx2=-0.037;Pvx1=0;Pvx2=0;Pex1=0.344;Pex2=0.095;Pex3=-0.02;Pex4=0;Pkx1=21.51;Pkx2=-0.163;Pkx3=0.245;Pcx1=1.685;Pcy1=1.193;Pdy1=-0.99;Pdy2=0.145;Pdy3=-11.23;Pey1=-1.003;Pey2=-0.537;Pey3=-0.083;Pey4=-4.787;Pey5=1;Pky1=-14.95;Pky2=2.13;Pky3=-0.028;Pky4=2;Pky5=0;Pky6=-0.92;Pky7=-0.24;Phy1=0.003;Phy2=-0.001;Phy3=0;Pvy1=0.045;Pvy2=-0.024;Pvy3=-0.532;Pvy4=0.039;同济大学汽车学院1133090qcz1=1.18;qhz1=0.007;qhz2=-0.002;qhz3=0.147;qhz4=0.004;qbz1=8.964;qbz2=-1.106;qbz3=-0.842;qbz5=-0.227;qbz6=0;qbz9=18.47;qbz10=0;qdz1=0.1;qdz2=-0.001;qdz3=0.007;qdz4=13.05;qdz6=-0.008;qdz7=0;qdz8=-0.296;qdz9=-0.009;qdz10=0;qdz11=0;qez1=-1.609;qez2=-0.359;qez3=0;qez4=0.174;qez5=-0.896;rcx1=1.092;rhx1=0.007;rex1=0.5;rex2=0.5;%setrbx1=12.35;rbx2=-10.77;rbx3=0;rcy1=1.081;rhy1=0.009;rhy2=0.1;%setrby1=6.461;rby2=4.196;rby3=-0.015;rby4=0;rvy1=0.053;rvy2=-0.073;rvy3=0.517;rvy4=35.44;rvy5=1.9;rvy6=-10.71;rey1=0.5;rey2=0.5;%setSsz1=0.043;Ssz2=0.001;Ssz3=0.731;Ssz4=-0.238;%LongitudinalForce(purelongitudinalslip)%input:slipratiokk=linspace(-1,1,200);dfz=(Fz-Fz0)/Fz0;ux=(Pdx1+Pdx2*dfz)*0.9;%取lam(ux*)=0.9Kxk=Fz*(Pkx1+Pkx2*dfz)*exp(Pkx3*dfz);%取lam(Ex*)=1Shx=Phx1+Phx2*dfz;%取lam(Hx*)=1kx=k+Shx;Ex=(Pex1+Pex2*dfz+Pex3*dfz^2)*(1-Pex4*sign(kx));%取lam(Ex*)=1Dx=ux*Fz;%caita(1-8)都取1Cx=Pcx1;%取lamCx=1Bx=Kxk/(Cx*Dx);%取伊布西隆=0Svx=Fz*(Pvx1+Pvx2*dfz);%后面的修正系数全部取1Fx0=Dx*sin(Cx*atan(Bx*kx-Ex.*(Bx*kx-atan(Bx*kx))))+Svx;subplot(2,3,1)plot(k,Fx0);gridset(gca,'xlim',[-11])%设置x轴范围set(gca,'xtick',[-1:0.2:1]);%设置x轴间隔set(gca,'ylim',[-50005000])%设置x轴范围set(gca,'ytick',[-5000:1000:5000]);%设置x轴间隔xlabel('纵向滑移率');ylabel('纵向力/（N）');title('纵向力(纯纵滑)');%LateralForce(puresideslip)%input:横向侧偏是侧偏角取tan后的值，即a*a=linspace(-10,10,200);ta=linspace(-0.17633,0.17633,200);%侧偏角取-10-10度，tan后变成-0.17633-0.17633，还是200个点同济大学汽车学院1133090sr=sin(r);%r是外倾角，sr表示r*uy=(Pdy1+Pdy2*dfz)/(1+Pdy3*sr^2);%取lam(uy*)=1Cy=Pcy1;%取lamCy=1Dy=uy*Fz;Kya=Pky1*Fz0*sin(Pky4*atan(Fz/((Pky2+Pky5*sr^2)*Fz0)))/(1+Pky3*sr^2);%取lam（Kya）=1By=Kya/(Cy*Dy);%取伊布西隆y=0；Svyr=Fz*(Pvy3+Pvy4*dfz)*sr;%后面的修正系数全部取1Svy=Fz*(Pvy1+Pvy2*dfz)+Svyr;%后面的修正系数全部取1Kyr0=Fz*(Pky6+Pky7*dfz);%取lam（Kyr)=1Shy=(Phy1+Phy2*dfz)+(Kyr0*sr-Svyr)/Kya;ay=ta+Shy;%tan(sa)表示横向侧偏角a*Ey=(Pey1+Pey2*dfz)*(1+Pey5*sr^2-(Pey3+Pey4*sr)*sign(ay));Fy0=Dy*sin(Cy*atan(By*ay-Ey.*(By*ay-atan(By*ay))))+Svy;subplot(2,3,2)plot(a,Fy0);gridset(gca,'xlim',[-1010])set(gca,'xtick',[-10:1:10]);set(gca,'ylim',[-40004000])set(gca,'ytick',[-4000:1000:4000]);xlabel('侧偏角/（度）');ylabel('侧偏力/（N）');title('侧偏力(纯侧滑)');%AligningTorque(puresideslip)%input:侧偏角Cr=1;%所有caita都取1Ct=qcz1;Bt=(qbz1+qbz2*dfz+qbz3*dfz^2)*(