Rayleigh-Ritz 法求解功能梯度梁的固有频率
别人发的程序，该程序按照参考文献的步骤和公式，使用Rayleigh-Ritz 法结合经典梁理论求解功能梯度梁的固有频率，但结果和参考文献的不符。烦请大神帮忙看看。参考文献 见附件，程序如下：%%% Rayleigh-Ritz 法求梁的屈曲温度，弹性模量各不相同,基于CBT梁理论
%%% 基于CBT梁理论
%%% 单位均用国际单位制
clc
clear
format short
%--------------------------------------------------
%梁的几何尺寸 a b h分别为长宽高
%--------------------------------------------------
L=5; %%梁轴向长度
b=1; %%梁的宽度
h=1; %%梁的高度
T0=300;
deltaT=0;%温升
%---------------------------------------------------
%材料属性
%FGM由metal和Ceramic组成(top表示上面，bottom表示下面)
n=10; %%%%%%%%%%功能材料的梯度指数
% EL
Et=380E9;
%---------------------------------------------------
% ER;
Eb=70E9;
%---------------------------------------------------
Dent=3800;vt=0.3;%top 的材料属性
Denb=2700;vb=0.3;%bottom的材料属性
alpha_t=18.591E-6;
alpha_b=6.941E-6;
syms z
%---------------------------------------------------
%得到FGM的等效材料参数
%---------------------------------------------------
Ez=(Et-Eb)*(z/h+1/2)^n+Eb;
Denz=(Dent-Denb)*(z/h+1/2)^n+Denb;
vz=(vt-vb)*(z/h+1/2)^n+vb;
alphaz=(alpha_t-alpha_b)*(z/h+1/2)^n+alpha_b;
%----------------------------------------------------
%定义积分常数
%--------------------------------------------------
Q11z=Ez/(1-vz^2);Q55z=Ez/(2*(1+vz));
A00=int(Ez,z,-h/2,h/2);%无量纲化时要用
A11=int(Q11z,z,-h/2,h/2);
B11=int(z*Q11z,z,-h/2,h/2);
D11=int(z^2*Q11z,z,-h/2,h/2);
I0=int(Denz,z,-h/2,h/2);
I1=int(Denz*z,z,-h/2,h/2);
I2=int(Denz*z^2,z,-h/2,h/2);
%小数化数据
A00=double(A00);
A11=double(A11);
B11=double(B11);
D11=double(D11);
I0=double(I0);
I1=double(I1);
I2=double(I2);
%---------------------------------------------------
%TSTB的位移和转角等关系，径向u0=a(t)*u(x);转角v0=b(t)*theta(x);横向w0=c(t)*phi(x)
%---------------------------------------------------
ul=1;ur=1;
pl=1;pr=1; %边界条件的选择
%%---------------------------------------------------
syms x z
m=10; %可求得的前m阶梁的频率值（试函数的项数）
i=1:m;
f_s(i)=(x/L).^(i-1);
u_x=(x/L)^ul*(1-x/L)^ur*f_s';
phi_x=(x/L)^pl*(1-x/L)^pr*f_s';
%以上均为矩阵表示
%------------------------------------------------------
%计算刚度矩阵和质量矩阵
%------------------------------------------------------
du1x=diff(u_x,x);
duphi1x=diff(phi_x,x);
duphi2x=diff(duphi1x,x);
K1=A11*int(du1x*du1x',x,0,L);
K2=B11*int(du1x*duphi2x',x,0,L);
K3=D11*int(duphi2x*duphi2x',x,0,L);
%-------------------------------------
M1=I0*int(u_x*u_x',x,0,L);
M2=I0*int(phi_x*phi_x',x,0,L);
M3=I1*int(u_x*diff(phi_x,x)',x,0,L);
M4=I2*int(diff(phi_x)*diff(phi_x'),x,0,L);
%-----------------------------------------------------
KE=;
ME=;
KE=double(KE);
ME=double(ME);
=eig(KE,ME);%求特征值
n1=1;
for kk=1:m
D(m+kk,m+kk);
p=sqrt(D(m+kk,m+kk))*L*L/h*sqrt(Denb/Eb);
pinlv(kk)=p; %无量纲频率
end
format short
pinlv'
ff1=V(m+1:2*m,1)'*phi_x;%第n1阶模态函数 横向振动
x=0:0.01:L;
f=eval(ff1);
figure(1)
fmax=max(abs(f));%最大位移值
plot(x,f./fmax,'-.r','Linewidth',3);
grid on
hold on
页:
[1]