[原创]Newmark-beta法求解动力系统函数

sogooda · 发表于 2007-11-5 10:33

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

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

x

本帖最后由 ChaChing 于 2010-8-25 22:30 编辑

function [x,v,a]=newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength)
% newmark-beta method
%          obtain the response of the dynamic system
%          [x,v,a]=newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength)
%          M - mass matrix
%          K - stiffness matrix
%          C - damping matrix
%          N - DOF
%          P - loads
%          x0 - initial displacement
%          v0 - initial velocity
%          a0 - initial acceleration
%          dt - interval
%          RecordLength - number of sampling points
x(:,1)=x0; v(:,1)=v0; a(:,1)=a0; deta=0.50; alpha=0.25;
a0=1/alpha/dt^2; a1=deta/alpha/dt; a2=1/alpha/dt; a3=1/2/alpha-1;
a4=deta/alpha-1; a5=dt*(deta/alpha-2)/2; a6=dt*(1-deta); a7=deta*dt;
K_=K+a0*M+a1*C; iK=inv(K_);

for i=1:RecordLength-1
P_(:,i+1)=P(:,i+1)+M*(a0*x(:,i)+a2*v(:,i)+a3*a(:,i))+C*(a1*x(:,i)+a4*v(:,i)+a5*a(:,i));
x(:,i+1)=iK*P_(:,i+1);
a(:,i+1)=a0*(x(:,i+1)-x(:,i))-a2*v(:,i)-a3*a(:,i);
v(:,i+1)=v(:,i)+a6*a(:,i)+a7*a(:,i+1);
end

eight · 发表于 2007-11-5 10:42

本帖最后由 ChaChing 于 2010-8-25 22:28 编辑

谢谢你的分享。有个建议：应对 x、a、v 这几个矩阵先预分配空间

sogooda · 发表于 2007-11-5 14:23

本帖最后由 ChaChing 于 2010-8-25 22:34 编辑

谢谢eight老兄指点啊，那是不是改成下面这样就好些了？

function [x,v,a]=newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength)
% newmark-beta method
%          obtain the response of the dynamic system
%          [x,v,a]=newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength)
%          M - mass matrix
%          K - stiffness matrix
%          C - damping matrix
%          N - DOF
%          P - loads
%          x0 - initial displacement
%          v0 - initial velocity
%          a0 - initial acceleration
%          dt - interval
%          RecordLength - number of sampling points
x=zeros(N,RecordLength); v=zeros(N,RecordLength); a=zeros(N,RecordLength);

x(:,1)=x0; v(:,1)=v0; a(:,1)=a0; deta=0.50; alpha=0.25;
a0=1/alpha/dt^2; a1=deta/alpha/dt; a2=1/alpha/dt; a3=1/2/alpha-1;
a4=deta/alpha-1; a5=dt*(deta/alpha-2)/2; a6=dt*(1-deta); a7=deta*dt;
K_=K+a0*M+a1*C; iK=inv(K_);

for i=1:RecordLength-1
P_(:,i+1)=P(:,i+1)+M*(a0*x(:,i)+a2*v(:,i)+a3*a(:,i))+C*(a1*x(:,i)+a4*v(:,i)+a5*a(:,i));
x(:,i+1)=iK*P_(:,i+1);
a(:,i+1)=a0*(x(:,i+1)-x(:,i))-a2*v(:,i)-a3*a(:,i);
v(:,i+1)=v(:,i)+a6*a(:,i)+a7*a(:,i+1);
end

eight · 发表于 2007-11-5 19:17

原帖由 sogooda 于 2007-11-5 14:23 发表

 登录/注册后可看大图

谢谢eight老兄指点啊，那是不是改成下面这样就好些了？

function [x,v,a]=newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength)
% newmark-beta method
% obtain the response of the dynamic system
...

P_ 也没有定义就使用了，类似这种问题其实你用高版本的 matlab （2006a以上）就可以在编辑窗口发现了

sogooda · 发表于 2007-11-6 10:41

这样啊？。。。呵呵。。。我用的是7.04，没有提示。不过又学到了一招，呵呵。。。谢谢eight版主

vib · 发表于 2007-11-16 12:37

最好再给个例子就更好了，呵呵

sogooda · 发表于 2007-11-17 14:27

例子1：三自由度系统的脉冲响应

clear;clc
N=3;
%Assemble mass matrix
m=ones(N,1)*4e5;
M=diag(m);
%Assemble stiffness matrix
k=ones(N,1)*2e8;
K(1,1)=k(1)+k(2);
K(1,2)=-k(2);
for i=2:N-1
K(i,i-1)=-k(i);
K(i,i)=k(i)+k(i+1);
K(i,i+1)=-k(i+1);
end
K(N,N-1)=-k(N);
K(N,N)=k(N);
%Assemble damping matrix
ac=0.7334;
bc=0.0026;
C=ac*M+bc*K;
x0(:,1)=0*ones(N,1);
v0(:,1)=0*ones(N,1);
a0(:,1)=0*ones(N,1);
RecordLength=1000;
dt=0.01;
P=zeros(N,RecordLength);
P(1,2)=1e5;
[x,v,a]=newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength);

复制代码

位移：

速度

加速度

例子2：三自由度系统的自由响应

clear;clc
N=3;
%Assemble mass matrix
m=ones(N,1)*4e5;
M=diag(m);
%Assemble stiffness matrix
k=ones(N,1)*2e8;
K(1,1)=k(1)+k(2);
K(1,2)=-k(2);
for i=2:N-1
K(i,i-1)=-k(i);
K(i,i)=k(i)+k(i+1);
K(i,i+1)=-k(i+1);
end
K(N,N-1)=-k(N);
K(N,N)=k(N);
%Assemble damping matrix
ac=0.7334;
bc=0.0026;
C=ac*M+bc*K;
x0(:,1)=0*ones(N,1);
v0(:,1)=0*ones(N,1);
a0(:,1)=0*ones(N,1);
v0(1,1)=0.05;
RecordLength=1000;
dt=0.01;
P=zeros(N,RecordLength);
[x,v,a]=newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength);

复制代码

结果图
位移：

速度

加速度

自由度1的位移fft变换

jupiter1281 · 发表于 2008-12-29 13:20

:@) :@) hao

leeking30 · 发表于 2008-12-30 04:55

初始加速度a 应该能自动算出来，不需要赋值吧！外力是知道的啊？

skyliu · 发表于 2013-10-16 16:30

不错的帖子，例子都有了

zcy09011391 · 发表于 2013-10-18 13:44

sogooda 发表于 2007-11-17 14:27

登录/注册后可看大图

例子1：三自由度系统的脉冲响应

位移：

好啊！！！！！！！！！！！

lihaitao123 · 发表于 2013-10-18 16:14

挺好的例子，如果刚度中出现非线性项的话不知道该怎么处理？

dgyezw007 · 发表于 2014-2-10 17:52

谢谢分享！

dollfish000 · 发表于 2014-6-13 14:59

7楼fft变换结果对称的有问题吗？

账号		自动登录	找回密码
密码			我要加入

[共享资源] [原创]Newmark-beta法求解动力系统函数

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

评分

本帖被以下淘专辑推荐:

两个例子：多自由度系统脉冲响应和自由响应

评分

评分