bubble-based EMMS/PFB模型的建立及在加压流化床浓相段的应用

宋素芳，郝振华，董立波，李俊国，房倚天



bubble-based EMMS/PFB模型的建立及在加压流化床浓相段的应用

宋素芳1,2,3，郝振华1，董立波1，李俊国1，房倚天1

（1中国科学院山西煤炭化学研究所，煤转化国家重点实验室，山西太原030001；2太原科技大学化学与生物工程学院，山西太原 030021；3中国科学院大学，北京100049）

基于多尺度分解和能量消耗分析方法，结合压力下锥形分布板射流床气泡直径关联式，建立了一个适于加压流化床（PFB）的能量最小多尺度模型——bubble-based EMMS/PFB模型。应用此模型模拟一个二维加压射流床，分析了操作压力、位置高度、空隙率及剩余速度对非均匀因子的影响。通过模拟结果与实验数据的对比，发现该模型相比于Gidaspow模型，能够更准确地模拟加压射流床内颗粒浓度的分布状态及颗粒靠近壁面处的速度变化；将这种曳力模型应用到流化床浓相段的模拟，预测了床内颗粒浓度瞬时分布及沿轴向的时均值分布、颗粒的速度分布等流动行为，使流化床浓相段的气固流动行为可视化，对流化床的设计、放大有一定的指导作用。

流化床；多尺度；数值模拟；气泡；两相流；流体动力学

引 言

加压流化床（PFB）通过增加气体密度，强化气体和颗粒之间的相互作用而改善气固流动行为，能够明显提高气固接触效率[1]。相比于传统常压流化床，加压流化床具有几个明显的优点：气固接触和混合效果好[2]，传热效率高[3-4]，化学反应速率高[5]，设备体积小[6]。因此，加压流化床在工业上具有良好的应用前景，尤其是在流化床煤气化领域。

研究气固流动特征对于优化流化床操作和进行工业放大都具有极其重要的意义。实验研究方法在加压情况下费用较高且测量效果不佳[7]，而数值模拟方法弥补了这一不足。气固流动模拟应用比较广泛的是双流体模型，该模型动量方程中固相应力一般采用颗粒动理学理论（KTGF）进行封闭；动量方程中气固相间曳力经验公式多数是建立在均匀流动前提下的，对气泡或颗粒聚团等非均匀流动结构的影响并没有考虑。针对此问题而发展起来的EMMS（energy-minimization multi-scale）方法在流化床气固非均匀流动分析方面早有应用[8-10]。

EMMS模型最早由Li等[11]提出时是基于颗粒聚团，即以颗粒聚团直径预测为求解方程之一，主要适用于高气速的、团聚物占主导的循环流化床，后来大多数的EMMS模型模拟研究[12-21]也都是以这个模型为基础。但是鼓泡流化床、湍动流化床等具有以气泡为主要不均匀结构的气固流动特征，在模型计算时可以将气泡直径作为模型参数之一，代替颗粒聚团的尺寸方程进行计算，而且相比于颗粒聚团，气泡直径的预测更具有丰富的理论依据，其预测结果也更准确，Shi等[22]的EMMS/bubbling 模型就是以这个思想为出发点提出的。

相关研究发现，压力和气泡所处的位置高度会直接影响到流化床气泡尺寸的大小[23]。传统EMMS模型并没有体现这些因素的影响，本文针对加压流化床气固流动特征，以气泡直径随操作压力及气泡位置变化的拟合公式为基础建立多尺度模型，从而使之适用于加压流化床气固两相流动系统，称之为bubble-based EMMS/PFB模型。

1 加压流化床气泡直径关联式

1.1 多尺度模型中气泡直径关联式分析

气泡直径是流化床气固流动中非常重要的一个参数，很多和速率相关的参数如气泡上升速度、相间传质传热速率等都取决于气泡的大小。可靠的气泡直径关联式对于流化床设计及精确模拟有着至关重要的决定作用。很多学者对气泡直径关联式做了大量的研究[24-25]，这些研究多是把气泡直径和剩余气速(g-mf)关联。但是，大量的研究表明，气泡直径不仅和上述因素有关，还和流化床操作压力以及气泡所处的位置高度有关。除此之外，分布板的结构也会对气泡直径有一定的影响。

1.2 压力下气泡尺寸的关联式

Cai等[23]在总结了大量的文献研究的基础上，指出除了操作气速非常低的情况下，不管是鼓泡区还是湍流区，气泡尺寸随操作压力的升高而减小，并拟合得到了关于压力和位置的气泡尺寸的关联式，该关联式压力适用范围为0.1～7.0 MPa。在Cai等的加压流化床气泡直径关联式基础上，考虑锥形分布板高度的影响，对该关联式进行局部修正，得到气泡直径预测表达式如下

式中，是气泡所处的垂直位置距离分布板的高度；0是常数，为射流入口与锥形分布板上端的垂直距离；是操作压力；(g-mf)是剩余气速；其余常数=-1.4×10-4,=0.25,=0.1。

在操作压力为0.7 MPa，剩余气速分别为0.3 m·s-1和0.35 m·s-1的实验条件下，不同床层高度下通过直径预测关联式得到的计算值与实验值吻合良好[7, 26]，基本能够满足计算要求。

2 模型的建立

2.1 压力下流化床气固流动系统分解

基于EMMS模型的思想，可以将加压流化床气固两相流动系统分为密相、稀相及相间相3个子系统，对应动力学参数分别为密相中表观气体速度(ge)、密相中颗粒速度(pe)、密相空隙率(ge)、气泡的体积分数(b)、气泡直径(b)、气泡速度(b)及气泡的加速度（b）。

稀相（气泡相）中假设只含气体[22]，即gb=1，气泡的密度和黏度等于其中气体的密度和黏度；密相（乳化相）可以看成一个连续的、均匀的混合相，视为拟流体，其密度、黏度及表观气速等参数根据所包含气固两相流体参数进行计算[27-28]。

2.2 平衡方程

单位体积内密相颗粒受力平衡方程：密相内气固曳力与相颗粒所受的有效重力相等，由于乳化相气体和颗粒的惯性相差很大，乳化相中颗粒加速度可以忽略[28]，即e=0。

单位体积内相间相受力平衡方程：曳力和有效重力平衡。

(3)

单位体积内气相和颗粒相质量守恒方程：简化起见，全体平均颗粒速度可以近似等于0(p=0)。

g=ge(1-b)+bb(4)

p=pe(1-b) (5)

平均空隙率通过密相和稀相的空隙率进行计算

g=ge(1-b)+bb(6)

2.3 稳定性方程

稳定性模型采用Shi等[22]提出的方程

其中b表示气泡中的气体占总的气体的比率，表达式为[28]

(8)

通过上面的分析可以看出考虑气泡非均匀结构的曳力系数可以表示为

将此曳力系数与不考虑非均匀结构的Wen &Yu曳力系数相除，得到非均匀因子如下

(10)

其中

将上述气泡直径的拟合公式、平衡方程以及稳定性条件联合起来，利用全局寻优方法来求解子系统的流动参数[26,29]，计算非均匀因子，从而对曳力系数进行修正。

3 模型的检验

3.1 模拟参数设置

为了验证模型的适用性，以一个二维加压射流床实验为基础进行模拟[7]，射流床结构如图1所示，操作条件、物性参数及模拟参数设置见表1。由于实验本身是一个二维床，为计算方便，采用二维网格模拟，网格尺寸为2 mm和10 mm两种规格，均为结构化网格。

表1 模拟参数设置

采用双流体模型进行模拟，用颗粒动理学理论(KTGF)封闭动量方程中固相应力，曳力系数采用基于气泡非均匀结构的bubble-based EMMS/PFB模型。

3.2 非均匀因子讨论

3.2.1 非均匀因子 为研究操作条件对非均匀因子的影响，分别计算了不同压力和不同剩余气速条件下的非均匀因子，将非均匀因子拟合成空隙率和位置高度的关系式见表2，它们有着共同的函数形式。

非均匀影响因子求出以后，参照Gidaspow模型模型，曳力系数的表达式可以写出如下

其中，d为非均匀影响因子d和1的交叉点的空隙率。

表2 不同操作条件下的非均匀影响因子

Note: Ifd≥1，then letd=1。

3.2.2 位置高度的影响 以=0.7 MPa，fmf=0.35 m·s-1为例分析非均因子的特点。首先，从图2可以看出非均因子是随位置升高而减小，气泡在上升过程中不断发展、长大、成熟，位置升高，气泡尺寸增大，流动局部非均匀性更加明显，非均匀影响因子就越小，对曳力的影响越显著。

3.2.3 空隙率的影响 非均因子在所研究空隙率范围内是一条有最小值的开口向上的曲线，越靠近空隙率两端，非均匀影响因子越大，流动也就越接近均匀状态。从图2中可以看出，不同位置高度，非均匀因子最低值对应的空隙率也不同。总体上，空隙率0.47以上，随空隙率增加，非均匀因子变大，说明气固两相整体流动趋于相对均匀。

3.2.4 操作压力的影响 在剩余气速相等的情况下，任取一个床层高度=0.4 m，分别计算操作压力为0.7、1.1、2.1 MPa时不同空隙率对应的非均匀影响因子d的值，如图3所示。可以看出，随压力升高，非均匀影响因子d值变得越大，原因可能在于剩余气速相同情况下，随压力升高，气泡尺寸减小，气泡数目增多，气泡分布变得相对均匀，因此，非均匀程度降低，非均匀影响因子升高。

3.2.5 剩余气速的影响 在操作压力=0.7 MPa时，任取一个床层高=0.5 m，分别计算剩余气速f-mf为0.20、0.35、0.50 m·s-1时非均匀影响因子d值，从图4可以看出剩余气速对d的影响主要表现在剩余气速增加会加强非均因程度，使非均匀因子减小，这与Shi等[22]得出的结论是一致的。

3.3 验证结果

采用bubble-based EMMS/PFB模型来模拟上述二维加压射流床，并与Gidaspow模型模拟结果及实验结果对比。图5为网格尺寸分别为10 mm和2 mm时瞬时颗粒相浓度分布。可以看出bubble-based EMMS/PFB模型模拟结果更接近实验值，由于Gidaspow模型没有考虑局部非均匀结构的影响，曳力被高估，计算的床层膨胀率也偏高。同时，相比于Gidaspow模型，bubble-based EMMS/PFB模型在10 mm和2 mm网格下的模拟结果相差较小，说明该模型对网格精度要求相对较低。

流化床在分布板上方靠近壁面的地方颗粒通常运动速度较慢，对于煤气化或煤燃烧流化床来说，此处近氧气入口，温度较高，热量如果不能被及时移走，易形成过热点导致结渣，加压操作下煤气化或煤燃烧流化床尤其如此。因此，对于加压操作流化床分布板上方靠近壁面的地方颗粒运动情况应该尤其关注。采用2 mm网格进行模拟，对靠近壁面处(分布板以上150 mm，距离壁面为15 mm)颗粒的流动速度进行对比，bubble-based EMMS/PFB模型预测值较接近实验值，相对精确，如图6所示。

通过上述对Gidaspow模型和bubble-based EMMS/PFB模型的模拟结果对比分析可以看出，在加压操作下，采用bubble-based EMMS/PFB模型能够更准确地模拟流化床内气固两相流动时颗粒浓度的分布状态及颗粒靠近壁面处的速度变化。

4 模型在流化床浓相段的应用

4.1 模拟对象及参数设置

将bubble-based EMMS/PFB模型用来模拟加压流化床下部浓相段气固两相流动，为简化问题，加压流化床浓相段几何模型如图7所示，采用三维模拟。划分网格时，将几何体分为上部直筒段和下部锥形段，网格尺寸均为20 mm。上部直筒段采用三棱柱网格，网格数目为637364个；下部锥形段采用金字塔网格，网格尺寸均为102899个。

浓相段气固两相的物性参数、操作参数及模拟设置见表3，采用双流体模型，耦合bubble-based EMMS/PFB曳力模型，湍流模型为标准模型，气相在壁面处为无滑移边界，固体颗粒采用Johnson等[30]的部分滑移条件来进行计算气固流动模拟计算。

表3 流化床浓相段模拟参数设置

4.2 瞬时固体浓度分布

图8为20 s时分别用Gidaspow模型与bubble-based EMMS/PFB模型预测得到的在不同位置的瞬时固体浓度切片图。可以看出两种模拟得到相似的流动分布，均可以观察到明显的气泡，颗粒浓度在分布板上方靠近壁面处较高，在接近床层表面时，颗粒浓度迅速降低；在床层下部靠近中心处则由于中心射流的缘故，颗粒浓度较低。由于Gidaspow模型没有考虑非均匀结构导致的气体携带颗粒能力的降低这一因素，计算的床层高度偏高，而且颗粒浓度较低，尤其是在分布板上方靠近壁面处，两种模型计算的结果差距明显。

4.3 时均颗粒浓度沿轴向分布

图9为15～25 s时间段两种模型计算出的时均颗粒浓度沿轴向的分布。由图可见，两种模拟得到的轴向颗粒浓度分布均呈现“S”形，即在床层底部浓度较高，随着到达床层表面，浓度迅速下降，在上方的自由空间，浓度接近于0。相对于bubble-based EMMS/PFB模型，Gidaspow模型由于高估了曳力，得到了偏高的床层高度，时均颗粒浓度随高度的变化较缓慢。

4.4 颗粒速度分布

分别采用Gidaspow模型和bubble-based EMMS/PFB模型模拟得到40 s时平面上颗粒瞬时速度分布，如图10所示。两种模型得到的颗粒速度分布类似，由于气泡的作用，气固流动存在强烈的返混，颗粒速度大小、方向均不均匀。由于气泡对颗粒的携带作用，颗粒在靠近床层中心的位置具有较高的上升速度，而在靠近壁面的地方多数向下运动。对于流化床煤气化来说，向下的颗粒流动有利于灰渣的排出。由于Gidaspow模型高估了气固之间的曳力，整体颗粒速度较bubble-based EMMS/PFB模型计算结果偏高。

5 结 论

（1）根据加压流化床气固流动特点，结合加压流化床的气泡预测关联式，建立了压力下基于气泡不均匀结构的多尺度模型，将其命名为bubble-based EMMS/PFB模型。

（2）采用一个二维加压射流床实验来验证bubble-based EMMS/PFB模型模拟的适用性，得到了瞬时颗粒浓度分布与近壁颗粒时均速度，模拟结果与实验数据吻合良好。

（3）将bubble-based EMMS/PFB模型应用到流化床浓相段的模拟，预测了不同位置的颗粒浓度分布、截面颗粒平均浓度轴向分布及颗粒速度分布变化等流动行为特征，对流化床煤气化的操作及放大有一定的指导意义。

符 号 说 明

ab, ae——分别为气泡、乳化相颗粒加速度，m·s-2 Cde,Cdi——分别为乳化相、气泡相有效曳力系数 db,dp——分别为气泡、颗粒直径，m Fde,Fdi, Fdb, Fe——分别为单位体积乳化相、相间相、气泡、乳化相颗粒的曳力，kg·m-2·s-2 fb——气泡气体体积占总气体比例 g——重力加速度，m·s-2 H——分布板上方垂直距离，m Hd——非均匀因子 me, mi——分别为单位体积乳化相内颗粒数目、单位体积内气泡数目 Nst——悬浮输送能，J·s-1·kg-1 P——操作压力，MPa Ub, Ue——分别为气泡相、乳化相平均表观气速，m·s-1 Ug——表观气速，m·s-1 Ugb,Uge——分别为气泡相、乳化相表观气速，m·s-1 Up,Upe——分别为表观颗粒速度，乳化相颗粒速度，m·s-1 Uslip——表观滑移速度，m·s-1 Uslip,e,Uslip,i——分别为乳化相、相间相表观滑移速度，m·s-1 be, bw——分别为控制体积内基于气泡曳力系数、Wen &Yu 曳力系数，kg·m-3·s-1 db——气泡相体积分数 eg——空隙率 下角标 g——气相 s——固相

References

[1] LI J, KUIPERS J A M. Effect of pressure on gas-solid flow behavior in dense gas-fluidized beds:a discrete particle simulation study[J]Powder Technology, 2002, 127: 173-184.

[2] Olowson P A, Almstedt A E . Influence of pressure and fluidization velocity on the bubble behaviour and gas flow distribution in a fluidized bed[J]Chemical Engineering Science, 1990, 45(7): 1733-1741.

[3] Li Z, ANNALAND M V S, KUIPERS J A M,. Effect of operating pressure on particle temperature distribution in a fluidized bed with heat production[J]Chemical Engineering Science, 2016, 140: 279-290.

[4] Nyoman S, Winaya P B. Effect of pressure and carbon dioxide concentration on heat transfer at high temperature in a pressurized circulating fluidized bed (PCFB) combustor[J]International Journal of Heat and Mass Transfer, 2001, 44: 2965-2971.

[5] Cao J T, Cheng Z H, Fang Y T. Simulation and experimental studies on fluidization properties in a pressurized jetting fluidized bed[J]Powder Technology, 2008, 183(1): 127-132.

[6] Valin S, Ravel S, Guillaudeau J. Comprehensive study of the influence of total pressure on products yields in fluidized bed gasification of wood sawdust[J]Fuel Processing Technology, 2010, 91(10): 1222-1228.

[7] 李俊国. 压力下射流流化床气固流动特性的研究[D]北京: 中国科学院大学, 2013. LI J G. Hydrodynamic characteristics in jetting fluidized bed under pressure[D]. Beijing: University of Chinese Academy of Sciences, 2013.

[8] Wang J W, Zhou Q, Hong K. An EMMS-based multi-fluid model (EFM) for heterogeneous gas-solid riser flows(Ⅱ): An alternative formulation from dominant mechanisms[J]Chemical Engineering Science, 2012, 75: 349-358.

[9] Wang W, Li J H. Simulation of gas-solid two-phase flow by a multi-scale CFD approach of the EMMS model to the sub-grid level[J]Chemical Engineering Science, 2007, 62(1/2): 208-231.

[10] Zhou Q, Wang J W. Coarse grid simulation of heterogeneous gas-solid flow in a CFB riser with EMMS drag model: effect of inputting drag correlations[J]Powder Technology, 2014, 253: 486-495.

[11] LI J H, Kwauk M. Particle-Fluid Two Phase Flow: Energy-minimization Multi-scale Method[M]. Beijing: Metallurgical Industry Press, 1994.

[12] Yang N, Wang W, Ge W. Choosing structure-dependent drag coefficient in modeling gas-solid two-phase flow[J]China Particuology, 2003, 1(1): 38-41.

[13] Yang N, Wang W, Ge W. CFD simulation of concurrent-up gas-solid flow in circulating fluidized beds with structure-dependent drag coefficient[J]Chemical Engineering Journal, 2003, 96(1/2/3): 71-80.

[14] Zeneli M, Nikolopoulos A, Nikolopoulos N. Application of an advanced coupled EMMS-TFM model to a pilot scale CFB carbonator[J]Chemical Engineering Science, 2015, 138: 482-498.

[15] ZHAO M, ZHANG N, LIU X H. Numerical simulation of loop seal for a circulating fluidized bed based on an improved EMMS bubbling model[J]Powder Technology, 2015, 284: 443-451.

[16] ZHOU Q, WANG J W. CFD study of mixing and segregation in CFB risers: extension of EMMS drag model to binary gas-solid flow[J]Chemical Engineering Science, 2015, 122: 637-651.

[17] 陈程, 祁海鹰. EMMS曳力模型及其颗粒团模型的构建和检验[J]化工学报, 2014, 65(6): 2003-2012. CHEN C, QI H Y. Development and validation of cluster and EMMS drag model[J]CIESC Journal, 2014, 65(6): 2003-2012.

[18] 鲁波娜, 程从礼, 鲁维民. 基于多尺度模型的MIP提升管反应历程数值模拟[J]化工学报, 2013, 64(6): 1983-1992. LU B N, CHENG C L, LU W M. Numerical simulation of reaction process in MIP riser[J]. CIESC Journal, 2013, 64(6): 1983-1992.

[19] 王维, 洪坤, 鲁波娜. 流态化模拟:基于介尺度结构的多尺度CFD[J]化工学报, 2013, 64(1): 95-106. WANG W, HONG K, LU B N.Fluidized bed simulation: structure-dependent multiscale CFD[J]CIESC Journal, 2013, 64(1): 95-106.

[20] 王雪瑶, 姜凡, 徐祥. 基于EMMS方法的曳力修正模型及实验验证[C]// 中国工程热物理学会多项流2009年学术会议. 乌鲁木齐:中国工程热物理学会, 2009: 1-7. WANG X Y, JIANG F, XU X. Modified drag model and experimental verification based on the EMMS method[C]// Proceedings of the Chinese Academy of Engineering Thermophysics 2009. Urumchi: CSET, 2009: 1-7.

[21] 周泉, 王军武, 李静海, 基于EMMS的双流体模型及其在高密度循环流化床中的应用[C]//第七届全国流态化会议. 北京:流态化专业委员会, 2013: 83-85. ZHOU Q, WANG J W, LI J H.EMMS-based two fluid model and its application to high-density circulating fluidized bed[C]//Proceedings of the 7th National Conference on Fluidization. Beijing: Fluidization Specialized Committee,2013: 83-85.

[22] SHI Z S, WANG W, LI J H. A bubble-based EMMS model for gas-solid bubbling fluidization[J]Chemical Engineering Science, 2011, 66(22): 5541-5555.

[23] Cai P, Schiavetti M, Michele G D. Quantitative estimation of bubble size in PFBC[J]Powder Technology, 1994, 80: 99-109.

[24] Horio M, NONAKA A. A generalized bubble diameter correlation for gas-solid fluidized beds[J]AIChE Journal, 1987, 33(11): 1865-1872.

[25] Mori S, Wen C Y. Estimation of bubble diameter in gaseous fluidized beds[J]AIChE Journal, 1975, 21(1): 109-115.

[26] SONG S F, HAO Z H, DONG L B. A bubble-based EMMS model for pressurized fluidization and its validations with data of a jetting fluidized bed[J]RSC Advances, 2016, 6: 111041-111051.

[27] Thomas D G. Transport characteristics of suspension(Ⅷ): A note on the viscosity of Newtonian suspensions of uniform spherical particles[J]Journal of Collioid Science, 1965, 20: 267-277.

[28] 石战胜. 基于气泡的EMMS模型及其在煤气化炉CFD模拟中的应用[D]北京: 中国科学院研究生院, 2011. SHI Z S. A bubble-based EMMS model for gas-solid fluidized beds and its application in CFD simulation of a coal gasifier[D]. Beijing: University of Chinese Academy of Sciences, 2011.

[29] 刘哲语, 李俊国, 李春玉. 多段分级转化流化床煤气化技术研发与进展[J]煤化工, 2016, (1): 3-6. LIU Z Y, LI J G, LI C Y. Research and development of multi-stage conversion fluidized bed coal gasification technology[J]Coal Chemical Industry, 2016, (1): 3-6.

[30] Johnson P C, Jackson R. Frictional-collisional constitutive relations for granular materials, with application to plane shearing[J]Journal of Fluid Mechanics, 1987, 176: 67-93.

Establishment of bubble-based EMMS/PFB model and its application on dense section of pressurized fluidized-bed

SONG Sufang1,2,3, HAO Zhenhua1, DONG Libo1,LI Junguo1, FANG Yitian1

(1State Key Laboratory of Coal Conversion, Institute of Coal Chemistry, Chinese Academy of Sciences, Taiyuan 030001, Shanxi, China;2School of Chemical and Biological Engineering, Taiyuan University of Science and Technology, Taiyuan 030021,Shanxi,China;3University of Chinese Academy of Sciences, Beijing 100049, China)

An energy minimization and multi-scale model, bubble-based EMMS/PFB, was established for pressurized fluidized-bed, based on methods of multi-scale decomposition and energy consumption analysis in combination with bubble diameter correlation to conical distributor. The model was used to simulate a 2D pressurized jetting fluidized-bed and to study influence of operating parameters such as operating pressure, bubble location, porosity and residual velocity on heterogeneous index. Comparison of simulation results to experimental data, it was showed that the new model provided better accuracy than the Gidaspow model in prediction of solid particle concentration distribution and change of particle velocity adjacent to wall above the distributor. Dense section simulation of a pressurized fluidized-bed by the new drag model yielded prediction of instantaneous and time-averaged axial distributions of particle concentration and particle velocity distribution. The simulation could make it available to visualize gas-solid flow behavior at dense section and provide guidance for design and industrial scale-up of pressurized fluidized-beds.

fluidized-bed; multi-scale; numerical simulation; bubble; two-phase flow; fluid dynamics

10.11949/j.issn.0438-1157.20170230

O 359；TQ 545

A

0438—1157（2017）08—2998—08

郝振华。第一作者：宋素芳（1982—），女，博士研究生，讲师。

国家自然科学基金项目（21606250）；中国科学院青年创新促进会基金项目（2016162）。

2017-03-09收到初稿，2017-05-17收到修改稿。

2017-03-09.

HAO Zhenhua, hzh@sxicc.ac.cn

supported by the National Natural Science Foundation of China (21606250) and the Youth Innovation Promotion Association, CAS(2016162).