摘要:以航天领域中研究再入飞行器热防护系统的感应耦合等离子体(inductively coupled plasma, ICP)风洞为研究对象, 通过流场-电磁场-化学场-热力场-湍流场多场耦合求解研究ICP风洞流场与电磁场的分布特性及其相互作用机理. 数值模拟中, 基于热化学非平衡等离子体磁流体动力学模型准确模拟了空气ICP的高频放电、焦耳加热、能量转化、粒子内能交换等过程, 通过多物理场耦合计算模拟得到了100 kW级ICP风洞内空气等离子体的电子温度、粒子数密度、洛伦兹力、焦耳加热率、速度、压强、电场强度的分布规律. 研究结果表明: 在感应线圈区靠近等离子体炬壁附近, 等离子体流动处于热力学非平衡状态; 洛伦兹力对感应线圈区空气粒子的动量传递和电子热运动起着控制作用.
关键词: 感应耦合等离子体/
流场与电磁场/
非平衡/
数值模拟

English Abstract

全文HTML

--> --> -->

再入飞行器完成太空任务重返地球大气层时, 由于其飞行速度很快, 在飞行器前方会形成很强的弓形激波, 激波层内气体被急剧压缩和加热形成温度高达几千摄氏度的等离子体气流, 飞行器头部材料将被等离子体气流急剧加热而发生分解反应, 因此, 若不在飞行器头部表面加载热防护材料, 飞行器很可能被等离子体气流损坏或烧毁.
为了开发耐高温且质量轻的飞行器热防护材料, 近年来, 很多国家建立了高温高焓风洞以开展飞行器防热材料烧蚀试验, 如电弧加热风洞、激波风洞、感应耦合等离子体(inductively coupled plasma, ICP)风洞. 由于ICP风洞能够产生连续、高纯度、接近真实飞行环境的高温、高焓气流, 因此它被广泛用于开发耐烧蚀的热防护材料^[1]、研究飞行器表面氮化反应机理^[2]、研制航天器新型薄膜减速伞材料^[3]等. ICP除了用于航天领域外, 也可以用于研究甲醛气体协同净化处理^[4]、超微金属粉末球化处理^[5]、微纳米颗粒制备^[6]等. 在进行上述工业应用研究时, 如果想对这些实验研究进行准确、细致地分析, 就必须先准确掌握ICP的流动特性与形成机理. 因此, 为了准确获得ICP的流动特性与形成机理, 我们需要对ICP的形成过程、高频放电特性以及电磁场与流场的作用机理等方面进行深入、细致的研究.
ICP风洞的系统结构通常由三部分组成: 进气部分、高频放电部分和试验部分, 这三部分的结构布局和组成单元如图1所示. ICP风洞的核心工作部分由高频电源、感应线圈和石英管组成, 由于等离子体流动率先形成于石英管内, 因此石英管通常又被称为等离子体炬. ICP的形成过程是一个复杂的物理、化学过程, 它涉及高频放电、电磁加热、趋肤效应、非平衡内能交换、气体电离/离解等过程. ICP的具体形成过程如下: 常温气体通过进气道(如图1右上角所示)进入到石英管前端进行充分预混合^[7], 当气体流经感应线圈下方时, 经电火花点火, 气体介质将在感应线圈产生的高频交变电磁场作用下发生感应放电, 此时由于趋肤效应和焦耳热效应的作用, 一部分电能将转化为分子热能. 随着气体温度不断升高, 气体分子将发生电离和离解反应并释放出大量的热, 进而形成由电子、原子、分子和阳离子组成的准电中性等离子体流动. 由于等离子体中存在大量的自由电子和活性阳离子, 它们在交变电磁场的作用下会产生感应磁场. 由于该感应磁场与感应线圈产生的外加磁场存在一定的相位差和耦合关系, 因此所生成的等离子体常被称为ICP. 由于ICP的形成过程较为复杂, 且流动温度也较高, 使用实验测量工具对其流动参数进行全面诊断难度很大, 因此开展数值模拟研究成为当前研究ICP的一种重要方法.
图 1 ICP风洞系统结构布局
Figure1. Schematic diagram of the ICP wind tunnel system.

自第一台ICP风洞问世以来^[8], 世界各国****纷纷开展了相关的实验与数值模拟研究工作^[9-18]. 1976年, Barnes和Nikdel^[19]使用一维电磁场模型和能量平衡方程联立求解的方法研究了氮气ICP的温度场和速度场的变化规律. Mostaghimi和Boulos ^[20]基于麦克斯韦方程组和毕奥-萨伐尔定律提出了一种二维磁矢势模型, 并成功地将该模型应用到了氩气ICP的数值模拟中. Punjabi等^[21,22]利用FLUENT软件对不同功率ICP炬内等离子体阻抗、吸收功率、耦合效率等进行了仿真分析, 研究发现: 随着ICP风洞输入功率的增大, 其等离子体火焰体积随之增大, 等离子体最高电子温度、最大速度和最大电子数密度均随之增大^[17,21]; 对于不同工作介质(如氮气、氧气、氩气和空气) ICP流动, 在相同的工况条件下(功率、频率、气压、体积流量均相同), 氧气ICP的等离子体火焰体积最小, 氩气、氮气和空气ICP的等离子体阻抗随输入功率增大而增大, 但氧气ICP的阻抗随功率增大而呈下降趋势^[21]; 此外, 无论以上何种气体介质, ICP流动的最大温度、速度和电子数密度均随电源放电频率的增大而减小^[23]. 但以上研究均是在局域热化学平衡假设条件下开展的, 只有当工作气压较高时, 该平衡假设才被认为是有效的. 因此, 当工作气压较低时, 将热化学非平衡模型引入数值模拟中是非常有必要的. 对于非平衡ICP研究, Lei等^[10]利用COMSOL软件对低压ICP流动进行了非平衡仿真, Stewart等^[24]使用二温度模型研究了氩气ICP的热非平衡放电过程, Munafò等 ^[25]、Zhang等^[26]和Lei等^[10]也构建了热非平衡和化学非平衡模型, 并对氩气ICP流动的非平衡特性进行了研究.
基于氩气热化学非平衡模型的发展, Degrez等^[27]采用空气化学反应动力学模型研究了空气电离、离解反应模型对流动速度、温度分布的影响, 研究发现, 由于气体粒子的扩散效应, 高频放电区氧原子和氮原子出现了分层现象. 然而, 该研究没有考虑热力学非平衡对其仿真结果的影响. Sumi等^[28]也采用双温度模型研究了ICP的热非平衡特性以及温度分布规律. Morsli和Proulx^[29]基于Dunn-Kang的32化学反应模型和双温度模型构建了空气ICP热化学非平衡磁流体力学模型, 并对超声速空气ICP的流动规律进行了研究, 研究表明空气化学反应模型对准确模拟空气粒子的组分浓度起着至关重要的作用. 此外, 本课题组之前通过对电磁场方程组进行简化、对电子输运系数计算方法进行完善、对四温度模型和粒子内能交换模型进行补充等, 也对非平衡ICP的仿真研究做了些许工作^[11,30-32]. 在前期研究中, 本课题组提出了一种收敛速度快、相对简单的焦耳热源模型, 该模型可以替代感应线圈区电磁场方程组的求解, 并且能较准确地模拟ICP风洞试验腔内等离子体气流的温度分布^[11]; 在此基础上, 我们还通过耦合求解远场电磁场和非平衡流体力学方程组研究了10 kW级氮气与空气ICP的流场差异及其热化学非平衡特性^[17], 并对10 kW ICP风洞内氮气离子体的超声速流动进行了瞬态模拟^[30]. 除此之外, 我们还对空气ICP的热化学非平衡特性与气压的非线性关系进行了研究, 研究发现当气压大于19 kPa时, 10 kW空气ICP流动趋于局域热力学平衡态^[31], 并发现当电导率的计算精度从一阶提高到三阶时, 仿真得到的气体最高温度会下降680 K, 所以最终整理提出了一种计算空气和氮气等离子体三阶精度电子输运系数的方法, 该方法在ICP仿真研究中具有很好的应用潜力^[32].
尽管上述研究工作对ICP数值模拟技术都起到了积极的促进作用, 然而目前尚缺少对高功率非平衡ICP电磁场与流场作用机理的深入研究, 关于洛伦兹力、焦耳加热率、电子数密度、电子温度等参数之间的相互影响关系目前尚不清晰. 此外, 我们前期研究主要是针对10 kW级中小功率ICP风洞展开的, 对于100 kW级高功率ICP风洞来说, 由于其输入功率较大, 其电磁场强度、气体温度和电子数密度必然将更高, 因此其电磁场、流场与热力场之间的耦合也将变得更加紧密、更加复杂. 而且, 由于100 kW级ICP风洞的结构参数、输入功率、放电频率与10 kW级ICP风洞存在诸多差异, 因此本文拟对高功率ICP风洞的流场与电磁场相互作用机理进行研究, 以揭示高功率ICP的高频放电规律和非平衡流动特性.
本文以100 kW级ICP风洞为研究对象, 基于流场-电磁场-热力场-化学场-湍流场多场耦合求解技术对ICP电磁场与流场的相互作用机理进行研究, 以揭示电磁场强度对流场参数的影响机理. 为了准确描述高温空气微观粒子的热运动现象, 本文采用含有最新电子振动弛豫时间算法的四温模型, 并使用三阶精度电导率和电子导热率进行仿真计算. 由于高频感应放电在线圈区发挥着重要作用, 因此通过求解麦克斯韦方程组来计算线圈电流产生的外加电场和ICP产生的感应磁场, 并使用11组分空气、32种化学反应来模拟高温空气粒子的电离、离解与电量交换等化学反应过程^[33,34]. 在进行空气ICP仿真时, Park等^[35]所提出的化学反应动力学模型常被用来模拟空气的电离和离解过程, 但Park等并未考虑N₂, O₂和NO分子的副电离反应^[33](如: O₂+N₂$\rightleftharpoons$NO+NO⁺+e^–, N₂+NO$\rightleftharpoons$N₂+NO⁺+e^–, O₂+NO$\rightleftharpoons$O₂+NO⁺+e^–), 由于空气在温度超过4000 K和9000 K时将很容易发生离解和电离反应, 故为了准确模拟高温空气的化学反应过程, 本文在数值模拟中将考虑这些副电离反应. 最终, 基于上述模型与二维黏性可压缩Navier-Stokes方程组的耦合求解来对ICP的形成过程与感应放电机理进行描述.

本文所研究的100 kW级ICP风洞几何形状与计算网格如图2所示. 图2(a)为远场电磁场和ICP炬流场的计算网格, 电磁场的计算域为–120 mm ≤ x ≤ 360 mm, 0 ≤ y ≤ 206 mm, 由101 × 66个网格节点组成, 其中ICP炬(ICP torch)的流场网格由101 × 38个网格节点组成, 为了准确计算线圈附近的电磁场与进气口处气流速度, 计算网格在感应线圈区和进气口处进行了加密处理. 图2(b)为ICP炬的结构示意图, 感应线圈围绕ICP炬管外壁面缠绕三圈, 线圈的内径和外径分别为47 mm和55 mm, 起始线圈的中心坐标为(x, y) = (51, 51) mm, 相邻感应线圈的中心线间隔为16.5 mm.
图 2 ICP炬计算网格和几何结构　(a) 电磁场与流场计算网格; (b) 等离子体炬几何结构
Figure2. Computational mesh and geometry of the inductively coupled plasma torch:　(a) Computational mesh of electromagnetic- and flow-field; (b) geometric construction of the ICP torch.

ICP风洞正常工作时, 感应线圈上通有高频交变电流, 电流频率为1.78 MHz, 工作气体如空气或氮气从壁面附近的狭窄开口注入, 经电火花高能点火后, 工作气体开始发生离解和电离反应并产生自由电子. 此后, 带电粒子在交变电流产生的交变电磁场作用下进一步发生剧烈的碰撞电离、副电离等反应, 并释放出大量能量而形成温度高达10000 K的ICP气流, 最终该等离子体气流从ICP炬出口流入大空间试验腔内进行飞行器气动性能测试、热防护材料烧蚀试验等实验研究.
本文所开展的仿真算例工况条件如下: 面板输入功率P = 90 kW, 质量流率m = 1.8 g/s, 工作气压p = 10 kPa, 电流频率f = 1.76 MHz, 工作气体为空气. 对于高功率ICP风洞, 由于存在较大的热能损失(比如, 冷却水引起的热损失, 电路上的电阻抗引起的能量损失以及辐射损失), 沉积到等离子体流动中的净输入功率通常与从供电系统读取的面板输入功率有所不同^[36], Fujita等^[28,37]对此进行了实验与仿真研究, 研究结果表明: 对于110 kWICP风洞其热效率通常为1/3. 因此, 本研究ICP风洞的热效率η也设定为1/3. 被等离子体实际吸收的沉积功率等于热效率η与面板输入功率的乘积. 仿真计算时, 为了确保能量守恒, 每一步流场迭代均使所有流体微元的焦耳加热率体积积分之和等于等离子体的沉积功率(即$\oint {{S_{{\rm{Joule}}}}} {\rm{d}}V = {P_{{\text{沉积}}}}$).
2

-->

2.1.流场方程组

为了简化计算, 本研究假设等离子体是电中性和二维轴对称的; 流体微元处于热化学非平衡态状态, 从微观上气体温度可分为平动温度${T_{{\rm{tr}}}}$、转动温度${T_{{\rm{rot}}}}$、振动温度${T_{{\rm{vib}}}}$和电子温度${T_{\rm{e}}}$. 仿真计算时, 考虑洛伦兹力、湍流以及等离子体产生的感应磁场. 最终所建立的非平衡磁流体力学方程组系统由总质量、动量、总能量守恒方程, 11组分粒子质量守恒方程, 分子振动、转动和电子能量守恒方程, 电磁场方程、湍流方程和32种空气化学反应动力学模型组成. 该方程组系统的矢量形式如下:

$\frac{{\partial {{Q}}}}{{\partial t}} + \frac{{\partial {{F}}}}{{\partial {x_j}}} = \frac{{\partial {{{F}}_{{\upsilon }}}}}{{\partial {x_j}}} + {{W}}.$

(1)式中, 守恒向量Q、通量矢量F、${{{F}}_{{\upsilon }}}$和源项矢量${{W}}$的表达式分别如下:

$\begin{split}&{{Q}} = \left[ \begin{gathered} \rho \\ \rho {u_i} \\ E \\ {\rho _s} \\ {E_{\rm rot}} \\ {E_{{\rm{vib}}}} \\ {E_{\rm{e}}} \\ \end{gathered} \right],\; {{F}} = \left[ \begin{gathered} \rho {u_j} \\ \rho {u_i}{u_j} + {\delta _{ij}}p \\ \left( {E + p} \right){u_j} \\ {\rho _s}{u_j} \\ {E_{{\rm{rot}}}}{u_j} \\ {E_{{\rm{vib}}}}{u_j} \\ {E_{\rm{e}}}{u_j} \\ \end{gathered} \right],\\&{{{F}}_{{\upsilon }}} = \left[ {\begin{array}{*{20}{c}} 0 \\ {{\tau _{ij}}} \\ {{q_j} + \rho \sum\limits_{s = 1}^{ns} {{h_s}{D_s}\dfrac{{\partial {X_s}}}{{\partial {x_j}}} + {u_i}{\tau _{ij}}} } \\ {\rho {D_s}\dfrac{{\partial {X_s}}}{{\partial {x_j}}}} \\ {{q_{{\rm{rot}},j}} + \rho \sum\limits_{s = M} {{e_{{\rm{rot}},s}}} {D_s}\dfrac{{\partial {X_s}}}{{\partial {x_j}}}} \\ {{q_{{\rm{vib}},j}} + \rho \sum\limits_{s = M} {{e_{{\rm{vib}},s}}} {D_s}\dfrac{{\partial {X_s}}}{{\partial {x_j}}}} \\ {{q_{{\rm{e}},j}} + \rho {h_{\rm{e}}}{D_{\rm{e}}}\dfrac{{\partial {X_{\rm{e}}}}}{{\partial {x_j}}} - {p_{\rm{e}}}{u_j}} \end{array}} \right], \\&{{W}} = \left[ \begin{gathered} 0 \\ {F_{{\rm{L}}i}} + {G_i} \\ {S_{{\rm{Joule}}}} \\ {{\dot \omega }_s} \\ {S_{{\rm{int,}}\;{\rm{rot}}}} \\ {S_{{\rm{int,}}\;{\rm{vib}}}} \\ {S_{{\rm{Joule}}}} + {S_{\operatorname{int},\;{\rm{e}}}}\end{gathered} \right].\end{split}$

(2)式中, ${\delta _{i,j}},\;{F_{{\rm{L}}i}},\;{G_i},\;{S_{{\rm{Joule}}}},\;{S_{{\mathop{\rm int}} }},\;\Delta h_s^0$和${\varTheta _{{\rm{vib,}}s}}$分别表示Kronecker函数、洛伦兹力、重力、焦耳加热效率、内能交换率、组分生成焓和振动特征温度; 下标s代表化学物种, M代表分子个数; ${\dot \omega _s}$为物质s的质量生成率, 并且是粒子质量守恒方程的源项(W的第4行); X_s 和D_s分别表示s的摩尔分数和扩散系数; ${u_i}$ (或者${u_j}$) 和 ${x_j}$代表速度和直角坐标系张量. (2)式中的应力张量${\tau _{ij}}$和热通量矢量${q_j}$分别由下式计算:

${\tau _{ij}} = \mu \left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}} - \frac{2}{3}\frac{{\partial {u_k}}}{{\partial {x_k}}}{\delta _{ij}}} \right),$

$ \begin{split}{q_j} &= {q_{{\rm{tr}},j}} + {q_{{\rm{rot}},j}} + {q_{{\rm{vib}},j}} + {q_{{\rm{e}},j}} \\&={\lambda _{{\rm{tr}}}}\dfrac{{\partial {T_{{\rm{tr}}}}}}{{\partial {x_j}}} + {\lambda _{{\rm{rot}}}}\dfrac{{\partial {T_{{\rm{rot}}}}}}{{\partial {x_j}}} + {\lambda _{{\rm{vib}}}}\dfrac{{\partial {T_{{\rm{vib}}}}}}{{\partial {x_j}}} + {\lambda _{\rm{e}}}\dfrac{{\partial {T_{\rm{e}}}}}{{\partial {x_j}}}.\end{split}$

高温气体的压强p由下式计算:

$p = \sum\limits_{s \ne {\rm{e}}}^{ns - 1} {{\rho _s}} {R_s}{T_{{\rm{tr}}}} + {\rho _{\rm{e}}}{R_{\rm{e}}}{T_{\rm{e}}} = \rho \hat R{T_{{\rm{tr}}}} + {p_{\rm{e}}}.$

气体的总内能E定义为

$ \begin{split}E =\, & {E_{{\rm{tr}}}} + {E_{{\rm{rot}}}} + {E_{{\rm{vib}}}} + {E_{\rm{e}}} \\&+\sum\limits_{s = 1}^{ns} {{\rho _s}} \Delta h_s^0 + \dfrac{1}{2}\rho {u_i}{u_i} + \dfrac{3}{2}\rho k.\end{split}$

气体粒子的平动能(${E_{{\rm{tr}}}}$)、转动能(${E_{{\rm{rot}}}}$)、振动能(${E_{{\rm{vib}}}}$)和电子内能(${E_{\rm{e}}}$)分别由下式计算:

${E_{{\rm{tr}}}} = \sum\limits_{s \ne {\rm{e}}} {\frac{3}{2}} {\rho _s}{R_s}{T_{{\rm{tr}}}},$

${E_{{\rm{rot}}}} = \sum\limits_{s = M} {{\rho _s}{R_s}{T_{{\rm{rot}}}}},$

${E_{{\rm{vib}}}} = \sum\limits_{s = M} {\dfrac{{{\rho _s}{R_s}{\varTheta _{{\rm{vib}},s}}}}{{{\rm{exp}}\left( {{{{\varTheta _{{\rm{vib}},s}}} / {{T_{{\rm{vib}},s}}}}} \right) - 1}}},$

${E_{\rm{e}}} = \sum\limits_{s \ne {\rm{e}}} {\frac{3}{2}} {\rho _{\rm{e}}}{R_{\rm{e}}}{T_{\rm{e}}}.$

(3)—(10)式中的气体输运系数(如气体黏性μ、平动导热系数${\lambda _{{\rm{tr}}}}$、振动导热系数${\lambda _{{\rm{vib}}}}$和转动导热系数${\lambda _{{\rm{rot}}}}$)均由Yos公式^[38,39]进行计算; 扩散系数${D_s}$由Curtiss和Hirschfelder^[40]给出的公式计算; 气体电导率σ和电子导热系数${\lambda _{\rm{e}}}$由三阶Sonine多项式和时新的电子-重粒子碰撞截面数据进行计算^[32,41,42]. 三阶精度电导率σ和电子导热系数${\lambda _{\rm{e}}}$的计算公式如下:

$\sigma \!=\! {\sigma ^{{\rm{3rd}}}} \!=\! \frac{{3{e^2}n_{\rm{e}}^2}}{{2k{T_{\rm{e}}}}}\sqrt {\frac{{2{\text{π}} k{T_{\rm{e}}}}}{{{m_{\rm{e}}}}}} {{\left| \begin{gathered} {q^{11}}\;{q^{12}} \\ {q^{21}}\;{q^{22}} \\ \end{gathered} \right|} \Bigg/ {\left| \begin{gathered} {q^{00}}\;{q^{01}}\;{q^{02}} \\ {q^{10}}\;{q^{11}}\;{q^{12}} \\ {q^{20}}\;{q^{21}}\;{q^{22}} \\ \end{gathered} \right|}},$

${\lambda _{\rm{e}}} = \lambda _{\rm{e}}^{{\rm{3rd}}} = \frac{{75n_{\rm{e}}^2k}}{{8k{T_{\rm{e}}}}}\sqrt {\frac{{2{\text{π}}k{T_{\rm{e}}}}}{{{m_{\rm{e}}}}}} \frac{{{q^{22}}}}{{{q^{11}}{q^{22}} - {{\left( {{q^{12}}} \right)}^2}}},$

式中${q^{mn}}$(m, n = 0—2) 为组分粒子数密度的函数, 它可由电子与其他重粒子的碰撞截面$\bar Q_{{\rm{ei}}}^{(l,s)}$(l = 1—2, s = 1—5)计算得到^[41]. 对于电子与带电粒子的碰撞截面$\bar Q_{{\rm{ei}}}^{(l,s)}$, 其可由库仑屏蔽电位计算^[42]; 电子和中性粒子的碰撞截面根据Laricchiuta等^[43]给出的方法计算.
2 -->

2.2.电磁场方程组

ICP的电磁场分布通常可通过求解以下磁矢势方程得到:

${\nabla ^2}{{A}} - \operatorname{i} {\mu _0}\omega \sigma {{A}} = - {\mu _0}{{{J}}_{\rm{c}}}.$

其中A为磁矢势, ω, ${J_{\rm{c}}}$, ${\mu _0}$和i分别表示线圈电流的角频率、电流密度、真空介电常数和复数单位(${\rm{i}} = \sqrt { - 1} $), 且$ \omega = 2{\text{π}}f$, f为线圈电流驱动频率. 此外, 由于磁矢势A与磁场强度H和电场强度E有如下关系^[20]:

${\mu _0}{{H}} = \nabla \times {{A}},$

${{E}} = \; - {\rm{i}}\omega {{A}},$

因此可通过求解磁矢势来得到电场和磁场分布, 进而得到等离子体的焦耳加热率和洛伦兹力.
对于圆柱形ICP炬, 线圈电流可以被认为是由若干平行环状电流微元组成, 因此磁矢势可以被认为只有切向分量^[20], 即${{A}} = ({A_r},\;{A_\theta },\;{A_z}) = (0,\;{A_\theta },\;0)$. 考虑到线圈电流产生的电磁场与等离子体产生的电磁场之间存在相位差, 故切向矢量势${A_\theta }$可表示为复数, 即${A_\theta } = {A_{\rm{R}}} + {\rm{i}}{A_{\rm{I}}}$. 最终, (13)式可写成如下形式:

${\nabla ^2}{A_{\rm{R}}} + \omega {\mu _0}\sigma {A_{\rm{I}}} = - {\mu _0}{J_{\rm{c}}},$

${\nabla ^2}{A_{\rm{I}}} - \omega {\mu _0}\sigma {A_{\rm{R}}} = 0,$

式中, 电流密度${J_{\rm{c}}}$可由${J_{\rm{c}}} = {I / {\left( {{\text{π}}R_{\rm{c}}^2} \right)}}$计算, 其中R_c为线圈半径, I为线圈电流. 待求解出磁矢势方程后, ICP的焦耳加热率S_Joule可由下式计算得到:

${S_{{\rm{Joule}}}} = \frac{1}{2}\sigma {\omega ^2}\left( {A_{\rm{R}}^{\rm{2}} + A_{\rm{I}}^{\rm{2}}} \right).$

轴向和径向洛伦兹力F_Lx和F_Ly可分别由下式进行计算:

${F_{{\rm{L}}x}} = - \frac{{\sigma \omega }}{2}\left( {{A_{\rm{R}}}\frac{{\partial {A_{\rm{I}}}}}{{\partial x}} - {A_{\rm{I}}}\frac{{\partial {A_{\rm{R}}}}}{{\partial x}}} \right),$

${F_{{\rm{L}}y}} = - \frac{{\sigma \omega }}{2}\left( {{A_{\rm{R}}}\frac{{\partial {A_{\rm{I}}}}}{{\partial y}} - {A_{\rm{I}}}\frac{{\partial {A_{\rm{R}}}}}{{\partial y}}} \right),$

式中, σ为等离子体的电导率, 焦耳加热率和洛伦兹力将被作为能量守恒和动量守恒方程的源项参与到流场方程组的迭代求解中. ICP风洞的电磁场与流场将通过焦耳加热率和洛伦兹力进行耦合联接.
2 -->

2.3.空气化学反应模型

为了准确模拟空气在高温条件下发生的离解、电离、电量交换、分子复合等化学反应, 本文将文献[33, 34]给出的由11种组分(N₂, O₂, NO, N₂⁺, O₂⁺, NO⁺, N, O, N⁺, O⁺, e^–)和32种化学反应组成的动力学模型应用到本研究的数值计算中. 该32种化学反应如表1所列, 化学反应r的正向反应速率${k_{{\rm{f}},r}}$由Arrhenius公式进行计算^[44]:

	r	反应物		生成物	Tf	Tb	Cr	n	θr	文献
离解/复合反应 (S1 = N2, O2, NO; S2 = N, O; S3 = N2, O2; S4 = NO, N, O)	1—3	N2 + S1	$ \rightleftharpoons $	N + N + S1	$\sqrt {{T_{{\rm{tr}}}}{T_{{\rm{vib}}}}} $	${T_{{\rm{tr}}}}$	7.0 × 1021	–1.60	113200	[35]
	4—5	N2 + S2	$ \rightleftharpoons $	N + N + S2	$\sqrt {{T_{{\rm{tr}}}}{T_{{\rm{vib}}}}} $	${T_{{\rm{tr}}}}$	3.0 × 1022	–1.60	113200	[35]
	6—8	O2 + S1	$ \rightleftharpoons $	O + O + S1	$\sqrt {{T_{{\rm{tr}}}}{T_{{\rm{vib}}}}} $	${T_{{\rm{tr}}}}$	2.0 × 1021	–1.50	59500	[35]
	9—10	O2 + S2	$ \rightleftharpoons $	O + O + S2	$\sqrt {{T_{{\rm{tr}}}}{T_{{\rm{vib}}}}} $	${T_{{\rm{tr}}}}$	1.0 × 1022	–1.50	59500	[35]
	11—12	NO + S3	$ \rightleftharpoons $	N + O + S3	$\sqrt {{T_{{\rm{tr}}}}{T_{{\rm{vib}}}}} $	${T_{{\rm{tr}}}}$	5.0 × 1015	0.00	75500	[35]
	13—15	NO + S4	$ \rightleftharpoons $	N + O + S4	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	1.1 × 1017	0.00	75500	[35]
泽尔多维奇反应	16	N2 + O	$ \rightleftharpoons $	NO + N	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	6.4 × 1017	–1.00	38400	[35]
	17	NO + O	$ \rightleftharpoons $	N + O2	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	8.4 × 1012	0.00	19450	[35]
电量交换反应	18	N2 + N+	$ \rightleftharpoons $	N2+ + N	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	1.0 × 1012	0.50	12200	[35]
	19	O2+ + O	$ \rightleftharpoons $	O+ + O2	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	4.0 × 1012	–0.09	18000	[35]
	20	NO+ + O	$ \rightleftharpoons $	NO + O+	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	3.63 × 1015	–0.60	13000	[33]
	21	O+ + N2	$ \rightleftharpoons $	N2+ + O	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	9.1 × 1011	0.36	22800	[35]
	22	NO+ + O2	$ \rightleftharpoons $	O2+ + NO	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	2.4 × 1013	0.41	32600	[35]
	23	NO+ + N	$ \rightleftharpoons $	NO + N+	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	1.0 × 1019	–0.93	61000	[33]
	24	NO+ + O	$ \rightleftharpoons $	N+ + O2	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	1.0 × 1012	0.50	77200	[35]
副电离反应	25	N + N	$ \rightleftharpoons $	N2+ + e–	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	4.4 × 107	1.50	67500	[35]
	26	O + O	$ \rightleftharpoons $	O2+ + e–	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	7.1 × 102	2.70	80600	[35]
	27	N + O	$ \rightleftharpoons $	NO+ + e–	${T_{{\rm{tr}}}}$	${T_{{\rm{tr}}}}$	8.8 × 108	1.00	31900	[35]
	28	O2 + N2	$ \rightleftharpoons $	NO + NO+ + e–	$\sqrt {{T_{\rm{e}}}{T_{{\rm{vib}}}}} $	${T_{\rm{e}}}$	1.38 × 1020	–1.84	141000	[33]
	29	N2 + NO	$ \rightleftharpoons $	N2 + NO+ + e–	$\sqrt {{T_{\rm{e}}}{T_{{\rm{vib}}}}} $	${T_{\rm{e}}}$	2.20 × 1015	–0.35	108000	[33]
	30	O2 + NO	$ \rightleftharpoons $	O2 + NO+ + e–	$\sqrt {{T_{\rm{e}}}{T_{{\rm{vib}}}}} $	${T_{\rm{e}}}$	8.80 × 1016	–0.35	108000	[33]
电子碰撞电离反应	31	N + e–	$ \rightleftharpoons $	N+ + e– + e–	${T_{\rm{e}}}$	${T_{\rm{e}}}$	2.5 × 1034	–3.82	168600	[35]
	32	O + e–	$ \rightleftharpoons $	O+ + e– + e–	${T_{\rm{e}}}$	${T_{\rm{e}}}$	3.9 × 1033	–3.78	158500	[35]

表1空气化学反应模型
Table1.Chemical reaction model of air.

${k_{{\rm{f}},r}}\left( {{T_{\rm{f}}}} \right) = {C_r}T_{\rm{f}}^n\exp \left( { - {\theta _r}/{T_{\rm{f}}}} \right),$

式中的化学反应系数${C_r}$, n 和 ${\theta _r}$取自文献[33, 35].
逆向化学反应速率${k_{{\rm{b}},r}}$由平衡常数K^eq进行计算:

${k_{{\rm{b}},r}}\left( {{T_{\rm{b}}}} \right) = {k_f}\left( {{T_{\rm{b}}}} \right)/{K^{{\rm{eq}}}}\left( {{T_{\rm{b}}}} \right).$

平衡常数由曲线拟合公式进行计算^[35]. 最终, 任一空气组分s 因发生化学反应而产生的质量变化率净质量生成率${\dot \omega _s}$可由下式计算:

$\begin{split}{\dot \omega _s} = & {M_s}\sum\limits_{r = 1}^{nr} {\left( {\nu _s^{{\rm{b}},r} - \nu _s^{{\rm{f}},r}} \right)} \\&\times\left[ {{k_{{\rm{f}},r}}{{\prod\limits_{j = 1}^{ns} {\left( {\dfrac{{{\rho _j}}}{{{M_j}}}} \right)} }^{{\nu _{j,r}}}} - {k_{{\rm{b}},r}}{{\prod\limits_{j = 1}^{ns} {\left( {\dfrac{{{\rho _j}}}{{{M_j}}}} \right)} }^{{{\nu'}\!_{j,r}}}}} \right],\end{split}$

式中, ν_s^b,r 和 ν_s^f,r是组分s在进行第r个化学反应前后的化学计量系数, M_s代表物质s的摩尔质量, 下标s和j代表物质组分编号.
2 -->

2.4.内部能量交换模型

由于电子、分子和原子之间存在弹性或非弹性碰撞, 碰撞时会产生平动能、转动能、振动能或电子能量的交换, 为了准确模拟这些内能交换过程, 根据相应的数学模型计算粒子内能交换率, 并将计算得到的转动能交换率S_int,rot、振动能交换率S_int,vib和电子能量交换率S_int,e作为能量源项添加到(2)式中的转动、振动和电子能量守恒方程中. 转动、振动和电子能量交换率S_int,rot, S_int,vib, S_int,e的具体计算方法如下:

${S_{{\rm{int,rot}}}} = {Q_{{\rm{T}}\text{-} {\rm{R}}}} - {Q_{{\rm{R}} \text{-}{\rm{V}}}} - {Q_{{\rm{R}} \text{-}{\rm{e}}}} + Q_{\rm{D}}^{{\rm{rot}}},$

${S_{{\rm{int,vib}}}} = {Q_{{\rm{T}} \text{-} {\rm{V}}}} + {Q_{{\rm{R}} \text{-} {\rm{V}}}} + {Q_{{\rm{e}} \text{-} {\rm{V}}}} + Q_{\rm{D}}^{{\rm{vib}}},$

${S_{{\rm{int,e}}}} = {Q_{{\rm{T}} \text{-}{\rm{e}}}} + {Q_{{\rm{R}} \text{-} {\rm{e}}}} - {Q_{{\rm{e}} \text{-}{\rm{V}}}} + Q_{\rm{D}}^{\rm{e}} + Q_{\rm{I}}^{\rm{e}},$

式中, 平动能-转动能交换率Q_T-R由Park^[45]给出的方法计算; 转动能-振动能量交换率Q_R-V由Millikan和White ^[46]给出的方法计算; 转动能-电子能交换率 Q_R-e由Lazdinis和Petrie^[47]给出的方法计算; 平动能-振动能量交换率Q_T-V 由Millikan和White ^[46] 和Park ^[48]给出的公式计算; 电子能-振动能交换率 Q_e-V 由 Landau-Teller 方程计算^[49], 其中氮分子和电子的振动能-电子能驰豫时间根据Kim等^[50]及Bourdon和Vervisch ^[51]给出的方法进行计算; 平动能-电子能交换率Q_T-e 由Appleton和Bray^[52]给出的公式计算; 空气分子因离解和电离反应产生的能量损失$Q_{\rm{D}}^{{\rm{rot}}},\;Q_{\rm{D}}^{{\rm{vib}}},\;Q_{\rm{D}}^{\rm{e}},\;Q_{\rm{I}}^{\rm{e}}$由 Gnoffo等^[53]给出的方法分别进行计算.
2 -->

2.5.湍流运输方程

本文采用低雷诺数湍流模型—Abe-Kondoh-Nagano k-ε模型^[54]来考虑湍流对ICP流动的影响, 包含湍流动能k和耗散率ε的湍流输运方程如下:

$\begin{split}& \dfrac{{\partial \rho k}}{{\partial t}} + \dfrac{\partial }{{\partial {x_j}}}\left( {\rho k{u_j}} \right) \\=\; &\dfrac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \dfrac{{{\mu _{\rm{t}}}}}{{{\sigma _k}}}} \right)\dfrac{{\partial k}}{{\partial {x_j}}}} \right] + \tau _{ij}^{\rm{t}}\dfrac{{\partial {u_j}}}{{\partial {x_j}}} - \rho \varepsilon,\end{split}$

$\begin{split}\dfrac{{\partial \rho \varepsilon }}{{\partial t}} + \dfrac{\partial }{{\partial {x_j}}}\left( {\rho \varepsilon {u_j}} \right) = \;& \dfrac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \dfrac{{{\mu _{\rm{t}}}}}{{{\sigma _\varepsilon }}}} \right)\dfrac{{\partial \varepsilon }}{{\partial {x_j}}}} \right] \\& +{C_{\varepsilon 1}}{f_1}\dfrac{\varepsilon }{k}\tau _{ij}^{\rm{t}}\dfrac{{\partial {u_j}}}{{\partial {x_j}}} - {C_{\varepsilon 2}}{f_2}\rho \dfrac{{{\varepsilon ^2}}}{k}.\end{split}$

湍流黏度由下式计算:

${\mu _{\rm{t}}} = {C_\mu }{f_\mu }\rho \frac{{{k^2}}}{\varepsilon }.$

上述方程中的模型常数如下:

$\begin{split}& {C_\mu } = 0.09,\;\;\;{\sigma _k} = 1.4,\;\;\;{\sigma _\varepsilon } = 1.4,\\& {C_{\varepsilon 1}} = 1.5,\;\;\;{C_{\varepsilon 2}} = 1.9.\end{split}$

模型函数表示为

${f_\mu }\!=\! {\left[ {1 \!- \!\exp \left( { - \frac{{{y^*}}}{{14}}} \right)} \right]^2}\left\{\! {1 \!+\! \frac{5}{{R_{\rm{t}}^{3/4}}}\exp \left[ \!{ - \!{{\left( {\frac{{{R_{\rm{t}}}}}{{200}}} \right)}^2}} \right]} \!\right\},$

${f_1} = 1,$

${f_2} \!=\! {\left[ {1 \!-\! \exp \left( { - \frac{{{y^*}}}{{3.1}}} \right)} \right]^2}\left\{\! {1 \!-\! 0.3\exp \left[ { - {{\left( {\frac{{{R_{\rm{t}}}}}{{200}}} \right)}^2}} \right]} \!\right\},$

式中, ${y^*} = {(\nu \varepsilon )^{1/4}}{y_{{\rm{wd}}}}/\nu $, ${R_{\rm{t}}} = {k^2}/\left( {\nu \varepsilon } \right)$, ν表示分子运动黏度, y_wd为与内壁面的距离.
2 -->

2.6.边界条件和数值求解方法

3 -->

2.6.1.边界条件

对于流场方程组: 1)在等离子体炬入口处, 工作气体由距壁面约2.4 mm的开口注入, 气体质量流率和初始温度T₀ = 300 K作为已知输入参数. 2)在ICP炬出口处, 工作气压设为定值p = 10 kPa, 其他参数由相邻的内部网格点线性插值得到. 3)在等离子体炬壁面处, 认为壁面上无滑移、无催化效应、无压力梯度, 壁面温度由导热方程${{\partial {T_{\rm{w}}}} / {\partial n}} = {q_{j{\rm{max}}}}$进行计算, 其中q_jmax为内壁面处的热流量密度. 4)在中心轴上, 采用轴对称边界条件, 即Q_i,j = Q_i,j+1, Q代表守恒变量. 对于磁矢势方程组, 电磁场计算区域的外边界(图2(a)中的x = –120 mm, x = 360 mm和y = 206 mm)设置的离感应线圈足够远以使磁矢势A_R和A_I在这些外边界上数值为零, 由于100 kW级高功率ICP风洞的电磁场强度较大, 其影响范围也更广, 故与10 kW ICP风洞的磁场外边界(电磁场趋近于0的边界)相比, 本文的远场电磁外边界向外侧进行了推移. 中心轴线上x = 0处我们采用轴对称边界条件: A_i,0 = A_i,1.
3 -->

2.6.2.数值求解方法

对于流场和湍流方程组, 本文采用有限体积法对其进行离散化, 通过时新的低耗散迎风差分法计算其对流项, 并通过MUSCL (monotonic upstream-centered scheme for conversation laws)格式将其计算精度提高到二阶精度; 黏性项采用二阶精度中心差分法进行计算. 对于非定常物理流动, 为获得其定常流场, 除了要在空间上进行推进求解外, 也需要在时间方向上进行推进求解, 故本文采用点隐式LU-SGS (lower and upper symmetric Gauss-Seidel)时间推进算法对流场方程组进行求解, 采用广义最小残差法(generalized minimal residualal method, GMRES)对湍流输运方程组进行求解.
对于电磁场方程组, 本文采用有限差分法对其进行离散化, 采用二阶中心差分法计算其数值通量, 采用亚松弛迭代法快速、稳定地求解离散后的方程, 松弛因子设定为0.2, 当磁矢势第n步和第n+1步数值解的相对误差小于10^–5时, 认为计算收敛. 对于化学反应方程组, 本文采用有限体积法对方程组进行离散化, 利用托马斯追赶法进行迭代求解.

下面将通过对100 kW级空气ICP风洞在典型工况下(P = 90 kW, m = 1.8 g/s, p = 10 kPa, f = 1.76 MHz)的速度、电子温度、压强、洛伦兹力、焦耳加热率、电场强度、电子数密度、粒子摩尔分数、平动温度的分布规律进行分析与总结, 明确ICP的流场与电磁场分布特性, 并揭示其流场与电磁场的相互作用机理.
ICP炬中气流的速度矢量、流线和电子温度分布云图如图3所示, 由于感应线圈区存在很强的外加电磁能, 从入口流入的常温新鲜空气在此区域经点火后迅速发生放电反应并吸收大量的焦耳热, 因此气体速度和温度在此区域都迅速升高. 由于空气中的氮分子与氧分子在此区域被迅速电离和离解并释放大量的热, 因此由电子、阳离子和原子组成的ICP温度在感应线圈区很高, 最大电子温度约10500 K. 此外, 从图3中的流线分布还可以发现, 在进气口附近形成了较强的气流回旋现象. 该回流的形成与感应线圈区的负压强梯度和电磁加热现象有很大关系.
图 3 等离子体炬内气体流线和速度矢量(上半部分)以及电子温度云图(下半部分)分布
Figure3. Distributions of streamlines and velocity vector (upper), and electron temperature (lower) in the torch.

图4为等离子体流动的流线与压力分布云图. 可见, 最大气压处在第二和第三圈线圈之间, 即在第二个线圈上游存在负压力梯度, 这是线圈上游靠近入口处产生涡流的原因之一. 此外, 在动量方程中, 速度和压力的分布都受到洛伦兹力的影响, 而且在能量守恒方程中焦耳加热率对流动速度也有影响. 因此, 线圈下方产生的大面积涡流是由感应线圈区域的负压力梯度、洛伦兹力和焦耳加热现象共同作用产生的.
图 4 等离子体流线(上半部分)和压力云图(下半部分)分布
Figure4. Distributions of streamlines (upper) and pressure contour (lower).

图5给出了轴向和径向洛伦兹力分布. 洛伦兹力作为能量守恒方程的惯性力源项参与到整个流场的迭代求解中. 由图5可见: 轴向洛伦兹力在感应线圈的第一与第二圈之间为正值, 即其方向向右. 然而在感应线圈第二圈与第三圈之间变为负值, 即轴向洛伦兹力的方向反向. 而对于径向洛伦兹力(图5下半部分), 其值始终为负值, 即其方向始终由壁面指向中心轴线, 这表明因趋肤效应在炬壁附近生成的自由电子将始终受到指向等离子体炬中轴线的惯性力, 且径向洛伦兹力的最大值比轴向洛伦兹的最大值高3.95倍, 说明等离子体的径向动量传递比轴向动量传递强烈的多.
图 5 轴向洛伦兹力(上半部分)和径向洛伦兹力(下半部分)的分布
Figure5. Distributions of axial Lorentz force (upper) and radial Lorentz force (lower).

图6给出了径向洛伦兹力(上半部分)和焦耳加热率(下半部分)分布. 可以看出, 焦耳加热率的分布形状和位置与径向洛伦兹力的非常相似. 这表明, 对于空气ICP流动, 进行焦耳加热现象的位置主要由径向洛伦兹力控制. 此外, 径向洛伦兹力还对等离子体最大温度和速度的位置有一定影响. 另一方面, 由于感应线圈上通有高频交变电流, 线圈区域会发生气体放电和焦耳加热现象, 由于焦耳加热率由气体电导率和电场E共同控制, 且等离子体炬壁上始终通有冷却水, 壁面温度被限制在1000 K以下, 导致壁面附近的流动温度较低, 即导电率很小, 因此, 尽管电场强度E在等离子体炬内壁面上很大, 但最大焦耳加热率并不出现在内壁面上, 而是出现在第二个感应线圈下方距离内壁面3.5 mm处.
图 6 径向洛伦兹力(上半部分)和焦耳加热率(下半部分)的分布
Figure6. Distributions of Joule heating rate(lower) and radial Lorentz force (upper)

图7给出了等离子体炬内电场强度虚部E_I (上半部分)和实部E_R (下半部分)的分布云图. 可知, 电场强度虚部的最大值比电场强度实部大2.9倍, 所以由磁矢势计算得到的电场强度虚部E_I是总电场强度的主要部分. 电场强度实部E_R在等离子体炬中始终为负值; 而虚部E_I在靠近壁面处为负值, 在第二个线圈下方靠近轴线位置处变为正值, 负的电场虚部主要是由外加的高频电流产生, 而正电场则是由等离子体内部电子和阳离子产生的感应电场.
图 7 电场强度分布(虚部E_I (上半部分)实部E_R(下半部分))
Figure7. Distribution of electric-field intensity (imaginary part E_I (upper) and real part E_R (lower)).

图8给出了等离子体炬中电场强度E_I(上半部分)和电子数密度n_e(下半部分)分布. 在强电场的作用下, 电子在径向洛伦兹力作用下向中轴线附近移动, 在第二个感应线圈下方距离等离子体炬内壁面5.5 mm处电子数密度达到最大值. 大量的自由电子聚集在感应线圈区, 电子数密度n_e = 1.0 × 10²⁰ m^–3所包络的区域与正的电场强度E_I所处位置近似相同, 这表明正的电场强度E_I由这些自由电子产生, 它是这些自由电子在等离子体炬中自由流动时形成的感应电流产生的感应电场.
图 8 电场强度E_i(上)和电子数密度n_e(下)的分布
Figure8. Distribution of electric field intensity E_I (upper) and electron number density n_e (lower).

图9为轴向位置x = 68 mm处11组分空气的摩尔分数分布. 如图9所示, 氧分子被完全解离和电离成氧原子和离子, 91%的氮分子在轴线附近被离解成N原子, 氮原子N和氧原子O在径向位置y ≤ 30 mm之前为最主要的化学组分. 因电离反应生成的自由电子和离子在本算例中仍然较少, 电子e^–和氧离子O⁺的摩尔分数约为10^–3, 二者电荷数趋于一致, 说明等离子体呈准电中性. 沿等离子体炬半径方向, 随着平动温度的降低, 氮分子N₂和NO分子的摩尔分数逐渐增大; 由于等离子体炬壁附近(30 mm < y < 37.5 mm)气流温度较低, 氮分子几乎不发生化学反应, 因此在靠近等离子体炬壁处N₂, O和O₂是最主要的空气粒子.
图 9 感应线圈中心(x = 68 mm)空气粒子径向摩尔分数分布
Figure9. Mole fraction of air species along the radial direction at the coil center x = 68 mm.

图10为等离子体平动温度和电子温度分布云图, 最大电子温度(9921 K)和最大平动温度(8507 K)均出现在第二个感应线圈下方靠近等离子体炬壁处. 在60 ≤ x ≤ 85 mm, 20 ≤ y ≤ 30 mm区域, 电子温度(8500 ≤ T ≤ 9500 K)比平动温度高约1000 K, 即在该区域等离子体处于热力学非平衡态. 然而, 在靠近中轴线附近, 电子温度与平动温度基本相等, 气体温度$T \approx$ 6500 K, 说明该区域空气重粒子(如N, O)与电子之间的能量交换基本达到了热力学平衡, 此时等离子体流动趋于局域热力学平衡态.
图 10 等离子体炬内平动温度(上半部分)和电子温度(下半部分)分布云图
Figure10. Distributions of translational (upper) and electronic temperatures (lower) in the torch.

本文基于多物理场耦合流动仿真研究了空气ICP流场与电磁场的分布特性及其相互作用机理, 通过数值模拟获得了100 kW级ICP炬内等离子体流动的电子温度、平动温度、粒子摩尔分数、速度、压强、洛伦兹力、电场强度、焦耳加热率等参数的分布规律, 研究发现:
1)在进气口与第二圈感应线圈之间出现了较强的涡流, 该涡流与感应线圈区的负压强梯度和电磁加热现象有很大关系. 径向洛伦兹力方向始终为负说明因趋肤效应在壁面附近生成的自由电子始终受到指向ICP炬中轴线的洛伦兹力作用; 且径向洛伦兹力的最大值比轴向洛伦兹的高3.95倍, 这表明动量传递以径向为主; 空气粒子的焦耳热效应也受径向洛伦兹力的影响.
2)电场强度虚部E_I的最大值比电场强度实部E_R的大2.9倍, 负的电场虚部E_I主要是由外加高频电流产生, 而正的电场虚部则是由等离子体内部的自由电子产生, 在第二圈感应线圈下方距离等离子体炬内壁面5.5 mm处自由电子的数密度达到最大值.
3)空气在感应线圈下方发生剧烈的离解和电离反应, 氧分子被完全解离和电离成氧原子或离子; 91%的氮分子在中轴线附近被离解成N原子, N和O原子在径向位置y ≤ 30 mm之前为最主要的化学组分. 本次模拟得到的空气ICP最大电子温度(9921 K)和最大平动温度(8507 K)均出现在第二圈感应线圈下方靠近等离子体炬壁面处. 在60 mm ≤ x ≤ 85 mm且 20 mm ≤ y ≤ 30 mm区域, 电子温度比平动温度高大约1000 K, 流动处于热力学非平衡状态. 然而在靠近中轴线附近, 电子温度与平动温度基本相等(约为6500 K), 在该区域空气重粒子与电子之间的能量交换处于局域热平衡.
感谢国家超级计算广州中心提供的天河二号超级计算机计算服务支持. 感谢日本宇宙航空研究开发机构Kazuhiko Yamada等的建议与支持.