朱益飛, 陳賢聰

(1. 空軍工程大學航空工程學院等離子體動力學重點實驗室, 陜西西安 710038; 2.西安交通大學機械工程學院, 陜西西安 710049)

引 言

大氣壓低溫等離子體由于能夠可控、 節(jié)能、 高效地產(chǎn)生活性粒子, 近年來越來越受到工業(yè)和科學界重視. 從臭氧產(chǎn)生、 聚合物處理、 激光制造、 污染控制、 生物除菌、 醫(yī)藥治療, 到流動控制、 點火助燃、 薄膜工藝, 低溫等離子體的應用領域正在迅速擴展[1-8].

在大氣壓下, 產(chǎn)生等離子體的基本方法之一是在由金屬制造的針板結(jié)構(gòu)的氣隙間施加高電壓. 高電壓在針尖形成局部強電場, 驅(qū)動電子加速并引發(fā)雪崩-流注放電, 形成一條細長的放電通道, 并在數(shù)十個納秒內(nèi)擊穿整個氣隙[9-10].

數(shù)值計算是測試理論假設、 驗證實驗結(jié)論和預測放電參數(shù)與行為的有力工具. 雪崩-流注放電是一個多尺度系統(tǒng), 包含了諸多過程, 包括化學反應、 組分輸運、 電磁場傳播以及流體耦合等, 高精度、 高效率地模擬這個系統(tǒng)是一個巨大的挑戰(zhàn).

首先, 流注的動力學過程是高度非線性的. 由于在流注頭部以及與介質(zhì)間隙, 電場與空間電荷分離層直接耦合, 對應區(qū)域的計算需要極高的空間精度. 另外, 電荷分離層通常具有弧度, 因而純一維模型不足以準確描述流注的傳播. 文獻[11]使用了1.5維模型, 通過將二維Poisson方程和一維輸運方程結(jié)合實現(xiàn)了徑向方向流注的模擬, 但是該方法無法模擬負流注和表面流注. 文獻[12-13]引入了SG方法模擬二維軸對稱幾何下的火花和輝光放電. 文獻[14-15]分析了正流注和負流注的動力學過程. SG方法經(jīng)文獻[16-17]的改良, 用于精確模擬針-板放電. 文獻[18]進一步對算法進行改良, 用于計算流注與介質(zhì)接觸問題. 其次, 流注傳播是一個多空間尺度問題. 大氣壓下, 流注頭部的電荷分離層厚度僅為幾十個電子自由程大小(幾十微米). 與此同時, 流注的傳播距離卻可達數(shù)厘米. 因此, 流注模型須解析跨越了4個量級的兩種空間尺度. 最常用的解決方案是使用非均勻網(wǎng)格, 并在等離子體區(qū)域加密. 文獻[19-20]等采用了最新的結(jié)構(gòu)化自適應網(wǎng)格, 使得計算效率得到了顯著提高. 流注也是一個多時間尺度問題. 流注的傳播時間一般為幾十到數(shù)百納秒, 而大氣壓下捕捉流注輸運過程所需的最小時間尺度為10-14s, 跨越了7～8個數(shù)量級. 文獻[21-23]使用了半隱式推進, 文獻[24]使用異步時間推進的方法克服了這些困難. 再次, 流注具有強非局域特性, 使得并行加速變得十分困難. 這種非局域性主要是由于電場和光電離的全局特性. 在有介質(zhì)的情況下, 邊界特性的變化會導致局部電場的陡增和間斷. 為了克服非局域性問題, 一種解決方案是使用局域能量近似, 將化學反應速率和輸運特性表達為局域平均電子能量的函數(shù). 另一種方法是引入對電離反應速率系數(shù)的修正項. 更加直接的方法是采用混合模擬, 考慮快速電子(逃逸電子)的形成. 文獻[25]將粒子模型用于模擬流注頭部, 將流體模型用于模擬流注主體. 這種“電子束-主體區(qū)”方法綜合了電子Monte Carlo模型和流體模型, 得到了大量研究[26-28]. 在一些專注于逃逸電子、 流注分叉現(xiàn)象的研究工作中, 純粒子-Monte Carlo模擬也得到了應用[29-32].

近年來, 很多課題組采用商用有限元軟件COMSOL Multiphysics開展流注模擬. 其在模擬等離子體放電中有一些獨特的優(yōu)勢: (1)具有友善的用戶界面, 便于快速入門并獲得一些初始結(jié)果; (2)有限元方法使其可以適應幾乎任何復雜幾何結(jié)構(gòu); (3)自帶的等離子體模塊能夠大大減少編程工作量; (4)自適應網(wǎng)格特性使其能夠在一定程度上克服有限元商用軟件的通病——網(wǎng)格數(shù)量太小; (5)自帶的后處理系統(tǒng)為計算結(jié)果分析提供了有力的工具; (6)擁有通用性的偏微分方程模塊, 可以較為自由地組織模型.

盡管COMSOL Multiphysics具有以上優(yōu)勢, 但是也存在很多非常顯著且難以解決的缺陷, 導致其無法作為研究在大氣壓下通常呈流注形態(tài)的納秒脈沖放電的理想工具: (1)COMSOL Multiphysics的等離子體模塊采用了標準Galerkin方法, 該方法無法捕捉輸運占優(yōu)的流注頭部區(qū)域的強間斷面, 導致計算結(jié)果嚴重失真, 呈偽擴散形體; (2)商用軟件的“黑盒子”模式使得無法對程序運行模式進行任何修改, 因而在模擬遇到問題時只能憑借猜測處理, 幾乎沒有辦法調(diào)試; (3)全隱式的數(shù)值格式使得時間積分對計算域任何一點擾動都十分敏感, 無關緊要處的一個奇異點就會導致整個模型計算效率下降, 收斂極其困難; (4)與基于質(zhì)量守恒的有限體積法不同, 有限元方法本質(zhì)是基于能量守恒原理的, 不具有保證質(zhì)量守恒的特性, 在等離子體計算中不可避免會出現(xiàn)負濃度情況; (5)COMSOL Multiphysics直接求解大規(guī)模矩陣而不充分利用并行計算機結(jié)構(gòu), 使其存在計算速度緩慢、 內(nèi)存消耗巨大的問題.

本文介紹一種以弱形式偏微分方程形式在COMSOL Multiphysics軟件中引入人工穩(wěn)定項, 用以克服等離子體模擬中出現(xiàn)的數(shù)值振蕩和偽擴散問題的解決方案. 文章分3個部分, 第1部分以弱形式的方式推導人工穩(wěn)定項的輸入形式, 第2部分給出基準驗證和實驗驗證, 最后給出結(jié)論.

1 弱形式與人工穩(wěn)定項

分析引言所述的COMSOL Multiphysics缺陷, 發(fā)現(xiàn)缺陷(3)和(4)是當?shù)入x子體與介質(zhì)、 金屬接觸時存在的最主要問題, 而對于針-板型的體放電而言, 這些缺陷不會對計算結(jié)果產(chǎn)生太大影響. 缺陷(5)無法克服, 而缺陷(1)和(2)可以通過引入弱形式偏微分方程, 添加人工穩(wěn)定項來解決. 本節(jié)將對在COMSOL中構(gòu)建弱形式PDE引入流線擴散人工穩(wěn)定性(streamline upwind Petrov Galerkin method, SUPG)進行闡述.

首先, 利用COMSOL Multiphysics建立一套最簡單的等離子體方程體系: Poisson方程、 Helmholtz方程和輸運方程. 其中, Poisson方程和Helmholtz方程是單純的橢圓形方程, 可以直接使用偏微分方程模塊以“黑盒子”的形式求解. 輸運方程則用弱形式偏微分方程表達.

在納秒放電時間尺度下, 可以假設離子和中性粒子不運動, 僅求解不包含輸運項的反應方程即可

式中,ni為離子密度,Di為離子擴散率,Si為化學反應主導的離子產(chǎn)生率,Sph為光電離主導的離子產(chǎn)生率.

電子的輸運方程為

式中,ne為電子密度,De為電子擴散率,μe為電子遷移率,Φ為電勢,Se為化學反應主導的電子產(chǎn)生率,Sph為光電離主導的電子產(chǎn)生率. 電子擴散率、 遷移率和化學反應電子產(chǎn)生率均可通過開源電子能量分布函數(shù)計算軟件BOLSIG+計算得到. 為了便于比對, 本文中的反應速率系數(shù)和輸運系數(shù)均采用了文獻[16]中給定的公式計算

以上方程乘以“試函數(shù)”u(x,y,z), 并對有限單元進行體積分, 可以得到

對上式前兩項進行分部積分, 并簡化可得

由此, 可以得到能夠直接輸入COMSOL Multiphysics的弱形式方程

(1)

式中, 第1行是等離子體計算域的弱形式偏微分方程, 第2行代表自然滿足的邊界條件. 式(1)可以寫成更明晰的形式

方程(1)所描述的漂移-擴散過程, 如果求解的區(qū)域被擴散過程主導, 那么傳統(tǒng)的COMSOL標準Galerkin差分方法可以獲得穩(wěn)定可靠的解. 如果計算區(qū)域是漂移主導, 那么方程(1)則表現(xiàn)出雙曲型方程特性, 傳統(tǒng)的數(shù)值方法變得非常不穩(wěn)定, 計算結(jié)果要么會在粗糙的網(wǎng)格下劇烈振蕩, 要么在精細的網(wǎng)格下形成偽擴散解. 在前述方程的基礎上, 利用COMSOL直接復現(xiàn)文獻[16]中的軸對稱針-板放電基準算例, 采用相同的方程(連續(xù)方程+Poisson方程), 分別使用50 μm和15 μm兩種尺寸網(wǎng)格開展計算, 結(jié)果如圖1(a)～(d)所示.

圖1(a)和(c)對比了使用COMSOL Multiphysics在兩種網(wǎng)格條件下計算的流注軸向電子密度分布(藍線)和基準算例中的軸向流注通道電子密度(紅線). 在粗網(wǎng)格下, COMSOL Multiphysics計算所得流注能夠復現(xiàn)基準算例流注傳播速度, 如圖1(b)所示, 但是整個放電通道內(nèi)出現(xiàn)了圖1(a)中劇烈的非物理數(shù)值振蕩, 計算結(jié)果無法分析理解. 如果僅僅加密網(wǎng)格, 如圖1(c)所示, 計算軸向電子密度不再產(chǎn)生數(shù)值振蕩, 但是流注傳播距離遠遠落后于基準算例, 圖1(d)表明, COMSOL Multiphysics標準有限元方法無法捕捉到基準算例中的電子雪崩-流注轉(zhuǎn)捩, 計算結(jié)果呈擴散流注狀. 這兩個案例清楚地表明, 使用COMSOL Multiphysics默認的數(shù)值格式和計算方法不可能正確計算大氣壓流注. 本文使用COMSOL Multiphysics內(nèi)置的等離子體模塊進行測試, 也得到了同樣的錯誤結(jié)果.

(a) Electron density in coarse meshes

(b) E/N in coarse meshes

(c) Electron density in fine meshes

(d) E/N in fine meshes

這種數(shù)值不穩(wěn)定性, 源于計算網(wǎng)格與算法不匹配, 有限元中的標準Galerkin方法與有限體積/有限差分中的標準SG方法一樣, 在高氣壓條件下, 計算等離子體流注的準確性會顯著下降. 假設電子擴散系數(shù)符合Einstein關系, 那么為了保證計算結(jié)果的準確性, 須保證以下關系[17]

(2)

式中,k表示第k個網(wǎng)格節(jié)點, ΔE為電場強度梯度,h為網(wǎng)格尺寸,Te為電子溫度. 上式表明, 兩個相鄰節(jié)點的電勢差須遠小于電子溫度. 對于高氣壓情況, 若要滿足該條件, 則需要極高的網(wǎng)格密度, 這使得高性能計算幾乎不可能. 方程(2)所規(guī)定的嚴苛的網(wǎng)格條件在COMSOL下基本不可能實現(xiàn)(除非計算域非常小). 為此, 可以通過引入人工穩(wěn)定性[33-37]提升計算精度. 最廣泛使用的人工穩(wěn)定性方法之一, 就是流線擴散穩(wěn)定SUPG方法[33-34]. 通過弱形式方程, 在COMSOL Multiphysics中對各個有限單元添加人工穩(wěn)定項

(3)

式中, 各項的含義如下

將方程(3)以弱形式方程的方式代入COMSOL中電子密度輸運方程弱形式, 即可得到穩(wěn)定的解. 針對該改進的COMSOL模型, 具體的驗證工作將在下節(jié)詳細介紹.

2 基準驗證與實驗驗證

針對改進的COMSOL Multiphysics流注放電模型, 基于經(jīng)典計算結(jié)果[16,38]和實驗結(jié)果[39]開展驗證.

2.1 計算基準驗證

針板放電是一種經(jīng)典和基礎的等離子體放電結(jié)構(gòu), 國外前期開展了大量研究, 主要采用有限體積方法編程求解連續(xù)方程、 Poisson方程與光電離. 其中, 俄羅斯科學院Kulikovsky等對針板放電開展了數(shù)值計算研究, 并構(gòu)建了一系列便于驗證的數(shù)值模擬基準算例. 文獻[16]計算了1 cm間隙條件下, 大氣壓針板放電構(gòu)型在13 kV恒定電壓下的擊穿過程, 考慮了電子碰撞電離反應、 三體和雙體吸附反應、 電子-離子復合和離子-離子復合反應, 成功復現(xiàn)了雪崩-流注轉(zhuǎn)捩現(xiàn)象和流注傳播特性. 文獻[16]的參數(shù)、 方程完整詳實, 得到了眾多課題組的復現(xiàn), 被國際上廣泛認可為大氣壓低溫等離子體數(shù)值計算的重要基準算例之一.

本文第1部分對比了使用標準有限元方法計算的結(jié)果與文獻[16]的結(jié)果, 在圖1中展示了數(shù)值振蕩和偽擴散問題. 本節(jié)使用修正后的COMSOL Multiphysics模型, 采用與文獻[16]相同的控制方程、 輸運參數(shù)和反應體系再次開展計算, 結(jié)果如圖2所示.

圖2(a)和(c)分別為文獻[16]計算得到的不同時刻軸向電子密度分布和電子密度空間分布云圖. 圖2(b)和(d)為采用修正COMSOL Multiphysics模型計算結(jié)果. 由圖可見, 添加人工穩(wěn)定性的修正模型計算結(jié)果與基準算例幾乎沒有區(qū)別: 電場強度、 電子密度與傳播過程均實現(xiàn)了準確的復現(xiàn). 注意到使用修正COMSOL Multiphysics模型計算的流注直徑比基準算例稍細, 引起這一區(qū)別的原因可能是光電離計算方式的不同: 文獻[16]采用了簡化的“積分環(huán)”方式, 僅計算流注頭部強電場附近的光電離源項, 并將全局積分簡化為對流注頭部附近的環(huán)形體積進行積分. 本文采用了近年來使用較為廣泛的Helmholtz三方程光電離模型[40], 該模型的準確性已經(jīng)得到了廣泛驗證[41], 本文不再贅述. 圖1中基準算例的紅線與圖2(a)(b)中左數(shù)第3條線(19 ns時刻)對應. 對比圖1， 2可見, 修正后的模型計算結(jié)果中已經(jīng)沒有了圖1藍線出現(xiàn)的數(shù)值振蕩和偽擴散的問題.

(a) Axial electron desnity from benchmark

(b) Axial electron desnity from corrected COMSOL model

(c) Plot of electron density contours from benchmark

(d) Electron density contours from corrected COMSOL model圖2 修正的COMSOL流注模型復現(xiàn)Kulikovsky基準流注算例[16]Fig. 2 Results comparison between corrected COMSOL model and Ref[16]

需要注意的是, 本節(jié)進行的計算基準驗證的核心目的是數(shù)值計算角度確保: 在相同的輸運參數(shù)、 幾何結(jié)構(gòu)和邊界條件下, 求解相同的控制方程能夠得到完全相同的解, 從而確保數(shù)值計算研究的一致性. 文獻[16]的結(jié)果作為數(shù)值計算的良好基準, 卻并未得到實驗的充分驗證. 在克服標準有限元方法引發(fā)的數(shù)值振蕩、 偽擴散問題的基礎上, 使用修正的模型與實驗結(jié)果進行了進一步比對驗證.

2.2 實驗結(jié)果驗證

利用修正后的COMSOL Multiphysics等離子體模型, 針對一例過電壓實驗進行直接比對[39]. 文獻[39]利用相增強高速相機ICCD研究了大氣壓下, 1.6 cm間隙針板結(jié)構(gòu)在納秒脈沖過電壓(86 kV)作用下的光學特性. 在實驗(見圖3(a))和后續(xù)的計算中[42]均發(fā)現(xiàn), 如果電壓上升時間明顯短于電場屏蔽效應的特征時間尺度, 局域場就會高于空氣電離閾值, 引起獨特的倒錐形放電結(jié)構(gòu).

(a) Measured discharge morphology

(b) Calculated distribution of secondary positive system(after direct Abel transformation)圖3 修正的COMSOL流注模型復現(xiàn)過電壓放電實驗Fig. 3 Overvoltage discharge experiment based on corrected COMSOL model

在與實驗完全相同的條件下, 本文試圖利用標準有限元模型計算該倒錐形放電結(jié)構(gòu), 模型因為未知(黑盒子)原因不收斂, 未能完成計算. 本文僅展示了利用修正后的COMSOL Multiphysics等離子體模型的計算結(jié)果.

在實驗中, 光學輻射的主要來源是氮氣分子的第二正帶系N2(C3Πu), 為此利用模型計算了該組分密度; 為了與實驗測量結(jié)果直接對照, 對計算結(jié)果進行了Abel變換以考慮放電區(qū)域的厚度. 將變換后的計算結(jié)果在0～1之間進行歸一化, 如圖3(b)所示. 實驗和計算均捕捉到了陽極和陰極區(qū)域的局域強輻射, 以及在大約Z=1 cm處的最大放電直徑1.2 cm, 實驗結(jié)果與計算結(jié)果對比吻合良好.

3 結(jié)論

針對商用有限元軟件模擬流注放電時存在的數(shù)值問題, 推導了基于流線擴散人工穩(wěn)定技術的修正項, 并在COMSOL Multiphysics軟件中進行了測試. 計算發(fā)現(xiàn), 采用標準有限元方法計算的針板流注結(jié)果存在顯著的數(shù)值振蕩和偽擴散流注問題, 而引入修正后的模型克服了這些數(shù)值問題, 能夠成功復現(xiàn)1 cm氣隙、 13 kV電壓條件下大氣壓空氣流注基準算例. 進一步計算了1.6 cm氣隙、 86 kV過電壓大氣壓空氣放電形態(tài)， 并與實驗結(jié)果進行了對比驗證, 證明了模型本身的可靠性. 計算與驗證結(jié)果表明, 以弱形式偏微分方程的方式將人工穩(wěn)定項引入基于有限元(標準Galerkin方法)的商用軟件中, 能夠有效克服放電計算中的數(shù)值振蕩和流注偽擴散問題, 確保等離子體計算準確性和可重復性.

致謝本文受到國家自然科學基金(51907204, 51790511, 91941105, 91941301)支持.