马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
如题,自己按照Newmark-beta法编制了一个单自由度的时程分析程序,但算出来的结果和sap2000的不符,而且反应加速度值很大,不知道问题出在哪里,源程序如下:
clear;clc;
% 结构模型初始参数----------------------------------------------------------
m=3e3; %质量(单位:kg)
k=1e6; %刚度((单位:N/m))
kesai=0.05; %阻尼比取0.05
c=2*kesai*sqrt(k*m); %阻尼系数
% 读取地震波数据------------------------------------------------------------
acc=textread('D:\研究生课题资料\科研课题资料\处理后的smc文件\51WCW_90_chnua370295.smc_090501.a','%f','headerlines',56);
PGA_Max=max(abs(acc)) %最大地面加速度绝对值
% Newmark-beta法的基本参数--------------------------------------------------
beta=1/6; gama=0.5; %按线性加速度法计算更接近真实结果,故取此组参数
dt=0.02; %地震加速度时程波记录时间间隔
b1=1/(beta*dt^2); b2=1/(beta*dt); b3=1-1/(2*beta); %计算参数
b4=gama/(beta*dt); b5=gama/beta-1; b6=(1-gama/(2*beta))*dt;
ke=k+m*b1+c*b4; %等效刚度
% 设定结构初始状态为零,生成向量空间存储计算值---------------------------------
u=zeros(100/dt,1); v=zeros(100/dt,1); a=zeros(100/dt,1);
% Newmark-beta法的主计算程序------------------------------------------------
for n=2:100/dt
fe=-m*acc(n)+[b1*u(n-1)+b2*v(n-1)-b3*a(n-1)]*m+[b4*u(n-1)+b5*v(n-1)-b6*a(n-1)]*c; %等效荷载
u(n)=fe/ke;
a(n)=b1*[u(n)-u(n-1)]-b2*v(n-1)+b3*a(n-1);
v(n)=b4*[u(n)-u(n-1)]-b5*v(n-1)+b6*a(n-1);
end
% 绘制结构在地震作用下的位移、速度、加速度时程曲线-----------------------------
subplot(3,1,1)
t=(0:length(a)-1)*dt;
plot(t,a) %加速度时程曲线
Acc_Max=max(abs(a))
title('Earthquake Response Curve Of Station 51WCW-90','fontsize',15)
ylabel('Acc(cm/s^2)','fontsize',12)
subplot(3,1,2)
plot(t,v) %速度时程曲线
Vel_Max=max(abs(v))
ylabel('Vel(cm/s)','fontsize',12)
subplot(3,1,3)
plot(t,u) %位移时程曲线
Dis_Max=max(abs(u))
xlabel('Time/s','fontsize',12)
ylabel('Dis/cm','fontsize',12)
% End---程序结束-------------
小弟初次发贴,恳请达人们帮分析一下,不胜感激!
|