超高纯氙去除氪低温精馏塔HYSYS模拟优化

doi:10.16183/j.cnki.jsjtu.2020.223

超高纯氙去除氪低温精馏塔HYSYS模拟优化

沙海东¹, 黄沛尧³, 王舟⁴, 崔祥仪⁵, 巨永林^,², 严锐², 李帅杰⁴, 王秀丽²

1.上海电力大学能源与机械工程学院, 上海 200090

2.上海交通大学制冷与低温工程研究所,上海 200240

3.上海交通大学巴黎高科卓越工程师学院, 上海 200240

4.上海交通大学粒子与核物理研究所, 上海 200240

5.上海交通大学李政道研究所, 上海 200240

Simulation and Optimization of Ultra-High Purity Krypton/Xenon Cryogenic Distillation System Based on HYSYS

SHA Haidong¹, HUANG Peiyao³, WANG Zhou⁴, CUI Xiangyi⁵, JU Yonglin^,², YAN Rui², LI Shuaijie⁴, WANG Xiuli²

1. College of Energy and Mechanical Engineering, Shanghai University of Electric Power, Shanghai 200090, China

2. Institute of Refrigeration and Cryogenics, Shanghai Jiao Tong University, Shanghai 200240, China

3. Paris Tech Elite Institute of Technology, Shanghai Jiao Tong University, Shanghai 200240, China

4. Institute of Particle and Nuclear Physics, Shanghai Jiao Tong University, Shanghai 200240, China

5. Tsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai 200240, China

通讯作者: 巨永林,男,教授,博士生导师,电话(Tel.):021-34206532;E-mail:yju@sjtu.edu.cn.

责任编辑: 陈晓燕

收稿日期: 2020-07-14

基金资助:

国家重点研发计划(2016YFA0400301)
上海市科学技术委员会重点基础研究项目(18JC1410200)

Received: 2020-07-14

作者简介 About authors

沙海东(1991-),男,河北省唐山市人,硕士生,研究方向为低温精馏与低温传热。

摘要

为满足大型暗物质探测器PandaX-4T低本底及高灵敏度的要求,探测介质氙(Xe)中氪(Kr)含量要求为n_Kr/n_Xe≤1×10^-13(n为物质的量).本文以超高纯氙去除氪低温精馏塔为例,首先简单介绍精馏塔及其结构设计.其次,对该塔的精馏塔参数进行HYSYS模拟及优化,得出现有运行工况下精馏塔内压力、温度以及氪浓度的分布,并研究了进料位置、回流比、进料流量、再沸器加热量、废品氙流量及进料压力等操作因素对精馏纯度的影响.提出优化提纯效果的精馏操作参数,即塔压为221 kPa、回流比为145时,原料氙中氪含量可从5×10^-7提纯到4.1×10^-14,这对超高纯氙去除氪低温精馏的产品纯度进一步提高有重要指导意义.

关键词： 低温精馏; 暗物质探测; HYSYS模拟

Abstract

To meet the requirements of the low background and high sensitivity of the large-scale dark-matter detector PandaX-4T, n_Kr/n_Xe≤1×10 ^-13 (n is the amount of substance) is required. The ultra-high purity krypton/xenon cryogenic distillation tower was used as an example in this paper. First, the distillation tower and the construction design of the distillation system were briefly introduced. Then, the operation parameters of the distillation tower were simulated and optimized by using HYSYS, and the distributions of the pressure, temperature and Kr concentration in the distillation tower were obtained. After that, the influences of those operational conditions including the feeding point position, the reflux ratio, the feeding flowrate, the heat of the reboiler, the flowrate of the off-gas, and the feeding pressure to the xenon purity of the distillation system were studied, and the optimized operation parameters were proposed. When the column pressure is set at 221 kPa and the reflux radio is set at 145, the krypton content in xenon can be purified from 5×10^-7 to 4.1×10^-14, which has an important guiding significance for further improvement of the purity of the xenon product of the xenon/krypton cryogenic distillation system.

Keywords： cryogenic distillation; dark matter detection; HYSYS simulation

PDF (1065KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

沙海东, 黄沛尧, 王舟, 崔祥仪, 巨永林, 严锐, 李帅杰, 王秀丽. 超高纯氙去除氪低温精馏塔HYSYS模拟优化[J]. 上海交通大学学报（自然版）, 2021, 55(7): 842-849 doi:10.16183/j.cnki.jsjtu.2020.223

SHA Haidong, HUANG Peiyao, WANG Zhou, CUI Xiangyi, JU Yonglin, YAN Rui, LI Shuaijie, WANG Xiuli. Simulation and Optimization of Ultra-High Purity Krypton/Xenon Cryogenic Distillation System Based on HYSYS[J]. Journal of shanghai Jiaotong University, 2021, 55(7): 842-849 doi:10.16183/j.cnki.jsjtu.2020.223

近百年来,天文观测从宇宙各个时间和距离尺度发现了很多证据,表明有一种神秘的、不发光的、全新未知物质充斥在宇宙里,质量是普通物质的5倍多.这类物质不参与电磁和强相互作用,被统称为暗物质.但是至今暗物质还未被实验探测发现,成为现代物理学中最难解释的谜团之一.

氙具有高原子序数(Z=54),且液氙的密度较高,能够阻止外界大部分γ射线,产生“自屏蔽”效应,同时液氙还具有高发光强度、低能阈及高能分辨率等优点,即氙本身具有较低的本底,故氙非常适合用作暗物质探测器的探测介质.PandaX实验是中国首个大型液氙暗物质直接探测实验^[1,2],探测灵敏度处于国际领先地位,并分别于2016年、2017年取得了当时国际对暗物质截面最领先的限制结果^[3,4].由于暗物质与普通物质相互作用非常微弱,需要不断的提高暗物质探测器的灵敏度.近年国际暗物质探测器规模不断升级,PandaX暗物质探测器中液氙的质量也由500 kg升级为4 t(PandaX-4T).

暗物质直接探测实验对探测介质的本底要求非常高,即氙中的杂质需要非常低以免干扰探测结果.其中氪-85(⁸⁵Kr)对探测结果影响尤其大,⁸⁵Kr属于放射性同位素,会产生β射线,干扰暗物质信号.当液氙中的n_Kr/n_Xe低于1×10^-12时才能保证暗物质探测器的灵敏度.而国内商业氙中氪的含量约5×10^-7,并不能满足暗物质探测实验要求,因此需要对商业氙进行精馏提纯,以降低其中的氪含量.

超高纯氪氙低温精馏系统^[5]于2019年在中国锦屏地下实验室完成了搭建工作,并于2019年1月完成调试运行,后经1.5个月稳定运行,共提纯5.75 t氙.该精馏塔可以99%的回收率,10 kg/h的速率将氙中氪含量从约5×10^-7降低到约1×10^-14.

1 超高纯氪氙精馏塔的结构与操作流程

1.1 精馏系统

超高纯氙除氪低温精馏^[6]是根据氪氙混合物沸点和挥发性的差异来实现的.在1 标准大气压下,氙(难挥发组分)的沸点是165 K,氪(易挥发组分)的沸点是120 K.氪氙低温精馏系统工艺流程如图1所示^[7],

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 氪氙低温精馏系统工艺流程图

1—原料氙,2—过滤器,3—气体纯化器,4—热交换器,5—PT60制冷机,6—精馏塔,7—AL300制冷机,8—冷凝器,9—液位计,10—再沸器,11—加热器,12—产品氙,13—真空罩,14—废品氙,15—真空抽口

Fig.1 Flowchart of Xe/Kr cryogenic distillation system

由精馏塔^[8]塔柱6、塔顶冷凝器8和塔底再沸器10构成,塔体外设置真空罩14防止热对流和热传导^[9].图中FM1、FM2及FM3为流量控制器,阀CV1、CV2为低温阀,阀SV1、SV2分别为内塔与外真空罩的安全阀.原料氙气由进口1进入,经过过滤器2去除其固体杂质,经气体纯化器3除去其中的水蒸气、氧气等气体杂质后进入换热器4,并由制冷机5预冷,随后从塔上部进入.氙气上升经过填料,与从塔顶冷凝器8中回流的液态氙以气液两相^[10]的方式进行传热传质,进行气液平衡交换,上升至塔顶冷凝器时,氙气成为含氪量较高的饱和气,称为废品氙,废品氙由塔顶引出后,进入废品氙储罐.在冷凝器中被制冷机液化的液氙回流后重新进入塔,在填料中与上升的氙气作用,两相进行传热传质和气液交换,使得液氙纯度进一步提升并流向塔底再沸器10,再沸器中的液氙一部分经加热器加热后汽化,重新进入塔内循环^[11,12],另一部分从底部被引出,作为产品氙进入产品氙储罐.

1.2 精馏塔运行操作流程

精馏塔运行分为5个阶段:预冷阶段→进料阶段→全回流阶段→提纯阶段→停机回收阶段^[9,13].预冷阶段是指在精馏塔正式运行前,向塔内充入一定量的原料氙,同时开启预冷制冷机和冷凝器上的制冷机,通过温度控制器控制制冷机的温度为178 K,为保证制冷和液化过程中精馏塔内压力维持在预设压力,需继续充入原料氙,直至塔底的再沸器开始出现液氙.此时进入进料阶段,在冷凝器制冷机维持178 K的情况下,向精馏塔继续补充原料氙,使再沸器内液位达到15 cm.随后进入全回流阶段,全回流阶段停止氙的进出,设定塔底再沸器的加热量(119 W),使精馏塔精馏段的部分液相氙汽化,塔顶冷凝器制冷机的温度仍为178 K,使得塔顶的氙气冷凝回流,塔内的氙气最终会达到比较稳定的气液平衡状态,纯度高的产品氙在再沸器中,纯度低的废品氙在冷凝器中.之后进入提纯阶段,以10 kg/h的速率提纯,设定原料氙、产品氙与废品氙的流量比值为100:99:1,这是个动态提纯过程,从塔顶和塔底不断流出废品氙和采出产品氙.最后是停机回收阶段,关闭制冷机和加热装置后将塔内剩余的气氙与液氙分别回收到液氮冷却的废品氙不锈钢气瓶和产品氙不锈钢气瓶中.

2 HYSYS稳态模拟计算及分析

2.1 数学模型

本文模拟的精馏塔塔径是125 mm,一般情况下,可假设认为塔板上液体与塔板间气体完全混合,即其具有均匀的组成、温度和压力,且精馏塔的设计计算是运用典型的M-T方法,如图2所示,图中Q_C为制冷机制冷量;Q_R为再沸器加热量.可假设每一块塔板汽、液两相接触后迅速达到平衡状态,离开塔板的汽液混合物也处于相平衡.故精馏塔中的每一块塔板都可以建立下列方程组^[14]:

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 精馏塔模型^[14]

Fig.2 Model of distillation tower

(1) 组分物料衡算方程

(1)s_ij+l_ij+g_ij+v_ij-v_i_,_j₊₁-l_i_,_j_-1-f_ij=0
(i=1,2,…,d;j=1,2,…,N)

式中:s_ij为液相侧线抽出S_j的摩尔流率;l_ij为组分i在塔中液流L_j的摩尔流率;g_ij为气相侧线抽出G_j的摩尔流率;v_ij为组分i在精馏塔中气流V_j的摩尔流率;f_ij为组分i在进料F_j中的摩尔流率;下角标j为塔板号,i为组分号.

(2) 组分相平衡方程

(2)v_ij=

\frac{K_{ij} V_{j}}{L_{j}}

l_ij

式中:K_ij为组分i在j板的平衡常数.

(3) 摩尔分率加和归一方程

(3)

\sum_{i} \frac{l_{ij}}{L_{j}}

\sum_{i} \frac{v_{ij}}{V_{j}}

\sum_{i} \frac{f_{ij}}{F_{j}}

(4) 热量衡算方程

(4) (S_j+L_j)h_j+(G_j+V_j)H_j-V_j₊₁H_j₊₁- L_j_-1h_j_-1-F_jh_Fj-Q_j=0

式中:H和h为汽、液相摩尔焓;h_F为加料的摩尔焓;Q_j为加入j板的热量.当j=1时,表示冷凝器,l_i₀=L₀=0;当j=N时,表示再沸器,V_N₊₁=v_i_,_N₊₁=0.

在计算过程中,除了上述方程组外,还需要添加K_ij、h_j及H_j的关系式,分别为

(5)K_ij=K_ij(p_j,T_j,v_ij,l_ij)

(6)h_j=h_j(p_j,T_j,l_ij)

(7)H_j=H_j(p_j,T_j,v_ij)

式中:T_j为第j块塔板处温度.

还需以下计算式:

(8)g_ij=

\frac{G_{j}}{V_{j}}

v_ij (i=1,2,…,d;j=2,3,…,N)

(9)s_ij=

\frac{S_{j}}{L_{j}}

l_ij (i=1,2,…,d;j=1,2,…,N-1)

因上述方程组与K_ij、h_j及H_j全为非线性关系,需要运用迭代法进行求解.Aspen Plus化工模拟软件通用于大型化工系统流程模拟,包含了近万种纯组分的物理数据以及多种专业的单元操作模块和特性方法.本文使用SRK特性方法进行模拟.该方法对所有热力学特性使用SRK状态方程^[15].该方程基于标准RKS状态方程和标准α函数:

(10)p=

\frac{RT}{V_{m} - b}

\frac{a}{V_{m} (V_{m} + b)}

(11)b=

\sum_{i}

x_ib_i

(12)b_i=

\frac{0.08664 R T_{c i}}{p_{ci}}

(13)a=a₀+a₁

(14)a₀=

\sum_{i = 1} \sum_{j = 1}

x_ix_j(a_ia_j)⁰^.⁵(1-k_ij)

(15)k_ij=

k_{ij}^{(1)}

k_{ij}^{(2)}

T+

\frac{k_{ij}^{(3)}}{T}

, k_ij=k_ji

(16)a₁=

\sum_{i = 1}

x_i

{\{\sum_{j = 1} x_{j} [(a_{i} a_{j})^{0.5} h_{ij}]^{1 / 3}\}}^{3}

(17)r_ij=

r_{ij}^{(1)}

r_{ij}^{(2)}

T+

\frac{r_{ij}^{(3)}}{T}

, r_ij≠r_ji

(18)a_i=0.427 47

\frac{R^{2} T_{ci}^{2}}{p_{ci}}

α_i(T)

(19)

\begin{array}{l} α_{i} (T) = {[1 + m_{i} (1 - T_{ri}^{1 / 2})]}^{2} \\ α_{i} (T_{ci}) = 1.0 \end{array}\}

(20)m_i=0.48508+1.5517ω_i-0.15613ω₂_i

式中:p为压力;R为气体常数;T为温度;V_m为原始SRK状态方程估算的摩尔体积;a₀为标准二次混合项;a₁为一个附加的不对称极性项;x_i为组分i的摩尔分数;a_i、a_j为计算过程中的求解变量;T_ci为组分i的临界温度;p_ci为组分i的临界压力;k_ij、r_ij分别为组分i、j的二元相互作用分数; $k_{ij}^{(1)}$ 、 $k_{ij}^{(2)}$ 、 $k_{ij}^{(3)}$ 及 $r_{ij}^{(1)}$ 、 $r_{ij}^{(2)}$ 、 $r_{ij}^{(3)}$ 是由数据回归获得的参数;T_ri为对应温度比,T_ri=T/T_ci;ω_i为偏心因子.

通过SRK状态方程,可以准确地计算出工质的物性,而精馏塔内的温度和压力的变化与精馏塔内氪氙组分的传热机理及理论塔板数的确定方法密切相关.汽液两相混合物平衡模型的温度、压力及浓度的关系为^[16]

(21)-X(ΔH/T)dT-XVdp+RTdX=0

式中:X为混合物浓度;ΔH为1 mol混合物在温度下从气态变为液态时的焓变;V为混合物体积.

当压力一定时,即在式(21)中dp=0时,则有:

(22)

\frac{dln X}{d T}

\frac{Δ H}{R T^{2}}

进一步求解式(22),可得混合物温度和浓度间的关系式为

(23)ln X=-

\frac{Δ H}{R T^{2}}

ΔT

当浓度一定时,即在式(21)中dX=0时,则有:

(24)

\frac{d p}{d T}

\frac{p Δ H}{R T^{2}}

对式(24)进行积分,可得混合物温度和压力的关系式为

(25)ln

\frac{p_{2}}{p_{1}}

\frac{Δ H}{R} (\frac{1}{T_{2}} - \frac{1}{T_{1}})

基于以上数学模型,进行边界条件设置及模拟分析^[17],与实验结果对比后优化了运行参数.

2.2 边界条件

根据设计,总塔板层数为17,从冷凝器记作第0块塔板,再沸器为第18块塔板,原料氙气进气位置为自上至下1/3处,即第5层塔板处进料.根据运行参数,选取冷凝器处压强为220.0 kPa,再沸器处压强为222.0 kPa,原料为气相,进气压强为221.0 kPa,进料流量设置为10 kg/h.原料氙为两元混合物,其中Kr的物质的量浓度为5×10^-7,Xe的物质的量浓度则为 0.999 999 5.最后设置其他参数:在参数设置窗口中设置回流比R'为145,回流液体流量为14.5 kg/h,产品氙流量为9.9 kg/h,即进料流量的99%,废品氙流量为0.1 kg/h,设置再沸器中的加热量为119 W,即428.4 kJ/h,精馏塔模型流程图如图3所示.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 精馏塔的模型流程图

Fig.3 Flowchart of model of distillation tower

2.3 结果与分析

应用HYSYS软件对实验操作参数进行模拟,并与实验测试结果^[18]对比,结果列于表1.

表1 实验结果与模拟结果对比

Tab.1 Comparisons of test and simulation results

项目	塔板数	理论塔板高度/cm	产品氙中氪的摩尔分数	废品氙中氪的摩尔分数	再沸器负荷/W	再沸器温度/K	塔内压力/ kPa	进料流量/ (kg·h^-1)	塔顶流量/ (kg·h^-1)
实验结果	17	35	<7.8×10^-12	1×10^-6	118	179.8	221	10	0.1
模拟结果	17	35	4.1×10^-14	8×10^-5	119	179.7	221	10	0.1

新窗口打开| 下载CSV

由表1可以看出,模拟结果与实验结果相差不大,基本上可以真实反映实验情况,模拟方法可靠.

本实验测量系统采用基于残余气体分析仪RGA对氪氙组分进行浓度测量的方法.目前RGA测量下限(1×10^-13)的限制导致测量结果可能无法准确反映产品氙的纯度,推测实际产品氙中的氪含量低于1×10^-13.

2.3.1 压力和温度分布精馏塔中压力和温度随塔板数n'的计算结果如图4所示.精馏塔中压力p_t和温度T呈梯度分布,随着塔板数的增大,压力由冷凝器的220 kPa线性增加至再沸器处的222 kPa,而温度也逐渐小幅升高,由冷凝器的179.48 K增加至再沸器处的179.67 K.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 精馏塔内温度与压力分布

Fig.4 Temperature and pressure distribution in distillation tower

2.3.2 进料位置对精馏纯度的影响保持加热量和回流比(数值为145)不变,中部进料(第8层塔板)和上部进料(第5层塔板)时各层塔板的氪气相摩尔分数的变化如图5所示.上部进料时氪气摩尔分数从第6块塔板开始明显低于中部进料,并且气相氪的摩尔分数小一个数量级甚至更多,同时上部进料在再沸器处(塔板18)得到的氪含量更低,产品氙更纯.从第5层板进料和第8层板进料时氪含量随塔板数的变化对比可以看出,第5层塔板进料为最优选,氪含量在塔底可达到约1×10^-14.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 上部进料和中部进料的各塔板气相氪摩尔分数对比

Fig.5 Comparisons of Kr concentration for the up and middle feeding in corresponding tray

2.3.3 回流比对精馏纯度的影响回流比R'为100、123、145、167及190时各塔板气相氪摩尔分数的变化如图6所示.回流比由100逐渐增大到190的过程中,再沸器(第18块塔板)处气相Kr的摩尔分数逐渐减小,说明当其它条件不变时,回流比越大,产品氙纯度越高,提纯效果更好.回流比为145时,氪含量可降低到约1×10^-14,回流比为190时氪含量可降低到约1×10^-16.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 不同回流比对应的各塔板气相Kr摩尔分数

Fig.6 Kr mole fraction in gas phase of each tray corresponding to different reflux ratios

2.3.4 进料量和再沸器加热量对精馏纯度的影响当加热量和回流比不变,进料流量q为1~14 kg/h,废品氙流量为0.1~1.4 kg/h(即进料流量的1%)时,再沸器处气相氪的摩尔分数如图7所示.随着进料流量由1 kg/h增大到14 kg/h,再沸器处气相氪的含量从2.35×10^-20上升到4.7×10^-8,增大幅度逐渐减弱.

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 进料流量和再沸器处气相氪摩尔分数的关系

Fig.7 Feeding flowrate versus Kr concentration in reboiler

再沸器加热量为119 W,进料流量q分别为6、8、10及14 kg/h时,回流比的变化如图8所示.进料流量为10 kg/h和废品氙流量为0.1 kg/h,再沸器加热量Q分别为100、110、119、130及140 W(即360、376、428.4、468 及504 kJ/h),R'的变化如图9所示.

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 加热量一定时进料流量和回流比的关系

Fig.8 Feeding flowrate versus reflux ratio at a constant heating amount

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 进料流量一定时再沸器加热量和回流比的关系

Fig.9 Heating capacity of reboiler versus reflux ratio at a constant feeding flowrate

由图8、9可知,加热量一定时,进料流量越大,回流比越小,产品氙中氪含量越大,产品氙纯度越差;而在进料流量一定时,再沸器加热量越大,回流比越大,产品氙中氪含量越小,产品氙纯度越好.

2.3.5 废品氙流量、进料压力和冷凝器压力加热量对精馏纯度的影响图10为进料流量为10 kg/h,回流比为145时,废品氙流量q、进料压力p₁和冷凝器压力p₂对产品氙中氪含量的影响.产品氙中氪的含量随废品氙的流量增大有明显的下降,在0.1~0.2 kg/h处下降明显,氪含量从约1×10^-14降低到约1×10^-18,之后随流量的增大下降趋于平缓,在废品氙流量达到1 kg/h时,氪含量已下降到10^-20量级以下,增大废品氙流量可以有效的提高产品纯度.进料压力和冷凝器压力均选取200~300 kPa进行模拟,产品氙中氪含量随进料压力和冷凝器压力的升高均增大,进料压力的升高使得产品氙中氪含量从 2.52×10^-14增大到2.97×10^-14,冷凝器压力的升高使得产品氙中氪含量从2.03×10^-14增大到 6.57×10^-14,但对其影响并不大且均在实验要求范围内.

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 废品氙流量、进料压力及冷凝器压力与产品氙中氪摩尔分数的关系

Fig.10 Kr concentration versus flowrate of off-gas Xe, feeding pressure, and pressure in flowrate condenser

3 结论

本文首先简单介绍了超高纯氪氙精馏塔的结构设计及操作流程,结合HYSYS模拟的方法对PandaX-4T暗物质探测实验的精馏塔进行了模拟计算和优化,其中包括数学模型、实验结果对比和进料位置、回流比、进料流量、再沸器加热量、废品氙流量、进料压力等操作因素对精馏纯度的影响.结果表明:再沸器加热量为119 W、塔压为221 kPa、回流比为145时可达到最优化操作条件.从上部(第5层塔板)进料可得到最优提纯效果.在实验条件允许的范围内,增大回流比明显有助于提纯效果,在确保进料流量即提纯速率为10 kg/h的情况下,建议用增大再沸器加热量的方法来增大回流比,以提高产品氙的纯度.当再沸器的功率受限于实验条件时,可考虑通过适当减小进料流量的方法提高产品氙纯度.与此同时,通过降低冷凝器压力、进料压力以及提高废品氙流量的方法也可以起到很好的提纯效果.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

季向东, 刘江来, 谌勋.

PandaX实验介绍和进展

[J]. 科学通报, 2016, 61(20): 2264-2272.