环丁砜萃取精馏过程模拟分析及工艺参数优化

汪勤，张冰剑，何畅，何昌春，陈清林



环丁砜萃取精馏过程模拟分析及工艺参数优化

汪勤，张冰剑，何畅，何昌春，陈清林

（中山大学化学工程与技术学院，广东省石化过程节能工程技术研究中心，广东广州 510275）

以环丁砜-烃类相平衡数据和NRTL-RK热力学方法为基础，对环丁砜萃取精馏过程进行了全流程模拟和工艺操作参数优化。综合考虑各个操作变量及其关联，提出了基于局部耦合参数迭代优化的整体协同优化的策略。通过文献数据回归和Aspen Plus物性估算系统相结合，补充修正了缺失的模型参数，并以此模拟分析了各关键操作参数对环丁砜萃取精馏过程能耗和分离效果的影响。结果表明：当萃取精馏塔操作压力为0.17 MPa，溶剂回收塔操作压力为0.05 MPa时，贫溶剂最佳温度为100℃，原料饱和气相进料的最佳进料位置为第50块塔板；溶剂回收塔最佳回流比为0.33；最佳进料位置为第6块塔板，汽提水量为2853 kg·h−1。优化后，装置最小热公用工程由1.158 GJ·t−1下降至0.802 GJ·t−1，节能效果显著。

萃取精馏；热力学；芳烃；计算机模拟；过程系统；优化

引 言

苯、甲苯、二甲苯（BTX）是重要的基础化工原料，通常由石脑油、裂解汽油及焦化石脑油通过加氢及连续重整得到。近年来通过轻质烃类芳构化及重芳烃轻质化生产BTX芳烃的技术也得到较快发展[1-3]。由于同碳数的芳烃和非芳烃的沸点接近，且能形成共沸物，难以通过普通精馏工艺实现高效分离。目前，从烃类混合物中分离芳烃主要有液-液抽提和抽提蒸馏两种工艺[4-6]。一般地，针对馏分较窄、芳烃含量较高的原料，宜采用抽提蒸馏工艺；而针对馏分较宽、芳烃含量低的原料宜采用液-液抽提工艺。用于萃取精馏回收芳烃的溶剂主要有环丁砜、-甲酰基吡咯烷酮､-甲酰基吗啉等。其中环丁砜具有溶解能力强、选择性好、分离容易、溶剂损耗少、无毒价廉等特点，在现有国内外萃取精馏装置中得到了广泛应用[7]。

由于装置设计和操作原因，芳烃萃取精馏装置普遍存在循环溶剂量大、装置低温热过剩、过程能耗偏高等问题。目前，国内针对环丁砜芳烃萃取精馏过程的工艺改进和节能降耗进行了广泛研究，如环丁砜-芳烃体系相平衡数据的测定[8-12]、芳烃回收过程的模拟与参数优化[13-18]、动态模拟在萃取精馏装置设计中的应用[19]等。Choi等[14]应用Aspen Hysys软件模拟了芳烃抽提过程，并基于模拟数据进行了操作参数优化，但由于缺少关键组分的物性数据，该模拟过程的实际应用受到了一定限制。霍月洋等[16]基于文献报道的环丁砜-芳烃抽提体系相平衡数据，采用UNIFAC-Dortmund基团贡献模型拟合了CH2-SULF、ACH-SULF、ACCH2-SULF 3对基团的交互作用参数，结合流程模拟软件Aspen Plus对其工艺过程进行模拟研究，模拟结果表明各塔操作参数与实际值基本吻合，但并没有涉及以降低能耗为目标的优化分析。李春利等[19]根据文献报道的环丁砜-芳烃抽提蒸馏体系相平衡数据，采用UNIFAC-Dortmund基团贡献模型估算了UNIQUAC方程中的二元交互作用参数，并将其应用于环丁砜萃取过程模拟，获得了与实际工况值较吻合的结果。以此，进一步考察了进料位置、回流比、溶剂温度、溶剂比等操作参数对分离效果及能耗的影响，获得了环丁砜萃取精馏单元过程的较优操作参数，但优化过程中将各单元过程独立考虑，全局性不足。

本研究以环丁砜-烃类的液液相平衡数据、汽液相平衡数据为基础，利用Aspen Plus对国内某石化企业环丁砜萃取精馏过程建立了全流程模拟模型。通过模拟各关键操作参数对环丁砜萃取精馏过程能耗和分离效果的影响，在保证分离要求的条件下，以能耗为目标，综合考虑整个装置各个操作参数的相互关联，提出了基于局部参数优化的全局协同优化策略，得到各参数的优化值。

1 环丁砜萃取精馏原理及工艺

1.1 环丁砜萃取精馏原理

环丁砜萃取精馏过程涉及的物系为高度非理想体系，表征其相平衡特征的汽液平衡常数K表示为[20-21]

式中,y和x分别表示组分在汽相和液相中的浓度；γ为组分在液相中的活度系数；ps为纯组分在温度为时的饱和蒸气压；Φs为纯组分在温度为，压力为ps时的逸度系数；为纯组分在温度为、压力为时的气相逸度系数；VL为纯组分的液相摩尔体积；、为系统压力、温度。对于芳烃（组分）和非芳烃（组分）分离，相对挥发度表示为

ɑ= K/K=(y/y)/(x/x) (2)

在汽液平衡中，随温度和压力变化不显著，在计算过程中可近似为常数。当被分离组分的相对挥发度接近于1或形成共沸物时，应用普通精馏并不可行。如芳烃萃取精馏装置进料，主要成分为C6～C8的非芳烃和芳烃，其中能形成66.68℃共沸物苯-正己烷（0.0721:0.9279，mol）、77.45℃共沸物苯-环己烷（0.5487:0.4513，mol）、110.57℃共沸物甲苯-乙基环（0.9006:0.0994，mol）等，同时还存在芳烃-非芳烃近沸点分离问题。此时可以通过加入第三种极性溶剂改变待分离组分间的相对挥发度[22-23]。

1.2 环丁砜萃取精馏工艺

以国内某石化企业环丁砜-水复合溶剂萃取精馏装置为例，其工艺流程如图1所示，主要包括萃取精馏塔T1、非芳烃精馏塔T2、溶剂回收塔T3及换热设备等。原料S1与贫溶剂S12在T1塔内逆向接触传质，进行萃取精馏。T1塔顶蒸出含少量溶剂的非芳烃蒸气S3进入T2塔精制非芳烃；溶解了芳烃的富溶剂被送入T3塔汽提回收溶剂并精制芳烃。T1、T3塔顶回流罐经油水分离后的水相经换热、减压后作为汽提水循环使用。

2 环丁砜萃取精馏工艺模拟

2.1 原料及操作参数

装置原料为重整石脑油F1及异构轻烃F2，其组成及流量见表1。装置主要参数及操作条件见表2。装置产品控制目标为芳烃产品中的非芳烃含量不大于1500 mg·kg-1（P, S8≤1500 mg·kg-1），溶剂含量不大于2 mg·kg-1（S, S8≤2 mg·kg-1）；非芳烃产品中的芳烃质量分数不大于1%（A, S6≤1%）。

表1 原料组成及其流量

表2 装置主要参数及操作条件

2.2 热力学模型及参数

流程模拟软件Aspen Plus提供了包括流体汽液平衡常数、液液平衡常数、焓、熵、密度和其他传递性能参数的多种热力学方法。针对环丁砜萃取精馏等非理想体系分离，热力学方法的选择较为关键。相关文献报道应用UNIQUAC及NRTL方程能获得较好的模拟结果[24-25]。本研究选用NRTL-RK方程，该方程以NRTL[26]方程计算式(1)中的液相活度系数，以RK[27]方程计算式(1)中的气相逸度系数。

NRTL方程：

其中

G=exp(-) (4)

=a+b/+eln+fT(5)

=c+d(-273.15K) (6)

=0 (7)

G=1 (8)

式中，a、b、c、d、e、f为与组成相关的常数，c、d是对称性常数；a、b、e、f为非对称性常数；G、G、τ是量纲1的二元交互作用参数。准确的模型参数能提高模型的准确性。对于环丁砜萃取精馏体系Aspen Plus内嵌数据库缺少部分关键组分二元交互作用参数。本研究针对数据库中缺失且能从相关文献查阅的相平衡数据的关键组分的参数进行拟合，见表3。对于数据库中缺失且无文献报道的体系则采用Aspen Plus的物性估算系统估算。

表3 文献拟合的NRTL模型二元交互作用参数

2.3 模拟结果

在Aspen Plus模拟系统中更新以上回归和估算的关键组分交互作用参数，并建立严格的模拟模型。流程中所有的塔设备采用RADFRAC模块，在T1塔内进行严格的汽液相平衡计算；在T2、T3塔内进行塔板上的汽液相平衡计算及冷凝器处的汽液水三相的相平衡计算[30]。模拟计算结果见表4，可以看出现场标定值与模拟值相对偏差在5%以内。表明所采用的模拟策略和物性计算方法能较准确地反映该分离操作过程，模拟方案可行，其结果可以作为后续分离及用能优化的基础数据。

表4 模拟结果与标定值对比

3 装置关键操作参数分析与优化

对于萃取精馏塔，降低塔的操作压力可减少再沸器负荷，但是过低的操作压力会使塔内形成两个液相，不利于操作。且T1塔顶汽相出料直接作为T2塔进料，T1塔的操作压力也受T2塔操作压力限制。另一方面，溶剂回收塔操作压力受限于塔顶循环水温度及环丁砜分解温度。控制整个装置的操作压力与现有操作压力相同，即T1：0.17 MPa，T3：0.05 MPa。萃取精馏塔顶汽相出料，贫溶剂S12必须在第1块塔板进料。综合考虑，装置的优化变量包括：（1）萃取精馏塔操作变量，包括原料进料温度及进料位置（S2、S2）；（2）溶剂回收塔操作参数，包括塔顶回流比（T3）、富溶剂进料位置（S5）；（3）溶剂操作参数，包括溶剂烃含量（CH,S12）及溶剂进料温度（S13）。操作变量包括连续操作变量（温度、压力、浓度）和非连续操作变量（进料位置），且包括两股循环物流，Aspen Plus自带的优化模块难以实现全局优化。综合考虑各个操作变量及其关联，提出了基于局部耦合参数迭代优化的整体协同优化的策略。在保证产品质量的前提下，通过单因子变化，调整关键操作参数，以装置能耗作为目标函数对相关影响因素进行综合分析，实现装置总体用能优化[29-30]。相关模拟优化过程如图2所示。

3.1 萃取精馏塔参数优化

以原料进料温度、进料位置作为变量分析，在满足分离要求的情况下，可得到如图3所示的装置能耗（R,T1+R,T3）与原料进料温度（S2）、进料位置（S2）的关系。由图3可知，当T1的操作压力为0.17 MPa，CH,S12=1500 mg·kg−1，S12=108℃时，提高原料进料温度，装置总再沸负荷减小，两者的变化与S2呈线性变化，S2的优化值应结合装置中可供回收的低能级的余热考虑。当S2≥113℃时，原料进料状态为部分汽相或汽相。当S2在95～115℃变化时，最佳原料进料位置为第39～43块塔板。

3.2 溶剂回收塔参数优化

溶剂回收塔T3的控制目标为塔顶芳烃产品中的溶剂含量（S,S8）不大于2 mg·kg−1。图4反映了当CH,S12=1500 mg·kg−1，S12=108℃时T3塔再沸器负荷随回流比及进料位置的变化关系。图5给出了当CH,S12=1500 mg·kg−1，S12=108℃时T3塔再汽提水量（S11）随回流比及进料位置的变化关系。图6为当CH,S12=1500 mg·kg−1，S12=108℃时芳烃产品中的溶剂含量（S,S8）随回流比及进料位置的变化关系。通过以上分析可知：减小T3塔顶回流比有利于降低塔底再沸器负荷，同时降低塔底汽提水的消耗。但是过小的回流比可能使得塔顶芳烃产品中的溶剂含量超标，当T3<0.3时，S,S8始终大于2 mg·kg−1。靠近塔顶的进料位置能降低塔再沸器负荷，同时减小汽提水量，但是进料位置应当和塔顶保持一定塔板数以回收溶剂。如当T3=0.4～0.7时，满足产品分离条件的最佳进料位置为第5块塔板，当T3=0.3时，满足产品分离条件的最佳进料位置为第6块塔板。

图6 芳烃产品的溶剂含量随回流比及进料位置的变化

Fig.6S, S8varying withT3andS5

3.3 溶剂操作参数优化

借助Aspen Plus，以循环贫溶剂烃含量（CH,S12）、萃取精馏塔进料温度（S12）为变量进行装置操作参数分析与优化。在满足产品分离要求的情况下，得到如图7所示的塔再沸器负荷（R,T1,R,T3）与贫溶剂烃含量、进料温度的关系。由图7可知，当T1的操作压力为0.17 MPa，S2=95℃，S2=39，S5=6，T3=0.3时，贫溶剂的烃含量越低、进料温度越高，T1、T3的总体能耗越小。但当S12大于某一值时，如CH,S12=1000 mg·kg−1，S12≥108℃时，T1塔难以达到其分离要求。T1再沸负荷随着S12的变化先减小再增大，T3塔的再沸负荷由于S12物流进料所携带的热量，导致T1操作温度升高，富溶剂进入T3塔温度升高，使T3塔再沸器负荷随S12的增加单调减小。溶剂的最佳操作温度可确定为T1再沸器负荷最小时对应的温度，通过优化原料进料温度实现减小T1、T3的再沸器负荷并有效回收装置低温热的双重目的。

图8为T3塔再沸器负荷与溶剂烃含量的关系。T3塔再沸器负荷随着溶剂烃含量的增加而减小，即循环贫溶剂中的烃含量越大，T3塔再沸器负荷越小。但是，从图8可以看出，过小的溶剂烃含量反而使装置总体能耗增加。如当S12=100℃，最佳贫溶剂烃含量为400 mg·kg−1；当S12=108℃，最佳贫溶剂烃含量为200 mg·kg−1。

3.4 优化结果

装置操作参数优化后的冷热物流数据汇总见表5，主要包括3股热物流、3股冷物流及3股再沸物流。规定最小传热温差Δmin=15℃，得到优化前后装置换热网络的总组合曲线GCC，如图9所示。优化后装置操作参数及产品质量汇总见表6。由表6可知，优化前后装置的非芳烃产品和芳烃产品均满足要求，通过装置操作参数与换热网络的集成优化，整个装置最小热公用工程由1.158 GJ·t−1下降到0.802 GJ·t−1，减小30.76%；最小冷公用工程由1.166 GJ·t−1下降到0.877 GJ·t−1，减小24.77%。

表5 装置冷热物流数据

表6 装置优化操作参数及产品质量

4 结 论

应用NRTL-RK方法对环丁砜萃取精馏装置进行了全流程模拟，通过文献数据回归和Aspen Plus物性估算系统相结合完善了模型参数，提高了模型准确度，结果表明模拟计算结果与实际操作参数吻合较好，NRTL-RK方法及回归和估算的模型参数能有效应用于环丁砜萃取精馏模拟计算。通过模拟各关键操作参数对环丁砜萃取精馏过程能耗和分离效果的影响，在保证分离要求的条件下，以能耗为目标，综合考虑整个装置各个操作参数的相互关联，对各操作参数进行优化；同时，基于局部参数优化进行全局协同优化。通过流程与操作参数优化，整个装置最小热公用工程由1.158 GJ·t−1下降到0.802 GJ·t−1，减小30.76%；最小冷公用工程由1.166 GJ·t−1下降到0.877 GJ·t−1，减小24.77%，节能效果显著。

References

[1] 范景新, 臧甲忠, 于海斌, 等. 重芳烃轻质化研究进展[J]. 工业催化, 2015, 23(9): 666-673. FAN J X, ZANG J Z, YU H B,. Research progress in conversion of heavy aromatics to light ones[J]. Industrial Catalysis, 2015, 23(9): 666-673.

[2] 段然, 巩雁军, 孔德嘉, 等. 轻烃芳构化催化剂的研究进展[J]. 石油学报(石油加工), 2013, 29(4): 726-737. DUAN R, GONG Y J, KONG D J,. Development of the catalyst for light paraffins aromatization[J]. Acta Petrolei Sinica(Petroleum Processing Section), 2013, 29(4): 726-737.

[3] 孔德金, 杨为民. 芳烃生产技术进展[J]. 化工进展, 2011, 30(1):16-25. KONG D J, YANG W M. Advance in technology for production of aromatic hydrocarbons[J]. Chemical Industry and Engineering Progress, 2011, 30(1): 16-25.

[4] 徐春明, 杨朝合, 林世雄. 石油炼制工程[M]. 北京: 石油工业出版社, 2009: 557-569. XU C M, YANG Z H, LIN S X. Petroleum Refinery Engineering[M]. Beijing: Petroleum Industry Press, 2009: 557-569.

[5] Díaz I, Palomar J, Rodríguez M,. Ionic liquids as entrainers for the separation of aromatic-aliphatic hydrocarbon mixtures by extractive distillation[J]. Chemical Engineering Research & Design, 2016, 115: 382-393.

[6] 刘璐, 全晓宇, 黄哲庆, 等. 芳烃抽提模拟与改造[J]. 化学工程, 2015, 43(1): 73-78. LU L, QUAN X Y, HUANG Z Q,. Simulation and reformation of aromatic extraction[J]. Chemical Engineering, 2015, 43(1): 73-78.

[7] Larriba M, Navarro P, Gonzalez-Miquel M,. Dicyanamide-based ionic liquids in the liquid-liquid extraction of aromatics from alkanes: experimental evaluation and computational predictions[J]. Chemical Engineering Research & Design, 2016, 109: 561-572.

[8] Wang Z R, Huang L, Xia S Q,. Isobaric (vapour + liquid) equilibria for sulfolane with toluene, ethylbenzene, and isopropylbenzene at 101.33 kPa[J]. Journal of Chemical Thermodynamics, 2011, 43(12): 1865-1869.

[9] SANTIAGO R S, AZNAR M. Liquid-liquid equilibria for quaternary mixtures of nonane + undecane + (benzene or toluene or-xylene) + sulfolane at 298.15 and 313.15 K[J]. Fluid Phase Equilibria, 2007, 253(2): 137-141.

[10] KO M, IM J, SUNG J Y,. Liquid-liquid equilibria for the binary systems of sulfolane with cycloalkanes[J]. Journal of Chemical & Engineering Data, 2007, 52(4): 636-638.

[11] Cele N P, Bahadur I, Redhi G G,. (Liquid + liquid) equilibria measurements for ternary systems (sulfolane + a carboxylic acid +-heptane) at= 303.15 K and at 0.1 MPa[J]. Journal of Chemical Thermodynamics, 2016, 96: 169-174.

[12] Wisniak J, Ortega J, Fernández L. A fresh look at the thermodynamic consistency of vapour-liquid equilibria data[J]. Journal of Chemical Thermodynamics, 2017, 105: 385-395.

[13] Lyu Z, Zhou T, Chen L,Simulation based ionic liquid screening for benzene-cyclohexane extractive separation[J]. Chemical Engineering Science, 2014, 113(13): 45-53.

[14] Choi Y J, Cho K W, Cho B W,. Optimization of the sulfolane extraction plant based on modeling and simulation[J]. Industrial & Engineering Chemistry Research, 2002, 41(22): 5504-5509.

[15] Lin W C. Studies on the extraction of aromatics from C9+oil[J]. Journal of Chemical Engineering of Japan, 2002, 35(12): 1257-1262.

[16] 霍月洋, 张树增. 环丁砜芳烃抽提蒸馏的计算机模拟[J]. 北京化工大学学报(自然科学版), 2012, 39(3): 19-23. HUO Y Y, ZHANG S Z. Computer simulation of a process for extractive distillation of sulfolane aromatics[J]. Journal of Beijing University of Chemical Technology (Natural Science Edition), 2012, 39(3): 19-23.

[17] 王凌燕, 孙晓岩, 李忠杰, 等. 基于回归二元交互参数的芳烃抽提过程模拟[J]. 化工学报, 2011, 62(12): 3452-3457. WANG L Y, SUN X Y, LI Z J,. Simulation of aromatic extraction process based on regressed binary interaction parameters[J]. CIESC Journal, 2011, 62(12): 3452-3457.

[18] QIN J W, YE Q, XIONG X J,. Control of benzene-cyclohexane separation systemextractive distillation using sulfolane as entrainer[J]. Industrial & Engineering Chemistry Research, 2013, 52(31): 10754-10766.

[19] 李春利, 伍武辉, 王志彦, 等. 环丁砜萃取精馏工艺模拟及优化[J]. 现代化工, 2014, 34(4): 147-153. LI C L, WU W H, WANG Z Y,. Simulation and optimization of extractive distillation process with sulfolane[J]. Modern Chemical Industry, 2014, 34(4): 147-153.

[20] 刘家祺. 传质分离过程[M]. 北京: 高等教育出版社, 2005: 8-21. LIU J Q. Mass Transfer and Separation Processes[M]. Beijing: Higher Education Press, 2005: 8-21.

[21] 徐承恩. 催化重整工艺与工程[M]. 北京: 中国石化出版社, 2009: 806. XU C E. Catalytic Reforming Process and Engineering[M]. Beijing: China Petrochemical Press, 2009: 806.

[22] Prausnitz J M, Anderson R. Thermodynamics of solvent selectivity in extractive distillation of hydrocarbons[J]. AIChE Journal, 1961, 7(1): 96-101.

[23] Aniya V, De D, Satyavathi B. A comprehensive approach towards dehydration of-butyl alcohol by extractive distillation: entrainer selection, thermodynamic modeling and process optimization[J]. Industrial & Engineering Chemistry Research, 2016, 55(25): 6982-6995.

[24] Mohsen-Nia M, Doulabi F S M. Liquid-liquid equilibria for mixtures of (ethylene carbonate + aromatic hydrocarbon + cyclohexane)[J]. Thermochimica Acta, 2006, 445(1): 82-85.

[25] Mesquita F M R, Pinheiro R S, Sant’Ana H B D,. Removal of aromatic hydrocarbons from hydrocarbon mixture using glycols at 303.15 K and 333.15 K and atmospheric pressure: experimental and calculated data by NRTL and UNIQUAC models[J]. Fluid Phase Equilibria, 2015, 387(15): 135-142.

[26] Renon H, Prausnitz J M. Local compositions in thermodynamic excess functions for liquid mixtures[J]. AIChE Journal, 1968, 14(1): 135-144.

[27] REDLICH O, KWONG J N S. On the thermodynamics of solutions (Ⅴ): An equation-of-state. fugacities of gaseous solutions[J]. Chemical Reviews, 1949, 44(1): 233-244.

[28] Im J, Lee H, Lee S,. Liquid-liquid equilibria for the binary systems of sulfolane with branched cycloalkanes[J]. Fluid Phase Equilibria, 2006, 246(246): 34-38.

[29] YU Y X, HE M Y, GAO G H,. Boiling points for five binary systems of sulfolane with aromatic hydrocarbons at 101.33 kPa[J]. Fluid Phase Equilibria, 2001, 190(1/2): 61-71.

[30] Hu S, Zhang B J, Hou X Q,. Design and simulation of an entrainer-enhanced ethyl acetate reactive distillation process[J]. Chemical Engineering & Processing2011, 50(11/12): 1252-1265.

Process simulation and optimization of sulfolane extractive distillation

WANG Qin, ZHANG Bingjian, HE Chang, HE Changchun, CHEN Qinglin

(School of Chemical Engineering and Technology, Guangdong Engineering Center for Petrochemical Energy Conservation, Sun Yat-Sen University, Guangzhou 510275, Guangdong, China)

Full process simulation and critical parameter optimization were performed on sulfolane extractive distillation (ED), based on sulfolane/arene phase equilibrium data and NRTL-RK thermodynamic method. With consideration of multiple variables and their interactions, a coordinative optimization strategy was further proposed from iterative optimization of local coupling parameters. In order to improve modeling accuracy, missing parameters in NRTL-RK database were augmented through literature data regression and Aspen Plus property estimation. At given separation specifications, energy consumption and separation efficiency were simulated through adjusting critical operating parameters. The results show that when extractive distillation column (EDC) and solvent recovery column (SRC) were operated at pressure of 0.17 MPa and 0.05 MPa, respectively, the optimum EDC operating conditions were lean solvent temperature at 100℃ and feed location of saturated crude vapor at 50thstage whereas the optimum SRC operating conditions were reflux ratio at 0.33, feed location at 6thstage, and stripping water load at 2853 kg·h−1. After optimization, energy saving is significant with minimum heat consumption reduced from 1.158 GJ·t−1to 0.802 GJ·t−1.

extractive distillation; thermodynamics; arene; computer simulation; process systems; optimization

10.11949/j.issn.0438-1157.20161593

TQ 028.3

A

0438—1157（2017）05—1969—08

陈清林。

汪勤（1990—），女，博士研究生。

国家自然科学基金项目（21276288，U1462113）。

2016-11-10收到初稿，2017-01-10收到修改稿。

2016-11-10.

Prof. CHEN Qinglin, chqlin@mail.sysu.edu.cn

supported by the National Natural Science Foundation of China (21276288, U1462113).