声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1533|回复: 0

[动力学和稳定性] 动力学方程数值解法:直接积分法(Newmark类)

[复制链接]
发表于 2018-5-11 16:14 | 显示全部楼层 |阅读模式

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

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

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

  tn+1时刻的位移、速度、加速度可以表示为:
1.png
  将式(3)代入到式(1),(2)可以得到:
2.png
  对上述式子进行进一步的整理:1. 可以直接略去高阶项;2. 用变权来调节。
3.png
  假设其在tn+1时刻近似满足运动方程:
4.png
  通过变换将速度和加速度用位移表示,代入运动方程,只剩tn+1时刻位移一个未知数,可得到:
5.png
  根据参数选择的不同,它包含了三个经典的算法。

  ① Newmark平均加速度法
6.png
  ② Newmark线加速度法
7.png
  ③ 中心差分法
8.png

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

  ② 定初始值 0.png

  ③ 选择时间步长△t,参数γ 、β。并计算积分常数:
9.png
  ④ 形成等效刚度矩阵
10.png
  ⑤ 矩阵进行三角分解
11.png

  (2) 对每一时间步
  ① 计算t+△t时刻的等效载荷
12.png
  ② 求解t+△t时刻的位移
13.png
  ③ 计算t+△t时刻的加速度和速度
14.png

  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

  ② 定初始值 15.png

  ③ 选择时间步长,并计算积分常数
16.png
  ④ 形成等效刚度
17.png
  ⑤ 将等效刚度进行三角分解
18.png

  (2) 对每一个时间步长
  ① 计算t+△t时刻的等效载荷
19.png
  ② 求解t+△t时刻的位移
20.png
  ③ 计算在t+△t时刻的加速度、速度和位移
21.png

  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

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

回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-3-28 19:33 , Processed in 0.164071 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表