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

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

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

function =newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength)
% newmark-beta method
%         obtain the response of the dynamic system
%         =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 =newmarkb(M,K,C,N,P,x0,v0,a0,dt,RecordLength)
% newmark-beta method
%         obtain the response of the dynamic system
%         =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 发表 http://www.chinavib.com/forum/images/common/back.gif
谢谢eight老兄指点啊，那是不是改成下面这样就好些了？

function =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;
=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);

=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

不错的帖子，例子都有了{:{05}:}

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

sogooda 发表于 2007-11-17 14:27 static/image/common/back.gif
例子1：三自由度系统的脉冲响应

位移：


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

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

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

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

谢谢分享！

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

7楼fft变换结果对称的有问题吗？
页: [1]
查看完整版本: [原创]Newmark-beta法求解动力系统函数