|
對于一個(gè)實(shí)際結(jié)構(gòu)，由有限元法離散化處理后，動(dòng)力學(xué)方程可寫為： 振型疊加法按照有限單元法的一般規(guī)則，經(jīng)過邊界條件的約束處理，結(jié)構(gòu)在強(qiáng)迫振動(dòng)時(shí)多自由度體系的運(yùn)動(dòng)平衡方程可以表示為： 其中，M是體系的質(zhì)量矩陣，C是體系的阻尼矩陣，而K則是剛度矩陣，R為外荷載向量。 上述動(dòng)力平衡方程實(shí)質(zhì)上是與加速度有關(guān)的慣性力和與速度有關(guān)的阻尼力及與位移有關(guān)的彈性力在時(shí)刻t與荷載的靜力平衡。 振型疊加法是把多自由度體系的結(jié)構(gòu)的整體振動(dòng)分解為與振型次數(shù)相對應(yīng)的單自由度體系，求得各個(gè)單自由度體系的動(dòng)力響應(yīng)后，再進(jìn)行疊加得出結(jié)構(gòu)整體響應(yīng)。 振型疊加法原理是利用結(jié)構(gòu)無阻尼自由振動(dòng)的振型矩陣作為變換矩陣，將結(jié)構(gòu)動(dòng)力方程式(1)式變換成一組非耦合的微分方程。逐個(gè)地求解這些方程后，將解疊加即可得到動(dòng)力方程的解。 將體系單元節(jié)點(diǎn)的位移向量表示為如下的變換形式： 式中的變換矩陣Φ是由動(dòng)力方程對應(yīng)的無阻尼自由振動(dòng)方程解出的前m階振型矩陣。即Φ=[φ1，φ2，...φm]，x(t)是與時(shí)間有關(guān)的m階向量，x的各分量稱為廣義位移。 將式(2)代入動(dòng)力方程(1)并左乘以Φt ，則可得廣義位移為未知數(shù)的方程： 式中 現(xiàn)在進(jìn)一步考察式(4)，考慮到特征向量的正交性，可得： 于是對應(yīng)于振型的廣義位移的平衡方程(3)可改寫為： 其中，Λ為特征值 將式(2)稍加運(yùn)算可得廣義位移用有限元位移表示的形式： 在(6)式中，當(dāng)忽略了阻尼的影響，平衡方程為互不耦合的，可以對每個(gè)方程逐個(gè)地進(jìn)行時(shí)間積分。出于相同的考慮，在對有阻尼的體系進(jìn)行分析時(shí)仍然希望采用相同的計(jì)算過程去求解互不耦合的平衡方程式。問題是式(6)中的阻尼陣C通常不能象體系的質(zhì)量陣和剛度陣那樣由單元的剛度陣和質(zhì)量陣裝配而成。但當(dāng)假定阻尼與固有頻率成比例嗎，即假定： 式中，ξ1是振型阻尼參數(shù)；δij是Kronecker符號(當(dāng)i=j時(shí)，δij=1。當(dāng)i≠j時(shí)，δij=0)。 這時(shí)式(6)可簡化為如下形式的若干個(gè)方程式：
其中xi(t)的初始條件為下式：
式(10)表示了一個(gè)具有單位質(zhì)量，剛度為ω12的自由度體系當(dāng)阻尼比為ξi時(shí)的運(yùn)動(dòng)平衡控制方程。這個(gè)平衡方程的求解可通過計(jì)算Duhamel積分求得。
式中
當(dāng)利用式(9)來考慮阻尼的影響時(shí)意味著假設(shè)結(jié)構(gòu)的總阻尼是每個(gè)振型的阻尼之和，而每個(gè)振型上的阻尼是能夠量測的，況且在大多數(shù)情況下結(jié)構(gòu)的阻尼比更易于量測。因而便于用來近似地反映結(jié)構(gòu)體系的阻尼特性。 同時(shí)在計(jì)算上也避免計(jì)算阻尼陣而只需計(jì)算剛度陣和質(zhì)量陣。 Duhamel積分遞推公式對以上方程式(10)，考慮某一模態(tài)的振動(dòng)，并略去下標(biāo)i可寫為：
考慮初始條件為：
此時(shí)其定解為：
式中：
將上式對時(shí)間求導(dǎo)，并將t，t+2△(其中2△為時(shí)間步長)帶入t0，最后通過拋物線法則計(jì)算式中的卷積則可得到：
記
則上述A矩陣元素為：
應(yīng)用上述遞推公式，以前一時(shí)刻來求后一時(shí)刻的結(jié)果。計(jì)算不重復(fù)。當(dāng)aij求出后，以后在時(shí)域中的步進(jìn)求解只是一些簡單的數(shù)組相乘。計(jì)算速度很快。 Duhamel積分matlab實(shí)例%%%%%%%%%%%%%%%%% %單自由度Duhamel積分 %內(nèi)部采用Simpson積分 %初始狀態(tài)為靜止?fàn)顟B(tài) %可計(jì)算無阻尼和阻尼強(qiáng)迫震動(dòng) %輸入可為數(shù)組或函數(shù)荷載 %只能得出位移時(shí)程圖 %%%%%%%%%%%%%%%%% clear all; %w=30 w=input('輸入自振頻率 ω:'); %n=10 n=input('輸入荷載步數(shù)n ：'); %m=96.6/32.3 m=input('輸入單元質(zhì)量m ：'); %T=0.05 T=input('輸入外荷載持續(xù)時(shí)間T：'); %si=0.0 si=input('輸入阻尼系數(shù)ξ：'); %k=2700 k=input('輸入剛度系數(shù)k：'); deta=T/n; wD=w*(1-si^2)^0.5; G=deta/(m*wD)/3; t1=[0:deta:T]; reply = input('輸入荷載數(shù)組直接回車或敲擊Y，輸入函數(shù)荷載輸入N： [Y]:','s'); if isempty(reply) reply = 'Y'; end if reply=='Y' p2=input('輸入荷載數(shù)組并用[]抱住為n＋1個(gè)元素：例如[019.3200 38.6400 57.9600 77.2800 96.6000 77.2800 57.9600 38.6400 19.32000]:\n');
%p2=[0:19.32:96.6,77.28:-19.32:0] exp1=exp(-2*si*w*deta); exp2=4*exp(-si*w*deta); A(1)=p2(1)*cos(wD*0); for i=2:n/2+1 A(i)=(A(i-1)+p2(2*i-3)*cos(wD*(2*deta*(i-2))))*exp1+p2(2*i-2)*cos(wD*(2*i*deta-3*deta))*exp2+p2(2*i-1)*cos(wD*2*(i-1)*deta); end
B(1)=p2(1)*sin(wD*0); for i=2:n/2+1 B(i)=(B(i-1)+p2(2*i-3)*sin(wD*(2*deta*(i-2))))*exp1+p2(2*i-2)*sin(wD*(2*i*deta-3*deta))*exp2+p2(2*i-1)*sin(wD*2*(i-1)*deta); end
t1=[0:2*deta:T]; disp('位移隨時(shí)間的變化'); V=G*(A.*sin(wD*t1)-B.*cos(wD*t1)); disp('彈性力隨時(shí)間的變化'); f=k*V; for i=1:n/2+1 p22(i)=p2(2*i-1); end disp('加速度隨時(shí)間的變化'); v11=(p22-k*V)/m; subplot(1,2,1),plot(t1,V); title('位移'),grid on,
t2=[T:deta:T*5]; sw=sin(w*t2); cw=cos(w*t2); disp('擴(kuò)大四倍后位移隨時(shí)間的變化'); y=G*A(n/2+1)*sw-G*B(n/2+1)*cw; subplot(1,2,2),
plot(t1,V,'r') ,hold on; plot(t2,y,'r'); title('加載時(shí)間增加四倍時(shí)的位移位移時(shí)程曲線'); grid on; else
syms tao t; y=input('函數(shù)荷載(關(guān)于未知數(shù)tao)：'); %y=-10000/8*tao+100; v1=1/m/wD*int(y*exp(-si*w*(t-tao))*sin(wD*(t-tao)),tao,0,t); v01=diff(v1,t); v02=diff(v01,t) t1=[0:T/n:T]; for i=1:n+1 v2(i)=subs(v1,t,t1(i)); v011(i)=subs(v01,t,t1(i)); v022(i)=subs(v02,t,t1(i)); end disp('荷載作用時(shí)間內(nèi)位移變化'); v2 disp('荷載作用時(shí)間內(nèi)位移最大值'); max(v2) disp('荷載作用時(shí)間內(nèi)速度變化'); v011 disp('荷載作用時(shí)間內(nèi)速度最大值'); max(v011) disp('荷載作用時(shí)間內(nèi)加速度變化'); v022 disp('荷載作用時(shí)間內(nèi)加速度最大值'); max(v022) subplot(2,3,1); plot(t1,v2),grid on; title('荷載作用時(shí)間內(nèi)位移時(shí)程曲線'); subplot(2,3,2); plot(t1,v011),grid on; title('荷載作用時(shí)間內(nèi)速度時(shí)程曲線'); subplot(2,3,3); plot(t1,v022),grid on; title('荷載作用時(shí)間內(nèi)加速度時(shí)程曲線'); A=1/m/wD*int(y*exp(si*w*tao-si*w*t)*cos(wD*tao),tao,0,T); B=1/m/wD*int(y*exp(si*w*tao-si*w*t)*sin(wD*tao),tao,0,T);
A1=subs(A,T); B1=subs(B,T); V33=A1*sin(wD*t)-B1*cos(wD*t); V033=diff(V33,t); V0033=diff(V033,t);
t2=[T:T/n:5*T]; v33=A1*sin(wD*t2)-B1*cos(wD*t2); for i=1:4*n+1 V03(i)=subs(V033,t,t2(i)); V003(i)=subs(V0033,t,t2(i)); end disp('擴(kuò)大四倍后位移隨時(shí)間的變化'); v33 disp('擴(kuò)大四倍后位移最大值'); max(v33) disp('擴(kuò)大四倍后速度隨時(shí)間的變化'); V03 disp('擴(kuò)大四倍后速度最大值'); max(V03) disp('擴(kuò)大四倍后加速度隨時(shí)間的變化'); V003 disp('擴(kuò)大四倍后加速度最大值'); max(V003) subplot(2,3,4); plot(t1,v2),hold on; title('加載時(shí)間增加四倍時(shí)的位移時(shí)程曲線'); plot(t2,v33),grid on; subplot(2,3,5); plot(t1,v011),hold on; title('加載時(shí)間增加四倍時(shí)的速度時(shí)程曲線'); plot(t2,V03),grid on; subplot(2,3,6); plot(t1,v022),hold on; title('加載時(shí)間增加四倍時(shí)的加速度時(shí)程曲線'); plot(t2,V003),grid on; end |
|
|
