|
楼主 |
发表于 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 |
|