【Matlab】结构在傅里叶展开下的周期荷载响应——文末附源码
生活随笔
收集整理的這篇文章主要介紹了
【Matlab】结构在傅里叶展开下的周期荷载响应——文末附源码
小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.
一、背景
研究振動體系對于周期荷載的結(jié)構(gòu)響應(yīng)狀態(tài),可借助傅里葉變換,將周期荷載 轉(zhuǎn)換為三角波的線性疊加后,運(yùn)用基本的結(jié)構(gòu)動力學(xué)知識即可解出體系位移情況。 本文重點(diǎn)內(nèi)容:- 分段函數(shù)的傅里葉展開
- 計算阻尼對振動的影響
二、理論推導(dǎo)
一般諧振荷載下的振動方程:?對周期荷載進(jìn)行傅里葉展開:
其中的傅里葉系數(shù):1.無阻尼情況
解出的通解:?
其中2.考慮阻尼情況
解出的通解: 其中?
?
?三、代碼實(shí)現(xiàn)
3.1物理模型
3.2實(shí)現(xiàn)分段函數(shù)的積分
? ? ? ? 此處默認(rèn)外力的分段函數(shù)為兩段,integrals文件的完整代碼
function [ output ] = integrals( ... % 此處默認(rèn)外力的分段函數(shù)為兩段t, ... % 積分變量fx_1, ... % 第一段函數(shù)表達(dá)式leftInterval_1, ... % 左區(qū)間rightInterval_1, ... % 右區(qū)間fx_2, ... % 第二段函數(shù)表達(dá)式leftInterval_2, ... % 左區(qū)間rightInterval_2 ... % 右區(qū)間) output=int(fx_1,t,leftInterval_1,rightInterval_1) + ... % 計算第一段函數(shù)的積分int(fx_2,t,leftInterval_2, rightInterval_2); % 計算第二段函數(shù)的積分 end?3.3傅里葉展開外力函數(shù)、無阻尼時的位移函數(shù)、有阻尼時的位移函數(shù)
????????method_fourier_unfolds文件的完整代碼
function [ ...external_force, ... % 經(jīng)傅里葉展開后的外力函數(shù)displacement_nodump, ... % 無阻尼時的位移函數(shù)displacement_dump ... % 存在阻尼時的位移函數(shù)] = method_fourier_unfolds( ...t1, ... % 前半周期時長t2, ... % 后半周期時長exact_external_force_1, ...exact_external_force_2, ...expand_num, ... % 傅里葉展開項(xiàng)數(shù)omega, ... % 體系固有頻率k, ... % 結(jié)構(gòu)剛度系數(shù)kesai ... % 阻尼比) T = t1 + t2; % 荷載總周期 theta=2*pi/T; % 外荷載頻率% 分段函數(shù)進(jìn)行傅里葉展開的核心代碼 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;% 以傳入的傅里葉展開項(xiàng)數(shù)進(jìn)行計算,最后將結(jié)果疊加即是傅里葉展開的近似 for n=1:expand_num% 以下運(yùn)算沒有復(fù)雜的邏輯,也未涉及復(fù)雜的循環(huán)迭代,只需要按數(shù)學(xué)表達(dá)式寫出即可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的表達(dá)式 displacement_nodump = (displacement_nodump+double(a0))/k; % 返回?zé)o阻尼時的位移表達(dá)式 displacement_dump = (displacement_dump+double(a0))/k; % 返回有阻尼時的位移表達(dá)式 end?3.4main文件中最終實(shí)現(xiàn)(完整代碼)
clear all%% 設(shè)置參數(shù) m = 2e4; % 小球質(zhì)量 EI = 1e6; % 剛度 L = 5; % 桿長 k = 3*EI/(L^3); % 計算體系剛度系數(shù) 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;% 進(jìn)行傅里葉展開 [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('傅里葉展開的外力函數(shù)圖像');????????結(jié)果
四、補(bǔ)充說明
? ? ? ? 因?yàn)檎駝臃匠讨锌紤]了結(jié)構(gòu)的眾多物理力學(xué)性質(zhì),如質(zhì)量,彈性模量,剛度、固有頻率等,而本文借助matlab實(shí)現(xiàn)時,未對公式做簡化,均采用變量形式進(jìn)行運(yùn)算,所以可以在此源碼基礎(chǔ)上,進(jìn)一步探討不同物理量對結(jié)構(gòu)振動的影響。(以下為阻尼比對結(jié)構(gòu)振動的影響)
如有任何疑問,請在評論區(qū)留言,
源文件github地址:https://github.com/darlingxyz/method_fourier_expansion_in_structural_mechanics
總結(jié)
以上是生活随笔為你收集整理的【Matlab】结构在傅里叶展开下的周期荷载响应——文末附源码的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: VC6 SDK 下载
- 下一篇: 京东商品详情API、通过商品ID获得京东