廣義特征值問題的快速傅里葉變換法

吳 鋒，徐小明，鐘萬勰

（大連理工大學工業(yè)裝備結構分析國家重點實驗室，工程力學系，遼寧大連 116024）

廣義特征值問題的快速傅里葉變換法

吳 鋒，徐小明，鐘萬勰

（大連理工大學工業(yè)裝備結構分析國家重點實驗室，工程力學系，遼寧大連 116024）

針對廣義特征值問題提出離散傅里葉變換法。該方法把結構的動力響應看作是一種信號，利用快速傅里葉變換進行分析，從而得到結構的振動頻率。該方法避免對剛度矩陣求逆，可同時計算出所有的特征值，是一種直接方法。數值算例驗證了該方法的正確性。

特征值；動態(tài)結構響應；快速傅里葉變換；采樣

在工程中，結構振動頻率分析和穩(wěn)定分析都涉及到廣義特征值問題。關于廣義特征值問題已經提出許多算法［1－9］，其中常用的算法有子空間迭代法、Ritz向量法、Lanczos算法等［4，6］。這些方法均為迭代方法，通常是利用一組基，把大規(guī)模的廣義特征值問題轉化成一個小規(guī)模的特征值問題，然后再用直接解法求解這一小規(guī)模特征值問題。對于小規(guī)模的特征值問題，通常采用Jacobi算法［2］。一般說來，迭代方法必須結合直接算法才可以進行，而且需要對剛度矩陣求逆，往往只能求解部分特征值。

對于一個結構，其結構振動響應中包含著關于結構振動頻率的所有信息，如果把結構的動力響應看作是一組信號，則可以利用信號處理中的快速傅里葉變換進行分析［10，11］，從而找到結構的振動頻率。在實際結構的模態(tài)分析中，常通過對測量得到的結構振動響應進行傅里葉變換，分析其結構的振動頻率［12］。實際上，這一思想也可用于計算廣義特征值問題。對于一組給定的剛度矩陣和質量矩陣，可以人為構造一種工況，采用計算機仿真計算其動力響應，然后利用快速傅里葉變換進行分析，從而得到特征值。關于結構的動力響應計算有許多成熟算法［13－16］?；谶@一認識，本文針對廣義特征值問題提出基于結構動力響應的離散傅里葉變換法。該方法把結構的動力響應看作是一種信號，而結構的動力響應可以利用成熟的算法如中心差分法、精細積分法等求解，也可以實際測量，從而避免對剛度矩陣求逆。得到結構動力響應后，利用快速傅里葉變換進行分析。本文方法可以同時計算出所有的特征值，是一種直接方法，也可以與迭代方法結合求解，處理迭代法中需要解決的小規(guī)模特征值問題。

1 基本思想

本文考慮結構的廣義特征值問題

式中：K，M為結構剛度、質量矩陣；x為特征向量；ωi為振動頻率。且有

求解式（1）的廣義特征值問題常用方法為子空間迭代、Lanczos迭代等。本文提出與此完全不同之方法，即由結構振動響應出發(fā)，對結構振動響應進行傅里葉變換，獲得結構振動頻率。

1.1 自由振動

設某工況下結構振動響應只含結構振動頻率，不含其它任何雜質，即無阻尼系統(tǒng)的自由振動。設結構位移向量a，其動力微分方程為

由式（9）可見結構位移響應中只含結構振動頻率，對結構位移響應進行傅里葉變換便可獲得結構的振動頻率。

1.2 時間步長及采樣長度

以上分析均基于連續(xù)情況，實際計算結構動力響應時只能得到各時間點位移值，故可用快速離散傅里葉變換（FFT）。計算結構動力響應應選精度較好的動力算法。對規(guī)模較大動力問題需避免求逆。獲得結構在各時間點位移值后須用FFT進行分析，其中時間步長取決于位移中所含最大頻率ωmax，而采樣長度則取決于位移中兩相鄰頻率之差Δω。設需分析信號為

根據式（15）知，如果頻率w為整數時，am＝δmw。因此根據am是否等于1，就可以找到頻率w。但是因為m≤0．5N，如果w＞0．5N，則w≠m，此時對于所有的m有am＝0，無法識別出信號的頻率，而其截斷N項的傅里葉級數為0，傅里葉信號失真。因此僅當w≤0．5N時才可以識別出信號的頻率，而其離散傅里葉變換在各個時間點上的值才不會失真。把w≤0．5N代入式（14）可得：

am是關于w的一個復函數，其模為：

因此，｜am（w）｜是一個關于w和m的函數，當固定w時，m越接近w，｜am｜越大，而遠離m時，｜am｜只有N－1的量級?，F在假設一個信號有兩個頻率組成，一個為ω，另一個為ω＋Δω，該信號離散傅里葉變換后的各個頻率的振幅即為：

只要函數｜bm（w）｜出現兩個峰值，一個在w處，一個在w＋Δw處即可以識別出這兩個峰值。當Δw很小時，由于w與w＋Δw相處太近，導致兩個峰值合成一個峰值，識別不出這兩個頻率。這是頻率分辨率的問題，解決的方法是增加采樣長度N，使得在［w］（［w］表示距離w最近的整數）和［w＋Δw］之間出現明顯的波谷，從而把兩個峰值明顯的隔開?，F在假設和［w＋Δw］的平均數的整數，要求為一種識別指標），因為

進一步，上式也可以寫成ΔωΔt≥8h／N，這一不等式就是測不準原理。當采樣點數N固定后，ΔωΔt大于一定值，因此，時域的分辨率Δt越小，則頻域的分辨率Δω就大，反之亦然，要想兩者都小，則必須增加采樣點數。

1.3 最大頻率

確定時間步長時需要用到結構的最大振動頻率，這一最大振動頻率可以根據結構的最小單元的最大特征值確定。當結構的質量矩陣是集中質量陣時，這里基于廣義蓋爾圓定理［17］，給出一個確定最大頻率的方法。在標準特征值問題中有個蓋爾圓定理，那么對于如式（1）所示的廣義特征值問題，我們可以很容易的建立一個廣義蓋爾圓定理，現在設：

剛度矩陣K中第i行第j列的元素為kij，質量矩陣M為對角矩陣，第i個對角元素為mi，x中第i個元素為xi，則矩陣形式的特征方程可以展開成如下所示的N個方程：

2 算法設計

根據式（38）和式（31）確定時間步長Δt和采樣長度N。然后根據式（3）和式（4）計算結構的響應，把響應進行離散傅里葉變換，既可以找到結構的各階振動頻率。在式（38）中，h是一種識別指標，描述了兩個相鄰頻率w和w＋Δw的分辨率。當h取得大時，｜bm（w）｜在這兩個頻率之間會出現一個明顯的波谷，從而可以識別出這兩個頻率。但是h取值很大時，采樣長度就會變大。經過計算發(fā)現，取h＝5較為合理。當通過傅里葉變換找到各個峰值時，此峰值所對應的頻率w為無量綱頻率，還需要通過式（14）化成圓頻率，也就是結構的振動頻率。下面講述具體計算過程。

（1）根據式（38）和式（31）確定時間步長和采樣長度N；

（2）選定初始位移M－1a0，根據式（3）計算結構的自由振動響應；

（3）計算出結構的動力響應之后，任意提取某節(jié)點的位移響應作為信號，對該信號進行快速傅里葉變換，得到頻域信號；

（4）對頻域信號取絕對值，找出頻域信號峰值所在的無量綱頻率w；

（5）根據式（14），把w化成圓頻率，即得到結構的振動頻率。

3 算例

計算圖1所示20層平面鋼框架結構的振動頻率。已知框架柱為450×450×25 mm箱型截面，彈性模量E＝2．06× 105MPa。按剛性樓板假設，采用層剪切模型，集中質量，每層的質量m＝168 750 kg。分別采用傳統(tǒng)求解特征值問題的QR法、Jacobi方法和本文方法計算。QR法和Jacobi方法屬于迭代法，把QR法相對誤差取10－13時的解作為標準解。本文方法計算時，初始位移中a0的各項元素均取為1，識別指標取h＝5，計算所用的時間步長按式（38）選取，而頻率分辨率取為Δω＝0．5。結構的位移響應分別采用精細積分法［11］、中心差分法和文獻［12］中所提出的的高精度攝動法計算。經過計算發(fā)現，把任一層的位移響應作為時域信號進行快速傅里葉變換，計算結果完全一樣，這是因為其時域信號中所包含的頻率分量完全相同，這里只給出采用第一層位移響應作為時域信號分析的結果。所有的計算結果列于表1。

圖1 平面框架Fig．1 Plane frame

表1 計算的頻率及相對誤差Tab.1 Com puted frequencies and relative errors

由表1可見，不管采用哪種方法計算結構的自由振動響應，計算得到的頻率完全一樣，這說明當時間步長和采用長度選定后，結構的自由振動響應計算手段對計算結果的影響不大。與標準解相比，各階頻率的誤差均很小。誤差最大的為第1個頻率，也只有0．42%，這表明本文方法的可靠性。為進一步比較本文方法和經典的QR法和Jacobi方法的性能，分別取QR法和Jacobi方法的迭代相對誤差為0．4%，比較計算時間，結果列于表2。表2中，本文方法計算結構動力響應時采用的是精細積分方法。根據表2可見，對于小規(guī)模的廣義特征值問題，當相對誤差均為0．4%時，本文方法僅需要很少的計算時間即可以給出十分滿意的結果。這說明本文方法比較適用于小規(guī)模廣義特征值問題的求解。根據表1的計算結果還可見，本文方法計算的高階頻率的相對誤差要小于低階頻率的相對誤差，這一特點不同于傳統(tǒng)的迭代方法，在傳統(tǒng)迭代方法中，往往高階頻率的相對誤差大于低階頻率的相對誤差。

表2 計算時間Tab.1 Computation times

4 結 論

本文針對廣義特征值問題提出基于結構動力響應的離散傅里葉變換法。該方法把結構的動力響應看作是一種信號，利用快速傅里葉變換進行分析，從而得到結構的振動頻率。在本文方法中，結構的動力響應可以利用成熟的算法如中心差分法、精細積分法等求解，也可以實際測量，從而避免對剛度矩陣求逆。本文方法可以同時計算出所有的特征值，是一種直接方法，也可以與迭代方法結合求解，處理迭代法中需要解決的小規(guī)模特征值問題。必須指出的是，本文方法是把工程中模態(tài)分析的做法進行數值化，從而用于計算廣義特征值問題。對于一個復雜的真實結構，其動力響應計算往往要比動力特性計算費時，從這一點上講，本文方法比較適用于小規(guī)模廣義特征值問題。但是從兩點上講本文方法仍然是有價值的。

（1）傳統(tǒng)的特征值計算方法如QR法、Jacobi方法、子空間迭代法等均為迭代方法，其收斂的理論依據均為反冪法，而本文方法則是基于傅里葉變換理論，因此本文方法是對現有的廣義特征值問題的數值算法的一種補充，從數值算法角度上講，本文方法是有價值的。

（2）本文算法把廣義特征值問題看成是結構動力響應問題和信號分析問題來處理，具有鮮明的工程背景。數值算例表明，對于小規(guī)模的廣義特征值問題，本文方法僅需要很少的計算時間即可以給出十分滿意的結果。在本文方法的求解過程中，結構動力響應的計算手段對特征值的計算影響不大。

［1］薛長峰，周樹荃．非對稱廣義特征值問題的并行連續(xù)同倫算法［J］．計算物理，1997，14（4／5）：619－621．

XUE Chang-feng，ZHOU Shu-quan．Parallel homotopy continuation algorithm for nonsymmetric generalized eigenvalue problem［J］．Chinese Journal of Computational Physics，1997，14（4／5）：619－621．

［2］王順緒，戴華．廣義特征值問題的并行塊Jacobi-Davidson方法及應用［J］．計算力學學報，2008，25（4）：428－433．

WANG Shun-xu，DAI Hua．Parallel block Jacobi-Davidson method for solving large generalized eigenvalue problems and it's application［J］．Chinese Journal of Computational Mechanics，2008，25（4）：428－433．

［3］劉寒冰，龔國慶，劉建設．一種有效的廣義特征值分析方法［J］．固體力學學報，2003，24（4）：419－428．

LIU Han-bing，GONG Guo-qing，LIU Jian-she．An effective algorithm for generalized eigenvalue problems［J］．ACTA Mechanica Solida Sinica，2003，24（4）：419－428．

［4］李秀梅，吳鋒．廣義特征值問題求解的改進Ritz向量法［J］．振動與沖擊，2012，31（7）：19－23．

LI Xiu-mei，WU Feng．Improved Ritz vector method for generalized eigenvalue problems［J］．Journal of Vibration and Shock，2012，31（7）：19－23．

［5］李海濱，黃洪鐘，趙明揚．有限元結構動力分析的廣義特征值的神經計算［J］．哈爾濱工業(yè)大學學報，2006，38（9）：1523－1527．

LI Hai-bin，HUANG Hong-zhong，ZHAO Ming-yang．Finite element approach for structural dynamic analysis using generalized eigenvalue-based neurocomputing［J］．Journal of Harbin Iistitute of Technology，2006，38（9）：1523－1527．

［6］黃吉鋒．求解廣義特征值問題的多重Ritz向量法［J］．力學學報，1999，31（5）：585－595．

HUANG Ji-feng．The multi-ritz-vector method in generalized eigenvalue problems［J］．Chinese Journal of Theoretical and Applied Mechanics，1999，31（5）：585－595．

［7］黃華娟，周永權，韋杏瓊，等．求解矩陣特征值的混合人工魚群算法［J］．計算機工程與應用，2010，46（6）：56－59．

HUANG Hua-juan，ZHOU Yong-quan，WEIXing-qiong，et al．Hybrid artificial fish school algorithm for solving matrix eigenvalues［J］．Computer Engineering and Applications，2010，46（6）：56－59．

［8］Nandy S，Sharma R，Bhattacharyya SP．Solving symmetric eigenvalue problem via genetic algorithms：serial versus parallel implementation［J］．Applied Soft Computing，2011，11（5）：3946－3961．

［9］Li H，Yang Y．The adaptive finite elementmethod based on multi-scale discretizations for eigenvalue problems［J］．Computers＆Mathematics with Applications，2013，65（7）：1086－1102．

［10］劉進明，應懷樵．FFT譜連續(xù)細化分析的富里葉變換法［J］．振動工程學報，1995，8（2）：162－166．

LIU Jin-ming，YING Huai-qiao．Zoom FFT spectrum by Fourier transform［J］．Journal of Vibration Engineering，1995，8（2）：162－166．

［11］易立強，鄺繼順．一種基于FFT的實時諧波分析算法［J］．電力系統(tǒng)及其自動化學報，2007，19（2）：98－102．

YILi-qiang，KUANG Ji-shun．FFT-based algorithm for realtime harmonic analysis［J］．Proceedings of The CSU－EPSA，2007，19（2）：98－102．

［12］Schwarz B，Richardson M．Experimentalmodal analysis［J］．CSIReliability Week，1999，36（6）：30－32．

［13］周杰梁，鄭百林．基于中心差分法的纖維結構碰撞動力學分析［J］．力學季刊，2011，32（3）：466－472．

ZHOU Jie-liang，ZHENG Bai-lin．Impact dynamics analysisof fiber-composed structure based on central difference method［J］．Chinese Quarterly of Mechanics，2011，32（3）：466－472．

［14］鐘萬勰．子域精細積分及偏微分方程數值解［J］．計算結構力學及其應用，1995，12（3）：253－260．

ZHONG Wan-xie．Subdomain precise integration method and numerical solution of partial differential equations［J］．Computational Structural Mechanics and Applications，1995，12（3）：253－260．［15］李秀梅，吳鋒，張克實．結構動力微分方程的一種高精度攝動解［J］．工程力學，2013，30（5）：8－12．

LI Xiu-mei，WU Feng，ZHANG Ke-shi．A high-precision perturbation solution for structural dynam ic differential equations［J］．Engineering Mechanics，2013，30（5）：8－12．

［16］李秀梅，吳鋒，黃哲華．結構動力方程一種新的級數形式的解析解［J］．振動與沖擊，2010，29（6）：219－222．

LIXiu-mei，WU Feng，HUANG Zhe-hua．New series form of analytical solution for the structural dynamic equations［J］．Journal of Vibration and Shock，2010，29（6）：219－222．

［17］李秀梅，吳鋒，黃哲華．PCG法的理論解釋及在結構分析中的應用［J］．廣西大學學報（自然科學版），2009（4）：463－468．

LI Xiu-mei，WU Feng，HUANG Zhe-hua．Theoretical explanation for PCG method and its application in structural analysis［J］．Journal of Guangxi University（Natural Science Edition），2009（4）：463－468．

Fast Fourier transform method for generalized eigenvalue problems

WU Feng，XU Xiao-ming，ZHONGWan-xie
（State Key Laboratory of Structural Analysis of Industrial Equipment，Department of Engineering Mechanics，Dalian University of Technology，Dalian 116024，China）

For generalized eigenvalue problems，a method based on the fast Fourier transform（FFT）was developed．In the proposed method，the dynamical structural response was viewed as a signal containing all information about the vibrational frequencies．Using FFT to the signal，the vibrational frequencies can be obtained．The method is a kind of direct solution method which can compute all eigenvalues without the matrix inversion．A numerical example manifests the correctness of the proposed method．

eigenvalue；dynamical structural response；fast Fourier transform；sampling

O327；O214．6

：A

10．13465／j．cnki．jvs．2014．22．012

國家基礎研究發(fā)展計劃973（2009CB918501）

2013－08－26 修改稿收到日期：2013－11－15

吳鋒男，博士生，1985年生