一種管道中的導波頻散計算方法

文立超，張應紅，劉文龍，張奧申

(桂林電子科技大學 機電工程學院, 桂林 541004)

基于導波的無損檢測方法具有傳播速度快、方向性好、應用范圍廣等特性。但是由于導波具有多模態、易頻散的特性，其頻率的選擇及頻散往往會影響檢測結果的準確性，因而研究管道中導波的頻散關系對于導波檢測的應用具有重要意義。

英國帝國理工大學的LOWE等[2-3]研究了充水埋地鋼管的頻散曲線及衰減特性，并基于解析法開發了頻散曲線軟件Disperse；SECO[4]在此方法基礎上通過MATLAB軟件開發了PC disp工具箱，該工具箱可用于分析管狀導波的頻散及傳播情況。北京工業大學的何存富[5]采用矩陣法對立柱埋地部分和未埋地部分的頻散特性進行了研究。近年來半解析有限元方法在頻散分析方面得到了廣泛應用。其中，具有代表性的有:BOCCHINI[6]基于半解析有限元法開發了GUWGUI軟件，可用于各種復雜截面的導波頻散關系的求解;浙江大學的胡劍虹[7]等對半解析有限元法的求解矩陣規模進行了改進，減少了計算量;海軍工程大學的王悅民[8-9]采用有限元的特征頻率法求解了管道的頻散曲線，該方法建立了三維模型，求解自由度大，需從可視化的角度統計周期數及模態數，進而求得頻散關系，但操作繁瑣。筆者提出一種基于有限元的模式分析法，通過推導導波傳播的波動方程，采用分離變量法得到亥姆霍茲方程，并用有限元法求得特征值解，從理論上證明了模式分析法的可行性。利用COMSOL Multiphysics軟件建立了二維管道截面的有限元模型。根據低頻下的導波頻散關系特點，進一步采用彈性波理論及鐵木辛柯梁理論對模式分析方法及半解析有限元法的求解結果進行驗證，最后對位移分量進行坐標變換及波結構特征分析。該方法與傳統半解析有限元法相比，需要設置模式階數，在有限元軟件中可以方便地得到頻散關系，為導波的頻散計算分析提供了一種新的思路。

1 理論基礎

假定管道材料為各向同性、均勻同質的線彈性體，內外徑分別為a,b，不考慮體力及慣性力的影響時，其一般的彈性動力學方程(Navier-stokes方程)如式(1)所示。

(1)

式中:u為時間諧振位移矢量;ρ為材料密度;λ和μ為材料的拉梅常量。

位移場通過亥姆霍茲分解可表示為標量φ的梯度和矢量H的旋度，即

(2)

將式(2)代入式(1)中可得

(3)

若要使式(3)成立，需要滿足下列關系式

(4)

由式(4)分析可得此方程為兩個形式一致的波動方程，而縱波波速cL可表示為

(5)

橫波波速cT為

(6)

因此，在波導中不僅存在以縱波傳播的伸縮擾動，而且存在以橫波傳播的切變擾動。模式分析法是對波動方程進行時空分離，并求解模式特征值解的方法。式(4)中由于求解的是同一類型的方程，所以可將φ,H統一記做ψ，cL、cT統一記做c，初始狀態下管道內外壁節點的位移及速度為0。往z方向添加入射的簡諧波，波的傳播函數可以用ψ=U0exp(iωt-ikz)進行描述,其中U0為波的幅值矢量。于是對ψ作變量分離ψ=U(x,y,z)T(t)=U(x,y,z)exp(iωt)，可得到下面兩個微分方程式。

(7)

式中：k為分離常數波數;U為波動位移幅值；t為波在波導里的傳播時間；T為僅與t有關的待定函數；x,y,z為空間坐標的笛卡爾坐標系;ω為波傳播的角頻率。

根據變分原理，亥姆霍茲方程相應的泛函可用下式表示。

(8)

式中:Ω為求解區域;J(U)為波動位移幅值的能量，min[J(U)]表示最小能量。

將其剖分為四邊形網格，設Ni為形函數，記B=Ni。則通過對泛函求極值可導出有限元方程,如式(9)所示。

AU=0

(9)

式中:A為總體剛度矩陣，其表達式如式(10)所述。

(10)

式中:Ne為網格數。

若給定參數f(波的特征頻率)，有k=2πf/c，對式(9)中矩陣A進行特征值求解，可得到一組k，將k代入式(9)，可求出U，即模態特征。

2 頻散分析

2.1 頻散關系求解

圖1 管道二維離散模型

以外徑a為30 mm，內徑b為25 mm的鋁材管道為例，建立幾何模型，添加材料為鋁，并劃分單元網格為四邊形的映射網格，管道二維離散模型如圖1所示。由于是理想情況下的管道，其不受外力作用，所以設置為自由邊界條件。鋁主要的材料參數如表1所示(表中E為彈性模量，ν為泊松比)。

表1 鋁的材料參數

將建立的幾何模型放置在固體力學場中求解，設定求解步驟為模式分析，并取求解模式數為30，選擇頻率步長為100 Hz，模式搜索基準值為2πf/c，這里根據經驗值設定c為330 m·s-1，并以頻率作為掃描參數。根據文獻[7]，通常導波檢測所用的檢測頻率較常規超聲波的頻率低，由于管道截面中的導波高階模態十分復雜，計算會引入較多的復雜模態，且對計算要求較高，因此選定求解頻率的最大值為100 kHz，通過求解得到了頻率-波數曲線，如圖2所示。為了驗證文章方法的合理性，將其和半解析有限元方法求解的結果進行了比較。

圖2 頻率-波數曲線

進一步地，由已求得的頻率f及波數可求得相速度cp，其滿足關系式

(11)

根據式(11)可得到其相速度頻散曲線，如圖3所示。

圖3 相速度頻散曲線

類似地，根據所求的頻率波數關系，可知其群速度cg滿足以下表達式

(12)

為求得群速度，可使Δf=1，則需要另計算f′=f+1的波數值，以得到兩組k值，求得Δk。所求得的群速度頻散曲線如圖4所示。

圖4 群速度頻散曲線

對圖2～4進行分析，可以得到以下結論。

(1) 基于有限元的模式分析方法不但與半解析有限元法所得到的結果相吻合，而且求解得到了由結構特征變化衍生出來的環狀模態[11-12](兩種方法所繪制曲線的未重合部分)，證明了該法的有效性。環狀模態是由于模式疊加產生的，在實際檢測中應該避免該類模態產生，以降低檢測信號的復雜性。

(2) 導波在管中傳播存在多模態現象，高頻情況下此現象更加突出。在0～100 kHz頻段中兩種方法所繪制曲線的重合部分，導波模態主要有縱向模態、扭轉模態及彎曲模態三類，每一類模態又細分為多種模式。而隨著頻率的增加，模式數相應增加。

(3) 0～20 kHz的低頻段與頻率超過20 kHz的頻段相比，導波模態較少，且頻率和波數之間的關系近似為線性關系。

2.2 頻散關系分析

(13)

式中:A為截面面積;I為截面慣性矩;ωb為彎曲位移;θb為轉角;κ為剪切系數；G為剪切模量。

根據參考文獻[11]可得剪切系數如下式所示。

(14)

式中:m為內外徑之比，即b/a。

對式(13)進行傅里葉變換，得到其解。對應的實數域下彎曲波數為

(15)

將ω以及表1中的相關材料參數代入上式即可求得彎曲模態下的理論解 。

綜上所述，可得頻率為10 010 Hz下的波數理論值，且可通過計算求得與其他兩種方法的相對誤差，如表2所示。

表2 不同方法得到的結果比較

由表2可得到，由理論所求得的波數解與半解析有限元法、模式分析法所求的波數基本一致，均可得到準確解。但是，針對L(0,1)和F(0,1)模態的波數解，模式分析法較半解析有限元法的精度有微弱的提高，而扭轉波模態T(0,1)波數的求解精度顯著提高，可見模式分析法更適用于扭轉模態的求解。綜合來看，半解析有限元法與模式分析法所得結果的誤差均較小，但是有限元的模式分析方法的精度更高，體現了該方法的優越性。

3 模態特征分析

模態特征可以表征導波在管道中傳播的波結構特征，通常選取的模態特征有位移變化、應力變化以及能量變化。由于導波檢測中較易獲得位移量[14]，所以文章采用位移變化來表征模態特征，以頻率為10 010 Hz時的各模態位移特征進行模態特征分析。為使得位移的表示更加合理，將笛卡爾坐標系下的位移變化量轉換為柱坐標系下的位移變化量，則各個方向的位移如下式所示。

(16)

式中:ux,uy,uz分別為模態在直角坐標系中的x,y,z三個方向上的位移分量；ur,uθ,uz分別為模態在柱坐標系中的徑向、切向及軸向上的位移分量；θ為旋轉角。

通過上面的轉換矩陣，在后處理結果中得到柱坐標下各個分量的位移，并對位移分量進行歸一化處理，從而得到各模態下的位移分布，如圖5～7所示。

圖5 10 010 Hz下L(0,1)模態的各分量位移圖

圖6 10 010Hz下T(0,1)模態的各分量位移圖

圖7 10 010 Hz下F(1,1)模態的各分量位移圖

由圖5可知，L(0,1)模態存在軸向位移，內外壁的位移分布幾乎相等，且遠大于徑向位移。而徑向位移和周向位移量接近于0。由圖6可看出，對于T(0,1)模態的導波，存在周向位移，且外表面壁的位移比內表面壁的位移更大。而F(1,1)模態在三個分量下均有位移，周向位移和軸向位移都較小，徑向位移最大，在管道中的傳播比較復雜。

4 結論

采用基于有限元的模式分析法求解管道中導波的頻散關系，通過亥姆霍茲方程得到波數和角頻率之間的關系式，并采用有限元方法進行了數值計算，利用COMSOL Multiphysics軟件建立模型，給定頻率參數，模式數及模式搜索基準值求解得到了頻散曲線及管道截面上的位移形變結果。與半解析有限元方法的求解結果相比，模式分析法得到了由模式疊加產生的環狀模態，說明了該法求解的全面性。同時，參照理論解，模式分析法的精度良好，且更適用于T(0,1)模態的求解。