一欠平衡气侵与重力置换气侵特征及判定方法
2015-03-24张兴全周英操郭庆丰
张兴全, 周英操, 刘 伟, 郭庆丰
(1.中国石油集团钻井工程技术研究院,北京 102206; 2.中国石油勘探开发研究院,北京 100083)
一欠平衡气侵与重力置换气侵特征及判定方法
张兴全1,2, 周英操1, 刘 伟1, 郭庆丰1
(1.中国石油集团钻井工程技术研究院,北京 102206; 2.中国石油勘探开发研究院,北京 100083)
在裂缝性碳酸盐岩地层中,地层气体以欠平衡或重力置换方式进入井筒形成气侵。对重力置换气侵和欠平衡气侵类型及进气量进行分析;基于气液两相流理论建立井筒气侵模型,改变边界条件计算精细控压钻井气侵期间井底压力和泥浆池增量的变化;通过改变井口压力观察泥浆池总量或出口流量的变化,分析井底气侵方式。结果表明:井口增加回压后,井底压力由欠平衡状态转变为过平衡状态,泥浆池增量略微增加,而重力置换气侵泥浆池增量保持原趋势增加,由此可以判断出井底的气侵方式;该判定方法在现场试验中得到验证,对精细控压钻井安全快速控制气侵有指导意义。
重力置换气侵; 欠平衡气侵; 井底压力; 泥浆池增量; 气液两相流
碳酸盐岩地层裂缝和缝洞发育,安全密度窗口窄,钻井过程中井漏、溢流时常发生[1-4]。精细控压钻井通过改变井口回压实现对井底压力的精确控制,从而解决碳酸盐岩地层窄安全密度窗口、易喷易漏地层钻井问题[5-7]。精细控压钻井允许一定量的气体进入井筒,但不同的气侵类型对井底压力的影响不同,因此不同的气侵类型允许的泥浆池增量不同,如欠平衡气侵泥浆池增量为1 m3,而重力置换气侵泥浆池增量为3 m3。目前国内外对欠平衡气侵和重力置换气侵的判定方法研究尚不够深入[8-13],现场难以根据钻井参数判断井底气侵类型。笔者建立气侵期间井筒流动参数计算模型,通过观察改变井口回压后泥浆池总量变化趋势可以快速判断井底气侵类型,并利用现场试验对该判定方法进行验证。
1 气侵类型
控压钻井气侵主要分为欠平衡气侵、重力置换气侵、岩屑破碎气气侵与浓度差气侵[14]。岩屑破碎气气侵与浓度差气侵井底进气量较小,对控压钻井井底压力控制产生的影响可以不予考虑。重力置换气侵和欠平衡气侵井底进气量较大,对控压钻井井筒参数及井底压力控制影响较大,应具体讨论这两种气侵形成的原因及进气量。
1.1 欠平衡气侵
欠平衡气侵是指井底压力小于地层压力,地层中的气体在压差作用下由地层渗流到井筒的过程。进入井筒中的气体体积量与井筒与地层压力的平方差成正比,欠压值越大,进入井筒中的气体越多。图1为欠平衡气侵示意图,地层气体在压差的作用下进入井筒,形成气侵。
图1 欠平衡气侵示意图Fig.1 Schematic diagram of underbalanced gas kick
欠平衡气侵井底进气还受储层特性、气体特性及钻开储层厚度等因素影响[15],但对于已发生气侵的井,对井底进气量影响较大且能够控制的因素为井底与地层的欠压值。井筒内进气量可以根据气井的产能方程确定,
(1)
式中,pwf为井底压力,MPa;pe为地层压力,MPa;Qsg为气井产能即气侵量,m3/d;k为储层渗透率,10-3μm2;hi为钻开储层厚度,m;T为气层温度,K;μ为天然气黏度,mPa·s;re和rw分别为气井控制的外边缘半径和井筒半径,m;S为表皮系数;Z为气体偏差系数。
1.2 重力置换气侵
重力置换气侵是指地层中的气体与井筒中的钻井液在密度差的作用下,地层气体进入井筒,井筒中的钻井液进入到地层中的过程,如图2所示。重力置换依靠气体和钻井液的密度差作为气侵的动力,气侵量与地层孔隙度及渗透率有关,地层孔隙度与渗透率越大,重力置换气侵量越大。对于裂缝及溶洞型地层,重力置换气侵明显,容易转变为恶性漏失或井喷事故。
图2 重力置换气侵示意图Fig.2 Schematic diagram of gravity displacement gas kick
重力置换气侵量与储层孔隙度、渗透率及井眼尺寸有关。对于特定地层,井底重力置换气侵可以看作进气量恒定,
Qsg=f(k,φ,D)=C.
(2)
式中,φ为孔隙度;D为井眼尺寸,mm。
地层压力或井底压力越大,则井底气体密度越大,重力置换作用越小。一定范围内改变井底压力对气体密度影响不大,因此改变井口回压对重力置换作用影响较小,可以忽略其影响。
2 气侵井筒计算模型
地层气体进入井筒后,井筒流动由纯钻井液的单相流动变为钻井液和地层气体的气液两相流动,由于须求取井底压力变化及井筒气液两相流分布,因此建立井筒气液两相流连续性方程、动量方程和气体状态方程。
2.1 井筒气液两相流动计算模型
液相连续性方程为
(3)
对于非产气层段,气相连续性方程为
(4)
对于产气层段,气相连续性方程为
(5)
气液两相流动量方程为
(6)
式中,ρm为钻井液密度,kg/m3;vm为钻井液流速,m/s;ρg为气体密度,kg/m3;vg为气体真实速度,m/s;λ为含气率;z为空间坐标;t为时间坐标;Qg为井底进气量,m3/s;p为节点压力,MPa;τ0为流体与管壁之间的剪切力,N/m2;A为环空面积,m2。
2.2 初始条件和边界条件
初始条件是指气侵发生前井筒的流动状态,根据单相流体的流动求取初始流动状态时的井底压力、钻井液流动速度,
(7)
式中,p(0,i)为初始时刻井底压力,MPa;pf为循环摩阻,MPa;vm(0,i)为初始时刻环空钻井液流速,m/s;Qm为钻井液排量,m3/s。
边界条件是气侵模拟的边界约束条件,也是判断连续性方程和动量方程的收敛条件。气侵发生后,井底进气量和井口回压为两个主要的边界条件。根据气侵方式不同,选择不同的计算方法,控压钻井井口回压为井口边界条件,
(8)
式中,vg(j,0)为j时刻井底气体流速,m/s;λ(j,0)为j时刻井底含气率;vsg为井底气体表观速度,m/s;V为气泡体积,m3;p(j,H/ΔH)为井口节点压力,MPa;pa为井口节流压力,MPa。
常规钻井中,井口回压为恒定大气压。若采用欠平衡钻井或控压钻井,在发生气侵后,随着气体的上返,井底压力下降,可增加井口回压控制气侵。施加井口回压时,模拟计算的边界条件及井筒节点参数会发生改变,计算公式为
(9)
式中,pa(j,H/ΔH)为j时刻井口节点压力,MPa;pa(j-1,H/ΔH)为j-1时刻井口节点压力,MPa;Δpa为井口回压的变化量,MPa;p(j,i)和p(j-1,i)分别为j时刻和j-1时刻i节点处的压力,MPa;V(j,i)和V(j-1,i)分别为j时刻和j-1时刻i节点处的气体体积,m3;Z(j,i)和Z(j-1,i)分别为j时刻和j-1时刻i节点处的偏差因子;T(i)为i节点处的温度,K;p(j,0)和p(j-1,0)分别为j时刻和j-1时刻井底压力,MPa;ρgh为井口回压压缩气体增加的静液柱压力,MPa。
2.3 求解步骤
(1)网格离散。对井筒按照空间和时间进行网格划分,将井筒深度划分为空间网格,溢流时间划分为时间网格。由于井口位置气体膨胀率较大,要求计算精度较高,井底位置气体膨胀率较小,要求计算精度较低,因此可以采用从井口到井底的等比数列空间网格。
时间网格的确定为气体从一个空间网格上升到另一个空间网格的时间,计算式为
(10)
式中,t(j)为j时刻的时间网格长度,s;z(i)为i空间网格长度,m;vg(i)为气体在i空间网格的真实速度,m/s。
(2)节点求解。根据初始条件和边界条件假设井底压力,计算井底进气速度,采用试算法计算节点含气率和压力,以井口压力为迭代约束条件验证假设的井底压力是否合适,若不合适再根据不同情况进行假设求解,直到假设的井底压力满足井口约束条件为止。
气侵井筒模拟计算步骤如图3所示。
图3 气侵井筒模拟计算步骤Fig.3 Gas kick wellbore simulation steps
3 不同气侵方式钻井参数变化特征
根据前述的井筒气侵计算模型对井筒发生重力置换气侵和欠平衡气侵进行模拟研究,模拟井深为5.5 km,钻井液密度为1.18 g/cm3,钻井液排量为10 L/s,储层渗透率为1×10-3μm2,钻开储层厚度为1 m,地层压力为64.5 MPa,钻头直径为15.24 cm,钻杆外径为8.89 cm,钻井液黏度计300转读数为58,600转读数为90,甲烷气体临界温度为191.05 K,甲烷气体临界压力为4.6 MPa,泥浆池报警值为0.3 m3。
3.1 井底含气率变化
当井底压力为64 MPa钻遇不同地层压力时,井底分别形成0.5、1.0和1.5 MPa的欠压值,形成欠平衡气侵,井底含气率随时间的变化如图4所示。
图4 欠平衡气侵井底含气率变化Fig.4 Downhole gas void change of underbalanced pressure gas kick
欠平衡气侵进气速度与井底欠压值有关,欠压值越大进气速度越大。当地层压力为64.5 MPa时,即井底欠压值为0.5 MPa,井底含气率为0.02;地层压力为65.5 MPa时,井底欠压值为1.5 MPa,井底含气率为0.06。
欠平衡气侵发生后,井底欠压值越大,井底进气速度越大,井底含气率随欠平衡气侵发生时间逐渐增加。当地层压力为64.5 MPa时,开始发生气侵时井底含气率约为0.02,在气侵40 min后,井底含气率达到0.2。
重力置换气侵井底含气率与进气量的关系如图5所示。重力置换气侵井底含气率不随时间发生变化,只与井底进气量有关,井底进气量越大,井底含气率越大。
图5 重力置换气侵井底含气率变化Fig.5 Downhole gas void change of gravity displacement gas kick
3.2 井筒含气率变化
井底压力为64 MPa,井底压力与地层压力负压差分别为0.5、1和1.5 MPa,气侵发生65 min时,不同井深处的含气率变化如图6所示。
图6 欠平衡气侵井筒含气率变化Fig.6 Wellbore gas void during underbalanced gas kick
气侵发生65 min时,欠压值越大进气速度越快,气体越先到达井口段。深度大于2.0 km井段含气率随着井深增加而增大,这是由于随着气侵的发生井底进气速度越来越快,越接近井底含气率越高。深度小于1.0 km井段由于气体的快速膨胀,井筒含气率快速增加。
不同井底进气量时,重力置换气侵发生60 min后井筒含气率的变化如图7所示。
由图7可以看出,井底进气量越大井筒含气率越大。深度大于2.0 km井段,井筒含气率变化较小,深度小于1.0 km井段,由于气体的膨胀,井筒含气率快速增加。
图7 重力置换气侵井筒含气率变化Fig.7 Wellbore gas void during gravity displacement gas kick
3.3 井底压力及泥浆池增量变化
当井底压力为64 MPa钻遇不同的地层压力时,井底分别形成0.5、1和1.5 MPa的欠压值,形成欠平衡气侵后,井底压力及泥浆池增量的变化如图8所示。
图8 欠平衡气侵特征Fig.8 Characteristics of underbalanced gas kick
由图8可知,进入井筒中的气体运移到井口需要一定时间,且气体在深度大于2.0 km井段膨胀性较小,因此气侵开始时刻井底压力下降较慢,泥浆池增量较小。在气侵30 min后,由于井底进气速度的加快和气体运移到井口位置的双重作用,井底压力快速下降,泥浆池增量快速增加。初始欠压值越大,井底压力下降越快,泥浆池增量越大。
不同井底进气量时,重力置换气侵的井底压力及泥浆池增量变化如图9所示。由图9可知,重力置换气侵进气速度恒定,井底压力的下降和泥浆池增量仅与气体的膨胀有关。在气侵发生后40 min内,由于气体的膨胀性较小,井底压力和泥浆池增量变化不大;40 min后,气体逐渐运移到井口,造成井底压力及泥浆池增量的快速变化。井底进气速度越大,井底压力下降越快,泥浆池增量越大。
图9 重力置换气侵特征Fig.9 Characteristics of gravity displacement gas kick
4 重力置换气侵与欠平衡气侵判定方法
欠平衡气侵和重力置换气侵都会造成井底压力下降,井口泥浆池总量增加,虽然两种气侵方式造成的井底压力和泥浆池增量变化速率不同,但钻井现场难以进行准确判断。两种气侵方式的井底含气率变化趋势不同,但无法通过仪器测出井底含气率的变化,因此现在钻井现场难以检测出井底发生气侵的类型。
控压钻井技术能够精确控制井底压力剖面,很好地解决窄密度窗口钻井问题,同时对气侵类型的判断提出了更高的要求,控压钻井井口回压的改变也为判断井底气侵类型提供了一种可行性方案。通过对两种气侵类型发生的原因及气侵特征分析可知,欠平衡气侵主要受井底压力的影响,而井底压力对重力置换气侵影响较小,因此可以通过改变井底压力的方式判断气侵类型。通过改变井口回压,井底压力由欠平衡气侵状态变为过平衡后,欠平衡气侵将不再发生,而重力置换气侵将继续维持原来的进气速度发生气侵。
井底发生欠平衡气侵时,根据表1中的模拟参数,当泥浆池增量超过0.3 m3,检测出溢流,增加井口回压,模拟改变回压后井底压力及泥浆池增量的变化,如图10所示。
由图10可知,钻遇地层压力为64.5 MPa,井底为欠平衡气侵,泥浆池增量超过0.3 m3时井口增加1 MPa回压,井底压力变为64.8 MPa,此时井底为过平衡状态。施加回压后,由于井筒内气体向上运移膨胀,井底压力略微下降,泥浆池增量少量增加。
当井口施加0.3 MPa回压,井底压力变为64.07 MPa,此时井底仍为欠平衡状态,如图11所示。在施加井口回压后,同一时刻下井底压力下降速度和泥浆池增量增加速度比未施加井口回压的气侵小。
图10 井口施加回压后井底为过平衡气侵特征Fig.10 Downhole overbalanced pressure after wellhead increased backpressure
图11 井口施加回压后井底欠平衡气侵特征Fig.11 Downhole underbalanced pressure after wellhead increased backpressure
井底发生重力置换气侵时,根据表1中的模拟参数,井底进气量为100 m3/d,泥浆池增量超过0.3 m3
后,检测出溢流,井口增加回压0.3 MPa,结果如图12所示。井底压力增加0.3 MPa后,继续下降,而泥浆池增量保持原来的趋势继续增加。
图12 井口施加回压后重力置换气侵特征Fig.12 Downhole gravity displacement gas kick after wellhead increased backpressure
在控压钻井现场试验与应用中,更多的是采用增加回压观察一个迟到时间后出口流量的变化。控压钻井井口装置安装有高精度的质量流量计,若出口流量大于入口流量,则泥浆池增量处于增加的趋势,由此可以更直观地判断井底气侵类型。不同气侵类型井口增加回压后气侵特征变化如表1所示。
井口施加回压后井底仍为欠平衡时,在井口设备允许的范围内,可以通过继续增加回压直到井底为过平衡的方式判断井底气侵类型。
表1 井口增加回压后不同气侵类型井底压力及泥浆池总量变化
5 实例应用
将通过改变井口回压观察泥浆池增量或出口流量变化判断井底气侵类型的方法在塔中地区进行试验。试验井目的层为上奥陶统良里塔格组,设计完钻井深为6.740 km,垂深5.005 km,靶点A井深为5.183 km,靶点B井深为6.740 km,井身结构如图13所示。
试验井段立压为18.7 MPa,排量为13.5 L/s,钻井液迟到时间为102 min,钻井液密度为1.10 g/cm3,钻井液黏度计600转读数为38,300转读数为24,黏度为11 mPa·s,正常钻进期间井口回压保持在2.5~3.0 MPa,接单根期间井口回压为4.5 MPa。在8∶40时,钻头钻达6.139 86 km,井口出口流量开始大于入口流量,检测出井底发生溢流,井口回压为2.8 MPa, 在8∶50时井口回压增加到3.4 MPa,一个迟到时间后观察井筒出口流量和入口流量的变化,结果如图14所示。
井口回压增加一个迟到时间后,在15∶40出口流量和入口流量相等,之后出口流量一直保持和入口流量相等。由此可知,在8∶50增加井口回压后,井底停止进气,从而判断在8∶50时通过井口增加回压使井底变为过平衡状态,井底停止进气。根据气侵判定方法可知,在8∶40发生的气侵为欠平衡气侵,泥浆池增量应该控制在1 m3内,井筒内气体循环完毕后,井口回压降低到2.5~3.0 MPa,进行正常钻进。
图13 井身结构Fig.13 Wellbore structure
图14 溢流检测试验Fig.14 Overflow testing
6 结 论
(1)欠平衡气侵初始欠压值越大,井底压力下降越快,泥浆池总量增加越快;重力置换气侵进气量越大,井底压力下降越快,泥浆池总量增加越快。
(2)精细控压钻井改变井口回压后观察泥浆池总量或出口流量变化可以判断井底发生气侵的类型。增加回压后,欠平衡气侵转变为过平衡,泥浆池总量略微增加或出口流量基本和入口流量相等;重力置换气侵泥浆池总量保持原趋势增加。
(3)工程应用实例证明通过增加井口回压观察泥浆池总量变化判断气侵类型的方法可行,为精细控压钻井判断气侵类型提供了理论依据。
[1] 周英操,杨雄文,方世良,等. 窄窗口钻井难点分析与技术对策[J]. 石油机械, 2010,38(4):1-7. ZHOU Yingcao, YANG Xiongwen, FANG Shiliang, et al. An analysis of the difficulties in narrow window drilling and technical solutions[J]. China Petroleum Machinery, 2010,38(4):1-7.
[2] 俞新永,周建东,滕学清.塔里木轮南奥陶系碳酸盐岩高压油气藏水平井及大斜度井欠平衡钻井技术[J].天然气工业,1999,19(2):63-68. YU Xinyong, ZHOU Jiandong, TENG Xueqing. Under balance pressure drilling techniques of the horizontal well and high angle deviated hole in Ordovician carbonate rock high-pressure reservoirs in Lunnan region of Tarim Basin[J].Natural Gas Industry, 1999,19(2):63-68.
[3] SCOTT Mcmillan, AN Wenhua, ZHANG Yaoming, et al. Overcoming drilling challenges in northwest China using an innovative MPD technique [R]. SPE 130317, 2010.
[4] URBIETA A, PEREZTELLEZ C, LUPO C, et al. Successful application of MPD technique in a HP/HT well focused on performance drilling in southern Mexico deep fractured carbonates reservoirs[R]. SPE 122200, 2009.
[5] MALIK H K, LITWIN R J, JONES BC, et al. Successful use of managed pressure drilling (MPD) to counteract loss circulation and HP gas zones: a case study from red rock south wapiti field in Grand Prairie Alberta[R]. SPE 137428, 2010.
[6] 石林,杨雄文, 周英操,等.国产精细控压钻井装备在塔里木盆地的应用[J]. 天然气工业, 2012,32(8):6-10. SHI Lin, YANG Xiongwen, ZHOU Yingcao, et al. Application of China-made precise managed pressure drilling equipment in the Tarim Basin[J].Natural Gas Industry,2012,32(8):6-10.
[7] RAMIREZ S, AGUILAR R, HERRERA R, et al. Managed pressure drilling (MPD) solves conventional drilling problems at the Xanab offshore field, which presents severe drilling fluid losses[R]. SPE 22713, 2011.
[8] 周英操, 高德利, 刘永贵. 欠平衡钻井环空多相流井底压力计算模型[J]. 石油学报, 2005,26(2):96-99. ZHOU Yingcao, GAO Deli, LIU Yonggui. New model for calculating bottom hole pressure of multiphase flow in annulus of underbalanced straight well[J]. Acta Petrolei Sinica, 2005,26(2):96-99.
[9] 刘永贵, 周英操, 王广新,等. 欠平衡钻井环空岩屑对井底负压的影响[J]. 石油学报, 2005,26(6):96-98. LIU Yonggui, ZHOU Yingcao, WANG Guangxin, et al. Effect of annular cuttings on negative pressure in under-balanced drilling process[J]. Acta Petrolei Sinica,2005,26(6):96-98.
[10] 周英操, 王广新, 翟洪军,等. 欠平衡钻井技术在大庆油田卫深5井中的应用[J]. 石油学报, 2003,24(6):90-93. ZHOU Yingcao, WANG Guangxin, ZHAI Hongjun, et al. Application of under-balanced drilling technology in Weishen-5 well of Daqing Oilfield[J]. Acta Petrolei Sinica, 2003,24(6):90-93.
[11] 杨虎, 鄢捷年. 欠平衡钻井井壁失稳和储层损害的最大负压差确定研究[J]. 天然气工业, 2005,25(5):53-55. YANG Hu, YAN Jienian. Determination of the maximum negative pressure difference causing wall destabilization and reservoir damagein underbalanced drilling[J]. Natural Gas Industry, 2005,25(5):53-55.
[12] 刘绘新, 李锋. 裂缝性储层井控技术体系探讨[J]. 天然气工业, 2011,31(6):77-80. LIU Huixin, LI Feng, Well control technologies for fractured gas reservoirs[J]. Natural Gas Industry,2011,31(6): 77-80.
[13] 史玉胜, 何世明, 林文秀,等. 微流量控压钻井技术研究现状[J]. 钻采工艺, 2011,34(3):12-14. SHI Yusheng, HE Shiming, LIN Wenxiu, et al. Research on micro-flux control pressure drilling technology[J]. Drilling & Production Technology, 2011,34(3):12-14.
[14] 隋秀香, 李相方, 齐明明,等. 高产气藏水平井钻井井喷潜力分析[J]. 石油钻探技术, 2004,32(3):34-35. SUI Xiuxiang, LI Xiangfang, QI Mingming, et al. Possibility analysis of blowout during horizontal drilling in gas reservoirs with high productivity[J]. Petroleum Drilling Techniques, 2004,32(3):34-35.
[15] 李相方, 庄湘琦. 关井压力恢复和读取时机分析[J]. 石油学报, 2002,23(5):110-112. LI Xiangfang, ZHUANG Xiangqi. Analysis on the pressure build-up and its recording time after shut-in[J]. Acta Petrolei Sinica,2002,23(5):110-112.
(编辑 李志芬)
A method for characterization and identification of gas kicks caused by underbalanced pressure and gravity displacement
ZHANG Xingquan1,2, ZHOU Yingcao1, LIU Wei1, GUO Qingfeng1
(1.CNPCDrillingResearchInstitute,Beijing102206,China;2.PetroChinaResearchInstituteofPetroleumExploration&Development,Beijing100083,China)
During well drilling in fracture-vuggy carbonate oil and gas reservoirs, gas may enter into wellbore in the forms of underbalanced pressure gas kick and gravity displacement gas kick. In order to analyze the gas kick velocity and identify the two different gas kick types, a wellbore gas kick model was proposed based on gas-liquid two phase flow to calculate downhole pressure and the incremental amount of drilling fluid in mud pit. It is possible to determine the types of gas kicks by observing the volume change in the mud pit via changing wellhead back pressure. The results show that, when the wellhead back pressure is increased, the downhole pressure may change and result in a transfer from underbalanced to overbalanced, in which only a slight volume increase in mud pit can be observed. On the other hand, the volume in mud pit may maintain its steady increase trend if a gravity displacement gas kick occurs downhole. This can be used to identify the gas kick types, and the method is verified in the field trials, which is useful for the managed pressure drilling to control gas kick quickly and safely.
gravity displacement gas kick; underbalanced pressure gas kick; downhole pressure; mud pit increment; gas-liquid two-phase flow
2014-04-10
国家自然科学基金项目(51274221);国家科技重大专项(2011ZX05021-003)
张兴全(1986-),男,博士研究生,主要从事控压钻井技术研究。E-mail:zxq336998@163.com。
1673-5005(2015)01-0095-08
10.3969/j.issn.1673-5005.2015.01.014
TE 21
A
张兴全,周英操,刘伟,等. 欠平衡气侵与重力置换气侵特征及判定方法[J]. 中国石油大学学报:自然科学版,2015,39(1):95-102.
ZHANG Xingquan, ZHOU Yingcao, LIU Wei, et al. A method for characterization and identification of gas kicks caused by underbalanced pressure and gravity displacement[J]. Journal of China University of Petroleum(Edition of Natural Science),2015,39(1):95-102.