汽車平順性仿真分析程序

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;% 計(jì)算頻率變化量Gq_n0=256*1e-6;%選取B級路面不平度等級n0=0.1;%給參考空間頻率賦值u=40/3.6;%速度換算md=2*pi*rand(1001,1);%得出1001個隨機(jī)變量矩陣for t=0:0.01:100%計(jì)算前輪路面輸入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;%將時(shí)間變量賦值給變量ap(j)=q;%將輸入變量賦值給pj=j+1;q=0;endfigure(1);plot(a,p); %作圖grid on;%給圖形加上柵格線xlim(0,100);ylim(auto);xlabel(時(shí)間(s);% 標(biāo)注ylabel(路面不平度位移(m);title(路面時(shí)域隨機(jī)數(shù)據(jù)(三角級數(shù)法);%對模擬隨機(jī)輸入數(shù)據(jù)進(jìn)行功率譜分析，采用矩形 窗%fs分析頻率,nf一頻域的采樣點(diǎn)數(shù)Py0,f=periodogr

3、am(p(1:500),hann(500),500,100,onesided); %利用周期圖法求前 500 個點(diǎn)(即前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個點(diǎn)(即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、對數(shù)坐標(biāo)作圖grid on; xlim(0.1,50); xlabel(頻率(Hz);% 標(biāo)注ylabel (路面位移功率譜密度(mA2/Hz);hold on %將下面的圖形一起繪制在figure(2)中 for i=1:500%計(jì)算路面隨機(jī)輸入理論上的功率譜密度fb(i)=0.1*i; Sq(i)=Gq_n0*0.01*u/(fb(i)A2);endloglog(fb,Sq,-);%采用對數(shù)坐標(biāo)作圖%根據(jù)軸距l(xiāng)w和車速V(m/s)生成后輪輸入時(shí)域數(shù)據(jù) lw=6.53;% 給軸距賦值=td=lw/u;%計(jì)算前后輪時(shí)間差l=round(td/0.01);%對 round 里面的數(shù)取整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); %產(chǎn)生(16*size)全零矩陣，其中 size 為二維矩陣 qf (1*9953)第二 維的尺寸(即列數(shù))，Q為(16*9553)的矩陣Q(10,:)=pf;%將qf的值賦給第9行，：指9953Q(11,:)=pr; %將qr的值賦給第1

11、0行，：指9953p=zeros(9);%p為8*8的矩陣，其中每個數(shù)都是0q=eye(9); %q為8*8的單位矩陣，其中對角線上都是1W=p q;-K -C; %p q k c均為9*9的矩陣，w為18*18的矩陣nt=size(pf,2);%數(shù)據(jù)個數(shù)，即 9953nf=500;h=1/100;%步長，時(shí)間間隔x(:,1)=zeros(18,1);%x為18*1的矩陣，且每個數(shù)都是0t(1)=0;%定義初始時(shí)間%改進(jìn)的歐拉算法for i=2:nt%nt=9953n=i-1;t(i)=h*i;xL=x(:,n)+h*(kt.*Q(:,n)+W*x(:,n);%當(dāng) n=1 時(shí)，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(時(shí)間(S);ylabel(垂直加速度rad/sA2);title(9DOF人體座椅垂直加速度時(shí)域變化);fs=100;%分析頻率=1/0.01，0.01為時(shí)間間隔nf=500;%每5s 一段，單個數(shù)據(jù)樣本個數(shù)nb=nf/2;%自譜驗(yàn)算

13、時(shí)的參數(shù) 重復(fù)點(diǎn)數(shù)n_n=19;%功率譜密度分析一一人體座椅俯系統(tǒng)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,循環(huán)結(jié)束%在

15、一個三分之一倍頻程內(nèi)積分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(