〔求助〕如何用MATLAB计算结构响应？

realyyy

根据动力学基本方程：M*ddX+C*dX+K*X=F
现已知M、C、K及F，如何用MATLAB计算结构响应呢？(X、dX、ddX分别为位移响应、速度响应、加速度响应)
先行谢过！
（补充一下：F为静力。）
（我本来想利用频响函数H（w）=inv(M-j*C/w-K/w^2)来算，但那样需要知道激励的频率，现在F只是静力啊。）

[ 本帖最后由 ChaChing 于 2009-10-19 16:09 编辑 ]

siyanger

静力是什么意思啊？你的意思就是求解这个微分方程吧？

feifeifool

可以尝试用Ｎemark方法或者中心差分法，威尔逊－theta法求解．
静力 只会产生一个特解，对应产生－静变形，关键还得看它的初始条件．

[ 本帖最后由 ChaChing 于 2009-10-19 16:35 编辑 ]

feifeifool

静力 只会产生一个特解，对应产生－静变形，关键还得看它的初始条件．

realyyy

感谢各位！
就是求解微分方程。
关键是其中的系数M、C、K以及X、F等都是矩阵啊！用通常的求解微分方程的方法不行的。(初始条件X0=dX0=ddX0=0)

feifeifool

对矩阵也是适用的，没有什么问题。初始向量为０的话，只作用一个静力是不会产生动态响应的。
你可以考虑用newmark方法或者用Runge Kutta（m，k维数不是很大的话）方法，这个不是很难的，我用newmark法算过二十几阶的。

realyyy

感谢feifeifool！
我试试。回头给大家结果。

realyyy

NEWMARK计算结构响应的MATLAB程序已经编制，出现了新问题！详见：http://forum.vibunion.com/forum/ ... ighlight=%2Brealyyy

feifeifool

clc;clear;close;
M=[1 1 1];M=diag(M);
C=diag([0.1,0.0,0.0]);
K=9*[1 1 1];K=diag(K);
u=[0 0 0]';
v=[0 0 0]';
a=[10 0 0]';
t(1)=0;               %时间
x(:,1)=u;             %位移
x1(:,1)=v;            %速度
x2(:,1)=a;            %加速度
%newmark参数
gama=0.5;
dt=1e-1;
delta=0.25;
a0=1/(delta*dt^2);
a1=gama/(delta*dt);
a2=1/(delta*dt);
a3=1/(2*delta)-1;
a4=gama/delta-1;
a5=dt*(gama/(2*delta)-1);
a6=dt*(1-gama);
a7=gama*dt;
q=zeros(5,1);
%等效刚度矩阵
K1=K+a0*M+a1*C;
% a0=zeros(5,1);
t_max=100;       %计算时间总长
i=1;
t(1)=0;
q=zeros(3,1);
while t(i)<t_max
      f=0*sin(3*t(i));
      Q=[f,0,0]';
      q(:,i+1)=Q+M*(a0*x(:,i)+a2*x1(:,i)+a3*x2(:,i))+C*(a1*x(:,i)+a4*x1(:,i)+a5*x2(:,i));
      x(:,i+1)=inv(K1)*q(:,i+1);
      x2(:,i+1)=a0*(x(:,i+1)-x(:,i))-a2*x1(:,i)-a3*x2(:,i);
      x1(:,i+1)=a1*(x(:,i+1)-x(:,i))-a4*x1(:,i)-a5*x2(:,i);
      i=i+1;
      t(i)=t(i-1)+dt;
end  
这个应该没有问题你可以试一下！你的程序好像newmark中的各个参数就有问题，其他的问题我也不是很清楚，你可以参考一下这个程序改一改。此外这种自由度比较小的拿matlab的ode45求解，很不错的方法。
你可以看一下，很简单。

realyyy

感谢feifeifool！
采用你的程序后运行正确通过！
（我原来是照《机械振动基础》(胡海岩主编)这书上的公式抄的啊，一个字母都不差，咋就不行呢？郁闷。。。）

[ 本帖最后由 realyyy 于 2006-8-14 15:32 编辑 ]

feifeifool

那你看看是不是写程序的时候写错了。我的程序是参照《结构动力学》邹经湘的，这个肯定没问题。
呵呵，不用那么客气

bettylong

这个程序很好,有好东西记得和大家一起分享哦!谢谢各位啦

ISHMII

谢谢，好用

gungungun

:handshake   3x！

pp786195

怎么判断求出的结果是否收敛？