结构动力学~求解悬臂梁的振动问题(振型及固有频率)~Matlab代码

利用有限元方法求解悬臂梁的振动问题,输入节点数可以算出前六阶振型以及计算固有频率

%% 利用有限元方法求解悬臂梁的振动问题 
clc
clear
syms x l rho b t E
A = b*t;
I = b*t^3/12;
n = input ('Please input the number of discrete elements n = '); % input the  total elements of the discrete element 
%% difine the shape function 
N1 = 1-3*(x/l)^2+2*(x/l)^3; 
N2 = (x/l-2*(x/l)^2+(x/l)^3)*l;
N3 = 3*(x/l)^2-2*(x/l)^3;
N4 = ((x/l)^3-(x/l)^2)*l;
%% solve the element mass matrix and stiffness matix, total mass matrix and stiffness matrix 
N = [N1,N2,N3,N4];
Me0 =  int(N'*N,x,0,l); % obtain the element mass matrix 
Me = rho*A*Me0;
Ke0 = int(diff(N.',x,2)*diff(N,x,2),x,0,l); % obtain the element stiffness matrix 
Ke = (E*I)*Ke0;
M0 = zeros(2*(n+1),2*(n+1));
K0 = zeros(2*(n+1),2*(n+1));
for i=1:1:n
    ae = zeros(4,2*(n+1)); 
    for j=1:1:4
        ae(j,2*i+j-2) = 1; % difine the coordinate transformation matrix
    end
    M0 = M0+ae.'* Me * ae; % solve the total mass stiffness matrix 
    K0 = K0+ae.'* Ke *ae; % solve the total stiffness matix 
end
disp('The element mass matrix Me = '); disp(Me); % show the element mass  matrix 
disp('The element stiffness matrix Ke =  ');disp(Ke); % show the element tiffness  matrix 
disp('The total mass matrix M = ');disp(M0); % show the total mass  matrix 
disp('The total stiffness matrix K = ');disp(K0); % show the total tiffness  matrix 
%%
Me1 = matlabFunction(Me);
Ke1 = matlabFunction(Ke);
M1 = matlabFunction(M0); % 将质量符号矩阵转换为代数矩阵
K1 = matlabFunction(K0); % 将刚度符号矩阵转换为代数矩阵
%%  input the geometry parameters and physical constants 
L = 0.4; % the total length of the beam, m
b = 0.02; % the width of the beam, m
t = 0.002; % the thickness of the beam, m 
E = 7*10^10; % the elastic modulus, GPa
rho = 2700; % the density of the beam, kg/m^3
l = L/n;
MNumeric = M1(b,l,t,rho);
KNumeric = K1(E,b,l,t);
disp('The numerical total mass matrix M = ');disp(MNumeric); % show the total mass  matrix 
disp('The numerical total stiffness matrix K = ');disp(KNumeric); % show the total tiffness  matrix
%% 1.悬臂梁的振动求解
%对于悬臂梁考虑约束条件节点1：挠度为零，转角为零；
%采用消元法：消去整体质量矩阵的前两行，前两列；消去整体刚度矩阵的前两行，前两列；
M = MNumeric([3:2*n+2],[3:2*n+2]);
K = KNumeric([3:2*n+2],[3:2*n+2]);
[X,lamda] = eig(K,M);
omega = sort(sqrt(diag(lamda)));
disp(omega);
XX = zeros(n,n);
XX(1,:) = 0;
for i = 2:n
    for j = 1:n
         XX(i,j) = X(2*(i-1)+1,j);
    end
end
%% 画出前六阶振型
r = 1:n;
for i = 1:n
    figure('Color',[1 1 1]);
    plot(r,XX(:,i));
    xlabel('\itx\rm/\itl','FontName','Times New Roman','FontSize',36);
    ylabel('u\rm(\itx\rm)','FontName','Times New Roman','FontSize',36);
    set(gca,'Linewidth',2); % setting up the linewidth of the coordinate
    set(gca,'FontName','Times New Roman','FontSize',36); 
    set(get(gca,'Children'),'linewidth',3); % setting up the linewidth of the curve
end
%% 根据振动力学方程得到悬臂梁固有频率
betal1 = 1.875; 
betal2 = 4.694;
betal = [];
omega_real = [];
error = [];
for i = 1:20
    betal(1) = betal1;
    betal(2) = betal2;
    if i>=3
        betal(i) = pi*(i-0.5);
    end
    omega_real(i) = (betal(i)/L)^2*sqrt(E*(b*t^3/12)/(rho*b*t));
    disp(omega_real(i));
end
if n>=20
    m = 20;
else 
    m = n;
end
for i = 1:m
      error (i) = (abs(omega_real(i)-omega(i))/omega_real(i))*100; % output the percentage error
      disp(error(i));
end