朱益飛, 陳賢聰

(1. 空軍工程大學(xué)航空工程學(xué)院等離子體動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室, 陜西西安 710038; 2.西安交通大學(xué)機(jī)械工程學(xué)院, 陜西西安 710049)

引 言

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

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

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

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

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

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

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

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

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

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

在納秒放電時(shí)間尺度下, 可以假設(shè)離子和中性粒子不運(yùn)動(dòng), 僅求解不包含輸運(yùn)項(xiàng)的反應(yīng)方程即可

式中,ni為離子密度,Di為離子擴(kuò)散率,Si為化學(xué)反應(yīng)主導(dǎo)的離子產(chǎn)生率,Sph為光電離主導(dǎo)的離子產(chǎn)生率.

電子的輸運(yùn)方程為

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

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

對(duì)上式前兩項(xiàng)進(jìn)行分部積分, 并簡(jiǎn)化可得

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

(1)

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

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

圖1(a)和(c)對(duì)比了使用COMSOL Multiphysics在兩種網(wǎng)格條件下計(jì)算的流注軸向電子密度分布(藍(lán)線)和基準(zhǔn)算例中的軸向流注通道電子密度(紅線). 在粗網(wǎng)格下, COMSOL Multiphysics計(jì)算所得流注能夠復(fù)現(xiàn)基準(zhǔn)算例流注傳播速度, 如圖1(b)所示, 但是整個(gè)放電通道內(nèi)出現(xiàn)了圖1(a)中劇烈的非物理數(shù)值振蕩, 計(jì)算結(jié)果無法分析理解. 如果僅僅加密網(wǎng)格, 如圖1(c)所示, 計(jì)算軸向電子密度不再產(chǎn)生數(shù)值振蕩, 但是流注傳播距離遠(yuǎn)遠(yuǎn)落后于基準(zhǔn)算例, 圖1(d)表明, COMSOL Multiphysics標(biāo)準(zhǔn)有限元方法無法捕捉到基準(zhǔn)算例中的電子雪崩-流注轉(zhuǎn)捩, 計(jì)算結(jié)果呈擴(kuò)散流注狀. 這兩個(gè)案例清楚地表明, 使用COMSOL Multiphysics默認(rèn)的數(shù)值格式和計(jì)算方法不可能正確計(jì)算大氣壓流注. 本文使用COMSOL Multiphysics內(nèi)置的等離子體模塊進(jìn)行測(cè)試, 也得到了同樣的錯(cuò)誤結(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)定性, 源于計(jì)算網(wǎng)格與算法不匹配, 有限元中的標(biāo)準(zhǔn)Galerkin方法與有限體積/有限差分中的標(biāo)準(zhǔn)SG方法一樣, 在高氣壓條件下, 計(jì)算等離子體流注的準(zhǔn)確性會(huì)顯著下降. 假設(shè)電子擴(kuò)散系數(shù)符合Einstein關(guān)系, 那么為了保證計(jì)算結(jié)果的準(zhǔn)確性, 須保證以下關(guān)系[17]

(2)

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

(3)

式中, 各項(xiàng)的含義如下

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

2 基準(zhǔn)驗(yàn)證與實(shí)驗(yàn)驗(yàn)證

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

2.1 計(jì)算基準(zhǔn)驗(yàn)證

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

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

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

(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流注模型復(fù)現(xiàn)Kulikovsky基準(zhǔn)流注算例[16]Fig. 2 Results comparison between corrected COMSOL model and Ref[16]

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

2.2 實(shí)驗(yàn)結(jié)果驗(yàn)證

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

(a) Measured discharge morphology

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

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

在實(shí)驗(yàn)中, 光學(xué)輻射的主要來源是氮?dú)夥肿拥牡诙龓礜2(C3Πu), 為此利用模型計(jì)算了該組分密度; 為了與實(shí)驗(yàn)測(cè)量結(jié)果直接對(duì)照, 對(duì)計(jì)算結(jié)果進(jìn)行了Abel變換以考慮放電區(qū)域的厚度. 將變換后的計(jì)算結(jié)果在0～1之間進(jìn)行歸一化, 如圖3(b)所示. 實(shí)驗(yàn)和計(jì)算均捕捉到了陽極和陰極區(qū)域的局域強(qiáng)輻射, 以及在大約Z=1 cm處的最大放電直徑1.2 cm, 實(shí)驗(yàn)結(jié)果與計(jì)算結(jié)果對(duì)比吻合良好.

3 結(jié)論

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

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