求解非线性方程-(cosh((4*((2*((361621451347209955078125*pi*om

华电机械

本人求受轴向力的悬臂梁的固有频率遇到个问题，用matlab编制了程序可是结果出不来，，，求高手指教；
clc
clear
syms E I a1 a2 m l s omiga  tr1 tr2  trr trr1 ;

%实验等截面梁参数
l=0.8;
d=0.01;
A=pi*d*d/4;
%梁材料参数（忽略悬臂梁的质量）
u=0.3;
rou=7800;
E=2e11;
I=pi*(d^4)/64;
m=rou*A;
s=1e4;
a1=sqrt((sqrt(s^2+4*m*omiga^2*E*I)-s)/(2*E*I));
a2=sqrt((sqrt(s^2+4*m*omiga^2*E*I)+s)/(2*E*I));

%计算总传递矩阵
tr1=[sin(a1*l) cos(a2*l) sinh(a2*l) cosh(a2*l);
    a1*cos(a1*l) -a1*sin(a1*l) a2*cosh(a2*l) a2*sinh(a2*l);
    -E*I*a1^2*sin(a1*l) -E*I*a1^2*cos(a1*l) E*I*a2^2*sinh(a2*l) E*I*a2^2*cosh(a2*l);
    (-E*I*a1^3)*cos(a1*l)  (E*I*a1^3)*sin(a1*l) (E*I*a2^3)*cosh(a2*l)  (E*I*a2^3)*sinh(a2*l)];

tr2 =[ -(a1*a2)/(a1^2 + a2^2),  1/a1,  -a2/(E*I*(a1^3 + a1*a2^2)),  0;
    a2^2/(a1^2 + a2^2),   0,    -1/(E*I*(a1^2 + a2^2)), 0;
-a1^4/(a2^2*(a1^2 + a2^2)), a1^2/a2^3, -a1^2/(E*I*a2^2*(a1^2 + a2^2)), 1/(E*I*a2^3);
%总传递矩阵
      tr=tr1.*tr2;
% 由悬臂梁的边界条件得特征方程
             trr=[ -(a1^2*sinh(a2*l))/(a1^2 + a2^2), cosh(a2*l)/a2;
                 -(cosh(a2*l)*(- E*I*a2^3 + s*a2))/(E*I*(a1^2 + a2^2)), 0 ];
             trr1=det(trr)  

得到的结果如下：
trr1 =

-(cosh((4*((2*((361621451347209955078125*pi*omiga^2)/4722366482869645213696 + 100000000)^(1/2))/125 + 160)^(1/2))/(5*pi^(1/2)))^2*(((361621451347209955078125*pi*omiga^2)/4722366482869645213696 + 100000000)^(1/2) - 10000))/(2*((361621451347209955078125*pi*omiga^2)/4722366482869645213696 + 100000000)^(1/2))


现在不知道怎么求trr1=0这个方程的解omiga。。。
试了用solve命令和fsolve命令好像都不行。。。。。。。哪位高手能指点迷津啊。。。实在是没办法啊。。

一路向前

trr1=0，omiga解为     0， -255.56703，  255.56703对不？对的话留言再讨论

华电机械

以下是悬臂梁不受轴向力的结果：
个人感觉轴向力对悬臂梁的固有频率不会影响很大，，而你算的是omiga=255.56703，换算成f=40.69，与第一阶不受轴向力的11.0688相差挺大的，，可能是我编写的程序有问题。。我再改进改进。。

华电机械

能发过求上述方程的源代码过来吗？？？感激。。。

华电机械

我又改进了下算的结果是这个：看看这个方程怎么解啊》》》
trr1 =

(390625*cosh((11*(((3779737860703086955078125*omiga^2)/12682136550675316736 + 100000000)^(1/2)/21856 + 625/1366)^(1/2))/25)^2*(625*(((3779737860703086955078125*omiga^2)/12682136550675316736 + 100000000)^(1/2)/21856 + 625/1366)^(1/2) - 683*(((3779737860703086955078125*omiga^2)/12682136550675316736 + 100000000)^(1/2)/21856 + 625/1366)^(3/2))*(((3779737860703086955078125*omiga^2)/12682136550675316736 + 100000000)^(1/2)/21856 + 625/1366)^(1/2)*(9676128923399902605*omiga^2 + 3246626956972881084416))/1034414265808893147420496494592

>>

一路向前

华电机械 发表于 2013-4-24 11:31

                               
登录/注册后可看大图

我又改进了下算的结果是这个：看看这个方程怎么解啊》》》
trr1 =

可能是你的数据或是程序有问题吧，一般轴向振动频率确实很小，比横向弯曲小很多，这样算出来是
528.23521135724722912687687058962
     -528.23521135724722912687687058962
  0.00035594562195853297287202288179877
-0.00035594562195853297287202288179877
   -18.317465070665901094120659558373*i
    18.317465070665901094120659558373*i，可能和结果差的更大，建议仔细查程序和理论公式，
，程序很简单，就是你的trr1精度太高了，matlab可能无法运算，你需要控制精度，可以用vpa(trr1,10),就把trr1截断到小数点后10位，即可运算。

令，记得点击回复，要不我看不到的

华电机械

一路向前 发表于 2013-4-24 12:34

                               
登录/注册后可看大图

可能是你的数据或是程序有问题吧，一般轴向振动频率确实很小，比横向弯曲小很多，这样算出来是
528.235 ...

好的，，谢谢你了，，，我仔细检查检查，感觉程序是有问题啊。

华电机械

一路向前 发表于 2013-4-24 12:34

                               
登录/注册后可看大图

可能是你的数据或是程序有问题吧，一般轴向振动频率确实很小，比横向弯曲小很多，这样算出来是
528.235 ...

控制精度又算了下得到这样的结果：
-(0.5*cosh(0.44*(4.58e-5*(2.98e5*omiga^2 + 1.0e8)^(1/2) + 0.458)^(1/2))^2*((2.98e5*omiga^2 + 1.0e8)^(1/2) - 1.0e4))/(2.98e5*omiga^2 + 1.0e8)^(1/2)
这个怎么解啊，刚才那么大可能是传递矩阵有问题

华电机械

下图方程的四个根（参考文献【利用Euler梁动模频型的计传算汽轮机 叶片静频和动频的传递矩形法】）是不是有错，我算的好像与它的不同，具体如word文档。。。请赐教。

一路向前

华电机械 发表于 2013-4-24 15:02

                               
登录/注册后可看大图

控制精度又算了下得到这样的结果：
-(0.5*cosh(0.44*(4.58e-5*(2.98e5*omiga^2 + 1.0e8)^(1/2) + 0.458) ...

这个怎么解？
syms omiga
solve(vpa(trr1,5))

悬臂梁的公式应当很多文献能查阅到，你自己看吧，我没做过

华电机械

恩，好吧，还是谢谢你