【悬架笔记零】半车垂向振动模型

模型如下
半车振动模型，被动悬架
力学推导

代码

%书第135页
clc
clear
close all
%% 1.车辆参数
ms=690;
Isy=1222; 
mwf=40.5; 
mwr=45.4;
Ksf=17000; 
Ksr=22000;
Csf=1500; 
Csr= 1500; 
Kwf= 192000; 
Kwr= 192000; 
a= 1.25; 
b= 1.51;
L=a+b;
%% 2.状态方程
A=[0 0 0 0 0 0 1 0 0 0;...
    0 0 0 0 0 0 0 1 0 0;...
    0 0 0 0 0 0 0 0 1 0;...
    0 0 0 0 0 0 0 0 0 1;...
    0 0 0 0 0 0 0 0 0 0;...
    0 0 0 0 0 0 0 0 0 0;...
    -1*(Ksf+Ksr)/ms (a* Ksf-b*Ksr)/ms Ksf/ms Ksr/ms 0 0 -1*(Csf+Csr)/ms (a*Csf-b* Csr)/ms Csf/ms Csr/ms;...
    (a* Ksf-b* Ksr)/Isy -1*(a*a* Ksf+b* b * Ksr)/Isy -1*a*Ksf/Isy b*Ksr/Isy 0 0 (a*Csf-b*Csr)/Isy -1*(a*a*Csf+b*b*Csr)/Isy -1*a*Csf/Isy b* Csr/Isy;...
    Ksf/mwf -1*a*Ksf/mwf -1*(Ksf+Kwf) 0 Kwf/mwf 0 Csf/mwf -1*a*Csf/mwf -1*Csf/mwf 0;...
    Ksr/mwr b*Ksr/mwr 0 -1*(Ksr+Kwr)/mwr 0 Kwr/mwr Csr/mwr b*Csr/mwr 0 -1*(Csr/mwr)];
    
B=[0 0 0 0 1 0 0 0 0 0;0 0 0 0 0 1 0 0 0 0]';

C=[-1*(Ksf+Ksr)/ms (a*Ksf-b*Ksr)/ms Ksf/ms Ksr/ms 0 0 -1*(Csf+Csr)/ms (a*Csf-b*Csr)/ms Csf/ms Csr/ms ; ...
(a*Ksf-b*Ksr)/Isy -1*(a*a*Ksf+b*b*Ksr)/Isy -1*a* Ksf/Isy b*Ksr/Isy 0 0 (a*Csf-b* Csr)/Isy -1*(a*a*Csf+b*b*Csr)/Isy -1*a*Csf/Isy b*Csr/Isy;...
1 -1*a -1 0 0 0 0 0 0 0;...
1 b -1 0 0 0 0 0 0 0;...
0 0 -1*Kwf 0 Kwf 0 0 0 0 0;...
0 0 0 -1*Kwr 0 Kwr 0 0 0 0];
%% 3.仿真设置
dt=0.01;%仿真间隔时间
t_end=10;%仿真总时长
t=0:dt:t_end;
u=60;%车速,km/h
u=u/3.6;%车速转化为m/s

n00=0.011;%下截止频率
f0=2*pi*n00*u;
Sq=256*10^-6;%路面等级系数
white_noise=1-2*1*rand(1,length(t)-1);%白噪声数据
for i=1:length(t)-1
    q(1)=0;
    q(i+1)=dt*(-2*pi*10*q(1)+2*pi*(sqrt(Sq*u))* white_noise(i))+q(1);
    dq=diff(q)/dt;
end
dq(2,:)=dq(1,:);
%% 4.仿真求解
for i=1:length(t)-1
    x(:,1)= [0;0;0;0;0;0;0;0;0;0];
    [G,H]=c2d(A,B,0.01);
    x(:,i+1)=G*x(:,i)+H*dq(:,i); 
    y(:,i)=C*x(:,i); 
end
%% 5.结果绘图
figure(1)
t1=0.01:0.01:10;
plot(t1,y(1,:))
xlabel('时间/s')
ylabel('车身垂直加速度/(m/s^2)')

figure(2)
t1=0.01:0.01:10;
plot(t1,y(2,:))
xlabel('时间/s');
ylabel("车身俯仰角加速度/(rad/s^2)")

figure(3)
tl=0.01:0.01:10;
plot(t1,y(3,:))
xlabel('时间/s');
ylabel('前悬架动挠度/m')

figure(4)
tl=0.01:0.01:10;
plot(t1,y(4,:))
xlabel('时间/s')
ylabel('后悬架动挠度/m')

figure(5)
t1=0.01:0.01:10;
plot(t1,y(5,:))
xlabel('时间/s')
ylabel('前轮胎动载荷/N')

figure(6)
t1=0.01:0.01:10;
plot(t1,y(6,:))
xlabel('时间/s')
ylabel('后轮胎动载荷/N')

效果（例）