查看: 3107|回复: 21

[随机振动] 虚拟激励法程序????

  [复制链接]
发表于 2010-11-10 13:27 | 显示全部楼层 |阅读模式

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

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

x
那位兄台有虚拟激励法的程序????最好是从有限元里提取信息的?可否分享一下。谢了!

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2010-11-12 06:33 | 显示全部楼层
  1. clear;
  2. clc;
  3. n=4;
  4. %虚数单位
  5. II=sqrt(-1);
  6. %结构质量矩阵
  7. M=eye(n)*1.0e+4;
  8. K=eye(n)*1.6*1.0e+7;
  9. %结构刚度矩阵聚合
  10. zk=zeros(n);
  11. for j=1:(n-1)
  12. zk(j,j)=K(j,j)+K(j+1,j+1);
  13. zk(j,j+1)=-K(j+1,j+1);
  14. zk(j+1,j)=-K(j+1,j+1);
  15. end
  16. zk(n,n)=K(n,n);
  17. K=zk;
  18. [tzxl,tzz]=eig(zk,M);
  19. d=diag(sqrt(tzz));

  20. %阻尼矩阵
  21. %振型规一化
  22. for i=1:n
  23.    [tzz1(i),j]=min(d);
  24.    tzxl1(:,i)=tzxl(:,j);
  25.    d(j)=max(d)+1;
  26. end
  27. %振型归一化取第一层的振型为1
  28. for j=1:n
  29.    tzxl1(:,j)=tzxl1(:,j)/tzxl1(1,j);
  30. end

  31. w0=tzz1;
  32. w=tzz1/(2*pi);
  33. %振型zhx
  34. zhx=tzxl1;

  35. %阻尼比
  36. ks0=0.05;
  37. %各阶模态阻尼比均为0.05
  38. ks=ones(n,1)*ks0;

  39. %求广义质量
  40. Mn=zhx'*M*zhx;

  41. %求广义阻尼矩阵
  42. C=zeros(n);
  43. for i=1:n
  44.    C(i,i)=2*ks(i)*w0(i)*Mn(i,i);
  45. end
  46. %通过振型和广义阻尼矩阵求解阻尼矩阵
  47. c=(zhx')\C/zhx;

  48. %过滤白噪声参数
  49. eg=0.6;
  50. wg=15.708;
  51. S0=0.001574;
  52. %频率间隔
  53. dw=0.3;
  54. %最大频率范围
  55. maxw=45;
  56. %最大时间值
  57. maxt=40;
  58. %时间间隔
  59. dt=0.2;
  60. %各层各时间点频率点的功率谱密度,循环变量:层数,时间点,频率点
  61. Pwt=zeros(n,maxt/dt,maxw/dw);
  62. %频率点数循环变量wn
  63. wn=1;

  64. %对频率进行循环,求解各频率点的时间历程
  65. for w=0:dw:maxw
  66. x1=1+4*eg^2*(w/wg)^2;
  67. x2=(1-(w/wg)^2)^2+4*eg^2*(w/wg)^2;
  68. Sgw=x1*S0/x2;
  69. s=sqrt(Sgw);
  70. %采用精细积分法进行求解时间历程,得到位移和速度时程
  71. [disp,velp]=JINGXI67(M,zk,c,dt,maxt,w,s,n);
  72. Ywt=disp;
  73. for kkk=1:maxt/dt
  74. %求确定频率下各时间点的功率谱
  75. Yw=Ywt(:,kkk);
  76. y1=conj(Yw);
  77. y2=transpose(Yw);
  78. %确定时间点确定频率下的功率谱
  79. Syyw=y1*y2;
  80.    for kk=1:n
  81.    Pwt(kk,kkk,wn)=Syyw(kk,kk);
  82.    end
  83. end
  84. wn=wn+1;
  85. end

  86. %求解完成后实际上wn为maxw/dw+1
  87. %各层的时变方差,循环变量为:层数,时间点
  88. Fangcha=zeros(n,maxt/dt);
  89. for tn=1:maxt/dt
  90. %求解各层的时变方差
  91. for kk=1:n
  92. xx1=zeros(wn-1,1);
  93. %每一个时刻的方差对各频率点进行积分
  94.    for wn0=1:wn-1
  95.    xx1(wn0)=Pwt(kk,tn,wn0);
  96.    end
  97. %采用复合梯形求积公式对功率谱进行积分得到方差
  98. Fangcha(kk,tn)=(xx1(1)+xx1(wn-1)+2*sum(xx1(2:wn-1-1)))*dw;
  99. end
  100. end

  101. %画图
  102. c1=(1:maxt/dt)*dt;
  103. d1=Fangcha(1,:)/S0;
  104. d2=Fangcha(2,:)/S0;
  105. d3=Fangcha(3,:)/S0;
  106. d4=Fangcha(4,:)/S0;
  107. figure(3)
  108. plot(c1,d1,'k',c1,d2,'r',c1,d3,'m',c1,d4,'r-')
复制代码
  1. function [disp,velp]=JINGXI67(m,k,c,dt,maxt,w,s,n)
  2. II=sqrt(-1);
  3. IIW=II*w;
  4. I=eye(n);
  5. Z=zeros(n);
  6. A=[Z,I;-m\k,-m\c];
  7. AF=-1*[zeros(n,1);ones(n,1)];

  8. %采用科茨公式求解非齐次项
  9. %求解积分点的矩阵指数
  10. EXPDT=EXPHDT(A,dt);
  11. EXP0_75DT=EXPHDT(A,0.75*dt);
  12. EXP0_50DT=EXPHDT(A,0.50*dt);
  13. EXP0_25DT=EXPHDT(A,0.25*dt);

  14. zeroq=zeros(n,maxt/dt);
  15. disp=zeroq;  
  16. velp=zeroq;  

  17. zerop=zeros(n,1);
  18. dispt=zerop;
  19. velpt=zerop;
  20. UK=[dispt;velpt];

  21. for i=1:maxt/dt
  22. %给出各积分点的X坐标,即时间
  23. t=(i-1)*dt;
  24. t0_25=t+0.25*dt;
  25. t0_50=t+0.50*dt;
  26. t0_75=t+0.75*dt;
  27. t1=i*dt;
  28. %调制函数
  29. gt=4*(exp(-0.0995*t)-exp(-0.199*t));
  30. gt0_25=4*(exp(-0.0995*t0_25)-exp(-0.199*t0_25));
  31. gt0_50=4*(exp(-0.0995*t0_50)-exp(-0.199*t0_50));
  32. gt0_75=4*(exp(-0.0995*t0_75)-exp(-0.199*t0_75));
  33. gt1=4*(exp(-0.0995*t1)-exp(-0.199*t1));

  34. FTK=AF*exp(IIW*t)*gt*s;
  35. FT0_25K=AF*exp(IIW*t0_25)*gt0_25*s;
  36. FT0_50K=AF*exp(IIW*t0_50)*gt0_50*s;
  37. FT0_75K=AF*exp(IIW*t0_75)*gt0_75*s;
  38. FTK1=AF*exp(IIW*t1)*gt1*s;

  39. TUK=EXPDT*UK;

  40. %采用科茨公式求解非齐次项
  41. X1=7*EXPDT*FTK;
  42. X2=32*EXP0_75DT*FT0_25K;
  43. X3=12*EXP0_50DT*FT0_50K;
  44. X4=32*EXP0_25DT*FT0_75K;
  45. X5=7*FTK1;
  46. TKTK1=1/90*dt*(X1+X2+X3+X4+X5);

  47. UK1=TUK+TKTK1;

  48. disp(:,i)=UK1(1:n,:);
  49. velp(:,i)=UK1(n+1:2*n,:);
  50. UK=UK1;
  51. end
复制代码
  1. function [T]=EXPHDT(A,dt)
  2. %计算矩阵指数
  3. N=20;
  4. m=2^N;
  5. t=dt/m;
  6. I=eye(size(A));
  7. Ta=A*t+(A*t)^2*(I+(A*t)/3+(A*t)^2/12)/2;
  8. %计算指数矩阵T
  9. for i=1:20
  10.   Ta=2*Ta+Ta^2;
  11. end
  12. T=I+Ta;
复制代码


程序详细说明.doc

184.5 KB, 下载次数: 155

精细直接积分法的积分方法选择.doc

54.22 KB, 下载次数: 130

虚拟激励法程序.doc

35.5 KB, 下载次数: 138

评分

2

查看全部评分

 楼主| 发表于 2010-11-18 11:17 | 显示全部楼层
谢谢你哥们
发表于 2011-1-7 11:09 | 显示全部楼层
怎么又是时程又是频域的,实现起来还挺复杂的。怀疑他是否真的像说的那样高效。一般的方法只是频域积分就可以了。还希望明白的人说说
发表于 2011-1-11 13:48 | 显示全部楼层
林家浩老师刚来我们所讲虚拟激励法啊
正缺这个啊
谢谢
发表于 2011-1-11 13:58 | 显示全部楼层
精细直接积分法的积分方法选择.doc
打开出现乱码,请楼主重新上传
谢谢
发表于 2011-1-16 18:26 | 显示全部楼层
就是,打开时乱码啊
发表于 2011-1-17 01:12 | 显示全部楼层
呵呵,明天早上林家浩正要来我们这儿开个讲座。。。
发表于 2011-5-9 10:17 | 显示全部楼层
正在学习中
楼主有没关于结构失效的程序啊
发表于 2011-6-17 16:44 | 显示全部楼层
非常感谢,正需要!
发表于 2011-6-23 12:46 | 显示全部楼层
但是貌似不准确啊!
发表于 2012-3-19 07:41 | 显示全部楼层
回复 6 # linzhidan 的帖子

看了一下,要用pdf打开的。
发表于 2012-4-13 16:48 | 显示全部楼层
谢谢了,我转化过来了
发表于 2012-5-15 17:13 | 显示全部楼层
好东西,谢谢,学习一下
发表于 2012-9-28 10:19 | 显示全部楼层
好东西,谢谢楼主分享!!!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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