newmark方法

newmark方法

方法介绍 http://en.wikilib.com/wiki/Newmark-beta_method

 

function [X,t] = newmark(M,C,K,f,dt,gam,beta,Xi,Xdi)
% % 输入参数说明：
% [M] = 质量 (nxn)
% [C] = 阻尼矩阵 (nxn)
% [K] = 刚度矩阵 (nxn)
% {f} = 激励力矩阵 (nxm), f的第k列是t(k)时刻的激励力，m-1为总的计算步数
% dt = 步长
% beta= Newmark 常数 (1/6 or 1/4 usually)
% gam = Newmark 常数 (1/2)
% Xi = 初始位移向量 (nx1)
% Xdi = 初始速度向量 (nx1)
% % 输出参数说明：
% {t} = 时间 (mx1)
% [X] = 响应 (nxm)，X的第k列为t(k)时刻的位移响应

%================================================================
n = size(M,1);
m = size(f,2);
t = zeros(m,1);
X = zeros(n,m);
Xd = zeros(n,m);
Xdd = zeros(n,m);

% 系数
c0 = 1/(beta*dt*dt) ;
c1 = gam/(beta*dt) ;
c2 = 1/(beta*dt) ;
c3 = 1/(beta*2) - 1 ;
c4 = gam/beta - 1 ;
c5 = 0.5*dt*(gam/beta - 2 ) ;
c6 = dt*(1 - gam ) ;
c7 = dt* gam;

% 初始刚度
Keff = c0*M + c1*C + K ;
Kinv = inv(Keff) ;

% 初始加速度
R0 = f(:,1);
Xddi= inv(M)*(R0 - K*Xi - C*Xdi);

% 第一步
f(:,1) = f(:,1) + M*(c0*Xi+c2*Xdi+c3*Xddi) ...
    +C*(c1*Xi+c4*Xdi+c5*Xddi);
X(:,1) = Kinv*f(:,1);
Xdd(:,1)= c0*(X(:,1)-Xi) - c2*Xdi - c3*Xddi ;
Xd(:,1) = Xdi + c6*Xddi + c7*Xdd(:,1);
t(1) = 0;

% 后续步
% ========================
for i=1:m-1
    t(i+1) = t(i) + dt;
    f(:,i+1) = f(:,i+1) + M * ...
        (c0*X(:,i)+c2*Xd(:,i)+c3*Xdd(:,i)) ...
        + C*(c1*X(:,i)+c4*Xd(:,i)+c5*Xdd(:,i)) ;
    X(:,i+1) = Kinv*f(:,i+1) ;
    Xdd(:,i+1)= c0*(X(:,i+1)-X(:,i))-c2*Xd(:,i)-c3*Xdd(:,i) ;
    Xd(:,i+1) = Xd(:,i)+c6*Xdd(:,i)+c7*Xdd(:,i+1) ;
end;

程序转自xyz999

newmark方法

方法介绍 http://en.wikilib.com/wiki/Newmark-beta_method

 

function [X,t] = newmark(M,C,K,f,dt,gam,beta,Xi,Xdi)
% % 输入参数说明：
% [M] = 质量 (nxn)
% [C] = 阻尼矩阵 (nxn)
% [K] = 刚度矩阵 (nxn)
% {f} = 激励力矩阵 (nxm), f的第k列是t(k)时刻的激励力，m-1为总的计算步数
% dt = 步长
% beta= Newmark 常数 (1/6 or 1/4 usually)
% gam = Newmark 常数 (1/2)
% Xi = 初始位移向量 (nx1)
% Xdi = 初始速度向量 (nx1)
% % 输出参数说明：
% {t} = 时间 (mx1)
% [X] = 响应 (nxm)，X的第k列为t(k)时刻的位移响应

%================================================================
n = size(M,1);
m = size(f,2);
t = zeros(m,1);
X = zeros(n,m);
Xd = zeros(n,m);
Xdd = zeros(n,m);

% 系数
c0 = 1/(beta*dt*dt) ;
c1 = gam/(beta*dt) ;
c2 = 1/(beta*dt) ;
c3 = 1/(beta*2) - 1 ;
c4 = gam/beta - 1 ;
c5 = 0.5*dt*(gam/beta - 2 ) ;
c6 = dt*(1 - gam ) ;
c7 = dt* gam;

% 初始刚度
Keff = c0*M + c1*C + K ;
Kinv = inv(Keff) ;

% 初始加速度
R0 = f(:,1);
Xddi= inv(M)*(R0 - K*Xi - C*Xdi);

% 第一步
f(:,1) = f(:,1) + M*(c0*Xi+c2*Xdi+c3*Xddi) ...
    +C*(c1*Xi+c4*Xdi+c5*Xddi);
X(:,1) = Kinv*f(:,1);
Xdd(:,1)= c0*(X(:,1)-Xi) - c2*Xdi - c3*Xddi ;
Xd(:,1) = Xdi + c6*Xddi + c7*Xdd(:,1);
t(1) = 0;

% 后续步
% ========================
for i=1:m-1
    t(i+1) = t(i) + dt;
    f(:,i+1) = f(:,i+1) + M * ...
        (c0*X(:,i)+c2*Xd(:,i)+c3*Xdd(:,i)) ...
        + C*(c1*X(:,i)+c4*Xd(:,i)+c5*Xdd(:,i)) ;
    X(:,i+1) = Kinv*f(:,i+1) ;
    Xdd(:,i+1)= c0*(X(:,i+1)-X(:,i))-c2*Xd(:,i)-c3*Xdd(:,i) ;
    Xd(:,i+1) = Xd(:,i)+c6*Xdd(:,i)+c7*Xdd(:,i+1) ;
end;

程序转自xyz999