任意边界条件下矩形薄板的固有频率计算

xiaofacai

如题，小弟想计算一个对边简支一边固支一边自由的矩形薄板的振动频率，通过ANSYS已得到，由于还要计算声辐射声功率，想通过双向梁函数用MATLAB来计算，但结果与ANSYS计算结果差距很大，

看了很多天也找不出问题，请求各位 路过看一下

xiaofacai

rhos=7800;
hs=0.05;
a=2.5;
b=2.5;
Es=21e10;
nus=0.3;
Nmod=6;

Is=hs^3/12;
Ds=Es*Is/(1-nus^2);
MN=Nmod;
f='cosh(x)*cos(x)+1';
for j=1:MN
a1(j)=fzero(f,(j-1/2)*pi); %方程f中找(j+1/2)*pi附近的零点赋值到a1(j)
end
b1(j)=(sinh(a1)+sin(a1))/(cosh(a1)+cos(a1));%求Betan
syms au bu xx L;
for m=1:Nmod
de1=sin(m*pi*xx/L);
de2=(cosh(au*xx/L)-cos(au*xx/L))-bu*(sinh(au*xx/L)-sin(au*xx/L));
df1=diff(de1,xx,4);
df2=diff(de1,xx,2);
df3=diff(de2,xx,2);
df4=diff(de2,xx,4); %对de求4阶导数

dd1=int(df1*de1,xx,0,L);
dd2=int(de2*de2,xx,0,L);
dd3=int(df2*de1,xx,0,L);
dd4=int(df3*de2,xx,0,L);
dd5=int(df4*de2,xx,0,L);
dd6=int(de1*de1,xx,0,L);

     I1(m)=subs(dd1,{au,bu,L},{a1(m),b1(m),a}); %利用{a1(m)/a,b1(m),a}的值代替符号表达式中{au,bu,L}的值,重新计算了表达式 I1。
     I2(m)=subs(dd2,{au,bu,L},{a1(m),b1(m),b});
     I3(m)=subs(dd3,{au,bu,L},{a1(m),b1(m),a});
     I4(m)=subs(dd4,{au,bu,L},{a1(m),b1(m),b});
     I5(m)=subs(dd5,{au,bu,L},{a1(m),b1(m),b});
     I6(m)=subs(dd6,{au,bu,L},{a1(m),b1(m),a});
end
%计算固有频率
for pm=1:Nmod
     for qm=1:Nmod
      wj(pm,qm)=sqrt(Ds/rhos/hs)*sqrt((I1(pm)*I2(qm)+2*I3(pm)*I4(qm)+I5(qm)*I6(pm))/I2(qm)/I6(pm));
       wj=vpa(wj,4);
      f=wj/2/pi;
      f=vpa(f,4);
     end
end

xiaofacai

xiaofacai 发表于 2019-9-5 11:46
rhos=7800;
hs=0.05;
a=2.5;

附上所求的程序，不知道哪里不对

xiaofacai

[img][img][/img]

欧阳中华

.
    将两种方法计算结果的频率和振型对照比较也许能发现些什么 .. . .

xiaofacai

谢谢欧阳老师的回复！
频率的截图上传上来了，第一张是ANSYS和MATLAB有限元程序计算的结果，第二张是双向梁函数的计算结果，2*pi的转换也计算过了，不知道错在哪。