汽车平顺性仿真分析程序

1、clear; %清空 clc;clf;q=0； %定义初始变量jT;f_l=0.05;%给最小频率赋值f_u=50; %给最大频率赋值ns=1000;%将0.5-50HZ的频率分为1000段delta_f=(f_u-f_l)/ns;% 计算频率变化量Gq_n0=256*1e-6;%选取B级路面不平度等级n0=0.1;%给参考空间频率赋值u=40/3.6;%速度换算md=2*pi*rand(1001,1);%得出1001个随机变量矩阵for t=0:0.01:100%计算前轮路面输入i=1;for n=1:1:1000f_k=f_l+(n-0.5)*delta_f;Gq_f=Gq_n0*n0A2

2、*u/f_kA2;a_k=sqrt(2*Gq_f*delta_f);q=q+a_k*sin(f_k*2*pi*(t+md(i);i=i+1;enda(j)=t;%将时间变量赋值给变量ap(j)=q;%将输入变量赋值给pj=j+1;q=0;endfigure(1);plot(a,p); %作图grid on;%给图形加上栅格线xlim(0,100);ylim(auto);xlabel(时间(s);% 标注ylabel(路面不平度位移(m);title(路面时域随机数据(三角级数法);%对模拟随机输入数据进行功率谱分析，采用矩形 窗%fs分析频率,nf一频域的采样点数Py0,f=periodogr

3、am(p(1:500),hann(500),500,100,onesided); %利用周期图法求前 500 个点(即前5s)的功率谱Pya=Py0;for i=2:(2*20-1)Py0,f=periodogram(p(i-1)*250+1:(i-1)*250+500),hann(500),500,100,onesided); % 依次每隔250个点(即2.5s)求一次功率谱Pya=Pya+Py0;endPya=Pya/(2*20-1); % 求平均 figure(2);%建立图形窗口 2% subplot(212); loglog(f(1:500/2),Pya(1:500/2),r);%用

4、对数坐标作图grid on; xlim(0.1,50); xlabel(频率(Hz);% 标注ylabel (路面位移功率谱密度(mA2/Hz);hold on %将下面的图形一起绘制在figure(2)中 for i=1:500%计算路面随机输入理论上的功率谱密度fb(i)=0.1*i; Sq(i)=Gq_n0*0.01*u/(fb(i)A2);endloglog(fb,Sq,-);%采用对数坐标作图%根据轴距lw和车速V(m/s)生成后轮输入时域数据 lw=6.53;% 给轴距赋值=td=lw/u;%计算前后轮时间差l=round(td/0.01);%对 round 里面的数取整pr=p(

5、1:(10000-l); pf=p(l+1):10000);m1=775; m2=2140; m3=19702; m5=1653; m7=730; m9=51; i4=238000; i6=270; i8=380; k1=1.96E+6; k2=7.84E+6; k3=3.82E+5; k4=1.78E+6; k5=3.14E+6; k6=2.45E+6; k7=1.96E+4; k8=1.96E+4; k9=18121; c3=3.18E+4; c4=1.17E+5; c5=5.00E+3; c6=3.64E+3; c7=7.5E+3; c8=7.3E+3; c9=362*0.75;11=5

6、.42;12=(6.53-5.42);13=(5.42+1.18);14=(5.42+0.32);15=(5.42-0.65);16=(5.42-0.64);a=(1.18-0.54);b=(0.64+0.54);c=(0.21+0.32);d=(0.65-0.21);m=zeros(9,9);m(1,1:9)=m1;m(2,1:9)=m2;m(3,1:9)=m3;m(4,1:9)=i4;m(5,1:9)=m5;m(6,1:9)=i6;m(7,1:9)=m7;m(8,1:9)=i8;m(9,1:9)=m9;k=zeros(9,9);k(1,1)=k1+k3;k(1,3)=-k3;k(1,4)=

7、11*k3;k(2,2)=k2+k4;k(2,3)=-k4;k(2,4)=-12*k4;k(3,3)=k3+k4+k5+k6+k7+k8;k(3,4)=-11*k3+12*k4-13*k7-14*k5-15*k6-16*k8;k(3,5)=-k5-k6;k(3,6)=c*k5-d*k6;k(3,7)=-k7-k8;k(3,8)=a*k7-b*k8;k(4,4)=11A2*k3+12A2*k4+13A2*k7+14A2*k5+15A2*k6+16A2*k8;k(4,5)=14*k5+15*k6;k(4,6)=-c*14*k5+d*15*k6;k(4,7)=13*k7+16*k8;k(4,8)=-

8、a*13*k7+b*16*k8;k(5,5)=k5+k6;k(5,6)=-c*k5+d*k6;k(6,6)=cA2*k5+dA2*k6;k(7,7)=k7+k8+k9;k(7,8)=-a*k7+b*k8;k(7,9)=-k9;k(8,8)=aA2*k7+bA2*k8;k(9,9)=k9;for i=1:9j=i：9；k(j,i)=k(i,j)endK=k./m;C=zeros(9,9);C(1,1)=c3;C(1,3)=-c3;C(1,4)=l1*c3;C(2,2)=c4;C(2,3)=-c4;C(2,4)=-l2*c4;C(3,3)=c3+c4+c5+c6+c7+c8;C(3,4)=-l1*

9、c3+l2*c4-l3*c7-l4*c5-l5*c6-l6*c8;C(3,5)=-c5-c6;C(3,6)=c*c5-d*c6;C(3,7)=-c7-c8;C(3,8)=a*c7-b*c8;C(4,4)=l1A2*c3+l2A2*c4+l3A2*c7+l4A2*c5+l5A2*c6+l6A2*c8;C(4,5)=l4*c5+l5*c6;C(4,6)=-c*l4*c5+d*l5*c6;C(4,7)=l3*c7+l6*c8;C(4,8)=-a*l3*c7+b*l6*c8;C(5,5)=c5+c6;C(5,6)=-c*c5+d*c6;C(6,6)=cA2*c5+dA2*c6;C(7,7)=c7+c8

10、+c9;C(7,8)=-a*c7+b*c8;C(7,9)=-c9;C(8,8)=aA2*c7+bA2*c8;C(9,9)=c9;C(9,9)=c9;for i=1:9j=i：9;C(j,i)=C(i,j)endC=C./m;kt=zeros(18,1);kt(10,1)=k1/m1;kt(11,1)=k2/m2;Q=zeros(18,size(pf,2); %产生(16*size)全零矩阵，其中 size 为二维矩阵 qf (1*9953)第二 维的尺寸(即列数)，Q为(16*9553)的矩阵Q(10,:)=pf;%将qf的值赋给第9行，：指9953Q(11,:)=pr; %将qr的值赋给第1

11、0行，：指9953p=zeros(9);%p为8*8的矩阵，其中每个数都是0q=eye(9); %q为8*8的单位矩阵，其中对角线上都是1W=p q;-K -C; %p q k c均为9*9的矩阵，w为18*18的矩阵nt=size(pf,2);%数据个数，即 9953nf=500;h=1/100;%步长，时间间隔x(:,1)=zeros(18,1);%x为18*1的矩阵，且每个数都是0t(1)=0;%定义初始时间%改进的欧拉算法for i=2:nt%nt=9953n=i-1;t(i)=h*i;xL=x(:,n)+h*(kt.*Q(:,n)+W*x(:,n);%当 n=1 时，x(：，n)表示

12、x 矩阵(18*9553)第n列的列矩阵，kt为16*16的矩阵，Q(；，n)表示矩阵(16*9553)的第n列的列矩阵；x(:,i)=x(:,n)+h/2*(kt.*Q(:,n)+W*x(:,n)+(kt.*Q(:,i)+W*xL);endxac9=W(18,:)*x;figure(3);plot(t,xac9,k),grid;axis auto;xlabel(时间(S);ylabel(垂直加速度rad/sA2);title(9DOF人体座椅垂直加速度时域变化);fs=100;%分析频率=1/0.01，0.01为时间间隔nf=500;%每5s 一段，单个数据样本个数nb=nf/2;%自谱验算

13、时的参数 重复点数n_n=19;%功率谱密度分析一一人体座椅俯系统Py09,f=periodogram(xac9(1:nf),hann(nf),nf,fs,onesided);Pya9=Py09;for i=2:(2*n_n-1)Py09,f=periodogram(xac9(i-1)*nb+1:(i-1)*nb+nf),hann(nf),nf,fs,onesided);Pya9=Pya9+Py09;endPya9=Pya9/(2*n_n-1);figure(4);subplot(2,1,1);loglog(f(1:nf/2),Pya9(1:nf/2),k,Linewidth,2);grid

14、on;set(gca,XMinorTick,off7YMinorTick,off7XMinorGrid,off7YMinorGrid,off)xlim(0.1 50)xlabel(频率Hz);ylabel (加速度功率谱密度(mA2/sA3);title(9DOF Model人体座椅垂直振动);hold on;delta_f=0.2;ai9=zeros(21,1);f_z9=zeros(21,1);% 中心频率for i=1:21f_l=0.5*0.89*1.26A(i-1);f_u=f_l*1.26;f_z9(i,1)=f_l*2A(1/6);% 求中心频率%若频率大于50HZ,循环结束%在

15、一个三分之一倍频程内积分for j=1:251f0=f(j);if(f0=f_l)&(f0=f_u)%频率值距1/3倍频带下限小于delta_f，但距上限大于delta_f if(abs(f0-f_l)delta_f) ai=Pya9(j)*abs(f0-f_l);%频率值距1/3倍频带上限小于delta_f，但距下限大于delta_f elseif(abs(f_u-f0)delta_f) ai=Pya9 (j) *(abs(f0-f_u)+delta_f);%频率值距1/3倍频带上限、下限均小于delta_felseif(abs(f0-f_l)delta_f)&(abs(f_u-f0)delta_f) ai=Pya9 (j) *(f_u-f_l);%频率值距1/3倍频带上限、下限均大于delta_felseai=Pya9 (j) *delta_f;end else ai=0;endai9(i,1)=ai9(i,1)+ai;endai9(i,1)=sqrt(ai9(i,1);end subplot(2,1,2);loglog(f_z9,ai9,k),grid;xlim(0.1 50)xlabel(