耦合电极直流等离子体物理场数值模拟研究

武劭恂, 司徒达志, 张子炜, 李 晴, 黄卫星

(1. 四川大学 化学工程学院, 四川 成都 610065; 2. 中国广核集团中广核研究院有限公司, 广东 深圳 518124)

1 引 言

国内近几年对固体垃圾问题，包括高危固体废弃物的重视，使探索出一种高效、清洁的处理方法变得日益重要。而等离子体技术具有能源利用率高、处理效果突出的特点，能够满足这一需求。目前在国外该技术已应用到包括城市生活、医院卫生、工业废置等诸多领域[1-2]。而该技术在应用中却存在一个问题，即由于等离子体温度高于电极金属的熔点，导致设备寿命减小，工作效率降低[3]。所以发展一种简便可靠的耦合电极-电弧的研究方法则成为关键。

由于等离子体发生装置尺寸小、结构复杂，在实验或生产操作中难以对其电极状态进行观测。除此之外，其具有的高温-电磁相互耦合的特点更使得数据测量工作难以进行。因此，数值模拟计算成为解决该问题的主要研究方法[4-6]。但同时，由于人们更多关注等离子体主体内的物理场分布，而忽略了电极本身热电子发射机制带来的热效应，因此在传统计算结果中对于电极温度的计算需要改进。有鉴于此，且考虑转移弧与非转移弧发生原理相似，因此本文针对自由燃烧电弧等离子体，将电极与电弧区耦合建立二维轴对称的磁流体动力学 (magnetohydrodynamic, MHD) 计算模型，并通过合理的设计方法，通过求解电子连续方程与附加面热源将电弧对电极的热效应、流体域的电导率修正包含在内，得到了与实验测量相近的温度场分布。该模型的建立可以为工业用直流电弧等离子体炬的设计改良提供方法基础。

2 数学模型

2.1 基本假设

研究宏观条件下发生装置内部的温度场分布，因此通过下面的合理假设对复杂的电弧模型进行简化，具体包括：

1) 在电弧区域采用局域热平衡 (local thermodynamic equilibrium, LTE) 假设；

2) 等离子体为光学薄；

3) 等离子体的密度、比热、速度、黏度、热导率、电导率、双极扩散系数、净辐射系数、气体粒子数密度分布等热力学参数和输运系数都仅为温度的函数；

4) 忽略重力和黏性耗散；

5) 由于等离子体温度高，对应的声速也高，使得当地流动马赫数并不高，非特别说明时等离子体流动是远低于声速的，视为不可压缩流动，可忽略压力功；

以上假设可以将耦合模型中控制方程的诸多项进行简化，并能够合理地选择需要的变物性函数关系。

2.2 控制方程

本文计算模型以自由燃烧电弧为基础，使用定常、二维轴对称条件，通过对纳维-斯托克斯方程组、Maxwell 方程组各项进行增减，并补充近电极区的电导率修正以及电极表面热源，得到以下耦合电极的磁流体动力学方程组(MHD)。

(1)纳维-斯托克斯方程组

式中zv 、rv 分别为速度的轴向与径向分量；Bθ为磁感应强度的周向分量；zj 、rj 分别为电流的轴向与径向分量；T 为温度；ek 为工质气体的导热系数；pc 为定压比热；σ 为电导率；Nε 为净辐射系数；kB为玻尔茲曼常数；e 为电子电荷。物理量单位均使用国际单位制。

式(2)、(3)为动量方程，添加了洛仑兹力项，以实现磁场中带电粒子受到洛仑兹力作用而产生的“箍缩”现象，这一现象发生在自由燃烧电弧中时，由于电弧受力向轴线方向压缩，导致电流流通面积变小，电流密度增大，其产生的焦耳热也会增大，因此不可忽略。而发生在非转移弧时，由于存在固壁约束，电流的流通面积在气流与洛仑兹力的共同作用下减少，此时洛仑兹力项可忽略。

(2) Maxwell 方程组

图1 氩等离子子体净辐射系数 Fig.1 Profile of net emissionn coefficient of argon plasma

其中 σeff为有效电导率， n0为中性性粒子数密度， nT为总粒子数， μe为电子迁移率率，σ 为热平衡状态下的电导率。值得注意的是，除了电电子迁移率 μe以外，粒子数密度等都可以看作温度度的函数，而电子迁移率，根据其定义

其中，em 为电子质量，ev 为粒子碰撞频率，thv 为电子热运动速度，eiQ 为电子-离子碰撞截面，eaQ 为电子-原子碰撞截面，这些参数并没有温度函数关系，意味电子迁移率eμ 需要用额外的方法去求解，本文使用了Bolsig 程序[12]求解玻尔兹曼方程得到了系列电子迁移率eμ ，用来修正电导率。

除了宏观等离子体对电极的热传导、热辐射外，电弧的行为也影响到电极的能量传递过程。阴极的过程具体表现为阴极金属发射的电子具有一定的电子焓热，使阴极表面温度降低。同时阴极表面的电离层产生的离子，将会在较大的电场作用下向阴极表面运动从而带来一定的离子焓热，使阴极表面温度升高，除此之外还有本身的热辐射，几种因素共同作用，作为电极表面的额外能量输入。阳极的过程表现为电子运动到阳极表面释放焓热，存在热辐射，而没有离子焓。电极的面热源形式如下[13]：

其中，ej 为近阴极区的电子电流，cφ 为阴极的逸出功，ij 为近阴极区的离子电流，iV 为工质气体一次电离电势，ε 为表面发射系数，α 为斯特藩-玻尔兹曼常数。

在具体计算中，由于离子电流无法提前得知，需要根据电流连续性条件，将阴极金属发射产生的热电子电流与弧柱区求解Maxwell 方程组得到的电弧电流做比较，继而通过离子电流补全守恒关系。热电子发射方程(Richardson 方程)为

其中， A 为热发射常数， kB为玻尔兹曼常数。离子电流 ji计算方法为，当弧区电流大于发射电流时( jn＞ jR)， je= jR， ji= jn- jR；当弧区电流小于发射电流时( jn＜ jR)， je= jn， ji= 0。

但由于在Fluent 软件中，无法进行面热源的设置，式(9)、(10)无法直接使用，因此采用另外一种，类似激光加热金属表面的处理方法[14]对其进行等效转化，从而能够正常加载。

2.3 模型归纳

2.2 节对该耦合模型的内容以及设置方法进行了描述，在此，为了方便与Fluent 软件程序相接，将上述方程组归纳为统一的控制方程形式，即基本的NS 方程形式。

表1 项类归纳 Table 1 Components of equations

根据以上归纳，即可使用户自定主标量(user-defined scalar，UDS)定义新的控制方程，使用用户自定义函数(user-defined function，UDF)写入物性函数、边界条件等。

3 计算模型

3.1 模型划分及边界条件

本文采用二维轴对称模型，图2 为模型结构，为了显示阴极平台，图示已作适当放缩，取半径35 mm，长40 mm 的计算域，其中包含阴极固体域、流体域、阳极固体域3 个部分。阴极固体域主体半径1.6 mm，尖端平台为30 度半锥角，平台半径为0.1 mm，材料为钍钨合金；流体域使用氩气工质气体；阳极固体域为大平板，材料为铜。由于本文只探究电极温度分布，对于温度超过金属熔点的情况不添加熔化-凝固模型，一律默认水冷材料，同时水冷材料不会影响电弧特性[9]。

网格建立过程中，采取先疏后密，并不断调节亚松弛因子，每阶段计算收敛后，使用自适应网格，将温度梯度较大区域以及电弧主要分布区域网格加密，不断调整计算最终达到要求。

图2 耦合电极自由燃烧电弧计算模型 Fig.2 Calculation model for coupling electrode free combustion arc

本文采用200 A 弧电流，电流连续方程与磁矢量势方程在全计算域计算，电子连续方程只在流体域计算，使用的工质气热力学参数与输运参数分别来自MURPHY[15]与张晓宁[16]。

所用程序使用Visual Studio2017 编写，Fluent18.1 求解，除了直接得到上述标量值，对计算域中总(分)电流值、焦耳热值、热辐射值、电场强度等黑箱值均单独分配存储空间UDM 以便分析。

表2 边界条件 Table 2 Boundary conditions

3.2 计算过程

该模型需要将电磁场与一般流场涉及到的物理模型耦合计算，需要使用Fluent 软件的UDF 与UDS功能，用以引入某些物性参数随温度变化的函数，并添加自定义标量电势φ 、磁矢量势的轴向与径向分量zA 、rA，以及电子数密度en 对应的标量方程。

通过上文对控制方程组的描述可以得知，除了电子连续方程较为特殊之外，其余方程求解中需要使用到的物性参数均与温度分布密切相关，因此从温度来设想一次迭代过程如下：由初始流场温度计算出局域热平衡电导率σ ，同时得到当前温度下的双极扩散系数 AD 、局域热平衡电子数密度eqn 用以求解式(8)电子连续方程，得到当前温度下的电子数密度en ，该温度下还可求解出1σ 参与的电流连续方程得到电势，即得到电场强度分布，以上计算结果可以得到有效电导率effσ ，于是其他过程均以有效电导率effσ 为基础继续计算。可以发现，虽然电磁场与流场耦合计算，但是电子方程和以电势为基础的方程组的求解相对独立，温度作为影响因素仅对电导率的计算有影响，而电子方程的求解反过来又会影响电势方程的求解，所以在大的耦合模型下，这两个方程又作为另一对耦合形式出现。此外，各方程添加的源项，极易引起迭代计算的震荡或者发散，因此为了得到稳定且较为准确的物理场分布，一方面需要保证高质量的网格以及合理的亚松弛因子，另一方面各项参数采集以及形式变化需要经过仔细斟酌。

4 计算结果分析

4.1 速度-压力分布

气体进口处速度较低，但由于阴极端部电流密度很大，工质气被迅速加热电离后，等离子体在洛仑兹力的作用下压缩，热量集中，从而使等离子体加速形成射流。如图3、4，在阴极附近，即温度最高的区域形成了一个由于热膨胀造成的高压区与阳极轴线附近形成了动能-势能转化高压区。

图3 速度分布 Fig.3 Map of speed distribution

图4 压力分布 Fig.4 Map of pressure distribution

4.2 温度分布

电极耦合模型的固体域温度分布是预估阴极寿命的重要依据，因此在自由燃烧电弧模型的基础上需要补充阴极金属发射电子对流体域电弧计算的影响，以及阴极弧根的形位对阴极表面的附近热源的影响。图5 中该模型计算得到流体域最高温度位于阴极尖端附近，约为23 150 K，阴极最高温度位于阴极尖端，约为4 000 K，阳极最高温度位于轴线附近，且呈近二次项分布，约为850 K。温度分布范围与图6 中前人研究结果基本一致。可以发现，在不采用电极水冷或添加剂熔融扩散的前提下[13]，采用200 A 电流强度势必会造成阴极的烧蚀。由于在计算过程中，气体逐渐电离，电极温度逐渐升高，可以设想，在此计算模型的基础上，如果采用非稳态模型计算，当阴极温度升高至熔点时所用计算时间即为电极的工作寿命，继而可以加入水冷与熔化-凝固模型对实际设备进行较为准确的设计优化，这也是本文的出发点。

图5 计算温度分布 Fig.5 Calculated temperature distribution

图6 文献中温度分布[11] Fig.6 Temperature distribution results from literature[11]

4.3 电子数密度分布

此处的电子数密度为求解电子连续性方程得到，由于该标量方程涉及到的物性参数较多，数量级普遍较大，采集时存在一定误差，此处不再展示其余粒子数密度等分布图。

图7 电子数密度分布Fig.7 Map of electronic number density distribution

由图7 中计算结果可知，电子数密度分布规律与温度分布规律近似，并且16 000～23 000 K 区域对应电子数密度最大的区域，这是因为一次电离的氩离子在14 000 K 左右会发生二次电离，但在25 000 K 以下由二次电离产生的电子数与一次电离提供的电子存在数量级的差别，直到高于25 000 K 时，一次电离的氩原子数才会呈现明显下降。所以，虽然本文中计算温度最高达到了23 000 K，但是使用的粒子数相关物性仍采用一次电离数据，使用的氩气电离能也采用了一次电离能，计算电子数密度分布规律基本合理。

4.4 电势-电流密度分布

与电势、电流密切相关的参数是电导率，由以上计算得到的电子数密度对局域热平衡电导率进行修正，最终得到电势为13.6 V(图8)，电流密度最大为7×108A⋅m-2位于阴极尖端附近(图9)。由于本文关于近阴极区的计算未包含鞘层的电位降计算，所以采用电势补偿的方法[17]，由ZHOU 等[18]的研究，当电流为200 A 大小时，鞘层内电势差约为11.7 V，则补偿后电势差为25.3 V，与实验测量值26 V 较为接近。相较于一般电导率σ 为温度的单调函数，有效电导率effσ 在修正过程中需要使用电子迁移率eμ ，这一参数在很大程度上与电场强度E 大小相关。由电弧的结构可知，在近电极区的较小轴向距离内存在很高的电势差，这意味着该区域的场强较大，因此，场强分布非线性，导致除了气体完全未电离和完全电离两种情况外，修正后的电导率分布关于温度并不是单调关系，由此导致计算结果里电势的数值会低于一般自由燃烧电弧，而电弧主体温度值接近，其中差别最明显的就是在近阴极区的温度分布。

图9 电流密度密度分布Fig.9 Map of density distribution of current density

5 结 论

针对耦合电极的直流等离子体进行了二维轴对称数值模拟。主要结论如下：

1) 速度-压力变化关系与一般自由燃烧电弧相似，速度高值位于电弧中心，压力高值位于电弧两端；

2) 流体域的温度-电子数密度分布规律基本一致，符合氩气正常电离过程的数量变化规律；固体域温度结果合理，可以近似预测电极寿命；

3) 提出了针对Fluent 软件功能的计算模型，以及在实验数据缺乏下的近似处理方法。