Rayleigh-Ritz 法求解功能梯度梁的固有频率

Apologize

别人发的程序，该程序按照参考文献的步骤和公式，使用Rayleigh-Ritz 法结合经典梁理论求解功能梯度梁的固有频率，但结果和参考文献的不符。烦请大神帮忙看看。参考文献 见附件，程序如下：

1. %%% Rayleigh-Ritz 法求梁的屈曲温度，弹性模量各不相同,基于CBT梁理论
   
2. %%% 基于CBT梁理论
   
3. %%% 单位均用国际单位制
   
4. clc
   
5. clear
   
6. format short
   
7. %--------------------------------------------------
   
8. %  梁的几何尺寸 a b h分别为长宽高
   
9. %--------------------------------------------------
   
10. L=5; %%梁轴向长度
    
11. b=1; %%梁的宽度
    
12. h=1; %%梁的高度
    
13. T0=300;
    
14. deltaT=0;%温升
    
15. %---------------------------------------------------
    
16. %材料属性
    
17. %FGM由metal和Ceramic组成(top表示上面，bottom表示下面)
    
18. n=10;                                    %%%%%%%%%%功能材料的梯度指数
    
19. % EL
    
20. Et=380E9;        
    
21. %---------------------------------------------------
    
22. % ER;
    
23. Eb=70E9;
    
24. %---------------------------------------------------
    
25. Dent=3800;vt=0.3;%top 的材料属性
    
26. Denb=2700;vb=0.3;%bottom的材料属性
    
27. alpha_t=18.591E-6;
    
28. alpha_b=6.941E-6;
    
29. syms z
    
30. %---------------------------------------------------
    
31. %得到FGM的等效材料参数
    
32. %---------------------------------------------------
    
33. Ez=(Et-Eb)*(z/h+1/2)^n+Eb;
    
34. Denz=(Dent-Denb)*(z/h+1/2)^n+Denb;
    
35. vz=(vt-vb)*(z/h+1/2)^n+vb;
    
36. alphaz=(alpha_t-alpha_b)*(z/h+1/2)^n+alpha_b;
    
37. %----------------------------------------------------
    
38. %定义积分常数
    
39. %--------------------------------------------------
    
40. Q11z=Ez/(1-vz^2);Q55z=Ez/(2*(1+vz));
    
41. A00=int(Ez,z,-h/2,h/2);%无量纲化时要用
    
42. A11=int(Q11z,z,-h/2,h/2);
    
43. B11=int(z*Q11z,z,-h/2,h/2);
    
44. D11=int(z^2*Q11z,z,-h/2,h/2);
    
45. I0=int(Denz,z,-h/2,h/2);
    
46. I1=int(Denz*z,z,-h/2,h/2);
    
47. I2=int(Denz*z^2,z,-h/2,h/2);
    
48. %小数化数据
    
49. A00=double(A00);
    
50. A11=double(A11);
    
51. B11=double(B11);
    
52. D11=double(D11);
    
53. I0=double(I0);
    
54. I1=double(I1);
    
55. I2=double(I2);
    
56. %---------------------------------------------------
    
57. %TSTB的位移和转角等关系，径向u0=a(t)*u(x);转角v0=b(t)*theta(x);横向w0=c(t)*phi(x)
    
58. %---------------------------------------------------
    
59. ul=1;ur=1;
    
60. pl=1;pr=1; %边界条件的选择
    
61. %%---------------------------------------------------
    
62. syms x z
    
63. m=10;                                 %可求得的前m阶梁的频率值（试函数的项数）
    
64. i=1:m;
    
65. f_s(i)=(x/L).^(i-1);
    
66. u_x=(x/L)^ul*(1-x/L)^ur*f_s';
    
67. phi_x=(x/L)^pl*(1-x/L)^pr*f_s';
    
68. %以上均为矩阵表示
    
69. %------------------------------------------------------
    
70. %计算刚度矩阵和质量矩阵
    
71. %------------------------------------------------------
    
72. du1x=diff(u_x,x);
    
73. duphi1x=diff(phi_x,x);
    
74. duphi2x=diff(duphi1x,x);
    
75. K1=A11*int(du1x*du1x',x,0,L);
    
76. K2=B11*int(du1x*duphi2x',x,0,L);
    
77. K3=D11*int(duphi2x*duphi2x',x,0,L);
    
78. %-------------------------------------
    
79. M1=I0*int(u_x*u_x',x,0,L);
    
80. M2=I0*int(phi_x*phi_x',x,0,L);
    
81. M3=I1*int(u_x*diff(phi_x,x)',x,0,L);
    
82. M4=I2*int(diff(phi_x)*diff(phi_x'),x,0,L);
    
83. %-----------------------------------------------------
    
84. KE=[K1+K1',-2*K2;-2*K2',K3+K3'];
    
85. ME=[M1+M1',-2*M3;-2*M3',M2+M2'+M4+M4'];
    
86. KE=double(KE);
    
87. ME=double(ME);
    
88. [V,D]=eig(KE,ME);%求特征值
    
89. n1=1;
    
90. for kk=1:m
    
91. D(m+kk,m+kk);
    
92. p=sqrt(D(m+kk,m+kk))*L*L/h*sqrt(Denb/Eb);
    
93. pinlv(kk)=p; %无量纲频率
    
94. end
    
95. format short
    
96. pinlv'
    
97. ff1=V(m+1:2*m,1)'*phi_x;%第n1阶模态函数 横向振动
    
98. x=0:0.01:L;
    
99. f=eval(ff1);
    
100. figure(1)
     
101. fmax=max(abs(f));%最大位移值
     
102. plot(x,f./fmax,'-.r','Linewidth',3);
     
103. grid on
     
104. hold on

复制代码