最新机械动力学第一次大作业
机械动力学大作业
曲柄连杆机构的建模与仿真
曲柄连杆机构如图1所示，曲柄AC长90mm,OC距离300mm。
计算β=300时曲柄的ｖ、θ、θ’。
1-1仿真测试
通过ADAMS/View建立的如图1-1-1所示的模型：
图1-1-1曲柄连杆模型
（1）在工具栏中单击仿真控制按钮。
（2）系统打开参数设置对话框，设置为Forever，Steps为0.01，如图1-1-2所示。
图1-2运动仿真
（3）点击开始仿真图标
（4）模型开始运动。
如果曲柄逆时针转动，对旋转运动（Motion_1）单击鼠标右键，选择Modify，在Function（time）一栏的数据前面加上“-”。
开始仿真。
1-2测试仿真结果
（1）在ADAMS工具菜单的Measure中选择，在To Point栏中输入PART3，即曲柄；
（2）在Characteristic栏中选择Translation velocity 以测其速度。
（3）选择Cylindrical坐标系，并选择R选项。
（4）设置完毕单击Apply按钮，弹出测量窗口如图1-2-1曲柄连杆运行一周的速度分量图。
图1-2-1 qubingliangan 测量曲线
（5）重复（（1）~（3）步再测量节点2、3沿Y轴方向的速度分量，测量结果如图1-2-2：
图1-2-2 qubingliangan2、3测量曲线。
机械动力学大作业
机械动力学大作业
一、问题及要求
建立单自由度杆机构（有无滑块均可）动力学模型，由静止启动，选择一个固定驱动力矩，绘制原动件在一周内的运动关系线图，具体
机构及参数自拟。
、建立模型
建立如图一所示机构：已知驱动力矩M=20Nm阻力矩M=10 Nm 杆1长120mm转动惯量为J i=0.16kgm2 ;杆3质量为1kg，转动惯量为
j3=0.16kgm2;杆2质量、转动惯量忽略，其他参数如图所示。
求：建立系统运动方程。
二、求解过程
M3图一
方法：利用等效力学模型法进行动力学分析
解：取杆1为等效件，有
f w3
M v= Mi - M3—（1）
jv=ji+02+j3）（—y ⑵
< XWj/
w3120 3
—二—=一=0.75 Wi 160 4
又有（2）可知J V为常数则可知：
Ki
M v=J v i
则错误！未找到引用源。
四、采用ADAMS软件或Matlab/Simulink 环境，建立机械系统的动力学模型，借助软件进行求解计算和结果分析。
（1）利用Adams软件，建模后如图：
图2
图3
(2)当杆1由图1所示位置开始运动一周，机构运动时间为0.03秒, 利用Adam漱件分析杆件1角加速度错误！未找到引用源。
随时间的变化关系图，如图4所示：
图 线系关动运的内周一在件动原
4
图
uoqB 」曰|总3。
\/」e[n6uv
0/SS3
y
s s C6唾
二 N 2
I 」善E。
机械动力学作业
仿真运动结束后的如图1-5所示。
图1-5仿真运动结束后的图
建立测量（滑块的位移、速度、加速度）
1）鼠标右键单击需要测量的部件，系统打开右键快捷菜单，选择Measure；
2）系统打开参数对话框，将Characteristic设为CM Position，Component设为X，测量X向位移；
图1-1曲柄滑块机构
本次所使用的运动仿真软件是ADAMS2012版，启动ADAMS/View程序，打开ADAMS软件，进行新建模型，然后在检查和设置建模基本环境，准备工作都完成后，就开始进行几何建模，建模完成图如下图1-2所示。
进行施加运动副和驱动，施加铰接副，在曲柄同地面框架连接，和曲柄连杆之间的连接，以及在连杆和滑块上处通过铰接副将连杆和滑块连接。滑块和地面也要进行棱柱副连接。运动副施加完成后的图，如图1-3所示。
3）点击Apply，出现空白的测量窗口；
4）重复上述步骤，将Characteristic设为CM Velocity，新建测量速度；
5）重复上述步骤，将Characteristic设为CM Acceleration，新建测量加速度；
测量后测得结果如图1-6,1-7,1-8所示。
图1-6滑块位移图
图1-7滑块速度图
图1-3运动副施加完成后的图
定义曲柄的运动
①在主工具箱的运动工具集，选择旋转运动工具图标 ，显示定义旋转运动对话框；
②在Set up栏，输入360；选择JOINT_1，完成转速设置。
定义曲柄运动后，如图1-4所示。
图1-4定义曲柄运动
对曲柄滑块机构进行仿真分析
1）仿真分析
机械原理课程大作业
机械原理课程大作业基于MATLAB平面连杆机构运动学和动力学分析指导老师：王玉丹目录作业一：平面连杆机构运动学分析第2页作业二：平面连杆机构动力学分析第15页作业一L(AE)=70mm,L(AB)=40mm,L(EF)=60mm,L(DE)=35mm,L(CD)=75m m,L(BC)=50mm,原动件以等角速度W1=10rad/s回转。
试以图解法求在θ1=50°时C点的速度和加速度.对机构进行运动分析，写出C点的位置、速度及加速度方程。
解题过程：令AB=r1, BC=r2, CD=r3, DE=r4，AE=r6，EF=r8, AF=r7，角EAF=θ1。
分析：对机构进行位置分析由封闭形ABCDEA可得：r1+r2=r6+r3+r4 （1）由封闭图形AEFA可得：r7=r6+r8 （2）将（1）（2）两式整理可得：r2-r3-r4=-r1+r6-r8+r7=r6【一】（1）位置方程：【二】速度方程：【三】加速度方程：【四】根据位置方程式编制如下函数：【五】进行数据输入，运行程序进行运算。
根据上面分析的θ1 的极限位置取θ1 的范围为40°-55°并均分成15个元素：输出的P、矩阵的第二列到第四列分别是θ2 、θ3 、4θ4 的值，第一列是AF杆的长度r1’。
【六】第二步根据速度方程式编写如下函数：根据第一步得到的数据进行数据输入，运行程序计算各速度值。
程序如下：程序运行得到q矩阵，第一行到第三行分别是a2、a3、a4 的值，第四行是杆AF上滑块运动的速度，即F点的速度。
【七】第三步编写加速度计算函数：【八】根据第一步和第二步输入数据，运行程序得到各加速度的值：【1】计算C点在θ1 =55°，w1 =10rad/s时的速度，加速度：总结数据绘出各构件的位置、速度和加速度的表格如下：【2】输出图像1）角位置程序及输出的图像：2）F点速度程序及输出的图像：3）角加速度程序及输出的图像：4）F点的加速度程序及输出图像：作业二在图示的正弦机构中，已知：L(AB)=100mm,h1=120mm,h2=80mm, W1=10rad/s(常数)，滑块2和构件3的重量分别为，G2 =40 N 和G3 =100 N,质心S2 和S3 的位置如图所示，加于构件3上的生产阻力Fr=400 N，构件1的重力和惯性力略去不计。
石大远程奥鹏-机械原理-第一次在线作业正确答案
错误
正确答案：正确
33.局部自由度是指有些机构中某些构件所产生的不影响其他构件的局部运动的自由度。
正确
错误
正确答案：正确
34.虚约束是指机构中某些运动副带入的对机构运动起重复约束作用的约束。虚约束虽对机构的运动并产生约束作用，但会导致机构自由度的计算结果与机构的实际自由度不相符。
正确
错误
正确答案：错误
A、无穷远
B、中心
C、法线
D、接触点
正确答案：D
25.三心定理是指三个彼此作平面运动的构件的（）个瞬心必位于同一直线上。
A、一
B、二
C、三
D、四
正确答案：C
26.机器都是由各种材料做成的制造单元零件经装配而成的各个运动单元零件的组合体；各个运动单元之间具有确定的相对运动。
正确
错误
正确答案：错误
27.构件和零件是两个不同的概念：构件是运动单元；零件是制造单元；机器中的构件可以是单一的零件，也可以是由几个零件装配成的刚性结构。
A、产生
B、不产生
C、主要
D、关键
正确答案：B
12.为了便于对含有高副的平面机构进行分析研究，可以将机构中的高副根据一定的条件虚拟的以低副加以代替，这种代替的方法就叫做（）。
A、低副高代
B、低副低代
C、高副高代
D、高副低代
正确答案：D
13.高副低代代替前后机构的（）完全相同。
A、低副数
B、高副数
C、构件数
D、自由度
正确答案：D
14.在平面机构中进行高副低代时，为了使得代替前后机构的自由度、瞬时速度和加速度都保持不变，只要用一个（）分别与两高副低代构件在过接触点的曲率中心处以转动副相联就行了。
机械原理大作业一 33题
机械原理大作业（一）作业名称：机械原理设计题目：连杆机构运动分析院系：机电工程学院班级： xxxxxx设计者： xxx学号： xxxxxx指导教师： xxxxx设计时间： 6.25---7.1哈尔滨工业大学机械设计一，运动分析题目二，建立以点G为原点的固定平面直角坐标系G-x,y三，对机构进行结构分析该机构由原动件AB（Ι级组),EFG(RRRⅡ级杆组），EHK（RRRⅡ级杆组),ECD(RRRⅡ级杆组)组成。
四，各基本杆组的运动运动分析数学模型（1）原动件AB（Ⅰ级组）已知原动件AB的转角ψ1=0~2π原动件AB的角速度ω1=10rad/s原动件AB的角加速度α1=0运动副A的位置坐标 XA=0 YA=0A点与机架相连，即该点速度和加速度均为0。
运动副A的速度 VxA=0 VyA=0运动副A的加速度 aXA=0 aYA=0原动件AB长度 lAB=200mm可求出运动副B的位置坐标 XA=XA+lABcosψ1 YB=YA+lABsinψ1 运动副B的速度 vXB=vXA-ω1lABsinψ1 vYB=vYA+ω1lABcosψ1运动副B的加速度aXB=aXA-ω12lABcosψ1-α1lABsinψ1aYB=aYA-ω12lABsinψ1+α1lABcosψ1（2）ECD （RRR Ⅱ级杆组）由（1）知B 点位置坐标、速度、加速度运动副D 点位置坐标 XD=XA+lADcos ψ2 YD=YA+lADsin ψ2D 点与机架相连，即该点速度和加速度均为0。
运动副D 的速度 v X D=0 v Y D=0运动副D 的加速度 a X D=0 a Y D=0杆BC 长 lBC=800mm杆CD 长 lCD=448mm可求得BC 杆相对于X 轴正方向转角Ψ3=2arctan （)00/()0202020(C A C B A B +-++）CD 杆相对于x 轴正方向转角Ψ4=arctan （（YC-YD ）/（XC-XD ）)其中，A0=2lBC(XD-XB)B0=2lBC(YD-YB)C0=l ２BC+l ２BD-l ２CDl ２BD=(XD-XB)２+(YD-YB)２求导可得BC 杆ω3、α3和CD 杆ω4、α4最后求导得vXC 、vYC 以及aXC 、aYCC 的轨迹即是E 的轨迹（3）EFG （RRR Ⅱ级杆组）运动副G点位置坐标 X G=XA+lAGcosψ5 YG=YA+lAGsinψ5G点与机架相连，即该点速度和加速度均为0。
机械动力学试验报告
机械动力学大作业院系：机械工程学院专业：机械设计及理论班级：4班学号：20404151姓名：王昌荣日期：2005.4.20机械动力学课程大作业（一）实验任务：1、观察激振试验；2、进行锤击试验，绘制特性曲线图及数据表，并利用所得数据对模态参数进行识别。
（二）仪器设备1、GF－10型功率放大器；2、XD-18A型超低频信号发生器；3、MO-1251型示波器；4、PDM-2000数据采集分析仪；5、带有传感器的力锤；6、加速度传感器。
（三）实验框图锤击实验框图：1．绘制实频、虚频、幅频、相频及Nyqiust圆源程序清单：%原始数据读入------------------------------------------------->clc;for j=1:4fid=0;while fid<1switch jcase 1filename='wcry01.txt';case 2filename='wcry02.txt';case 3filename='wcry03.txt';case 4filename='wcry04.txt';end[fid,j]=fopen(filename,'r');if fid==-1disp(message)endenddata_all=fscanf(fid,'%g');N=length(data_all); %数据长度fs=2560; %采样频率f=1000; %分析频率data_Re=data_all(1:N/2); %前半段数据为实频data_Im=data_all(N/2+1:N); %后半段数据为虚频for i=1:N/2Aw(i)=sqrt(data_Re(i)*data_Re(i)+data_Im(i)*data_Im(i)); %幅频Qw(i)=atan(data_Im(i)/data_Re(i)); %相频endfiguresubplot(221);plot(data_Re(1:400),'b');title('—实频特性—')subplot(222);plot(data_Im(1:400),'g');title('—虚频特性—')subplot(223);plot(Aw(1:400),'c');title('—幅频特性—')subplot(224);plot(Qw(1:400),'r');title('—相频特性—')figureplot(data_Re(1:400),data_Im(1:400),'k');title('—Nyquist特性—')end第一组数据绘图：第二组数据绘图：第三组数据绘图第四组数据绘图:2．用五点求拟合圆，并计算各阶固有频率，阻尼比：采用第一组数据计算：源程序清单：clc;fprintf(‘----------采用第1组数据计算----------\n')fid=0;while fid<1filename='wcry01.txt'; %原始数据文件名称 [fid,message]=fopen(filename,'r'); %读入数据if fid==-1disp(message);endenddata_all=fscanf(fid,'%g');N=length(data_all); %数据长度%fs=2560; %采样频率%f=1000; %分析频率data_Re=data_all(1:N/2); %前半段数据为实频 data_Im=data_all(N/2+1:N); %后半段数据为虚频 for i=1:N/2Aw(i)=sqrt(data_Re(i)*data_Re(i)+data_Im(i)*data_Im(i)); %幅频endpeak_value=zeros(1,3);subpoint=[1 40 91 130 131 170];[max_data1,peak_value(1)]=max(Aw(subpoint(1):subpoint(2)));[max_data2,peak_value(2)]=max(Aw(subpoint(3):subpoint(4)));[max_data3,peak_value(3)]=max(Aw(subpoint(5):subpoint(6)));peak_value(2)=peak_value(2)+subpoint(3)-1;peak_value(3)=peak_value(3)+subpoint(5)-1;fprintf('\n 三阶峰值位置分别是：');fprintf('%d %d %d\n',peak_value(1),peak_value(2),peak_value(3));m=5;k=(m+1)/2;dataX_curve=zeros(1,m);dataY_curve=zeros(1,m);for j=1:3peakva=peak_value(j);[center_X,center_Y,curve_R,A_mainarray]=calc_curve(data_Re,data_Im,peakva); ring_X=zeros(1,628);ring_Y=zeros(1,628);for i=1:628ring_X(i)=curve_R*cos(i/100)+center_X;ring_Y(i)=curve_R*sin(i/100)+center_Y;endfigureplot(ring_X,ring_Y,'m');hold onfor i=1:mdataX_curve(i)=data_Re(peak_value(j)-k+i);dataY_curve(i)=data_Im(peak_value(j)-k+i);endplot(dataX_curve,dataY_curve,'b-*'); %绘制拟合圆所取各点switch jcase 1title('一阶五点拟合园法');case 2title('二阶五点拟合园法');case 3title('三阶五点拟合园法');endx1=ring_X(1)+curve_R/10; %x轴绘制y1=0;x2=ring_X(314)-curve_R/10;y2=0;point_X=[x1,x2];point_Y=[y1,y2];plot(point_X,point_Y,'k','LineWidth',1);y1=ring_Y(157)+curve_R/10; %y轴绘制x1=0;y2=ring_Y(471)-curve_R/10;x2=0;point_X=[x1,x2];point_Y=[y1,y2];plot(point_X,point_Y,'k','LineWidth',1);plot(center_X,center_Y,'mo'); %圆心坐标x1=ring_X(1)+curve_R/10; %圆对称轴xy1=center_Y;x2=ring_X(314)-curve_R/10;y2=y1;point_X=[x1,x2];point_Y=[y1,y2];plot(point_X,point_Y,'r--','LineWidth',0.8);y1=ring_Y(157)+curve_R/10; %圆对称轴yx1=center_X;y2=ring_Y(471)-curve_R/10;x2=x1;point_X=[x1,x2];point_Y=[y1,y2];plot(point_X,point_Y,'r--','LineWidth',0.8);hold off;axis('image');chord_angle=zeros(1,6);chord_anglem=0;for i=2:m+1chord_angle1=atan(data_Im(peak_value(j)-k+i-1)/data_Re(peak_value(j)-k+i-1));chord_angle2=atan(data_Im(peak_value(j)-k+i )/data_Re(peak_value(j)-k+i ));chord_angle(i) =abs(chord_angle2-chord_angle1);end[chord_anglem,peak_num]=max(chord_angle(1:5)); %取得最大角度即为最大弦长位置chord_num=peak_value(j)-k+peak_num;hold onif peak_num==5 %如果固有频率为第5点则绘制下一点 x1=data_Re(peak_value(j)+k-1);y1=data_Im(peak_value(j)+k-1);x2=data_Re(peak_value(j)+k);y2=data_Im(peak_value(j)+k);point_X=[x1,x2];point_Y=[y1,y2];plot(point_X,point_Y,'b-*');endfor i=1:3 %绘制固有频率与圆心连线x1=center_X;y1=center_Y;x2=data_Re(peak_value(j)-k+peak_num-i+2);y2=data_Im(peak_value(j)-k+peak_num-i+2);point_X=[x1,x2];point_Y=[y1,y2];plot(point_X,point_Y,'b','LineWidth',0.8);endnatural_FREQ(j)=(peak_value(j)-k+peak_num)*2.56; %得到各阶固有频率pkvl_plus =(peak_value(j)+1)+(peak_value(j)-1);pkvl_divid=(peak_value(j)+1)-(peak_value(j)-1);%得到各阶阻尼比damping_ratio(j)=(pkvl_plus-2*k+2*peak_num)*pkvl_divid/2/...(peak_value(j)-k+peak_num)^2*(tan(chord_angle(peak_num)/2)+tan(chord_angle(peak_num+1)/ 2));endfprintf('\n 各阶固有频率分别是：'); %输出各阶固有频率disp(natural_FREQ);fprintf('\n 各阶阻尼比分别是：'); %输出各阶阻尼比disp(damping_ratio);程序清单：function [center_X,center_Y,curve_R,A_mainarray]=calc_curve(data_Re,data_Im,peakva) m=5;k=(m+1)/2;for i=1:mdataX_curve(i)=data_Re(peakva-k+i);dataY_curve(i)=data_Im(peakva-k+i);endcurve_A=zeros(3,3);curve_B=zeros(3,1);curve_C=zeros(3,1);curve_A(1,1)=dataX_curve*dataX_curve';curve_A(1,2)=dataX_curve*dataY_curve';curve_A(2,2)=dataY_curve*dataY_curve';for i=1:mcurve_A(1,3)=curve_A(1,3)+dataX_curve(i);curve_A(2,3)=curve_A(2,3)+dataY_curve(i);curve_B(1 )=curve_B(1)-dataX_curve(i)^3-dataX_curve(i)*dataY_curve(i)^2;curve_B(2 )=curve_B(2)-dataX_curve(i)^2*dataY_curve(i)-dataY_curve(i)^3;curve_B(3 )=curve_B(3)-dataX_curve(i)^2 -dataY_curve(i)^2;endcurve_A(2,1)=curve_A(1,2);curve_A(3,1)=curve_A(1,3);curve_A(3,2)=curve_A(2,3);curve_A(3,3)=m;A_mainarray=curve_A;curve_C=curve_A\curve_B;center_X=-curve_C(1)/2;center_Y=-curve_C(2)/2;curve_R=sqrt(center_X^2+center_Y^2-curve_C(3));%子程序结束结果输出----------采用第1组数据计算---------- 三阶峰值位置分别是： 14 326 361各阶固有频率分别是： 35.8400 837.1200 921.6000各阶阻尼比分别是： 1.2261 0.0025 0.0017三阶拟合圆图形：2．求各阶主振型，仍然调用上述求圆心子程序程序清单：clc;file_num=4;main_array=zeros(file_num,3);for ks=1:file_numfid=0;while fid<1switch kscase 1filename='wcry01.txt';case 2filename='wcry02.txt';case 3filename='wcry03.txt';case 4filename='wcry04.txt';end[fid,w]=fopen(filename,'r');if fid==-1disp(message)endenddata_all=fscanf(fid,'%g');N=length(data_all); %数据长度data_Re=data_all(1:N/2); %前半段数据为实频 data_Im=data_all(N/2+1:N); %后半段数据为虚频 for i=1:N/2Aw(i)=sqrt(data_Re(i)*data_Re(i)+data_Im(i)*data_Im(i)); %幅频endpeak_value=zeros(1,3);subpoint=[1 40 91 130 131 170];[max_data1,peak_value(1)]=max(Aw(subpoint(1):subpoint(2)));[max_data2,peak_value(2)]=max(Aw(subpoint(3):subpoint(4)));[max_data3,peak_value(3)]=max(Aw(subpoint(5):subpoint(6)));peak_value(2)=peak_value(2)+subpoint(3)-1;peak_value(3)=peak_value(3)+subpoint(5)-1;m=5;k=(m+1)/2;for j=1:3peakva=peak_value(j);[center_X,center_Y,curve_R,A_mainarray]=calc_curve(data_Re,data_Im,peakva);peak_flag=1;if data_Im(peakva)<0;peak_flag=-1;endmain_array(ks,j)=curve_R*peak_flag;endendarray_max=zeros(1,3);for i=1:3array_max(i)=main_array(4*i-2);endfor i=1:3for j=1:4main_array(j,i)=main_array(j,i)/array_max(i);endendfor i=1:3figureplot(main_array((4*i-3):(4*i)),'m');hold offswitch icase 1title('一阶主阵型');case 2title('二阶主阵型');case 3title('三阶主阵型');endend%程序结束绘图及计算结果：的现代化方法和手段。
2021年机械系统动力学大作业
《机械系统动力学》课程作业小组组员: 王凌飞 0702081t王毅 0702041指导老师: 何荇兮学院: 机械工程学院专业: 机械工程重庆大学机械工程学院十一月机械系统动力学大作业一、 问题描述图1为汽车结构简化模型:图1 汽车结构简化模型图2为汽车结构受力分析:图2 受力分析图已知22120.64m 4000kg 2000N s/m r m c c ====⋅121220000N/m0.9m1.4m k k l l ====r : 车辆回转半径。
初始条件为: 0x x θθ====。
外部冲击力矩: )(10t δ。
试用MATLAB 中ode45函数求解并画出0-5s 内位移x 和转角θ响应。
单位冲击函数()t δ定义:1,0()0,t t t δ=⎧=⎨≠⎩, 其图像如图3所表示。
00.511.5tδ图3 单位冲击函数图像二、 求解过程1.系统运动方程不考虑冲击力矩, 由图2机构受力分析得到系统运动方程以下:0)()()()(112221112221=-+++-+++θθk l k l x k k c l c l x c c x m (1) 0)()()()(222121112212122211222=++-+++-+θθk l k l x l k l k c l c l x l c l c x mr (2) 考虑t=0时刻, 系统受到一个冲击力矩)(10t δ, 此时运动学方程表示为:0)()()()(112221112221=-+++-+++θθk l k l x k k c l c l x c c x m (3) 222222211221122112211()()()()10()mr c l c l x l c l c k l k l x l k l k t θθθδ+-+++-++= (4)2．运动方程一阶常微分方程组形式令，，，，4321x x x x x x ====θθ 则 t=0时:12212222114121221133422222422112221142211111223()/()/()/()/[10()/()()()]/x x x c c x m l c l c x m k k x m l k l k x m x x x l c l c x l c l c x l k l k x l k l k x mr =⎧⎪=-+---+--⎪⎨=⎪⎪=---+---+⎩ t>0时:12212222114121221133422222422112221142211111223()/()/()/()/[()/()()()]/x x x c c x m l c l c x m k k x m l k l k x m x x x l c l c x l c l c x l k l k x l k l k x mr =⎧⎪=-+---+--⎪⎨=⎪⎪=---+---+⎩ 3. MATLAB 程序求解运动方程ode45函数主程序文件solution.m, 以下所表示: clc;clear;r=0.8; m=4000; c1=; c2=; k1=0; k2=0; l1=0.9; l2=1.4; t0=0; tf=5;x0=[0,0,0,0];[t,x]=ode45('f1',[t0 tf],x0); subplot(3,1,1) plot(t,x);subplot(3,1,2) plot(t,x(:,1)); subplot(3,1,3) plot(t,x(:,3));ode45函数微分关系函数文件f1.m, 以下所表示:function xdot=f1(t,x) r=0.8; m=4000; c1=; c2=; k1=0; k2=0; l1=0.9; l2=1.4; t0=0; tf=5;xdot=zeros(4,1); xdot(1)=x(2);xdot(2)=-(c1+c2)/m*x(2)-(l2*c2-l1*c1)/m*x(4)-(k1+k2)/m*x(1)-(l2*k2-l1*k1)/m*x(3);xdot(3)=x(4); if t<=0xdot(4)=10/(m*r*r)-(c2*l2-c1*l1)/(m*r*r)*x(2)-(c2*l2*l2+c1*l1*l1)/(m*r*r)*x(4)-(k2*l2-k1*l1)/(m*r*r)*x(1)-(l1*l1*k1+l2*l2*k2)/(m*r*r)*x(3); elsexdot(4)=-(c2*l2-c1*l1)/(m*r*r)*x(2)-(c2*l2*l2+c1*l1*l1)/(m*r*r)*x(4)-(k2*l2-k 1*l1)/(m*r*r)*x(1)-(l1*l1*k1+l2*l2*k2)/(m*r*r)*x(3); end在MATLAB 中运行以上两个程序, 即可得0～5s 内汽车系统位移x 与转角θ响应。