1、基于 R-K 法的 1/4 汽车模型非线性分析 摘 要 :文章采用 RK法 4阶方法解微分方程 ,在求解当中 ,加入非线性判断语句 ,分析车轮离地等造成的非线性工况。给出了线性系统和非线性系统图形的差别。简要叙述了频谱图 ,得出了解决此类非线性问题的一种简单可行的办法。 下载 关键词 :微分方程组 ;MATLAB;非线性分析 ;1/4 汽车模型 中图分类号 :U461.4 文献标识码 :A 文章编号 :1006-8937(2010)08-0076-02 matlab 能轻易的建立 LTI 对象 ,能够满足绝大部分的连续线性时不变系统仿真的需求。但是由于汽车悬架的阻尼 ,刚度等不是一个常数 ,甚
2、至车轮离地等 ,所以要考虑非线性问题。文章使用龙格 -库塔法解微分方程组 ,加入非线性判断语句 ,从而更真实的模拟悬架垂线运动。程序输出了轮胎离地情况下的图形 ,与线性系统输出图像做个比较。 1 龙格库塔法编程 1.1 R-K编程实现 1/4 双 质量振动数学模型为以下方程组 MY+DY+KY=R Y=M-1R-KY-DY 其中 M 为质量矩阵 ,D 为阻尼矩阵 ,K 为刚度矩阵。 利用四阶龙格 -库塔解微分方程 (经典 R-K),直接用数值方法来解微分方程。 matlab 自带 ode45 函数 ,但是其系数 ,不能根据情况随时变化 ,不能体现解决非线性的微分方程组。下面是编程实现 : (参
3、数设置 :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; %在仿真中 ,如果满
4、足非线性条件 ,那么它们会改变。 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
5、(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
6、*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+
7、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线性模型得出的图形是一致的。在程序的 “ 非线性判断 ” 部分可以添加对非线性问题处理的语句。例如阻尼与速度 ,位移等状态量的函数关系。还有轮胎刚度非线性变化 ,
8、及轮胎跳离地面等。可以是数值型对应关系 ,可见灵活性极大。 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。现在随着变
9、形量减小 ,弹性力 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
10、离地 ,此结果与速度加速度时域图形很吻合。 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)的功率谱非常分散。说明 考虑轮胎离地情况下
11、,输出加速度不仅与受迫振动的频率有关 ,而且还有输入的幅值有关。 FFT 处理 ,对非线性系统不是很管用了。此外 ,可以看到在固有频率处引起了共振 ,在 1.5Hz 处 ,簧载质量的加速度比非簧载质量还要大。 在 5Hz(图 2)处非簧载加速度大约是簧载加速度的 5倍 ,在 15Hz处非簧载加速度是簧载的 20 倍 ,可见 ,远离一阶固有频率可以极大的降低簧载加速度。并且在固有频率附近 ,簧载质量还会产生很多高频振荡 ,对人身体造成伤害。 此外令 INY=1,0;1,0就可以得出汽车从 1 米高跌落 ,可以得出垂 向加速度和悬架挠度图。此外还可以轻易的画出轮胎附着力图 ,可以用来分析汽车在某路
12、面的安全性。 3 结 论 可以看出线性系统和非线性系统的差别。这种差别在汽车遇到恶劣条件的时候 ,尤其明显。线性系统输出与输入的谐波频率对应 ,而非线性系统则没有这样的性质。所以非线性系统没有如线性系统那样的传递函数 ,bode 图等 ,用解析方法研究很困难 ,而数值方法比较适合。非线性系统的输出不仅与输入的频率有关 ,还与其幅值大小等等有关。 4 拓 展 汽车是一个非常复杂的系统。各个环节相互关联 ,这里 只是对 1/4模型仿真 ,为了更加详细的仿真模型 ,可以将以上的代码制作为 simulink的 s函数 ,那么在 simulink 里面很容易实现代码的重用。经过一定的完善就可以建立1/2 模型等。可以单独的研究某些问题 ,也可以为 adams 等大型仿真软件提供基本的参考。 参考文献 : (德 )H-P 威鲁麦特 .车辆动力学 :模拟及其方法 M.北京 :北京理工大学出版社 ,1998. 颜庆津 .数值分析 M.北京 :北京航空航天大学出版社 ,2006. 陈立平 .机械系统动力学分析及 ADAMS应用教 程 M.北 京 :清华大学出版社 ,2005.
