張 淼， 于 瀾 ， 鞠 偉

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

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

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

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

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

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

(1)

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

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

(2)

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

(3)

并取坐標變換

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

(4)

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

(5)

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

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

(6)

或其逆變換:

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

(7)

則:

(8)

和

(9)

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

(10)

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

(11)

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

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

對一般動力系統(tǒng)式(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)

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

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

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

(18)

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

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

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

(19)

即由式(16)可有:

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

(20)

(r=1,2,…,n)

(21)

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

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

其中Qr為比例換算因子。

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

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

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

(22)

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

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

(23)

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

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

(24)

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

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

(25)

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

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

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

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

與(20)式類似地可有

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

這是關于矩陣eDt(D2-Λ)e-Dt的特征問題表達式。模態(tài)阻尼矩陣D雖不是對角陣，但它為對稱陣，因此可正交相似于對角陣，即存在正交陣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)的質量陣為單位陣，因此由定理3系統(tǒng)式(11)的頻響函數(shù)矩陣可簡化為對角陣：

(28)

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

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

(29)

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

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

(30)

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

5 數(shù)值算例

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

圖1 兩自由度阻尼振動系統(tǒ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)的無阻尼正則振型矩陣為:

(32)

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

(33)

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

(34)

其中:

由本文提出的方法，

(35)

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

(36)

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

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

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

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

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

Fig.2 Spring-mass system

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

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

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

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

{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內(nèi)應用本文方法求解響應，并與Newmark方法進行比較，結果見表2。

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

6 結 論

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

參 考 文 獻

[1]徐醫(yī)培，李素有，吳立言.結構動態(tài)響應的求解方法分析[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]張森文，曹開彬.計算結構動力響應的狀態(tài)方程直接積分法[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]劉寒冰，張 淼，魏 健.結構動力響應分析的多重網(wǎng)格方法[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]郭興旺，鄒家祥.對機械振動系統(tǒng)的六種動態(tài)響應分析方法的評述[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īng)典阻尼系統(tǒng)的求解[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]楊彥芳，宋玉普，紀衛(wèi)紅.基于實測頻響函數(shù)主成分的在役網(wǎng)架損傷識別方法[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]楊海峰，吳子燕，吳 丹.基于加速度頻率響應函數(shù)的結構損傷測量方法研究[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.