《杜哈梅积分》由会员分享，可在线阅读，更多相关《杜哈梅积分（13页珍藏版）》请在金锄头文库上搜索。
1、杜哈梅积分程序（1）任务：1、试编写单自由度杜哈梅数值积分的计算程序。（2）要求：可以使用至少2种积分方法；可以处理的荷载种类：a.荷载数组；（？）b.脉冲荷载；c.函数荷载。中的两 种。（3）绘出系统响应时程曲线，求出最大响应（位移、速度、加速度）。（4）算例分析：图示质量m受爆炸荷载作用，P=100kN, m=6.4kg,结构刚度k=34847.77N/m,阻尼比为5%, 荷载如下图示：777解：所做MATLAB程序有6个文件（真正包含5个有用文件，还有1个为组合梯形积分函数，共 参考，未利用）。文件名：example.m, IFunction.m, FuncLoad.m, simpson
2、.m, trapezoid.m, stepFuncTion.m。各文件所实现的功能：example.m :运行程序，调用FuncLoad函数，传入lFunction函数句柄及时间向量、结构刚 度、集中质量、阻尼比，然后由FuncLoad函数生成的位移、速度、加速度向量，绘制“位 移一时间” “速度一时间” “加速度一时间”曲线，同时找到最大响应值，并做标注。 lFunction.m :荷载的函数，荷载关于时间t的函数，作为句柄导入FuncLoad函数中。在A（t）、 B（t）的计算中，它与ege*sin（T）及ee*cos（T）的乘积形成新的被积函数。选择需要DD计算的代码，并将其他荷载函数注
3、释。simpson.m :辛普森积分函数，输入参数为函数句柄、积分上下限（为标量），输出为对应 的积分值。trapezoid.m :梯形组合积分函数，输入、输出同辛普森积分函数。（功能与simpson.m相同， 所以只用一种即可。）FuncLoad.m :连续荷载处理函数，输出为对应时间向量的位移、速度、加速度向量，输入 参数为荷载函数句柄、时间向量、结构刚度、集中质量、阻尼比。调用它，可计算连续荷载 作用下结构的响应，如三角函数荷载、阶跃荷载及其他连续函数荷载。stepFuncTion.m :脉冲荷载处理函数，输出为对应时间向量的位移、速度、加速度向量，输 入参数为荷载函数句柄、时间向量、荷
4、载终止时刻、结构刚度、集中质量、阻尼比。调用它， 可计算脉冲荷载作用下结构的响应。注：对于不同的荷载函数lFunction.m中已写出，计算不同荷载的响应，只需将要计算 的荷载函数取消注释，运行example.m文件。程序中积分数值算法为辛普森积分，速度、加 速度的数值算法为边界节点为向前或向后差分。荷载函数作为函数句柄导入FuncTion、stepFuncTion函数中，与其他函数相乘形成新的函数。对于位移响应计算，杜哈梅积分的数值法在理论上应该对于任意荷载都成立，但是本程序 在载入跳跃荷载函数后，位移荷载响应图形貌似符合要求，但根据位移响应生成的速度、加 速度响应出现问题(出现无数次突变)
5、，因此可判断位移响应有问题。程序文件如下:simpson.mfunction fv =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;for j=1:10000;s=h/3*(B(2*j-1)+4*B(2*j)+B(2*j+1);fv=fv+s;endtrapezoid.mfunction fv=trapezoid(f,a,b)%组合梯形积分
6、法函数%input%outputf 为被积函数为积分下限为积分上限 为积分值%A为(a,b)的10000等分横坐标向量 %B 为 A 对应的函数值abfvA=a:(b-a)/10000:b;B=f(A);fv=0; h=(b-a)/10000;for n=1:10000;s=h/2*(B(n)+B(n+1);fv=fv+s;%s 为小梯形面积end stepFuncLoadfunction u,v,a=stepFuncLoad（lFunction,t,t1,k,m,Xi）%函数荷载作用下的单自由度体系动力响应位移U、速度V、加速度a%input%output% w=sqrt(k/m);- ff
7、2 为何在函数-t为时间向量（单位：S）- k 为结构刚度（单位： N/m）- m 为集中质量（单位： kg）-Xi为阻尼比E- u 为位移列向量- v 为速度列向量- a 为加速度列向量- t1 为脉冲荷载突变时刻%圆频率 ww1=w*sqrt(1-XS2);%阻尼振动频率 w1,w1=sqrt(1- E A2)n1=nUmel(t);h=(t(n1)-t(1)/n1;function fv1=f1(x)%对应 A(t)的被积函数fV1=lFUnction(x).*exp(Xi*w*x).*cos(w1*x);end function fv2=f2(x)%对应 B(t)的被积函数 fv2=l
8、Function(x).*exp(Xi*w*x).*sin(w1*x);endA1=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);for i=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
9、*t(i)/(m*w1);endu0=u(n2);%自由振动v0=(3*u(n2)-4*u(n2-1)+u(n2-2)/(2*h);for i=(n2+1):n1u(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;% 求速度列向量 for i=2:(n1-1)v(i)=(u(i+1)-u(i-1)/(2*h);endv(n1)=(u(n1)-u(n1-1)/h;a(1)=(v(2)-v(1)/h;%求加速度列向量for i=2:(numel(t)
10、-1)a(i)=(v(i+1)-v(i-1)/(2*h);enda(numel(t)=(v(numel(t)-v(numel(t)-1)/h;endFuncLoad.m function u,v,a=FuncLoad(lFunction,t,k,m,Xi)阻尼振动频率w1,w1=sqrt（1- E八2）%函数荷载作用下的单自由度体系动力响应位移u、速度V、加速度%input- ff2为何在函数%-为时间向量（单位：S）%-为结构刚度（单位：N/m）%-为集中质量（单位：kg）%- X为阻尼比E%output- u为位移列向量%-为速度列向量%-为加速度列向量w=sqrt(k/m);圆%频率 w
11、w1=w*sqrt(1-Xi入2);n1=numel(t);h=(t(n1)-t(1)/n1;function fv1=f1(x)对应A(t)的被积函数fv1=lFunction(x).*exp(Xi*w*x).*cos(w1*x);end function fv2=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(
12、t),1);for i=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*w 1);endv(1) = (u(2)-u(1)/h;% 求速度列向量for i=2:(n1-1) v(i)=(u(i+1)-u(i-1)/(2*h);endv(n1)=(u(n1)-u(n1-1)/h;a(1) = (v(2)-v(1)/h;%求加速度列向量for i=2:(n1-1)
13、a(i)=(v(i+1)-v(i-1)/(2*h);enda(n1)=(v(n1)-v(n1-1)/h; endlFunction.mfunction pv=lFunction(t)%荷载函数(load function)%pv=(-100*t/0.08+100)*1000.*(t=0.08);%爆破(脉冲)荷载作用 pv=100*1000.*(t0&t=0.08); %脉冲荷载(恒载，t1=0.08s) %pv=100*1000.*(t0);%恒载作用%pv=100*1000*sin(73.79*t).*(t0);%简谐荷载作用(共振) example.m%运行示例t=eps:0.0001:
14、0.5;u,v,a=stepFuncLoad(lFunction,t,0.08,34847.77,6.4,0.05);%调用 Function 函数，传入荷载函数句柄及其他参数，得到位移、速度、加速度向量plot(t,lFunction(t)/1000);grid on %绘制“荷载时间”曲线 axis(0 0.5 -120 120);xlabel(t(s);ylabel(P(kN);title(荷载一时间曲线);figure(2)%绘制“位移时间”曲线并标出位移最大点plot(t,u,-r),grid onxlabel(t(s);ylabel(u(m);title(位移一时间曲线);C1,I1=max(abs(u); u_max=C1(1);t_max1=I1(1)*0.0001;hold on;text(t_max1,u(l1(1),leftarrow 位移最大值点,FontSize,15); disp(位移最大值出现在：t=,num2str(t_max1),s);disp(位移最大值：u_max=,num2str(u_max),m);figure(3)%绘制“速度一时间”曲线并标出速度最大点plot(t,v,-.b),grid onxlabel(t(s);ylabel(v(m/s)
