分治法在管道渦流檢測阻抗解析中的應用

王志剛，羅清旺，師奕兵

（電子科技大學自動化工程學院 成都 611731）

·自動化技術·

分治法在管道渦流檢測阻抗解析中的應用

王志剛，羅清旺，師奕兵

（電子科技大學自動化工程學院 成都 611731）

介紹了一種鐵磁性管道渦流阻抗模型的數值解析方法。該方法利用阻抗模型中貝塞爾函數在大宗量（分治法產生的一個子區間）時的漸近性，在整個廣義積分區間采用分治法簡化積分函數，降低了對阻抗模型廣義積分的計算難度?；趯τ嬎懔颗c計算準確度的折中，討論了該解析方法中分治點的選擇原則。將該方法的阻抗模型計算結果與通過物理檢測設備測試實際管道的值進行比較，驗證了分治法解阻抗模型的可行性。該方法對于解析柱坐標系下管道渦流檢測阻抗模型的應用，具有簡單、快速和高精度的優點。

貝塞爾函數； 大宗量； 分治法； 阻抗模型； 數值解析

在管道渦流無損檢測研究中，常利用Maxwell方程及某種邊界條件，列出電磁波在鐵磁性管道中傳播的模式方程，然后根據該模式方程，研究電磁波在鐵磁性管道中的傳播阻抗特性，并從傳播阻抗特性中提取出管道物理信息［1-3］，最后基于這些物理信息判別鐵磁性管道的損傷情況。上述針對鐵磁性管道的檢測方法，涉及到對電磁場中傳播阻抗的數值解析。由文獻［4-5］可知，管道渦流阻抗模型是關于修正貝塞爾函數的一種復雜廣義積分形式，該積分形式決定了它難于直接利用貝塞爾函數的積分性質進行計算。雖然也有一些關于貝塞爾廣義積分函數的研究性論文發表，但是它們多是針對特定的貝塞爾函數形式［6-9］，如較常見且采用的有復雜遞推法［7］以及針對特定形式的貝塞爾函數積分［8］方法。前者只能用于特定的遞推形式，后者應用于有限的積分區間，文獻［9］提出了修正貝塞爾函數的數值積分方法的數學演繹。對于本文所需解析的渦流阻抗模型中具有的函數復雜性、函數高速振蕩衰減及在廣義區間積分等特征，它們均不適用。

另一方面，函數的漸近性被廣泛應用于解析復雜積分函數，并具有很好的解析效果［10-13］?？梢岳秘惾麪柡瘮禎u進性，采用分治法處理復雜積分函數。分治法的思想是將一個難以直接解決的大問題，分割成一些規模較小的相同問題，以便分而治之。依此將阻抗模型的廣義積分區間分治為小宗量積分區間和大宗量積分區間，并在大宗量積分區間應用貝塞爾函數的漸近性，選定逼近函數簡化阻抗函數。一般的數列排序分治法［14］，通常將最后的數值點作為分治點，它的選取不影響計算結果?？紤]到需對渦流阻抗模型中積分區間［0，∞］分治，所以不能取最后數點作為分治點，這樣對分治點的選取時需要考慮逼近函數對原函數的逼近精度，以及阻抗函數的數值計算準確度。

本文以分治法求解阻抗模型的過程包含：分治積分區間及阻抗模型積分函數的衰減特性和逼近特性分析；積分函數計算和結果準確度分析兩個階段。依靠分治法和貝塞爾函數逼近的特性的簡化，可在保證計算結果精度的條件下大大減少計算量。

1 鐵磁性管道渦流阻抗模型

圖1所示，激勵及接收線圈均處于管道內部（管道無限長）。a、b、c分別為激勵、接收線圈、管道的內半徑，線圈軸向距離為s，s

圖1 管道電磁傳播阻抗模型

在柱坐標系中，各區域中任一點的磁位量設為A（r，z；r′，z′）其為包含復變量的貝塞爾方程：

式中，I1（xnr）,K1（xnr）分別為1階第一類和第二類修正貝塞爾函數；Cn（x）、Dn（x）分別代表在區域n的函數系數； xn為介質傳播參數，可表示為，

式中，下標n=1，2，3，4分別對應區域I，II，III，IV，區域I和II為空氣；x=x1=x2,μ0=μ1=μ2。對式（1）求解，并設a=b＜c，經復雜推導可得鐵磁性管道電磁渦流檢測阻抗模型為［15］：

其中，式（3）為空氣中線圈組傳輸阻抗，式（4）表示由于鐵磁性管道存在所導致阻抗函數變化的增量。因此，阻抗函數可寫為：Z=Zd+Λ。由式（3）、式（4）的積分特性可知，被積函數是衰減函數，且計算難度較大。另外，由管道趨膚特性可知，|x3c|＞＞1，根據修正貝塞爾函數的性質，式（6）成立：

將式（6）代入式（5）簡化，有：

由于管道具有較大的電導率，可對式（2）化簡，得到：

觀察式（4）～式（8），由于貝塞爾函數特性及廣義的積分區間，化簡之后仍然很難直接積分計算或利用貝塞爾函數性質進行計算。

2 分治法解阻抗模型

分治法解阻抗模型，是將廣義積分區間分治為小宗量和大宗量區間；在大宗量區間，以簡化近似函數代替原函數進行計算；而在小宗量區間，則進行定積分或求和計算。

2.1 貝塞爾函數的Melin-Barnes逼近

由分治法思想可知，在積分區間采用分治法，若要對積分函數化簡，其在大宗量區間時必須有高精度逼近的近似函數。由文獻［13］可知，漢克函數在大宗量時具有良好的逼近函數，并且，修正貝塞爾函數可由漢克函數線性表示。所以，如果已知漢克函數在大宗量時的高精度逼近函數，就可以通過兩類貝塞爾函數的線性關系，求得其在大宗量時修正貝塞爾函數的高精度逼近式，該線性關系為：

Borel求和法和Melin-Barnes數值法對貝塞爾函數均可實現良好逼近［13］，而兩者比較，Borel求和法計算量大，本文采用Melin-Barnes數值法對大宗量區間的修正貝塞爾函數進行逼近。

由于柯西留數定理在一些特殊實積分，如反常積分、廣義積分是很好的分析手段，而Melin-Barnes數值法的基本思想是：在形如式（11）的一般冪級數中應用柯西留數定理，有：

在一定條件下［13］，在式（11）中應用柯西留數定理，可將 S（N，z）表示成Melin-Barnes積分，有：

式中，N-1＜c＜N 。由于漢克函數的完全等式是形如式（11）所示的無窮級數的冪級數，不利于直接數值計算，所以應用柯西留數定理將漢克函數的冪級數形式轉換成形如式（12）的Melin-Barnes積分，得到漢克函數的高精度逼近式，以便于達到快速數值計算漢克函數的目地，具體過程參見文獻［13］。

在大宗量（λ＞＞1）時，可由Melin-Barnes數值法［13］逼近的漢克函數簡化為：

式中，α表示階數。式（13）、式（14）提供了大宗量區間精確的漢克函數逼近公式，再結合式（9）、式（10）可以得到修正貝塞爾函數在大宗量區間內高精度的逼近式。然而，這種逼近不可避免存在一定誤差，它與分治點的選取有直接關系。

2.2 選取分治點

在分治法解析阻抗模型的過程中，分治點λ的選取決定了漸近函數引入起點。理論上，分治點λ越大，大宗量積分區間越小，計算結果越準確，但是計算量就會越大；反之，分治點λ越小，計算量越小，大宗量區間的逼近誤差就會越大。合理選取分治點，使得計算量和結果準確性得以折中。

為了便于分治點的選取，首先需要對阻抗模型被積函數的變化趨勢進行分析，它是一個穩定系統的輸出（該穩定系統的輸出是一個電壓信號）。由穩定系統特性可知，其輸出信號必然收斂，阻抗模型在［0，∞）內廣義積分被積函數必然衰減，這為大宗量區間內被積函數的逼近提供了可能，函數的衰減特性導致大宗量區間內的積分量占積分總量較少。圖2所示為阻抗Zd、Λ的被積函數在有限區間內的分布情況。

由圖2可知， Zd的被積函數呈震蕩衰減的趨勢，Λ的被積函數的實部、虛部均呈快速衰減的趨勢，其值在λ＞4時基本趨于0。阻抗模型的被積函數變化趨勢為分治算法提供了可行性，亦為分治算法中分治點的起點提供一個大致范圍，如λZd≥10，λΛ≥4。為獲得誤差小的漸近貝塞爾函數分治起點，結合文獻［16］，修正貝塞爾函數漸近曲線的誤差為：

式中，α為修正貝塞爾函數的階數；λ為分治點。由式（15）可知λ值越大，修正貝塞爾函數的漸近誤差越小，Melin-Barnes數值法的逼近度越高，但是會增加數值計算量。

圖2 渦流阻抗模型的被積函數分布

分治點的選取通過觀察函數計算結果的穩定度并進一步確定阻抗模型數值結果的準確度為原則。結合圖2和式（15）就可為分治算法提供誤差小，數值大小合適的分治點。

2.3 分治法解渦流阻抗模型

在阻抗模型中，物理參數滿足如下關系式，0

實際測量中通常取 a=b，而a＞b和a＜b情況下的阻抗函數衰減速度均大于 a=b情況下函數的衰減速度。因此，若該方法對 a=b的情況適用，那么對a＞b或a＜b也適用。在式（3）中，選取分治點 λ0分治積分區間，在大宗量區間進行漸近性分析，將式（9）、式（10）、式（13）、式（14）代入式（3）的被積函數，當xa=xb≥λ0，a=b＜c時，有：

所以，式（3）可化為：

同理，在式（4）中選取分治點 λ1，在大宗量區間進行漸近分析。當1xaλ≥時，有：

將式（18）代入式（4），得：

由式（1）、式（2）可知，dx表示磁位量 A（r， z；r′， z′）在介質傳播點x的微積分， x3可由式（8）化簡。式（17）所示的積分式在數值計算軟件中可直接采用函數命令調用的方式得到數值結果。而式（19）在積分區間［0，λ1/a］的表達式仍很難直接函數調用得到積分結果。由圖2可知，Λ的被積函數衰減很快，故在［0，λ1/a］分治區間內，可用積分定義將式（19）變化為：

式（20）中，在分治區間［0，λ1/a］中累加和的步進xΔ在一定程度上影響結果的精度，本文選取Δx=0.005，以保證計算結果的準確性。利用MATLAB軟件函數庫及高效的數值計算和符號計算功能，對于可直接計算的貝塞爾函數數值積分，采用函數命令調用的形式實現，如式（17）。對于難于積分的復雜貝塞爾函數形式，采用積分定義的方法積分，如式（19）在［0，λ1/a］的積分。計算中，約定在計算量和計算結果精度滿足要求前提下，當3次結果相對誤差低于0.5%時，即采用第2次的計算結果，其對應的截斷點λ判為合適。算法過程如下：

1） 參數初始化。包括：c ＞ a = b ＞ 0，ω ＞ 0，μ0＞ 0，σ3＞0,ε3＞0,λ=1,λstep=2,x=0，xstep=0.005等。

2） 若x≤λ / a ，λ=λ+λstep

3 結果驗證

為了驗證分治法解析阻抗模型的正確性，選取合適的管道參數值，進行數值計算。依實際情況取參數：a=b=0.018 9 m，s=0.063 5m，c=0.067 2 m ，f =20kHz ，μ0=4π ×10-7H/m，μ3=80×4π ×10-7H/m，σ3=5.95×106S/m。根據圖2和式（15）?。害?≥10， λ1≥4 。由分治點λ與參量的關系式可知，在大宗量區域：λ≤xa。應用分治法計算Zd和Λ的數值結果如表1、表2所示。

表1 λ0不同時　Zd的值

表2 λ1不同時Λ的值

由表1、表2可知，在分治點較小時（如 λ0=1，λ1=1時）Zd與Λ的分治結果與后面穩定的分治結果有一定的差距，正如前文分析，高精度的貝塞爾逼近函數須在大宗量區間引入，即分治點較大時引入。當分治點較小時引入，會造成分治結果不正確。而在 λ0≥10和 λ1≥4的范圍，Zd與Λ的分治結果基本趨于穩定，其值的穩定性體現了應用分治法解阻抗模型的正確性。由表1可知 λ0=15作為Zd的分治點；由表2知 λ1=5作為Λ的分治點最合適，既減小了計算量，又保證了結果值的準確度。為了體現分治點對積分結果的影響（因為在 λ0≥10的范圍， Zd的值趨于穩定，呈現水平的特性，所以分治點對分治結果的影響主要體現在0≤λ0≤10范圍內），選取 Zd～λ0（λ0≤12）繪圖，如圖3所示。

由圖3可知，分治點值越大，Zd積分結果越穩定。在分治點 λ0≤8的區域內，其對Zd積分結果影響很大， λ0越小，影響效果越大。表1、表2和圖3的積分結果的準確性是通過積分結果漸近誤差的連續穩定性來保證的，即

圖3 Zd積分值與分治點取值關系

為了更好地驗證分治法解阻抗模型的正確性，本文通過實際管道阻抗檢測值來驗證其一致性。

對于式（7）所示的算子，在電磁渦流近場的阻抗模型中反映了管道帶來的耦合影響，對其可做如下變形，以獲得更加簡潔的表述：

式中，

首先，運用分治法計算獲得式（3）、式（4）的L值，再通過管道測試試驗，獲得實際渦流阻抗的L值，將上述結果與去耦度ζ做變化曲線進行比較，曲線如圖4所示。圖4從電學試驗角度展示了分治算法解阻抗模型的可行性和正確性。由圖4可知運用分治算法解電磁渦流近場阻抗模型的結果與實際物理檢測值具有很好的一致性（個別實測異常點除外，如2c=0.153 7 m，ζ=13.1579時）。另外，由于分治算法的快速、資源占用少的特點，其在解阻抗模型中更具有可操作性。

圖4 分治法L值與實測L值比較圖

4 結 論

通過對鐵磁性管道近場渦流傳輸阻抗模型的分析，利用分治法將該阻抗模型的廣義積分區間分治為小宗量積分區間和大宗量積分區間，并在大宗量積分區間應用貝塞爾函數的漸近性，選定漢克函數逼近簡化阻抗函數，完成復雜阻抗函數的解析，同時，通過計算結果的逼近誤差穩定性來確定結果的準確性。電學對比實驗的結果表明，該方法不僅可以快速有效的解析輸出，且與實際測試結果具有很好的一致性，充分證明了該方法的可行性。

［1］ DODD C V， DEEDS W E. Analytical solutions to eddy -current probe-coil problems［J］. J Appl Phys， 1968， 39（6）：2829-2838.

［2］ DODD C V， CHENG C C， DEEDS W E. Induction coils coaxial with an arbitrary number of cylindrical conductors［J］. J Appl Phys， 1974， 45（2）： 638-647.

［3］ DESJARDAINS D P R， KRAUSE T W， GAUTHIER N. Analytical modeling of the transient response of a coil encircling a ferromagnetic conducting rod in pulsed eddy current testing［J］. NDT and E International， 2013， 60：127-131.

［4］ VASIC D， BILAS V. Model-based stochastic inversion of coil impedance for determination of tube inner radius and electromagnetic properties［J］. AUTOMATIKA： Journal for Control， Measurement， Electronics， Computing and Communications， 2009， 50（3）： 227-233.

［5］ VASIC D， BILAS V. Validation of a coil impedance model for simultaneous measurement of electromagnetic properties and inner diameter of a conductive tube［J］. IEEE Trans Instrum Meas， 2006， 55（5）： 337-341.

［6］ CHEN Ru-yun， AN Cong-pei. On the evaluation of infinite integrals involving Bessel functions［J］. Applied Mathematics and Computation， 2014， 235： 212-220.

［7］ LIN Q G. Infinite integrals involving Bessel functions by contour integration［J］. Integral Transforms and Special Functions， 2013， 24（10）： 783-795.

［8］ CHEN Ru-yun， AN Cong-pei. On evaluation of Bessel transforms with oscillatory and algebraic singular integrands［J］. Journal of Computational and Applied Mathematics， 2014， 264： 71-81.

［9］ BARICZ A， POGANY T. On a sum of modified Bessel functions［J］. Mediterranean Journal of Mathematics， 2014，11（2）： 349-360.

［10］ NICULESCU C， POPOVICI F. The asymptotic behavior of integrable functions［J］. Real Analysis Exchange， 2012，38（1）： 157-168.

［11］ XIANG Shu-huang. Laplace transforms for approximation of highly oscillatory Volterra integral equations of the first kind［J］. Applied Mathematics and Computation， 2014， 232：944-954.

［12］ GIULIANO R， MACCI C. Large deviation principles for sequences of maxima and minima［J］. Communications in Statistics-Theory and Methods， 2014， 43（6）： 1077-1098.

［13］ KOWALENKO V. Exactification of the asymptotics for Bessel and Hankel functions［J］. Applied Mathematics and Computation， 2002， 133（2）： 487-518.

［14］ CORMEN T H， LEISERSON C E， RIVEST R L， et al. Introduction to algorithms［M］. 3rd ed. Beijing： China Machine Press， 2013： 16-21.

［15］ 王志剛， 羅清旺， 師奕兵. 鐵磁性管道內渦流線圈耦合分析與管道參數檢測［J］. 儀器儀表學報， 2014， l35（12）：2843-2851. WANG Zhi-gang， LUO Qing-wang， SHI Yi-bing. Analysis of eddy current coil coupling in ferromagnetic pipe and parameter detection［J］. Chinese Journal of Scientific Instrument， 35（12）： 2843-2851.

［16］ 仇海亮. 復宗量貝塞爾函數的數值解法及在測井解釋中的應用［J］. 中國石油和化工標準與質量， 2013（17）： 128. QIU Hai-liang. Numerical solution of complex argument Bessel function and its application in the log interpretation［J］. China Petroleum and Chemical Standard and Quality， 2013（17）： 128.

編 輯 漆 蓉

Divide-and-Conquer PAippp e rEodadchy ACpuprrlieendt tToe Ismtpedance Analysis of

WANG Zhi-gang， LUO Qing-wang， and SHI Yi-bing
（School of Automation engineering， University of Electronic Science and Technology of China Chengdu 611731）

This paper proposes a unique analytical method for eddy current impedance model of ferromagnetic pipe. Based on the asymptotic behavior of Bessel function in model under bulk quantity occasion，the method exploits the divide-and-conquer approach in the generalized integral interval to simplify model calculation. The calculation results of the model are compared with the results of practical pipe eddy current test，verifying the correctness and feasibility of calculation on the divide-and-conquer approach. The analytical method has the advantages of conciseness and high precision for impedance model applications under cylindrical coordinate systems.

Bessel function； bulk quantity； divide-and-conquer approach； impedance model；numerical analysis

TE1

A

10.3969/j.issn.1001-0548.2016.02.011

2014 - 07 - 21；

2015 - 11 - 25

“十二五”國家科技重大專項（2011ZX05020-006-005）

王志剛（1967 - ），男，副教授，主要從事檢測技術方面的研究.