动力学方程数值解法：直接积分法（Newmark类）
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 =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 pointsx=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;clcN=3;%Assemble mass matrixm=ones(N,1)*4e5;M=diag(m);%Assemble stiffness matrixk=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);endK(N,N-1)=-k(N);K(N,N)=k(N);%Assemble damping matrixac=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);
Wilson-θ法 Wilson- θ法的一般步骤：
(1) 初始值计算
① 形成系统矩阵K，M和C
② 定初始值
③ 选择时间步长，并计算积分常数
④ 形成等效刚度
⑤ 将等效刚度进行三角分解
(2) 对每一个时间步长
① 计算t+△t时刻的等效载荷
② 求解t+△t时刻的位移
③ 计算在t+△t时刻的加速度、速度和位移
Wilson-θ法程序代码 仅供参考：
function =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; =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
本文摘录自百度文库《结构动力学方程常用数值解法》一文，程序来源于声振论坛。
页:
[1]