声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 7622|回复: 22

[共享资源] MATLAB编程求单跨梁振型及频率(两端任意约束)

[复制链接]
发表于 2006-4-12 16:24 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
  1. b=input('please input b.' );
  2. h=input('please input h.' );
  3. Rho=input('please input Rho.' );
  4. E=input('please input E.' );
  5. ds=input('How many elments do you want to devide?');
  6. l=input('Pleaes input the length');
  7. Lres=input('Pleaes input the left restraint,Fix->3,Hinged2->2,Hinged1->1,Free->0 ');
  8. Rres=input('Pleaes input the left restraint,Fix->3,Hinged2->2,Hinged1->1,Free->0 ');
  9. I=b*(h^3)/12;
  10. EI=E*I;

  11. x1=0;x2=sym('L');
  12. x=sym('x');
  13. j=0:3;v=x.^j;
  14. a=[1,x1,x1^2,x1^3;
  15. 0,1,2*x1,3*x1^2;
  16. 1,x2,x2^2,x2^3;
  17. 0,1,2*x2,3*x2^2];
  18. d=v/a;
  19. dt=d';
  20. m=dt*d;
  21. M=Rho*b*h*int(m,x,0,'L');
  22. M=subs(M,'L',l/ds);
  23. M=double(M);
  24. Ni=diff(d,x,2);
  25. Nt=Ni';
  26. k=Nt*Ni;
  27. KL=EI*int(k,x,0,'L');
  28. KL=subs(KL,'L',l/ds);
  29. KL=double(KL);

  30. KGC=cell(1,ds);
  31. for n=1:ds
  32. KGC{n}=zeros((ds+1)*2);
  33. end
  34. for n=1:ds;
  35. for i=1:4
  36. for j=1:4
  37. KGC{n}(i+2*n-2,j+2*n-2)=KL(i,j);
  38. end
  39. end
  40. end
  41. KG=zeros((ds+1)*2);
  42. for n=1:ds
  43. KG=KGC{n}+KG;
  44. end

  45. MGC=cell(1,ds);
  46. for n=1:ds
  47. MGC{n}=zeros((ds+1)*2);
  48. end
  49. for n=1:ds;
  50. for i=1:4
  51. for j=1:4
  52. MGC{n}(i+2*n-2,j+2*n-2)=M(i,j);
  53. end
  54. end
  55. end
  56. MG=zeros((ds+1)*2);
  57. for n=1:ds
  58. MG=MGC{n}+MG;
  59. end

  60. if Lres==3;
  61. MG=MG(3:end,3:end);KG=KG(3:end,3:end);
  62. for n=1:ds;
  63. KGC{n}=KGC{n}(3:end,3:end);
  64. end
  65. elseif Lres==2|Lres==1;
  66. MG=MG(2:end,2:end);KG=KG(2:end,2:end);
  67. for n=1:ds;
  68. KGC{n}=KGC{n}(2:end,2:end);
  69. end
  70. else MG=MG;KG=KG;
  71. end

  72. if Rres==3;
  73. MG=MG(1:end-2,1:end-2);KG=KG(1:end-2,1:end-2);
  74. for n=1:ds;
  75. KGC{n}=KGC{n}(1:end-2,1:end-2);
  76. end
  77. elseif Rres==2|Rres==1;
  78. MG=MG([1:end-2,end],[1:end-2,end]);KG=KG([1:end-2,end],[1:end-2,end]);
  79. for n=1:ds;
  80. KGC{n}=KGC{n}([1:end-2,end],[1:end-2,end]);
  81. end
  82. else MG=MG;KG=KG;
  83. end

  84. [EVector,EValue]=eig(KG/MG);
复制代码

[ 本帖最后由 suffer 于 2006-10-9 20:00 编辑 ]

点评

赞成: 3.0
赞成: 3
  发表于 2014-3-27 18:50
回复
分享到:

使用道具 举报

发表于 2006-4-30 23:13 | 显示全部楼层
<P>验算了一下,发现频率值偏大很多!振型也需要校核下!<BR>大家的看法呢?</P>
发表于 2006-5-1 01:26 | 显示全部楼层
我算了简支梁:b=0.05m,h=0.008m,l=0.68m,E=2.06e11,mm=0.0314kg/cm,这是我们实验做的对象,实测数据证实了理论值完全可靠。开始只算了ds=2,只有一阶频率符合地很好,后来算了ds=8,发现前三阶频率都符合地很好。理论值f=pi/2*n^2*sqrt(EI/mm/l^4)=40.19*n^2;取ds=8,从算得的EValue中取出最小的三值:63764,1.0207e+006,5.178e+006,f=sqrt(63764)/2/pi=40.19,...,等等,符合得很好。<BR>楼主真是强人!<BR>我开始只用前面部分算质量矩阵和刚度矩阵,结果误差很大。《有限元方法编程》(ISBN7-5053-9115-1)一书中说分布质量矩阵用集中质量矩阵近似表示会在预测梁的频率方面产生较大的误差,可能是这个原因吧。<BR>
 楼主| 发表于 2006-5-18 17:05 | 显示全部楼层
<P>谢谢各位学友的验证,这是我论文的一部分,但导师说他没有经费让我们做实验,同学们都学ANSYS了,本想让他们帮我算算的,但无赖他们只会静力分析,何况ANSYS算出后我也不知道怎样把ANSYS的数据导入MATLAB里,因为把振型和频率以及质量阵和刚度阵算出来后,用MATLAB编程诊断结构损伤。所以干脆我就用MATLAB算频率和振型了。第一次编程,还没有来得及验算理论值是不是对的就发上来了,谢谢学友的指正!</P>
发表于 2006-5-30 14:52 | 显示全部楼层
很好的程序,怎样提高运算精确度?<BR>
发表于 2006-5-30 14:59 | 显示全部楼层
有运行错无啊
发表于 2006-5-30 15:00 | 显示全部楼层
有点小错误,怎么调节?
发表于 2006-5-30 16:10 | 显示全部楼层
确实不错,比只会用ANSYS的强多了
 楼主| 发表于 2006-6-7 15:39 | 显示全部楼层

的确会有点错误。

<P>是的,就是那个input语句的存在而导致的。我用MATLAB6.5运行,有时可以通过,有时MATLAB却假死。非得强制关闭。input语句是想做得更加人性化一点。但无赖却有点小问题。希望学友们指正。想每次运行都通过,不用那个语句直接赋值就可以了。</P>
发表于 2006-8-22 15:30 | 显示全部楼层
强!请问有没有做过求薄板频率和振型阿?
发表于 2006-9-14 12:50 | 显示全部楼层
大哥,还有没有其他关于这方面的东西
发表于 2006-10-23 20:19 | 显示全部楼层
运行的时候出现   Warning: Explicit integral could not be found.    这是怎么回事呢?
发表于 2006-10-24 15:51 | 显示全部楼层
原帖由 haiyangnvshen 于 2006-10-23 20:19 发表
运行的时候出现   Warning: Explicit integral could not be found.    这是怎么回事呢?


估计你的数据是随便输入的
导致积分不存在
发表于 2006-10-24 15:52 | 显示全部楼层
原帖由 xinyuxf 于 2006-8-22 15:30 发表
强!请问有没有做过求薄板频率和振型阿?


呵呵,你去写一个发出来就有了
发表于 2006-10-24 15:53 | 显示全部楼层
原帖由 wsh6759 于 2006-5-30 14:52 发表
很好的程序,怎样提高运算精确度?<BR>


改变建模方法,把集总质量法改为有限元法估计精度就能提高很多
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-12 13:37 , Processed in 0.071590 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表