汪泽涛,郭凯伦,王成龙,张大林,田文喜,秋穗正,苏光辉
(西安交通大学 核科学与技术学院,陕西 西安 710049)
固态热管反应堆具有紧凑性、简易性、模块化等多种优点,可广泛应用于多种场景,在未来的能源市场有着巨大潜力[1]。高温碱金属(钠、钾、锂等)热管是热管堆内的关键部件,依靠内部工质的相变及循环,可实现堆内非能动热量的传输。
分子动力学作为一种微观模拟手段,在研究液态薄膜气液交界处的蒸发与冷凝中已有诸多应用[10-12],但研究对象主要为水、氩以及一些液态有机物,有关液态碱金属的研究十分稀少[13]。
因此,为深入理解高温热管内部的工质的蒸发与冷凝机理,本研究使用分子动力学软件LAMMPS,基于热管堆中常用的高温钠热管的启动与运行工况[14-15],模拟4组不同工况下的液态钠薄膜的平衡态蒸发,并求解MAC。随后加入氩气原子,作为非凝结性气体进行模拟,考察非凝结性气体对液态钠薄膜蒸发和MAC的影响。
如图1所示,利用LAMMPS建立了尺寸为8.4 nm×8.4 nm×58.2 nm的模拟区域,区域内部上下侧均设定了9层金原子作为衬底壁面,在壁面上部放置了8.4 nm×8.4 nm×8 nm的固态钠块,模拟区域的x、y方向均采用周期性边界,z方向采用固定边界。在考察非凝结性气体的影响时,除模拟区域中央设定了3 811个氩原子(在模拟中蒸发为气态原子作为非凝结性气体)外,其余设置均保持不变。
图1 模拟体系
原子晶格形式均采用FCC(面心立方),钠、金、氩原子的晶格常数分别为5.94、4.08、5.2[16]。原子间的相互作用采用12-6 Lennard-Jones势[17],其分布形式如图2所示。
图2 12-6 Lennard-Jones势函数
12-6 Lennard-Jones势的表达式如下:
(1)
式中,ε、σ分别为势阱深度、特征长度。势阱深度是指两原子间相互吸引作用的最大值,也是原子吸引作用和排斥作用的界限;特征长度是指两原子间不存在相互作用时的距离。各原子间的势参数[13,16,18-19]列于表1。
表1 原子相互作用势参数
根据钠热管的运行温度范围,共设定了600、700、800、900 K 4组工况,进行平衡态蒸发模拟。其中,前两组工况主要针对钠热管的启动阶段,后两组工况主要针对钠热管的正常运行阶段。
模拟步长设为1 fs,截断半径为3.5σNa-Na。除最外侧固定壁面外,使体系内其余所有原子在NVT系综下运行6 ns,达到每组工况下所指定的温度范围,实现平衡态蒸发。随后在NVE系综下运行1 ns,用于输出有关后处理信息。用OVTIO和MATLAB进行可视化演示和数据后处理。在分子动力学中,系综是指具有相同宏观状态,不同微观状态的热力学体系的集合[17]。NVT为定容定温且体系原子数目不变的系综,NVE为定容且体系总能与原子数目不变的系综。
如图3所示,分子运动理论[20]认为,气体原子在到达气液交界后,一般会发生3种行为:被交界面反弹;冷凝为液体原子;先冷凝为液体原子,然后又蒸发为气体原子。
图3 气体原子在气液交界的3种行为
依据Liang等[19]的方法,在距离气液交界3.5σNa-Na处设定虚拟面(图4),利用式(2)~(4)先后确定了气体原子的平均速度um和特征时间Δt。其中,Δt指气体原子穿入虚拟面,又穿出虚拟面所要经历的平均时间。从统计起始时刻开始,经历Δt后,穿入虚拟面的气体原子中被冷凝的数目比例即为MAC(式(5)中记为α,Nref特指穿出虚拟面的原子数目)。4组工况下的特征时间分别为12.2、11.3、10.6、9.9 ps。
图4 统计原理示意图
d=3.5σNa-Na
(2)
(3)
Δt=2d/um
(4)
(5)
分子动力学模拟的有效性,在于模拟中所选取的势函数及其参数能否准确反映出流体的物性特征。如图5所示,单独建立了8.4 nm×8.4 nm×20 nm的模拟区域,区域中央设定了12 nm的液态钠薄膜,利用NVT系综实现了其在600~910 K的平衡态蒸发。
图5 势函数验证模型
获取了饱和钠蒸气在不同温度下的压力。利用式(6)、(7)对温度、压力进行无量纲处理(式(6)中kB指玻尔兹曼常数),与NIST专设的12-6 Lennard-Jones流体物性库[21-22]中的压力进行了比较(图6),误差在3%~25%。氩的势函数及参数已被广泛应用于研究液态氩蒸发与冷凝的分子模拟中[10,19,23-24],有效性可以保证。
图6 验证结果对比
以上结果及论述表明,本研究所选用的势函数及其参数具有准确性,所进行的模拟具备有效性。
T*=kBT/ε
(6)
P*=Pσ3/ε
(7)
图7为无非凝结性气体和含非凝结性气体影响下的平衡态蒸发图像。随着温度的升高,钠的平衡态蒸发愈发剧烈,模拟区域主体逐渐被气体原子占据。除少量氩原子在钠液膜内部,非凝结性气体主要集中在钠蒸气区。
4组工况下,无非凝结性气体的MAC和非凝结性气体影响下的MAC汇总如图8所示。总体来看,从工况1到工况4,随着温度的升高,MAC出现了下降趋势。除工况2外,非凝结性气体的存在会使MAC降低。不存在非凝结性气体时,工况3的MAC相比工况2出现了较小回升。交界面气态钠原子的反弹比例如图9所示,可看出,从工况1到工况4,随着温度的升高,交界面气态钠原子的反弹比例进一步升高。当含有非凝结气体时,这种反弹作用则进一步加剧。
图9 4组工况下气液交界位置气体钠原子的反弹比例
如图10、11所示,获取4组工况下的液膜和邻近蒸汽区域的温度、势能分布,进行分析。温度分布反映原子的动能大小(统计热力学中,ke=1.5kBT)。势能分布反映原子间的相互作用,其绝对值越大,表明作用越强烈(负值表明吸引作用)。从微观角度来看,蒸发是液体内部原子动能增强,克服周围原子的吸引作用,跃出气液交界成为气体原子的过程。从工况1到工况4,平衡态蒸发的温度升高,原子的动能增加。同时,液膜内部的相互吸引作用逐步减弱(-0.4~-0.3 eV变化至-0.2~-0.1 eV),且交界面的气相钠原子的反弹逐步加剧。这三者的综合作用,使得气体原子不宜凝结,导致MAC呈现出下降趋势。工况3中,随着液膜的减薄,壁面的吸引作用开始凸显(-0.7~-0.6 eV),使工况3的MAC出现较小回升,但不影响总的下降趋势。
图10 液膜及邻近钠蒸气区域温度分布
图11 液膜及邻近钠蒸气区域势能分布
由于非凝结性气体主要集中在液膜外的钠蒸气区域,仅有少部分在液膜内部,交界面反弹加剧,穿入虚拟面的气相钠原子中被反弹的比例显著增加。此外,在前两组工况下,液膜中存在的非凝结性气体也使液膜内部的吸引作用分别减弱至-0.3~-0.2 eV和-0.2~-0.1 eV。后两组工况下,由于液膜十分稀薄,非凝结性气体几乎全部集中在气相区,对液膜内部原子间的相互吸引作用影响不大,主要起到的是进一步加深界面反弹的作用,气相钠原子反弹比例相比前两组工况显著提升,分别为0.806 8和0.825 3。因此,在工况1、3、4下,非凝结性气体的存在,使得冷凝不易发生,降低了MAC。
工况2中,以上两种影响虽然存在,但MAC出现了较小增加。非凝结性气体的存在,使工况2下液膜的厚度降低了2 nm。这使得更多的气态钠原子充斥在蒸气空间,样本量发生轻微改变,进而导致非凝结性气体存在条件下统计所得MAC的较小回升。
上述结果为钠热管的相关数值模拟提供了新的参数依据,也表明了在启动初期(600 K)和钠热管正常运行阶段(800 K和900 K),非凝结性气体会阻碍气液界面的冷凝,而在启动中后期(700 K),非凝结性气体的存在却起到了相反作用。非凝结性气体对热管性能有着重要影响。大量宏观角度的研究[25-27]认为,一定量的非凝结性气体一方面会促进钠热管的启动进程;另一方面,当钠热管处于正常运行阶段,非凝结性气体滞留在气液界面和冷凝段末端,阻碍热管内的相变换热。本文的结果与这些宏观角度的观点吻合。同时,本工作也表明了非凝结性气体在钠热管启动中后期对气液交界相变的特殊影响作用,该作用需进一步深入研究。
本研究利用分子动力学软件LAMMPS,以热管堆中常用的高温钠热管的启动和运行工况为基准,模拟了4组不同温度下液态钠薄膜的平衡态蒸发,求解了每组工况下的MAC。并以气体氩原子作为非凝结性气体,考察了非凝结性气体对MAC的影响,得到主要结论如下。
1)在600、700、800、900 K下的液态钠薄膜平衡态蒸发中,MAC分别为0.388 6、0.211 9、0.261 5、0.241 6。存在非凝结性气体时,MAC分别为0.282 9、0.254 3、0.129 5、0.107 2。这些结果为钠热管的数值模拟工作提供了新的参数借鉴。
2)温度升高,使液膜内部原子动能增加,液膜内原子的相互吸引作用减弱,使液体原子更易蒸发。同时,也使得气液交界面的反弹加剧,冷凝不易发生。
3)非凝结性气体的存在,一方面会进一步加深交界面原子间的反弹作用,另一方面也会使液膜内部的原子相互吸引作用减弱。这使得冷凝不易发生,造成MAC降低。600、700 K下,这两方面的影响都很显著。而在较高温度下(800、900 K),主要体现为第1种影响。
4)微观角度下,非凝结性气体对气液界面的相变影响将会随温度范围而改变。在钠热管启动初期(600 K)和正式运行(800 K和900 K)范围,非凝结性气体起到阻碍作用。其在启动中后期(700 K),起到促进作用。这一结果与许多宏观实验研究相吻合,也表明非凝结性气体在钠热管启动中对后期的影响需进一步深入研究。