温稠密铝等离子体物态方程及其电离平衡研究*

王天浩 王坤† 张阅 姜林村

1) (河北工业大学, 省部共建电工装备可靠性与智能化国家重点实验室, 天津 300130)

2) (河北工业大学, 河北省电磁场与电器可靠性重点实验室, 天津 300130)

温稠密物质的物性参数在惯性约束聚变能源、Z箍缩等高能量密度物理领域的实验结果分析和物理过程数值模拟等方面有着重要的应用价值.本文应用部分电离等离子体模型, 在理想自由能的基础上考虑了库仑相互作用、排斥体积作用和极化作用等非理想特性, 开展了温稠密等离子体物态方程和电离平衡的研究.计算了温稠密铝等离子体的压强等物态方程数据和在密度为 1.0 × 10–4—3.0 g/cm3, 温度为 1.0 × 104—3.0 ×104 K范围内的粒子组分.计算结果显示, 铝等离子体的平均电离度在临界密度区域内随着密度的增加而突然增大.根据非理想Saha方程中有效电离能这一关键参数, 分析了铝等离子体平均电离度在临界密度区域内随密度迅速增大的现象.

1 引 言

温稠密物质是指介于固态和理想等离子体态之间的物质形态[1], 常见于丝阵Z箍缩[2]、X射线源[3]、惯性约束聚变能源[4]等高能量密度物理领域.物态方程和电离平衡模型描述一定热力学条件下物质的宏观性质和粒子组分[5,6], 是高能量密度物理领域相关实验结果分析和物理过程数值模拟的重要物性参数.

物质在固态和理想等离子体态都有较为精确的物态方程模型.基于晶格点阵运动理论的Grüneisen物态方程是比较完善的高压固体物态方程[7], 而理想等离子体满足经典Boltzmann统计,可以采用理想气体物态方程进行描述.温稠密物质处于部分简并、强耦合和非理想等离子体状态[8],这给准确描述温稠密等离子体特性带来了极大的困难和挑战.研究者通常将不同相区的物态方程模型进行衔接处理, 或采用半经验的方法来建立适用温度-密度范围宽的多项物态方程[9].Thomas-Fermi (TF)模型及相关的改进模型以其适用温度-密度范围宽、数学处理简洁等优点常用于描述温稠密等离子体的性质[10], 然而TF模型只能估算等离子体的平均电离度, 不能得到不同价态的离子分布.为了获得更为精确的物性参数, 从头算方法、化学模型等也应用于温稠密物质物态方程和输运参数研究[11].汤文辉等[12]根据离子耦合系数和电子简并度参数对等离子体进行了划分, 并对不同区域内等离子体物态方程理论进行了总结.张颖等[13]结合了微扰理论和线性响应理论, 在位形相互作用中考虑了原子间的微扰自由能和量子力学修正自由能, 库仑相互作用中考虑了电子间的交换作用,能够自洽地对电离能进行修正, 建立了稠密氦等离子体的物态方程, 计算了电离度、压强等物性参数,得到了粒子间的非理想相互作用引起的电离能降低是出现压致电离现象的主要原因的重要结论, 并将模型进一步改进, 推广至氩、疝等惰性元素等离子体物态方程的计算中[14−16].俄罗斯高温联合研究院的Apfelbaum[17,18]在提出的广义化学模型中考虑了电子间交换相互作用的影响, 计算了镍、钛、锌等离子体在温稠密区域的物态方程.罗斯托克大学的Kuhlbrodt等[19]利用考虑了库仑相互作用和极化作用的自由能模型计算了等离子体组分及输运参数.考虑非理想特性的温稠密等离子体特性研究显示平均电离度在临界密度区域内随着密度增加突然增大, 而解释这一现象需要针对温稠密等离子体物态方程和电离平衡开展更为深入的研究.

本文基于部分电离等离子体模型, 在理想自由能模型的基础上考虑了库仑相互作用、排斥体积作用和极化作用等非理想特性, 对温稠密铝等离子体物态方程和电离平衡开展了研究.获得了宽温度-密度范围内的物态方程和等离子体组分数据, 分析了非理想特性对等离子体组分内不同价态离子电离能的影响.通过自由能和电离能的变化对临界密度区域内温稠密铝等离子体平均电离度随密度增加突然上升的现象进行了分析.

2 温稠密等离子体物态方程模型

温稠密等离子体处于部分电离等离子体状态,由电子、原子和不同价态的离子组成.对于给定体积V、温度T的热力学系统, 其自由能为F, 本文采用自由能密度f = F/V来建立温稠密等离子体物态方程模型.等离子体的自由能密度f可以分为理想自由能部分fid和非理想自由能部分fint, 理想自由能密度可以表示为, 其中上标i,e分别表示重粒子和电子.重粒子的理想自由能密度可以表示成如下形式[16,18]

其中kB是玻尔兹曼常数; nk是粒子数密度; Qk是粒子的内部电子配分函数; lk= h/(2πmkkBT)1/2是德布罗意波长, h是普朗克常数, mk是粒子质量.电子的理想自由能密度可以表示为[16]

其中 ne是电子数密度;是理想电子化学势; In是Fermi-Dirac积分.

当等离子体处于温稠密区域时, 粒子间的相互作用不可忽略[20].在本文中, 非理想自由能密度fint主要考虑带电粒子之间的库仑相互作用fcoul、重粒子之间的排斥体积作用fhs以及中性原子与带电粒子之间的极化作用fpol的贡献, 如下式所示:

带电粒子之间的库仑相互作用可以分为四部分[21]

排斥体积作用对自由能密度的贡献可以表示为[19,23]

带电粒子和中性原子间的极化作用可以表示为[24,25]

其中B0j是第二维里系数, 可以由极化势UPP计算,极化势表达式为

其中e0是真空介电常数, aD是偶极化率, k是德拜长度, rc是截断半径.本文采用部分电离等离子体模型来描述温稠密等离子体, 因此极化作用项中只计算了电子-中性原子相互作用.通过上述物态方程模型可以获得温稠密等离子体压强、粒子化学势等热力学量.

3 电离平衡模型

电离平衡模型对准确计算物态方程和表征等离子体粒子组分有着重要作用.铝等离子体内电离过程如下式所示:

式中k = 0表示中性原子.若等离子体处于电离平衡状态, 不同电离态离子、中性原子和电子也处于相平衡状态, 由化学势表示的多粒子体系电离平衡对应的相平衡表达式为[26]

其中 Ek + 1是电离能, µk是重粒子化学势, µe是电子化学势.

温稠密等离子体物态方程中自由能密度划分为理想部分和非理想部分.与之相对应, 等离子体中粒子化学势也可划分为理想部分和非理想部分,即µ = µid+ µint, 将该式代入 (9)式可以得到

等离子体中粒子化学势的理想部分及非理想部分可由自由能密度求解, 如下式所示:

将(1)式和(11)式代入(10)式, 化简可得到非理想Saha方程

非理想Saha方程结合质量守恒定律和电荷守恒定律形成的封闭方程组, 通过与物态方程模型迭代求解可以得到等离子体的组分.质量守恒定律和电荷守恒定律如(13)式和(14)式所示:

本文应用上述物态方程和电离平衡模型针对温稠密铝等离子体物性参数开展了研究, 模型中相关参数(里德伯单位制)详见表1.

表1 物态方程及电离平衡模型参数列表Table 1.Parameters of equation of state and ionization equilibrium model.

4 计算结果与讨论

温稠密等离子体压强等热力学量可通过自由能求得 P = –( ∂ F/ ∂ V)|T.本文模型计算的铝等离子体压强在0.1 g/cm3时随温度的变化曲线见图1,由 SCAALP-FT, QEOS 以及 VASP 等模型[27−29]计算的结果也列于图中对比.本文得到的压强曲线与VASP, SCAALP-FT模型计算的结果较为接近, 而QEOS模型的结果明显高于其他模型, 并且差别随着温度升高而变大.

图1 在密度为 0.1 g/cm3 时, 不同模型计算的铝等离子体压强随温度的变化Fig.1.The pressure of aluminum plasma calculated by different models as a function of temperature at density of 0.1 g/cm3.

描述多粒子体系的物态方程与非理想Saha方程迭代求解可获得等离子体的粒子组分.本文计算了铝等离子体在密度 1.0 × 10–4g/cm3< r <3.0 g/cm3, 温度 1.0 × 104K < T < 3.0 × 104K范围内平均电离度的等值线分布, 结果如图2所示.等离子体的电离机制在低密度、高温度区域由热电离主导, 而在高密度、低温度区域由压致电离主导.临界点是气-液相界线的终点, 已有大量研究工作对铝的临界点参数进行了估算[30−32], 一般认为金属铝的临界温度处于 8.0 × 103—1.0 × 104K范围内, 临界密度处于 0.4—1.0 g/cm3范围内.在本文中, 将上述密度区域称为临界密度区域.在临界密度区域内, 平均电离度呈现随密度增大而迅速上升趋势.在高密度区域, 当温度小于 1.2 × 104K时, 铝等离子体主要发生一级电离和二级电离, 体系主要由Al1+, Al2+和电子组成, 平均电离度的最大值为 2; 而当温度大于 1.2 × 104K 时, 二级电离随密度增大而减弱, 铝等离子体发生三级电离, 等离子体中带电粒子主要由Al3+和电子组成, 平均电离度的最大值为3.

图2 铝等离子体的平均电离度随密度与温度变化等值线分布Fig.2.Contour map of average ionization degree of aluminum plasma as a function of density and temperature.

等离子体中不同粒子的相对粒子分数在密度为0.1 g/cm3时随温度的变化如图3所示.从图3可以看出, 铝等离子体在温度为 1.0 × 104—2.0 ×106K的范围内电离形成的离子最高, 为六价铝离子.铝原子相对粒子分数随着温度的升高不断降低, 当温度达到 3.0 × 104K 时, 铝原子相对粒子分数已低于0.1.等离子体中一价到五价离子的相对粒子分数随着温度的增加都呈现出先增大后减小的变化规律.各价态离子的相对粒子分数峰值所对应的温度随着电离态的增大而逐渐升高.在温度为 1.0 × 104—3.0 × 104K、密度为 1.0 × 10–4—3.0 g/cm3的范围内, 等离子体中重粒子组分主要由 Al0, Al1+, Al2+和 Al3+组成, 更高电离态离子可忽略不计.

图3 等离子体中不同粒子的相对粒子分数在密度为0.1 g/cm3时随温度的变化Fig.3.Dependence of relative particle fraction of different particles on temperature at density of 0.1 g/cm3.

采用非理想Saha方程计算得到的铝等离子体在温度为10000和15000 K时平均电离度随密度的变化如图4所示.德拜理论修正的Saha方程、Ecker-Kröll (DmEK)模型[33,34]、FEOS 模型[35]和COMPTRA04程序[19]等计算的结果以及实验测量的平均电离度数据[36]也列于图4中进行对比.从图4(a)可以看出, 本文计算结果、德拜理论修正的Saha方程、DmEK模型和COMPTRA04程序计算的平均电离度分别在密度为0.63, 0.02, 1.71和0.21 g/cm3时出现突然上升的现象.德拜理论修正的Saha方程在耦合系数小于1的范围内计算结果处于比较合理的范围内, 但在耦合系数大于1的高密度区域, 其计算结果存在比较大的偏差.DmEK模型在小于临界密度区域内能够描述电离机制由热电离向压致电离的转变, 但平均电离度突变时的密度较高.COMPTRA04程序计算的平均电离度在高密度区域最高达到2, 本文计算结果显示平均电离度在10000 K时最高可达到2, 而在15000 K时可达到3.平均电离度随密度的突然增加导致采用过小的体积步长计算压强等热力学量时出现波动, 因此在该温度、密度区域内采用稍大的体积步长来获得平滑的压强曲线.

物态方程模型中非理想特性对非理想Saha方程准确描述等离子体电离机制随密度的增加由热电离到压致电离的转变有非常重要的影响.图5是不同温度下铝等离子体物态方程中非理想部分的库仑相互作用、排斥体积作用、极化作用和等离子体中电子以及原子的相对粒子分数随密度的变化.随着密度的增加, 等离子体由弱耦合区域过渡到强耦合区域, 非理想效应增强.

图4 不同模型计算的铝等离子体平均电离度在不同温度下随密度的变化 (a) 10000 K; (b) 15000 KFig.4.Average ionization degree of aluminum plasma calculated by different models as a function of density at different temperatures: (a) 10000 K; (b) 15000 K.

图5 不同温度下非理想效应自由能密度以及电子、原子相对粒子分数随密度的变化 (a) 库仑相互作用; (b) 排斥体积作用;(c) 极化作用; (d) 电子、原子相对粒子分数Fig.5.Free energy density of different non-ideal effects and relative particle fraction for electrons and atoms as a function of density at different temperatures: (a) Coulomb interaction; (b) excluded volume effect; (c) polarization effect; (d) relative particle fraction for electrons and atoms.

从图5可以看出, 排斥体积作用和极化作用等非理想部分远小于库仑相互作用.因此, 温稠密铝等离子体的非理想特性由库仑相互作用主导, 这与Apfelbaum[18]在其文章中所述结论一致.由(11)式可知, 非理想自由能密度的急剧变化会导致粒子化学势的变化, 进而影响有效电离能和等离子体组分.各非理想特性随密度、温度的变化趋势并不一致.库仑相互作用随着密度的增加而增强, 但随温度的变化并不明显.如图5(b)所示, 排斥体积作用在密度小于1 g/cm3的区域内随密度的增加呈现先增强后减弱的变化规律.当温度为15000 K 时, 排斥体积作用在 0.54 g/cm3时达到最大值.此时, 铝原子相对粒子分数已很低, 等离子体平均电离度接近1, 处于平均电离度迅速增大的密度区域.极化作用在密度为0.42 g/cm3时最强,随着密度进一步增加, 铝原子相对粒子分数降低导致极化作用减弱.排斥体积作用和极化作用的峰值所对应的密度随着温度的升高稍有降低.图5(d)显示在温度15000 K时铝原子相对粒子分数在0.03 g/cm3时达到最大值 0.27.随着温度升高, 铝原子相对粒子分数不断降低.当温度为30000 K时, 铝原子相对粒子分数随密度变化的峰值低于0.1.在平均电离度随密度迅速增大的高密度区域,铝原子几乎全部被电离.

稠密等离子体处于部分电离等离子体状态, 很高的重粒子数密度使得电离产生的自由电子可占据的自由空间减少, 而带电粒子也会使原子和离子自身的化学势发生改变, 重粒子之间的空间也因此具有较低的势能, 这会导致离子的能级结构发生变化, 出现电离能下降的现象[34,37].为了获得电离能的分布, 本文计算了15000 K时各粒子相互作用部分化学势随密度的分布, 如图6所示.从图6可以看出, 粒子的非理想化学势随密度增加而增大.当密度达到 0.5 g/cm3时, 由于离子密度增加, 非理想效应增强, 价态较高的离子非理想化学势的数值增大趋势明显.

图6 温度为 15000 K 时等离子体中不同粒子化学势的非理想部分随密度的变化Fig.6.Dependence of non-ideal chemical potential of particles on density at temperature of 15000 K.

粒子的非理想部分化学势决定了电离能下降,本文模型计算的铝等离子体在温度为15000 K时不同价态离子电离能降低随密度的变化如图7所示.DmEK模型计算结果也画于图中对比.从图7可以看出, 不同价态离子电离能的降低均随密度增加而增大, 且价态越高的离子电离能下降程度越大, 这与Son等[38]在其文章中所述结论是一致的.非理想效应的增强导致有效电离能减小, Al1+,Al2+和Al3+的有效电离能分别在密度为0.11, 0.53和0.54 g/cm3时接近0, 过小的电离能阈值会使原子和低价态离子进一步电离.当密度高于0.54 g/cm3时, Al2+全部电离为Al3+, 而这导致铝等离子体的平均电离度在临界密度区域出现突然上升的趋势.DmEK模型计算的电离能在密度小于0.1 g/cm3的区域与本文计算结果符合较好, 但低估了各级粒子电离能的降低, 从而平均电离度在密度约为3 g/cm3时才出现增大的现象.

图7 温度为 15000 K 时不同模型计算的电离能的降低∆E随密度的变化(黑色曲线, 非理想Saha方程; 红色曲线,DmEK模型)Fig.7.Depression of ionization potential calculated by different models as a function of density at 15000 K.Black lines correspond to nonideal Saha equation; red lines correspond to DmEK model.

5 结 论

温稠密等离子体态是高能量密度物理领域中最为常见的物理状态, 准确的描述温稠密态物质的特性是非常复杂的物理问题.本文在理想自由能的基础上考虑了库仑相互作用、排斥体积作用和极化作用等重要的非理想特性, 结合非理想Saha方程开展了温稠密铝等离子体物态方程和电离平衡研究, 获得了压强、等离子体粒子组分等数据, 并与其他模型计算结果及相关实验数据作了对比、分析.不同温度、密度下的非理想自由能曲线表明,库仑相互作用随着密度增大而增强, 但对温度并不敏感; 在小于 1 g/cm3的密度范围内, 排斥体积作用和极化作用都随着密度增加呈现先增强后减弱的变化规律.多粒子体系物态方程与非理想Saha方程迭代求解获得了等离子体中不同电离态离子分布.在温度为 10000 K 时, 高密度等离子体平均电离度为2, 而在温度为15000 K时, 高密度等离子体平均电离度为3.通过等离子体中粒子非理想化学势分析了由于非理想特性导致的不同价态离子的电离能下降.计算结果显示, 在某些温度、密度下, 离子的电离能下降与其电离能相当,从而导致有效电离能接近 0.在 1.0 × 104K < T <3.0 × 104K 的温度范围内, 当等离子体密度达到0.54 g/cm3时, Al3+的有效电离能降低至接近 0,从而导致Al2+电离成为Al3+, 等离子体平均电离度在这一密度区域内随着密度增加迅速地增大.本文模型能够计算密度为 1.0 × 10–4—3.0 g/cm3范围内的铝等离子体物态方程和粒子组分, 当密度进一步增加时, 迭代求解计算量增加, 迭代收敛性逐渐变差, 电离平衡方程变得难以求解, 并且过小的体积步长也会在临界密度区域内对压强的计算结果造成影响.今后将对上述问题开展进一步深入的研究.

附录A 等离子体物态方程中库仑相互作用项各部分表达式

本文采用Pade近似公式描述库仑相互作用, 电子-电子相互作用、电子-离子相互作用和离子-离子相互作用的表达式如下[39]:

其中 f0—2, e1—2, q0—2以及 Q0—5的表达式为

其中 z4= 0.827702, z5= 1.131, z6= 0.55045, 函数 qi(t) 可以通过(A17)—(A20)式计算