基于頻響函數矩陣計算阻尼系統動力響應的新方法

張 淼， 于 瀾 ， 鞠 偉

(1.長春工程學院 理學院，長春 130012；2.中國第一汽車股份有限公司 技術中心，長春 130012)

計算結構動力響應時若采用直接積分法，對于每一個時間步長，其運算次數與半帶寬、自由度數的乘積成正比。當半帶寬較大且時間歷程遠大于系統的最小固有振動周期時，結構動力響應的計算將是很耗時的。而振型迭加法在一定條件下可以取得比直接積分法高的計算效率，它主要是利用系統自由振動的模態振型將多自由度的動力方程組轉換成為相互解耦的獨立方程，對每一個方程可以求解其響應的解析解和數值解。在對每個方程求解時常常采用Duhamel積分法，在一般情況下，它也需要數值積分來計算，只有極少數簡單情況才可以得到解析解。直接積分法和振型迭加法各有優勢，目前被工程界廣泛應用。

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

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

1 多自由度阻尼系統的轉化

對自由度為n的阻尼系統，設M、C、K分別是對稱的質量、阻尼和剛度矩陣，相應的運動微分方程為:

(1)

若取V=[{v1},…,{vn}]為無阻尼正則振型矩陣，則模態質量和模態剛度矩陣變為:

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

(2)

記模態阻尼矩陣為:

(3)

并取坐標變換

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

(4)

在模態坐標{q(t)}下，式(1)等價于:

(5)

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

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

(6)

或其逆變換:

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

(7)

則:

(8)

和

(9)

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

(10)

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

(11)

這時式(11)變為無阻尼振動系統的運動微分方程。為了討論其特征問題，引入狀態方程。

2 阻尼系統的復模態分析

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

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

(12)

令:

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

則:

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

即:

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

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

|Z(λ)|=0

(13)

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

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

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

(14)

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

(15)

其中:

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

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

(16)

其特征值滿足方程:

|A-λE|=0

(17)

(i=1,2,…,n)

基于上面關于一般動力系統的復模態理論，具體考慮無阻尼系統式(11)的極點和模態振型，并利用它們求系統式(11)的響應。

3 經典阻尼系統的精確響應求解

對經典阻尼系統，系統(1)的模態阻尼矩陣D為對角陣，所以矩陣D2也為對角陣，根據矩陣代數理論，矩陣函數f(D)=diag(f(d1),…,f(dn))也是對角陣，即

(18)

從而無阻尼系統式(11)的系數矩陣便約化為常數對角陣,即

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

相應地，在坐標變換式(7)下，式(5)等價于:

(19)

即由式(16)可有:

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

(20)

(r=1,2,…,n)

(21)

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

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

其中Qr為比例換算因子。

定理2[13]對阻尼系統(1)，其頻率響應矩陣為:

其中2jωrmrQr=1，mr為模態質量。

由于系統式(19)的模態向量為單位坐標向量，即VTV=I，且該系統的質量陣為單位陣，因此由定理3可得系統式(19)的頻響函數矩陣為對角陣:

(22)

由此，在簡諧激勵力作用下，當輸入力為{f}={F}eiωt時，系統式(19)的穩態響應為:

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

(23)

再根據式(4)和式(7)，可得原阻尼系統式(1)的響應為:

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

(24)

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

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

(25)

由式(25)可知只需利用無阻尼正則振型、無阻尼固有頻率及模態阻尼矩陣等信息即可求得簡諧激勵下的任一經典阻尼動力系統的精確響應。

4 非經典阻尼系統的精確響應求解

對非經典阻尼系統，無法直接使用振型迭加法求解，為此本文提出一種基于模態參數求解非經典阻尼精確響應的新方法。

對非經典阻尼系統，模態阻尼矩陣D不是對角陣，對應系統式(11)的狀態矩陣A為

與(20)式類似地可有

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

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

D=Pdiag[k1,…,kn]PT

(26)

根據矩陣代數理論，矩陣函數:

eDt=Pdiag[ek1t,…,eknt]PT

(27)

(r=1,2,…,n)

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

(28)

由此，在簡諧激勵條件下，當輸入力為{f}={F}eiωt時，系統式(11)的穩態響應為：

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

(29)

再由式(4)和式(7)可得系統式(1)的穩態響應為：

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

(30)

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

5 數值算例

作為算例1，考慮如圖所示的動力系統。

圖1 兩自由度阻尼振動系統

對上面兩自由度阻尼振動系統，令:

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

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

M1=M2=m=2 kg

(31)

解得系統的無阻尼正則振型矩陣為:

(32)

由于模態阻尼陣:

(33)

為對角陣，所以該系統為經典阻尼系統。由振型迭加法可推得響應的解析表達式為:

(34)

其中:

由本文提出的方法，

(35)

根據式(25)得:

(36)

表1 本文算法與振型迭加法計算響應解的比較

下面在時間6 s內，分別以時間步長Δt=0.25 s和Δt=1 s應用本文方法求解響應，并與振型迭加法的計算結果進行了比較，具體結果見表1。

從表1的結果看，對于經典阻尼系統的響應問題，本文方法的計算結果與振型迭加法的計算結果基本一致，而本文的響應表達式是通過解析方法得到的，這說明本文方法與振型迭加法的表達式是等價的。二者的偏差存在于振型迭加法需反復計算兩倍于系統自由度數的反正切和三角正弦值，在計算機運算過程中不可避免地產生了誤差的累積。

作為算例2，考慮如下的彈簧質量系統，如圖2所示。

圖2 彈簧質量系統

Fig.2 Spring-mass system

圖中k1=0.95，k2=k3=…=k10=0.03，k11=1.05，m1=m2=…=m10=1.0，每個頻率的自然阻尼系數都為0.05。用有限元方法提取系統的性質矩陣后，計算其無阻尼正則振型矩陣為:

由式(3)計算模態阻尼陣為:

表2 本文算法及Newmark法計算響應解的比較

D并非純對角矩陣，因此該系統為非經典阻尼系統。設簡諧激勵為:

{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。下面在時間2 s內應用本文方法求解響應，并與Newmark方法進行比較，結果見表2。

從表2的結果看，對于非經典阻尼系統的響應問題，本文方法與Newmark方法求得的響應值相近，即驗證了本文方法的正確及可行性。二者的偏差存在于本文方法為解析解，而Newmark方法為數值解。而在計算效率方面，Newmark方法在每個時間步長上都需要求解一個未知數個數與自由度數相同的方程組，為了提高精度還需要縮小時間步長，勢必帶來計算量的增長，因此本文方法在計算效率及精度方面均優于Newmark方法。

6 結 論

本文提出了一種計算經典及非經典阻尼系統響應的新方法。只需利用任意動力系統所對應的無阻尼正則振型、無阻尼固有頻率及模態阻尼陣，即可求得任意激勵下經典及非經典阻尼動力系統的精確響應。本文提出的算法公式簡潔、緊湊和精確，不存在計算誤差，不僅在計算響應時使效率及精度均得到提高，而且還可作為優化目標函數應用于模型修正、結構設計優化及損傷識別等工程領域，有著良好的前景。

參 考 文 獻

[1]徐醫培，李素有，吳立言.結構動態響應的求解方法分析[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]張森文，曹開彬.計算結構動力響應的狀態方程直接積分法[J].計算力學學報，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]劉寒冰，張 淼，魏 健.結構動力響應分析的多重網格方法[J].吉林大學學報：工學版，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]蒲軍平，劉 巖，王元豐，等.基于精細時程積分的結構動力響應降維分析[J].清華大學學報：自然科學版，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]郭興旺，鄒家祥.對機械振動系統的六種動態響應分析方法的評述[J].振動與沖擊，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]周錫元，馬東輝，俞瑞芳.工程結構中的阻尼與復振型地震響應的完全平方組合[J].土木工程學報，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].吉林大學學報：工學版，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].振動與沖擊，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]楊海峰，吳子燕，吳 丹.基于加速度頻率響應函數的結構損傷測量方法研究[J].振動與沖擊，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.