%此段程序可用来数值/理论计算、比较二维悬臂梁的固有频率和振型
clear %清空数据，否则有可能会因使用先前的数据而造成向量维数不一样画图错误
yyy=inline('cos(h)*cosh(h)+1','h'); %产生方程cos(h)*cosh(h)=-1 的数值解
i=1;
for a=1:0.01:100
b=a+0.01;
aa=(cos(a)*cosh(a)+1);
bb=(cos(b)*cosh(b)+1);
if aa*bb<=0
jie(i)=fzero(yyy,[a b]);
i=i+1;
end
end
jie; %jie为数值解的向量，此处为32个数值解
%-----------------------------
fprintf('此段程序可用来计算悬臂梁的固有频率和振型\n\n');
f=input('输入梁的单元数即将梁等分成几部分:'); %梁被分成f份
js=input('输入要查看的前js(不可大于32，也不可大于梁的单元数)阶振型:') ; %可查看1--js阶的图形
E=7e10; %梁的弹性模量
I=1e-2; %梁的截面惯性矩
l=1; %梁的长度
p=3e3; %梁的密度
A=1e-2; %梁的截面面积
L=l/f; %单位梁的长度
num=2*(f+1); %单位梁的总自由度
ke=E*I/(L^3)*[12 -6*L -12 -6*L; %梁单元的刚度矩阵，机械振动基础P137
-6*L 4*L^2 6*L 2*L^2;
-12 6*L 12 6*L;
-6*L 2*L^2 6*L 4*L^2];
me=p*A*L/420*[156 -22*L 54 13*L; %梁单元的一致质量矩阵，机械振动基础P138
-22*L 4*L^2 -13*L -3*L^2;
54 -13*L 156 22*L;
13*L -3*L^2 22*L 4*L^2;];
k=zeros(num); %初始化 总的刚度矩阵元素全为0
m=zeros(num); %初始化 总的质量矩阵元素全为0
for n=0:f-1;
for i=1:4;
for j=1:4;
k(i+2*n,j+2*n)=k(i+2*n,j+2*n)+ke(i,j); %单元刚度矩阵元素放入总的刚度矩阵里
m(i+2*n,j+2*n)=m(i+2*n,j+2*n)+me(i,j); %单元质量矩阵元素放入总的质量矩阵里
end
end
end
kk=k(3:num,3:num); %因为第一个节点为固定端，转角和挠度为0，故删去第一行第一列，第二行第二列
mm=m(3:num,3:num); %因为第一个节点为固定端，转角和挠度为0，故删去第一行第一列，第二行第二列
[v d]=eig(kk,mm); %求特征向量v 和特征值d
for j=1:2*f;
w(j)=sqrt(d(j,j)); % w（j）为第j个固有频率
v(:,j)=v(:,j)/v(2*f-1,j); %以每列最后一个数值为基准，进行归一化
for i=1:f
vv(i,j)=v(2*i-1,j); %将每列奇数位的值赋给新建的矩阵vv，该矩阵行数等于可动节点f*2
end
end
x=linspace(1/f,1,f); %建立一个x向量
xl=0:0.01:l; %建立一个xl向量
sL=jie; %进行赋值，为计算理论振型提供服务
s=sL./l;
wl=s.^2*sqrt(E*I/p/A); %理论固有频率
vn=-(sinh(sL)-sin(sL))./(cosh(sL)+cos(sL)); %为计算理论振型提供服务
for i=1:js
vl(:,i)=cosh(s(i)*xl)-cos(s(i)*xl)+vn(i)*(sinh(s(i)*xl)-sin(s(
i)*xl));
%理论固有振型函数，机械振动基础P105，原书中有错误，丢失了cosh（sn）后的x，已修正
vljs(:,i)=vl(:,i)/vl(101,i); %以每列最后一个数值为基准，进行归一化
end
for i=1:js
figure %显示前js阶振型
plot(x,vv(:,i),'b*') %数值计算所得振型
hold on
plot(xl,vljs(:,i),'r--') %理论计算所得振型
r=strcat('第',num2str(i),'阶固有振型');
title(r)
legend('数值计算值','理论计算值')
end
figure
plot(w,'b-') %下面进行固有频率比较
hold on
plot(wl,'rx')
title('数值计算和理论计算固有频率比较');
legend('数值计算所得固有频率','理论计算所得固有频率')