废弃物尤其是放射性高危废弃物的处理已成为可持续发展的重要课题。相较于传统的处理方法，等离子体技术因能源利用率高、处理效果突出成为了近年的研究和应用热点^[1]，其中尤以无害化处理放射性固体核废料为核心^[2]。该技术目前在国外的成功应用范例是等离子体熔融玻璃化（PM/V），其中，利用等离子体的高温特性，通过加入合适的添加剂，将硅酸盐、二氧化硅等介质加热至熔融状态，从而将有害的核固废固定在晶体材料中，形成低浸出率的玻璃体，实现无害化处理^[3]。

由于等离子发生装置的结构特殊，高温–电磁相互耦合等特点，其实验研究尤其是电弧区域的数据测试等技术要求高、难度大，甚至有些参数难以测量，因此,模拟计算技术成为研究发生装置内部温度–电磁分布行为的主要手段。例如：在非转移弧领域较早是由1989年Scott等^[4]进行的数值模拟研究，通过环路积分求解电磁场方程，进而获得流场分布；随后，Murphy等^[5]采用2维模型进行了等离子体射流及有机物裂解的数值模拟；李和平等^[6]也分析了非转移弧的3维计算结果。又由于其中涉及的学科跨度较大、专业知识较多，针对不同的模型，采用的假设不同，得到的计算方法并不通用，例如：Trelles^[7]采用基于亚网格尺度的有限元方法模拟计算了3维非定常等离子炬，以及关于洛仑兹力对流动影响程度的讨论^[8]等。鉴于此，针对直流电弧层流等离子体炬，将流场与电磁场结合建立2维轴对称的磁流体动力学（MHD）计算模型，对发生器内部电弧–等离子体的分布特征进行了模拟，得到了与实际相符的物理场分布，为等离子炬装置的工业应用设计优化提供了方法基础。

1 数学模型 1.1 基本假设

主要探究在宏观条件下，装置内部的热–电分布特性，因此对模拟计算的电弧–等离子体的发生过程做出了合理假设与简化，具体包括：

1）在弧柱区采用局域热力学平衡假设，默认该区域内各种质量的粒子碰撞频率足够高，使其速度分布满足Maxwell速率分布，对应的特征温度相等。

2）等离子体为光学薄。

3）等离子体的密度、比热、黏度、热导率和电导率等热力学参数和输运系数都仅为温度的函数。

4）忽略重力和黏性耗散。

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

相比较目前采用有限元、双温模型及非平衡状态下的物理数学模型等方法来延伸输运方程计算范围以及提高计算准确度，其相应的计算成本要高得多。因此，根据以上假设，将磁流体动力学方程组中的诸多控制方程进行大幅简化，并更加合理地选择部分实验难以获得或精度难以保障的参数项，从而更加高效地服务于工业设计优化。

1.2 控制方程

使用定常、2维轴对称条件，通过对纳维–斯托克斯方程组、Maxwell方程组各项进行增减，得到以下磁流体动力学方程组（MHD）。

1）纳维–斯托克斯方程组

$ \frac{\partial }{{\partial {\textit{z}}}}\left( {\rho {v_\textit{z}}} \right) + \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\rho {v_{{r}}}} \right) = 0 $

(1)

$ \begin{aligned}[b] \dfrac{\partial }{{\partial {\textit{z}}}}\left( {\rho {v_{{{{r}}}}}{v_{\textit{z}}}} \right) +& \dfrac{1}{r}\dfrac{\partial }{{\partial r}}\left( {r\rho {v^2_{{{{r}}}}}} \right) = - \dfrac{{\partial p}}{{\partial r}} + \dfrac{1}{r}\dfrac{\partial }{{\partial r}}\bigg(2r\mu \dfrac{{\partial {v_{{r}}}}}{{\partial r}}\bigg) +\\ & \dfrac{\partial }{{\partial {\textit{z}}}}\left[ {\mu \left( {\dfrac{{\partial {v_{\textit{z}}}}}{{\partial r}} + \dfrac{{\partial {v_r}}}{{\partial {\textit{z}}}}} \right)} \right]- 2\mu \dfrac{{{v_r}}}{{{r^2}}} - {{j}_r}{{B}_\theta } \end{aligned} $

(2)

$ \begin{aligned}[b] \dfrac{\partial }{{\partial {\textit{z}}}}\left( {\rho {v^2_{\textit{z}}}} \right) +& \dfrac{1}{r}\dfrac{\partial }{{\partial r}}\left( {r\rho {v_{{{{r}}}}}{v_{\rm{{\textit{z}}}}}} \right) = - \dfrac{{\partial p}}{{\partial {\textit{z}}}} + \dfrac{\partial }{{\partial {\textit{z}}}}\bigg(2\mu \dfrac{{\partial {v_{\textit{z}}}}}{{\partial {\textit{z}}}}\bigg) + \\ & \dfrac{1}{r}\dfrac{\partial }{{\partial r}}\left[ {r\mu \left( {\dfrac{{\partial {v_{\textit{z}}}}}{{\partial r}} + \dfrac{{\partial {v_{{r}}}}}{{\partial {\textit{z}}}}} \right)} \right] + {{j}_{{r}}}{{B}_\theta } \end{aligned} $

(3)

$ \begin{aligned}[b] &\dfrac{\partial }{{\partial {\textit{z}}}}\left( {\rho {v_{\textit{z}}}{C_{\rm{p}}}T} \right) + \dfrac{1}{r}\dfrac{\partial }{{\partial r}}\left( {r\rho {v_{{r}}}{C_{\rm{p}}}T} \right) = \dfrac{\partial }{{\partial {\textit{z}}}}\left( {k\dfrac{{\partial T}}{{\partial {\textit{z}}}}} \right) + \\ &\dfrac{1}{r}\dfrac{\partial }{{\partial r}}\left( {rk\dfrac{{\partial T}}{{\partial r}}} \right) - 4\text{π} {\varepsilon _{\rm{N}}} + \dfrac{{{{{j}}_{{r}}}^2 + {{{j}}_{\textit{z}}}^2}}{\sigma } + \dfrac{5}{2}\dfrac{{{k_{\rm{e}}}}}{e}\left( {\dfrac{{{{{j}}_{\textit{z}}}}}{{{C_{\rm{p}}}}}\dfrac{{\partial h}}{{\partial {\textit{z}}}} + \dfrac{{{{{j}}_{{r}}}}}{{{C_{\rm{p}}}}}\dfrac{{\partial h}}{{\partial r}}} \right) \end{aligned} $

(4)

式中， ${{{B}}_\theta }$ 为磁感应强度的周向分量， ${{{j}}_{\textit{z}}}$ 、 ${{{j}}_{{r}}}$ 分别为电流的轴向与径向分量， ${k_{\rm{e}}}$ 为工质气体的导热系数， ${C_{\rm{p}}}$ 为定压比热， $\sigma $ 为电导率， ${\varepsilon _{\rm{N}}}$ 为净辐射系数。

式（2）、（3）添加了洛伦兹力源项，以满足带电粒子在自磁场中受到洛伦兹力的“箍缩”对动量方程求解的影响；式（4）添加了焦耳热、热辐射、电子焓输运项，以满足等离子体的电流热效应、光学薄且辐射不可忽略及碰撞支配能量传递的特殊性质。

2）Maxwell方程组

电流连续性方程：

$ \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\sigma \frac{{\partial \varphi }}{{\partial r}}} \right) + \frac{\partial }{{\partial {\textit{z}}}}\left( {\sigma \frac{{\partial \varphi }}{{\partial {\textit{z}}}}} \right) = 0 $

(5)

径向磁矢量势方程：

$ \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {{{A}}_{\textit{z}}}}}{{\partial r}}} \right) + \frac{\partial }{{\partial {\textit{z}}}}\left( {\frac{{\partial {{{A}}_{\textit{z}}}}}{{\partial {\textit{z}}}}} \right) = - {\mu _0}{{{j}}_{\textit{z}}} $

(6)

轴向磁矢量势方程：

$ \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {{{A}}_r}}}{{\partial r}}} \right) + \frac{\partial }{{\partial {\textit{z}}}}\left( {\frac{{\partial {{{A}}_r}}}{{\partial {\textit{z}}}}} \right) = - {\mu _0}{{{j}}_r} $

(7)

其中，欧姆定律：

$ {{{j}}_{\textit{z}}} = - \sigma \frac{{\partial \varphi }}{{\partial {\textit{z}}}},\;{{{j}}_r} = - \sigma \frac{{\partial \varphi }}{{\partial r}}{\text{。}} $

磁矢量势与磁感应强度转化关系：

${{{B}}_\theta } = \frac{{\partial {{{A}}_r}}}{{\partial {\textit{z}}}} - \frac{{\partial {{{A}}_{\textit{z}}}}}{{\partial r}}{\text{。}} $

式中， $\varphi $ 为电势， ${{{A}}_{\textit{z}}}$ 、 ${{{A}}_r}$ 分别为磁矢量势的轴向与径向分量。

该方程组中电势 $\varphi $ 及磁矢量势 ${{{A}}_i}$ 为自定义标量，需要使用UDS功能添加方程。

以氩气单质为工质气体，用到的输运性质和热力学性质由Murphy等^[9]给出。

2 计算模型 2.1 网格划分及边界条件

使用的发生器结构尺寸参考中科院力学所潘文霞课题组设计的等离子体实验装置^[10-11]，并做了适当简化，整体尺寸相对较大，而采用的气体流量较小（入口氩气流量为1.7 g/s），故假设流场处于层流状态，并进行模拟计算， $\gamma = \varphi \partial \varphi /\partial {{n}}{\text{。}}$

简化模型及网格如图1所示。

在Fluent软件中采用2维轴对称模型，ABCDEFGHIJKA为计算域，其中，DE为对称轴，AK为工质气体入口，EF为等离子体出口，ABC为阴极铜座边壁，CD为嵌有阴极金属的电子发射端，FGHI为阳极边壁，IJK为电势非零结构边界。边界条件（表1）中通量为零项均使用如下符号。

采用200 A的弧电流，根据文献[12]可知，阴极附着区域半径为0.9 mm。同时由文献[13]，阴极表面的电流分布形式为抛物线型或是指数型，对计算域的流场分布没有影响，只是对近阴极区的电势电流分布有影响，因此，在阴极端面采用抛物线型分布（式（8）），此边界条件也作为UDF写入。由此开始计算等离子体炬内物理场分布。

$ {{{j}}_{\rm{c}}} = {{{j}}_{\max }}\left( {1 - b{r^2}} \right) $

(8)

2.2 计算过程

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

通过第1节能量方程的形式可以看出，随温度变化的物理量同时也会随着每次迭代的温度场分布不同而变化，进而对其余各方程的求解造成不同精度的影响。因此从温度来设想1次迭代过程如下：由初始流场温度计算出电导率，从而得到此次迭代中的电势分布；由电势和电导率计算得到各方向电流密度；电流密度再作为源项计算得到磁感应强度分布；磁感强再作为源项计算动量方程得到速度场；同时，电流密度还作为能量方程源项，与得到的速度场（对流项）再次计算温度场，从而形成1次迭代过程。由该过程可以发现：虽然电磁场与流场耦合计算，但是以电势为基础的方程组的求解相对独立，温度作为影响因素仅对电导率的计算有影响；此外，各方程源项及标量方程的添加，容易引起迭代计算的震荡或者发散，其中关系最为密切的是能量方程中所添加的3个源项。因此，为了得到稳定且较为准确的物理场分布，一方面，需要保证高质量的网格及合理的方程松弛因子；另一方面，与能量方程相关的各项参数采集及形式变化需要经过仔细斟酌。

3 计算结果分析

使用Fluent软件对以上模型进行计算，得到的物理场分布如下。

3.1 速度分布

由于直流等离子体发生器应用的领域越来越多，配合不同的工况，需要调节不同的工艺参数以达到合适的等离子体状态。例如，在给定弧电流的条件下，随着输入工质气体质量流量的增大，流动会由层流向湍流转变。层流等离子体射流平稳的输出特点降低了与环境气体的卷吸混合，进而保障了射流质量；而湍流等离子体的强卷吸特性可以有效地增强处理有机废物时的解离效果。

图2为速度–压强分布。由图2可知，速度最大值为800 m/s，发生位置为阳极筒体中压力最小的区域。由于阴极附近的电流密度最大，热效应最明显，输入的氩气被急速加热，在洛仑兹力及气体的热膨胀作用下形成加速射流。发生器存在阶梯扩张结构，因此在输出等离子体前会经历一次膨胀作用，更有利于形成层流等离子体。

3.2 温度分布

关于等离子体熔融玻璃化技术，最关注的是出口射流的温度分布，因为后续工段使用的是硅酸盐或是二氧化硅，要满足使其融化而不至于气化的温度条件，就要合理控制发生器内部的温度分布。图3为本次模拟的温度场分布。

由图3可知，等离子体温度最高的部位集中于轴向近阴极区域，中心温度可达24 000 K以上，发生器内氩气流量小，本身近似层流流动，因此电弧行为较为平缓，等离子体整体呈现平稳均匀的流动状态，因此，高温区域在轴线方向上分布范围较广。

将轴向距离为14 mm处截面M–M的温度分布值提取出，绘制温度与径向距离的对数关系曲线，如图4所示。由图4可知，温度梯度较大的半径位置刚好是在等离子体与冷气层的交界面处。这个特点对应了等离子体内部能量传递是以粒子间的碰撞为主，其中，主要是电子与重粒子的碰撞，而电子在其中的运动规律宏观表现即为电流在发生器内部的分布。

3.3 电流密度分布

图5为发生器内部的电流密度分布情况。显然，阴极附近的电流密度最为集中，沿轴线方向逐渐减低，最后终至止于某处。在本文计算条件下，电流密度最后与阳极壁面接触于轴线23 mm处，可以理解为阳极弧根附着位置为此处。由于将阳极整个壁面的电势值都设定为0，区别于传统固定位置的虚拟阳极法，随着使用弧电流大小及氩气流量的不同，弧根的轴线距离将发生变化，进而改变了等离子体高温区域的延伸范围。为了进一步说明阳极弧根位置与等离子体高温区域的关系，将两者的分布图放在一起，如图6所示。

由图6红色标注部分可知，高温区域终止位置明显与弧根分布区域重合，这一结论也证明了前文在构建能量方程源项，需要着重考虑电流的热效应，即焦耳热项。

3.4 电势分布

在描述发生器性能优劣时经常会使用伏安特性曲线，数值模拟中一般给定总的弧电流大小来计算电势，而考虑到电势分布，就不得不提到鞘层效应。通常来说，鞘层具有几何尺寸小、压降高的特点，并且由于其本身处于近电极区，自身为非局域热力学平衡状态（NLTE），不仅在工质气体热力学性质、输运性质需要重新选择，还要在磁流体力学方程组中对能量方程强调电子–重粒子的能量传递过程。除此之外，为了计算到稳定且准确的耦合鞘层效应的流体模型，还需要将网格细化到小于鞘层尺寸的程度。显然，仅为了得到宏观物理场分布，用以调节等离子体射流特性时，只需要对鞘层模型进行适当简化，根据文献[14]，使用简化鞘层相关参数的设定，并采用电势补偿的方法。本次计算电势分布如图7所示。

通过图7可知，电势最大值为35 V，集中在阴极附近，并随着轴向距离的增大，电势迅速下降，在23 mm处基本稳定在较低电势值。整个腔体弧压降为35 V左右，根据Pan等^[15]的实验测定，在氩气质量流量为1.7 g/s左右时，弧压降在50 V左右，即补偿压降为15 V左右。在排除计算误差外，该差值体现在阴极鞘层内的弧压降，根据文献[16]的研究结果，当弧电流为200 A时，鞘层的电势差约为12 V，与计算结果较为接近。

提取轴线上电势分布数据，并与电流分布对比，如图8所示。由图8可知：电流密度与电势梯度趋近平缓的轴向位置相同；与电流–温度分布图比较可以发现，电流密度分布长度，即电弧长度影响着工质气的加热距离，进而影响射流的温度分布；同时，弧根附着位置也会影响到整个发生器的弧压降，从而影响到设备效率。因此，不难发现，弧根附着位置与停留时间的长短会涉及到阳极烧蚀的问题，所以，寻找合适的工艺参数，使得兼顾电弧的脉动行为与射流热效率是工业化高效生产的关键。

4 结　论

使用Fluent软件对直流非转移弧等离子体发生器进行了2维轴对称数值模拟。得到以下结论：

1) 最高速度产生于阳极筒体内压力最小区域，并随着轴线距离平稳而缓慢地降低。

2) 层流条件下，高温区域分布均匀且较长，与阳极弧根附着轴向尺度有关，等离子主体与冷气界面上温度梯度较大。

3) 电流电势基本符合阳极弧根弥散型电弧分布特点，结合国内外实验研究，鞘层结构压降与未计算压降大小基本相近。

因此，文中所涉及到的基础数据（气体热力学、输运性质、净辐射系数等）与计算模型（MHD模型、简化鞘层方法等）可以在其他相似结构的直流非转移弧等离子炬的数值模拟中使用。此外，由结果分析可知，电流密度分布尤其是阳极弧根附着位置会在很大程度上影响温度场及电势分布。后续将会开展探究影响电流分布及热效率因素的工作。