4汽车模型非线性分析.docx

基于R-K法的1/4汽车模型非线性分析 摘 要:文章采用RK法4阶方法解微分方程,在求解当中,加入非线性判断语句,分析车轮离地等造成的非线性工况。给出了线性系统和非线性系统图形的差别。简要叙述了频谱图,得出了解决此类非线性问题的一种简单可行的办法。 下载 关键词:微分方程组;MATLAB;非线性分析;1/4汽车模型 中图分类号:U461.4文献标识码:A 文章编号:1006-8937(2010)08-0076-02 matlab能轻易的建立LTI对象,能够满足绝大部分的连续线性时不变系统仿真的需求。但是由于汽车悬架的阻尼,刚度等不是一个常数,甚至车轮离地等,所以要考虑非线性问题。文章使用龙格-库塔法解微分方程组,加入非线性判断语句,从而更真实的模拟悬架垂线运动。程序输出了轮胎离地情况下的图形,与线性系统输出图像做个比较。 1 龙格库塔法编程 1.1 R-K编程实现 1/4双质量振动数学模型为以下方程组 MY+DY+KY=R Y=M-1R-KY-DY 其中M为质量矩阵,D为阻尼矩阵,K为刚度矩阵。 利用四阶龙格-库塔解微分方程(经典R-K),直接用数值方法来解微分方程。matlab自带ode45函数,但是其系数,不能根据情况随时变化,不能体现解决非线性的微分方程组。下面是编程实现: (参数设置:k1=604685,k2=45482,m1=98, m2=500,r=2500) clear 输入参数区域 n=10000;%仿真步总数 h=linspace(0,1,n); %时间序列 INY=0,0;0,0; %初始状态设置。1行为非簧载质量位移和速度,2行为簧载质量 U=(x)0.1; %仿真输入函数,模拟时域的路面文件。 m1=98;c1=604685; %模型的参数 m2=500;c2=45482;r=2500; 初始化参数 K0=c1+c2,-c2;-c2,c2; %K0为初始刚度矩阵,D0为 初始阻尼矩阵 D0=r,-r;-r,r; %在仿真中,如果满足非线性条件,那么它们会改变。 M=m1,0;0,m2; %质量矩阵,也可变化 dR=(m1+m2)*9.8/c1;%轮胎静平衡时的形变量。 z=zeros(n,1); %参数用于提取轮胎离地高度 y=zeros(2,n);v=zeros(2,n);a=zeros(2,n);%分别是位移,速度和加速度。 dh=h(2)-h(1);%步长 for i=1:n 进入循环,载入初值 K=K0;D=D0; if i=1 Y=INY(:,1);d1=INY(:,2);X=Y; else d1=V; end 非线性的判断语句 % if casestatement ; end if Y(1)-U(h(i)dR %非簧载质量位移大于输入位移,即轮胎离地 c1=(m1+m2)*9.8/(Y(1)-U(h(i);%轮胎离地时转换的等效的刚度 K(1,1)=c1+c2; z(i)=Y(1)-U(h(i)-dR;%轮胎离地高度 end R=c1*U(h(i);0;%根据c1更改输入 以下是解微分方程组 %RK需要4点的斜率,第1点的速度和位移,已经在初始化载入。 dd1=M(R-K*X-D*d1); %(1点加速度) d2=d1+dh/2*dd1; %(2点速度) X=Y+dh/2*(d1+d2)/2; %(2点位移)改进EULER法? dd2=M(R-K*X-D*d2); %(2点加速度) d3=d1+dh/2*(dd1+dd2)/2; %(3点速度) X=Y+dh/2*(d1+4*d2+d3)/6;%(3点位移)3阶KUTTA法 dd3=M(R-K*X-D*d3); %(3点加速度) d4=d1+dh*(dd1+4*dd2+dd3)/6;%(4点速度)3阶KUTTA X=Y+(d1+2*d2+2*d3+d4)/6*dh;%(n+1时刻位移状态量)4阶KUTTA dd4=M(R-K*X-D*d4); %(4点加速度) A=(dd1+2*dd2+2*dd3+dd4)/6;% (n+1时刻加速度) V=d1+A*dh;% (n+1时刻速度)4阶KUTTA Y=X; %赋值Y,作为下一个循环的初值 下面是提取该循环的状态量 a(:,i)=A;v(:,i)=V;y(:,i)=Y; end plot(h,y,h,z+0.1); legend(非簧载位移,簧载位移,轮胎离地高度); grid on 可以输出加速度a,速度v,位移y,轮胎离地距离z,以及悬架绕度,轮胎附着力等等。程序屏蔽非线性判断语句时,与simulink线性模型得出的图形是一致的。在程序的“非线性判断”部分可以添加对非线性问题处理的语句。例如阻尼与速度,位移等状态量的函数关系。还有轮胎刚度非线性变化,及轮胎跳离地面等。可以是数值型对应关系,可见灵活性极大。 1.2 阻尼非线性 阻尼在被压缩时一般比伸长时小,比如r1=2500,r2=5000。那么程序选择部分就可以这样写:Y(2)-Y(1)0,悬架被拉伸,r1=5000。此外可以用数值方法输入任何想要的阻尼曲线。 1.3 轮胎离地引起的非线性 关于轮胎离地问题,就是Y(1)-U(t)dR。R0为轮胎不受力(包括不受重力)时的半径,dR为R0受mg作用时径向形变量,dR=mg/k1。由于原来系统平衡时,系统的重力是由轮胎的弹性力Ft来平衡。我们假象的弹簧力F=Ft-mg,向上为正。平衡位置为0。现在随着变形量减小,弹性力Ft减小,相当于我们假想的弹簧被拉伸,F减小。如果Y(1)-U(t)dR关系式满足,Ft=0那么假象的弹簧力F=-(m1+m2)g。之后,就算Y(1)-U(t)再增大,F也不变了。可以把轮胎的刚度k1分成两段,在Y(1)-U(t)dR,k1=k0,当Y(1)-U(t)dR阶段,k1=mg/(Y(1)-U(t)。则可满足F不再变化的需求。 2 图形比较 2.1 瞬态响应比较 图1对结果有一个感性的认识。该工况为汽车以极快的速度驶上一个为0.1 m高的台阶的情况。轮胎离地最高达0.08 m左右,大概在0.03 s,0.18s离地,此结果与速度加速度时域图形很吻合。 2.2 频率特性比较 将r设置为0,U=0,Y=1,0;1,0,屏蔽非线性判断语句,经过FFT处理可以得出系统的固有频率大概是1.5Hz。那么我们取幅值0.1,频率1.5 Hz的谐波来观察下两个系统的频率响应的差别。改输入为: U=(t)0.1*sin(2*pi*1.5*t);。 在固有频率1.5Hz处,非线性系统与线性系统有很大的差别。线性系统的谱密度图集中在1.5Hz处。非线性系统除了簧载质量(sprung)的功率谱密度受1.5Hz影响大,非簧载质量(unsprung)的功率谱非常分散。说明考虑轮胎离地情况下,输出加速度不仅与受迫振动的频率有关,而且还有输入的幅值有关。FFT处理,对非线性系统不是很管用了。此外,可以看到在固有频率处引起了共振,在1.5Hz处,簧载质量的加速度比非簧载质量还要大。 在5Hz(图2)处非簧载加速度大约是簧载加速度的5倍,在15Hz处非簧载加速度是簧载的20倍,可见,远离一阶固有频率可以极大的降低簧载加速度。并且在固有频率附近,簧载质量还会产生很多高频振荡,对人身体造成伤害。 此外令INY=1,0;1,0就可以得出汽车从1米高跌落,可以得出垂向加速度和悬架挠度图。此外还可以轻易的画出轮胎附着力图,可以用来分析汽车在某路面的安全性。 3 结 论 可以看出线性系统和非线性系统的差别。这种差别在汽车遇到恶劣条件的时候,尤其明显。线性系统输出与输入的谐波频率对应,而非线性系统则没有这样的性质。所以非线性系统没有如线性系统那样的传递函数,bode图等,用解析方法研究很困难,而数值方法比较适合。非线性系统的输出不仅与输入的频率有关,还与其幅值大小等等有关。 4 拓 展 汽车是一个非常复杂的系统。各个环节相互关联,这里只是对1/4模型仿真,为了更加详细的仿真模型,可以将以上的代码制作为simulink的s函数,那么在simulink里面很容易实现代码的重用。经过一定的