杜哈梅积分.docx
《杜哈梅积分.docx》由会员分享，可在线阅读，更多相关《杜哈梅积分.docx（14页珍藏版）》请在冰豆网上搜索。
杜哈梅积分
杜哈梅积分程序
（1）任务：
1、试编写单自由度杜哈梅数值积分的计算程序。
（2）要求：
①可以使用至少2种积分方法；
②可以处理的荷载种类：
a.荷载数组；（？
）b.脉冲荷载；c.函数荷载。
中的两种。
（3）绘出系统响应时程曲线，求出最大响应（位移、速度、加速度）。
（4）算例分析：
图示质量m受爆炸荷载作用，P=100kN，m=6.4kg，结构刚度k=34847.77N/m，阻尼比为5%，荷载如下图示：
解：
所做MATLAB程序有6个文件（真正包含5个有用文件，还有1个为组合梯形积分函数，共参考，未利用）。
文件名：
example.m，lFunction.m，FuncLoad.m，simpson.m，trapezoid.m，stepFuncTion.m。
各文件所实现的功能：
example.m：
运行程序，调用FuncLoad函数，传入lFunction函数句柄及时间向量、结构刚度、集中质量、阻尼比，然后由FuncLoad函数生成的位移、速度、加速度向量，绘制“位移—时间”、“速度—时间”、“加速度—时间”曲线，同时找到最大响应值，并做标注。
lFunction.m：
荷载的函数，荷载关于时间t的函数，作为句柄导入FuncLoad函数中。
在A（t）、B（t）的计算中，它与
及
的乘积形成新的被积函数。
选择需要计算的代码，并将其他荷载函数注释。
simpson.m：
辛普森积分函数，输入参数为函数句柄、积分上下限（为标量），输出为对应的积分值。
trapezoid.m：
梯形组合积分函数，输入、输出同辛普森积分函数。
（功能与simpson.m相同，所以只用一种即可。
）
FuncLoad.m：
连续荷载处理函数，输出为对应时间向量的位移、速度、加速度向量，输入参数为荷载函数句柄、时间向量、结构刚度、集中质量、阻尼比。
调用它，可计算连续荷载作用下结构的响应，如三角函数荷载、阶跃荷载及其他连续函数荷载。
stepFuncTion.m：
脉冲荷载处理函数，输出为对应时间向量的位移、速度、加速度向量，输入参数为荷载函数句柄、时间向量、荷载终止时刻、结构刚度、集中质量、阻尼比。
调用它，可计算脉冲荷载作用下结构的响应。
注：
①对于不同的荷载函数，lFunction.m中已写出，计算不同荷载的响应，只需将要计算的荷载函数取消注释，运行example.m文件。
程序中积分数值算法为辛普森积分，速度、加速度的数值算法为边界节点为向前或向后差分。
荷载函数作为函数句柄导入FuncTion、stepFuncTion函数中，与其他函数相乘形成新的函数。
②对于位移响应计算，杜哈梅积分的数值法在理论上应该对于任意荷载都成立，但是本程序在载入跳跃荷载函数后，位移荷载响应图形貌似符合要求，但根据位移响应生成的速度、加速度响应出现问题（出现无数次突变），因此可判断位移响应有问题。
程序文件如下:
simpson.m
functionfv=simpson（f,a,b）
%辛普森积分法函数
%input-@f为被积函数的句柄
%-a为积分下限
%-b为积分上限
%output-fv为积分值
A=a:
（b-a）/20000:
b;%A为（a,b）的20000等分横坐标向量
B=f（A）;%B为A对应的函数值向量
fv=0;
h=（b-a）/20000;
forj=1:
10000;
s=h/3*（B（2*j-1）+4*B（2*j）+B（2*j+1））;
fv=fv+s;
end
trapezoid.m
functionfv=trapezoid（f,a,b）
%组合梯形积分法函数
%input-f为被积函数
%-a为积分下限
%-b为积分上限
%output-fv为积分值
A=a:
（b-a）/10000:
b;%A为（a,b）的10000等分横坐标向量
B=f（A）;%B为A对应的函数值
fv=0;
h=（b-a）/10000;
forn=1:
10000;
s=h/2*（B（n）+B（n+1））;%s为小梯形面积
fv=fv+s;
end
stepFuncLoad
function[u,v,a]=stepFuncLoad（lFunction,t,t1,k,m,Xi）
%函数荷载作用下的单自由度体系动力响应——位移u、速度v、加速度a
%input-ff2为何在函数
%-t为时间向量（单位：
s）
%-k为结构刚度（单位：
N/m）
%-m为集中质量（单位：
kg）
%-Xi为阻尼比ξ
%output-u为位移列向量
%-v为速度列向量
%-a为加速度列向量
%-t1为脉冲荷载突变时刻
w=sqrt（k/m）;%圆频率w
w1=w*sqrt（1-Xi^2）;%阻尼振动频率w1,w1=sqrt（1-ξ^2）
n1=numel（t）;
h=（t（n1）-t
（1））/n1;
functionfv1=f1（x）%对应A（t）的被积函数
fv1=lFunction（x）.*exp（Xi*w*x）.*cos（w1*x）;end
functionfv2=f2（x）%对应B（t）的被积函数
fv2=lFunction（x）.*exp（Xi*w*x）.*sin（w1*x）;end
A1=zeros（length（t）,1）;
B1=zeros（length（t）,1）;
u=zeros（length（t）,1）;
v=zeros（length（t）,1）;
a=zeros（length（t）,1）;
%求位移列向量
n2=int64（t1/h+1）;
fori=1:
n2%强迫振动
A1（i）=simpson（@f1,0,t（i））;%对应A（t）
B1（i）=simpson（@f2,0,t（i））;%对应B（t）
u（i）=（A1（i）*sin（w1*t（i））-B1（i）*cos（w1*t（i）））*exp（-Xi*w*t（i））/（m*w1）;
end
u0=u（n2）;%自由振动
v0=（3*u（n2）-4*u（n2-1）+u（n2-2））/（2*h）;
fori=（n2+1）:
n1
u（i）=exp（-Xi*w*t（i-n2+1））*（u0*cos（w1*t（i-n2+1））+（v0+Xi*w*u0）/w1*sin（w1*t（i-n2+1）））;
end
v
（1）=（u
（2）-u
（1））/h;%求速度列向量
fori=2:
（n1-1）
v（i）=（u（i+1）-u（i-1））/（2*h）;
end
v（n1）=（u（n1）-u（n1-1））/h;
a
（1）=（v
（2）-v
（1））/h;%求加速度列向量
fori=2:
（numel（t）-1）
a（i）=（v（i+1）-v（i-1））/（2*h）;
end
a（numel（t））=（v（numel（t））-v（numel（t）-1））/h;
end
FuncLoad.m
function[u,v,a]=FuncLoad（lFunction,t,k,m,Xi）
%函数荷载作用下的单自由度体系动力响应——位移u、速度v、加速度a
%input-ff2为何在函数
%-t为时间向量（单位：
s）
%-k为结构刚度（单位：
N/m）
%-m为集中质量（单位：
kg）
%-Xi为阻尼比ξ
%output-u为位移列向量
%-v为速度列向量
%-a为加速度列向量
w=sqrt（k/m）;%圆频率w
w1=w*sqrt（1-Xi^2）;%阻尼振动频率w1,w1=sqrt（1-ξ^2）
n1=numel（t）;
h=（t（n1）-t
（1））/n1;
functionfv1=f1（x）%对应A（t）的被积函数
fv1=lFunction（x）.*exp（Xi*w*x）.*cos（w1*x）;end
functionfv2=f2（x）%对应B（t）的被积函数
fv2=lFunction（x）.*exp（Xi*w*x）.*sin（w1*x）;end
A1=zeros（length（t）,1）;
B1=zeros（length（t）,1）;
u=zeros（length（t）,1）;
v=zeros（length（t）,1）;
a=zeros（length（t）,1）;
fori=1:
numel（t）%求位移列向量
A1（i）=simpson（@f1,0,t（i））;%对应A（t）
B1（i）=simpson（@f2,0,t（i））;%对应B（t）
u（i）=（A1（i）*sin（w1*t（i））-B1（i）*cos（w1*t（i）））*exp（-Xi*w*t（i））/（m*w1）;
end
v
（1）=（u
（2）-u
（1））/h;%求速度列向量
fori=2:
（n1-1）
v（i）=（u（i+1）-u（i-1））/（2*h）;
end
v（n1）=（u（n1）-u（n1-1））/h;
a
（1）=（v
（2）-v
（1））/h;%求加速度列向量
fori=2:
（n1-1）
a（i）=（v（i+1）-v（i-1））/（2*h）;
end
a（n1）=（v（n1）-v（n1-1））/h;
end
lFunction.m
functionpv=lFunction（t）
%荷载函数（loadfunction）
%pv=（-100*t/0.08+100）*1000.*（t<0.08）+0.*（t>=0.08）;%爆破（脉冲）荷载作用
pv=100*1000.*（t>0&t<0.08）+0.*（t>=0.08）;%脉冲荷载（恒载，t1=0.08s）
%pv=100*1000.*（t>0）;%恒载作用
%pv=100*1000*sin（73.79*t）.*（t>0）;%简谐荷载作用（共振）
example.m
%运行示例
t=eps:
0.0001:
0.5;
[u,v,a]=stepFuncLoad（@lFunction,t,0.08,34847.77,6.4,0.05）;
%调用Function函数，传入荷载函数句柄及其他参数，得到位移、速度、加速度向量
plot（t,lFunction（t）/1000）;gridon%绘制“荷载—时间”曲线
axis（[00.5-120120]）;
xlabel（'t（s）'）;ylabel（'P（kN）'）;
title（'荷载—时间曲线'）;
figure
（2）%绘制“位移—时间”曲线并标出位移最大点
plot（t,u,'-r'）,gridon
xlabel（'t（s）'）;ylabel（'u（m）'）;
title（'位移—时间曲线'）;
[C1,I1]=max（abs（u））;
u_max=C1
（1）;t_max1=I1
（1）*0.0001;
holdon;
text（t_max1,u（I1
（1））,'\leftarrow位移最大值点','FontSize',15）;
disp（['位移最大值出现在：
t=',num2str（t_max1）,'s']）;
disp（['位移最大值：
u_max=',num2str（u_max）,'m']）;
figure（3）%绘制“速度—时间”曲线并标出速度最大点
plot（t,v,'-.b'）,gridon
xlabel（'t（s）'）;ylabel（'v（m/s）'）;
title（'速度—时间曲线'）;
[C2,I2]=max（abs（v））;
v_max=C2
（1）;t_max2=I2
（1）*0.0001;
holdon;
text（t_max2,v（I2
（1））,'\leftarrow速度最大值点','FontSize',15）;
disp（['速度最大值出现在：
t=',num2str（t_max2）,'s']）;
disp（['速度最大值：
v_max=',num2str（v_max）,'m/s']）;
figure（4）%绘制“加速度—时间”曲线并标出加速度最大点
plot（t,a,':
k'）,gridoff
xlabel（'t（s）'）;ylabel（'a（m/s^2）'）;
title（'加速度—时间曲线'）;
[C3,I3]=max（abs（a））;
a_max=C3
（1）;t_max3=I3
（1）*0.0001;
holdon;
text（t_max3,a（I3
（1））,'\leftarrow加速度最大值点','FontSize',15）;
disp（['加速度最大值出现在：
t=',num2str（t_max3）,'s']）;
disp（['加速度最大值：
a_max=',num2str（a_max）,'m/s^2']）;
①运行结果：
命令行显示：
位移最大值出现在：
t=0.0427s
位移最大值：
u_max=5.3216m
速度最大值出现在：
t=0.0207s
速度最大值：
v_max=196.2601m/s
加速度最大值出现在：
t=0.0003s
加速度最大值：
a_max=15605.6988m/s^2
图片窗口：
②修改程序，将爆破荷载函数载入，其他荷载函数注释，example.m选择调用FuncTion函数，同时调整参数输入，得运行结果：
命令行：
位移最大值出现在：
t=0.0379s
位移最大值：
u_max=4.0449m
速度最大值出现在：
t=0.0611s
速度最大值：
v_max=208.7957m/s
加速度最大值出现在：
t=0.0003s
加速度最大值：
a_max=15567.1761m/s^2
图片窗口：
升级会员 