液态锂在α-Fe （100）表面滑移特性的分子动力学研究

Molecular dynamics study of slip properties of liquid lithium on α-Fe （100） surface

摘要

采用分子动力学方法模拟液态Li在α-铁（Fe）纳米通道中的Couette流动，分析三角形表面粗糙单元对液态Li在α-Fe （100）表面滑移特性的影响。结果表明三角形两侧角度对滑移长度的影响并不显著，但是存在一个临界角度，当三角形两侧角度超过该临界值时，液态Li会在三角形前侧形成微尺度旋涡，显著增加流动阻力，进而引起滑移长度的突然大幅下降；三角形高度对滑移长度具有更加显著的影响，且同样存在一个临界高度，只有当三角形高度超过临界值时才会形成旋涡，否则无论如何改变三角形的角度都不会在三角形附近形成旋涡；温度对于临界角度和临界高度都有一定的影响，随着温度的升高，临界角度和高度都呈下降趋势。

Abstract

Liquid lithium （Li） has a wide range of applications in nuclear fusion test devices， but the slip characteristics of liquid Li on the surface of solid materials have not been fully studied. In this paper， the Couette flow of liquid Li in α-iron （Fe） nanochannels was simulated by the molecular dynamics method， and the effect of triangular roughness units on the slip characteristics of liquid Li on the α-Fe （100） surface was analyzed. The results show that the influence of the angle on both sides of the triangle on the slip length is not significant， but there is a critical angle， when the angle on both sides of the triangle exceeds the critical value， the liquid Li will form a microscale vortex on the front side of the triangle， which significantly increases the flow resistance， and then causes a sudden and large decrease in the slip length， the critical angle and the height show a downward trend. The height of the triangles has a more pronounced impact on slip length， and a similar critical height exists. Vortices only form when the height exceeds this critical value， regardless of changes to the angles. Temperature has a certain influence on these parameters. As temperature rises， both the critical angle and critical height decrease.

Graphical abstract

关键词

分子动力学 / 液态锂 / Couette流动 / 滑移长度

Key words

molecular dynamics / liquid lithium / Couette flow / slip length

随着能源危机问题日趋严重，磁约束受控核聚变技术受到越来越广泛的关注。在当前的磁约束核聚变试验装置中，直接面向等离子体的第一壁材料在运行过程中需承受多种极端物理场的作用，包括超高热流密度的热负荷、高能等离子体轰击、中子辐照等^［1-6］。这些极端条件使得第一壁材料的选择非常有限，目前最好的选择是金属钨（W），但是在长时间工作状态下仍然面临熔化、腐蚀、开裂等问题^［7-10］。除固体材料外，采用流动的液态金属作为面向等离子体材料是另一种选择，其基本原理是在固态材料表面覆盖一层液体金属薄膜，通过液膜的持续流动，避免固体材料与等离子体的直接接触，从而有效地解决固态第一壁材料所面临的大部分问题^［11］。同时，实验证实采用液态锂（Li）作为面向等离子体材料还可以显著提高等离子体的稳定性和约束性能^［12-15］，荷兰的研究团队采用3D打印技术制作了多孔毛细结构的面向等离子体部件，利用毛细力将液态Li吸附到材料表面，结果表明该部件在16 MW·m^-2热流密度的热负荷作用下仍可以保持结构完整性^［16］。

目前，已有不少研究人员针对液态Li第一壁开展了相关数值模拟研究，在宏观尺度上，通常是通过求解Navier-Stokes方程来模拟液态Li在聚变环境下的流动行为^［17-20］。在这些模拟中，通常会采用无滑移边界条件。而实际上，真实的液固界面滑移状态可能存在不同程度的正滑移和负滑移^［21-23］。具体滑移特性与原子间的相互作用、温度、表面粗糙度、晶格匹配度等多种因素有关^［23］。在微观尺度上，目前的研究焦点大多集中于液态Li在金属表面的润湿性和铺展性方面^［24-26］，有关液态Li在固体材料表面的滑移特性，还需要开展进一步的研究。

我们前期已经开展了液态Li在金属W表面的滑移特性研究^［27］，重点讨论了矩形粗糙单元、壁面应力状态、壁面温度等因素的影响。考虑到真实材料表面的粗糙单元通常具有复杂的几何形状，而三角形更具一般性，通过改变三角形两侧的角度，可以更好地模拟真实的粗糙表面。因此，本文采用分子动力学方法，模拟液态Li在α-Fe通道内的Couette流动，讨论三角形粗糙单元对液态Li在α-Fe （100）表面滑移特性的影响规律。研究结果可以为更好地理解液态Li在不锈钢微通道内的流动提供科学参考。

1 模型和方法

本文中的分子动力学模拟全部采用开源软件LAMMPS^［28］完成。Li-Li、Fe-Fe和Li-Fe原子之间的相互作用采用Liu等^［29］专门针对Li-Fe系统开发的嵌入原子势进行描述，时间步长取为0.001 ps。原子构型如图1所示，其中红色原子为Fe，蓝色原子为Li。Fe壁面的晶格取向为：x-［100］，y-［010］，z-［001］，上壁面的厚度为5a，其中a=2.83 Å（1 Å=0.1 nm，下同）^［29］，为体心立方α-Fe的晶格常数。下壁面最小厚度同样为5a，在壁面的中心部位加入三角形的粗糙单元，并采用3个参数进行表征，分别为迎风面角度α、背风面角度β和高度k。如图1所示，通过控制这3个参数，模拟不同粗糙单元的壁面。整个模拟区域沿z方向的高度为50a，在上下壁面之间根据液态Li在不同温度下的数密度随机填充Li原子。模型沿y方向的尺寸为20a，沿x方向的尺寸则根据三角形粗糙单元的几何尺寸进行调整，同时保证三角形左右两侧角点之间的距离为40a。图1中所示的原子模型沿x方向的总长度为76a，共包含30 400个Fe原子、70 198个Li原子。

为了构建Couette流动，在x和y方向上采用周期性边界条件，在z方向采用自由边界条件。所有的Fe原子固定在晶格点位，不产生热振动，即处于0 K状态，仅提供粗糙壁面条件。在液态Li开始流动之前，系统在正则系综下弛豫0.5 ns，一方面使得液态Li达到需要的稳态温度，另一方面通过上壁面在z方向的自由移动，消除液态Li的初始压力。此后，保持下壁面不动，上壁面作为刚体沿x正向以U=\mathrm{3}\mathrm{\AA }·\mathrm{p}{\mathrm{s}}^{-\mathrm{1}}的速度运动，经过4 ns的模拟以后，形成稳定的Couette流动。

对于Couette流而言，如果采用无滑移边界条件，那么流体的速度将沿着通道厚度方向呈线性分布，而在微观尺度上，由于流体与壁面原子之间的相互作用，流速在壁面附近通常呈现非线性变化。如图2所示，根据流速的非线性特征，可以分为3种情况，即：负滑移、无滑移和正滑移。本文采用滑移长度表征这3种界面滑移特性，具体计算方法为：对通道中心部位的流速进行线性插值并外推到与壁面速度相同的点，该点与壁面的距离即为滑移长度。为计算液态Li的流速，将流体区域划分为厚度为2.83 Å的等厚度薄片，将薄片内所有原子的速度进行算术平均，作为薄片中心点处的流动速度，同时为消除热振动带来的随机误差，再对薄片速度进行4 ns的时间平均。

2 结果与讨论

2.1 迎风面角度对滑移特性的影响

为考察三角形粗糙单元的迎风面角度α对滑移特性的影响，模拟计算不同角度条件下液态Li在α-Fe （100）晶面上的滑移长度。三角形粗糙单元的高度k=41.035 Å，背风面角度β=35°，α的变化范围为26.57°～90°，对应的模型沿x方向的长度为249.04～141.5 Å。

图3所示为典型的液态Li流速随高度的分布曲线（α=45°）。可以看出，在壁面附近，液态Li的流速几乎为0，这是由于Li-Fe之间具有较强的相互作用，距离Fe表面最近的一层Li原子被紧紧地束缚在Fe表面，无法流动。从图1所示的原子构型图中也可以看出，在壁面附近Li原子呈现分层状态，表现出类固态特征。这种现象与壁面Fe原子的0 K条件并没有直接的关系，即使当壁面存在有限温度时，壁面Fe原子的热振动会引起Li原子势能和受力的变化，但是这种变化很小，Li原子仍然会呈现类固态的分层现象。随着与Fe表面距离的增加，液态Li的速度逐渐增加，在远离壁面的区域，液态Li流速与高度具有良好的线性关系，符合Couette流动的特征，将线性区（53.77 Å z 138.67 Å）插值到壁面流速，即可计算出滑移长度为-36.62 Å。此外，从图中还可以看出，上表面附近的非线性区域并不明显，这是由于下表面具有三角形粗糙单元，而上表面为光滑表面。

图4所示为滑移长度随迎风面角度α的变化曲线。可以看出，滑移长度全部为负值，表明液态Li在α-Fe （100）表面具有负滑移特征，Fe表面对液态Li的流动具有一定的阻碍作用。随着角度α的增加，滑移长度呈现出分段线性的规律：在\mathrm{26.57}°\le \alpha \le \mathrm{45}°角度范围内，拟合可得斜率为-0.043；在\mathrm{49.40}°\le \alpha \le \mathrm{90}°角度范围内，斜率为-0.084。显然在大角度区域内，滑移长度对角度变化更加敏感，即随着角度的增加下降更快。然而，实际上在两段线性区域内，滑移长度的变化并不显著，当α从26.57°增加到45°时，滑移长度仅降低了0.8 Å。而在两段线性区域之间，滑移长度却有明显的变化，当α从45°增加到49.40°时，滑移长度降低了3.3 Å。仔细观察不同工况下的原子构型，我们发现产生这一显著变化的原因是由于液态Li的流场结构发生了变化。

为直观展示流场结构，绘制了模拟区域的流线图。具体做法为：在x-z平面内将模拟区域划分为32×32的均匀网格，每个单元的尺寸为2.83 Å×2.83 Å，在每个单元内对原子速度进行空间和时间（4 ns）的平均，作为该单元的速度，得到速度场分布，进而可绘制出液态Li的流线分布，如图5所示。图5（a）为α=45°时的流场，可以看出液态Li平滑地流过三角形表面，表现为完全层流状态。实际上，根据我们的模拟结果，当\alpha \le \mathrm{45}°时，液态Li一直处于完全层流状态，这个阶段滑移长度随着角度α的增加而轻微下降，是因为下壁面的长度随着α的增加而增加，因此与液态Li接触的Fe面积增大，受到Fe表面吸附作用的Li原子数量也相应增多，从而导致流动阻力增加。当α=49.40°时，流线如图5 （b）所示，可以看出，此时在三角形粗糙单元的底部左侧出现了微尺度旋涡，正是该旋涡的产生引起了流动阻力的显著增加，对应图4中滑移长度的突然下降。此后，随着α角度的继续增加，旋涡不断增大（如图5（c）和5（d）所示），受影响的液态Li范围扩大，使得滑移长度以更快的速度降低。显然，在当前工况下，产生旋涡的临界角度为α=49.40°。

2.2 背风面角度对滑移特性的影响

与迎风面角度α相同，我们也模拟分析了背风面角度β对液态Li在α-Fe （100）晶面上滑移特性的影响，计算模型与前一小节所表述的类似，三角形粗糙单元的高度k=41.035 Å，迎风面角度α=35°，背风面角度的取值范围为\mathrm{26.57}°\le \beta \le \mathrm{90}°，对应的模拟区域长度为249.04～141.5 Å。

滑移长度随背风面角度β变化的曲线如图6所示。从图中可以看出，背风面角度对滑移长度的影响规律与迎风面角度非常相似，呈现分段线性的关系。在\mathrm{26.57}°\le \beta \le \mathrm{54.46}°角度范围内，滑移长度随着β的增加缓慢线性降低，当β=60.26°时，滑移长度突然显著下降，随后，滑移长度继续随着β的增加缓慢线性降低。通过观察流线图（见图7），可以确认在前一个线性阶段，液态Li处于完全层流状态，在后一个线性阶段，在三角形粗糙单元尾流区出现了不断长大的旋涡，而滑移长度的突变正是由于尾流区旋涡的突然出现所引起的。与迎风面角度不同的是，引起滑移长度突变的临界背风面角度为β=60.26°，比临界迎风面角度更大，滑移长度的突变值为3.89 Å，比迎风面角度引起的突变略大。

2.3 三角形高度的影响

已有研究表明，粗糙单元的高度通常对滑移长度有比较显著的影响^［30-31］。为揭示三角形粗糙单元高度对液态Li在α-Fe （100）晶面上滑移特性的影响，我们模拟了多种工况（不同迎风和背风面角度）下滑移长度随三角形高度的变化曲线，2种典型的结果如图8所示。图中黑色实线对应的工况为：α=45°，β=54.46°。在这种工况下，液态Li处于完全层流状态，三角形粗糙单元前后始终没有出现旋涡，也正因如此，滑移长度与三角形高度全程保持线性关系。图8中红色实线对应的工况为：α=49.40°，β=60.26°。此时在三角形的迎风面和背风面均应出现旋涡，然而，模拟结果显示，当三角形高度k\le \mathrm{29.715}\mathrm{\AA }时，并未出现旋涡，只有当k\ge \mathrm{35.375}\mathrm{\AA }时，才形成旋涡。表现在图8中的曲线上，在高度范围7.075\mathrm{\AA }\le k\le \mathrm{29.715}\mathrm{\AA }内，滑移长度随高度的增加线性降低，且与α=45°、β=54.46°的结果重合，当高度增大到35.375 Å时，滑移长度突然下降，表明此时旋涡开始形成，此后随着高度的继续增加，旋涡不断增大，滑移长度以更大的斜率与三角形高度保持线性关系。

除图8所示的2种工况以外，本文还模拟计算了其他多种迎风面角度和背风面角度的组合工况，例如当α=60.26°、β=49.40°时仅在三角形前方形成旋涡，当α=30.26°、β=74.05°时仅在三角形后方形成旋涡。而模拟结果表明，对于这些工况，产生旋涡的临界高度始终保持不变，即k=35.375 Å。

2.4 温度效应

在微观尺度上，温度反映了原子随机振动的强弱，当温度升高时，原子的动能增加，从而更容易逃离壁面的束缚。因此，可以预见，温度对于液态Li在α-Fe表面的滑移特性也应具有一定的影响。本文模拟计算了不同温度条件下滑移长度随α、β和k的变化关系，结果显示温度对于产生旋涡的临界角度和临界高度都有一定的影响。

图9所示为不同温度下滑移长度随迎风面角度α的变化关系。可以看出，在每一个固定温度下，曲线都具有分段线性的特征，而且，随着温度的升高，在相同的线性阶段内，滑移长度随着温度的升高而增大。这是由于温度升高意味着原子动能增加，流体黏性降低，流动阻力降低。然而，随着温度的升高，滑移长度发生突变的临界α角度不断减小，意味着在高温条件下，较小的迎风面角度也可以促使液态Li在三角形前方形成旋涡，从而增大流动阻力，减小滑移长度。例如当迎风面角度α=41.19°时，根据图9中的结果，当温度T=600和700 K时液态Li为完全层流，没有旋涡产生，因此随着温度的升高，流动阻力减小，滑移长度增加，而当温度T=800 K时，旋涡开始出现，滑移长度大幅下降，甚至低于600 K时的结果，此后随着温度的继续升高（900 K），流动阻力再次减小，滑移长度增加。

对于三角形粗糙单元的背风面角度β，温度也有类似的影响规律，即随着温度的升高临界角度β逐渐减小。图10展示了临界α和β随温度的变化关系，可以看出，2个临界角度都与温度呈近似线性关系。此外，模拟结果显示，产生旋涡的临界高度同样随着温度的升高而降低，如图11所示，其主要原因是原子随机热振动速度的增加使得旋涡更容易形成。

3 结论

本文采用分子动力学方法模拟液态Li在α-Fe纳米通道中的Couette流动，探讨三角形表面粗糙单元对液态Li在α-Fe （100）表面滑移特性的影响，主要结论如下：

1）滑移长度随着三角形迎风面角度的增大而减小，但是存在一个临界角度。当角度达到临界值49.40°时，液态Li会在三角形粗糙单元的前侧形成微尺度旋涡，显著增加流动阻力，进而引起滑移长度的突然下降，当角度大于或者小于临界值时，滑移长度都随着角度的增大线性减小，但是减小的幅度较小。三角形背风面角度对滑移长度有类似的影响规律，只是临界角度较大，为60.26°。

2）三角形高度对滑移长度具有更加显著的影响。当三角形两侧角度均小于产生旋涡的临界角度时，滑移长度随着高度的增加线性降低。当角度超过临界值时，滑移长度随高度的增加呈现分段线性的变化规律，即存在一个临界高度35.375 Å。低于该高度时，即使三角形角度超过临界值也不会产生微尺度旋涡；高于该高度时，由于旋涡的存在，滑移长度首先发生突降，随后随着高度的增加线性降低。

3）温度对于临界角度和临界高度均有一定影响，随温度升高，临界角度和高度均呈下降趋势。

参考文献

[1] Federici G， Biel W， Gilbert M R， et al. European DEMO design strategy and consequences for materials［J］. Nuclear Fusion， 2017， 57（9）： 092002. DOI： 10.1088/1741-4326/57/9/092002 .

[2] Eich T， Sieglin B， Thornton A J， et al. ELM divertor peak energy fluence scaling to ITER with data from JET， MAST and ASDEX upgrade［J］. Nuclear Materials and Energy， 2017， 12： 84-90. DOI： 10.1016/j.nme.2017.04.014 .

[3] You J H， Visca E， Bachmann C， et al. European DEMO divertor target： operational requirements and material-design interface［J］. Nuclear Materials and Energy， 2016， 9： 171-176. DOI： 10.1016/j.nme.2016.02.005 .

[4] Hasegawa A， Fukuda M， Nogami S， et al. Neutron irradiation effects on tungsten materials［J］. Fusion Engineering and Design， 2014， 89（7/8）： 1568-1572. DOI： 10.1016/j.fusengdes.2014.04.035 .

[5] Yu X G， Gou F J， Tian X. Molecular dynamics study of the effect of hydrogen on the mechanical properties of tungsten［J］. Journal of Nuclear Materials， 2013， 441（1/2/3）： 324-330. DOI： 10.1016/j.jnucmat.2013.06.018 .

[6] 刘松畅， 余新刚. 单晶钨辐照损伤的分子动力学研究［J］.中国科学院大学学报，2024，41（4）：452-460.DOI： 10.7523/j.ucas.2022.087 .

[7] Pitts R A， Bardin S， Bazylev B， et al. Physics conclusions in support of ITER W divertor monoblock shaping［J］. Nuclear Materials and Energy， 2017， 12： 60-74. DOI： 10.1016/j.nme.2017.03.005 .

[8] You J H， Visca E， Barrett T， et al. European divertor target concepts for DEMO： design rationales and high heat flux performance［J］. Nuclear Materials and Energy， 2018， 16： 1-11. DOI： 10.1016/j.nme.2018.05.012 .

[9] Kallenbach A， Bernert M， Dux R， et al. Impurity seeding for tokamak power exhaust： from present devices via ITER to DEMO［J］. Plasma Physics and Controlled Fusion， 2013， 55（12）： 124041. DOI： 10.1088/0741-3335/55/12/124041 .

[10] Li M Y， You J H. Design options to mitigate deep cracking of tungsten armor［J］. Fusion Engineering and Design， 2017， 124： 468-472. DOI： 10.1016/j.fusengdes.2017.01.015 .

[11] 陈然， 王增辉， 倪明玖. 强磁场对导电流体热毛细流动和换热影响的实验研究［J］. 中国科学院大学学报， 2019， 36（1）： 25-30. DOI： 10.7523/j.issn.2095-6134.2019.01.005 .

[12] Wang H Y， Hu J S， Gao X， et al. Influence of Li and B coatings of metal walls on deuterium retention and plasma confinement in HT-7［J］. Nuclear Fusion， 2012， 52（10）： 103002. DOI： 10.1088/0029-5515/52/10/103002 .

[13] Zuo G Z， Hu J S， Zhen S， et al. Comparison of various wall conditionings on the reduction of H content and particle recycling in EAST［J］. Plasma Physics and Controlled Fusion， 2012， 54（1）： 015014. DOI： 10.1088/0741-3335/54/1/015014 .

[14] Kugel H W， Mansfield D， Maingi R， et al. Evaporated lithium surface coatings in NSTX［J］. Journal of Nuclear Materials， 2009， 390/391： 1000-1004. DOI： 10.1016/j.jnucmat.2009.01.262 .

[15] Bell M G， Kugel H W， Kaita R， et al. Plasma response to lithium-coated plasma-facing components in the National Spherical Torus Experiment［J］. Plasma Physics and Controlled Fusion， 2009， 51（12）： 124054. DOI： 10.1088/0741-3335/51/12/124054 .

[16] Rindt P， Korving S Q， Morgan T W， et al. Performance of liquid-lithium-filled 3D-printed tungsten divertor targets under deuterium loading with ELM-like pulses in Magnum-PSI［J］. Nuclear Fusion， 2021， 61（6）： 066026. DOI： 10.1088/1741-4326/abf854 .

[17] Yamamoto Y， Kunugi T. Direct numerical simulation of liquid metal free-surface turbulent flows imposed on wall-normal magnetic field［J］. Fusion Engineering and Design， 2018， 136： 925-930. DOI： 10.1016/j.fusengdes.2018.04.041 .

[18] Hoashi E， Yoshihashi-Suzuki S， Nanba H， et al. Numerical study on free surface flow of liquid metal lithium for IFMIF［J］. Fusion Engineering and Design， 2013， 88（9/10）： 2515-2519. DOI： 10.1016/j.fusengdes.2013.05.059 .

[19] Bühler L， Mistrangelo C， Najuch T. Magnetohydrodynamic flows in model porous structures［J］. Fusion Engineering and Design， 2015， 98： 1239-1243. DOI： 10.1016/j.fusengdes.2015.01.018 .

[20] Khodak A， Maingi R. Modeling of liquid lithium flow in porous plasma facing material［J］. Nuclear Materials and Energy， 2021， 26： 100935. DOI： 10.1016/j.nme.2021.100935 .

[21] Watts E T， Krim J， Widom A. Experimental observation of interfacial slippage at the boundary of molecularly thin films with gold substrates［J］. Physical Review B， 1990， 41（6）： 3466-3472. DOI： 10.1103/physrevb.41.3466 .

[22] Israelachvili J N. Measurement of the viscosity of liquids in very thin films［J］. Journal of Colloid and Interface Science， 1986， 110（1）： 263-271. DOI： 10.1016/0021-9797（86）90376-0 .

[23] Thompson P A， Robbins M O. Shear flow near solids： epitaxial order and flow boundary conditions［J］. Physical Review. A， Atomic， Molecular， and Optical Physics， 1990， 41（12）： 6830-6837. DOI： 10.1103/physreva.41.6830 .

[24] Zou C X， Sun X G， Xu C， et al. Wetting characteristics of lithium droplet on iron surfaces in atomic scale： a molecular dynamics simulation［J］. Computational Materials Science， 2018， 149： 435-441. DOI： 10.1016/j.commatsci.2018.03.058 .

[25] Sun X G， Xiao S F， Deng H Q， et al. Molecular dynamics simulation of wetting behaviors of Li on W surfaces［J］. Fusion Engineering and Design， 2017， 117： 188-193. DOI： 10.1016/j.fusengdes.2016.06.037 .

[26] Xu S， Fan X F， Gu C Z， et al. First principles and molecular dynamics study of Li wetting and diffusion on W surfaces［J］. Journal of Nuclear Materials， 2020， 539： 152345. DOI： 10.1016/j.jnucmat.2020.152345 .

[27] Liu S C， Yu X G. Molecular dynamics study on the slippage of liquid lithium flow in tungsten nanochannels［J］. Nuclear Fusion， 2023， 63（3）： 036007. DOI： 10.1088/1741-4326/acb27c .

[28] Plimpton S. Fast parallel algorithms for short-range molecular dynamics［J］. Journal of Computational Physics， 1995， 117（1）： 1-19. DOI： 10.1006/jcph.1995.1039 .

[29] Liu S C， Yu X G， Zhang N M. An embedded-atom method interatomic potential for lithium-iron binary system and its applications in the liquid first wall system［J］. Journal of Nuclear Materials， 2024， 589： 154867. DOI： 10.1016/j.jnucmat.2023.154867 .

[30] Ushigusa K， Seki M， Suganuma K， et al. Electron cyclotron resonance discharge cleaning by using LHRF system on JT-60U［J］. Fusion Engineering and Design， 1999， 45（2）： 137-144. DOI： 10.1016/S0920-3796（99）00002-2 .

[31] Ning R H， Li Y G， Zhou W H， et al. Modeling D retention in W under D ions and neutrons irradiation［J］. Journal of Nuclear Materials， 2012， 430（1/2/3）： 20-26. DOI： 10.1016/j.jnucmat.2012.06.029 .

基金资助

国家自然科学基金(12172356)

中国科学院青年交叉团队项目资助