请教一个三自由度质量－弹簧振动问题

tianlei2005

一个常见的三自由度质量－弹簧系统，其动力学方程为：
[M]{x''}+[K]{x}={F}
质量、刚度和激励矩阵分别为：
M=diag([1;1;1]);k=[3 -1 0;-1 2 -1;0 -1 3];F={sin(3*t);0;0};
我分别用模态叠加法和Runge-Kutta算法求解，但是两种解法得到的结果却不相同，请问这是什么原因，何种方法才是正确的。
这是我用模态叠加法的代码：
%利用模态叠加法分析扭转振动的实例
%定义符号,a、b为阻尼的比例系数
syms t;
%齿轮系统参数
m1=1;m2=1;m3=1;k1=2;k2=2;k3=1;k4=2;
M=diag([m1;m2;m3]);
% K=[k1+k2 -k2 0;
%     -k2 k2+k3 -k3;
%     0 -k3 k3+k4];
K=[3 -1 0;-1 2 -1;0 -1 3];
F1=1;F2=0;F3=0;w=3;
w0(1)=w;w0(2)=0*w;w0(3)=0*w;
f1=F1*sin(w0(1)*t);f2=F2*sin(w0(2)*t);f3=F3*sin(w0(3)*t);
%求振动系统的正则振型矩阵和固有频率
[N,V]=eig(K,M);
ModeValue=sqrt(diag(V));
w1=ModeValue(1);
w2=ModeValue(2);
w3=ModeValue(3);
%定义振动系统的初始条件
x0=[0 0 0]';
X0=N'*M*x0;
    a1=X0(1);
    a2=X0(2);
    a3=X0(3);
%正则变换，求出正则坐标下的激励力
Q=[f1;f2;f3];
Q_Fai=N'*[f1;f2;f3];
Qn=N'*diag([F1;F2;F3]);
% Qn1=N'*[f1;0;0];
% Qn2=N'*[0;f2;0];
% Qn3=N'*[0;0;f3];
Q_Fai1=N'*[f1;0;0];
Q_Fai2=N'*[0;f2;0];
Q_Fai3=N'*[0;0;f3];
%求出在正则坐标下的响应
n=length(M);
%求系统的无阻尼受迫振动解
for i=1:3
    b1(i)=1/(ModeValue(i)^2-w^2);
    b2(i)=1/(ModeValue(i)^2-(3*w)^2);
    b3(i)=1/(ModeValue(i)^2-(8*w)^2);
end
Xn1=b1'.*Q_Fai1;
Xn2=b2'.*Q_Fai2;
Xn3=b3'.*Q_Fai3;
x1=N*Xn1;
x2=N*Xn2;
x3=N*Xn3;
x=x1+x2+x3;
% x=C*Xn;
% Qn=N'*[F1;F2;F3];
x=-0.088.*[1;2;1].*sin(3*t)-2.63.*[-1.732;0;1.732].*sin(3*t)+0.21.*[1.414;-1.414;1.414].*sin(3*t)
%定义仿真时间和采样点数
t=0:0.01:100;
%对结果进行fft变换
u1=eval(x(1));
u2=eval(x(2));
u3=eval(x(3));（省去后面的画图和FFt变换部分），


下面是调用ode45函数的代码
%test4.m
function f=test4(t,y);
m1=1;m2=1;m3=2;k=1;
U=[0 1 0 0 0 0;
    -3 0 1 0 0 0;
    0 0 0 1 0 0;
    1 0 -2 0 1 0;
    0 0 0 0 0 1;
    0 0 1 0 -3 0];
f=U*y+[0 sin(3*t) 0 0 0 0]';

%test4Result.m
[t,y]=ode45('test4',[0:0.01:100],[0 0 0 0 0 0]');
u1=y(:,1);
u2=y(:,3);
u3=y(:,5);
(后面省去画图和fft变换部分)

vib

你研究得很仔细，向你学习，我们学机械振动时我都没思考过这些问题，

xzf198347

这跟我要研究的内容挺相关的，请问是用matlab仿真的吗？我一定要好好学习:lol

VibrationMaster

模态叠加法中一般不考虑初条件，而龙个哭他法的解包含初始条件的贡献

hyl2323

我也用matlab编写过，我用的是模态叠加法，已经验证程序正确，我想你的不同应该是龙哥窟他法的初始条件引起的，模态叠加也要考虑初始条件的。有时间我用我程序来帮你算下。

VibrationMaster

谢谢你，模态叠加法中可以考虑进初条件。但是龙个哭他法肯定要有初条件。
比较2。jpg和4。jpg就可以知道对应与瞬态响应的三个固有频率谱峰没有出现在模态叠加法中，所以我推测楼主在模态叠加法中没有考虑初条件。
将来有问题，一定找你要程序。

keyi1224

下面是调用ode45函数的代码
%test4.m
function f=test4(t,y);
m1=1;m2=1;m3=2;k=1;
U=[0 1 0 0 0 0;
    -3 0 1 0 0 0;
    0 0 0 1 0 0;
    1 0 -2 0 1 0;
    0 0 0 0 0 1;
    0 0 1 0 -3 0];
f=U*y+[0 sin(3*t) 0 0 0 0]';

这里的U什么意思啊？？？？

whecust

k1=2;k2=2;k3=1;k4=2;

% K=[k1+k2 -k2 0;
%     -k2 k2+k3 -k3;
%     0 -k3 k3+k4];
K=[3 -1 0;-1 2 -1;0 -1 3];
好像k算错了？

vibrater

你的图象没有纵横坐标，也没有说明求出来的是什莫？

huaideggs

学习，学习:@)

cumtjack

你好，能不能把你程序分享一下呢，正在做轮文，希望能够借鉴一下

cumtjack

楼主，能不能把程序不全呢？
我正在做毕业论文，想参考一下

cumtjack

请问楼主，你在用eig求解固有频率和振型时，振型满足正交性吗
我也是用这个函数求得，结果振型不满足正交性，不知道对不对
希望大家能够讨论一下

WSYcxl

怎么大家都说Runge-Kutta法需要初始条件那？
你们仔细的分析一下龙格库塔积分形式(四阶为例)：
dX=EX，X为构造的状态量；X(i+1)=X(1)+[K(1)+2K(2)+2K(3)+K(4)]/6，(小括号内为下标，下同)，由于K(1）包含了X，而X是只与M、C、K有关的。所以龙格库塔法是自启动的，是不需要初始条件的。
ode45的内核是什么我不知道，你可以把初始条件设为[0]试试看，看看结果是不是一样的。
关于这个我印象比较深刻。我硕士答辩时，为了凑论文字数，用龙格库塔法时过程写的很详细。然后一个老师问：你这个我怎么没看到初始条件那，难道这个系统无缘无故的自己就振动起来了？然后我解释了一番，但因为他是搞信号检测、探伤方面的，年纪也很大了，什么数值模拟他不知道。幸亏大家平时都熟，最后还是老板给我打了圆场。

[ 本帖最后由 WSYcxl 于 2009-3-16 21:45 编辑 ]

a043588

% x=C*Xn;
% Qn=N'*[F1;F2;F3];
x=-0.088.*[1;2;1].*sin(3*t)-2.63.*[-1.732;0;1.732].*sin(3*t)+0.21.*[1.414;-1.414;1.414].*sin(3*t)X是怎么算出来的请详细说明下呗