一種求解蘭姆波特征方程的改進方法

何呂龍，尚柏林，張亞豪

（空軍工程大學航空航天工程學院，西安710038）

一種求解蘭姆波特征方程的改進方法

何呂龍，尚柏林，張亞豪

（空軍工程大學航空航天工程學院，西安710038）

針對Rayleigh-Lamb波動方程求解過程中存在分母為零和正切函數(shù)間斷點的問題，對原有求解方法和解的數(shù)據(jù)存儲結構進行了改進，并使用該改進方法繪制了單層鋁合金板中蘭姆波的相速度和群速度頻散曲線。通過與實驗結果進行對比，驗證了方法的合理性和正確性。仿真分析結果表明，改進的求解方法在不增加其他約束的情況下，減小了數(shù)學分析的難度，降低了程序設計的復雜程度，節(jié)約了程序運行時間；改進的數(shù)據(jù)存儲結構可以根據(jù)工程需要方便地繪制指定模式的頻散曲線，對基于蘭姆波的主動結構健康監(jiān)測系統(tǒng)研究中激勵信號中心頻率的選取也有一定的指導作用。

蘭姆波，頻散曲線，健康監(jiān)測，損傷定位

0　引言

傳統(tǒng)的結構健康監(jiān)測系統(tǒng)（Structural Health Monitoring System，SHMS）大多基于被動無損檢測/評估技術（NDT/NDE），需要大量專業(yè)儀器。由于這些儀器大多體積較大，結構復雜，重量較大，一般不能實現(xiàn)對結構健康狀態(tài)的實時在線監(jiān)測，尤其是針對飛機、飛船等對結構重量要求嚴格的航空、航天飛行器。因此，發(fā)展一種體積小、重量輕的在線結構健康監(jiān)測系統(tǒng)是航空航天技術的發(fā)展趨勢和要求。目前，基于蘭姆波（Lamb Waves）的主動結構健康監(jiān)測被認為是一種最有潛力的在線結構健康監(jiān)測技術。蘭姆波在結構中傳播時，結構內部的各種損傷都會引起應力集中，導致裂紋擴展加速。這些損傷以及損傷周圍的邊界都會引起在結構中傳播的蘭姆波信號的散射和能量吸收，通過對傳感器接收到的損傷信號與健康信號進行對比分析，便可對結構的健康狀況做出評定。正是基于這種現(xiàn)象，蘭姆波可以被用來對結構中的損傷進行檢測和定位。與基于聲發(fā)射的結構健康監(jiān)測原理類似，基于蘭姆波的損傷定位也是通過不同位置的傳感器接收到的損傷反射信號與健康信號的時間延遲進行損傷定位的［1］。由于蘭姆波具有多模式和頻散特性，使得對損傷信號的分析處理變得困難。因此，想要有效利用蘭姆波進行損傷定位，必須首先解決蘭姆波的頻散曲線繪制問題。

蘭姆波頻散曲線的求解主要有數(shù)值計算和試驗兩種方法，國內外學者對此也做了大量的研究工作。劉增華等［2］借鑒文獻［3］中求解多層板結構頻散方程的方法，提出采用子波法對自由板的Navier運動方程進行求解來繪制蘭姆波頻散曲線。徐西寧等在文獻［4-5］中分別采用半解析有限元法和圖解法得到了蘭姆波的頻散曲線。數(shù)值計算法具有成本低，速度快，精度可控等優(yōu)點，是求解蘭姆波頻散方程的常用方法。由于繪制完整的滿足工程需要的蘭姆波頻散曲線需要大量的數(shù)據(jù)點，采用實驗方法周期長，成本高，誤差大，一般只能繪制低階模式的頻散曲線，因此，試驗方法主要用于檢驗確認數(shù)值計算結果的準確性。雖然以上方法都可以得到蘭姆波的頻散曲線，但是這些方法都較為復雜，很難被一般工程技術人員理解應用。解的數(shù)據(jù)存儲結構也不夠合理，不能根據(jù)實際需要繪制指定模式的頻散曲線。

針對以上問題，本文對原有求解方法和解的數(shù)據(jù)存儲結構進行了改進，使之更適合一般工程技術人員使用。通過與實驗結果進行對比，驗證了本文方法的合理性和正確性，對基于蘭姆波的主動結構健康監(jiān)測系統(tǒng)研究中激勵信號中心頻率的選取也有一定的指導作用。

1　蘭姆波特征方程分析

1.1Rayleigh-Lamb波動方程

蘭姆波，也即板波，是在作為波導的具有表面邊界的板中傳播的彈性應力波［6］。由于蘭姆波是介質中橫波和縱波遇到板的上、下表面時發(fā)生波形轉換形成的不同波形相互疊加生成的，因此，蘭姆波具有多模式和頻散特性。在均質薄板中，根據(jù)質點振動位移特點，蘭姆波可分為對稱和反對稱兩種模式，通常用S0，A0，S1，A1…Sn，An表示，其中Sn表示對稱模式，An表示放對稱模式，其波動特性采用Rayleigh-Lamb波動方程描述如式（1）、式（2）所示［7］。

對稱模式：

反對稱模式：

式中，

可以看出，式（1）、式（2）是兩個復雜的超越方程。k0與ω的關系是非線性的，而且對不同的模式也對應不同的非線性關系，這也從數(shù)學角度證明了蘭姆波的多模式和頻散特性。因此，對蘭姆波在板結構中的傳播特性的完整描述通常是一系列分散曲線形式，作為頻厚積的函數(shù)來表示板結構中蘭姆波傳播模式的速度，每條曲線代表一種特定的模式，這就是蘭姆波的頻散曲線［8］。

蘭姆波在材料中傳播時會發(fā)生頻散，是一種多模式復合波，相速度即是指各個分波的傳播速度，用cp表示。不同模式、相位相同的振動會相互疊加，這種相同相位振動的傳播速度就是波群傳播的速度，即蘭姆波的群速度，用cg表示。群速度可以通過相速度來求解，相速度與群速度的關系如式（5）所示：

將ω=2πf帶入上式得到：

式中，d為薄板厚度，fd表示頻率與厚度的乘積，為了便于處理，通?？勺鳛橐粋€參數(shù)處理，稱為頻厚積。頻厚積的單位為赫茲·米（Hz·m），但大多數(shù)情況，所研究的薄板厚度都在毫米量級，因此，通常用兆赫茲·毫米（MHz·mm）作為fd單位。

1.2Rayleigh-Lamb波動方程分析

觀察式（1）、式（2）不難發(fā)現(xiàn)，隨著cp取值不同，方程兩邊都可能會出現(xiàn)復數(shù)、分母為零、正切函數(shù)間斷點等問題。在考慮實際物理意義繪制頻散曲線時，只需在實數(shù)范圍內求解方程。本文以對稱模式方程為例進行分析，反對稱模式同理可得。

將式（3）、式（4）分別帶入式（1），得到：

對于均勻薄板，有式（9）成立:

式中，μ為泊松比，ρ為密度。容易證明cl>cs恒成立，因此，在方程求解時，只需對cp進行分段處理。

①0＜cp＜cs

當cp在（0，cs）內取值時，由式（8）可得：

化簡可得：

②cs＜cp＜cl

當cp在（cs，cl）內取值時，由式（8）可得：

方程兩邊同乘復數(shù)單位i，化簡可得：

③cp＞cl

當cp在（cl，+∞）內取值時，原方程不變，即：

2　改進算法

2.1改進求解方法

觀察式（10）～式（12）并結合正切函數(shù)的性質，不難發(fā)現(xiàn)：

和

都是方程分母的零點。

和

都是方程左邊的正切函數(shù)的間斷點。

這兩個問題都是數(shù)值求解蘭姆波頻散曲線時必須考慮的，否則會導致求解結果不嚴謹，甚至出現(xiàn)錯誤。

式中，系數(shù)A、B、C、D隨cp的取值范圍不同而變化。

這種處理方法雖然避免了分母為零的情況，通過合理設置迭代步長，也得到了蘭姆波頻散曲線，但是并沒有消除正切函數(shù)間斷點的影響，求解過程并不嚴謹，必須合理設置求解步長，以避免正切函數(shù)出現(xiàn)間斷點。

文獻［10］中根據(jù)正切函數(shù)的性質，首先利用間斷點進一步將求解區(qū)間劃分為若干子區(qū)間，然后再在每一個子區(qū)間內進行求解。但是，由于蘭姆波的多模式和頻散特性，使用這種方法需要經(jīng)過復雜的數(shù)學分析，討論函數(shù)間斷點，得出多個求解區(qū)間，并且區(qū)間長度隨模式不同而變化，使編程求解變得更加困難，不適合一般工程技術人員使用。如果數(shù)學分析時漏掉某一區(qū)間，還可能導致繪制的頻散曲線不完整。

針對上述問題，本文利用三角函數(shù)的性質，將式（10）～式（12）進一步轉化為如下形式：

式中，系數(shù)A、B、C、D隨cp的取值范圍不同而變化。

這種處理方法同時避免了分母為零和正切函數(shù)存在間斷點的問題。編程求解時，不再需要經(jīng)過復雜的數(shù)學分析確定分段點和求解子區(qū)間，除了實際物理意義，也沒有增加其他約束條件，使得數(shù)學分析和程序設計的難度都明顯降低。

在繪制群速度頻散曲線時，由于數(shù)值計算得到的相速度頻散曲線是一系列離散點構成，無法求導。如果將這些點進行曲線擬合，不僅需要大量的時間，還會帶來較大的誤差。為了得到群速度頻散曲線，從導數(shù)的數(shù)學定義出發(fā)，采用如下所示近似方法。

顯然，當Δcp、Δfd足夠小時，這種近似方法是合理的。

2.2改進解的存儲結構

文獻［8-10］都通過數(shù)值計算方法得到了蘭姆波的頻散曲線，但是在對這些數(shù)值解的存儲時都采用簡單的數(shù)組存儲，這使得必須一次繪制出所有的數(shù)據(jù)點才能得出各個模式蘭姆波的頻散曲線，但是在實際工程應用中，人們往往只需要特定模式的頻散曲線，因此，這種存儲方式不能滿足工程需要。

為解決這個問題，本文將數(shù)據(jù)存入矩陣C中：

其中，C的第1列為頻厚積fdn，從C的第2列開始，每一列都是一種模式的蘭姆波在相應頻厚積下的速度（相速度或群速度）。以頻厚積fdn為橫坐標，相速度為縱坐標，把特定一列的速度值繪制出來就可以得出相應模式的頻散曲線了，這樣就不必再繪制出所有模式的頻散曲線或者使用其他數(shù)據(jù)歸類算法，使頻散曲線的繪制更加簡便、靈活。圖1為通過這種方法選擇繪制的對稱模式S0和S1。

圖1　對稱模式頻散曲線

通過這種數(shù)據(jù)存儲結構，與式（13）的近似方法相結合，可以非常方便地由相速度頻散曲線得到群速度頻散曲線。如果采用文獻［9-10］中簡單數(shù)組存儲方式，則這種近似方法就不能使用了。

3　改進效果分析

3.1算例分析

為了驗證本文所用方法的合理性和改進效果，以飛機上大量使用的鋁合金板（cl=6 450 m/s，cs=3 090 m/s）為例，通過仿真計算，繪制了鋁合金板中的相速度和群速度頻散曲線。

將ω=2πf帶入式（10）～式（12），三式均是cp和fd的函數(shù)。在程序設計時，有兩種循環(huán)方式可以選擇，一種是以cp為外層循環(huán)，另一種是以fd為外層循環(huán)。如圖1所示，若任取cp，對fd由小到大取值，當cp=cp1時，所求得的第1個解對應對稱模式S0，當cp=cp2時，求得的第1個解對應對稱模式S1，兩者對應不同模式。反之，若任取fd，對cp由小到大取值，求得的第1個解始終對應對稱模式S0，第2個解始終對應對稱模式S1，以此類推。因此，為方便數(shù)據(jù)存儲和曲線繪制，選擇以fd為外層循環(huán)，cp為內層循環(huán)，并利用根的存在判定定理和二分法快速求解方程的根，將最終的計算結果存入矩陣C中。方程解的個數(shù)可以通過步長Δcp、Δfd來控制。為節(jié)約程序運行時間，程序中設置了最大迭代次數(shù)和求解精度兩個迭代終止條件，當?shù)螖?shù)或者求解精度達到程序設定值，迭代便停止，進入下一步。完整的程序框圖如下頁圖2所示。

3.2效果分析

文獻［9］中的處理方法雖然避免了分母為零的情況，但是沒有消除正切函數(shù)間斷點的影響。在編程求解時，對迭代步長較為敏感，需要合理選擇步長，迭代前需判定是否為間斷點。

文獻［10］利用間斷點進一步將求解區(qū)間劃分為若干子區(qū)間，在每一個子區(qū)間內進行求解。但是，由于蘭姆波的多模式和頻散特性，這種方法需要經(jīng)過復雜的數(shù)學分析，得出多個求解區(qū)間，并且區(qū)間長度隨模式變化而變化，使編程求解變得困難，如果分析錯誤還可能導致繪制的頻散曲線不完整。

圖2　程序設計框圖

本文提出的改進方法與文獻［8-10］相比，不僅同時避免了分母為零和正切函數(shù)存在間斷點的問題，而且在編程求解之前，也不需要經(jīng)過復雜的數(shù)學分析確定分段點和求解子區(qū)間，除了實際物理意義，沒有其他約束條件，使得數(shù)學分析和程序設計的難度都明顯降低。

在程序運行效率方面，通過編程驗證，在相同運行環(huán)境和求解精度下，文獻［9］中所用算法的運行時間為8.45 s，文獻［10］中所用算法的運行時間為12.36 s。雖然改進后的方法方程的長度增加，但是方程中不含tan和tanh，程序運行時間較文獻［9］所述方法運行時間減少，降至7.33 s。與文獻［10］所述方法相比，因為求解前不需要計算間斷點和求解區(qū)間，程序運行時間減少更加明顯。

通過以上分析可知，與原方法相比，改進后的方法數(shù)學分析難度和程序設計復雜性都得到降低，程序運行時間也得到減少。最終繪制的完整的相速度和群速度頻散曲線分別如圖3（a）、圖3（b）所示。圖中，S、A分別表示對稱和反對稱模式，下標0、1、2…分別表示對應模式階數(shù)。

4　結束語

根據(jù)Rayleigh-Lamb波動方程的特點，對其求解方法和解的數(shù)據(jù)存儲結構進行了改進，消除了原方程存在的分母為零和正切函數(shù)不連續(xù)帶來的問題，使頻散曲線繪制更加簡便、靈活，通過實例分析，取得了如下結論：

圖3　蘭姆波頻散曲線

①通過與參考文獻［8-10］結果進行對比，可知本文提出的改進算法是合理可行的，計算結果是正確無誤的，能夠正確求解并繪出蘭姆波頻散曲線，且本文使用的方法在數(shù)學分析和編程時都更為簡單；②觀察圖3（a）可知，除對稱模式S0和反對稱模式A0外，其他模式的蘭姆波都存在比較明顯的截止頻率，頻率低于截止頻率的蘭姆波將不能在介質中傳播。這對于主動結構健康監(jiān)測系統(tǒng)中選擇激勵信號中心頻率，減少多模式信號干擾具有指導作用；③觀察圖3（b）可知，不同模式蘭姆波群速度隨頻厚積的變化規(guī)律不同。應將激勵信號中心頻率選在頻散曲線隨頻厚積變化平緩的區(qū)間，可以減小損傷定位的誤差，提高損傷定位精度和對細小裂紋的識別能力。從圖中看，當頻厚積大于1 MHz·mm時，反對稱模式A0群速度幾乎不再變化，不考慮其他因素時是最理想的激勵信號。

［1］袁慎芳.結構健康監(jiān)控［M］.北京：國防工業(yè)出版社，2007: 221.

［2］劉增華，何存富，吳斌.單層板頻散曲線的矩陣表示及Matlab實現(xiàn)［J］.無損檢測，2005，27（5）:225-227.

［3］Lowe M J S.Matrix Techniques for Modeling Ultrasonic Waves in Multilayered Media［J］.IEEE Transactions on Ultrasonics，F(xiàn)erroelectrics and Frequency Control，1995，42（4）:525-542.

［4］徐西寧，余祖俊，朱力強，等.半解析有限元法分析蘭姆波頻散特性［J］.儀器儀表學報，2013，34（2）：247-253.

［5］徐西寧，余祖俊，朱力強.圖解法求解Lamb波頻散方程［J］.電子測量與儀器學報，2012，26（11）：966-971.

［6］張海燕，他得安，劉鎮(zhèn)清.層狀各向異性復合板中的蘭姆波［M］.北京:科學出版社，2008:3.

［7］閻石，張海鳳，蒙彥宇.Lamb波頻散曲線的數(shù)值計算及實驗驗證［J］.華中科技大學學報（城市科學版），2010，27（1）:1-4.

［8］邱雷.基于壓電陣列的飛機結構健康監(jiān)測與管理系統(tǒng)研究［D］.南京:南京航空航天大學，2011.

［9］艾春安，李劍.蘭姆波頻率方程的數(shù)值解法［J］.無損檢測，2005，27（6）:294-296.

［10］鄭祥明，趙玉珍，史耀武.蘭姆波頻散曲線的計算［J］.無損檢測，2003，25（2）:66-68.

An Improved Method for Solving Lamb Frequency Dispersion Equation

HE Lü-long，SHANG Bo-lin，ZHANG Ya-hao
（School of Aeronautics and Astronautics Engineering，Air Force Engineering University，Xi'an 710038，China）

In order to avoid the zero denominator and existing of discontinuous points in solving Rayleigh-Lamb wave equation，an advanced calculation method and data saving structure are promoted. The frequency dispersion curves of phase velocity and group velocity in an monolayer aluminum plate are plotted respectively，contrastive analysis shows the results are in good agreement with experimental data，which verifies the rationality and accuracy of the improved method.The analysis result demonstrates the improved method can decrease the difficulty of mathematics analysis，lower the complexity of program design，save the calculating time with no additional constraint.The new data saving structure can plot the specified mode of dispersion curves easily，which helps to choose the central frequency of drive signal in an active structural health monitoring system based on Lamb wave.

lamb waves，dispersion curves，structural healthy monitoring，damage location

TG115.28

A

1002-0640（2015）08-0101-05

2014-07-05

2014-08-06

何呂龍（1990-），男，四川綿陽人，碩士研究生。研究方向：飛機結構健康監(jiān)控。