二自由度振动仿真：matlab直接解微分方程&&virtual.lab motion仿真

一、二自由度振动解析法与多种数值算法见：https://blog.csdn.net/lijil168/article/details/67647924

 

二、用matlab符号建立微分方程，并用龙格库塔法 ode45()直接进行数据计算

%SH760小轿车空载主要参数
clear;
m=1340;
a=1.54;
b=1.29;
l=a+b;
Ic=2395; %绕质心的转动惯量
rou=sqrt(Ic/m);
k1=40*1000;
k2=44*1000;
M=[m*(b^2+rou^2)/l^2,m*(a*b-rou^2)/l^2;m*(a*b-rou^2)/l^2,m*(a^2+rou^2)/l^2];
K=[k1,0;0,k2];
 
%用matlab ode45数值解------------------------------------------------
 
A=[zeros(2,2),eye(2);-M\K,zeros(2,2)]; %四阶参数4X4矩阵X'=AX-->X=expm(A*t)*X0 X=[x1;x2;x1';x2']
syms x1 x2 dx1 dx2
df_sym=A*[x1;x2;dx1;dx2];
df_sym=subs(df_sym,{'x1','x2','dx1','dx2'},{'x(1)','x(2)','x(3)','x(4)'});
n=length(df_sym);
i=1;
ss='[';   %先定义好很重要，否则再循环体中定义时，每一循环ss不累加。
while i<n
    ss=strcat(ss,char(df_sym(i)),';');
    i=i+1;
end
ss=strcat(ss,char(df_sym(i)));
ss=strcat(ss,']');
f=inline(ss,'t','x');
[t,x]=ode45(f,[0 10],[1,0,0,0]);%初始
%subplot(2,2,1)
plot(t,[x(:,1),x(:,2)]) %时间状态系列

三、virtual.lab motion仿真

注意，当设置初始位移x1=1m,x2=0m时，由方程x1=x-l1*theta,x2=x+l2*theta,计算小车重心的初始位移x=0.455830和小车的初始转角theta=0.353357时，有误差（线性近似计算）。

由motion制图容易取得真实的值：x_3d=0.456000，theta_3d=0.361039（20.686度）。

另外，力平衡方程也是线性近似的（只有在小位移时成立）。这是造成理论计算和motion仿真结果有部分误差的主要原因！

特别需要注意的是motion用欧拉四元素来定义小车的初始姿态，欧拉四元素可参考维基百科：https://en.wikibooks.org/wiki/Multibody_Mechanics/Euler_Parameters

u为旋转轴的单位矢量，phai 为相对旋转轴的转角。欧拉四元素组成四维单位矢量，即四元素平方和=1。

为了方便得到各参数，用计算器软件先计算：

仿真结果为：

可见，小车的运动趋势和理论解很接近，而且理论微分方程是近似的，在大位移时误差会偏大，用motion仿真的可行度更高。