化学势在钨表面氢平衡浓度计算中的应用

吕广宏，李宇浩，张雪松，金 硕，周 苗，程 龙

(北京航空航天大学 物理学院，北京 100191)

热力学与统计物理是高等学校物理专业本科生的一门核心课程，目标旨在研究热运动规律以及与热运动相关的物性和宏观物质系统演化行为[1,2].相较于物理专业的其它核心课程，如理论力学、电动力学和量子力学等，热力学统计物理的研究对象为近乎于无穷多个粒子(～1023)所组成的系统.这就导致热力学统计物理的课程学习不能仅仅依赖数学公式的推导和简单的物理直觉，必须在理解其物理本质的基础上，建立完备的知识体系，并采用适当的数学方法近似求解.其中，化学势作为热力学与统计物理课程中一个重要而又抽象的物理概念，不同于传统的能量和温度，往往难以直接理解其物理含义[3].因此，在课程中如何通过引入实际相关算例，加深对化学势概念和物理本质的理解，这对于热力学统计物理的学习是十分必要且有益的.

算例引入：在未来核聚变反应堆中，钨(W)基面对等离子体材料会受到低能、大束流的氢(H)同位素等离子体辐照[4].在长时间辐照后，H同位素在W表面的吸附和脱附将达到动态平衡，此时表面H同位素的浓度即为平衡浓度.那么如何求解不同条件(温度、压强)下H在W表面的平衡浓度？这一算例源于聚变材料领域的实际问题，但并不要求很高的理论物理基础和专业背景，可以仅从基本的热力学统计物理概念出发，按照预先设计的关键问题和已知物理量，结合简单的微积分知识来研究.同时，在研究过程中会涉及计算物理方面的知识，教学的同时开拓学生视野，提高其自主学习的兴趣.

动态平衡是指H在W表面的吸附和脱附达到平衡状态，W材料本身仅作为基底，不参与平衡过程.因此，整个等离子体辐照过程可以近似为单元两相系所构成的孤立系统，即H在W表面的H-W系统和等离子体辐照H脱附的环境H2气体系统.根据单元系的复相平衡条件，两相系统达到平衡时其对应的温度(T)、压强(P)和化学势(μ)必须分别相等.由于环境中H2分子数量巨大，其压强和温度几乎保持不变，与W-H系统等同，两相系统达到平衡的关键在于化学势相同，而化学势又与体系的吉布斯自由能直接相关.因此，为了研究H在W表面的平衡浓度，必须首先给出环境和表面系统的吉布斯自由能，进而获得对应的H化学势.

1 环境系统中H原子的化学势

显而易见，在整个双相系统中总H原子个数保持不变，即当W表面吸附1个H原子时，环境中减少1/2个H2分子.根据开系(物质的量发生变化)热力学基本方程，环境H2分子系统的吉布斯函数可以写成

(1)

(2)

即在温度和压强保持不变的情况下，化学势等于增加1个H原子(或1/2个H2分子)时吉布斯函数的变化.

根据热力学基本方程，环境H2分子体系的吉布斯自由能由下式得到

GE=FE+PV

(3)

其中F是体系的亥姆霍兹自由能，P和V分别为压强和体积.在标准大气压和100 K以上时，H2分子近似看作理想气体.此时，体系的亥姆霍兹自由能可以写成

FE=F0+Ft+Fr+Fvib

(4)

其中F0是0 K下体系的亥姆霍兹自由能，可作为常数项，仅与参考点的选取有关.若以真空中单个H原子作为参考点，则有F0=-458.1 kJ/mol，对应于H2分子的结合能；若以H2分子作为参考点(这也是本文选取的参考点)，则有F0=0.上式右边的第2、3和4项分别为H2分子平动、转动和振动对于亥姆霍兹自由能的贡献.

对于平动项，其能量表达式可以写为

(5)

据此，平动的配分函数为

(6)

其中V为H2分子的体积，m为H2分子质量，h为普朗克常量数，k和T分别为玻尔兹曼常数和温度.因此，平动对亥姆霍兹自由能的贡献可以由下式得到，即

Ft=-NkTlnZt+kTlnNH2!

(7)

(8)

根据微观粒子全同性原理，H2分子的转动状态与两个H核的自旋状态有关.如果H核自旋平行，则转动量子数l取奇数，称为正氢(o-H2)；如果H核自旋反平行，l取偶数，称仲氢(p-H2).一般而言，正氢和仲氢含量之比为3:1，即环境H2分子系统为正氢和仲氢的混合物.此时，H2分子的转动能级可以写为

(9)

其中I为H2分子绕质心的转动惯量.对于o-H2和p-H2，其配分函数分别写为

(10)

(11)

(12)

根据之前的文献[8]，我们将配分函数取前8项，即l≤7.

对于H2分子中原子的相对振动(此项对自由能的贡献远低于平动和转动项的贡献)，可近似作为简谐振动.H2分子的振动能级可以写为

(13)

其中ω为振动圆频率.振动配分函数可以写为

(14)

因此，振动对于亥姆霍兹自由能的贡献为

(15)

(16)

(17)

表1 标准大气压下不同温度时单个H原子的化学势[6,7].

2 表面系统中H原子的化学势

类似的，对于W-H表面系统，单个H原子的化学势同样可写为

⑨张丹丹:《社区自治的特征:偏态自治和无序自治——社区自治空间有限性的原因》，《华东理工大学学报》(社会科学版)2015年第2期。

(18)

(19)

其中内能U和熵S均会随表面吸附H原子个数的变化而改变.因此，H在W表面的化学势是

(20)

其中第1项是平均1个H原子从环境中到吸附在W表面所导致的内能变化，可通过计算机模拟的方法得到.W作为一种体心立方金属(立方晶胞8个顶角和体心均有原子占据)，具有周期性和对称性的晶格结构.由于(110)和(100)表面分别是W中最稳定和最简单的表面，在此仅以(110)和(100)表面为例进行讨论.最近，Piazza等人[8,9]通过基于密度泛函理论的第一性原理计算，给出了H原子在W表面的平均吸附能(Ead).吸附能的定义为

(21)

其中Esurf和Esurf+H分别是纯W表面和含H原子W表面的体系总能量，EH2是真空中单个H2分子的能量.需要指出的是，在此我们仅考虑H原子吸附在W最表层上的情况，这是由于H原子难以溶解或吸附在W的近表层[10].吸附能为负意味着H原子倾向于从环境中吸附到W表面.根据Piazza等人的计算结果[8,9]，单个H原子吸附在W(110)和(100)表面的平均吸附能为-0.71 eV和-0.75 eV，此即为单个H原子吸附到W表面系统所导致的内能变化.

图1给出了H原子在W(110)和(100)表面吸附达到饱和时，H原子的占位.需要指出的是，我们在此采用了周期性边界条件，且只考虑第1层W原子.图1(a)中4个H原子位于棱上，为相邻2个晶胞共享，即每个相当于1/2个H原子，因此1个晶胞内H原子数为2个；第1层含5个W原子，其中顶角W原子由4个晶胞共享(每个占1/4)，而中心W原子完全属于该晶胞，因此1个晶胞内W原子数也为2个.据此，当H在W(110)表面达到饱和时，对应的H和W原子比为1:1.类似地，图1(b)中4个H原子也位于棱上，即与(a)图一样，1个晶胞内含有2个H原子；第1层4个W原子均占据顶角，单个晶胞占比为1/4，所以1个晶胞内W原子数仅为1个. 因此，当H在W(100)表面达到饱和时，对应的H和W原子比为2:1.

熵为系统在某一宏观状态下可能出现的微观分布状态数，反应了系统的混乱程度.吸附系统中熵的计算只考虑被吸附原子与第1层表面的相互作用.一般而言，金属体系的熵包括构型熵和振动熵两个部分.其中构型熵反映了材料内部原子的无序度，即混乱程度的度量，与系统状态数的对数成正比.振动熵反映了晶格振动频率的无序度，晶格振动可以通过计算声子谱振动频率来获得.计算结果表明，对于H在W表面的吸附过程，振动熵与构型熵相比可忽略不计，因而本文只考虑构型熵的贡献.同时，相较于振动熵，构型熵也更适合热力学统计物理的案例教学.

图1 H原子占据(a) W(110)表面和(b) W(100)表面达到饱和时的原子占位示意图，其中灰色和白色的大圆分别表示表面第1层和第2层的W原子，黑色小圆表示吸附在W表面的H原子.在(110)表面，H和最表层W原子之比为1:1，而在(100)表面，H和最表层W原子之比为2:1.

(22)

(23)

此时采用斯特林公式可将上式化简为

(24)

可以看到，体系的构型熵仅取决于占据W表面的H原子数与W表面所有可能的占位数之比，而与总占位数的绝对值无关.基于式(20)、(21)和(24)，可以得到

(25)

3 基于两相的化学势计算H平衡浓度

根据复相平衡条件(温度、压强和化学势均相等)以及式(17)、(25)，我们可以得到

(26)

图2和图3展示了不同条件(温度和压强)下，H在W(110)和(100)表面吸附的平衡浓度.在此，平衡浓度是指W表面所有可能的H原子占位被占据的可能性，即x.可以看到，在300～500 K温度和1～10000 Pa压强范围内，W表面H平衡浓度基本不随温度和压强的变化而改变，始终接近于1(趋于定值)，即基本W表面所有可能的占位均被H原子占据.这是由于在低温下(<500 K)吸附能主导了H在W表面的吸附.根据公式(26)，H平衡浓度强烈依赖于H在W表面的吸附能和环境H2分子系统中H原子的化学势.其中，前者近似为常数项[(110)表面：-0.71 eV；(001)表面：-0.75 eV]，不随温度和压强的变化而改变.根据式(17)，环境系统中H原子的化学势始终为负值，且随温度的增加(或压强的下降)单调下降.而在300～500 K和1～10000 Pa的范围内，环境系统中H原子化学势绝对值明显低于H在W表面吸附能的绝对值，即后者决定了氢在钨表面的平衡浓度.由于吸附能为负值，意味着H原子倾向于从环境系统吸附到W表面系统，对应于W表面所有可能的占位均被H原子占据.随着温度的进一步上升，W表面H平衡浓度随之下降；在1000 K以上时，H平衡浓度均低于0.2，且随压强的下降逐渐降低，如图2和图3所示.这说明温度升高和压强下降均会抑制H在W表面的吸附.这是由于温度升高和压强降低都将导致环境系统中单个H原子化学势的降低，使其逼近甚至低于H在W表面的吸附能，进而显著促进H原子的脱附，降低H在W表面的平衡浓度.进一步对比H在(110)面和(100)面的平衡浓度可以发现，在相同条件下，H在(100)面的平衡浓度略高于其在(110)面的平衡浓度，表明H更倾向于吸附在W(100)表面.这对应于H在(100)面的吸附能低于其在(110)面的吸附能.

图2 W(110)表面H吸附的平衡浓度.

图3 W(100)表面H吸附的平衡浓度.

需要指出的是，即便H在W(110)和(100)面的平衡浓度均为1，其对应的H原子密度却并不相同，这是由于两个表面的W原子密度和H可能占位的不同(见图1).其中，对于W(110)表面，当H吸附浓度达到饱和时，对应的H原子密度为～1.39×1019m-2；对于W(100)表面，饱和时对应的H原子密度为～1.97×1019m-2，与实验上300 K下测量的结果基本一致[11].

有趣的是，即便不从单元系的复相平衡条件出发，而是将整个吸附/解吸附看成一个等温等压的孤立系统，基于热动平衡判据(δG=0)，我们同样可以得到H在W表面的平衡浓度.此时，整个体系的吉布斯自由能可以表示为

G(n,P,T)=G(0,P,T)+ΔGH2(n,P,T)+

ΔGW-H(n,P,T)

(27)

其中第一项为没有H原子吸附在W表面时，整个体系的吉布斯自由能.第二项和第三项分别为有n个H原子吸附在W表面时环境系统和W-H表面系统吉布斯自由能的变化.接下来，类似于步骤一和步骤二的处理，根据热动平衡判据，同样可以得到W表面H平衡浓度的表达式.

另一种思路则是从玻尔兹曼分布函数出发，同样将整个吸附/解吸附看成一个等温等压的孤立系统.定义两种不同的H原子状态，即H在W表面吸附的状态(na)和H从W表面脱附的状态(nd).两种状态的H原子数目之和保持不变，即为W表面H原子所有可能的占位数.据此，基于玻尔兹曼分布函数，我们可以知道两种分布状态的H原子个数之比取决于两种状态下H原子的能量之差，即有

(28)

其中Ea和Ed分别是吸附和脱附时H原子的能量状态.根据定义，我们可以得到

(29)

(30)

关于H平衡浓度的计算可同样应用于材料体内.材料体内H可能的占位有2类：间隙位置和空位位置.一般而言，材料中间隙位置很多，但化学势较高；而空位位置较少，但化学势很低.因此，在一定温度和压强条件下，材料体内同样存在H平衡浓度，即H在材料间隙和空位位置的吸附和脱附达到动态平衡.这部分拓展已在作者前期的研究论文中进行了详细的介绍[12]，在此不再具体展开.

综上所述，本文针对化学势这一热力学统计物理课程中重要物理量，介绍了W材料表面H平衡浓度这一典型算例.通过分别考察环境H2分子系统和W表面系统中H原子的化学势，利用所学过的热力学统计物理知识，最终给出了不同温度和压强条件下W材料表面H平衡浓度.这不仅能帮助学生对化学势等相关重要知识点的深入理解，强化了热力学统计物理知识体系的系统构建，并且非常有助于学生解决实际问题能力的提升.