动力学方程数值解法：直接积分法（Newmark类）

weixin

　　Newmark类直接积分法原理

       下面首先来讨论一下Newmark类直接积分法的基本原理。

　　tn+1时刻的位移、速度、加速度可以表示为：

　　将式(3)代入到式(1)，(2)可以得到：

　　对上述式子进行进一步的整理：1. 可以直接略去高阶项;2. 用变权来调节。

　　假设其在tn+1时刻近似满足运动方程：

　　通过变换将速度和加速度用位移表示，代入运动方程，只剩tn+1时刻位移一个未知数，可得到：

　　根据参数选择的不同，它包含了三个经典的算法。

　　① Newmark平均加速度法

　　② Newmark线加速度法

　　③ 中心差分法

　　Newmark法

       Newmark法的一般步骤：
　　(1) 初始值计算
　　① 形成系统矩阵K，M和C

　　② 定初始值

　　③ 选择时间步长△t，参数γ 、β。并计算积分常数：

　　④ 形成等效刚度矩阵

　　⑤ 矩阵进行三角分解

　　(2) 对每一时间步
　　① 计算t+△t时刻的等效载荷

　　② 求解t+△t时刻的位移

　　③ 计算t+△t时刻的加速度和速度

　　Newmark法应用实例

       newmarkb子程序：

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

　　程序调用实例：三自由度系统的脉冲响应

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

　　Wilson-θ法

       Wilson- θ法的一般步骤：
　　(1) 初始值计算
　　① 形成系统矩阵K，M和C

　　② 定初始值

　　③ 选择时间步长，并计算积分常数

　　④ 形成等效刚度

　　⑤ 将等效刚度进行三角分解

　　(2) 对每一个时间步长
　　① 计算t+△t时刻的等效载荷

　　② 求解t+△t时刻的位移

　　③ 计算在t+△t时刻的加速度、速度和位移

　　Wilson-θ法程序代码

       仅供参考：

function [ua,ua1,ua2]=wilson(m,k,c,pt,u,u1,u2,n,theta,dertat)

% by swell2005

%---------------

a0=6/((theta*dertat)^2);

a1=3/(theta*dertat);

a2=2*a1;

a3=theta*dertat/2;

a4=a0/theta;

a5=-a2/theta;

a6=1-3/theta;

a7=dertat/2;

a8=dertat^2/6;

k=k+a0*m+a1*c;

[r,s]=ldl(k);

%----------------

for i=0:n,

  pta=pt+m*(a0*u+a2*u1+2*u2)+c*(a1*u+2*u1+a3*u2);

  uaa=r\pta;

  uaa=s\uaa;

  uaa=r'\uaa;

  ua2(:,i+1)=a4*(uaa-u)+a5*u1+a6*u2;

  ua1(:,i+1)=u1+a7*(ua2(:,i+1)+u2);

  ua(:,i+1)=u+dertat*u1+a8*(ua2(:,i+1)+2*u2);

  u=ua(:,i+1);

  u1=ua1(:,i+1);

  u2=ua2(:,i+1);

end

　　本文摘录自百度文库《结构动力学方程常用数值解法》一文，程序来源于声振论坛。