張思全，胡盛斌，陸文華

(上海工程技術(shù)大學(xué)航空運(yùn)輸學(xué)院，上海 201620)

在石油、化工、航空等領(lǐng)域，許多關(guān)鍵設(shè)備和部件由多層導(dǎo)電材料制成，對其內(nèi)部缺陷和結(jié)構(gòu)完整性的檢測是一個較困難的問題。超聲檢測法雖具有可檢測較厚材料的優(yōu)勢，但不適于檢測相鄰層間含有空氣的多層導(dǎo)電結(jié)構(gòu)。電渦流檢測技術(shù)利用材料中缺陷與渦電流、線圈激勵電磁場之間相互作用引起線圈阻抗變化的原理來檢測缺陷并對結(jié)構(gòu)完整性進(jìn)行評價，并且具有不需耦合劑、可檢測多層導(dǎo)電結(jié)構(gòu)隱藏缺陷的優(yōu)勢［1－4］。

如何快速、精確求解線圈阻抗是渦流檢測的重要問題。數(shù)十年來，許多學(xué)者對此進(jìn)行了研究。Dodd和Deeds研究了正弦電流激勵線圈位于有均勻厚度導(dǎo)電材料覆層的半無限大導(dǎo)電平板上方的情況，得到了積分形式的矢量磁位閉合解析解［5］;其后，Luquire等人又研究了線圈置于任意層導(dǎo)電結(jié)構(gòu)中間時的阻抗變化模型［6］。在國內(nèi)，雷銀照研究了半無限大導(dǎo)體上方覆蓋一定厚度金屬涂層構(gòu)成的兩層平板導(dǎo)電結(jié)構(gòu)的軸對稱渦流解析解［7］;黃平捷研究了置于多層導(dǎo)電結(jié)構(gòu)上方線圈響應(yīng)與多層導(dǎo)電結(jié)構(gòu)厚度變化之間的關(guān)系模型［8］;范孟豹推導(dǎo)了位于任意層導(dǎo)電結(jié)構(gòu)上方渦流探頭的阻抗解析模型，在阻抗計算中引入符號運(yùn)算法求解矢量磁位表達(dá)式系數(shù)，減小了程序計算量［9］。但這些研究采用徑向無限大導(dǎo)體，所得解為積分形式，存在計算困難、時間花費(fèi)大的問題。

本文在前人研究基礎(chǔ)上，采用解域截取法分析矩形截面圓柱線圈置于任意層導(dǎo)電平板上方的正向問題，推導(dǎo)渦流檢測探頭線圈與任意層導(dǎo)電結(jié)構(gòu)相互作用的級數(shù)形式阻抗變化計算模型。為驗(yàn)證模型的正確性，針對單層導(dǎo)電結(jié)構(gòu)的特例，計算了線圈置于其上方時的阻抗變化，分析了影響計算精度的主要因素，并與有限元計算結(jié)果進(jìn)行了對比。

1 任意層導(dǎo)電板上方線圈阻抗計算

1.1 任意層導(dǎo)電板上方單匝線圈阻抗計算

1.1.1 任意層導(dǎo)電板上方單匝線圈的矢量磁位

如圖1a)所示，一個單匝通電圓形線圈置于任意層(設(shè)為n層)導(dǎo)電平板上方，其軸線與任意層平板表面垂直;設(shè)各層平板的厚度、電導(dǎo)率、介電常數(shù)都是任意的，各層介質(zhì)都是線性、各向同性、均勻的。設(shè)線圈中為正弦激勵電流。采用圓柱坐標(biāo)系O－(r，φ，z)，各單位矢量分別為er，eφ，ez。設(shè)線圈軸與Z軸重合，與任意導(dǎo)電平板表面的交點(diǎn)為原點(diǎn)O，線圈的半徑為r0，線圈在平板上方的提離距離為z0。任意層導(dǎo)電結(jié)構(gòu)最下方第1層，記為L1，設(shè)其相對磁導(dǎo)率為μr1，電導(dǎo)率為σ1，在實(shí)際運(yùn)算中可假設(shè)其厚度為無限大或?yàn)?。自下往上，各層結(jié)構(gòu)編號、電磁參數(shù)與厚度表達(dá)式依次增加，即第i層記為Li，相對磁導(dǎo)率為μri，電導(dǎo)率為σi。最上方第n層導(dǎo)電平板記為Ln。通電圓環(huán)線圈下方空氣層(0＜z＜z0)記為Ln+1，線圈上方(z＞z0)空氣層記為Ln+2。

圖1 線圈位于徑向被圓柱邊界截取的多層導(dǎo)電結(jié)構(gòu)上方

許多文獻(xiàn)在研究多層導(dǎo)電結(jié)構(gòu)時，研究區(qū)域?yàn)闊o窮大，電磁場以傅立葉－貝塞爾積分項(xiàng)的形式表示，在數(shù)值計算時收斂很慢。將積分變?yōu)榧墧?shù)形式，并通過調(diào)整級數(shù)的求和項(xiàng)數(shù)來控制計算的精度和速度可達(dá)到提高計算速度的目的。為此，在距離線圈適當(dāng)距離處施加一個圓柱形邊界并令邊界處的解為零，就得到一個被截取的圓柱形解域，這樣雖然會增加解的誤差，但只要包含足夠多項(xiàng)，用級數(shù)解代替積分解就不會影響計算精度［10－11］。

1.1.2 矢量磁位的邊值問題

麥克斯韋方程組是渦流檢測的理論基礎(chǔ)。為了簡化電磁場的計算及分析，引入矢量磁位A。

圖1中所有區(qū)域被與線圈同軸、半徑為b的圓柱面截取，從下至上，共有n+2個區(qū)域，有n+1個相鄰層之間界面zk，其中k=1，2，…，n，n+1。zk也表示各界面的z軸坐標(biāo)。

半徑為r0的單線圈位于內(nèi)邊界面z=z0平面內(nèi)，其電流密度可以用δ函數(shù)寫為

式(1)中I為單匝線圈激勵電流的幅值。由于介質(zhì)和場源關(guān)于對稱軸r=0呈軸對稱分布，所以圖1中場的矢量磁位僅含有周向分量Aφ。

忽略位移電流，在一個各向同性、線性均勻媒質(zhì)中，單匝線圈在(r0，z0)處產(chǎn)生的矢量磁位A(r，z)滿足下列約束方程及其邊界條件［5］:

(1)約束方程

式(2)中i=1，2，3，…，n+2，表示相應(yīng)的層數(shù);=－jωμriμ0σi;當(dāng)I≠0 時，各區(qū)域矢量磁位非零且有界。

(2)邊界條件

在所有界面zi，i=1，2，…，n，n+1，滿足:

在界面zi，i=1，2，…，n，滿足:

圓形單線圈的上下區(qū)域，即第Ln+1區(qū)域和Ln+2區(qū)域之間界面zn+1滿足:

1.1.3 矢量磁位約束方程的通解

在圓柱坐標(biāo)系下，約束方程式(2)可寫為:

用分離變量法求解式(6)得:

式(7)中J1(λjr)和Y1()分別是第一類和第二類一階貝塞爾函數(shù)。λj稱為分離常數(shù)，它與r和z無關(guān)，由于在r=b處施加了截取邊界，所以引入離散本征值λj來代替Dodd和Deeds模型中的連續(xù)值，由J1(λjb)=0可求得本征值λj。式(7)中

第L1層

第Li層，i=2，…，n+1

第Ln+2層

根據(jù)邊界條件式(3)，在界面z=zi，i=1，2，…，n，n+1，滿足

根據(jù)邊界條件式(4)，在界面z=zi，i=1，2，…，n，滿足

根據(jù)式(5)、式(10)和式(11)，并利用貝塞爾函數(shù)的正交性

線圈上下區(qū)域交界面z=zn+1=z0，滿足

這樣，共得到2n+2個方程，需要求解2n+2個未知系數(shù)。

1.1.4 矢量磁位的求解

在導(dǎo)電結(jié)構(gòu)層數(shù)較多時，各層矢量磁位的系數(shù)較復(fù)雜，將式(13)重寫為

式(16)中定義

利用式(12)、式(16)，對界面z=zi，i=1，2，…，n，因?yàn)镃12=0，可以求解得到以下關(guān)系

在界面z=zn+1，因?yàn)镃(n+2)1=0，并且第Ln+1、Ln+2層為空氣，αn+1=αn+2=λj，由式(12)可得

由式(15)、式(20)可得

式(18)、式(19)可寫成矩陣形式

令Ti為式(23)中的轉(zhuǎn)換矩陣，則各系數(shù)之間關(guān)系可寫為

令

因此

將式(21)代入式(29)得

由式(22)、式(30)可得

由式(21)、式(28)得

這樣可以得到任意區(qū)域矢量磁位的系數(shù)，進(jìn)而可以寫出任意區(qū)域的矢量磁位。其中解析計算中要利用線圈上、下方區(qū)域的矢量磁位。

線圈上方區(qū)域的矢量磁位

線圈下方區(qū)域的矢量磁位

1.2 任意層導(dǎo)電平板上方多匝線圈阻抗計算

確定了由單匝線圈在Li(i=1，2，…，n+1，n+2)區(qū)域的矢量磁位Ai(r，z)后，對于圖1b)所示具有矩形截面的多匝線圈，可以通過將單匝線圈導(dǎo)出的矢量磁位沿線圈橫截面積分獲得各區(qū)域總的矢量磁位:

在獲得矩形截面多匝線圈上方Ln+2區(qū)域總的矢量磁位和下方Ln+1區(qū)域總的矢量磁位后，令中的l2為z，令中的l1為z，并將兩式相加，可以得到線圈區(qū)域任一點(diǎn)的矢量磁位Atotal。利用公式V(r，z)=jω2πrAtotal(r，z)可以獲得單線圈上的感應(yīng)電動勢，通過疊加可以得到m匝線圈的自感電壓為

從而得出任意層導(dǎo)電結(jié)構(gòu)上方m匝線圈阻抗為:

式(37)中，I為單匝線圈內(nèi)的電流幅值。

2 單層導(dǎo)電結(jié)構(gòu)上方線圈阻抗計算

為了驗(yàn)證所導(dǎo)出任意層導(dǎo)電結(jié)構(gòu)上方m匝線圈阻抗表達(dá)式(37)的正確性，下面計算單層導(dǎo)電結(jié)構(gòu)上方線圈阻抗的特例。

2.1 單層導(dǎo)電結(jié)構(gòu)上方單匝線圈的矢量磁位

則置于單層導(dǎo)電結(jié)構(gòu)上方線圈的阻抗變化為

令式(23)中的轉(zhuǎn)換矩陣Ti中的i取1，可得單層導(dǎo)電結(jié)構(gòu)中轉(zhuǎn)換矩陣T1;令式(33)、式(34)中的n=1，可分別得到單層導(dǎo)電結(jié)構(gòu)線圈上部區(qū)域L3和下方區(qū)域L2的矢量磁位A3、A2

圖2 多線圈位于單層導(dǎo)電結(jié)構(gòu)上方

2.2 單層導(dǎo)電結(jié)構(gòu)上方多匝線圈的阻抗計算

確定了單匝線圈上下區(qū)域的矢量磁位A3、A2后，對于圖2所示具有矩形截面的圓柱線圈，可以利用式(35)通過將單匝線圈的矢量磁位沿線圈橫截面積分獲得總的矢量磁位。

計算中采用線圈電流密度

可分別得到m匝線圈上方區(qū)域矢量磁位和m匝線圈下方區(qū)域矢量磁位。令中的l2=z中的l1=z，并將兩式相加，可得m匝線圈區(qū)域的矢量磁位A(2－3)為

利用式(36)，通過疊加可以得到探頭線圈的自感電壓為

從而得出單層導(dǎo)電結(jié)構(gòu)上方圓柱形線圈的阻抗表達(dá)式為:

式(42)中

m為矩形線圈的匝數(shù)，I為單匝線圈內(nèi)的電流幅值。m匝線圈的電流密度如式(39)所示。

式(38)中的單層導(dǎo)體上方線圈阻抗Zw即為式(42)。令Zw中α1=λj、μr1=1即得到無導(dǎo)體時的線圈阻抗Zh

利用式(38)即可計算出單層導(dǎo)電結(jié)構(gòu)中渦電流引起的線圈阻抗變化為

2.3 級數(shù)展開計算及有限元仿真驗(yàn)證

為了檢驗(yàn)級數(shù)計算法的正確性，采用Mathematica軟件，根據(jù)式(45)，在不同線圈激勵頻率下，對線圈阻抗變化情況進(jìn)行了計算。計算中所采用線圈及導(dǎo)電結(jié)構(gòu)參數(shù)如表1所示。

表1 線圈及導(dǎo)電結(jié)構(gòu)相關(guān)參數(shù)

為驗(yàn)證級數(shù)計算法的正確性，同時采用有限元法對置于單層導(dǎo)電板上方線圈的阻抗變化進(jìn)行仿真計算。一般文獻(xiàn)通常采用ANSYS軟件［12］，這里采用Maxwell 2D電磁場計算軟件進(jìn)行仿真。首先構(gòu)造如圖3所示的2D模型，由于模型關(guān)于線圈軸對稱，所以只需畫出一半;然后指定模型的材料屬性，并設(shè)定邊界條件和激勵源;在模型左邊界施加奇對稱邊界條件，在線圈上施加電流激勵源;設(shè)定求解參數(shù)，將線圈上的阻抗與電感值包含在計算過程中。

如圖3所示進(jìn)行網(wǎng)格劃分，并設(shè)定求解器殘差，最后啟動求解過程。

圖3 2D有限元模型及網(wǎng)格劃分

在不同激勵頻率下，將級數(shù)展開法與有限元仿真計算所得的阻抗變化結(jié)果進(jìn)行了對比，其中電阻和電抗的對比分別如圖4、圖5所示。在級數(shù)法計算中，取解域截取半徑b=30r2，計算級數(shù)前150項(xiàng)的和。

圖4 線圈電阻變化與激勵頻率關(guān)系

圖5 線圈電抗變化與激勵頻率關(guān)系

由圖4、圖5可見，級數(shù)法計算所得線圈阻抗變化值與有限元仿真計算結(jié)果能夠很好吻合，驗(yàn)證了級數(shù)展開計算模型的正確性。

3 討論

為進(jìn)一步研究解域截取半徑與求和項(xiàng)數(shù)的不同選擇對級數(shù)法計算線圈阻抗變化值精度的影響，在激勵頻率為1 kHz、10 kHz、50 kHz和100 kHz下，分別取解域截取半徑b=10r2和b=30r2，分別計算級數(shù)前50項(xiàng)和前150項(xiàng)的和，計算中所采用線圈及導(dǎo)電結(jié)構(gòu)參數(shù)如表1所示。計算各參數(shù)下級數(shù)法與仿真計算法所得線圈阻抗變化值之間的相對誤差，結(jié)果如表2所示。

表2 級數(shù)計算與有限元仿真結(jié)果之間的相對誤差單位:%

由表2中可見，在較大的截取半徑b=30r2情況下，如果級數(shù)相加項(xiàng)數(shù)為較少的30項(xiàng)，則所得計算結(jié)果誤差大于較小截取半徑b=10r2的情況;但如果將級數(shù)求和項(xiàng)數(shù)增大到150項(xiàng)時，計算精度就大于b=10r2的情況。這說明一個較大的截取半徑需要更多的求和項(xiàng)來保證較好的計算精度。

級數(shù)計算式(42)、式(45)中只有截取半徑b和求和項(xiàng)數(shù)需要確定，積分項(xiàng)Int(x1，x2)僅依賴于線圈內(nèi)外半徑r1和r2，而與線圈的提離和材料特性無關(guān)，可以提前計算。因此級數(shù)計算法可獲得較快的計算速度。

4 結(jié)論

從Maxwell方程組出發(fā)，引入矢量磁位，采用解域截取和變量分離的渦流檢測線圈阻抗計算方法，將有限層導(dǎo)電平板推廣到任意層，導(dǎo)出了位于任意層有限區(qū)域?qū)щ娖桨迳戏骄€圈阻抗變化的模型，用易于計算的級數(shù)解代替積分解，拓寬了經(jīng)典線圈阻抗計算公式的應(yīng)用范圍，提高了計算速度，這在需要計算大量阻抗值的情況下非常有用。對單層導(dǎo)電結(jié)構(gòu)的特例進(jìn)行了計算并用有限元仿真方法進(jìn)行了驗(yàn)證，證明了級數(shù)展開計算方法的正確性，可以應(yīng)用在多層導(dǎo)電結(jié)構(gòu)缺陷渦流檢測法定量檢測正向模型中。

［1］Theodoulidis T，Poulakis N，Dragogias A.Rapid Computation of Eddy Current Signals from Narrow Cracks［J］.NDT＆E International，2010，43:13－9.

［2］Guohou Li，Pingjie Huang，Peihua Chen，et al.Quantitative Nondestructive Estimation of Deep Defects in Conductive Structures［J］.International Journal of Applied Electromagnetics and Mechanics，2010，33(3－4):1273－1278.

［3］Sun Y S，Ouyang T，Udpa S.Remote Field Eddy Current Testing:One of the Potential Solution for Detecting Deeply Embedded Discontinuities in Thick and Multilayer Metallic Structures［J］.Materials Evaluation，2001，59(5):632－637.

［4］張玉華，孫慧賢，羅飛路，等.一種用于盤孔裂紋檢測的差動式渦流探頭的設(shè)計與實(shí)現(xiàn)［J］.傳感技術(shù)學(xué)報，2008，21(6):1079－1083.

［5］Dodd C V，Deeds W E.Analytical Solutions to Eddy Current Probe-Coil Problems［J］.Journal of applied physics，1968，39(6):2829－2838.

［6］Luquire J W，Deeds W E，Dodd C V.Alternating Current Distribution Between Planar Conductors［J］.Journal of Applied Physics，1970，41(10):3983－3991.

［7］雷銀照，馬信山.渦流線圈的阻抗計算［J］.電工技術(shù)學(xué)報，1996，11(1):17－20.

［8］黃平捷，吳昭同，鄭建才，等.多層厚度電渦流檢測反演算法及實(shí)驗(yàn)研究［J］.儀器儀表學(xué)報，2005，26(4):428－432.

［9］范孟豹，黃平捷，葉波，等.多層導(dǎo)電結(jié)構(gòu)厚度電渦流檢測解析模型及實(shí)驗(yàn)驗(yàn)證［J］.中國電機(jī)工程學(xué)報，2008，28(27):142－147.

［10］Theodoulidis T，Kriezis E.Series Expansions in Eddy Current Nondestructive Evaluation Models［J］.Journal of Materials Processing Technology，2005，161(5):343－347.

［11］于亞婷，杜平安，李代生.電渦流傳感器線圈阻抗計算方法［J］.機(jī)械工程學(xué)報，2007，43(2):210－214.

［12］任吉林，刁海波，唐繼紅，等.渦流傳感器提離效應(yīng)的ANSYS模擬［J］.傳感技術(shù)學(xué)報，2008，21(6):967－971.