%%单自由度系统杜哈梅数值积分的实现及考题验证
%输入参数:
% 系统参数:
% 质量、刚度、阻尼比
% 选择输入荷载类型:
% 谐荷载、脉冲荷载、爆炸荷载
% 选择积分方法:
% 梯形积分、辛普森积分
%梯形积分:V = H (S_1 + S_2 ) /2.
%辛普森积分:V = H (S_1 + 4S_0 + S_2) /6.
M=input('质量M(kg):');
K=input('刚度K(N/m):');
n=input('阻尼比η:');
Wn=sqrt(K/M); %无阻尼圆频率
Wd=Wn*sqrt(1-n*n); %有阻尼圆频率
loadType=input('选择荷载种类:(谐荷载:1;爆炸荷载:2) 请输入数字:');
switch loadType
case 1 %谐荷载参数定义
p0=input('简谐荷载幅值(N):');
w=input('简谐荷载圆频率(rad/s):');
t=input('简谐荷载作用时间(s):');
dt=input('时间步长(s):');
N=2*t/dt;
P=zeros(1,N+1); %扩充荷载矩阵
for i=1:t/dt;
P(i)=p0*sin(w*(i-1)*dt);%每步谐荷载大小
end
integration=input('选择数值积分方法:(梯形积分:11;辛普森积分:12),请输入数字:');
switch integration
case 11 %梯形积分方法
v=zeros(1,N);%速度
a=zeros(1,N-1);%加速度
Y=zeros(N+1,N+2);%位移
for i=1:N;
for j=1:i;%对每一荷载步进行计算
Y(i+1,j)=P(j)*(exp((-n)*Wn*((i*dt-(j-1)*dt))))*sin(Wd*(i*dt-(j-1)*dt));
Y(i+1,j+1)=P(j+1)*(exp((-n)*Wn*((i*dt-j*dt))))*sin(Wd*(i*dt-j*dt));
Y(i+1,N+2)=Y(i+1,N+2)+(dt/(2*M*Wd))*((Y(i+1,j)+Y(i+1,j+1)));%取Y中间值并乘以区间长
end
end
y=Y(:,N+2);
for i=1:N; %通过位移求速度
v(i)=(y(i+1)-y(i))/dt;
end
for i=1:N-1; %通过速度求加速度
a(i)=(v(i+1)-v(i))/dt;
end
t1=(0:1:N)*dt;%定义各时间点
subplot(2,2,1);%在左上角作荷载-时间图像
plot(t1,P);
title('简谐荷载图(采用梯形法)');
xlabel('时间t(s)');ylabel('荷载P(N)');
subplot(2,2,2);plot(t1,y);%右上角作位移-时间图像
title('竖向位移-时间图');
xlabel('时间t(s)');ylabel('位移y(m)');
t2=0:dt:(N-1)*dt;
subplot(2,2,3);plot(t2,v);%左下角作速度-时间图像
title('速度-时间图');
xlabel('时间t(s)');ylabel('速度v(m/s)');
t3=0:dt:(N-2)*dt;
subplot(2,2,4);plot(t3,a);%右下角作加速度-时间图像
title('加速度-时间图');
xlabel('时间t(s)');ylabel('加速度a(m/s^2)')
hold off;
case 12 %辛普森积分方法
v=zeros(1,N/2);%速度
a=zeros(1,(N/2)-1);%加速度
Y=zeros((N/2)+1,2*N+2);%位移
for i=1:(N/2);
for j=1:2:2*i-1;
Y(i+1,j)=P(j)*exp(-1*n*Wn*(2*i*dt-(j-1)*dt))*sin(Wd*(2*i*dt-(j-1)*dt));
Y(i+1,j+1)=P(j+1)*exp(-1*n*Wn*(2*i*dt-j*dt))*sin(Wd*(2*i*dt-j*dt));
Y(i+1,j+2)=P(j+2)*exp(-1*n*Wn*(2*i*dt-(j+1)*dt))*sin(Wd*(2*i*dt-(j+1)*dt));
Y(i+1,2*N+2)=Y(i+1,2*N+2)+(dt/(3*M*Wd))*(Y(i+1,j)+4*Y(i+1,j+1)+Y(i+1,j+2));
%辛普森公式:V = H/6*(S_1 + 4S_0 + S_2).注意本例中区间长为2dt.
end
end
y=Y(:,2*N+2);
for i=1:N/2; %通过位移求速度
v(i)=(y(i+1)-y(i))/(2*dt);
end
for i=1:N/2-1; %通过速度求加速度
a(i)=(v(i+1)-v(i))/(2*dt);
end
t1=(0:1:N)*dt;
subplot(2,2,1);plot(t1,P);%左上角作荷载图
title('简谐荷载图(辛普森法)');
xlabel('时间t(s)');ylabel('荷载P(N)');
t2=(0:1:N/2)*2*dt;
subplot(2,2,2);plot(t2,y);
title('竖向位移-时间图');
xlabel('时间t(s)');ylabel('位移y(m)');
t3=(0:1:(N/2)-1)*2*dt;
subplot(2,2,3);plot(t3,v);
title('速度-时间图');
xlabel('时间t(s)');ylabel('速度v(m/s)');
t4=(0:1:(N/2)-2)*2*dt;
subplot(2,2,4);plot(t4,a);
title('加速度-时间图');
xlabel('时间t(s)');ylabel('加速度a(m/s^2)');
hold off;
end
case 2 %爆炸荷载参数定义,设爆炸荷载函数为P=A*t+B
A=input('爆炸荷载函数斜率A(N/s):');B=input('爆炸荷载纵轴截距B(N):');
t=input('爆炸荷载横轴截距t(s):');
dt=input('时间步长(s):');
N=2*t/dt;P=zeros(1,N+1);
for i=1:t/dt;
P(i)=A*(dt*i)+B;%从爆炸荷载函数中提取瞬时荷载
end
integration=input('选择数值积分方法:梯形积分:21;辛普森积分:22),请输入数字:');
switch integration
case 21 %梯形积分法
v=zeros(1,N);%速度
a=zeros(1,N-1);%加速度
Y=zeros(N+1,N+2);%位移
for i=1:N;
for j=1:i;
Y(i+1,j)=P(j)*(exp((-n)*Wn*((i*dt-(j-1)*dt))))*sin(Wd*(i*dt-(j-1)*dt));
Y(i+1,j+1)=P(j+1)*(exp((-n)*Wn*((i*dt-j*dt))))*sin(Wd*(i*dt-j*dt));
Y(i+1,N+2)=Y(i+1,N+2)+(dt/(2*M*Wd))*((Y(i+1,j)+Y(i+1,j+1)));
end
end
y=Y(:,N+2);
for i=1:N;
v(i)=(y(i+1)-y(i))/dt;
end
for i=1:N-1;
a(i)=(v(i+1)-v(i))/dt;
end
t1=(0:1:N)*dt;
subplot(2,2,1);plot(t1,P);
title('爆炸荷载图(梯形法)');
xlabel('时间t(s)');ylabel('荷载P(N)');
subplot(2,2,2);plot(t1,y);
title('竖向位移-时间图');
xlabel('时间t(s)');ylabel('位移y(m)');
t2=0:dt:(N-1)*dt;
subplot(2,2,3);plot(t2,v);
title('速度-时间');
xlabel('时间t(s)');ylabel('速度v(m/s)');
t3=0:dt:(N-2)*dt;
subplot(2,2,4);plot(t3,a);
title('加速度-时间图');
xlabel('时间t(s)');ylabel('加速度a(m/s^2)');
hold off;
case 22 %辛普森积分方法
v=zeros(1,N/2);
a=zeros(1,(N/2)-1);
Y=zeros((N/2)+1,2*N+2);
for i=1:(N/2);
for j=1:2:2*i-1;
Y(i+1,j)=P(j)*exp(-1*n*Wn*(2*i*dt-(j-1)*dt))*sin(Wd*(2*i*dt-(j-1)*dt));
Y(i+1,j+1)=P(j+1)*exp(-1*n*Wn*(2*i*dt-j*dt))*sin(Wd*(2*i*dt-j*dt));
Y(i+1,j+2)=P(j+2)*exp(-1*n*Wn*(2*i*dt-(j+1)*dt))*sin(Wd*(2*i*dt-(j+1)*dt));
Y(i+1,2*N+2)=Y(i+1,2*N+2)+(dt/(3*M*Wd))*(Y(i+1,j)+4*Y(i+1,j+1)+Y(i+1,j+2));
end
end
y=Y(:,2*N+2);
for i=1:N/2;
v(i)=(y(i+1)-y(i))/(2*dt);
end
for i=1:N/2-1;
a(i)=(v(i+1)-v(i))/(2*dt);
end
t1=(0:1:N)*dt;
subplot(2,2,1);plot(t1,P);
title('爆炸荷载图(辛普森法)');
xlabel('时间t(s)');ylabel('荷载P(N)');
t2=(0:1:N/2)*2*dt;
subplot(2,2,2);plot(t2,y);
title('竖向位移-时间图');
xlabel('时间t(s)');ylabel('位移y(m)');
t3=(0:1:(N/2)-1)*2*dt;
subplot(2,2,3);plot(t3,v);
title('速度-时间图');
xlabel('时间t(s)');ylabel('速度v(m/s)');
评论1