冯国瑞，朱卫兵，李 竹，白锦文，罗泽强

(1.太原理工大学 矿业工程学院，山西 太原 030024；2.山西省绿色采矿工程技术研究中心，山西 太原 030024；3.中国矿业大学 矿业工程学院，江苏 徐州 221116)

我国中西部矿区煤炭资源储量丰富，开采初期限于开采技术和采掘设备水平，采煤方法多为房柱式开采，且呈现出采厚弃薄、采肥丢瘦、采优弃劣和采易弃难特点。随着矿井煤炭资源日益枯竭，残煤复采成为提升矿井资源采出率和实现“精细化采煤”的重要手段，解放受老采空区遗留煤柱影响的整层弃煤逐步成为研究重点。研究层间岩层采动稳定性判别方法、煤柱失稳机理、煤柱群中各煤柱“多米诺骨牌式”链式失稳特征及其动态失稳时序，探索科学开采方法并构建安全复采技术体系是残煤复采的关键所在。

针对该领域，国内外学者对单个煤柱的力学性能开展研究，提出了评价煤柱强度的经验公式和评估煤柱稳定性安全系数的确定方法，促进了对煤柱失稳模式和加固方法领域的研究与实践。事实上，煤柱作为多孔且裂隙发育介质，因胶结弱面、矿井水侵蚀以及采动影响的作用，其承载能力逐步劣化并最终破坏失稳，国内外许多煤矿深受煤柱群下如何安全开采的困扰，一旦煤柱群大面积失稳，易造成下煤层工作面压架和地面大范围塌陷。周子龙等分析了循环载荷作用下煤柱群体系的承载特性及疲劳损伤机理，提出了煤柱群稳定性判别方法及降低煤柱群失稳风险的有效措施。朱德福等等采用重整化群的方法，研究揭示了采动影响下煤柱群应力集中区和释放区的位置及范围。朱卫兵等通过高清摄像机取帧技术，分析了相似模拟中近距离下煤层开采上覆煤柱群的动态失稳过程，再现了煤柱群中各个煤柱破坏形式及其“多米诺骨牌式”链式失稳时序特征。笔者等提出了遗留煤柱群链式失稳的“关键柱”理论，认为遗留煤柱群的整体失稳的根本原因是遗留煤柱群中“关键柱”的失稳，进而引起覆岩载荷转移诱发邻近煤柱联动失稳。白锦文等提出了关键柱柱旁充填岩层控制的技术方法,并揭示了关键柱柱旁充填岩层控制的核心机理。

纵观已有研究成果，有关底板煤柱群安全复采工程技术方面的研究鲜有报道。笔者结合山西元宝湾煤矿浅埋深工作面蹬空复采案例，综合实验室强度测试实验、物理模拟、相似模拟及现场实测方法，研究底板煤柱群应力分布及其破坏特征，提出柱式采空区封闭灌浆充填方法，确定底板蹬空安全复采的临界充填高度，为类似条件下残留煤柱资源的安全开采提供借鉴。

1 工作面概况

山西省元宝湾煤矿6107工作面宽度240 m，推进长度500 m，主采6号煤层，平均厚度3 m，倾角4°～8°，埋深150 m。6107工作面下方为9号煤层房采区，实测底板蹬空区长度普遍为25～50 m，宽度5～15 m，高度5～8 m，9号煤层巷道宽度2.5～5.0 m，高度2.2～4.0 m，6107工作面上方的4号煤层同样为房采区，遗留柱采区宽度7～30 m，高度3～6 m，长度30～60 m，物探发现遗留柱采区顶板较为完整，仅部分区域出现垮落，6107工作面顶底板采掘情况，如图1所示，如何实现6107工作面蹬空阶段的安全复采是亟待解决的重大技术难题。

图1 6107工作面上下柱采区对照Fig.1 Pillar and goaf area overlying and underlying 6107 working face comparison diagram

2 蹬空开采底板煤柱承载及变形特征

2.1 力学模型

根据6107工作面开采现状，将蹬空阶段地质条件简化，如图2所示。其中，为上煤层工作面采高，m；为下煤层工作面采高，m；为底板煤柱宽度，m；为煤柱间距，m。层间岩层厚度取15 m，下煤层残留工作面采高为8 m，各煤柱间距30 m，煤柱宽度为10 m。根据类似开采条件下工作面超前影响范围的研究成果，结合6107工作面采高及岩层柱状，推测该工作面超前影响最大范围为60 m。因此，图2所示工作面超前应力增高区波及底板2个煤柱，即2,3号煤柱处于超前应力增高区，4号处于过渡区，而5～7号煤柱处于原岩应力区。

图2 试验工作面蹬空开采示意Fig.2 Schematic diagram of kick-off mining in test working face

鉴于图2试验工作面蹬空开采模型属于“多跨连续梁”的超静定结构，包含了7个残采房式煤柱和2个半无限边界煤柱。为了便于分析计算，将图2中的蹬空开采覆岩结构，简化为左右边界为固支，中间均匀布设7个弹簧支撑结构，如图3所示。其中，为煤壁处支承应力，MPa；为工作面支承应力峰值，MPa；为原岩应力，MPa。每个弹簧代表图2中对应位置的底板煤柱。

图3 多跨连续梁示意Fig.3 Schematic diagram of multi-span continuous beam

认清底板煤柱的载荷分布规律是评价底板煤柱稳定的先决条件。通过结构力学求解器SMSolver，首先，求解各个弹簧支撑的压缩量；其次，弹簧支撑点作为煤柱的中心支撑点，考虑到煤柱宽度为10 m，间隔1 m提取中心支撑点左、右两侧各5 m的下沉数据。在此基础上，考虑到6107工作面开采前底板煤柱弹性核区占比较大且仍处于稳定状态，将底板煤柱顶界面下沉数据除以煤柱高度获得煤柱压缩应变，该应变乘以弹性模量，获得煤柱顶界面垂直压应力。结果表明，底板各煤柱顶界面垂直压应力分布近似呈线性非均匀分布状态，以6107工作面下伏柱采区1，2，3号煤柱为例，3者顶界面垂直载荷分布形式如图4所示。

图4 底板煤柱顶界面垂直压应力分布形态Fig.4 Stress distribution of floor coal pillar top interface

基于此，建立顶界面载荷非均匀分布条件下煤柱内应力计算力学模型，依据各煤柱内应力计算结果，结合第四强度理论，深入探究各煤柱内应力分布规律。如图5所示，其中，为底板煤柱顶界面左端应力，MPa；为底板煤柱顶界面右端应力，MPa；为煤柱内任一点正应力，MPa；为煤柱内任一点切应力，MPa。

图5 煤柱顶界面应力分布示意Fig.5 Stress distribution diagram of coal pillar top interface

设图5中煤柱内Airy应力函数为

=+++

(1)

可得各应力分量为

(2)

式中，，，，为待定系数。

根据上述边界条件，可得待定系数，，，计算结果为

(3)

2.2 内应力分布规律

房采煤柱因上覆顶板岩层不均匀下沉而呈现偏心压缩现象，斜切变形是导致煤柱发生屈服及破坏的主要形式。依据第四强度理论，由主应力求解煤柱内部各点Mises等效应力可更为准确地表征煤柱是否发生塑性破坏，相关参数取值中超前应力集中系数为2，工作面埋深为150 m，地层岩石容重取25 000 N/m，故=3.75 MPa，=7.00 MPa；峰值距煤壁距离为10 m。煤柱弹性模量2 GPa，层间岩层弹性模量25 GPa，层间岩梁厚度取平均值15 m且为单位宽度，煤柱宽度10 m，煤柱高度8 m，以2，3号煤柱和5，6号煤柱为例，其内应力分布形态如图6所示。

(4)

式中，为等效应力，MPa；，，分别为第1、第2、第3主应力，MPa。

图6 底板各煤柱Mises等效应力分布Fig.6 Distribution of Mises equivalent stress of coal pillars in floor and the division of dangerous areas

如图6所示，在超前支承压力作用下，2，3，5，6号各煤柱的左上角和右上角均出现不同程度的应力集中，但2，3号煤柱应力集中差异程度更大，2号煤柱的右上角以及3号煤柱的左上角是易发生破坏的危险区域，其斜切破坏趋势更为明显；而5，6号煤柱左上角和右上角应力集中程度基本一致，表明5，6号煤柱趋近于均匀压缩状态，煤柱内部各点等效应力基本保持一致，处于尚未破坏的状态。

3 煤柱破坏形式及柱群失稳扩展规律

3.1 相似材料物理力学参数测试

在掌握底板煤柱因顶板非均匀下沉和载荷非均匀分布导致斜切破坏的基础上，通过相似模拟研究，进一步探究底板煤柱群的动态失稳次序。根据试验工作面覆岩关键层判别结果，将全部地层简化为关键层1，2，3、煤层及软岩5种岩层。针对这5种岩性进行实验，获得符合要求的实验材料配比。制作的标准试件如图7(a)所示。不同配比对应的应力应变曲线如图7所示，根据模拟实验1∶160的应力相似比要求，配比1最为合适。

3.2 物理模拟方案

通过物理模拟研究近距离煤层开采上、下煤柱采区条件下层间岩层运动特征、柱群破坏失稳先后次序及其扩展规律，模型长度130 cm，宽度12 cm，高度63 cm。基于地层柱状及相似材料强度测试结果，各岩层厚度及材料配比见表1。模型的几何相似比1∶100，密度相似比1.0∶1.6，应力相似比1∶160。模型顶部采用铁块加载的方式对模型进行载荷补偿，用以模拟150 m埋深。

图7 实验材料配比试件参数测定Fig.7 Experimental material proportioning specimen and parameter measurement

表1 相似模拟材料基本物理力学参数

6107工作面煤层开采前，先采4号煤层，煤柱宽度为5 cm，间距为5 cm；而后开采9号煤层，煤柱宽度为2 cm，间距为7 cm，如图8所示。6号煤层开采时，两侧各留宽度为10 cm的煤柱，工作面每间隔10 min向前开采5 cm，模型开采过程中通过高速摄影技术进行实时记录，以准确捕捉和再现煤柱破裂、失稳及岩层破断运动的全过程。

图8 物理模拟实验方案Fig.8 Scheme diagram of physical simulation experiment

3.3 底板煤柱破坏形式及失稳时序

6107工作面回采过程中，物理相似模拟实验很好地揭示了底板煤柱的破坏形式及底板煤柱群中“关键柱”的位置，对文献[18]提出的关键柱理论起到重要支撑。当6107工作面推进距小于30 cm时，顶、底板煤柱群均能够保持稳定状态。当6107工作面推进60 cm时，基本顶悬空段中部产生张拉裂隙，此时煤柱出现较大变形，但仍能保持完好和稳定。通过高速摄像机逐帧分析，可以发现顶板岩层产生了明显的非均匀下沉特征，且层间岩层的非均匀下沉对9号煤层柱采区中遗留煤柱产生了显著的斜切作用，使底板煤柱剥落体产生旋转下落的现象，如图9所示。6107工作面开采过程中，上、下煤柱群失稳及覆岩破断运动过程如图10所示。

结果表明，煤柱裂隙发育、剥落至整个煤柱群失稳可分为5个阶段：

(1)煤柱均匀承载压密阶段。中层遗煤开采前，层间岩层基本处于未采动状态，呈现整体弯曲下沉的现象，各煤柱顶界面载荷近似于均匀分布，如图10(a)所示。

(2)煤柱非均匀承载阶段。伴随着中层遗煤工作面向前推进，底板煤柱逐步过渡为非均匀承载状态，且非均匀承载随工作面推进愈发显著，如图10(b)所示。

(3)关键柱破坏失稳阶段。当工作面推进至60 cm时，超前支承压力影响下底板煤柱非均匀承载达最大化，此时距工作面最近的底板煤柱率先失稳，如图10(c)，(d)所示。

(4)失稳载荷双向传递阶段。关键柱失稳后，上覆载荷快速向前、后邻近煤柱传递，诱发邻近煤柱的破坏及失稳，覆岩载荷双向传递特征如图10(e)所示。

(5)煤柱群整体失稳阶段。伴随着覆岩载荷急速地在底板煤柱群中双向传递，整个煤柱群形成以“关键柱”为中心“多米诺骨牌式”链式失稳特征，如图10(f)所示。

图9 物理模拟中破坏失稳煤柱剥落体运动特征Fig.9 Movement characteristics of failure pillar fractured bodies in physical simulation

图10 中层遗煤顶底板柱采区煤柱动态失稳过程Fig.10 Dynamic collapse process of coal pillars in middle residual coal roof and floor

4 充填率对煤柱稳定性的影响规律

上述研究结果表明，6107工作面蹬空开采过程中，若不对底板煤柱进行加固，势必导致工作面及超前煤岩体塌陷事故。为此，提出对底板柱采区进行注浆充填，注浆固结体对底板煤柱形成侧向夹持作用，抑制底板煤柱塑性区扩展并提高煤柱承载能力，以此确保6107工作面安全回采。

4.1 充填体对煤柱弹塑性扩展的抑制效应

根据文献[17]统计结果，当煤柱弹性核区占比达35%～40%时，可认为煤柱处于稳定状态。为确保6107工作面安全复采，以40%作为煤柱临界失稳判定标准。借助9号煤层煤柱煤样应力应变测试结果对模拟参数进行标定，如图11(a)，(b)所示。以4号煤柱为例，充填体高度对底板煤柱弹塑性扩展的影响规律如图11(c)～(h)所示。

图11 煤体模拟参数标定及充填体高度对煤柱塑性区扩展的抑制效应Fig.11 Inhibition effect on plastic zone expansion due to calibration of coal simulation parameters and backfilling body height

由图11可知，在9号煤层柱采区未充填时，4号底板煤柱塑性区完全贯通，其破坏失稳势必导致整个煤柱群失稳及工作面人员、采掘设备及超前煤岩体的大面积塌陷。随着充填体高度的增加，受充填体的侧向夹持作用，底板煤柱弹性核区占比快速提升。当充填体高度达到6 m时，煤柱弹性核区占比达60%。表明随着充填体高度的增加，充填体对煤柱塑性区扩展形成强力的抑制效应，提高了底板煤柱的承载能力。

4.2 临界充填高度及蹬空开采安全评价

当充填体高度为6 m，即不接顶高度为2 m时，可确保6107工作面安全复采。为此，通过相似模拟研究不接顶高度2 m时底板煤柱稳定性，以佐证数值模拟结果。为确保充填效果，充填材料被制作成类牙膏稠度状态，以避免浆液漏失及浸泡弱化底板煤柱，在柱式采空区前后安装挡板，随之在柱采区中进行了注浆充填，待充填材料固结后撤去挡板，如图12所示。

结果表明，虽然仍有2 m不接顶高度，但下煤层煤柱与充填体胶结良好，充填部分的煤柱从单轴受力变成二轴受力，减小了底板煤柱承载高度，极大提升了煤柱的稳定性。这意味着当充填体高度为6 m时，可有效确保6107工作面中层遗煤的安全复采。

图12 下煤层柱采区充填后岩层运动及煤柱稳定性Fig.12 Strata movement and pillar stability after backfilling of underlying pillar mining area

5 蹬空段充填效果实测及评价

5.1 充填效果观测

基于前述研究成果，在开采前对9号煤层柱采区进行地面钻孔注浆充填，注浆完成后钻孔并未封闭，利用钻孔电视设备在未封闭的5，13号钻孔内对9煤层柱采区的充填情况进行窥视，钻孔窥视方案如图13所示。

钻孔窥视发现：5号钻孔电视下探至-161 m到达9号煤层柱采区，如图14(a)所示，继续下探至-163 m见充填区积水面，如图14(b)所示，此处柱采区充填未接顶高度约为2 m，充填效果较好；13号钻孔电视下探至-164 m见钻孔坍塌位置，如图14(c)所示，继续下探至-168 m见充填区积水面，如图14(d)所示，表明此处尚存在高度约4 m的未充实区域，推测13号孔附近柱采区充填效果一般，近4 m的充填未接顶空洞将对6107工作面的安全回采造成影响。

5.2 柱采区充填未接顶防治对策

在6107工作面推进至5号钻孔对应位置时，鉴于充填体高度已达6 m，未接顶高度仅为2 m，底板煤柱承载能力较强，6107工作面未采取任何措施得以安全推过。而在6107工作面推进至13号钻孔对应位置时，在工作面前方发现由于9号煤层柱采区充填不接顶及顶板垮落形成的空洞区域，其长度约10 m，宽度5～6 m，高度1～3 m，空洞垮塌范围贯穿9号煤层与6号煤层的层间岩层，且6号煤层也发生局部垮落，迫使6107工作面短暂停止回采。针对上述情况，采用工作面采煤机割落的煤块，再填上少量废旧工字钢对工作面前方空洞区域进行充填，最终确保工作面液压支架安全移架，安全采出后续煤炭资源。

图13 窥视孔位置及钻孔窥视方案Fig.13 Borehole location and borehole peeping scheme diagram

图14 钻孔窥视结果Fig.14 Borehole peeping result diagram

6 结 论

(1)揭示了底板煤柱偏心压缩条件下破坏形式及运动规律，9号煤层柱采区若在未充填情况下，蹬空开采阶段的底板煤柱将成为承载地基，导致其内应力急剧增大，且由于层间岩层的非均匀下沉，易导致底板煤柱产生偏心压缩现象，使得煤柱剥落体呈现出旋转下沉的运动特征，最终会造成底板煤柱发生斜切破坏。

(2)揭示了底板煤柱群破坏失稳次序，底板煤柱群失稳次序呈现出“关键柱”率先失稳，覆岩载荷双向传递，煤柱群“多米诺骨牌式”双向链式失稳特征。整个底板房柱式采空区中，则呈现出以“关键柱”为辐射中心，由内向外多向辐射的波状失稳现象。

(3)研究确定6107工作面蹬空段安全复采底板煤柱群临界充填高度为6 m，并得到了现场实测的验证。5号钻孔区域充填体达到6 m，不接顶高度为2 m，工作面实现了安全回采；13号钻孔区域充填不接顶高度为4 m，工作面蹬空开采存在塌陷风险。现场采用局部强化充填底板空洞的方法，确保了6107工作面的安全回采。

[1] 冯国瑞,张玉江,戚庭野,等. 中国遗煤开采现状及研究进展[J].煤炭学报,2020,45(1):151-159.

FENG Guorui,ZHANG Yujiang,QI Tingye,et al. Status and research progress for residual coal mining in China[J]. Journal of China Coal Society,2020,45(1):151-159.

[2] 冯国瑞,侯水云,梁春豪,等. 复杂条件下遗煤开采岩层控制理论与关键技术研究[J]. 煤炭科学技术,2020,48(1):144-149.

FENG Guorui,HOU Shuiyun,LIANG Chunhao,et al. Basic theories and key technologies of rock strata control for residual coal resources mining under complex conditions[J]. Coal Science and Technology,2020,48(1):144-149.

[3] 周小建,冯国瑞,白锦文. “上垮落-下刀柱”复合残采区中层弃煤开采对下部采空区煤柱稳定性的分析[J]. 煤炭技术,2017,36(6):47-49.

ZHOU Xiaojian,FENG Guorui,BAI Jinwen. Coal-pillar stability analysis of “collapse on-column under knife” complex residues abandoned coal mining area be mined in middle[J]. Coal Technology,2017,36(6):47-49.

[4] 冯国瑞,张绪言,李建军,等. 刀柱采空区上方遗弃煤层上行开采可行性判定[J]. 煤炭学报,2009,34(6):726-730.

FENG Guorui,ZHANG Xuyan,LI Jianjun,et al. Feasibility on the upward mining of the left-over coal above goaf with pillar supporting method [J]. Journal of China Coal Society,2009,34(6):726-730.

[5] SINGH R,SINGH A K,MAITI J,et al. An observational approach for assessment of dynamic loading during underground coal pillar extraction[J]. International Journal of Rock Mechanics and Mining Sciences,2011,48(5):794-804.

[6] POULSEN B A,SHEN B. Subsidence risk assessment of decommissioned bord-and-pillar collieries[J]. International Journal of Rock Mechanics and Mining Sciences,2013,60:312-320.

[7] 祁和刚,于健浩. 深部高应力区段煤柱留设合理性及综合卸荷技术[J]. 煤炭学报,2018,43(12):3257-3264.

QI Hegang,YU Jianhao. Rationality and comprehensive unloading technology of deep high stress section coal pillars[J]. Journal of China Coal Society,2018,43(12):3257-3264.

[8] 谭毅,郭文兵,白二虎,等. 条带式Wongawilli煤柱特征及作用机理分析[J]. 煤炭学报,2019,44(4):1003-1010.

TAN Yi,GUO Wenbing,BAI Erhu,et al. Coal pillar characteristics and its mechanism of strip wongawilli mining[J]. Journal of China Coal Society,2019,44(4):1003-1010.

[9] 张杰,王斌. 浅埋间隔采空区隔离煤柱稳定性及覆岩失稳特征研究[J]. 采矿与安全工程学报,2020,37(5):936-942.

ZHANG Jie,WANG Bin. Stability of isolated coal pillar and overburden instability in shallow-buried interval gob[J]. Journal of Mining and Safety Engineering,2020,37(5):936-942.

[10] 张明,成云海,王磊,等. 浅埋复采工作面厚硬岩层-煤柱结构模型及其稳定性研究[J]. 岩石力学与工程学报,2019,38(1):87-100.

ZHANG Ming,CHENG Yunhai,WANG Lei,et al. Struct-ure model and stability research of thick hard strata-coal pillar in shallow-buried re-mined panels[J]. Chinese Journal of Rock Mechanics and Engineering,2019,38(1):87-100.

[11] 徐敬民,朱卫兵,鞠金峰. 浅埋房采区下近距离煤层开采动载矿压机理[J]. 煤炭学报,2017,42(2):500-509.

XU Jingmin,ZHU Weibing,JU Jinfeng. Mechanism of dynamic mine pressure occurring below adjacent upper chamber mining goaf with shallow cover depth[J]. Journal of China Coal Society,2017,42(2):500-509.

[12] 顾大钊,颜永国,张勇,等. 煤矿地下水库煤柱动力响应与稳定性分析[J]. 煤炭学报,2016,41(7):1589-1597.

GU Dazhao,YAN Yongguo,ZHANG Yong,et al. Experimental study and numerical simulation for dynamic response of coal pillars in coal mine underground reservoir[J]. Journal of China Coal Society,2016,41(7):1589-1597.

[13] ZHOU Zilong,WANG Haiquan,CAI Xin,et al. Bearing characteristics and fatigue damage mechanism of multi-pillar system subjected to different cyclic loads[J]. Journal of Central South University,2020,27(2):542-553.

[14] CHEN Lu,ZHOU Zilong,ZANG Chuanwei,et al. Failure pattern of large-scale goaf collapse and a controlled roof caving method used in gypsum mine[J]. Geomechanics and Engineering,2019,18(4):449-457.

[15] ZHU Defu,TU Shihao,TU Hongsheng,et al. Mechanisms of support failure and prevention measures under double-layer room mining gobs-A case study:Shigetai coal mine[J]. International Journal of Mining Science and Technology,2019,29(6):955-962.

[16] ZHU Weibing,CHEN Lu,ZHOU Zilong,et al. Failure propag-ation of pillars and roof in a room and pillar mine induced by longwall mining in the lower seam[J]. Rock Mechanics and Rock Engineering,2019,52(4):1193-1209.

[17] 朱卫兵,许家林,陈璐,等. 浅埋近距离煤层开采房式煤柱群动态失稳致灾机制[J]. 煤炭学报,2019,44(2):358-366.

ZHU Weibing,XU Jialin,CHEN Lu,et al. Mechanism of disaster induced by dynamic instability of coal pillar group in room-and-pillar mining of shallow and close coal seams[J]. Journal of China Coal Society,2019,44(2):358-366.

[18] 冯国瑞,白锦文,史旭东,等. 遗留煤柱群链式失稳的关键柱理论及其应用展望[J]. 煤炭学报,2021,46(1):164-179.

FENG Guorui,BAI Jinwen,SHI Xudong,et al. Key pillar theory in the chain failure of residual coal pillars and its application prospect[J]. Journal of China Coal Society,2021,46(1):164-179.

[19] 白锦文,崔博强,戚庭野,等. 关键柱柱旁充填岩层控制基础理论[J]. 煤炭学报,2021,46(2):424-438.

BAI Jinwen,CUI Boqiang,QI Tingye,et al. Fundamental theory for rock strata control of key pillar-side backfilling[J]. Journal of China Coal Society,2021,46(2):424-438.

[20] 潘岳,顾士坦,李文帅. 煤层弹性、硬化和软化区对顶板弯矩特性影响分析[J]. 岩石力学与工程学报,2016,35(2):3846-3857.

PAN Yue,GU Shitan,LI Wenshuai. Analysis on bending moment property of hard roof influenced by elastic,hardening and softening zone of coal seam[J]. Chinese Journal of Rock Mechanics and Engineering,2016,35(2):3846-3857.