文立超，張應(yīng)紅，劉文龍，張奧申

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

基于導(dǎo)波的無損檢測方法具有傳播速度快、方向性好、應(yīng)用范圍廣等特性。但是由于導(dǎo)波具有多模態(tài)、易頻散的特性，其頻率的選擇及頻散往往會影響檢測結(jié)果的準(zhǔn)確性，因而研究管道中導(dǎo)波的頻散關(guān)系對于導(dǎo)波檢測的應(yīng)用具有重要意義。

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

1 理論基礎(chǔ)

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

(1)

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

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

(2)

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

(3)

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

(4)

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

(5)

橫波波速cT為

(6)

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

(7)

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

根據(jù)變分原理，亥姆霍茲方程相應(yīng)的泛函可用下式表示。

(8)

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

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

AU=0

(9)

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

(10)

式中:Ne為網(wǎng)格數(shù)。

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

2 頻散分析

2.1 頻散關(guān)系求解

圖1 管道二維離散模型

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

表1 鋁的材料參數(shù)

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

圖2 頻率-波數(shù)曲線

進(jìn)一步地，由已求得的頻率f及波數(shù)可求得相速度cp，其滿足關(guān)系式

(11)

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

圖3 相速度頻散曲線

類似地，根據(jù)所求的頻率波數(shù)關(guān)系，可知其群速度cg滿足以下表達(dá)式

(12)

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

圖4 群速度頻散曲線

對圖2～4進(jìn)行分析，可以得到以下結(jié)論。

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

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

(3) 0～20 kHz的低頻段與頻率超過20 kHz的頻段相比，導(dǎo)波模態(tài)較少，且頻率和波數(shù)之間的關(guān)系近似為線性關(guān)系。

2.2 頻散關(guān)系分析

(13)

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

根據(jù)參考文獻(xiàn)[11]可得剪切系數(shù)如下式所示。

(14)

式中:m為內(nèi)外徑之比，即b/a。

對式(13)進(jìn)行傅里葉變換，得到其解。對應(yīng)的實(shí)數(shù)域下彎曲波數(shù)為

(15)

將ω以及表1中的相關(guān)材料參數(shù)代入上式即可求得彎曲模態(tài)下的理論解 。

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

表2 不同方法得到的結(jié)果比較

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

3 模態(tài)特征分析

模態(tài)特征可以表征導(dǎo)波在管道中傳播的波結(jié)構(gòu)特征，通常選取的模態(tài)特征有位移變化、應(yīng)力變化以及能量變化。由于導(dǎo)波檢測中較易獲得位移量[14]，所以文章采用位移變化來表征模態(tài)特征，以頻率為10 010 Hz時(shí)的各模態(tài)位移特征進(jìn)行模態(tài)特征分析。為使得位移的表示更加合理，將笛卡爾坐標(biāo)系下的位移變化量轉(zhuǎn)換為柱坐標(biāo)系下的位移變化量，則各個(gè)方向的位移如下式所示。

(16)

式中:ux,uy,uz分別為模態(tài)在直角坐標(biāo)系中的x,y,z三個(gè)方向上的位移分量；ur,uθ,uz分別為模態(tài)在柱坐標(biāo)系中的徑向、切向及軸向上的位移分量；θ為旋轉(zhuǎn)角。

通過上面的轉(zhuǎn)換矩陣，在后處理結(jié)果中得到柱坐標(biāo)下各個(gè)分量的位移，并對位移分量進(jìn)行歸一化處理，從而得到各模態(tài)下的位移分布，如圖5～7所示。

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

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

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

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

4 結(jié)論

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