【Matlab】结构在傅里叶展开下的周期荷载响应——文末附源码

一、背景

研究振动体系对于周期荷载的结构响应状态，可借助傅里叶变换，将周期荷载 转换为三角波的线性叠加后，运用基本的结构动力学知识即可解出体系位移情况。 本文重点内容：

• 分段函数的傅里叶展开
• 计算阻尼对振动的影响

二、理论推导

一般谐振荷载下的振动方程：

 对周期荷载进行傅里叶展开：

其中的傅里叶系数：

1.无阻尼情况

解出的通解：

 

其中

2.考虑阻尼情况

解出的通解： 其中

 

 

 三、代码实现

3.1物理模型

3.2实现分段函数的积分

        此处默认外力的分段函数为两段，integrals文件的完整代码

function [ output ] = integrals( ... % 此处默认外力的分段函数为两段t, ... % 积分变量fx_1, ... % 第一段函数表达式leftInterval_1, ... % 左区间rightInterval_1, ... % 右区间fx_2, ... % 第二段函数表达式leftInterval_2, ... % 左区间rightInterval_2 ... % 右区间) output=int(fx_1,t,leftInterval_1,rightInterval_1) + ... % 计算第一段函数的积分int(fx_2,t,leftInterval_2, rightInterval_2); % 计算第二段函数的积分 end

 3.3傅里叶展开外力函数、无阻尼时的位移函数、有阻尼时的位移函数

        method_fourier_unfolds文件的完整代码

function [ ...external_force, ... % 经傅里叶展开后的外力函数displacement_nodump, ... % 无阻尼时的位移函数displacement_dump ... % 存在阻尼时的位移函数] = method_fourier_unfolds( ...t1, ... % 前半周期时长t2, ... % 后半周期时长exact_external_force_1, ...exact_external_force_2, ...expand_num, ... % 傅里叶展开项数omega, ... % 体系固有频率k, ... % 结构刚度系数kesai ... % 阻尼比) T = t1 + t2; % 荷载总周期 theta=2*pi/T; % 外荷载频率% 分段函数进行傅里叶展开的核心代码 syms t n; a0=1/T * integrals('t',exact_external_force_1,0,t1,exact_external_force_2,t1,T); an=2/T * integrals('t',exact_external_force_1*cos(n*theta*t),0,t1,exact_external_force_2*cos(n*theta*t),t1,T); bn=2/T * integrals('t',exact_external_force_1*sin(n*theta*t),0,t1,exact_external_force_2*sin(n*theta*t),t1,T);% 初始化必要的值 displacement_nodump = 0; displacement_dump = 0; beta = 1:expand_num; miu = 1:expand_num; external_force = 0;% 以传入的傅里叶展开项数进行计算，最后将结果叠加即是傅里叶展开的近似 for n=1:expand_num% 以下运算没有复杂的逻辑，也未涉及复杂的循环迭代，只需要按数学表达式写出即可beta(n) = n*theta/omega;miu(n) = abs(1./(1-beta(n).^2));epsilon = atan(2*kesai*beta(n)/(1-beta(n).^2));miu_dump = 1./((1-beta(n).^2).^2+(2*kesai*beta(n)).^2).^0.5;displacement_nodump=displacement_nodump + miu(n) .* vpa(eval(an*cos(n*theta.*t)+bn*sin(n*theta.*t)));if (nargin > 7 ) % 若输入值大于7，说明需要同时计算有阻尼和无阻尼情况，否则只计算无阻尼情况displacement_dump = displacement_dump + miu_dump .* vpa(eval(an*cos(n*theta.*t-epsilon)+bn*sin(n*theta.*t-epsilon)));endexternal_force = external_force+vpa(eval(an)*cos(n*theta.*t)+eval(bn)*sin(n*theta.*t)); endexternal_force = external_force+a0; % 返回作用力P的表达式 displacement_nodump = (displacement_nodump+double(a0))/k; % 返回无阻尼时的位移表达式 displacement_dump = (displacement_dump+double(a0))/k; % 返回有阻尼时的位移表达式 end

 3.4main文件中最终实现(完整代码)

clear all%% 设置参数 m = 2e4; % 小球质量 EI = 1e6; % 刚度 L = 5; % 杆长 k = 3*EI/(L^3); % 计算体系刚度系数 w = (k/m)^0.5; % 计算体系固有频率 p_max = 1000; % 最大外力值t1 = 1; % 前半周期 t2 = 1; % 后半周期% 定义符合变量t syms t % 方波荷载 force_1 = 0*t+p_max; force_2 = 0*t;% 进行傅里叶展开 [external_force, displacement_nodump, displacement_dump] = ...method_fourier_unfolds(t1, t2, force_1, force_2, 50, w, k, 0.6);len = 200; y1 = 1:len; y2 = 1:len; p = 1:len; time = 1:len; t = 0; for i=1:lent = t + 0.05;time(i) = t;y1(i) = eval(displacement_nodump);y2(i) = eval(displacement_dump);p(i) = eval(external_force); endsubplot(2, 2, 1); plot(time, y1); title('无阻尼时位移与时间曲线');subplot(2, 2, 2); plot(time, y2); title('有阻尼时位移与时间曲线');subplot(2, 2, [3,4]); plot(time, p); title('傅里叶展开的外力函数图像');

        结果

---

四、补充说明

        因为振动方程中考虑了结构的众多物理力学性质，如质量，弹性模量，刚度、固有频率等，而本文借助matlab实现时，未对公式做简化，均采用变量形式进行运算，所以可以在此源码基础上，进一步探讨不同物理量对结构振动的影响。(以下为阻尼比对结构振动的影响)

---

如有任何疑问，请在评论区留言，

源文件github地址：https://github.com/darlingxyz/method_fourier_expansion_in_structural_mechanics

总结

以上是生活随笔为你收集整理的【Matlab】结构在傅里叶展开下的周期荷载响应——文末附源码的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。