楊旭輝, 余厚全, 陳強(qiáng), 程峰

(長江大學(xué), 湖北 荊州 434023)

0 引 言

Lamb波是最常見的一種平面導(dǎo)波形式。對Lamb波在缺陷板中反射、散射、模式改變等現(xiàn)象的研究,為超聲無損檢測提供了重要的理論依據(jù)。但這些研究多集中于自由薄板的情況,即僅考慮薄板處于真空或空氣中且薄板為完全彈性介質(zhì),從而不需要考慮能量損失,也就是不需要考慮幅值衰減。目前這種自由薄板中的Lamb波從模型的建立到頻散曲線的繪制已有比較成熟的方法[1-4]。

隨著工程應(yīng)用的發(fā)展,Lamb波被引入套管井固井水泥膠結(jié)質(zhì)量的評價中[5]。若滿足一定條件,Lamb波在圓柱層狀套管井中的傳播可近似為在層狀介質(zhì)中的傳播。Lamb波在上述介質(zhì)中向前傳播時由于能量損失會出現(xiàn)幅值衰減的現(xiàn)象。這種波稱為泄漏Lamb波,層狀介質(zhì)中泄漏Lamb波的情況更為復(fù)雜。不僅要考慮傳播過程中的頻散性還要考慮幅值衰減。實(shí)踐表明,對這種多層介質(zhì)建模所得到的泄漏Lamb波的模型比自由薄板中的Lamb波模型要復(fù)雜得多。若依舊采用自由薄板中Lamb波的數(shù)值求解方法將會出現(xiàn)求解困難甚至完全失效的情況。因此,為在固井質(zhì)量評價中更加有效地應(yīng)用泄漏Lamb波,本文研究了層狀介質(zhì)中泄漏Lamb波的建模以及頻散方程的求解方法,繪制了泄漏Lamb波頻散及衰減曲線,以期為泄漏Lamb波在工程中的應(yīng)用提供一種實(shí)用的數(shù)值計算方法。

1 層狀介質(zhì)泄漏Lamb波的建模

1950年Thomson[6]首次采用傳遞矩陣方法對任意層狀介質(zhì)中Lamb波傳播進(jìn)行建模,隨后Haskell[7]進(jìn)一步發(fā)展了該方法,并使之在地震學(xué)中得到了應(yīng)用。在此基礎(chǔ)上,文獻(xiàn)[8-10]給出了針對黏彈性層狀介質(zhì)以及浸液薄板中泄漏Lamb波的研究思路。目前大部分文獻(xiàn)對于層狀介質(zhì)中泄漏Lamb波傳播的建模仍沿用傳遞矩陣方法。

圖1 層狀介質(zhì)模型

圖1為泄漏Lamb波傳播的層狀介質(zhì)模型?？蓪⒃撃Ｐ椭械目v波和橫波勢函數(shù)分別設(shè)為

(1)

(2)

u(n)=φ(n)+×ψ(n)

(3)

(4)

根據(jù)Snell定理以及波數(shù)的定義,可將D(n)化簡為僅與各層參數(shù)及x1方向的波數(shù)有關(guān)的矩陣。其元素可參考文獻(xiàn)[13]。式(4)適用于圖1層狀介質(zhì)模型中的所有層以及層內(nèi)的任意位置。將式(4)分別應(yīng)用于各層的頂部和底部,并考慮到位移及應(yīng)力的連續(xù)性,可得到關(guān)系式

(5)

式中,Lm,top為層m的頂部;[L]L2為各層的傳遞矩陣[13]。以下將矩陣積[L]L(m-1)…[L]L3·[L]L2記為 [S]。與自由薄板中的Lamb波不同,為了得到泄漏Lamb波的模型,必須考慮圖1中模型的上下半空間為非真空的情況?？紤]式(4)、式(5)可得到

(6)

(7)

可得

(8)

f(K)=T22·T44-T42·T24=0

(9)

K為前面提到的x1方向的波數(shù)。當(dāng)考慮波在傳播過程中能量損失時,常用的處理方法是引入復(fù)頻率或波數(shù)。若引入復(fù)波數(shù),則復(fù)波數(shù)的實(shí)部為頻率與相速度的比值,虛部對應(yīng)泄漏Lamb波在x1方向傳播時的衰減。求取泄漏Lamb波的頻散曲線和衰減曲線即為使頻率、相速度和衰減值同時滿足式(9),也即求解復(fù)頻散方程f(K)的0點(diǎn)。

該泄漏Lamb波的頻散方程是復(fù)超越方程,不存在解析形式的根。必須通過數(shù)值方法進(jìn)行求解,而且該方程的復(fù)雜性隨著層數(shù)的增加而增加。因此,對數(shù)值計算方法的效率提出了更高的要求。

2 常用方法及分析

與泄漏Lamb波相關(guān)的大部分文獻(xiàn)和專著中多采用文獻(xiàn)[13]提出的方法求取復(fù)頻散方程f(K)的0點(diǎn)。其核心步驟:粗略尋根、最優(yōu)尋根、曲線追蹤。①粗略尋根如圖2所示,產(chǎn)生如圖3所示的頻散方程f(K)的幅值變化曲線。當(dāng)頻散方程的幅值處于局部極小值時,將對應(yīng)一相速度,該相速度和固定頻率即構(gòu)成頻散曲線上的一近似點(diǎn)。②最優(yōu)尋根主要思路是保持某一參量不變而另2項參量以較小步長交替搜索,找到頻散方程幅值新的極小點(diǎn)。當(dāng)極小點(diǎn)處頻散方程的幅值與0值的差滿足精度要求時,停止搜索,該極小點(diǎn)對應(yīng)的頻率、相速度、衰減值即可作為頻散方程的最優(yōu)根。③曲線追蹤時首先需確定2至3個初始最優(yōu)根,根據(jù)這些最優(yōu)根進(jìn)行外推以獲取新的最優(yōu)根的預(yù)測值,得到預(yù)測值后對預(yù)測值再進(jìn)行最優(yōu)尋根。最終可以得到整個選定模式的頻散曲線或衰減曲線。

圖2 相速度掃描示意圖

圖3 頻散方程幅值變化曲線示意圖

圖4 頻散方程幅值變化曲線示意圖

對上述算法進(jìn)行驗(yàn)證。在求解自由或浸水鋼板Lamb波頻散曲線和衰減曲線時可得到滿意的結(jié)果。當(dāng)鋼板的上下介質(zhì)聲參數(shù)為非對稱時,模型如表1所示。在粗略尋根后進(jìn)行最優(yōu)尋根時,部分極小點(diǎn)對應(yīng)的頻散方程的幅值沒有出現(xiàn)趨近于0值的現(xiàn)象,視作無法收斂。為分析原因,固定一個頻率進(jìn)行相速度掃描。得到了如圖4的幅值隨相速度變化的曲線。顯然曲線中的局部極小值點(diǎn)多于圖3中的局部極小值點(diǎn)。進(jìn)一步驗(yàn)證,局部極小值點(diǎn)A、D在最優(yōu)尋根過程中,頻散方程的幅值將逐漸趨近于0,視作可以收斂。局部極小值點(diǎn)B、C進(jìn)行最優(yōu)尋根時無法收斂,將導(dǎo)致根的取舍的抉擇難度增加,從而增加誤判。這種誤判是致命的,將直接導(dǎo)致第3步曲線追蹤時出現(xiàn)大量錯誤,以至無法繪制與頻散曲線對應(yīng)的衰減曲線,而且也使得大量時間耗費(fèi)在無法收斂點(diǎn)的計算上。圖5是保留所有誤判時a0模式的頻散曲線圖,其對應(yīng)參數(shù)為表1?？捎^察到a0模式出現(xiàn)大量間斷點(diǎn),明顯不符合模式波的傳播規(guī)律。表1中將水的橫波速度設(shè)為較小值(相比縱波速度小2個數(shù)量級)。這種處理方式在含流體層的層狀介質(zhì)中較為常見,可避免傳播矩陣推導(dǎo)過程中出現(xiàn)矩陣奇異的現(xiàn)象?？梢宰C明,這樣做流體層不必特別處理就能得到正確的結(jié)果。

表1 層狀介質(zhì)模型參數(shù)

圖5 保留誤判時a0模式頻散曲線示意圖

3 改進(jìn)算法及計算結(jié)果

由前面的分析,計算泄漏Lamb波頻散方程f(K)的0點(diǎn)是獲取頻散曲線和衰減曲線的關(guān)鍵,其實(shí)質(zhì)是求解復(fù)超越方程f(K)的復(fù)根。文獻(xiàn)[12]中Dayal和Kinral應(yīng)用二分法和牛頓-弦截法分別求出了頻散方程復(fù)根的實(shí)部和虛部。但在2條頻散曲線的交叉點(diǎn)處,求根有可能出現(xiàn)異常[13]。其他方法還有圍數(shù)積分法、Muller法、Newton-Raphson法等[14]。圍數(shù)積分法利用了復(fù)變函數(shù)的廣義對數(shù)留數(shù)定理[15],可以穩(wěn)定地給出高精度復(fù)數(shù)根,但由于在解的搜索過程中需要對頻散方程進(jìn)行求導(dǎo)和積分計算,當(dāng)層數(shù)較多時,耗時較長。Muller法也稱作拋物線法,同樣可求解頻散方程的復(fù)根,但需要要構(gòu)造二次插值多項式導(dǎo)致計算略顯繁瑣。Newton-Raphson法思路簡潔適合編程實(shí)現(xiàn)且在單根附近具有較高的收斂速度,但對初始近似值的精度要求較高,如果初始近似值的誤差較大,則可能給出發(fā)散的結(jié)果。

綜合考慮后的算法在粗略尋根和最優(yōu)尋根方面做了改進(jìn),曲線追蹤方法不變。改進(jìn)后粗略尋根的目的是盡可能以較小的誤差得到到頻散方程0點(diǎn)的初始近似值,然后在最優(yōu)尋根中利用初始近似值和Newton-Raphson方法進(jìn)行迭代計算得到0點(diǎn)的最優(yōu)值。搜索初始0點(diǎn)近似值時,首先仍然是采用前面提到的文獻(xiàn)[13]的方法,得到如圖4所示的曲線。然后在圖4中的局部極小值點(diǎn)及其附近進(jìn)行局部峰的搜索,最后判斷局部峰內(nèi)是否確定存在0點(diǎn)。若存在0點(diǎn),則局部峰所處位置即為所求0點(diǎn)的近似值。

圖6為某一固定頻率點(diǎn)形成的局部峰形狀示意圖。圖6中水平軸分別對應(yīng)離散化了的相速度和衰減。離散間距需足夠小,以使得局部峰凸顯出來,垂直軸對應(yīng)H(K)的幅值。在固定頻率處相速度和衰減同時參與了掃描,形成了有關(guān)H(K)幅值的三維圖,避免了圖4中粗略尋根時對Lamb波衰減預(yù)估中的人為因素。

圖6 頻散方程在某固定頻率點(diǎn)處的局部峰形狀圖

得到H(K)的幅值,通過直接搜索的方法確定局部峰值點(diǎn)。圖7對應(yīng)的是某一固定頻率處離散的相速度—衰減平面示意圖。圖7中矩形所選區(qū)域?qū)?yīng)的H(K)已離散化,離散化值為H(m,n),m,n=1,2,3,…。若H(m,n)為該區(qū)域的最大值,則H(m,n)即為所求局部峰。相應(yīng)的頻散方程f(K)在該處就可能存在0點(diǎn)。

為進(jìn)一步降低誤判率,可引入函數(shù)

(10)

該函數(shù)用來量化離開局部峰頂點(diǎn)時H(m,n)的下降速度。下降速度越快,局部峰處存在0點(diǎn)的可能性越大。因此可將Γ>ε作為判斷標(biāo)準(zhǔn),ε為小于1的正數(shù)需根據(jù)經(jīng)驗(yàn)和實(shí)際情況選取。當(dāng)計算得到的Γ滿足Γ>ε判定為局部峰,即可能存在0點(diǎn)的區(qū)域,可用于下一步繼續(xù)分析,否則認(rèn)為是無0點(diǎn)存在的普通小突起。

利用上述方法獲取局部峰,即得到了頻散方程f(K)存在0點(diǎn)的可能區(qū)間。再判斷在已求得的可能區(qū)間是否確定存在0點(diǎn)。由復(fù)變函數(shù)理論可知,如果f(K)在簡單閉曲線C上與C內(nèi)解析,且在C上不等于0,那么f(K)在C內(nèi)0點(diǎn)的個數(shù)等于1/2π乘以當(dāng)K沿C的正向(逆時針)繞行一周f(K)的輻角的改變量?？杀硎緸?/p>

(11)

式中,N表示0點(diǎn)個數(shù);ΔC+Argf(K)表示K沿C正向繞行1周時f(K)輻角的改變量。在圖7中如果H(m,n)已確認(rèn)為局部峰,則利用輻角原理判斷0點(diǎn)個數(shù)時繞行路徑C即指H(m,n)處的虛線方框。如圖7中箭頭所示,顯然它是簡單閉曲線。由頻散方程f(K)知道,在路徑C內(nèi)無極點(diǎn)存在,且對于實(shí)際的物理系統(tǒng)總認(rèn)為是連續(xù)的可導(dǎo)的。因此,只考慮f(K)在C內(nèi)與C上是解析的且在C上不為0(C上為0則該點(diǎn)即為所求0點(diǎn))的情況。

圖7 相速度—衰減平面示意圖

計算機(jī)在計算輻角改變量之前需要先對路徑離散化,離散化后相鄰2點(diǎn)輻角的改變量表示為

Δi=Δ[Ki,Ki+1]Argf(K)i=1,2,3,…M

(12)

式中,M為路徑C上總的離散點(diǎn)數(shù),且Km+1與K1點(diǎn)重合。那么由繞行路徑C計算得到的0點(diǎn)個數(shù)為

(13)

同樣,由復(fù)變函數(shù)理論知道,若路徑C在離散化過程中滿足|Δi|≤π,i=1,2,3,…M,則

Δi=Arg [f(Ki+1)/f(Ki)]i=1,2,3,…M

(14)

文獻(xiàn)[16]中給出了詳細(xì)的針對繞行路徑C的離散化方法,應(yīng)用該方法可以使離散化的Δi滿足|Δi|≤π。文獻(xiàn)[16]還給出了完整的基于FORTRAN的源代碼。

將式(14)帶入式(13)后可得

(15)

式(15)用于判斷局部峰內(nèi)是否確定存在0點(diǎn)。如果存在多個0點(diǎn),可以將圖7中離散點(diǎn)的間距進(jìn)一步減小,總可以得到只包含1個0點(diǎn)的局部峰。

當(dāng)確定局部峰內(nèi)包含0點(diǎn),可以采用文獻(xiàn)[13]的最優(yōu)尋根方法計算最優(yōu)根。但這種交替式的搜索方式收斂速度較慢?？紤]到此時可以取離局部峰最近的離散點(diǎn)作為0點(diǎn)的近似值,Newton-Raphson方法可以利用該近似值作為初始近似值進(jìn)行迭代計算。此時,迭代公式為

(16)

式中,x0、y0、x1、y1為實(shí)數(shù);K0、K1為不相等的初值,均選自局部峰附近的離散點(diǎn),且滿足|f(K1)|<|f(K0)|。具體計算步驟

(1) 局部峰附近選定2個初始值K0、K1,計算相應(yīng)的函數(shù)值f(K0)、f(K1);

(2) 代入式(12)的迭代公式計算Ki+1的值;

(3) 如果出現(xiàn)|Ki+1-Ki|<ξ的情況,則終止迭代計算,Ki+1即為所求,否則轉(zhuǎn)(4),ξ為給定精度值;

(4) 如果迭代計算超過最大指定次數(shù)N,則認(rèn)為迭代過程不收斂,計算失敗,否則轉(zhuǎn)(2)繼續(xù)進(jìn)行迭代計算。

通過以上步驟可以順利得到頻散方程f(K)0點(diǎn)的最優(yōu)根,為曲線追蹤提供有力保障。多次計算后發(fā)現(xiàn),在求取局部峰的過程中,當(dāng)ε選擇合適時,局部峰內(nèi)無0點(diǎn)存在的幾率很小。所以,對于局部峰內(nèi)是否存在0點(diǎn)的判斷只在頻散曲線出現(xiàn)異常時引入。例如曲線追蹤過程中在預(yù)測值附近并未搜索到局部峰的存在,原因可能是ε設(shè)置得過大或是頻散曲線波動過大。此時可適當(dāng)調(diào)低ε設(shè)置值并擴(kuò)大搜索局部峰的范圍。這種情況可能導(dǎo)致多個局部峰的出現(xiàn),此時可引入0點(diǎn)判斷的算法,排除無0點(diǎn)存在的普通小突起的干擾。采用這種可剝離式的算法可大大提高計算速度。

圖8 鋼板位于慢水泥與水之間時的a0模式頻散曲線圖

圖9 鋼板位于慢水泥與水之間時的a0模式衰減曲線圖

圖10 鋼板位于快水泥與水之間時的a0模式頻散曲線圖

圖11 鋼板位于快水泥與水之間時的a0模式衰減曲線圖

圖8、圖9采用改進(jìn)的算法分別得到了鋼板位于慢水泥與水之間時a0模式的頻散曲線和衰減曲線。參數(shù)見表1。由曲線可知相速度被慢水泥的縱波速度截為2支,即圖8中在200 kHz附近出現(xiàn)了2個模式波。圖10、圖11中鋼板位于快水泥與水之間,此時快水泥的縱波速度大于a0模式的最大相速度,橫波速度仍然將衰減曲線截為2支,但相速度所受影響并不明顯。圖8至圖10與文獻(xiàn)[17]中利用商業(yè)軟件DISPERSE所計算曲線完全一致。

4 結(jié) 論

(1) 泄漏Lamb波頻散和衰減曲線的繪制在工程應(yīng)用中具有重要意義。已有文獻(xiàn)的相關(guān)算法多適用于自由或浸水的層狀介質(zhì)。當(dāng)層狀介質(zhì)上下半無限介質(zhì)的聲學(xué)參數(shù)為非對稱時,在尋找最優(yōu)根的過程中,出現(xiàn)誤判的幾率增大。

(2) 對層狀介質(zhì)中泄漏Lamb波進(jìn)行建模并對相應(yīng)的頻散方程進(jìn)行分析,設(shè)計了一種針對該頻散方程的數(shù)值計算方法。該方法在局部峰附近利用復(fù)變函數(shù)理論,對根的存在性予以明確的判斷,從而極大地降低了最優(yōu)根的誤判率。

(3) 利用該方法可準(zhǔn)確計算并繪制層狀介質(zhì)中泄漏Lamb波頻散和衰減曲線,具備一定的通用性。

參考文獻(xiàn):

[1] 曹正敏, 林莉, 李喜孟. 蘭姆波頻散曲線的繪制與試驗(yàn)驗(yàn)證 [J]. 理化檢驗(yàn): 物理分冊, 2008, 44(9): 193-199.

[2] 艾春安, 李劍. 蘭姆波頻率方程的數(shù)值解法 [J]. 無損檢測, 2005, 27(6): 294-296.

[3] 艾春安, 吳安法. 蘭姆波頻散方程的分析及數(shù)值迭代算法 [J]. 上海航天, 2008, 28(5): 42-44.

[4] 閻石, 張海鳳, 蒙彥宇. Lamb波頻散曲線的數(shù)值計算及試驗(yàn)驗(yàn)證 [J]. 華中科技大學(xué)學(xué)報: 城市科學(xué)版, 2010, 27(1): 1-4.

[5] 王華, 陶果, 尚學(xué)峰, 等. 輕質(zhì)水泥固井質(zhì)量聲波測井評述與方法研究 [J]. 測井技術(shù), 2014, 38(2): 165-173.

[6] Thomson W T. Transmission of Elastic Waves Through a Strified Solid [J]. J Appl Phys, 1950, 21: 89-93.

[7] Haskell N A. Dispersion of Surface Waves in Multilayered Media [J]. Bull Seim Soc Am, 1953, 43: 17-34.

[8] Nayfeh A H, Chimenti D E. Elastic Wave Propagation Influidloaded Multiaxial Anisotropic Media [J]. Acoust Soc Am, 1991, 89: 542-549.

[9] Cervenka P, Challande P. A New Efficient Algorithm to Compute the Exact Reflection and Transmission Factors for Plane Waves in Layered Absorbing Media (Liquids and Solids) [J]. Acoust Soc Am, 1991, 89: 1579-1589.

[10] Hosten B. Bulk Heterogeneous Plane Waves Propagation Through Viscoelastic Plates and Stratified Media with Large Values of Frequency Domain [J]. Ultrason., 1991, 29: 445-450.

[11] Chimenti D E, Nayfeh A H, Butler D L. Leaky Rayleigh Waves on a Layered Halfspace [J]. J Appl Phys, 1982, 53: 170-176.

[12] Dayal V, Kinra V K. Leaky Lamb Waves in an Anisotropic Plate I: An Exact Solution and Experiments [J]. J Acoust Soc Am, 1989, 85: 2268-2276.

[13] Lowe M J S. Plate Waves for the NDT of Diffusion Bonded Titanium [D]. London: University of London, 1993.

[14] 易大義, 沈云寶, 李有法. 計算方法 [M]. 2版. 杭州: 浙江大學(xué)出版社, 2003: 150-154.

[15] 張忠誠, 王成. 對數(shù)留數(shù)定理的推廣及應(yīng)用 [J]. 洛陽大學(xué)學(xué)報, 2006, 21(2): 20-22.

[16] Xingren Ying, Katz I N. A Reliable Argument Principle Algorithm to Find the Number of Zeros of an Analytic Function in a Bounded Domain [J]. Numer Math, 1988, 53: 143-163.

[17] Froelich B. Multimode Evaluation of Cement Behind Steel Pipe [J]. J Acoust Soc Am, 2008, 123: 3648.