一、背景

研究振動體系對于周期荷載的結構響應狀態，可借助傅里葉變換，將周期荷載 轉換為三角波的線性疊加后，運用基本的結構動力學知識即可解出體系位移情況。 本文重點內容：

• 分段函數的傅里葉展開
• 計算阻尼對振動的影響

二、理論推導

一般諧振荷載下的振動方程：

?對周期荷載進行傅里葉展開：

其中的傅里葉系數：

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】结构在傅里叶展开下的周期荷载响应——文末附源码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。