按照书上给的矩阵编程，求解出来的固有频率不对，求大佬帮忙看看






clear;clc;
%初始参数：质量，转矩，刚度
k_s1=5e8;k_s2=5e8;k_s3=5e8;k_r1=5e8;k_r2=5e8;k_r3=5e8;
k_c=1e8;k_r=1e8;k_s=1e8;k_p=1e8;k_ct=0;k_st=0;k_pt=0;k_rt=1e9;
f_1=0;f_2=2*pi/3;f_3=4*pi/3;m_c=5.43;m_r=2.35;m_s=0.4;m_p=0.66;
a_s=24.6*pi/180;a_r=a_s;
f_s1=f_1-a_s;f_s2=f_2-a_s;f_s3=f_3-a_s;
f_r1=f_1+a_r;f_r2=f_2+a_r;f_r3=f_3+a_r;
%1质量矩阵
M_c=diag([m_c,m_c,6.29]);
M_r=diag([m_r,m_r,3]);
M_s=diag([m_s,m_s,0.39]);
M_p=diag([m_p,m_p,0.61]);
M=blkdiag(M_c,M_r,M_s,M_p,M_p,M_p);%质量矩阵M
%2螺旋矩阵
G=zeros(18,18);
G(1,2)=-2*m_c;
G(4,5)=-2*m_r;
G(7,8)=-2*m_s;
G(10,11)=-2*m_p;
G(13,14)=-2*m_p;
G(16,17)=-2*m_p;
G;%质量矩阵M
k_bc=diag([k_c k_c k_ct]);
k_br=diag([k_r k_r k_rt]);
k_bs=diag([k_s k_s k_st]);
k_c1_1=k_p*[1 0 -sin(f_1);0 1 cos(f_1);-sin(f_1) cos(f_1) 1 ];
k_c1_2=k_p*[1 0 -sin(f_2);0 1 cos(f_2);-sin(f_2) cos(f_2) 1 ];
k_c1_3=k_p*[1 0 -sin(f_3);0 1 cos(f_3);-sin(f_3) cos(f_3) 1 ];
k_c2_1=k_p*[-1 0 0;0 -1 0;sin(f_1) -cos(f_1) 0];
k_c2_2=k_p*[-1 0 0;0 -1 0;sin(f_2) -cos(f_2) 0];
k_c2_3=k_p*[-1 0 0;0 -1 0;sin(f_3) -cos(f_3) 0];
k_c3_1=k_p*[1 0 0;0 1 0;0 0 0];k_c3_2=k_c3_1;k_c3_3=k_c3_1;
K_b=[k_bc+k_c1_1+k_c1_2+k_c1_3,zeros(3,3),zeros(3,3),k_c2_1,k_c2_2,k_c2_3;
    zeros(3,3),k_br,zeros(3,3),zeros(3,3),zeros(3,3),zeros(3,3);
    zeros(3,3),zeros(3,3),k_bs,zeros(3,3),zeros(3,3),zeros(3,3);
    k_c2_1,zeros(3,3),zeros(3,3),k_c3_1,zeros(3,3),zeros(3,3);
    k_c2_2,zeros(3,3),zeros(3,3),zeros(3,3),k_c3_2,zeros(3,3);
    k_c2_3,zeros(3,3),zeros(3,3),zeros(3,3),zeros(3,3),k_c3_3];%支撑刚度矩阵Kb
%向心刚度矩阵
K_omg=diag([m_c m_c 0 m_r m_r 0 m_s m_s 0 m_p m_p 0 m_p m_p 0 m_p m_p 0]);
%6啮合刚度矩阵
K_s1_1=k_s1*([sin(f_s1)^2, -cos(f_s1)*sin(f_s1),-sin(f_s1);0,cos(f_s1)^2,cos(f_s1);0,0,1]+...
    +triu([sin(f_s1)^2, -cos(f_s1)*sin(f_s1),-sin(f_s1);0,cos(f_s1)^2,cos(f_s1);0,0,1],1)');
K_s1_2=k_s2*([sin(f_s2)^2, -cos(f_s2)*sin(f_s2),-sin(f_s2);0,cos(f_s2)^2,cos(f_s2);0,0,1]+...
    +triu([sin(f_s2)^2, -cos(f_s2)*sin(f_s2),-sin(f_s2);0,cos(f_s2)^2,cos(f_s2);0,0,1],1)');
K_s1_3=k_s3*([sin(f_s3)^2, -cos(f_s3)*sin(f_s3),-sin(f_s3);0,cos(f_s3)^2,cos(f_s3);0,0,1]+...
    +triu([sin(f_s3)^2, -cos(f_s3)*sin(f_s3),-sin(f_s3);0,cos(f_s3)^2,cos(f_s3);0,0,1],1)');
K_s2_1=k_s1*[-sin(f_s1)^2 cos(f_s1)*sin(f_s1) -sin(f_s1);cos(f_s1)*sin(f_s1) -cos(f_s1)^2 cos(f_s1);...
    sin(f_s1) -cos(f_s1),1];
K_s2_2=k_s2*[-sin(f_s2)^2 cos(f_s2)*sin(f_s2) -sin(f_s2);cos(f_s2)*sin(f_s2) -cos(f_s2)^2 cos(f_s2);...
    sin(f_s2) -cos(f_s2),1];
K_s2_3=k_s3*[-sin(f_s3)^2 cos(f_s3)*sin(f_s3) -sin(f_s3);cos(f_s3)*sin(f_s3) -cos(f_s3)^2 cos(f_s3);...
    sin(f_s3) -cos(f_s3),1];
K_s3_1=k_s1*([sin(f_s1)^2,-cos(f_s1)*sin(f_s1),sin(f_s1);0,cos(f_s1)^2 -cos(f_s1);0 0 1]+...
+triu([sin(f_s1)^2,-cos(f_s1)*sin(f_s1),sin(f_s1);0,cos(f_s1)^2 -cos(f_s1);0 0 1],1)');
K_s3_2=k_s2*([sin(f_s2)^2,-cos(f_s2)*sin(f_s2),sin(f_s2);0,cos(f_s2)^2 -cos(f_s2);0 0 1]+...
    +triu([sin(f_s2)^2,-cos(f_s2)*sin(f_s2),sin(f_s2);0,cos(f_s2)^2 -cos(f_s2);0 0 1],1)');
K_s3_3=k_s3*([sin(f_s3)^2,-cos(f_s3)*sin(f_s3),sin(f_s3);0,cos(f_s3)^2 -cos(f_s3);0 0 1]+...
    +triu([sin(f_s3)^2,-cos(f_s3)*sin(f_s3),sin(f_s3);0,cos(f_s3)^2 -cos(f_s3);0 0 1],1)');
K_r1_1=k_r1*([sin(f_r1)^2,-sin(f_r1)*cos(f_r1),-sin(f_r1);0,cos(f_r1)^2,cos(f_r1);0,0,1]+...
    +triu([sin(f_r1)^2,-sin(f_r1)*cos(f_r1),-sin(f_r1);0,cos(f_r1)^2,cos(f_r1);0,0,1],1)');
K_r1_2=k_r2*([sin(f_r2)^2,-sin(f_r2)*cos(f_r2),-sin(f_r2);0,cos(f_r2)^2,cos(f_r2);0,0,1]+...
    +triu([sin(f_r2)^2,-sin(f_r2)*cos(f_r2),-sin(f_r2);0,cos(f_r2)^2,cos(f_r2);0,0,1],1)');
K_r1_3=k_r3*([sin(f_r3)^2,-sin(f_r3)*cos(f_r3),-sin(f_r3);0,cos(f_r3)^2,cos(f_r3);0,0,1]+...
    +triu([sin(f_r3)^2,-sin(f_r3)*cos(f_r3),-sin(f_r3);0,cos(f_r3)^2,cos(f_r3);0,0,1],1)');
K_r2_1=k_r1*([-sin(f_r1)^2,cos(f_r1)*sin(f_r1),sin(f_r1);0,-cos(f_r1)^2,-cos(f_r1);0 0 -1]+...
    triu([-sin(f_r1)^2,cos(f_r1)*sin(f_r1),sin(f_r1);0,-cos(f_r1)^2,-cos(f_r1);0 0 -1],1)');
K_r2_2=k_r2*([-sin(f_r2)^2,cos(f_r2)*sin(f_r2),sin(f_r2);0,-cos(f_r2)^2,-cos(f_r2);0 0 -1]+...
    triu([-sin(f_r2)^2,cos(f_r2)*sin(f_r2),sin(f_r2);0,-cos(f_r2)^2,-cos(f_r2);0 0 -1],1)');
K_r2_3=k_r3*([-sin(f_r3)^2,cos(f_r3)*sin(f_r3),sin(f_r3);0,-cos(f_r3)^2,-cos(f_r3);0 0 -1]+...
    triu([-sin(f_r3)^2,cos(f_r3)*sin(f_r3),sin(f_r3);0,-cos(f_r3)^2,-cos(f_r3);0 0 -1],1)');
K_r3_1=k_r1*([sin(f_r1)^2,-cos(f_r1)*sin(f_r1),-sin(f_r1);0,cos(f_r1)^2,cos(f_r1);0,0,1]+...
    triu([sin(f_r1)^2,-cos(f_r1)*sin(f_r1),-sin(f_r1);0,cos(f_r1)^2,cos(f_r1);0,0,1],1)');
K_r3_2=k_r2*([sin(f_r2)^2,-cos(f_r2)*sin(f_r2),-sin(f_r2);0,cos(f_r2)^2,cos(f_r2);0,0,1]+...
    triu([sin(f_r2)^2,-cos(f_r2)*sin(f_r2),-sin(f_r2);0,cos(f_r2)^2,cos(f_r2);0,0,1],1)');
K_r3_3=k_r3*([sin(f_r3)^2,-cos(f_r3)*sin(f_r3),-sin(f_r3);0,cos(f_r3)^2,cos(f_r3);0,0,1]+...
    triu([sin(f_r3)^2,-cos(f_r3)*sin(f_r3),-sin(f_r3);0,cos(f_r3)^2,cos(f_r3);0,0,1],1)');
K_pp_1=K_r3_1+K_s3_1;K_pp_2=K_r3_2+K_s3_2;K_pp_3=K_r3_3+K_s3_3;
skm_K_r1_n=K_r1_1+K_r1_2+K_r1_3;
skm_K_s1_n=K_s1_1+K_s1_2+K_s1_3;
K_m=[zeros(3,18);
     zeros(3,3),skm_K_r1_n,zeros(3,3),K_r2_1,K_r2_2,K_r2_3;
     zeros(3,3),zeros(3,3),skm_K_s1_n,K_s2_1,K_s2_2,K_s2_3;
     zeros(3,3),K_r2_1,K_s2_1,K_pp_1,zeros(3,3),zeros(3,3);
     zeros(3,3),K_r2_2,K_s2_2,zeros(3,3),K_pp_2,zeros(3,3);
     zeros(3,3),K_r2_3,K_s2_3,zeros(3,3),zeros(3,3),K_pp_3];
K=K_m+K_b
A=inv(M)*K;  %% 函数inv()是求矩阵的逆的函数；

[V,D]=eig(A);  %% 特征值和特征向量；

%% 计算固有频率并按大小排序

ww=diag(D);  %% 提出特征值

w1=sqrt(ww);  %%  计算出固有频率

w=sort(w1/2/pi);  %%  固有频率排序


1.jpg



2.jpg



3.jpg



4.jpg 返回小木虫查看更多