基于COMSOL的流注頭部分叉過程仿真與分析

朱毅佳 朱 武* 張佳民

1(上海電力大學電子與信息工程學院 上海 200090) 2(上海電力大學自動化工程學院 上海 200090)

0 引 言

隨著我國特高壓輸電工程的深入開展，電壓等級不斷升高，高壓設備空氣間隙放電極易發生，對電網的安全、穩定運行構成潛在的危害。如何優化輸變電工程的空氣間隙結構一直是高電壓工程亟待解決的核心問題之一[1]。因此，對空氣間隙放電特性的研究具有重要的理論價值和實際應用價值。20世紀以來，國內外學者對氣體放電的研究不斷深入，一般認為，流注是一種以弱電離、低溫等離子體絲的形式傳播的放電現象，其頭部由強電場驅動。在許多的氣體放電實驗中，放電通道在發展過程中分裂成若干條樹枝一樣外形的放電通道，這種現象被稱為流注分叉[2]。流注在唯象層面的分叉圖樣可以通過分形生長模型(DBM)生成，然而這種唯象層面的生成方法并沒有反映出流注分叉出現的物理原因。關于流注分叉現象的原因，文獻[3-5]認為流注頭部分叉是流注頭部空間凈電荷薄層在發展過程中出現Laplace不穩定性的表現。凈電荷薄層與流注前方電場形成正反饋，其越接近平坦，不穩定性就越高。文獻[6-7]給出了由Laplace不穩定性導致流注分叉的物理圖像。Li等[8]采用粒子模型研究大氣中光電離對流注形成和傳播的影響，獲得了流注的分叉結構。

縱觀流注放電的研究現狀，描述放電過程的流體模型仍是流注數值仿真的常用手段。張赟等[9]采用伽遼金有限元法求解氣體放電流體模型，結合FCT方法抑制計算誤差，提高了計算精度和效率。Li等[10]與Teunissen博士和Ebert教授合作，開發出一套流注放電的三維粒子仿真模型，反映真實的放電結構特征。

本文通過對流擴散方程建立流體力學模型，將初始電離區域形狀簡化為兩個具備高斯分布的初始電離區域的疊加，使用COMSOL Multiphysics軟件在二維笛卡爾坐標系中對標準大氣壓下板-板電極結構中流注放電的過程進行模擬，目的是研究初始電離區域形狀和電壓對流注分叉情況的影響，觀察流注通道的形成發展過程并對流注放電的頭部分叉情況做出預測。

1 仿真模型

1.1 物理模型

搭建板-板電極間均勻電場的模型框架如圖1所示。其中：陰極和陽極為兩平行矩形板，間距d為10 mm，背景氣體被簡化為N2-O2體積比為8 ∶2的氮氧混合氣體，氣壓p為101 kPa，環境溫度設為300 K。

圖1 平行板模型x-y截面圖

在二維笛卡爾坐標系中，陽極(電壓為+U)位于y=10 mm處，陰極(接地)在y=0處。

陰極表面的初始電離種子沿徑向和軸向均滿足高斯分布[11]如下：

(1)

為了更好地描述初始電離區域形狀，將兩個滿足高斯分布的電離種子疊加構成一個整體。

描述該整體的高斯公式為：

(2)

式中：ne和np分別是放電過程中的電子數密度和正離子數密度；n0=2×1018m-3表示初始種子團的密度，即初始預電離水平；±x0是兩個高斯中心的橫軸坐標，令D=2x0；σy和σx是初始電子分別在y軸和x軸方向的特征尺寸，表示其分布的空間寬度。

1.2 化學反應模型

為進一步探討流注分叉情況和相關電子密度的關系，本文以電子和正負離子在電場中的漂移擴散反應為流注模型的控制方程,以標準大氣壓下的流注為研究目標，采用仿真軟件COMSOL Multiphysics的等離子體模塊來模擬流注的形態。

雖然文中的空氣已做簡化處理，但氮氧混合氣體中生成的化學反應種類繁多且復雜，分析全部的反應意義不大，本文選取了具有代表性的反應模型，展示了包括彈性碰撞、電離碰撞、相應附著反應和粒子間反應在內的21種化學反應，如表1所示。

表1 計算中考慮的化學反應式

續表1

彈性碰撞和電離反應速率k1-k4及有效電子溫度Te(單位為K)由BOLSIG+軟件求解玻爾茲曼動力學方程計算得到。反應方程5-21的反應系數為常數來自文獻[13]。

1.3 方程的邊界條件及光電離

本文采用傳統流體動力學模型，其由粒子的連續性方程(又稱對流擴散方程)耦合電勢的泊松方程組成，光電離則作為源項加在電子和正離子連續性方程上[14]。

由Hagelaar等[15]提出的電極表面處的電子通量邊界條件為：

(3)

(4)

在實際應用中可簡化為：

(5)

式中：E為場強；μe為電子磁矩；Γe為電子通量；n為電極表面的法向量；re為經典二次電子反射系數；vth是電子的熱力學速度；kB為玻爾茲曼常數；αe是方向函數(當電子通量指向電極時為1，相反時則為0)。

電極處的離子通量可設置為[16]：

(6)

式中：Γk為離子通量；nk和mk為第k類離子的離子密度和離子質量；Tk為第k類離子的溫度；γk為離子在電極表面衰變成中性的反應系數，表示反應的概率。

在非電極邊界上則應用零梯度密度邊界條件：

n·(?ni)=0

(7)

本文中光電離項的計算是選取了文獻[17]中的三項Helmholtz模型以及文獻[18]中使用過的模型。三項Helmholtz模型計算結果準確，在遠離表面的區域，它的計算效率更高。在不考慮邊界表面光電離作用的情況下，利用Dirichle邊界條件為0來降低計算成本[19]。

在該仿真計算中最細網格(靠近對稱軸x=0)設置為0.1 μm，元素總數大約為107。采用MUMPS直接線性求解器進行求解，盡管該求解器的內存開銷大，但其魯棒性強。本文所有仿真結果均使用Intel(R) Xeon(R) E5-2637 v3處理器，且需要10～20個小時來模擬流注的發展過程。

2 仿真結果與討論

2.1 初始電離區域形狀對流注發展的影響

為探討距離D其臨界值的意義和相關電子密度的影響，并進一步探究種子解的包絡面和流注頭部分叉之間的聯系，本文根據式(2)繪制了特征尺寸σx=σy=0.45 mm時不同距離D下的初始種子包絡面圖，如圖2所示。由圖2(d)可明顯看出：隨著電子密度ne的增加(輪廓線由外及內表示電子密度依次增大)，包絡面的輪廓逐漸由一個較大的類半圓縮小為兩個較規整的小半圓。圖2(a)-圖2(d)中最內側的線直觀地顯示了電子密度ne=1.5×1019m-3在距離D取不同值時的初始種子包絡面示圖。

圖2 初始種子包絡面圖

經觀察發現，對于某一固定的電子密度值ne，只要適當提高距離D的值，包絡面的輪廓線就會重復相似的轉變過程(以圖2所示的情景為例)：當電子密度ne=1.5×1019m-3時距離D從0開始緩慢增加時，包絡面的輪廓中心開始向下凹陷，逐漸形成類駝峰狀；當距離D大致介于0.88 mm到0.90 mm之間時，包絡面的輪廓完全分裂成兩個對稱的半圓片，并且在剛好完成完全分裂(包絡面的輪廓線經過原點)時，D達到其臨界值。

這種通過識別與特定電子密度相關的指示輪廓線，利用輪廓線的幾何形狀以判斷流注發展過程中頭部是否分叉的方法即為包絡面輪廓指示法。在模型基礎參數固定不變的基礎上，為獲得距離D對流注發展動態的具體影響，本文模擬計算了一組不同距離的情況下電場背景均為38 kV·cm-1的電子密度分布圖，如圖3所示。

圖3 不同距離下流注放電電子密度分布圖

從圖3(a)-圖3(d)中可以明顯看出，距離D在小于0.82 mm時，只有一條頭部未發生分叉的流注由疊加的高斯種子產生，這表明距離D應存在一個臨界值，低于該臨界值時流注頭部便不發生分叉。其中圖3(a)中距離D=0的情況可視為實際只有一個初始高斯種子在發展，作為參考項。當距離D較小時(圖3(b))，流注的頭部未發生分叉的情況，而是在x軸方向上產生了橫向擴展；當距離D到達其臨界值(圖3(d))時，流注頭部開始出現分叉；當距離D增加到1 mm(圖3(f))時，流注頭部的分叉情況得到更加充分的展示。進一步增加距離D到1.4 mm時，這兩個由于距離較遠而完整分開的初始高斯種子將一條流注的頭部分叉問題轉化為兩條流注相互作用影響的問題。由文獻[20]中的結論可知，由于靜電排斥的作用，兩條流注在前進過程中會形成一定的角度，在以后的發展過程中不會產生融合的情況，仿真結果與該結論是一致的。而后當距離D再繼續增加至更大值時，靜電排斥效應將不再明顯，兩條流注將保持相互平行的狀態向陽極發展，此時兩條流注之間的相互影響可以忽略不計。

流注發展速率v也是描述流注時空特性的重要參數之一。根據仿真所得數據，可以確定流注頭部的坐標位置，進而計算出流注的發展速率。

在38 kV·cm-1的電場背景下，距離D對流注的發展速率vave(前9 ns的平均速率)的影響如圖4所示。通過分析可以發現，當距離D逐漸增大并超過0.4 mm時，其對vave的影響由微弱轉變為迅速增大，并且vave的下降加速度在距離D介于0.8～0.9 mm間開始減緩，而流注頭部分叉的臨界值剛好處于這個范圍內。但由于本文仿真數據的局限性，無法準確定位這個點。為此，本文探究出另一種指示方法。

圖4 流注平均速率示意圖

2.2 電壓對流注發展的影響

流注通道的發展特性與電場背景緊密相關，為獲得電場對流注形成和發展的影響，本文基于上述模型模擬計算了3組不同參數的情況：1) 場強為36 kV·cm-1，時間設置為10 ns；2) 場強為38 kV·cm-1，時間設置為9 ns；3) 場強為40 kV·cm-1，時間設置為7 ns。仿真結果如圖5所示。從仿真結果可以看出，在相同電壓下增大距離D的值，流注在發展過程中仍存在著3種狀態：1) 距離D低于臨界值時，只有一條頭部未發生分叉的流注由疊加的高斯種子產生；2) 距離D等于臨界值時，流注的頭部開始出現分叉；3) 距離D超過臨界值時，流注頭部的分叉更加明顯。

圖5 不同電壓下流注放電電子密度分布圖

但在更高的電壓背景下，流注的發展速率明顯提高，流注頭部的分叉現象出現得更晚，即距離D的頭部分叉臨界值與電壓背景成正相關，隨著電壓U的增高而增大(U=36 kV時為0.76 mm，U=38 kV時為0.81 mm，U=40 kV時為0.83 mm)。這表明更高的電壓或更大的場強對流注頭部的分叉起抑制作用，更加促進了頭部的融合(兩條流注由兩個疊加的高斯種子發展而來)。較為合理的解釋是，在更高的電壓背景下，光電離則更為突出，它抑制了初始種子頭部包絡面的不規則性。同時，更高的電壓加速了流注的傳播，這也不利于流注頭部分叉的出現。

2.3 包絡面輪廓指示法

基于上述對流注發展規律的總結，本文提出初始電離區域的包絡面輪廓指示法，通過與特定電子密度相關的指示輪廓線的識別，判斷流注發展頭部分叉過程：當輪廓線的幾何形狀恰好分開為兩個小半圓(指示輪廓線經過原點)時，對應的距離D即為流注頭部分叉的距離臨界值。

在電壓U分別為36 kV、38 kV和40 kV的背景條件下，剛好落在橫軸上完成分裂的曲線其電子密度分別為1.96×1019m-3、1.78×1019m-3和1.71×1019m-3。

根據包絡面輪廓指示法，利用式(2)，可以解出對應的距離D，從而確定流注頭部分叉臨界距離：U=36 kV時為0.76 mm，U=38 kV時為0.81 mm，U=40 kV時為0.83 mm。計算結果與圖5展示的仿真結果一致。

這種方法將初始電離區域形狀的幾何不規則性與流注放電的頭部分叉現象聯系起來，可以作為一種預測流注頭部分叉的簡單工具。

3 結 語

本文使用仿真軟件COMSOL Multiphysics對空氣中的負流注進行了二維模擬，研究了電壓背景及初始電離區域的形狀對流注發展的影響。主要結論如下：

1) 相同電壓背景下改變初始電離區域形狀會對流注的時空發展特性產生影響，隨著距離的增大，流注的發展速率會減緩。

2) 隨著電壓降低，流注頭部分叉的趨勢更加明顯，分叉的臨界值更小。

3) 初步探究了初始種子電離區域形狀的不規則性和流注放電頭部分叉之間的關系：當高斯中心的距離較近、高斯種子重疊較多時，通過識別與特定電子密度相關的指示性輪廓線，利用輪廓線的幾何形狀可以可靠預測流注頭部的分叉情況。