clear;%清除内存变量 % 动力特征参数赋值++++++++++++++++++++
m=1;%质量
period=(0:0.01:6);%
% 读入地震波++++++++++++++++++++
xg=xlsread('2.0.xlsx');%读取地震波
t=xg(:,1);%给时间向量提出来
dt=t(2)-t(1);%地震波的时间间隔
% 积分常数+++++++++++++++++++++++++
gama=0.505;beta=0.2525;
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/beta/2)*dt;
for jj=1:length(period)
Tn=period(jj);%周期
omega=2*pi/Tn;%圆频率
k=m*omega^2;%刚度
zeta=0.2;%阻尼比
c=2*zeta*omega*m;%阻尼系数
ke=k+b1*m+b4*c;%等效刚度矩阵
% 循环+++++++++++++++++++++++++
x1=0;%初值x1
v1=0;%初值速度
a1=0;%初值加速度
u=zeros(length(t),1);%用来放位移结果的向量
for jjj=1:length(t)
x2=(-m*xg(jjj,2)+...
m*(b1*x1+b2*v1-b3*a1)+...
c*(b4*x1+b5*v1-b6*a1)...
)/ke;%第二时刻的位移
v2=b4*(x2-x1)-b5*v1+b6*a1;%速度
a2=b1*(x2-x1)-b2*v1+b3*a1;%相对加速度
u(jjj)=x2;
v(jjj)=v2;
a(jjj)=a2;
x1=x2; v1=v2; a1=a2;
end
x_max(jj)=max(abs(u));
v_max(jj)=max(abs(v));
ag_max(jj)=max(abs(a'+xg(:,2)));
iii=max(ag_max)
end
%%
figure(1);%指定图号
plot(t,xg(:,2))%画二维图,加速度时程
xlabel('时间(s)');%横轴标签
ylabel('加速度(m/s^2)');%纵轴标签
grid on;%显示网格
figure(2);%指定图号
plot(t,u)%画二维图,位移时程
xlabel('时间(s)');%横轴标签
ylabel('位移(m)');%纵轴标签
grid on;%显示网格
figure(3);%指定图号
plot(t,v)%画二维图,速度时程
xlabel('时间(s)');%横轴标签
ylabel('速度(m/s)');%纵轴标签
grid on;%显示网格
%%
figure(4)
plot(period,ag_max)%画二维图,加速度反应谱
xlabel('周期(s)');%横轴标签
ylabel('加速度(m/s^2)');%纵轴标签
grid on;%显示网格
figure(5)
plot(period,x_max)%画二维图,位移反应谱
xlabel('周期(s)');%横轴标签
ylabel('位移(m)');%纵轴标签
grid on;%显示网格
figure(6)
plot(period,v_max)%画二维图,速度反应谱
xlabel('周期(s)');%横轴标签
ylabel('速度(m/s)');%纵轴标签
grid on;%显示网格
%%
评论14