题目：在Matlab 环境下，推导半车模型的运动方程并进行响应分析。
假设车身是刚体，并具有垂直和俯仰两个自由度，要求得出半车模型的系统响应功率谱密度和响应函数。
推导过程：
设车身质心垂向振动量为b Z ，俯仰角为b θ，则有如下运动方程组：
r
f b hp r f b hb r tr wr f f tf wf bF aF I F F z m Fr z z K z
m F z z K z m +-=+=--=--=θ )()(303101
将b b b b b z z a z z θθ+=-=42,以及
)()()()(43432121z z
C z z K F z z C z z K F sr sr r sf sf f -+-=-+-=
代入上面运动方程组，将4321z z z z 、、、作为状态变量，即][43
2
1
z z z z Z =，
]0*0*[00r
tr f tf z K z K B =，可得到状态方程：
''B AZ =
矩阵B 中含有前轮输入和后轮输入的耦合，假设后轮所行驶的轨迹为前轮在时间段t ∆前行驶过的轨迹，且
u L t /=∆，既有如下关系：
t i f e Z t Z z ωω000)sin(==
u L i f u L i t i r e z e Z u L t Z t t Z z /0/0000)/sin()](sin[ωωωωωω--==-=∆-=
因此，可将B 改写为
f u
L i tr tf
z e K K B 0/*]00[ω-=
而矩阵A 通过化简，得到其为：
⎥⎥
⎥⎥⎥⎦⎤
⎢⎢
⎢⎢⎢
⎣
⎡++-+-++---+++-++-++-+---+++-)()
()
()(00)()
()()(0
323222221212sr sr sr sr sf sf sf sf sr
sr sr
sr tr wr sr sr sr sr sf sf sf sf sf sf sf
sf tf wf C i K n C i K n C i K n C i K n C i K C i K K m C i K n C i K n C i K n C i K n C i K C i K K m ωωωωωωωωωωωωωωωω
其中，hp
hb hp hb hp hb I b m n I ab m n I a m n 2
32211,1,1+
=-=+= 感想： 这次作业比较难，所以很多地方借鉴了别人的想法。
我觉得主要难点在于把时域模型转换成频域模型，十分抽象不好理解，这涉及到傅里叶变换，采用现代
控制理论的状态空间表达式，建立状态方程进行求解。
在MATLAB 中编程调试了很久，才和书上的图例基本保持一致，验证了我的推导和分析是正确的。
示意图来拍于书上，同时从规范性考虑，上述推导过程中所列公式均采用公式编辑器完成。
附件： 含有Matlab 代码的M 文件
半车模型的系统响应功率谱密度和响应函数（共5张图）
Clc
Clear all
close all
mhb=690;
Ihp=1222;
mwf=40.5;
mwr=45.4;
Ktf=192e3;
Ktr=192e3;
Ksf=17e3;
Ksr=22e3;
Csf=1.5e3;
Csr=1.5e3;
a=1.25;b=1.51;L=2.76;v=20;
a1=1/mhb+a^2/Ihp ;b1=1/mhb-a *b/Ihp ;c1=1/mhb+b^2/Ihp;
for i=1:150
f(i)=i/10;%频率
omega=2*pi*f(i);
s=(4.47e-4)/(f(i)^2.5);%路面输入
A=[Ktf+Ksf-omega^2*mwf+1j*Csf*omega,-Ksf-1j*Csf*omega,0,0;
a1*(Ksf+1j*Csf*omega),omega^2+a1*(-Ksf-1j*Csf*omega),b1*(Ksr+1j*Csr*omega),b1*(-Ksr-1j*Csr*omega);
0,0,-omega^2*mwr+Ktr+Ksr+1j*omega*Csr,-Ksr-1j*Csr*omega;
b1*(Ksf+1j*Csf*omega),b1*(-Ksf-1j*Csf*omega),c1*(Ksr+1j*Csr*omega),c1*(-Ksr-1j*Csr*omega)+omega^2];
B=[Ktf;0; Ktr *exp(-1j* omega* L/v);0];
Z=A\B;
f_susp_pmd(i)=s*abs(Z(2)-Z(1))^2;%前悬架动行程功率谱密度
r_susp_pmd(i)=s*abs(Z(4)-Z(3))^2;%后悬架动行程功率谱密度
f_a(i)=(-omega)^2*abs(Z(2));%求前轮簧载质量加速度增益
r_a(i)=(-omega)^2*abs(Z(4));%求后轮簧载质量加速度增益
cg_a(i)=abs(-omega^2*Z(2) *b-omega^2*Z(4)*a)/L;%求质心处垂向加速度增益
f_a_pmd(i)=s*f_a(i)^2;%前轮簧载质量垂向加速度功率谱密度
r_a_pmd(i)=s*r_a(i)^2;%后轮簧戟质量垂向加速度功率谱密度
cg_a_pmd(i)=s*abs(cg_a(i))^2;%质心簧载质量垂向加速度功率谱密度
f_tire_pmd(i)=s*(Ktf/1000*abs(Z(1)-1))^2;%前轮轮胎动载荷功率谱密度
r_tire_pmd(i)=s*(Ktr/1000*abs(Z(3)-exp(-1j* omega* L/v)))^2;%后轮轮胎动载荷功率谱密度thb(i)=abs(atan(Z(4)-Z(2))/L)*180/pi;%求俯仰角增益
if f(i)<1
jq=10^(-0.6);
elseif f(i)>=1&&f(i)<=4
jq=10^((3*log2(f(i))-6)/10);
elseif f(i)>4&&f(i)<8
jq=1;
else
jq=10^((-6*log2(f(i))+18)/10);%加权系数
end
jqf_a_pmd(i)=jq* f_a_pmd(i);%经加权后的前轮垂向载荷功率谱密度
jqr_a_pmd(i)=jq*r_a_pmd(i);
jqcg_a_pmd(i)=jq* cg_a_pmd(i);
end
figure(1)
plot(f,f_susp_pmd,'b',f,r_susp_pmd,'--r');legend('前轮','后轮');
xlabel('频/Hz');
ylabel('悬架动行程功率谱密度(m^2/Hz');
figure(2)
plot(f,jqf_a_pmd,'b',f,jqr_a_pmd,'--r',f,jqcg_a_pmd,':k');legend('前轮','后轮','质心'); xlabel('频率Hz');
ylabel('垂向加速度功率谱密度/[(m/s^2)Hz]');
figure(3)
plot(f, f_tire_pmd,'b', f, r_tire_pmd,'--r');legend('前轮','后轮');
xlabel('频率Hz');
ylabel('轮胎动载荷功率谱密度/(kN^2/Hz)');
figure(4)
plot(f,cg_a);
xlabel('频率Hz');
ylabel('质心处垂向加速度增益/[(m/s^2)/m]');
figure(5)
plot(f,thb);
xlabel('频率Hz');
ylabel('俯仰角增益/[(°)/m]');。