基于頻響函數(shù)矩陣計(jì)算阻尼系統(tǒng)動(dòng)力響應(yīng)的新方法

張 淼， 于 瀾 ， 鞠 偉

(1.長(zhǎng)春工程學(xué)院 理學(xué)院，長(zhǎng)春 130012；2.中國(guó)第一汽車股份有限公司 技術(shù)中心，長(zhǎng)春 130012)

計(jì)算結(jié)構(gòu)動(dòng)力響應(yīng)時(shí)若采用直接積分法，對(duì)于每一個(gè)時(shí)間步長(zhǎng)，其運(yùn)算次數(shù)與半帶寬、自由度數(shù)的乘積成正比。當(dāng)半帶寬較大且時(shí)間歷程遠(yuǎn)大于系統(tǒng)的最小固有振動(dòng)周期時(shí)，結(jié)構(gòu)動(dòng)力響應(yīng)的計(jì)算將是很耗時(shí)的。而振型迭加法在一定條件下可以取得比直接積分法高的計(jì)算效率，它主要是利用系統(tǒng)自由振動(dòng)的模態(tài)振型將多自由度的動(dòng)力方程組轉(zhuǎn)換成為相互解耦的獨(dú)立方程，對(duì)每一個(gè)方程可以求解其響應(yīng)的解析解和數(shù)值解。在對(duì)每個(gè)方程求解時(shí)常常采用Duhamel積分法，在一般情況下，它也需要數(shù)值積分來(lái)計(jì)算，只有極少數(shù)簡(jiǎn)單情況才可以得到解析解。直接積分法和振型迭加法各有優(yōu)勢(shì)，目前被工程界廣泛應(yīng)用。

對(duì)于經(jīng)典阻尼系統(tǒng)，由于其阻尼矩陣能被系統(tǒng)的無(wú)阻尼固有振型對(duì)角化，因此可以實(shí)現(xiàn)系統(tǒng)解耦，在理論上其響應(yīng)求解不存在困難，只是在求解每一個(gè)解耦方程的過(guò)程中，需要考慮所使用的算法的穩(wěn)定性及計(jì)算代價(jià)[1]。若將系統(tǒng)轉(zhuǎn)入狀態(tài)間格式，在狀態(tài)方程中使用直接積分法，迭代產(chǎn)生響應(yīng)的近似解序列[2～4]，與其伴隨的則是計(jì)算精度及計(jì)算效率的平衡問(wèn)題[5]。對(duì)于非經(jīng)典阻尼引起的耦合問(wèn)題，傳統(tǒng)的實(shí)模態(tài)理論無(wú)法使方程解耦。目前常用的計(jì)算響應(yīng)的方法是近似對(duì)角化法，但若結(jié)構(gòu)的整體模態(tài)阻尼矩陣的某些非對(duì)角元數(shù)值較大時(shí)，直接忽略非對(duì)角項(xiàng)所引起的誤差的可控性還有待研究[6]；再者是采用復(fù)模態(tài)解耦方法來(lái)計(jì)算響應(yīng)[7]，但計(jì)算過(guò)程相當(dāng)復(fù)雜。

文獻(xiàn)[8]提出了一種函數(shù)變換方法，轉(zhuǎn)移系統(tǒng)阻尼項(xiàng)的影響，使用數(shù)值技術(shù)求解時(shí)變矩陣的特征向量，給出了自由振動(dòng)的響應(yīng)近似解。文獻(xiàn)[9]中利用矩陣分析原理及變換前后系統(tǒng)的振型向量之間的內(nèi)在聯(lián)系，給出了系統(tǒng)相應(yīng)的模態(tài)計(jì)算方法。由于頻率響應(yīng)矩陣包含的信息豐富，實(shí)測(cè)方便，測(cè)量的精度較高，因此使用頻率響應(yīng)矩陣解決許多工程應(yīng)用的實(shí)際問(wèn)題，諸如模型修正[10]及結(jié)構(gòu)損傷識(shí)別[11-12]等。在這些研究的基礎(chǔ)上，本文擬使用矩陣函數(shù)變換將阻尼系統(tǒng)轉(zhuǎn)化為無(wú)阻尼時(shí)變系統(tǒng)，從而既克服了阻尼項(xiàng)所帶來(lái)的困擾，又更容易求得其模態(tài)參數(shù)，進(jìn)而推導(dǎo)并證明了計(jì)算其頻率響應(yīng)矩陣的公式，并以此為基礎(chǔ)提出了一種基于頻響函數(shù)矩陣求解經(jīng)典和非經(jīng)典阻尼系統(tǒng)精確響應(yīng)的新方法，結(jié)論的公式精確簡(jiǎn)單，易于實(shí)施。本文方法在解決經(jīng)典阻尼問(wèn)題時(shí)，與振型迭加法同為解析解；對(duì)于非經(jīng)典阻尼系統(tǒng)，無(wú)法直接使用振型迭加法解決，而相較于Newmark法，本文方法的最大優(yōu)勢(shì)在于它為精確解，而非數(shù)值解。

1 多自由度阻尼系統(tǒng)的轉(zhuǎn)化

對(duì)自由度為n的阻尼系統(tǒng)，設(shè)M、C、K分別是對(duì)稱的質(zhì)量、阻尼和剛度矩陣，相應(yīng)的運(yùn)動(dòng)微分方程為:

(1)

若取V=[{v1},…,{vn}]為無(wú)阻尼正則振型矩陣，則模態(tài)質(zhì)量和模態(tài)剛度矩陣變?yōu)?

VTMV=diag(m1,m2,…,mn)=I

(2)

記模態(tài)阻尼矩陣為:

(3)

并取坐標(biāo)變換

{x(t)}=V{q(t)}

(4)

在模態(tài)坐標(biāo){q(t)}下，式(1)等價(jià)于:

(5)

當(dāng)式(5)中的模態(tài)阻尼陣為對(duì)角矩陣D=diag(d1,…,dn)時(shí)，式(1)稱為經(jīng)典阻尼系統(tǒng)，否則稱為非經(jīng)典阻尼系統(tǒng)。事實(shí)上，對(duì)經(jīng)典阻尼系統(tǒng)，式(5)表明，原動(dòng)力系統(tǒng)已經(jīng)實(shí)現(xiàn)解耦，正如引言中所敘述的那樣，可以采用振型迭加法，求解出每個(gè)節(jié)點(diǎn)處的動(dòng)力響應(yīng)。為了得到更加精確的結(jié)論，本文提出一種計(jì)算阻尼系統(tǒng)精確響應(yīng)的新方法。為此，定義一個(gè)與模態(tài)阻尼矩陣有關(guān)的矩陣函數(shù)變換:

{y(t)}=eDt{q(t)}

(6)

或其逆變換:

{q(t)}=eDt{y(t)}

(7)

則:

(8)

和

(9)

將式(7)～式(9)代入式(5)得:

(10)

用eDt左乘式(10)兩端，有:

(11)

這時(shí)式(11)變?yōu)闊o(wú)阻尼振動(dòng)系統(tǒng)的運(yùn)動(dòng)微分方程。為了討論其特征問(wèn)題，引入狀態(tài)方程。

2 阻尼系統(tǒng)的復(fù)模態(tài)分析

對(duì)一般動(dòng)力系統(tǒng)式(1)，把時(shí)間域上的矩陣方程變換到以λ為變量的拉氏域中，并假定初始位移和初始速度均為零，則得:

(λ2M+λC+K){X(λ)}={F(λ)}

(12)

令:

Z(λ)=λ2M+λC+K

則:

Z(λ){X(λ)}={F(λ)}

即:

{X(λ)}=H(λ){F(λ)}

式中H(λ)稱為傳遞矩陣。沿頻率軸jω計(jì)算的傳遞矩陣稱為頻率響應(yīng)矩陣。且:

|Z(λ)|=0

(13)

即為式(1)的特征方程，特征方程的根為式(1)的極點(diǎn)，決定式(1)的共振頻率。

為了把式(12)轉(zhuǎn)化為一般特征問(wèn)題，我們引入恒等式：

(λE-λE){X(λ)}=0

(14)

將式(12)與式(14)相結(jié)合得:

(15)

其中:

如果右端向量為零，式(17)就成了關(guān)于實(shí)值矩陣A的一般特征問(wèn)題:

(A-λE){Y}={0}

(16)

其特征值滿足方程:

|A-λE|=0

(17)

(i=1,2,…,n)

基于上面關(guān)于一般動(dòng)力系統(tǒng)的復(fù)模態(tài)理論，具體考慮無(wú)阻尼系統(tǒng)式(11)的極點(diǎn)和模態(tài)振型，并利用它們求系統(tǒng)式(11)的響應(yīng)。

3 經(jīng)典阻尼系統(tǒng)的精確響應(yīng)求解

對(duì)經(jīng)典阻尼系統(tǒng)，系統(tǒng)(1)的模態(tài)阻尼矩陣D為對(duì)角陣，所以矩陣D2也為對(duì)角陣，根據(jù)矩陣代數(shù)理論，矩陣函數(shù)f(D)=diag(f(d1),…,f(dn))也是對(duì)角陣，即

(18)

從而無(wú)阻尼系統(tǒng)式(11)的系數(shù)矩陣便約化為常數(shù)對(duì)角陣,即

eDt(Λ-D2)e-Dt=Λ-D2

相應(yīng)地，在坐標(biāo)變換式(7)下，式(5)等價(jià)于:

(19)

即由式(16)可有:

[(D2-Λ)-λ2E]{X}={0}

(20)

(r=1,2,…,n)

(21)

其中{er}為單位向量。至此求出了系統(tǒng)式(19)的模態(tài)參數(shù)，就可以利用它們來(lái)表示其頻響函數(shù)。

定理1[13]對(duì)阻尼系統(tǒng)(1)，其傳遞矩陣為:

其中Qr為比例換算因子。

定理2[13]對(duì)阻尼系統(tǒng)(1)，其頻率響應(yīng)矩陣為:

其中2jωrmrQr=1，mr為模態(tài)質(zhì)量。

由于系統(tǒng)式(19)的模態(tài)向量為單位坐標(biāo)向量，即VTV=I，且該系統(tǒng)的質(zhì)量陣為單位陣，因此由定理3可得系統(tǒng)式(19)的頻響函數(shù)矩陣為對(duì)角陣:

(22)

由此，在簡(jiǎn)諧激勵(lì)力作用下，當(dāng)輸入力為{f}={F}eiωt時(shí)，系統(tǒng)式(19)的穩(wěn)態(tài)響應(yīng)為:

{y(t)}=H(jω)eDtVT{f}

(23)

再根據(jù)式(4)和式(7)，可得原阻尼系統(tǒng)式(1)的響應(yīng)為:

{x(t)}=Ve-DtH(jω)eDtVT{f}

(24)

根據(jù)式(18)及式(22)可知，eDt,e-Dt及H(ω)均為對(duì)角陣，所以上式簡(jiǎn)化為:

{x(t)}=VH(jω)VT{f}

(25)

由式(25)可知只需利用無(wú)阻尼正則振型、無(wú)阻尼固有頻率及模態(tài)阻尼矩陣等信息即可求得簡(jiǎn)諧激勵(lì)下的任一經(jīng)典阻尼動(dòng)力系統(tǒng)的精確響應(yīng)。

4 非經(jīng)典阻尼系統(tǒng)的精確響應(yīng)求解

對(duì)非經(jīng)典阻尼系統(tǒng)，無(wú)法直接使用振型迭加法求解，為此本文提出一種基于模態(tài)參數(shù)求解非經(jīng)典阻尼精確響應(yīng)的新方法。

對(duì)非經(jīng)典阻尼系統(tǒng)，模態(tài)阻尼矩陣D不是對(duì)角陣，對(duì)應(yīng)系統(tǒng)式(11)的狀態(tài)矩陣A為

與(20)式類似地可有

[eDt(D2-Λ)e-Dt-λ2E]{X}={0}

這是關(guān)于矩陣eDt(D2-Λ)e-Dt的特征問(wèn)題表達(dá)式。模態(tài)阻尼矩陣D雖不是對(duì)角陣，但它為對(duì)稱陣，因此可正交相似于對(duì)角陣，即存在正交陣P，使:

D=Pdiag[k1,…,kn]PT

(26)

根據(jù)矩陣代數(shù)理論，矩陣函數(shù):

eDt=Pdiag[ek1t,…,eknt]PT

(27)

(r=1,2,…,n)

其中{er}為單位向量。由于系統(tǒng)式(11)的模態(tài)向量為eDt{er}(r=1,2,…,n)，且該系統(tǒng)的質(zhì)量陣為單位陣，因此由定理3系統(tǒng)式(11)的頻響函數(shù)矩陣可簡(jiǎn)化為對(duì)角陣：

(28)

由此，在簡(jiǎn)諧激勵(lì)條件下，當(dāng)輸入力為{f}={F}eiωt時(shí)，系統(tǒng)式(11)的穩(wěn)態(tài)響應(yīng)為：

{y(t)}=[H(jω)]eDtVT{f}

(29)

再由式(4)和式(7)可得系統(tǒng)式(1)的穩(wěn)態(tài)響應(yīng)為：

{x(t)}=Ve-Dt[H(jω)]eDtVT{f}

(30)

其中H(jω)的算法見(jiàn)式(28)， eDt的算法見(jiàn)式(27)。

5 數(shù)值算例

作為算例1，考慮如圖所示的動(dòng)力系統(tǒng)。

圖1 兩自由度阻尼振動(dòng)系統(tǒng)

對(duì)上面兩自由度阻尼振動(dòng)系統(tǒng)，令:

k1=k2=k3=k=2 000 N/m

C1=C2=C3=c=3 N/(m·s-1)

M1=M2=m=2 kg

(31)

解得系統(tǒng)的無(wú)阻尼正則振型矩陣為:

(32)

由于模態(tài)阻尼陣:

(33)

為對(duì)角陣，所以該系統(tǒng)為經(jīng)典阻尼系統(tǒng)。由振型迭加法可推得響應(yīng)的解析表達(dá)式為:

(34)

其中:

由本文提出的方法，

(35)

根據(jù)式(25)得:

(36)

表1 本文算法與振型迭加法計(jì)算響應(yīng)解的比較

下面在時(shí)間6 s內(nèi)，分別以時(shí)間步長(zhǎng)Δt=0.25 s和Δt=1 s應(yīng)用本文方法求解響應(yīng)，并與振型迭加法的計(jì)算結(jié)果進(jìn)行了比較，具體結(jié)果見(jiàn)表1。

從表1的結(jié)果看，對(duì)于經(jīng)典阻尼系統(tǒng)的響應(yīng)問(wèn)題，本文方法的計(jì)算結(jié)果與振型迭加法的計(jì)算結(jié)果基本一致，而本文的響應(yīng)表達(dá)式是通過(guò)解析方法得到的，這說(shuō)明本文方法與振型迭加法的表達(dá)式是等價(jià)的。二者的偏差存在于振型迭加法需反復(fù)計(jì)算兩倍于系統(tǒng)自由度數(shù)的反正切和三角正弦值，在計(jì)算機(jī)運(yùn)算過(guò)程中不可避免地產(chǎn)生了誤差的累積。

作為算例2，考慮如下的彈簧質(zhì)量系統(tǒng)，如圖2所示。

圖2 彈簧質(zhì)量系統(tǒng)

Fig.2 Spring-mass system

圖中k1=0.95，k2=k3=…=k10=0.03，k11=1.05，m1=m2=…=m10=1.0，每個(gè)頻率的自然阻尼系數(shù)都為0.05。用有限元方法提取系統(tǒng)的性質(zhì)矩陣后，計(jì)算其無(wú)阻尼正則振型矩陣為:

由式(3)計(jì)算模態(tài)阻尼陣為:

表2 本文算法及Newmark法計(jì)算響應(yīng)解的比較

D并非純對(duì)角矩陣，因此該系統(tǒng)為非經(jīng)典阻尼系統(tǒng)。設(shè)簡(jiǎn)諧激勵(lì)為:

{f}=(F1,F2,F3,F4,F5,F6,F7,F8,F9,F10)Tsinωt

(37)

其中:F1=F3=F5=F7=F9=1,F2=F4=F6=F8=F10=2,ω=2。下面在時(shí)間2 s內(nèi)應(yīng)用本文方法求解響應(yīng)，并與Newmark方法進(jìn)行比較，結(jié)果見(jiàn)表2。

從表2的結(jié)果看，對(duì)于非經(jīng)典阻尼系統(tǒng)的響應(yīng)問(wèn)題，本文方法與Newmark方法求得的響應(yīng)值相近，即驗(yàn)證了本文方法的正確及可行性。二者的偏差存在于本文方法為解析解，而Newmark方法為數(shù)值解。而在計(jì)算效率方面，Newmark方法在每個(gè)時(shí)間步長(zhǎng)上都需要求解一個(gè)未知數(shù)個(gè)數(shù)與自由度數(shù)相同的方程組，為了提高精度還需要縮小時(shí)間步長(zhǎng)，勢(shì)必帶來(lái)計(jì)算量的增長(zhǎng)，因此本文方法在計(jì)算效率及精度方面均優(yōu)于Newmark方法。

6 結(jié) 論

本文提出了一種計(jì)算經(jīng)典及非經(jīng)典阻尼系統(tǒng)響應(yīng)的新方法。只需利用任意動(dòng)力系統(tǒng)所對(duì)應(yīng)的無(wú)阻尼正則振型、無(wú)阻尼固有頻率及模態(tài)阻尼陣，即可求得任意激勵(lì)下經(jīng)典及非經(jīng)典阻尼動(dòng)力系統(tǒng)的精確響應(yīng)。本文提出的算法公式簡(jiǎn)潔、緊湊和精確，不存在計(jì)算誤差，不僅在計(jì)算響應(yīng)時(shí)使效率及精度均得到提高，而且還可作為優(yōu)化目標(biāo)函數(shù)應(yīng)用于模型修正、結(jié)構(gòu)設(shè)計(jì)優(yōu)化及損傷識(shí)別等工程領(lǐng)域，有著良好的前景。

參 考 文 獻(xiàn)

[1]徐醫(yī)培，李素有，吳立言.結(jié)構(gòu)動(dòng)態(tài)響應(yīng)的求解方法分析[J].機(jī)械設(shè)計(jì)與制造，2009,6:12-14.

XU Yi-pei,LI Su-you,WU Li-yan.Analysis about several methods of solving dynamic response of structures[J].Machinery Design & Manufacture,2009,6: 12-14.

[2]張森文，曹開(kāi)彬.計(jì)算結(jié)構(gòu)動(dòng)力響應(yīng)的狀態(tài)方程直接積分法[J].計(jì)算力學(xué)學(xué)報(bào)，2000,17(1): 94-97.

ZHANG Sen-wen,CAO Kai-bin.Direct integration of state equation method for dynamic response of structure [J].Chinese Journal of Computation Mechanics,2000,17(1): 94-97.

[3]劉寒冰，張 淼，魏 健.結(jié)構(gòu)動(dòng)力響應(yīng)分析的多重網(wǎng)格方法[J].吉林大學(xué)學(xué)報(bào)：工學(xué)版，2008，38(3)：619-623.

LIU Han-bing,ZHANG Miao,WEI Jian.Multigrid method for structure dynamic response analysis [J].Journal of Jilin University: Engineering and Technology Edition,2008,38(3): 619-623.

[4]蒲軍平，劉 巖，王元豐，等.基于精細(xì)時(shí)程積分的結(jié)構(gòu)動(dòng)力響應(yīng)降維分析[J].清華大學(xué)學(xué)報(bào)：自然科學(xué)版，2002，42(12)：1681-1683.

PU Jun-ping,LIU Yan,WANG Yuan-feng,et al.Structure dynamic response analysis by reduced dimensions based on a precise time-integration method [J].Journal of Tsinghua University: Science and Technology,2002,42(12): 1681-1683.

[5]郭興旺，鄒家祥.對(duì)機(jī)械振動(dòng)系統(tǒng)的六種動(dòng)態(tài)響應(yīng)分析方法的評(píng)述[J].振動(dòng)與沖擊，1996,15(2): 43-46.

GUO Xing-wang,ZOU Jia-xiang.A review of six dynamic response analysis methods of mechanical vibration systems [J].Journal of Vibration and Shock,1996,15(2): 43-46.

[6]周錫元，馬東輝，俞瑞芳.工程結(jié)構(gòu)中的阻尼與復(fù)振型地震響應(yīng)的完全平方組合[J].土木工程學(xué)報(bào)，2005,38(1): 31-39.

ZHOU Xi-yuan,MA Dong-hui,YU Rui-fang.Damping in structures and complete qudratic combination (CCQC)of complex mode seismic responses[J].China Civil Engineering Journal,2005,38(1): 31-39.

[7]Greco A，Santini A.Comparative study on dynamic analysis of non-classically dymped linear system[J].Structural Engineering and Mechanics,2002,14(6):679-698.

[8]Liu Z S，Song D T，Huang C.Vibration analysis of non-classically damped linear systems [J].Journal of Vibration and Acoustics,2004,126:456-458.

[9]徐 濤，程 飛，于 瀾,等.非經(jīng)典阻尼系統(tǒng)的求解[J].吉林大學(xué)學(xué)報(bào)：工學(xué)版，2006，36：49-52.

XU Tao,CHENG Fei,YU Lan,et al.Solution of non-classical damped system [J].Journal of Jilin University: Engineering and Technology Edition,2006,36: 49-52.

[10]Li W M，Hong J Z.Research on the iterative method for model updating based on the frequency response function [J].Acta Mechanica Sinica,2012,28(2):450-457.

[11]楊彥芳，宋玉普，紀(jì)衛(wèi)紅.基于實(shí)測(cè)頻響函數(shù)主成分的在役網(wǎng)架損傷識(shí)別方法[J].振動(dòng)與沖擊，2007,26(9): 128-132.

YANG Yan-fang,SONG Yu-pu,JI Wei-hong.An identification method of existing truss structural damage bases on principal component analysis measured frequency response functions[J].Journal of Vibration and Shock,2007,26(9): 128-132.

[12]楊海峰，吳子燕，吳 丹.基于加速度頻率響應(yīng)函數(shù)的結(jié)構(gòu)損傷測(cè)量方法研究[J].振動(dòng)與沖擊，2007,26(2): 90-92.

YANG Hai-feng,WU Zi-yan,WU Dan.Structural damage detection method based on acceleration frequency response function [J].Journal of Vibration and Shock,2007,26(2): 90-92.

[13]Ward H，Stefan L，Paul S.Modal analysis theory and testing[M].Katholieke Universiteit Levven,Brussel,Belgium,1997.