杜哈梅积分Duhamel's integral Matlab编程案例讲解

在马来西亚理工大学UTM的结构抗风和抗震设计课程所布置的课后作业，如图：

对此，编程的源代码为

（建议仅采用部分，或仅作参考，因为未进行验算，可能会出错）：

%定义有具体数值的参数：构件的参数：截面高度为1m，宽度为0.3m，长度为10m，弹性模量E=20GPa；引擎：重量10kN；阻尼比zeta，即ζ：0.01；

h=1;b=0.3;L=10;M=10;zeta=0.01;E=20.*10^9;

%利用上述参数进行计算，得出数值：截面模量I，构件动态刚度K，引擎质量m； m=M.*10^3/9.8; I=b.*h^3/12; K=48.*E.*I./L^3; %利用上述参数进行计算，得出数值：结构的固有角速度wd（单位：1/s），自振频率fn（单位：Hz） %结构动力学原理可知：外力振动频率与结构自振频率接近时，结构的位移最大。 %根据维基百科上的定义，严格来说，wd=wn(1-zata^2)^0.5，wn才是固有角速度，但此处zeta较小，近似取wd=wn； wd=(K./m)^0.5; fn=wd/(2.*pi); %定义时间t，即杜哈梅积分表达式中的积分上限t；因为本结果要用数值解替代参数解，所以需要用t迭代。 %此处，为t建立一个矩阵，设t从1开始，以步长0.1的间隔，进行迭代 t=1:0.1:20; %建立一个零矩阵，用于输出数据：后续求出了结果之后，可以将矩阵中的元素替换。为画图方便，其大小与t矩阵保持一致； yp = zeros(size(t)); %{ 不太清楚原理，是从参考的豆丁网的代码上抄的；没有这个表达式的话，图像会有锯齿 是：引擎的动力方程，P=P0sin(f*tau)；定义函数F代表P，P其变量为杜哈梅积分中的τ（即tau），f为外力震动频率 而在f与tau的乘积小于等于pi的时候，P有取值；否则，取P=0； 在本案例中，P0为结构的自重M，即10kN；f则在0.1Hz到1000Hz之间变化。本代码试取f=24Hz以示范。 %} f=24; F = @(tau) M.*10^3.*sin(f.*tau).*(f*tau<=pi)+0.*(f.*tau>pi); % 第一个循环，求位移：遍历t，t在所有值上，杜哈梅积分的解： for i = 1:length(t) %{ 杜哈梅积分Duhamel Integral的方程式示意： yp(t) = (1/(m*wd)) * integral(F(tau) * exp(-zeta*wd*(t-tau)) * sin(wd*sqrt(1-zeta^2)*(t-tau)), tau=0 to t) 在此基础上设计此函数。此函数的含义为：定义积分函数integrand，在t=i的时候，该函数的表达式为： %} integrand = @(tau) (1/(m*fn)) *F(tau) .* exp(-zeta*fn*(t(i)-tau)) .* sin(fn*sqrt(1-zeta^2)*(t(i)-tau)); % 将输出矩阵yp的第i个元素替换为积分值： yp(i) = integral(integrand, 0, t(i)); end %此时，得到的yp矩阵，即在频率f下，杜哈梅积分的t时间内的所有值的矩阵。 %本程序不进行封装了，方便校验参数的值 %画图： %与此同时，结构自振时最大位移的值Δ为： delta=M.*10^3.*L^3/(3*E*I); %w为引擎的外部动力的范围，题目给出的范围是0.1Hz到1000Hz，但是100Hz到1000Hz期间图像变化较小，无讨论价值；设每次变化幅度为0.1Hz； w=0.1:0.1:100; %同上，设置0矩阵以储存运算结果，D储存位移最大值结果，V储存剪力最大值结果,Stress储存截面压强的最大值的结果 D=zeros(size(w)); V=zeros(size(D)); Stress=zeros(size(D)); %设置beta即β值，其为引擎施加的外部动力震动频率与结构自振频率之比； beta=w./f; %循环历遍，计算最大位移值Δd=Δ for i = 1:length(w) %位移最大值公式：Δd=Δ*1/(((1-β^2)^2+(2*ζ*β)^2)^0.5) D(i)=delta.*1/((1-beta(i)^2)^2+(2*zeta*beta(i))^2)^0.5; %剪力最大值公式：V=kD V(i)=K.*D(i); %截面所受的压强 Stress(i)=V(i)./(h*b); end %查看位移最大值，单位mm DMax=max(D) %查看最大位移状态下，剪力的最大值，单位kN VMax=max(V) %最大剪力作用下，截面所受压强最大值； StressMax=max(V)/(h*b) %判断构件是否会破坏：当最大剪力除以截面面积大于10MPa时，损坏。 if StressMax<=10^7 fprintf('suitable') else fprintf('Structure will destory') end figure(1); %f=24Hz时杜哈梅方程解，位移与时间的关系 subplot(2,2,1) plot(t,yp) xlabel('t') title('Engine Natural frequency f=24Hz') ylabel('D') %引擎的自振频率f变化时，最大位移值的变化 subplot(2,2,2) plot(w,D); xlabel('Engine Natural frequency') title('Max Deflection') ylabel('D') %引擎的自振频率f变化时，最大剪力的变化 subplot(2,2,3) plot(w,V); xlabel('Engine Natural frequency') title('Max Shear Force') ylabel('V') %引擎的自振频率f变化时，截面所受最大压强的变化 subplot(2,2,4) plot(w,Stress); yline(10^7,'--r'); legend('Stress','allowable tension') xlabel('Engine Natural frequency') title('Max Stress') ylabel('Stress')

最后得到的结果为（可能有错，恳请指正）：

DMax =0.3333

VMax =  8.0000e+06

StressMax =  2.6667e+07

Structure will destory

与该课程相关的课件，我会在学期结束之后，打包好上传云盘，贴于本文末尾评论区处。届时读者可以学习其原理。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删