声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2963|回复: 7

[线性振动] newmark法和wilson_cita法为什么计算结果差很大

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

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

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

x
这是NewMark程序
  1. %NewMark Function
  2. clc,clear
  3. m1=341;                                       % Quality of the former landing gear
  4. m2=363;                                       % Quality of the back landing gear
  5. mc=24900;                                     % Quality of the body
  6. Jc=424378;                                    % Moment of inertia of the body
  7. Lz=5.805;                                     % Distance between the former and the back tyre(m
  8. a=0.91*Lz;                                    % Distance of the formwe tyre spread
  9. b=0.09*Lz;                                    % Distance of the back tyre spread
  10. k1=8.1123e5;                                  % Rigidity coefficient of the former tyre
  11. c1=3*1e4;                                     % Damp coefficient of the former tyre
  12. k3=1.6285e5;                                  % Rigidity coefficient of the former Amortine shoring
  13. c3=5*1e4;                                     % Damp coefficient of the former Amortine shoring
  14. k2=4.41046;                                   % Rigidity coefficient of the back tyre
  15. c2=3*1e4;                                     % Damp coefficient of the back tyre
  16. k4=7.3361e5;                                  % Rigidity coefficent of the back Amortine shoring
  17. c4=5*1e4;                                     % Damp coefficient of the back Amortine shoring
  18. g=9.81;
  19. M =  [(mc*b^2+Jc)/Lz^2  (mc*a*b-Jc)/Lz^2  0   0
  20.       (mc*a*b-Jc)/Lz^2  (mc*a^2+Jc)/Lz^2  0   0
  21.         0         0    m1  0
  22.         0         0    0   m2];
  23. C =  [c4   0  -c4    0
  24.       0    c3  0  -c3
  25.       -c4  0  c2+c4  0
  26.       -c3  0  0    c1+c3];
  27. K =  [k4   0  -k4    0
  28.       0    k3  0  -k3
  29.       -k4  0  k2+k4  0
  30.       -k3  0  0    k1+k3];
  31. u=[0 0 0 0]';
  32. v=[0 0 0 0]';
  33. a=[g 0 g g]';
  34. t(1)=0;               %time
  35. x(:,1)=u;             %distance
  36. x1(:,1)=v;            %speed
  37. x2(:,1)=a;            %accelerrate
  38. %newmar kcoefficient
  39. gama=0.5;
  40. dt=1e-2;
  41. delta=0.25;
  42. b0=1/(delta*dt^2);
  43. b1=gama/(delta*dt);
  44. b2=1/(delta*dt);
  45. b3=1/(2*delta)-1;
  46. b4=gama/delta-1;
  47. b5=dt*(gama/(2*delta)-1);
  48. b6=dt*(1-gama);
  49. b7=gama*dt;

  50. %等效刚度矩阵
  51. K1=K+b0*M+b1*C;
  52. t_max=10;       %计算时间总长
  53. i=1;
  54. t(1)=0.01;
  55. q=zeros(4,1);
  56. while t(i)<t_max
  57.       Q(:,i)=[0,0,0,0]';
  58.       q(:,i+1)=Q(:,i)+M*(b0*x(:,i)+b2*x1(:,i)+b3*x2(:,i))+C*(b1*x(:,i)+b4*x1(:,i)+b5*x2(:,i));
  59.       x(:,i+1)=inv(K1)*q(:,i+1);
  60.       x2(:,i+1)=b0*(x(:,i+1)-x(:,i))-b2*x1(:,i)-b3*x2(:,i);
  61.       x1(:,i+1)=b1*(x(:,i+1)-x(:,i))-b4*x1(:,i)-b5*x2(:,i);
  62.       i=i+1;
  63.       t(i)=t(i-1)+dt;
  64. end
复制代码


这是wilson_cita程序
  1. %wilson_cita Function
  2. clc,clear
  3. m1=341;                                       % Quality of the former landing gear
  4. m2=363;                                       % Quality of the back landing gear
  5. mc=24900;                                     % Quality of the body
  6. Jc=424378;                                    % Moment of inertia of the body
  7. Lz=5.805;                                     % Distance between the former and the back tyre(m
  8. a=0.91*Lz;                                    % Distance of the formwe tyre spread
  9. b=0.09*Lz;                                    % Distance of the back tyre spread
  10. k1=8.1123e5;                                  % Rigidity coefficient of the former tyre
  11. c1=3*1e4;                                     % Damp coefficient of the former tyre
  12. k3=1.6285e5;                                  % Rigidity coefficient of the former Amortine shoring
  13. c3=5*1e4;                                     % Damp coefficient of the former Amortine shoring
  14. k2=4.41046;                                   % Rigidity coefficient of the back tyre
  15. c2=3*1e4;                                     % Damp coefficient of the back tyre
  16. k4=7.3361e5;                                  % Rigidity coefficent of the back Amortine shoring
  17. c4=5*1e4;                                     % Damp coefficient of the back Amortine shoring
  18. g=9.81;
  19. M =  [(mc*b^2+Jc)/Lz^2  (mc*a*b-Jc)/Lz^2  0   0
  20.       (mc*a*b-Jc)/Lz^2  (mc*a^2+Jc)/Lz^2  0   0
  21.         0         0    m1  0
  22.         0         0    0   m2];
  23. C =  [c4   0  -c4    0
  24.       0    c3  0  -c3
  25.       -c4  0  c2+c4  0
  26.       -c3  0  0    c1+c3];
  27. K =  [k4   0  -k4    0
  28.       0    k3  0  -k3
  29.       -k4  0  k2+k4  0
  30.       -k3  0  0    k1+k3];
  31. u=[0 0 0 0]';
  32. v=[0 0 0 0]';
  33. a1=[g 0 g g]';
  34. t(1)=0;               %time
  35. x(:,1)=u;             %distance
  36. x1(:,1)=v;            %speed
  37. x2(:,1)=a1;            %accelerate
  38. %cita coeffi
  39. cita=1.42085;
  40. dt=1e-2;
  41. s=cita*dt;
  42. %等效刚度矩阵
  43. K1=K+6/s^2*M+3/s*C;
  44. t_max=10;       %full time
  45. i=1;
  46. t(1)=0.01;
  47. while t(i)<t_max
  48.       Q=[mc*g,0,m1*g,m2*g]';
  49.       gt=M*(6/s^2*x(:,i)+6/s*x1(:,i)+2*x2(:,i))+C*(3/s*x(:,i)+2*x1(:,i)+s/2*x2(:,i))+Q;
  50.       us=inv(K1)*gt;
  51.       u2s=6/s^2*(us-x(:,i))-6/s*x1(:,i)-2*x2(:,i);
  52.       u1s=3/s*(us-x(:,i))-2*x1(:,i)-s/2*x2(:,i);
  53.       x2(:,i+1)=6*dt/s^3*(us-x(:,i))-6*dt/s^2*x1(:,i)+(1-3/cita)*x2(:,i);
  54.       x1(:,i+1)=x1(:,i)+dt/2*(u2s-x2(:,i));
  55.       x(:,i+1)=x(:,i)+dt*x1(:,i)+dt^2/6*(u2s+2*x2(:,i));
  56.       i=i+1;
  57.       t(i)=t(i-1)+dt;
  58. end
复制代码


两个程序计算结果相差很大,希望大家帮忙看看。
其中NewMark程序算出稳定解
而wilson_cita算出的竟然不稳定

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2006-11-15 13:27 | 显示全部楼层

各有优劣

 楼主| 发表于 2006-11-15 15:02 | 显示全部楼层
那么你有没有wilson_cita的Matlab程序啊
我觉的我的程序好像没有错
发表于 2006-11-15 22:00 | 显示全部楼层
q(:,i+1)=Q(:,i)+M*(b0*x(:,i)+b2*x1(:,i)+b3*x2(:,i))+C*(b1*x(:,i)+b4*x1(:,i)+b5*x2(:,i));
好像这个语句有点问题,b0*x(:,i)应该是b6*x(:,i)吧
发表于 2006-11-15 22:45 | 显示全部楼层
将你的程序改正如下。
      x2(:,i+1)=6*dt/s^3*(us-x(:,i))-6*dt/s^2*x1(:,i)+(1-3/cita)*x2(:,i);
      x1(:,i+1)=x1(:,i)+dt/2*(x2(:,i+1)+x2(:,i))
觉得
  u1s=3/s*(us-x(:,i))-2*x1(:,i)-s/2*x2(:,i);
此行程序没有作用

评分

1

查看全部评分

发表于 2006-11-17 20:11 | 显示全部楼层
威尔逊法:plot(x2(2,:))得到如下 图象,好象收敛
发表于 2006-11-17 20:24 | 显示全部楼层

图象

回复
1.JPG
 楼主| 发表于 2006-11-18 18:47 | 显示全部楼层
谢谢了,heyong 兄
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-3 05:39 , Processed in 0.083442 second(s), 27 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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