基于通量校正傳輸算法的局部放電光學信號機理仿真

張 悅， 錢 勇， 劉 偉， 臧奕茗， 盛戈皞， 江秀臣

(1. 上海交通大學 電氣工程系， 上海 200240； 2. 國網安徽省電力有限公司電力科學研究院， 六氟化硫氣體特性分析與凈化處理技術國網公司重點實驗室， 合肥 230022)

局部放電在線監測對于電力設備絕緣缺陷故障診斷具有重要的意義.目前常用的檢測手段包括特高頻檢測法、超聲波法及光測法等[1]，其中光測法具有絕緣性能優異、響應時間短及抗干擾能力強等優點，因而受到廣泛關注.

認知光信號產生的機理、各種光信號的特征及其變化規律是實現電力設備局部放電光學檢測及可靠診斷的基礎.絕緣缺陷導致的不均勻電場會形成電子崩、電離和復合的過程，進而產生大量光子.在光子的作用下，產生光致電離，形成較為復雜的流注放電.流注放電發展速度快、放電量較大，是局部放電檢測的主要對象.因此，在研究局部放電光信號的產生機理時，可以重點對流注放電進行仿真研究.

流注放電過程一般等效為流柱模型進行求解，將等離子體視為流體，建立電子、陽離子以及陰離子的對流擴散方程，并與電場控制的泊松方程相互耦合.在氣體放電過程中，除了電子碰撞電離之外，光電離對放電發展起主要作用.因此，作為流體模型的源項，光電離項的大小可以在一定程度上表示局部放電產生的光信號的輻射強度，分析光電離項的變化規律和特性有助于了解局部放電光信號的產生機理.

目前國內外學者對光電離的研究主要集中在快速求解上.早期光電離項通過對空間中每一個點進行全域積分得到[2].1995年，Pancheshnyi等[3]提出了 Eddington近似，進而提升了運算效率.2000年，Kulikovsky等[4]提出將光電離輻射的發射量限制在放電主軸周圍的一個小圓柱體上，并將其劃分為小環.2006年，Ségur等[5]進一步分析光電離模型，將光電離項與光子分布函數聯系在一起.2007年，Luque等[6]先后提出Helmholtz方程進行近似研究，也是目前光電離項的主要研究方法.2011年， Bonaventura 等人在Ségur基礎上研究了N2電離過程中光子通量的變化[7].國內學者蔡新景等[8]對氮氣、氧氣及其混合氣體進行流注仿真.汪沨等[9]基于歐拉-泰勒-伽遼金 (ETG)-通量校正傳輸法對短間隙SF6/N2混合氣體流注放電進行數值仿真，其重點均為分析不同氣體下的流注放電特性，鮮有文章對光子通量的變化進行研究.

本文采用通量校正傳輸(FCT)算法，通過COMSOL LiveLink with MATLAB對二維平面針板結構下的局部放電模型進行仿真，模擬尖端缺陷時的放電過程并對流注運動時光子通量的變化規律以及特性進行研究.

1 數學模型的建立

1.1 流體模型

局部放電過程中，由于電子的碰撞電離，原子或分子躍遷到高能級態，而當被激發的原子或分子回到低能級態時會向外產生電磁輻射，通常以紫外或可見光的形式散發，這部分能量又會進一步激發新的原子或分子產生光電離.文獻[10]中普遍認為空氣中的光電離是由氧分子吸收激發的氮分子發出的光子引起的，光的波長范圍為98～102.5 nm.放電過程中，電子和離子的擴散和遷移通常當作流體進行求解，根據Zheleznyak的光電離模型，其過程表示為[4]

Sph+αne|μeE|-ηne|μeE|-kepnenp

(1)

Sph+αne|μeE|-kepnenp-knpnnnp

(2)

ηnenp-knpnnnp

(3)

(4)

式中：ne、np及nn分別為未知的電子、陽離子及陰離子的數密度；μe、μp、μn、De、Dp及Dn分別為電子、陽離子及陰離子的遷移率和擴散系數；E為電場強度；Sph為光電離項；α、η分別為Townsend電離系數與吸附系數；kep、knp分別為電子與陽離子、陽離子與陰離子的復合系數；V為電勢；e為電子電量；ε0為真空介電常數.式(1)～(3)為描述等離子體的對流擴散方程，式(4)為電場的泊松方程，涉及系數的具體數值參考文獻[11].

1.2 光電離項的計算

關于空間中某點的Sph的計算，可以由整個空間中每一個點對該點積分得到.光子分布由激發態物質密度得到，光子密度的分布函數為

(5)

式中：V′為體積.光電離源項的計算式為

(6)

(7)

(8)

(9)

g(R)=

(10)

考慮到分布函數的計算復雜性，Sph通常采用3組Helmholtz方程近似得到[8]：

(11)

(12)

式中：ps為物質s的分壓；系數Aj、λj分別采用表1中的數據[8]，1Torr=133.3 Pa；I(r) 為單位體積內光子的產生數量.忽略光子在邊界上的光電發射效應，在邊界上設置Sph= 0，則可求得Sph(r)，即單位時間r處由于光電離作用而產生的電子數目.

表1 Helmholtz方程系數

1.3 光子通量的計算

激發態的原子分子是由電離產生的高速電子碰撞產生，而單位體積內光子的產生數量通常也正比于碰撞電離產生的電子數量：

(13)

式中：ξ為光電離效率，取值范圍為0.01～0.1[10].

2 數學模型的求解

本文通過超松弛迭代法(SOR)算法求解1.1節中的泊松方程.采用FCT算法對對流擴散方程即式(1)～(3)進行求解.這種高精度算法的基本思想是加入非線性的平衡反擴散來抑制數值發散，同時有效地避免計算過程中產生非物理負值[12].FCT算法分為高階格式與低階格式兩部分，高階格式為數值算法提供足夠高的數值計算精度，低階算法保證數值計算的穩定性和非負性.本文中高階格式采用BDF4 (Fourth-order Euler backward difference)[13]，這種隱式格式具有很好的收斂性和穩定性，可以采用較大的時間步長進行計算，從而減小總體運算量，低階格式通過在高階格式基礎上添加調制耗散得到：

MLun+1=(ML+D)uH

(14)

為了驗證FCT算法的有效性，本文采用一維及二維對流擴散方程進行方波檢驗，結果如圖1～2所示.圖中：A為幅值，t為時間.可以看出，算法得到的方波在運行過程中幅值基本保持不變，波形也沒有發生較大的改變，精度滿足要求[14].

圖1 一維平流方程FCT算法檢驗

圖2 二維平流方程FCT算法檢驗

3 仿真模型的建立

明確數學模型及其求解方法后，本文采用COMSOL建立仿真模型，對二維平面針板結構的流注模型進行仿真，模擬空氣中尖端缺陷局部放電發生時的放電過程.

針尖采用雙曲型結構，設置陽極電壓為15 kV，陰極接地，針板之間距離為6 mm，初始等離子云分布在針尖前0.15 mm處，滿足高斯分布，其輪廓如圖3所示.針尖以及初始等離子體的具體表示方程如下：

(15)

(16)

式中：x、y為針尖的橫、縱坐標，用于確定針尖的幾何輪廓，針尖為雙曲線結構.圓錐曲線參數a=1 mm，b=0.07 mm；參數n0=1020m-3；概率分布函數的參數y0=0.15 mm，σx=0.05 mm，σy=0.05 mm.幾何模型和網格剖分如圖4所示.為了提升計算精度， 在流注通過的路徑上將網格剖分尺度設置為最小尺寸3.0 μm.針尖處的曲率半徑為0.02，網格增長率為1.01，共剖分 41 398 個單元，時間步進為0.005 ns，對 3 249 個步長進行求解，共耗時5 h 32 min.

圖3 初始等離子體的輪廓

圖4 針板的幾何模型與網格剖分

4 結果與分析

通過FCT算法對上一章中建立的放電仿真模型進行求解，即可得到等離子體對流擴散方程中各粒子的分布.本節分析在流注發展過程中光電離項與光子通量的大小，并得到局部放電發生時光信號輻射強度的變化規律.將流注開始后的0～5 ns記為流注發展前期，5～7 ns記為流注發展后期.

4.1 流注發展初期

在流注發展的初始階段，空間電荷效應還不明顯，此時場強和光子通量(Φ)的最大值都集中在陽極頭部附近，如圖5所示.由于電離效應大于吸附效應，陽極附近會產生大量的電子和離子，導致電場強度和光子通量不斷增強.

圖5 初始電場和光子通量分布

4.2 流注發展前期

在光電離效應的作用下，流注頭部前方產生大量的光電子，新產生的電子在向流注頭部移動過程中又持續電離出新的離子，這部分的電子和離子使得流注通道向前延長.如圖6(a)所示，流注頭部附近的場強最大且流注體內的場強遠低于周圍電場強度.這是由于流注體外部的等離子體鞘層中陽離子數目大于電子和陰離子數目，使得空間電荷產生的電場方向與電極的電場方向相反，導致流注體內的電場被屏蔽.圖6(b)為軸線處電場強度在0～5 ns內的變化，從左至右依次表示不同時間的軸線電場強度分布，圖中L為弧長.可以看到在流注發展前期，流注發展速度較慢，隨著流注向前移動，頭部場強先是逐漸減小隨后穩定在8×106V/m，這是由于背景的拉普拉斯電場隨著遠離電極不斷減小導致的.

如圖7所示，由于光子數與電離項近似呈正比，光子通量最大值主要集中在流注頭部附近且變化趨勢基本與場強保持一致.隨著流注的向前發展，發光半徑不斷變大.這里由于忽略了激發態原子分子的壽命，流注通道內幾乎不產生光子[10]. 圖7(b)中從左至右依次表示不同時間的軸線光子通量分布，最右處藍色線表示t=0時刻軸線光子通量分布，最左側為t=5 ns時刻軸線光子通量分布，每條線的時間間隔為0.2 ns.

圖6 0～5 ns電場強度變化

圖7 0～5 ns光子通量變化(從左到右t分別為5, 4.8,…,0 ns)

4.3 流注發展后期

(1) 場強變化.隨著流注不斷靠近陰極，由于流注前方到陰極之間的場強不斷增強，使得流注的運動速度變快，約為1.8×106m/s.從圖8中可以看到，流注頭部的半寬不斷變長，相對流注發展初期增大2.5～3倍，約為0.5 mm.圖9為軸線處電場強度的變化，從左至右依次表示不同時間的軸線電場強度分布，最右藍色線表示t=5 ns時刻軸線電場強度分布，最左側為t=7 ns時刻軸線處電場強度分布，每條線的時間間隔為0.2 ns.可以看出，在靠近陰極時流注前方的場強迅速增強.

圖8 5～7 ns電場強度分布

(2) 光子通量變化.如圖10所示，流注發展后期，發光點的發光半徑明顯增加，由0.1 mm增長到0.5 mm，且隨著流注前方的電場強度的增強，光子通量幅值呈指數的增長，如圖11所示.圖中從左至右依次表示不同時間的軸線光子通量分布，最右處藍色線表示t=5 ns時刻軸線光子通量分布，最左側為t=7 ns時刻軸線光子通量分布，每條線的實際間隔為0.05 ns.可以看出，當流注到達陰極后，光子通量約為1.4×1028s-1·m-3.

圖9 5～7 ns電場強度變化(從左到右分別為t=7，6.8，…，5 ns)

圖10 5～7 ns光子通量分布

圖11 5～7 ns光子通量變化(從左到右分別為t=7，6.95，…，5 ns)

5 結語

通過數值研究針板結構下的正流注模型，模擬尖端缺陷局部放電過程.計算模型以3種基本粒子的漂移擴散方程和泊松方程為基礎，得到電子密度、電場以及光子通量隨時間的變化規律，解釋了局部放電光信號的產生機理及發展過程.結果表明：空氣放電發生時流注的移動速度約為8.6×105m/s，流注通道的最大半寬為0.5 mm，相對初始寬度增大2.5～3倍.正流注的發光位置往往集中在快速移動的流注頭部，隨著流注等離子體的發展向陰極運動，光子通量的最大值約為1.4×1028s-1·m-3，隨著流注向陰極移動緩慢上升，發光的半寬由0.1 mm增長到0.5 mm，當電壓足夠高時，產生的光子通量與流注頭部到陽極的距離近似呈指數關系.由于空間電荷層的寬度(或等效的電離區寬度)是由光電離輻射的特征吸收長度決定的.所以，可以通過改變氣體的光電離長度，控制等離子體的形狀、密度，并在一定程度上控制流注的傳播速度，進而可以控制局部放電的放電特性.