考虑黏结特性的水电地下洞室块体地震响应分析

李文倩，佟大威，王 振，鲁文妍 2，朱晓斌，钟登华

（1. 天津大学水利工程仿真与安全国家重点实验室，天津 300072；2. 南京水利科学研究院，南京 210029）



考虑黏结特性的水电地下洞室块体地震响应分析

李文倩1，佟大威1，王 振1，鲁文妍1 2，朱晓斌1，钟登华1

（1. 天津大学水利工程仿真与安全国家重点实验室，天津 300072；2. 南京水利科学研究院，南京 210029）

摘要：在水电工程地下洞室的施工及运行过程中，地震作用下块体可能出现滑移失稳甚至破坏，将会给工程带来严重威胁，因此研究地震作用下地下洞室块体的抗滑稳定性具有十分重要的意义.实际工程中块体与周边岩体界面上存在黏结作用，而ABAQUS的接触属性中自带的抗滑模型未考虑该特点，为此建立了考虑黏结强度影响效应的地下洞室块体修正抗滑模型，完善了块体与基岩间的静态作用与动态相对滑动的影响机理.通过验证模型与工程算例，研究了块体在静态及动态荷载激励下的响应特征与稳定性，实现了抗滑模型合理性与有效性的验证.对比经典抗滑模型与修正抗滑模型的分析结果，在重力阶段最大竖向位移分别为1.04 cm和0.05 cm，在地震过程中的平均竖向位移与最大竖向位移分别为2.24 cm、3.91 cm和2.07 cm、3.38 cm.上述结果反映了块体在黏结特性的作用下具有更好的稳定性.依次改变抗滑接触属性中的关键参数，讨论了摩擦系数及黏聚力对块体响应与接触面动态特性的影响.

关键词：岩石块体；地下洞室；抗滑模型；动力响应分析；ABAQUS

我国大型水力发电工程多位于高山峻岭之间，且地下洞室往往处于极其复杂的地质条件下，其周边分布的由地质构造面、随机裂隙及临空面等切割而成的曲面块体是影响地下洞室稳定性的主要安全隐患，在强震作用下块体一旦发生滑移失稳将对工程产生极其严重的影响，因此对大型地下洞室中局部块体进行动力稳定性的模拟研究具有重要的工程意义.

大型地下洞室块体结构的稳定性不仅受到工程施工过程与外界激励作用的影响，同时还受到块体自身具有的动态稳定机理与接触特性的影响.国内较早的有关块体稳定分析方法由石根华［1］提出，他从块体稳定分析的几何原理与岩体的影响因素出发建立滑动体的数学模型.在国外，Goodman等［2-3］更早地对块体理论与接触问题在岩体工程中的应用进行了系统的论述，并广泛应用于实际工程中［4］.地下洞室块体结构动力分析方法主要有连续数值方法［5］、非连续数值方法［6-8］与耦合数值方法［9］.目前在块体领域的研究中常基于块体理论并采用离散元法将岩石和节理分开考虑，以反映岩体在裂隙面上的滑动、张开等特性［10］.郑银河等［11］基于块体理论识别块体，并在抗滑力（包括摩擦力和黏聚力）的作用下研究了考虑岩桥破坏的块体稳定性.付晓东等［12］基于非连续变形分析（DDA）在接触面采用刚塑性模型处理节理黏聚力的作用，研究了二维滑坡块体的安全系数［12］.钟登华等［13］在曲面块体的快速识别方法的基础上，运用有限元理论并采用Coulomb摩擦模型研究了复杂地质条件下的曲面块体的地震响应.综上所述，目前基于有限元理论的三维块体动态响应特性的模拟较少，且相关研究忽视了岩体间黏结强度对接触面动态特性的影响效应，无法有效地体现实际工程中的接触面动力情况.

针对以上研究中存在的不足，本文考虑了实际工程中节理裂隙的不完整性与黏结特性，基于ABAQUS的二次开发平台，编制了摩擦接触子程序，在ABAQUS自带抗滑模型的基础上建立并实现了适用于大型地下洞室稳定性分析的“块体周边接触面抗滑模型”.该修正抗滑模型充分考虑了块体与基岩间的相互影响与作用机理，将离散元中关于裂隙面抗滑特性研究的理论在有限元方法中更好的实现，模拟实际工程中块体周边结构面上的摩擦与黏结特性，以期更加合理地反映局部块体及基岩在强震作用下的动态变化过程.首先在理论上给出了接触面修正抗滑模型的基本原理与适用条件；其次，分别从验证模型与工程算例出发研究了“块体周边接触面抗滑模型”的适用性与有效性；最后，研究了修正抗滑模型中相关参数变化对块体响应与接触面动态特性的影响规律.

1　块体周边接触面抗滑模型

大型水电工程地下洞室中块体与基岩结构间的接触状态可以分为脱开、黏结和滑动3种情况［14］.对于地下洞室中块体与结构间抗滑特性的模拟的关键在于，其一判断块体与基岩结构间的接触状态；其二为在已知存在接触的条件下采用合理的抗滑模型.

对于块体与结构间接触状态的判断，主要基于块体与基岩结构间接触面上的法向变形（位移）的大小及方向来确定，接触面间法向模型表达式为

式中：σcn为接触面法向应力，Pa；kn为接触面法向刚度，Pa/m；un为接触面上法向相对变形（位移），m.

块体在滑动过程中受到的抗滑力包括摩擦力及黏聚力.结合广义库伦模型可以将块体与基岩结构间接触面上的抗滑力描述为［14-15］

式中：τct为接触面切向应力，Pa；ut为接触面上切向相对变形（位移），m；kt为接触面切向刚度，Pa/m；f为接触面上的摩擦系数；c为块体与岩体结构接触面上的黏聚力，Pa；sgn（·）为符号函数，给定接触面抗滑力的符号与方向.

ABAQUS有限元分析中针对接触面间摩擦模型的极限剪应力计算默认采用Coulomb定律［16］，其抗滑状态的判断公式中不存在接触面间的黏聚力项c，而在实际工程中地下洞室岩体与由裂隙形成的块体间应存在黏结作用，本文对此进行了完善.基于ABAQUS提供的开放的二次开发平台Subroutine，运用Fortran语言编写代码，根据“块体周边接触面抗滑模型”，自定义接触面摩擦模型子程序FRIC，实现与ABAQUS中参数的对接，以反映剪应力与接触面作用的压力之间的关系.

2　抗滑模型验证

2.1验证模型

建立图1中的数值模型进行ABAQUS摩擦子程序的合理性与有效性的验证.上部滑块尺寸为10 m×10 m×10 m（长×宽×高），用于模拟实际地下洞室中的块体结构，下部基础结构模拟地下洞室岩体，尺寸为100 m×50 m×10 m.两部分结构所赋材料属性相同，即密度为2 500 kg/m3，弹性模量为30.0 GPa，泊松比为0.20.下部基础底面采用全约束处理，侧面采用法向约束处理.施加荷载时先对上部滑块施加一个方向竖直向下、大小为10 m/s2的重力加速度，重力平衡后在水平向再施加一个时长为8 s且最大幅值为10 mN的动态推力，在1～6 s时由0 mN线性增加至最大推力，并保持最大.

图1　数值模型Fig.1 Numerical model

2.2验证结果分析

基于上述数值验证模型，分别按照经典抗滑模型与本文采用的修正抗滑模型对滑块与基础接触面上的抗滑特性进行模拟，设定接触面摩擦系数为0.2，黏聚力为10 kPa，剪切刚度为10 mPa/m.根据上述数值模型与材料属性，基于理论公式（3）可以计算获得考虑黏结强度效应的块体抗滑模型条件下的解析解，其中上部滑块将在总加载时间为5.0 s时开始滑动，滑块最终恒定加速度a为1.6 m/s2，而接触面上最大剪应力τ 为60 kPa.

式中：Fmax为加载的水平推力最大值，N；u为接触面摩擦系数；P为接触面法向压力，N；A为接触面间接触面积，m2.

图2和图3分别给出了经典抗滑模型、修正抗滑模型与理论解析解3种条件下，滑动块体加速度和接触面上的剪应力曲线对比.由图2可知，在经典抗滑模型条件下，滑块在总加载时间的第4.5 s开始出现滑动，而本文采用的修正抗滑模型考虑了块体与基础结构间的黏聚力作用，从而增加了块体滑动的极限抗滑力，滑块在总加载时间的第5.0 s发生滑动，这与理论分析保持一致.而在总加载时间的第7.0 s时水平推力幅值达到最大，此时采用修正抗滑模型计算的滑块最终恒定运动加速度为1.6 m/s2，与解析解相同，而经典抗滑模型的计算结果为2.0 m/s2，从计算结果可知修正的抗滑模型中考虑黏聚力可使滑块具有更好的抗滑稳定性.

由图3可知，采用两种抗滑模型计算接触面剪应力均是在总加载时间的第2.0 s开始随着水平推力的增加而逐渐加大，直至经典模型的第4.5 s和修正模型的第5.0 s时滑块开始滑动，其界面剪应力基本保持不变.考虑到修正模型中的黏聚力作用，其对应的接触面剪应力为59.6 kPa，与解析解60.0 kPa基本相同，且大于经典抗滑模型对应的49.8 kPa，因此，本文所采用的修正抗滑模型及摩擦子程序的有效性得以验证.

图2　滑动块体加速度曲线Fig.2 Acceleration curve of sliding block

图3　滑动块体接触面剪应力曲线Fig.3 Shear stress on contact surface of sliding block

3　工程应用

3.1工程概况与模型建立

所选工程位于我国西南地区澜沧江下游，引水发电系统位于山体左岸，地下式厂房结构，包括主副厂房、主变室、尾水闸门室及尾水调压室等.建立整体计算模型（见图4），模型尺寸为837 m×446 m× （299～551）m（长×宽×高），主变室尺寸为348 m× 19 m×（23.6～38.6）m.其中，采用NURBS-TINBRep混合数据结构，建立复杂的地形体，结合几何建模技术建立地下洞室群，并在此基础上基于块体理论耦合随机结构面网络模拟技术，采用识别曲面块体的三大定理：封闭性、完备性及唯一性，快速建立地下洞室区域的曲面块体［17］.基于搜索到的曲面块体，选取位于主变室拱顶部位的块体为研究对象，建立的地形模型、洞室模型及块体模型见图4.

图4　块体与整体模型位置的注释图Fig.4 Position annotation of block and whole model

模型的本构采用理想弹塑性模型，采用Drucker-Prager屈服准则，岩体的弹性模量为16 GPa，密度为2.71 g/cm3，泊松比为0.23，抗拉强度为1.6 mPa，内摩擦角为52°，围岩阻尼采用Rayleigh阻尼，接触属性中摩擦系数为0.3，黏聚力为0.5 mPa.主要计算工况包括重力工况、开挖工况及地震工况.计算输入的地震波采用唐山余震天津医院地震波，峰值加速度为145.8 cm/s2，地震持续时间为19.20 s，经过调幅、地震波反演、基线校正及滤波等处理，基岩处的地震波加速度时程曲线见图5，并在基岩模型底面的全部节点处施加x向地震激励.

图5　基岩处地震加速度时程曲线Fig.5 Seismic acceleration time history in rock bottom

3.2对比分析结果

图6和图7分别给出了两种抗滑模型条件下洞顶水平x向加速度及块体水平x向相对基岩底部加速度.由图7可知，经典模型条件下块体水平x向相对基岩底部加速度的均方根值为0.437 m/s2，而修正模型条件下对应响应为0.427 m/s2，也体现出修正模型考虑黏结特性后对块体动态响应的影响特性.

图6　洞顶水平x向加速度Fig.6 Horizontal acceleration of roof in x direction

图7　块体水平x向相对加速度Fig.7 Horizontal relative acceleration of block in x direction

图8给出了两种抗滑模型条件下块体运动特征点处的竖向位移对比.前2.0 s为重力与开挖的静力作用阶段，经典抗滑模型条件下两阶段的最大竖向位移分别为1.04 cm和1.19 cm，而修正的抗滑模型条件下对应值分别为0.05 cm和0.49 cm.在地震阶段，经典抗滑模型条件下平均竖向位移与最大竖向位移分别为2.24 cm和3.91 cm，而修正抗滑模型对应值分别为2.07 cm和3.38 cm，说明在静、动力荷载激励下，块体在修正抗滑模型时发生的运动程度要小于经典抗滑模型的响应结果，在一般情况下，上述结果反映了块体在考虑黏结特性的修正抗滑模型作用下具有更好的稳定性.块体在地震荷载激励作用下相对周边岩体会发生明显错动，见图9.

图8　块体竖向位移Fig.8 Vertical displacement of block

图9　块体相对岩体错动示意Fig.9 Dislocation of block relative to rock

图10和图11分别给出了两种抗滑模型条件下块体水平x向绝对与相对位移时程.由图中数值可知，基于修正抗滑模型计算结果整体均小于经典抗滑模型计算结果，在整个地震阶段经典与修正抗滑模型最大水平x向位移分别为7.01 cm和6.78 cm，平均水平x向位移分别为2.26 cm和2.08 cm.从块体的不同响应结果可以看出，由于考虑了块体与周边基岩间的黏结特性，使得块体在修正抗滑模型条件下的动态响应整体上要明显小于在经典抗滑模型条件下的动态响应.

图10　块体水平x向绝对位移Fig.10 Horizontal absolute displacement of block in x direction

图11　块体水平x向相对位移Fig.11 Horizontal relative displacement of block in x direction

图12和图13从块体与周边岩体接触面的动态特性出发，分别给出了图中节点位置在两种抗滑模型条件下切向剪应力与界面脱开情况对比.由图12可知，基于修正抗滑模型计算获得的剪切应力曲线在整体趋势中大于经典抗滑模型的计算结果，其中前者最大剪应力值达到0.22 mPa，而后者为0.16 mPa，前者在整个计算过程中的平均剪应力为0.12 mPa，而后者为0.10 mPa，其进一步说明了在动力计算中，块体与基岩间的黏聚力对结构动态特性起到了有效且明显的作用.

图12　块体与周边岩体界面剪应力Fig.12 Shear stress on surface between block and rock

由图13可知，在整个动力计算过程中所选取节点处接触面上出现脱开的程度在修正抗滑模型作用下要更小，这说明采用考虑接触面上黏结特性的抗滑模型时块体与周边岩体具有更好的稳定性.

图13　块体与周边岩体界面脱开Fig.13 Disengage on surface between block and rock

4　参数变化对响应结果的影响

通过改变修正抗滑模型中摩擦系数与黏聚力关键参数，研究其变化对响应结果的影响.重点研究了块体与岩体摩擦系数取值范围为0.1～0.7、界面黏聚力取值范围为0.1～0.5 mPa时块体与基岩间接触面的动态特性变化规律.图14和图15分别给出了块体在重力阶段与地震阶段的竖向位移随参数变化的响应结果.由图14可知，在重力阶段，不同黏聚力条件下随着摩擦系数的增加，块体的竖向位移呈现减小的变化趋势，而在选取相同的摩擦系数时，块体竖向位移随着黏聚力的增加也逐渐变小，以上各现象特别是在摩擦系数小于0.4且黏聚力值小于0.3 mPa时尤为明显.然而当摩擦系数达到0.7及黏聚力达到0.5 mPa时，改变对应参数对块体竖向位移的变化影响已经较小，说明该参数在此范围时对块体抗滑起到了决定性作用，而改变其余参数对块体响应影响不大.

图14　重力阶段块体竖向位移Fig.14 Vertical displacement of block in gravity stage

图15中模型在地震作用下块体的竖向位移响应结果随着修正抗滑模型中摩擦系数与黏聚力的增加而呈现整体上的减小趋势，各工况下块体响应幅度（均方根值）最大为2.22 cm，最小为1.69 cm，在摩擦系数与黏聚力分别处于较低水平时，改变对应参数对块体整体动态响应具有较为明显的影响.而当摩擦系数为0.4且黏聚力达到0.5 mPa这一较大量级时，摩擦系数的改变对块体的竖向位移影响效果基本不变，说明在此时的动力条件下块体整体稳定性主要受到抗滑模型中黏聚力项的影响，与摩擦系数的关联较弱.

图15　地震阶段块体竖向位移Fig.15 Vertical displacement of block in seismic stage

5　结 论

（1）本文采用的修正抗滑模型考虑了块体与基岩间的相互影响与作用机理，将离散元中关于裂隙面抗滑特性研究的理论在有限元理论中更好的实现，模拟实际工程中块体接触面上的摩擦与黏结特性，更合理地反映了局部块体及基岩在强震作用下的动态稳定与变化过程.

（2）通过验证模型结果可知，在动态荷载作用下基于修正抗滑模型的计算结果要明显小于通过经典抗滑模型计算的动力响应结果，并基于理论公式得到修正抗滑模型条件下的解析解，与修正抗滑模型条件下验证模型的响应结果基本吻合，说明本文采用修正抗滑模型的有效性与正确性.

（3）基于某工程算例，全面分析了计算模型在初始阶段、开挖阶段及地震阶段，两种抗滑模型条件下块体结构动态响应及块体周边接触面上的动态接触特性，通过对比经典抗滑模型与修正抗滑模型的响应结果，进一步分析了本文采用的修正抗滑模型能够合理地模拟地下洞室块体的动态响应特性.

（4）通过改变修正抗滑模型中关键参数的取值，研究了其对块体结构在静、动力作用下响应结果的影响.研究结果表明，地下洞室块体结构在摩擦系数与黏聚力增大的情况下，结构的响应呈现逐渐减小的趋势；而在摩擦系数或黏聚力较大时，另一参数的改变已经不能显著影响块体的整体稳定性.

参考文献：

［1］ 石根华. 岩体稳定分析的几何方法［J］. 中国科学，1981（4）：487-495. Shi Genhua. Geometric method of stability analysis of rockmass［J］. Science China，1981（4）：487-495（in Chinese）.

［2］ Goodman R E，Shi G H. Block Theory and Its Application to Rock Engineering［M］. London：Prentice-Hall Inc，1985.

［3］ Goodman R E，Taylor R L，Brekke T L. A model for the mechanics of jointed rock［J］. Journal of Soil Mechanics and Foundations，1968，94（3）：637-659.

［4］ 黄正加，乌爱清，盛 谦. 块体理论在三峡工程中的应用［J］. 岩石力学与工程学报，2001，20（5）：648-652. Huang Zhengjia，Wu Aiqing，Sheng Qian. Usage of block theory in the Three Gorges projects［J］. Chinese Journal of Rock Mechanics and Engineering，2001，20（5）：648-652（in Chinese）.

［5］ 李海波，马行东，李俊如，等. 地震荷载作用下地下岩体洞室位移特征的影响因素分析［J］. 岩土工程学报，2006，28（3）：358-362. Li Haibo，Ma Xingdong，Li Junru，et al. Study on influence factors of rock cavern displacement under earthquake［J］. Chinese Journal of Geotechnical Engineering，2006，28（3）：358-362（in Chinese）.

［6］ Beyabanaki S A R，Mikola R G，Hatami K. Threedimensional discontinuous deformation analysis（3-D DDA）using a new contact resolution algorithm［J］. Computers and Geotechnics，2008，35：346-356.

［7］ 乌爱清，丁秀丽，卢 波，等. DDA方法块体稳定性验证及其在岩质边坡稳定性分析中的应用［J］. 岩石力学与工程学报，2008，27（4）：664-672. Wu Aiqing，Ding Xiuli，Lu Bo，et al. Validation for rock block stability and its application to rock slope stability evaluation using DDA method ［J］. Chinese Journal of Rock Mechanics and Engineering，2008，27（4）：664-672（in Chinese）.

［8］ Oliveira D A F，Indraratna B. Comparison between models of rock discontinuity strength and deformation ［J］. Journal of Geotechnical and Geoenvironmental Engineering，2010，136：864-874.

［9］ Kaidi S，Rouainia M，Ouahsine A. Stability of breakwaters under hydrodynamic loading using a coupled DDA-FEM approach ［J］. Ocean Engineering，2012（5）：62-70.

［10］ 王 帅，盛 谦，朱泽奇，等. 地震荷载作用下地下洞室不利地质结构塌落机制研究［J］. 岩土力学，2012，33（10）：2897-2902，2966. Wang Shuai，Sheng Qian，Zhu Zeqi，et al. Study of collapse mechanism of underground caverns with unfavorable geological structures under seismic loading［J］. Rock and Soil Mechanics，2012，33（10）：2897-2902，2966（in Chinese）.

［11］ 郑银河，夏 露，于青春. 考虑岩桥破坏的块体稳定性分析方法［J］. 岩土力学，2013，34（增1）：197-203. Zheng Yinhe，Xia Lu，Yu Qingchun. Stability analysis method of block considering cracking of rock bridge ［J］. Rock and Soil Mechanics，2013，34（Suppl 1）：197-203（in Chinese）.

［12］ 付晓东，盛 谦，张勇慧. 基于矢量和一非连续变形分析的滑坡安全系数计算方法研究［J］. 岩石力学与工程学报，2014，33（增2）：4122-4128. Fu Xiaodong，Sheng Qian，Zhang Yonghui. Research on safety factor computation of landslide based on vector sum-discontinuous deformation analysis method［J］. Chinese Journal of Rock Mechanics and Engineering，2014，33（Suppl 2）：4122-4128（in Chinese）.

［13］ 钟登华，鲁文妍，刘 杰，等. 复杂地质条件下地下洞室曲面块体地震响应分析［J］. 天津大学学报：自然科学与工程技术版，2014，47（6）：471-478. Zhong Denghua，Lu Wenyan，Liu Jie，et al. Surfaceblock identification and seismic response analysis of underground structures under complicated geological conditions［J］. Journal of Tianjin University：Science and Technology，2014，47（6）：471-478（in Chinese）.

［14］ Su Huaizhi，Hu Jiang，Li Jinyou，et al. Deep stability evaluation of high-gravity dam under combining action of powerhouse and dam［J］. International Journal of Geomechanics，2013，13：257-272.

［15］ Jiang Yusong，Su Chao. Mixed finite element method for contact problems of multibody［C］//Proceedings of the 12th Biernial International Conference on Engineering，Science，Construction ，and Operations in Challenging Environments. Honolulu，Hawaii，USA，2010：606-620.

［16］ 费 康，张建伟. ABAQUS在岩土工程中的应用［M］.北京：中国水利水电出版社，2010. Fei Kang，Zhang Jianwei. Application of ABAQUS in Geotechnical Engineering［M］. Beijing：China Water & Power Press，2010（in Chinese）.

［17］ 钟登华，鲁文妍，刘 杰，等. 基于三维地质模型的地下洞室曲面块体分析［J］. 岩石力学与工程学报，2011，30（增2）：3696-3702. Zhong Denghua，Lu Wenyan，Liu Jie，et al. Surfaceblock analysis of underground cavern based on 3D geological modeling［J］. Chinese Journal of Rock Mechanics and Engineering，2011，30（Suppl 2）：3696-3702（in Chinese）.

（责任编辑：樊素英）

Seismic Responses of Underground Cavern Block in Hydraulic Engineering Considering Bond Characteristics

Li Wenqian1，Tong Dawei1，Wang Zhen1，Lu Wenyan1 2，Zhu Xiaobin1，Zhong Denghua1
（1.State Key Laboratory of Hydraulic Engineering Simulation and Safety，Tianjin University，Tianjin 300072，China；2.Nanjing Research Institute of Water Conservancy，Nanjing 210029，China）

Abstract：During the construction and operational process of underground cavern in hydraulic engineering，it is of great significance to study the anti-sliding stability of underground cavern blocks excited by earthquake because the occurrence of block’s instability and damage induced by seismic loads may bring a serious threat to the projects.In view of the bond characteristics on the contact surface between block and rock in the practical projects，an antisliding model for the underground cavern block considering the effect of bond strength was proposed to improve the contact simulation performance of original Coulomb model in ABAQUS and perfect the influencing mechanism of the static behavior and dynamic relative sliding between block and rock.Further，not only the response behavior and stability of the blocks under the static and dynamic excitations was researched，but also the rationality and validity of the anti-sliding model was verified.Under the separate simulation by classical and modified models，the maximum vertical displacements of block in gravity stage are 1.04 cm and 0.05 cm，respectively，and the average vertical displacements and the maximum vertical displacements of block in seismic stage are 2.24 cm，3.91 cm and 2.07 cm，3.38 cm，respectively.The comparison result indicates that the anti-sliding model considering the bond effect has better stability.In addition，the influences of friction coefficient and cohesion in the anti-sliding model on block responses and dynamic characteristics on the contact surface were also discussed comprehensively.

Keywords：rock block；underground cavern；anti-sliding model；dynamic response analysis；ABAQUS

中图分类号：TV312

文献标志码：A

文章编号：0493-2137（2016）04-0369-07

DOI:10.11784/tdxbz201507002

收稿日期：2015-07-01；修回日期：2015-10-22.

基金项目：国家创新研究群体基金资助项目（51321065）；国家自然科学基金重点资助项目（51439005）；天津市自然科学基金青年基金资助项目（13JCQNJC08900）.

作者简介：李文倩（1986— ），女，博士研究生，bxqe@163.com.

通讯作者：鲁文妍，wylu@nhri.cn.