佛顶宫矿坑边坡三维静动力有限元模拟分析研究

张菊连，梁志荣，李 伟

(上海申元岩土工程有限公司，上海 200040)

佛顶宫矿坑边坡三维静动力有限元模拟分析研究

张菊连，梁志荣，李 伟

(上海申元岩土工程有限公司，上海 200040)

对国内首个既有矿坑地下空间开发项目的佛顶宫边坡进行了三维有限元研究，分别对静力和地震这两种工况进行了模拟。结果表明：削坡或锚索(杆)加固后，最大位移均发生在坡顶、坡面填土、残积土和强风化岩处；最大应力均发生在坡内，最大剪应力发生在各坡脚和岩土交界面处；加速度响应随着边坡高度方向呈放大趋势，表层土体加速度响应明显放大；极限状态下，广义塑性应变区集中在坡顶、坡面处。静力工况下，锚索(杆)加固工况相对于无支护工况，边坡最大位移减小11.9%～28.4%，最大应力减小5.6%～67.1%，安全系数提高了34.5%，说明锚索(杆)加固效果显著。地震相对静力支护工况，最大位移增大3.03～19.73倍，主应力增大6.5%～146.7%，最大剪应力增大168.4%～474.6%，安全系数减小0.6%～23.1%。大震相对于小震支护工况下，坡面加速度峰值及放大倍数分别增大15.04～27.00倍、1.44～1.93倍。对矿坑南北、东西方向施加地震计算的结果显示：前者动位移、动应力均较后者大，安全系数较后者小，说明矿坑长轴方向施加地震较短轴方向不利。

隧道工程；既有深坑；地震工况；位移；应力；地震加速度响应

0 引 言

边坡稳定性方面的研究有200多年的历史，计算方法在理论、数值分析方面都有非常丰富的成果[1-7]。地震工况下的边坡动力响应包括加速度、速度、位移、动应力和动应变响应等。比如：毕忠伟等[3]对其进行了分析研究，但多数集中在简单的均质假想边坡，研究结论很难代表现场实际情况。

地震荷载作用下边坡稳定性研究常采用的分析方法有：拟静力法[4]、滑块分析法(Newmark法)、有限元方法。拟静力法和Newmark法的假设与实际相差较大，难以真实地反映地震荷载作用下岩体边坡的动力学行为；数值模拟法能够考虑振幅、频谱、持时等地震荷载的特性，能够真实模拟边坡在地震荷载作用过程中的动力特征和破坏机理。

反应谱法和时程分析方法则是建筑结构抗震设计中常用的两种方法。对一些重要工程和特殊结构，为了解其在整个地震过程中的动力行为，需要使用时程分析方法。时程分析法直接从运动微分方程出发，用数值方法逐步积分，可以获得结构在整个地震作用中的运动状态变化情况。时程分析法包括确定性动力分析与非确定性动力分析两大类[5-7]。

利用有限元法对无限域或半无限域问题进行数值模拟时，通常采取的办法是人为地截取一定宽度的地基范围，将无限区域变成有限区域，再将区域离散为有限个单元，并在截取边界上施加人为的约束来近似处理。由于截取的模型往往具有更大的范围和更多的自由度，离散范围很大，划分的单元及结点数成倍增加，使计算工作量增大，若涉及到非线性问题求解将更加困难。无限元是几何上趋于无穷的单元，是对有限元求解无限域问题的有效补充，有限元无限元联合法成为解决该类问题的有效方法。有限元/无限元法，在有限元研究法的基础上引进无限元理论，形成有限元/无限元联合计算。无限元的引入，能合理地反映真实边界条件，使截取的岩土体范围减少，既提高了计算精度，也减少了计算工作量。由于无限元多布置在有限区域外围，在弹塑性等非线性分析中，无限元按弹性单元处理即可，可免去许多烦琐的重复计算。

笔者结合国内第一个既有矿坑地下空间开发项目—南京牛首山佛顶宫项目，采用有限元-无限元单元、时程分析结合强度折减法，对复杂岩土体矿坑边坡进行了三维有限元模拟。研究结果有助于进一步揭示边坡在地震作用下的破坏机理。

1 项目概况

南京牛首山文化旅游区一期工程位于牛首山遗址公园核心区域，其场地原貌见图1。

图1 矿坑全景Fig.1 Mine panorama

根据文献[8]，南京市抗震设防烈度为7度，设计地震加速度峰值为0.10 g，设计地震分组为第一组，场地类别为二类，特征周期为0.35 s。支护抗震设计按小震(50年超越概率63%)设计，变形控制按大震(50年超越概率2%)验算。拟建场地属对建筑抗震不利地段，部分建筑位于滑坡顶部的危险地段。场地内从上至下分布有：①杂填土、②残坡积土、③1强风化凝灰岩、③2中风化凝灰岩[9]、④1强风化蚀变安山质凝灰岩、④2中风化蚀变安山质凝灰岩；局部为：③2a中风化(破碎)凝灰岩、④2a中风化(破碎)蚀变安山质凝灰岩。

2 计算模型和参数

2.1 计算模型

根据勘察资料[9]、佛顶宫削坡加固图纸及建筑方案等相关资料建立边坡三维有限元分析模型。边坡尺寸选取如下：模型选取宽度665 m，长786 m，高85 m(图2)。

图2 削坡状态下三维模型几何模型Fig.2 Three dimensional geometric model of cutting slope

模型中的土层和风化岩层均采用的是摩尔-库伦模型。锚杆、锚索采用二结点的梁单元模拟。对梁单元施加预压力从而产生预应力效果。计算中考虑了锚杆、锚索的和轴向拉压作用及抗弯作用，忽略了锚杆、锚索与岩体接触面上的剪切滑移作用。在静力计算中，将重力作为主要作用力施加在锚固坡体上。在动力计算中，预先将锚杆锚索单元设为空单元，对模型施加重力求出边坡内的初始应力分布，然后激活空单元，导入初始应力场，施加地震动荷载，求出在地震动荷载作用下基坑的基本受力变形特征以及相应的安全系数值。实际施工后，由于各种原因带来的预应力损失，砂浆的剪切滑移表现在预应力松弛上，在计算中取设计预应力值的80%作为锚杆锚索的预应力值。削坡锚索(杆)加固后的模型见图3。

2.2 物理、力学参数

根据文献[9]的建议值和文献[10]的规定，确定了材料参数，如表1。

图3 锚固状态下三维模型网格划分Fig.3 Three dimensional mesh model of reinforced slope表1 有限元强度折减法计算采用的材料参数Table 1 The material parameters for strength reduction FEM calculation

参数岩土层①杂填土②残积土③1强风化凝灰岩③2中风化凝灰岩④1强风化安山质凝灰岩④2中风化安山质凝灰岩③2a中风化(破碎)凝灰岩④2a中风化(破碎)蚀变安山质凝灰岩容重γ/(kN·m-3)19.319.920.424.920.524.823.723.1弹性模量E×102/MPa167171235250158399932832672泊松比γ0.430.430.330.250.340.240.260.25动弹性模量Ed×102/MPa117117170039991838691126293843动泊松比γd0.420.420.400.370.390.380.400.40抗剪强度黏聚力c/kPa010(18)(600)(18)(500)(42)(42)内摩擦角φ/(°)3230(35)(45)(35)(44)(37)(37)基底摩擦因数μ0.200.280.400.500.400.50地基承载力特征值f/kPa9015026040002203800岩土体与锚固体黏结强度特征值frb/kPa100260100300180200

注：括号中的参数为经验值。

2.3 边界条件

静力计算时，模型除边界采用无限元外，内部均采用四节点单元进行划分。节点总数为67 947个，单元总数为362 029个。其中坑周边单元分布密集，密集范围内的单元最大尺寸不超过2.0×2.0 m2，在远离坡体处，单元稀疏分布，网格尺寸是5～20 m。坐标取向沿高度方向为z正方向，沿边坡方向为x，y正方向，长度单位为m。模型中所有单元均采用摩尔-库伦屈服准则。边坡边界采用普通边界，底部为固定边界。

动力计算时，除左右两边边界为无限元边界外，其余网格均与静力计算相同。底部边界将x，y方向放开，变为自由向，z方向仍然固定约束，在模型底部x，y方向输入南京人工合成波作为地震动荷载的输入。

2.4 工况和加载情况

数值计算主要考虑以下几种工况：①天然工况；②开挖过程及开挖后无支护状态工况；③开挖后支护状态下自然工况；④开挖后无支护状态下地震荷载工况；⑤支护状态下地震荷载工况。

其中，输入的地震动为南京地区人工合成波[11]，分别代表大震状态，50年超越概率为2%，输入最大加速度值为180 gal(gal=10-2m/s2)；中震状态，50年超越概率为10%，输入最大加速度值为103 gal；小震状态，50年超越概率为63%，输入最大加速度值为34 gal(图4)。在地震工况下，对x(东西方向)和y(南北方向)分别施加荷载进行计算。

图4 南京牛首山人工合成地震波Fig.4 Artificial seismic waves of Nanjing Niushou Mountain

3 静力有限元分析

3.1 位移场分布

削坡及锚索(杆)加固后，最大位移均发生在坡顶和坡面为填土、残积土层和强风化岩中。边坡最大水平向位移(U1)减小11.9%、沉降(U2)增大13.8%、坡顶与坡底的最大水平相对位移(U3)减小28.4%。说明锚索(杆)加固后，边坡的水平位移、坡顶与坡底水平相对位移减小显著，这体现了锚索(杆)提供水平抗力、限制水平位移的特点。

3.2 应力场分布

对削坡及锚索(杆)加固后的边坡应力(Mises应力σm、最大主应力σ1、最小主应力σ3及剪应力τ)进行了分析。受自重作用，坡面处应力较小，最大应力均发生在坡内如图5(a)，最大剪应力发生在各坡脚和岩土交界面处如图5(b)。锚索加固后，应力分布规律没有显著的变化，但数值减小5.6%～67.1%，说明锚索(杆)加固后边坡的应力得到改善。

图5 静力锚索加固工况下的云图Fig.5 Cloud chart under anchor cable reinforcement condition

3.3 稳定系数

将模型材料参数进行折减，计算不收敛时或位移发生突变时，边坡整体处于极限稳定状态，此时的折减系数即为安全系数。加固前后安全系数分别为1.19,1.60，安全系数提高了34.5%，加固效果显著如图6。

图6 边坡等效塑性应变云图Fig.6 Equivalent plastic strain cloud chart

极限稳定状态下，从等效塑性应变图中可以看到：支护前，边坡的潜在的破坏位置为矿坑的正南侧、正北侧和东北侧区域，具体潜在滑动坡体大小如图6(a)。支护后，其潜在滑移面位置仍为正南、正北侧区域，具体潜在滑动坡体大小如图6(b)。与支护前相比，塑性变形区明显缩小，滑动块体大幅减小，支护加固效果明显。

4 动力有限元分析

4.1 位移响应

跟踪边坡的位移变化可以发现：位移不断随地震迟时发生变化。x方向(东西方向)施加地震荷载时，在10 s的时刻，位移达到最大值；y方向(南北方向)施加地震荷载时，在15 s的时刻，位移达到最大值。位移场分布显示最大的水平位移、沉降与静力下的分布情况一致，均发生在坡顶或坡面的填土、残积土及强风化岩层中。将水平位移(U1)、坡顶与坡底的最大水平相对位移(U3)分别绘制于图7。

图7 支护工况下曲线Fig.7 Chart of slops with reinforcement

由图7可以看出：支护状态下，位移随地震幅值的增大而增大。小震工况下U1增大3.03～5.45倍，U3增大4.17～5.93倍；中震工况下U1增大8.29～8.72倍，U3增大9.00～9.34倍；大震工况下U1增大9.17～19.73倍，U3增大11.10～19.72倍，加速度峰值越大，位移响应越显著；y方向地震引起的最大位移量比x方向地震引起的大，说明矿坑长轴方向的动位移响应大于短轴方向。

4.2 应力分布

动力响应下的最大应力基本上发生在坡内部如图8(a)，最大剪应力τmax发生在坡脚如图8(b)，与静力情况下应力分布相似。分别将静力、地震工况下的应力大小反映在图9中，可以发现主应力和剪应力随地震加速度峰值的增大而增大，主应力相对于静力工况分别增加6.5%～47.2%(x方向地震)、57.8%～146.7%(y方向地震)，剪应力相对静力工况分别增加168.4%～328.4%(x方向地震)、270.0%～474.6%(y方向地震)，剪应力相对增大的更明显一些，反映了锚索(杆)起了显著的加固效果；y方向地震下的应力比x方向地下的应力明显要大，说明矿坑长轴的动应力响应比短轴方向要强烈。

图8 应力云图Fig.8 Stress distribution of pit slope

图9 各种地震工况下的应力变化情况Fig.9 Stress changes under various seismic conditions

4.3 加速度响应

加速度响应随着边坡高度方向呈放大趋势，表层土体加速度响应明显放大。最大加速度响应均发生在坡顶边缘处。将坡顶边缘、坡面加速度响应峰值amax及放大倍数Ta分别绘于图10。

图10 加速度响应曲线Fig.10 Maximum acceleration curve

图10(a)显示坡顶边缘处加速度与支护与否没有负的相关关系，也不随地震加速度峰值的增大而增大，可能原因是坡顶边缘处的土体没有受到约束的作用。图10(b)、图10(c)表明坡面加速度峰值及放大倍数在支护后小震工况下均减小20 %，并随地震加速度峰值的增大而增大，说明坡面加固效果良好。相对小震，中震情况下加速度峰值增大5.00～14.82倍，放大倍数增大0.75～1.86倍；大震情况下加速度峰值增大15.04～27.00倍，放大倍数增大1.44～1.93倍。

4.4 稳定系数

将模型材料参数进行折减，然后输入地震波进行计算。当结果中某时刻出现塑性区贯通或者局部贯通状态，或者计算不收敛时，此刻是极限稳定状态，此时的折减系数即为安全系数。3种地震工况下的安全系数(表2)，能满足小震不坏、中震可修、大震不倒的设计思想。y方向地震下的安全系数较x方向小，说明y方向施加地震对矿坑安全性更加不利。

表2 边坡安全系数Table 2 Slope safety coefficients

塑性应变区发生在坡顶和坡面为填土、残积土层和强风化岩中(图11)，塑性区向坡体深部转移，说明锚索加固区的坡体的安全系数均大于最小安全系数，加固后坡体得到较好的稳固。塑性区范围随地震强度的增大而增大：小震情况下，塑性区仅发生在矿坑的北侧顶部；中震情况下，北侧塑性区范围增大，西侧有新的塑性开展区；大震情况下，原有塑性区范围进一步扩大，矿坑的多处出现新的塑性开展区。相对静力无支护工况，安全系数在小震工况下减小了10.1%～15.1%；相对静力有支护工况，安全系数在小震工况下减小0.6%～4.4%、中震工况下减小3.8%～5.6%、大震工况下减小了20.0%～23.1%，地震对边坡的稳定性影响显著。

图11 边坡等效塑性应变云图(x方向)Fig.11 Equivalent plastic strain cloud chart(x direction)

5 结 论

采用三维有限元-无限元法模拟，分析了国内首个矿坑边坡在静力和地震工况下的位移场、应力场、加速度响应和极限状态下的安全系数，得到如下结论。

1)静力和地震工况下，削坡或锚索(杆)加固后，最大位移均发生在坡顶和坡面为填土、残积土层和强风化岩中。加固后的静力工况，边坡最大位移减小11.9%～28.4%，说明加固效果显著。相对于静力支护工况，小震支护工况下位移增大3.03～5.93倍，中震支护工况下位移增大8.29～9.34倍，大震支护工况下位移增大9.17～19.73倍，位移增大的幅度比较大，设计应考虑以大震工况下的变形作为控制变形。

2)静力和地震工况下，最大应力均发生在坡内，最大剪应力发生在各坡脚和岩土交界面处。静力工况下，锚索加固后，应力减小5.6%～67.1%，说明锚索(杆)加固后边坡的应力得到改善。地震工况下，相对于静力工况，主应力增大6.5%～146.7%，剪应力增大168.4%～474.6%，边坡应力随地震加速度峰值的增大而增大。

3)加速度响应随着边坡高度方向呈放大趋势，表层土体加速度响应明显放大，最大加速度响应均发生在坡顶边缘处。坡面加速度峰值及放大倍数在支护后小震工况下均减小20%，说明支护对坡面岩土体起到良好的加固效果。

4)极限稳定状态下，广义塑性应变区集中在坡顶和坡面为填土、残积土层和强风化岩中。静力加固后，安全系数提高了34.5%，加固效果显著。相对静力支护工况下的安全系数，小震工况减少0.6 %～4.4 %、中震工况减小3.8%～5.6%、大震工况减小了20.0%～23.1%，地震对边坡的稳定性影响显著。

5)通过y方向和x方向动力响应的比较发现，矿坑长轴方向较短轴方向的动位移、动应力大，而安全系数小，说明矿坑长轴方向的地震稳定性较短轴方向较差。

[1] 陈寅春,胡文亮.基于有限元强度折减法多级边坡稳定性研究[J].重庆交通大学学报(自然科学版),2014,33(6):100-103. CHEN Yinchun, HU Wenliang. Stability study of multi-step slope based on strength reduction method[J].JournalofChongqingJiaotongUniversity(ScienceEdition),2014,33(6):100-103.

[2] 赵彪,陈洪凯.边坡地震稳定性完全动力分析法[J].重庆交通大学学报(自然科学版),2011,30(增刊1):580-583. ZHAO Biao,CHEN Hongkai. Fully dynamic analysis method for seismic stability of slope[J].JournalofChongqingJiaotongUniversity(ScienceEdition), 2011,30(Sup1):580-583.

[3] 毕忠伟,张明,金峰,等.地震作用下边坡的动态响应规律研究[J].岩土力学,2009,30(增刊1):180-183. BI Zhongwei, ZHANG Ming, JIN Feng, et al. Dynamic response of slopes under earthquakes[J].RockandSoilMechanics, 2009, 30 (Sup1): 180-183.

[4] 陈林杰,郑晓卫.基于有限元强度折减法的地震区三维边坡稳定性分析[J]. 重庆交通大学学报(自然科学版), 2013,32(3):415-418. CHEN Linjie, ZHENG Xiaowei. Three-dimensional slope stability of earthquake zone based on strength reduction FEM[J].JournalofChongoingJiaotongUniversity(NaturalScience),2013,32(3):415-418.

[5] 王思敬,薛守义.岩体边坡楔形体动力学分析[J].地质科学,1992(2):177-182. WANG Sijing, XUE Shouyi. Dynamic analysis of wedge sliding on rock slopes[J].ScientiaGeologicaSinica, 1992(2): 177-182.

[6] 张菊明,王思敬.层状边坡岩体滑动稳定的三维动力学分析[J].工程地质学报,1994,2(3):1-12. ZHANG Juming, WANG Sijing. 3-D dynamic analysis of rock sliding on the layered rock slope[J].JournalofEngineeringGeology, 1994, 2(3): 1-12.

[7] 薛守义,王思敬,刘建中.块状岩体边坡地震滑动位移分析[J].工程地质学报,1997,5(2):131-136. XUE Shouyi, WANG Sijing, LIU Jianzhong. Analysis of sliding displacements of rockmass on slopes earthquake[J].JournalofEngineeringGeology, 1997, 5(2): 131-136.

[8] 中国建筑科学研究院.建筑抗震设计规范:GB 50011—2010[S].北京:中国建筑工业出版社,2010. China Academy of Building Research.CodeforSeismicDesignofBuilding: GB 50011—2010[S]. Beijing: China Building Industry Press, 2010.

[9] 江苏南京地质工程勘察院.南京牛首山文化旅游区-核心区岩土工程详细勘察报告[R].南京:江苏南京地质工程勘察院,2013. Jiangsu Nanjing Institute of Geo-Engineering Investigation.DetailedGeotechnicalEngineeringSurveyReportfortheCoreAreaofNanjingNiushouCulturalTourismZone[R]. Nanjing: Jiangsu Nanjing Institute of Geo-Engineering Investigation, 2013.

[10] 重庆市设计院,中国建筑技术集团有限公司.建筑边坡工程技术规范:GB 50330—2013[S].北京:中国建筑工业出版社,2014. Chongqing Architectural Design Institute, China Building Technique Group Co.,Ltd..TechnicalCodeforBuildingSlopeEngineering: GB 50330—2013[S].Beijing: China Building Industry Press, 2014.

[11] 江苏省地震工程研究院.南京江宁牛首山胜境一期工程场地地震安全性评价报告[R].南京:江苏省地震工程研究院,2012. Jiangsu Earthquake Engineering Institute.SiteSeismicSafetyEvaluationReportforOnePhaseofNanjingJiangningNiushouMountainScenicSpotProject[R].Nanjing: Jiangsu Earthquake Engineering Institute, 2012.

Three-dimensional Static and Seismic Finite Element Simulation and Analysis of Foding Palace Mine Pit

ZHANG Julian, LIANG Zhirong, LI Wei

(Shanghai Shenyuan Geotechnical Engineering Co.,Ltd., Shanghai 200040, P.R.China)

The first domestic existing pit-Foding palace mine pit slope under static state and seismic condition was simulated and analyzed by three dimensional finite element. Some results were obtained that after which cutting or reinforcement with anchor cable(rod), the maximum displacement all occurred on the slope top and the surface was filled with soil, residual soil and highly weathered rock; the maximum stress all occurred in the depth inside slope and the maximum shear stress occurred at interface between slope toe and ground. The acceleration response is in trend of expanding with increase of slope height and the surface soil acceleration response expanded significantly. Under limit state condition, the generalized plastic strain area concentrated on slope top and surface. Under static loading case, with anchor cable(rod) reinforcement, the maximum side slope displacement was reduced by 11.9%～28.4%, stress decrease by 11.9%～28.4%, 5.6%～67.1% respectively while slope safety factor increased by 34.5% compared with the condition without support, all of which showed strong strengthening effect of anchor cable(rod). Under seismic condition, compared with static condition, maximum displacement was increased by 3.03～19.73 times, main stress was increased by 6.5%～146.7%, maximum shear stress was increased by 168.4%～474.6%, safety factor was decreased by 0.6%～23.1%. Under major earthquake condition, compared with small earthquake condition, slope surface acceleration peak and magnification was increased by 15.04～27.00 folds and 1.44～1.93 folds respectively. The calculation results from applying earthquake on the mine pit from initial north-south and later east-west directions respectively show that the former generated greater seismic displacement and dynamic stress than the latter and the former safety factor was less than the latter indicating that the earthquake striking from pit longer axis direction is worse than that striking from shorter axis to slope stability.

tunnel engineering; existing pit; earthquake; displacement; stress; seismic acceleration response

10.3969/j.issn.1674-0696.2016.04.04

2015-10-12；

2015-10-26

上海市科技人才计划基金项目(14XD1420400)

张菊连(1984—)，女，江西丰城人，高级工程师，博士，主要从事基坑工程及边坡治理方面的研究。E-mail：327804253@qq.com。

TU457；U456.2

A

1674-0696(2016)04-013-07