Newmark数值分析法matlab程序

clc

clear

M=[1,2,3,4];

m=diag(M); %质量对角线矩阵

k=[ 800 -800 0 0 ;

-800 2400 -1600 0 ;

0 -1600 4800 -3200 ;

0 0 -3200 8000 ]; %计算刚度矩阵

c=0.05m+0.02k;%c=2m0.05*sqrt(k/m) ;

f0= 100;

t1=5; %力的作用时间

nt=2000; %分2000步完成

dt=0.01; %时间步长

alfa=0.25; %γ=0.25

beta=0.5; %β=0.5

a0=1/alfa/dt/dt; %

a1=beta/alfa/dt; %

a2=1/alfa/dt; %

a3=1/2/alfa-1; %

a4=beta/alfa-1; %

a5=dt/2*(beta/alfa-2); %

a6=dt*(1-beta); %

a7=dtbeta; %

d=zeros(4,nt); %初位移为0

v=zeros(4,nt); % 初速度为0

a=zeros(4,nt); % 初加速度为0

for i=2:nt

t=(i-1)dt;

if (t<t1),

f=[f0sin(4pit/t1);0;0;0]; %f=[300sin(6pit)-50cos(3pit);0;0;0]; %力作用时间内对结构进行加载

else

f=[0;0;0;0]; %力作用时间外结构不受力

end

ke=k+a0m+a1c; % 有效刚度矩阵

fe=f+m(a0d(:,i-1)+a2v(:,i-1)+a3a(:,i-1))+c(a1d(:,i-1)+a4v(:,i-1)+a5a(:,i-1));

% t+dt时刻的有效荷载

d(:,i)=ke\fe;%d(:,i)=inv(ke)fe; %求解t+dt时刻的位移

a(:,i)=a0(d(:,i)-d(:,i-1))-a2v(:,i-1)-a3a(:,i-1); %计算t+dt时刻的加速度

v(:,i)=v(:,i-1)+a6a(:,i-1)+a7*a(:,i); %计算t+dt时刻的速度

end

T=(0:dt:19.99); %离散系统dt为采样周期 19.99为终端时间

close all

figure %控制窗口数量

plot(T,d(1,:)) %绘制位移函数图像

title(‘各质点位移总图’) %添加标题为各质点位移总图

xlabel(‘s’) %对x轴进行标注为时间（s）

ylabel(‘m’) %对y轴进行标注为位移（m）

grid %显示画图中的个网线

figure

plot(T,a) %绘制加速度函数图像

title(‘加速度总图’) %添加标题为各质点加速度总图

legend(‘a1’,‘a2’,‘a3’,‘a4’)

xlabel(‘s’) %对x轴进行标注为时间（s）

ylabel(‘m/s^2’) %对y轴进行标注为加速度（m/s2）

grid

figure

plot(T,v) %绘制速度函数图像

title(‘速度总图’) %添加标题为各质点速度总图

legend(‘v1’,‘v2’,‘v3’,‘v4’)

xlabel(‘s’) %对x轴进行标注为时间（s）

ylabel(‘m/s)’) %对y轴进行标注为加速度（m/s）

grid