xiaofacai 发表于 2019-9-5 11:46

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

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

看了很多天也找不出问题,请求各位 路过看一下{:4_74:}

xiaofacai 发表于 2019-9-5 11:46

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 发表于 2019-9-5 11:47

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


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

xiaofacai 发表于 2019-9-5 11:59

欧阳中华 发表于 2019-9-5 13:56

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

xiaofacai 发表于 2019-9-5 14:20

谢谢欧阳老师的回复!
频率的截图上传上来了,第一张是ANSYS和MATLAB有限元程序计算的结果,第二张是双向梁函数的计算结果,2*pi的转换也计算过了,不知道错在哪。
页: [1]
查看完整版本: 任意边界条件下矩形薄板的固有频率计算