杜哈梅积分

杜哈梅积分 (1)任务:1、试编写单自由度杜哈梅数值积分的计算程序。 (2)要求:?可以使用至少2种积分方法; ?可以处理的荷载种类:a.荷载数组;(,)b.脉冲荷载;c.函数荷载。中的两种。 (3)绘出系统响应时程曲线，求出最大响应(位移、速度、加速度)。 (4)算例 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 : 图示质量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)的计算中，它与e???*sin(?D?)及——————————————————————————————————————————————— e???*cos(?D?)的乘积形成新的被积函数。选择需要计算的代码，并将其他荷载函数注释。 simpson.m :辛普森积分函数，输入参数为函数句柄、积分上下限(为标量)，输出为对应的积分值。 trapezoid.m :梯形组合积分函数，输入、输出同辛普森积分函数。(功能与simpson.m相同，所以只用一种即可。) FuncLoad.m :连续荷载处理函数，输出为对应时间向量的位移、速度、加速度向量，输入参数为荷载函数句柄、时间向量、结构刚度、集中质量、阻尼比。调用它，可计算连续荷载作用下结构的响应，如三角函数荷载、阶跃荷载及其他连续函数荷载。 stepFuncTion.m :脉冲荷载处理函数，输出为对应时间向量的位移、速度、加速度向量，输入参数为荷载函数句柄、时间向量、荷载终止时刻、结构刚度、集中质量、阻尼比。调用它，可计算脉冲荷载作用下结构的响应。 注:?对于不同的荷载函数，lFunction.m中已写出，计算不同荷载的响应，只需将要计算的荷载函数取消注释，运行example.m文件。程序中积分数值算法为辛普森积分，速度、加速度的数值算法为边界节点为向前或向后差分。荷载函数作为函数句柄导入FuncTion、stepFuncTion函数中，与其他函数相乘形成新的函数。 1 / 11 程序文件如下: simpson.m ——————————————————————————————————————————————— function 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; end trapezoid.m function fv=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; for n=1:10000; s=h/2*(B(n)+B(n+1)); %s 为小梯形面积 fv=fv+s; 2 / 11 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; 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(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*t(i))/(m*w1); end u0=u(n2); %自由振动 v0=(3*u(n2)-4*u(n2-1)+u(n2-2))/(2*h); for i=(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))); 3 / 11 end v(1)=(u(2)-u(1))/h;%求速度列向量 for i=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;%求加速度列向量 for i=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; 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(t),1); 4 / 11 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); end v(1)=(u(2)-u(1))/h;%求速度列向量 for i=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;%求加速度列向量 for i=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 function pv=lFunction(t) %荷载函数(load function) %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); 5 / 11 %调用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 on xlabel('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(I1(1)),'\leftarrow 位移最大值点 ','FontSize',15); ——————————————————————————————————————————————— disp(['位移最大值出现在: t=',num2str(t_max1),'s']); disp(['位移最大值: u_max=',num2str(u_max),'m']); figure(3) %绘制“速度—时间”曲线并标出速度最大点 plot(t,v,'-.b'),grid on xlabel('t(s)');ylabel('v(m/s)'); title('速度—时间曲线'); [C2,I2]=max(abs(v)); v_max=C2(1);t_max2=I2(1)*0.0001; hold on; 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'),grid off xlabel('t(s)');ylabel('a(m/s^2)'); title('加速度—时间曲线'); ——————————————————————————————————————————————— [C3,I3]=max(abs(a)); a_max=C3(1);t_max3=I3(1)*0.0001; hold on; text(t_max3,a(I3(1)),'\leftarrow 加速度最大值点 ','FontSize',15); disp(['加速度最大值出现在: t=',num2str(t_max3),'s']); disp(['加速度最大值: a_max=',num2str(a_max),'m/s^2']); ?运行结果: 6 / 11 命令行显示:位移最大值出现在:t=0.0427s 位移最大值:u_max=5.3216m 速度最大值出现在:t=0.0207s 速度最大值:v_max=196.2601m/s 加速度最大值出现在:t=0.0003s 加速度最大值:a_max=15605.6988m/s^2 图片窗口: 7 / 11 8 / 11 ?修改程序，将爆破荷载函数载入，其他荷载函数注释，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 图片窗口: 9 / 11 10 / 11 11 / 11 ———————————————————————————————————————————————