free12.m

汽车振动的频谱分析。很有意思

function [fz1,fx1,fz4,fz5,fz6,fz7,fz8,fy1,fy4,fy5,fy6,fy7,fy8,fx5]=free12(fm,c)
%mass matrix
mp=184.1;
mo=8.2;
jpx=6.77;
jpy=8.871;
jpz=6.584;
jox=0.62;
joy=0.2;
joz=1.27;
m_diag=[mp,mp,mp,mo,mo,mo,jpx,jpy,jpz,jox,joy,joz];
m=diag(m_diag);

%stiffness matrix
a=0.25;%variable
b=0.254;%variable
c1=0.334;
d=0.466;
l=0.2;
la=0.24;
lb=0.36;
lc=0.45;
ld=0.12;

F1=[1 0 0 0 0 0 0 -a 0 0 0 0
    0 1 0 0 0 0 a  0 b 0 0 0
    0 0 1 0 0 0 0 -b 0 0 0 0];
F2=[1 0 0 -1 0 0 0 -a -l 0 0 lb
    0 1 0 0 -1 0 a 0 c1 0 0 -la
    0 0 1 0 0 -1 l -c1 0 -lb la 0];
F3=[1 0 0 -1 0 0 0 -a -l 0 0 -lb
    0 1 0 0 -1 0 a 0 c1 0 0 -la
    0 0 1 0 0 -1 -l -c1 0 lb la 0];
F4=[1 0 0 0 0 0 0 -a 0 0 0 0
    0 1 0 0 0 0 a 0 -d 0 0 0
    0 0 1 0 0 0 0 d 0 0 0 0];
F5=[0 0 0 1 0 0 0 0 0 0 0 -lb
    0 0 0 0 1 0 0 0 0 0 0 ld
    0 0 0 0 0 1 0 0 0 lb -ld 0];
F6=[0 0 0 1 0 0 0 0 0 0 0 lb
    0 0 0 0 1 0 0 0 0 0 0 ld
    0 0 0 0 0 1 0 0 0 -lb -ld 0];
F7=[0 0 0 1 0 0 0 0 0 0 0 lb
    0 0 0 0 1 0 0 0 0 0 0 -lc
    0 0 0 0 0 1 0 0 0 -lb -lc 0];
F8=[0 0 0 1 0 0 0 0 0 0 0 -lb
    0 0 0 0 1 0 0 0 0 0 0 -lc
    0 0 0 0 0 1 0 0 0 lb -lc 0];
k11=79.11805e3;%variable
k12=79.11805e3;%variable
k13=79.11805e3;%variable
k1=[k11 0 0
    0 k12 0
    0 0 k13];
k22=1.137396e5;%variable
k23=1.137396e5;
k2=[0 0 0
    0 k22 0
    0 0 k23];
k32=1.358418e5;%variable
k33=1.358418e5;
k3=[0 0 0
    0 k32 0
    0 0 k33];
k41=79.11805e3;
k42=79.11805e3;
k43=79.11805e3;%variable
k4=[k41 0 0
    0 k42 0
    0 0 k43];
k51=273.6905e3;
k52=273.6905e3;
k53=273.6905e3;%variable
k5=[k51 0 0
    0 k52 0
    0 0 k53];
k61=273.6905e3;
k62=273.6905e3;
k63=273.6905e3;%variable
k6=[k61 0 0
    0 k62 0
    0 0 k63];
k71=259.8413e3;
k72=259.8413e3;
k73=259.8413e3;%variable
k7=[k71 0 0
    0 k72 0
    0 0 k73];
k81=259.8413e3;
k82=259.8413e3;
k83=259.8413e3;%variable
k8=[k81 0 0
   0 k82 0
   0 0 k83];
K=F1'*k1*F1+F2'*k2*F2+F3'*k3*F3+F4'*k4*F4+F5'*k5*F5...
  +F6'*k6*F6+F7'*k7*F7+F8'*k8*F8;
%damp matrix
C1=c*eye(3);
C2=[0 0 0
    0 c 0
    0 0 c];
C3=C2;
C4=C1;
C5=C1;
C6=C1;
C7=C1;
C8=C1;

C=F1'*C1*F1+F2'*C2*F2+F3'*C3*F3+F4'*C4*F4+F5'*C5*F5...
  +F6'*C6*F6+F7'*C7*F7+F8'*C8*F8;
y=0.05/0.135;
mp=0.33;
m1=0.3;
mj=mp+m1;
r=0.05;
delta=0.0865;
lm=1;

t=0.0025:0.0025:10;
f1=-0.04667*(2*pi*fm)^2*cos(4*pi*fm*t);
f2=-0.004*(2*pi*fm)^2*cos(4*pi*fm*t);
F1=fft(f1);
F2=fft(f2);

for n=1:4000
    f=0.1*n;
    sf3(n)=pi*(81e-3)^2*60*fm*lm/f/f;
    H=inv(K-(2*pi*f)^2*m+i*2*pi*f*C);
    sf1(n)=10*conj(F1(n))*F1(n);
    sf2(n)=10*conj(F2(n))*F2(n);  
%%%%%%%%%%%%%%%%%%%%%          Sf         %%%%%%%%%%%%%%%%%%%%%%
    i_Sf=[3,7,8];
    j_Sf=i_Sf;
    s_Sf=[sf1(n),sf3(n),sf2(n)];
    Sf=sparse(i_Sf,j_Sf,s_Sf,12,12); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    X=conj(H)*Sf*H.';
    x=sqrt(X);

    %compute the mount force(z)
    fz1(n)=x(3,3)-b*x(8,8);
    fz1(n)=sqrt(2)*sqrt(abs(k13*fz1(n))^2+abs(2*pi*f*c*fz1(n))^2);
    fz4(n)=x(3,3)+d*x(8,8);
    fz4(n)=sqrt(2)*sqrt(abs(k43*fz4(n))^2+abs(2*pi*f*c*fz4(n))^2);
    fz5(n)=x(6,6)+lb*x(10,10)-ld*x(11,11);
    fz5(n)=sqrt(2)*sqrt(abs(k53*fz5(n))^2+abs(2*pi*f*c*fz5(n))^2);
    fz6(n)=x(6,6)-lb*x(10,10)-ld*x(11,11);
    fz6(n)=sqrt(2)*sqrt(abs(k63*fz6(n))^2+abs(2*pi*f*c*fz6(n))^2);
    fz7(n)=x(6,6)-lb*x(10,10)+lc*x(11,11);
    fz7(n)=sqrt(2)*sqrt(abs(k73*fz7(n))^2+abs(2*pi*f*c*fz7(n))^2);
    fz8(n)=x(6,6)+lb*x(10,10)+lc*x(11,11);
    fz8(n)=sqrt(2)*sqrt(abs(k83*fz8(n))^2+abs(2*pi*f*c*fz8(n))^2);


    %compute the mount force(y)
    fy1(n)=x(2,2)+a*x(7,7)+b*x(9,9);
    fy1(n)=sqrt(2)*sqrt(abs(k12*fy1(n))^2+abs(2*pi*f*c*fy1(n))^2);
    fy4(n)=x(2,2)+a*x(7,7)-d*x(9,9);
    fy4(n)=sqrt(2)*sqrt(abs(k42*fy4(n))^2+abs(2*pi*f*c*fy4(n))^2);
    fy5(n)=x(5,5)+ld*x(12,12);
    fy5(n)=sqrt(2)*sqrt(abs(k52*fy5(n))^2+abs(2*pi*f*c*fy5(n))^2);
    fy6(n)=x(5,5)+ld*x(12,12);
    fy6(n)=sqrt(2)*sqrt(abs(k62*fy6(n))^2+abs(2*pi*f*c*fy6(n))^2);
    fy7(n)=x(5,5)-lc*x(12,12);
    fy7(n)=sqrt(2)*sqrt(abs(k72*fy7(n))^2+abs(2*pi*f*c*fy7(n))^2);
    fy8(n)=x(5,5)-lc*x(12,12);
    fy8(n)=sqrt(2)*sqrt(abs(k82*fy8(n))^2+abs(2*pi*f*c*fy8(n))^2);

    %compute the mount force(x)
    fx1(n)=x(1,1)-a*x(8,8);
    fx1(n)=sqrt(2)*sqrt(abs(k11*fx1(n))^2+abs(2*pi*f*c*fx1(n))^2);
    fx4(n)=x(1,1)-a*x(8,8);
    fx4(n)=sqrt(2)*sqrt(abs(k41*fx4(n))^2+abs(2*pi*f*c*fx4(n))^2);
    fx5(n)=x(4,4)-lb*x(12,12);
    fx5(n)=sqrt(2)*sqrt(abs(k51*fx5(n))^2+abs(2*pi*f*c*fx5(n))^2);

end