曹艷君 王 皓

（復(fù)旦大學(xué)航空航天系，上海 200433）

基于快速Fourier變換法的廣義特征值問(wèn)題重根辨識(shí)方法?

曹艷君 王 皓?

（復(fù)旦大學(xué)航空航天系，上海 200433）

對(duì)具有重根的廣義特征值問(wèn)題，采用基于快速Fourier變換的方法進(jìn)行求解，實(shí)現(xiàn)重根辨識(shí).文章中采用多次單點(diǎn)初始激勵(lì)的方式，仿真計(jì)算測(cè)點(diǎn)上的自由振動(dòng)響應(yīng)，對(duì)響應(yīng)進(jìn)行快速Fourier變換后得到頻域數(shù)據(jù).而后對(duì)頻域數(shù)據(jù)分析，得到固有頻率和多組測(cè)點(diǎn)振型數(shù)據(jù).根據(jù)單頻和重頻處的振型特性，引入振型的余弦相似度為判別參數(shù)，辨識(shí)重根.數(shù)值算例表明，該方法可有效實(shí)現(xiàn)重根辨識(shí)，同時(shí)特征值的計(jì)算能達(dá)到較高精度.

廣義特征值問(wèn)題， 重根辨識(shí)， 快速Fourier變換法， 固有頻率， 動(dòng)力學(xué)響應(yīng)

引言

廣義特征值問(wèn)題的求解是結(jié)構(gòu)動(dòng)力分析和穩(wěn)定性分析的關(guān)鍵.求解廣義特征值問(wèn)題的傳統(tǒng)算法有Lanczos 法、子空間迭代法、Ritz 向量法［1，2］和 Jacobi算法等.近年來(lái)也涌現(xiàn)出一些新興算法，如混合人工魚(yú)群算法［3］、神經(jīng)網(wǎng)絡(luò)算法［4，5］、瀑布型多重網(wǎng)格法［6］等.

2014年，吳鋒、徐小明和鐘萬(wàn)勰提出將快速Fourier變換法用于求解廣義特征值問(wèn)題［7］.該方法的主要思想是，采用計(jì)算機(jī)仿真計(jì)算結(jié)構(gòu)的自由振動(dòng)響應(yīng)，對(duì)響應(yīng)進(jìn)行快速Fourier變換，從頻域數(shù)據(jù)可直觀地得到特征值信息.但對(duì)于具有重頻的廣義特征值問(wèn)題，作者并未敘及.

重頻是多自由度系統(tǒng)中固有頻率重疊的現(xiàn)象，在動(dòng)力學(xué)分析中較為常見(jiàn).在求解廣義特征值問(wèn)題過(guò)程中若沒(méi)有考慮重頻存在的可能性，求解時(shí)會(huì)遺漏重頻處的振型.這會(huì)導(dǎo)致分析結(jié)果誤差增大，甚至得到完全錯(cuò)誤的結(jié)果［8］.因此，在廣義特征值問(wèn)題的求解中，有必要對(duì)重頻進(jìn)行辨識(shí).

本文主要對(duì)具有重頻的廣義特征值問(wèn)題進(jìn)行研究，采用多次單點(diǎn)初始激勵(lì)，仿真計(jì)算測(cè)點(diǎn)上的動(dòng)力學(xué)響應(yīng).在對(duì)頻域數(shù)據(jù)進(jìn)行分析時(shí)，除提取固有頻率信息以外，還同時(shí)得到其對(duì)應(yīng)的振型.得到固有頻率和多組振型數(shù)據(jù)后，根據(jù)單頻和重頻處振型的不同特性，以振型向量的余弦相似度為判別參數(shù)辨識(shí)重頻.

1 重頻結(jié)構(gòu)與仿真計(jì)算

1.1 重頻結(jié)構(gòu)的自由振動(dòng)問(wèn)題

無(wú)阻尼系統(tǒng)的自由振動(dòng)方程為:

其中x為結(jié)構(gòu)的位移向量，M為質(zhì)量矩陣，K為剛度矩陣，n為系統(tǒng)的總自由度數(shù)．

本文僅討論K為正定矩陣的情形，此時(shí)ωi大于零．

若系統(tǒng)固有頻率中存在重頻，假設(shè)ωm為k重頻（k＞1），即:

對(duì)應(yīng)的 k 個(gè)主振型記為 φm，φm+1，…，φm+k-1．系統(tǒng)的自由振動(dòng)的解寫(xiě)成以下形式:x（t） ＝

對(duì)于非零位移、零速度初始條件，系統(tǒng)自由振動(dòng)的響應(yīng)式為:

對(duì)于第j個(gè)自由度，有響應(yīng):

式（6）中，單頻 ωs前系數(shù)為而重頻ωm前系數(shù)為在仿真計(jì)算中，這些系數(shù)可通過(guò)模態(tài)參數(shù)辨識(shí)方法得到．

第二次計(jì)算 ωs處的振型為＝，而ωm處的振型為:

由此可知，對(duì)于具有二重頻的問(wèn)題，單頻處兩次計(jì)算得到的振型歸一化后是一致的．而重頻處，只要:

得到的兩個(gè)振型就線性獨(dú)立．于是，根據(jù)此特性可以對(duì)重頻進(jìn)行判別．

自由振動(dòng)響應(yīng)算法可采用精細(xì)積分法［9］、Runge-Kutta法［10］等．實(shí)際計(jì)算中，無(wú)需在所有的自由度上得到完整的振型．只需在結(jié)構(gòu)中選取若干個(gè)測(cè)點(diǎn)，得到這些測(cè)點(diǎn)上的振型，便可識(shí)別出重頻．

1.2 快速Fourier變換及振型提取

在利用快速Fourier變換進(jìn)行仿真計(jì)算時(shí)，要求數(shù)據(jù)在時(shí)域和頻率均是離散的，且為有限長(zhǎng)［11］．若采樣間隔為Δt，采樣點(diǎn)數(shù)為Nt，采樣得到的響應(yīng)數(shù)據(jù)序列為 xj（r）（r＝ 0，1，2，…，Nt-1），xj（r）做快速Fourier變換后得到的頻域數(shù)據(jù)為Xj（l），有:

其中 Aj（l）為 Xj（l）的幅值，φj（l）為相位．

對(duì)Xj（l）提取峰值便可得到結(jié)構(gòu)的固有頻率．同時(shí)，分析頻域數(shù)據(jù)可得到振型．此處的振型指在選取的測(cè)點(diǎn)上的振型，而非完整振型．

1.3 采樣時(shí)間間隔Δt和采樣點(diǎn)數(shù)Nt

仿真計(jì)算時(shí)，采樣時(shí)間間隔Δt和采樣點(diǎn)數(shù)Nt的選取需要滿足一定要求．

選取采樣間隔Δt之前，要對(duì)結(jié)構(gòu)的最大頻率作估計(jì)．由采樣定理可知，采樣頻率應(yīng)至少是最大頻率的兩倍．基于廣義蓋爾圓定理［12］，有最大特征值:

仿真計(jì)算時(shí)同樣需要考慮頻率分辨率．頻率分辨率指算法能將信號(hào)中兩個(gè)相近的譜峰保持分開(kāi)的能力．若信號(hào)由兩個(gè)相近頻率組成，分別為ω0和ω0+Δω0．當(dāng) Δω0很小時(shí)，可能導(dǎo)致兩個(gè)峰合在一起．此時(shí)需要增加采樣點(diǎn)數(shù)Nt，使得頻域數(shù)據(jù)兩峰之間出現(xiàn)明顯的波谷．根據(jù)文獻(xiàn)［7］，采樣點(diǎn)數(shù)Nt可通過(guò)下式確定:

1.4 振型向量相似度度量

兩次計(jì)算在不同的點(diǎn)上給定初始條件，得到兩組頻率和測(cè)點(diǎn)上振型的數(shù)據(jù)．滿足一定條件時(shí)，可望重頻處得到的兩個(gè)振型向量線性無(wú)關(guān)．此處引入振型向量的余弦相似度刻畫(huà)振型的相似度．若有兩個(gè) n 維向量 μ＝｛μ1，μ2，…，μn｝T和 η＝｛η1，η2，…，ηn｝T，則有μ和η夾角的余弦值為:

cosθ的值越接近于1，兩向量的相似度越高．| cosθ＝1時(shí)， μ 和 η 可以互相線性表出；cosθ＝0時(shí)，μ和η相互正交．

理論上，單頻處兩個(gè)振型向量的夾角余弦值|cosθ|＝1，而重頻處，夾角余弦值|cosθ|＜1．由此，以振型向量的余弦相似度為判別參數(shù)，可實(shí)現(xiàn)重頻的判別．

2 算例驗(yàn)證與分析

2.1 算法

以具有二重根的廣義特征值問(wèn)題為例，具體的算法如下:

（1）選定s個(gè)測(cè)點(diǎn)，確定識(shí)別指標(biāo)h和頻率分辨率Δω0的值，由式（13）和式（14）確定時(shí)間間隔Δt和采樣點(diǎn)數(shù)Nt．

（2）取一個(gè)初始激勵(lì)點(diǎn)，給定初始條件．

（3）計(jì)算給定初始條件下測(cè)點(diǎn)的位移響應(yīng)，對(duì)響應(yīng)數(shù)據(jù)進(jìn)行快速Fourier變換，得到頻域數(shù)據(jù)．

（4）從頻域數(shù)據(jù)中提取固有頻率和振型信息，振型做歸一化處理．

（5）另取一個(gè)初始激勵(lì)點(diǎn)，給定初始條件，執(zhí)行步驟（3）和（4），得到另一組頻率和振型數(shù)據(jù)．

（6）計(jì)算每一階頻率兩振型向量的夾角余弦值，判別重頻．

2.2 算例驗(yàn)證

如圖1所示，有一四邊固支的對(duì)稱鋼板結(jié)構(gòu)，尺寸為 1m×1m×0.002m，彈性模量 E＝2×1011Pa，泊松比 μ＝0．3，密度 ρ＝ 7800kg／m3．采用矩形板單元進(jìn)行網(wǎng)格劃分（板單元的每個(gè)節(jié)點(diǎn)上有3個(gè)自由度，分別為撓度和x方向、y方向的轉(zhuǎn)角）．如圖，劃分網(wǎng)格后共有8×8＝64個(gè)單元，81個(gè)節(jié)點(diǎn)．節(jié)點(diǎn)編號(hào)方式為:左下角節(jié)點(diǎn)為1號(hào)節(jié)點(diǎn)，沿左側(cè)邊向上編號(hào)，至左上角點(diǎn)為9號(hào)節(jié)點(diǎn)，而后以底邊左面第二個(gè)節(jié)點(diǎn)為10號(hào)節(jié)點(diǎn)，依次向上編號(hào)，以此類(lèi)推，直至編至右上角點(diǎn)為81號(hào)節(jié)點(diǎn)．

圖1 鋼板結(jié)構(gòu)Fig.1 Steel plate

選擇5個(gè)節(jié)點(diǎn)計(jì)算響應(yīng)，節(jié)點(diǎn)編號(hào)為20、29、51、62和71號(hào)．第一次計(jì)算時(shí)，在22號(hào)節(jié)點(diǎn)撓度方向給定大小為0．05m初位移，第二次在65號(hào)節(jié)點(diǎn)撓度方向給定同樣大小的初位移．兩次計(jì)算初速度均為零．

取 Δω0＝0．5，h＝5，相應(yīng)的采樣時(shí)間間隔 Δt＝2×10-4s，計(jì)算響應(yīng)的數(shù)值方法為四階Runge-Kutta法．

提取頻域數(shù)據(jù)的峰值得到固有頻率，綜合多個(gè)測(cè)點(diǎn)的數(shù)據(jù)，得到如表1第二列和第四列所示的固有頻率結(jié)果．

根據(jù)式（10）和式（11）從頻域數(shù)據(jù)中提取振型，對(duì)于每個(gè)測(cè)點(diǎn)，其撓度方向計(jì)入振型．

其振型向量夾角余弦值為 cosθ1＝1．000．

夾角余弦值cosθ2＝-0．235．由此，可以判定處存在重頻．

由于篇幅原因，其余各階頻率處的振型不再一一列出，振型向量的夾角余弦值如表1第三列和第六列所示．

表1 固有頻率和振型夾角余弦值Table 1 Natural frequency and cosine values

表2 計(jì)算結(jié)果與精確解Table 2 Calculated and exact results

由表 1 可知，頻率 35.138Hz、76.550Hz、101.663Hz、133.225Hz和160.050Hz是重頻.用Matlab的eig函數(shù)計(jì)算精確解，將計(jì)算結(jié)果和精確解進(jìn)行對(duì)比.表2中列出了前20階固有頻率的結(jié)果（帶?的為重頻）.由此可見(jiàn)，本方法能較好地實(shí)現(xiàn)重頻辨識(shí)，同時(shí)固有頻率計(jì)算達(dá)到較高的精度.

2.3 分析與討論

對(duì)于具有二重頻的問(wèn)題，以上算法需要進(jìn)行兩次響應(yīng)計(jì)算．理論上，只計(jì)算一次響應(yīng)也可實(shí)現(xiàn)重頻的辨識(shí)．計(jì)算響應(yīng)的初始條件為:在一個(gè)點(diǎn)上給定初位移，同時(shí)另一個(gè)點(diǎn)上給定初速度．此時(shí)重頻處實(shí)部分量得到的振型?φreal和虛部分量得到的振型?φimag線性無(wú)關(guān)．理論上可據(jù)此實(shí)現(xiàn)重頻的判別，但實(shí)際操作時(shí)，重頻辨識(shí)結(jié)果受初始條件的影響較大，效果不理想．因此，對(duì)于具有二重頻的廣義特征值問(wèn)題，本文仍選擇求兩次響應(yīng)的方式．

若研究具有k重頻（k≥2）的問(wèn)題，首先，可按照二重頻的方法，計(jì)算固有頻率并辨識(shí)重頻．在每階重頻處，可得兩個(gè)線性無(wú)關(guān)的振型．而后繼續(xù)選用新的初始激勵(lì)點(diǎn)，計(jì)算動(dòng)力學(xué)響應(yīng)，并得到新的一組振型數(shù)據(jù)．若重頻處的振型不可用之前的振型線性表出，則表示該階頻率處又找到一個(gè)新的振型．繼續(xù)選用新的初始激勵(lì)點(diǎn)，計(jì)算響應(yīng)，直到重頻處的振型均可用之前的振型線性表出．此時(shí)，重頻處得到的線性無(wú)關(guān)的振型的數(shù)目即是該階頻率的重?cái)?shù)．

3 結(jié)論

本文采用多次單點(diǎn)初始激勵(lì)的方式，從響應(yīng)的頻域數(shù)據(jù)中得到多組測(cè)點(diǎn)上的振型，引入振型的余弦相似度為判別參數(shù)辨識(shí)重根.以具有二重頻的問(wèn)題為例進(jìn)行特征值計(jì)算和重頻辨識(shí).計(jì)算結(jié)果表明，振型向量的余弦值在單頻和重頻處有明顯差異，能夠根據(jù)余弦值有效辨識(shí)重頻，同時(shí)特征值計(jì)算可達(dá)到較高精度.

1 Li R C，Zhang L H.Convergence of the block Lanczos method for eigenvalue clusters.Numerische Mathematik，2015，131（1）:83～113

2 宮玉才，周洪偉，陳璞等.快速子空間迭代法、迭代Ritz向量法與迭代Lanczos法的比較.振動(dòng)工程學(xué)報(bào)，2005，18（2）:227～232 （Gong Y C， Zhou H W， Chen P， et al.Comparative study of fast subspace iteration method，Ritz method and Lanczos method.Journal of Vibration Engineering， 2005，18（2）:227～232 （in Chinese））

3 黃華娟，周永權(quán)，韋杏瓊等.求解矩陣特征值的混合人工魚(yú)群算法.計(jì)算機(jī)工程與應(yīng)用，2010，46（6）:56～59（Huang H J， Zhou Y Q， Wei X Q， et al.Artificial fish school method for eigenvalue problems.Computer Engineering and Applications， 2010，46（6）:56～59 （in Chinese））

4 Nandy S，Sharma R，Bhattacharyya S P.Solving symmetric eigenvalue problem via genetic algorithms:Serial versus parallel implementation.Applied Soft Computing， 2011，11（5）:3946～3961

5 Tang Y，Li J P.Another neural network based approach for computing eigenvalues and eigenvectors of real skew-symmetric matrices.Computers＆Mathematics with Applications， 2010，60（5）:1385～1392

6 Chen H T，He Y H，Li Y，et al.A multigrid method for eigenvalue problems based on shifted-inverse power technique.European Journal of Mathematics， 2015，1（1）:207～228

7 吳鋒，徐小明，鐘萬(wàn)勰.廣義特征值問(wèn)題的快速傅里葉變換法.振動(dòng)與沖擊， 2014，33（22）:67～71 （Wu F， Xu X M，Zhong W X.Fast Fourier transform method for generalized eigenvalue problems.Journal of Vibration and Shock，2014，33（22）:67～71 （in Chinese））

8 劉福林，閆維明，何浩祥.重頻結(jié)構(gòu)的振型識(shí)別及處理方法.工業(yè)建筑， 2011，41（S1）:368～372 （Liu F L， Yan W M，He H X.Mode shape identification and disposal in the structure of overlap frequencies.Industrial Construction， 2011，41（S1）:368～372 （in Chinese））

9 張繼鋒，鄧子辰，徐方暖等.一種新的改進(jìn)精細(xì)直接積分法.動(dòng)力學(xué)與控制學(xué)報(bào)， 2015，13（4）:241～245 （Zhang J F，Deng Z C，Xu F N，et al.A new improved precise direct integration method.Journal of Dynamics and Control，2015，13（4）:241～245 （in Chinese））

10 吳志橋，高普云，任鈞國(guó).Runge-Kutta方法求解結(jié)構(gòu)動(dòng)力學(xué)方程.系統(tǒng)仿真學(xué)報(bào)，2010，22（9）:2085～2090（Wu Z Q，Gao P Y，Ren J G.Runge-Kutta methods for time integration in structural dynamics.Journal of System Simulation， 2010，22（9）:2085～2090 （in Chinese））

11 陳恒亮，蔣勇.基于DSP的實(shí)數(shù)FFT算法研究與實(shí)現(xiàn).動(dòng)力學(xué)與控制學(xué)報(bào)， 2005，3（2）:52～55 （Chen H L，Jiang Y.Design and realization of real FFT based on DSP.Journal of Dynamics and Control， 2005，3（2）:52～55 （in Chinese））

12 Nakatsukasa Y.Gerschgorin′s theorem for generalized ei-genvalue problems in the Euclidean metric.Mathematics of Computation， 2011，80（276）:2127～2142

A METHOD OF MULTIPLE-FREQUENCY IDENTIFICATION FOR GENERALIZED EIGENVALUE PROBLEMS BASED ON FAST FOURIER TRANSFORM

Cao Yanjun Wang Hao?
（Department of Aeronautics and Astronautics，F(xiàn)udan University，Shanghai 200433，China）

A method based on the fast Fourier transform was proposed to solve the generalized eigenvalue problems with multiple roots.This paper studied the dynamic responses of the measure nodes with nonzero initial condition on one point.Both natural frequencies and mode shapes were extracted from the data in the frequency domain.Responses under different initial conditions were calculated to get several sets of mode shapes.Taking the cosine similarity of mode shapes as the discriminant parameter，multiple roots were then identified.The numerical example shows that this method can identify the multiple roots efficiently，and the result reaches a high accuracy.

generalized eigenvalue problem， multi-frequency identification， fast Fourier transform， natural frequency， dynamic response

13 September 2016，revised 21 December 2016.

10.6052／1672-6553-2017-015

2016-09-13收到第1稿，2016-12-21收到修改稿.

?國(guó)家自然科學(xué)基金資助項(xiàng)目（11572089）

?通訊作者 E-mail:wanghao＠fudan.edu.cn

?The project supported by the National Natural Science Foundation of China（11572089）

? Corresponding author E-mail:wanghao＠fudan.edu.cn