落石冲击荷载作用下的桩板拦石墙结构动力响应*

0 引 言

崩塌落石是高山峡谷带常发的斜坡地质灾害之一，落石突发、高能的动力特性对人民生命财产、房建、交通及环境带来巨大威胁。中国、Nepal、Italy、Canada、Japan、Peruvian Andes(Jr，1998；Yamagishi，2000；张路青等，2005；Paronuzzi，2009；Procter et al.，2015；韩振华，2018；王学良，2018；祝介旺，2018)等地均有大量落石灾害实例报道。鉴于落石的危害性，各国学者开展了广泛的防治研究。Ritchie(1963)是第一个系统性研究落石问题的学者。Peila et al.(2007)开展了大比例加筋土拦石墙结构模型的抗冲击性能现场试验研究。Ronco et al.(2009)通过ABAQUS对Peila et al.(2007)的试验进行数值模型验证。Mougin et al.(2005)采用1/3缩尺模型试验，研究钢筋混凝土板在冲击荷载下的动力响应。王东坡等(2016)以Olsson冲击动力学控制方程为基础，提出了落石冲击钢筋混凝土板的计算方法并与有限元进行了对比分析。王爽等(2016)采用LS-DYNA动力有限元进行了棚洞抗冲击行为研究。Zineddin(2008)基于落锤冲击加载试验研究了不同配筋率钢筋混凝土板冲击作用下的动力响应。

图1 孔玉镇高位崩塌桩板拦石墙
Fig.1 Pile-plate rock retaining wall in Kongyu Town

目前国内外关于崩塌落石防治措施主要分为主动和被动防护。除了针对原位危岩能够实施主动加固外，受危岩崩塌区地形地质条件影响，至少有70%的危岩崩塌采用被动拦挡方案。大量的被动网、圬工拦石墙被采用，但是被动网质量控制目前尚无国家统一标准，实施效果一般或较差，而圬工拦石墙则由于自身基础宽度、落石槽侧的缓冲层厚度，在实施过程中往往占地过大。此外拦挡部位地形坡度较陡(25°～35°)时，墙后落石空间严重受限。针对这种情况，在2008年 “5·12” 汶川震区和2013年“4·20”芦山震区的高位危岩崩塌落石，一些学者提出了桩板拦石墙结构的被动拦挡方案，即在滑坡抗滑桩板结构基础上增加缓冲层而成的复式结构，具有占地面积小、有效拦截能力高、易于施工等特点，对于防护高冲击能量的落石存在无可比拟的优势，其中缓冲层土体可以有效减小结构体积，增大抗冲击能力，在实际应用中取得了良好效果(图1)。但这种桩板拦石墙结构中的桩间距、桩间板厚度、缓冲层材料及厚度等设计参数取值均是基于经验，缺乏令人信服的结构受力依据，而目前大多数研究主要集中于不同结构类型的棚洞结构(张群利等，2015)。作为一种相对较新的轻型防治结构，桩板拦石墙在高位崩塌落石冲击下的结构受力研究比较匮乏。为此本文依托四川省康定县孔玉乡寸达河坝后山危岩带落石防治工程实例，釆用有限元与无限元耦合的办法，开展桩板拦石墙动荷载作用下的动力响应数值模型研究，对揭示其力学特性的动态演变机理及设计具有现实意义。

1 落石冲击的桩板拦石墙动力响应数值模拟

1.1 工程概况

研究区位于康定县孔玉乡寸达河坝后山，受强降雨及2013年4.20芦山强烈地震影响，崩塌区危岩体变形加剧。2005年7月、2012年8月和2014年6～8月，发生多次落石灾害，共造成3人死亡，3人受伤及12辆车辆被毁，严重威胁坡下河坝村民及S211省道车辆和行人的生命财产安全。现场调查危岩带分布及地形地貌特征表明，边坡坡度介于35°～55°，无明显的缓坡平台，故采用桩板拦石墙进行治理(图2)。拦石结构设于居民点及S211省道以内靠山侧，主要包括混凝土桩-板及土工格栅组成的缓冲层。桩长20m，地面以上部分为10.6m，截面尺寸为1.5m×1.2m。缓冲层顶部宽度0.8m，底部2.7m，每0.6m铺设1层土工格栅。抗滑桩内部配有d=25mm纵向钢筋及d=16mm@200mm箍筋。桩间板厚0.4m，宽5.0m，板内以正交垂直方式上下铺设两层d=16mm@250mm钢筋网(图3、图4、图5)。

图2 孔玉镇典型危岩与保护对象剖面
Fig.2 Sectional view of typical perilous rock and protected object in Kongyu Town

图3 桩板拦石墙模型图
Fig.3 FEM model of structure

图4 结构侧立面图
Fig.4 Structure side elevation

图5 混凝土板及配筋
Fig.5 Reinforcement of slab

1.2 数值计算模型

利用ABAQUS/Explicit动力有限元对桩板拦石墙结构进行数值计算，落石简化为半径等于0.8m的刚性球体，质量为5361.7ikg，冲击作用点距墙顶部3.54m(图4)。结构形式与尺寸(图3)。由于结构延伸较大，累计总长255m，属于大纵横比结构，如果将整体作为计算模型，计算规模将以级数倍增加，因此本研究取其中21.2m作为建模长度。但是选取有限长度会使两端产生人工截断边界，进而导致边界效应。为了消除人工边界处引起波的反射、振动带来的误差，采用了有限元与无限元耦合的方式进行计算。“无限元”边界具有物理尺寸，但几何上趋于无穷远(图6)。通过相应的映射函数，使得无限元与有限元结合且在公共边界上的位移协调得到保证，而在远场位移趋于0，从而合理反映真实的边界条件(戚玉亮等，2014)。边界处使用无限元可以在不增加计算规模的前提下，通过吸收截断边界处应力波，实现波的逸散，从而增加计算精度，减小计算量。模型中，钢筋及土工格栅与基体的相互作用为“embedded”。考虑到混凝土桩与板在变形过程中始终贴合，因此采用“Tie”模拟桩体与板的连接作用。落石与缓冲层的接触方式采用有限滑移。

图6 三维单向映射无限元
Fig.6 Three-dimensional mapping infinite elementa.C3D8R单元；b.CIN3D8无限元

1.3 材料参数及本构模型

本文所涉及的材料共4种，其中混凝土材料在冲击作用下表现为微裂纹累计损伤、扩展及最终断裂动态破坏过程。为了表征这种弱化效应，用损伤因子D定义其损伤程度。当微裂纹积累达到破坏阀值时，混凝土发生脆性或粉碎性破坏。其损伤因子定义为：

(1)

式中，W0为无损混凝土的应变能密度；Wε为损伤混凝土的应变能密度。对于无损混凝土，Wε=W0，则D=0；对于损伤混凝土材料，0<Wε<W0，于是D≠0；在极限状态，W0≫Wε，则D→1。

本文采用混凝土CDP损伤模型，可以较好地模拟混凝土材料在动态加载时的主要性能(李敏等，2011)。根据王中强等(2004)的研究，通过混凝土本构方程，并引入损伤因子变量，建立损伤演变方程，进而可得到拉压损伤-应变关系曲线(图7，图8)。

按《混凝土结构设计规范》，拉压本构方程分别为：

(2)

式中，αa、αd为单轴受压应力-应变曲线上升段与下降段参数，按规范选取；fc为混凝土单轴抗压强度；εc为与fc对应的混凝土峰值压应变。

图7 混凝土受压损伤因子
Fig.7 Concrete compression damage factor

图8 混凝土受拉损伤因子
Fig.8 Concrete tensile damage factor

受拉本构方程为：

(3)

式中，αt为单轴受拉应力-应变曲线的上升段与下降段参数值，按规范选取；ft为混凝土单轴抗拉强度；εt为与ft对应的混凝土峰值拉应变。

桩板拦石墙使用的混凝土桩、板纵向钢筋与箍筋为率相关材料，为了考虑其动力增强效应，使用三折线模型，根据林峰等(2008)的研究，在不同的应变率下，钢筋的杨氏模量基本一致，钢筋的屈服强度和极限强抗拉变化规律为：

(4)

式中，为钢筋应变率；分别为动、静态屈服强度；fu，d、 fu，st分别为动、静态极限强度；D1、 D2为线性回归系数。一般落石冲击速度小于25m·s-1，因此本文对于钢筋HRB400选用 时的参数，钢筋HRB335选用 时的参数，具体见表1。

表1 钢筋材料参数(林峰等，2008)
Table 1 Parameters of reinforcement material(Lin et al.，2008)

缓冲层材料一般为现场开挖土体装袋砌筑而成，实际应用中有些采用土工格栅作为缓冲层加筋材料，但目前关于这类复合材料的冲击特性研究较少，根据Dong et al.(2011)的研究，格栅应力-应变满足理想弹塑性，相关参数(表3)。缓冲层材料使用Drucker-Prager较合理，混凝土材料参数(表2)。

表2 混凝土及缓冲层土体材料参数(周珉等，2017)
Table 2 Parameters of concrete and cushion soil material(Zhou et al.，2017)

图9 典型土工格栅本构(Dong，2011)
Fig.9 Typical geogrid constitutive(Dong，2011)

表3 土工格栅材料参数
Table 3 Parameters of geogrid material

2 计算结果与分析

2.1 系统能量与冲击力分析

Muraishi et al.(2005)研究表明，冲击能量小于100ikJ的落石事件所占比例为68%，冲击能量小于1000ikJ的落石所占比例达到了90%。Chau et al.(2002)的研究表明，落石冲击过程中旋转动能仅为平动能量的10%左右，因此本文分析中忽略角速度，冲击能量均是由落石动能转化。选取滚石冲击速度V分别为10m·s-1、 15m·s-1、 20m·s-1、 25m·s-1，对应冲击平动动能分别为268ikJ、 603ikJ、1072ikJ、1675ikJ。冲击角度θ分别取30°、45°、60°、75°、90°，通过计算每种工况下动力响应，进行横纵向对比分析。通过分析能量间的相互转化关系，量化伪应变能所占系统总能量百分比是评价计算结果合理性的主要方法，一般而言，伪应变能占总能量的10%以下，计算结果可靠。以冲击速度15m·s-1，冲击角度30°工况为例(图10)，系统总能量为563kJ，冲击过程中最大伪应变能8.506kJ，仅占总能量的1.45%，其他各种工况下，伪应变能所占比例均在5%以下。此外，系统总能量在整个分析过程中保持不变，系统内能与弹性应变能先增大后保持稳定，动能由最大逐渐减小，在回弹阶段有所增大，碰撞结束后趋于稳定，总体上系统各部分能量的变化趋势均符合相关理论，本文计算是合理的。

图10 系统能量历程曲线
Fig.10 Duration curves of system energy

图11表明，冲击总时长介于0.08～0.11s，在相同冲击角度下，接触时间随速度的增大而有所增加，冲击力随冲击速度的增加而增大，在0.04s冲击力达到最大，在不同的冲击速度下，峰值冲击力分别为1.2 MN、1.7 MN、2.1 MN、2.5 MN。冲击力时程曲线均为单峰值曲线，整体波动较小。曲线趋势表明，本文考虑的“无限元”边界，可以有效吸收截断边界处应力波，避免应力波反射造成的多峰或数据波动剧烈所造成的误差。

图11 不同冲击速度下冲击力时程曲线
Fig.11 Duration curves of impact force at different speeds

为了研究落石与冲击角度的关系，计算冲击速度为15m·s-1，冲击角度分别为30°、45°、60°、75°、90°情况下冲击力的变化规律(图12)，法向冲击力随冲击角度的增加而增大，而切向冲击力随冲击角度的增加而减小，在90°的时候，切向冲击力趋于0。各个冲击角度下最大法向冲击力分别为1.16MN、1.70MN、2.14MN、2.41MN、2.51MN，随冲击角度增加，法向冲击力的增幅有所减小。峰值切向冲击力在各个角度下分别为0.40MN、0.46MN、0.44MN、0.28MN、0.05MN，相同冲击速度下，冲击角度为45°时，切向冲击力最大。根据力的合成原理，得到每个冲击角度下合冲击力的最大值，分别为1.21MN、1.73MN、2.15MN、2.41MN、2.51MN。相对于法向冲击力而言，随着冲击角度的增加，切向冲击力占合冲击力的比例愈来愈小。此外，法向峰值冲击力在0.037s时刻出现，而切向冲击力在0.049s达到最大，法向冲击力与切向冲击力同时开始并结束，但非同时达到最大，峰值切向冲击力滞后于峰值法向冲击力，这与波动理论结果相似(Zhang et al.，2015)。

图12 不同冲击角度下冲击力时程曲线
Fig.12 Duration curves of impact force at different angles

图13得到了冲击角度与峰值冲击力之间的关系，各冲击速度对应的最大合冲击力随冲击角度的增加而增大，但75°后冲击力增加的趋势较平缓。

图13 冲击角度与峰值冲击力关系
Fig.13 Relationship between impact angle and impact force

2.2 冲击力计算模型结果比较

为了分析本文计算的合理性，选取Labiouse et al.(1966)、Kawahara et al.(2006)、Pichler et al.(2006)及Hertz(1997)弹性理论与本文计算结果进行对比分析(图14)。结果表明，冲击力与速度基本为线性相关，而 Hertz弹性理论计算结果整体最小。相同条件下Pichler模型冲击力随速度变化率最大，在冲击速度分别为10m·s-1、15m·s-1、20m·s-1、25m·s-1 的情况下，本文计算最大冲击力为1.9MN、2.5MN、3.1MN、3.7MN，Kawahara模型计算结果分别为1.3MN、2.0MN、2.9MN、3.7MN。整体上本文的计算结果与其他模型均处于相同数量级，当速度大于 15m·s-1 时，本文FEM解与Kawahara模型解趋于一致。

图14 冲击力与冲击速度关系
Fig.14 Relationship between impact speed and impact force

2.3 冲击位移分析

图15 不同冲击角度下冲击深度时程曲线
Fig.15 Duration curves of depth at different impact angles

图15得到了不同冲击角度下法向最大冲击深度时程曲线。冲击深度曲线整体分为上升段、下降段与平稳段。在弹塑性压入第1阶段，冲击速度减小、冲击力增大，进而冲击深度增大。当冲击速度减小到0时，此时冲击力与冲击深度均达到最大。在卸荷回弹阶段，土体发生弹性恢复，冲击深度有所减小，除一定的弹性恢复外其余均为永久塑性变形。最大冲击深度与最大冲击力均在0.04s左右所达到最大，两者的峰值时刻基本对应，随着冲击角度的增加，可恢复的弹性变形有小幅增大。

按经典弹塑性理论(Johnson，1992)可知，冲击过程可分为典型弹塑性压入与卸荷回弹两阶段，第1阶段又可分为弹性压入与塑性压入阶段。分析不同冲击速度下的加载-卸载曲线(图16)表明，压入阶段，冲击深度与最大冲击力呈正相关，冲击初始阶段，土体未发生屈服，加载路径表现为直线。此后塑性屈服区发生在土体的某一点，并随着冲击力的增大而扩展，当冲击速度减小至0时，发生卸荷回弹，回弹路径为直线且与加载路径不同，这表明卸荷段仍然满足弹性理论，卸荷结束除小部分回弹变形，其余均为永久塑性变形。在冲击速度为10m·s-1、 15m·s-1、 20m·s-1 的情况下，最大冲击深度分别为0.24m、0.33m、0.42m，最终塑性变形分别为0.19m、0.26m、0.37m，回弹变形量为0.03m、0.07m、0.05m，回弹变形占总塑性变形比例很小。

图16 不同冲击速度下加载-卸载曲线
Fig.16 Load-unload curves at different impact speeds

图17表明，冲击深度随角度的变化趋势较冲击角度和冲击力关系曲线相似，10m·s-1、15m·s-1、 20m·s-1、25m·s-1 均在冲击角度为90°时达到最大。

由图18可知，以冲击时间为中间变量，根据速度的分解法则表明，得到冲击过程中法向冲击深度随法向冲击速度减小，在法向速度为0时，动能全部转化为结构的应变能，此时达到最大冲击深度。当冲击角度为90°时，法向速度最大，冲击深度也最大。

图17 冲击角度与最大冲击深度关系
Fig.17 Relationship between impact angle and maximum depth

图18 最大冲击深度随速度变化曲线
Fig.18 Variation curves of maximum depth with speeds

图19 质心运动轨迹曲线
Fig.19 Motion tracks of centroid

观察质心运动轨迹曲线(图19)，随冲击角度的增大，法向位移增加，而切向位移减小，当冲击角为90°时，法向位移最大，切向位移基本为0。

2.4 结构耗能分析

图20 不同部件吸收能量曲线
Fig.20 Absorption energy curve of different components

本文所采用的缓冲层是由土工格栅与土体组成的复合材料，为了充分了解这种新型复合缓冲层的缓冲耗能效果，分析不同冲击角度下结构整体、土体材料及土工格栅在冲击过程中非弹性耗能情况得到图20曲线。结果表明，在20m·s-1 的高速冲击作用下，结构整体产生的非弹性耗能随冲击角度的增加而增大最终趋于稳定，在30°、45°、60°冲击角度下整体非弹性耗能曲线与土体曲线重合，表明结构所有的塑性应变能均被土体所吸收，格栅与其他部件基本没有吸收落石产生的能量，在75°与90°冲击角度下，结构整体非弹性耗能大幅增大，除土体与仅格栅吸收了部分塑性变形能外，剩余的部分被结构的其他部件吸收，产生如混凝土桩、板的损伤及钢筋的应变能等。相比于土体耗能作用，土工格栅整体产生的非弹性耗能较小，最大仅为5300J，仅占总非弹性耗能的6%左右，因此对于桩板拦石墙结构而言，缓冲层土体是最为有效的耗能部件。

2.5 混凝土损伤分析

图21 20m·s-1 冲击速度下板损伤情况
Fig.21 Plate damage at the impact speed of 20m·s-1

图22 20m·s-1 冲击速度下桩损伤情况
Fig.22 Pile damage at the impact speed of 20m·s-1

图23 不同冲击速度下板损伤情况
Fig.23 Plate damage at the impact speed of 20m·s-1

图24 不同冲击速度下桩损伤情况
Fig.24 Pile damage at different impact speeds

研究及实践表明，由于混凝土初始缺陷，冲击荷载作用下拦石墙中混凝土构件的响应主要表现为损伤与劣化，根据前面研究可知，桩板拦石墙在落石冲击作用下，其中75%以上的冲击能量均被缓冲层吸收，其余部分由格栅、钢筋、混凝土吸收。其中部分能量产生钢筋的塑性变形及混凝土的不可逆损伤。如果冲击速度过大或缓冲层后面的混凝土桩、板设计不合理，易造成结构的损伤甚至倒塌丧失承载力。如何定量描述混凝土结构在冲击作用下的劣化程度，合理定义有效的损伤准则，需要考虑：(1)钢筋混凝土结构的整体特性。(2)评价指标易获取且具有一般适用性。本文按照损伤因子的定义，将钢筋混凝土破坏程度分为4个等级：

D∈(0.0，0.25] 轻度损伤

D∈(0.25，0.5] 中度损伤

D∈(0.5，0.75] 重度损伤

D∈(0.75，1.0] 完全损伤

根据以上定义，分别计算了在不同工况下混凝土桩、板在冲击作用下产生的损伤体积，在冲击速度为10m·s-1 及15m·s-1 工况下，混凝土桩-板无论是受压损伤还是受拉损伤均比较小，缓冲层的存在能够吸收大部分能量，从而有效地保护桩板结构完整性(图23，图24)，同时验证了前文图20的结论。以冲击速度20m·s-1 为例，在30°入射角情况下，混凝土板开始发生轻度受拉损伤(图21)，而桩结构在45°时开始发生受拉损伤(图22)。桩和板的损伤程度及损伤百分比均随着冲击角度的增大而增加，总体上看，同一部件受拉损伤程度与损伤体积百分比均要远大于受压损伤，靠近碰撞中心部位损伤明显大于远离碰撞中心区域(桩、板相互位置关系见图2)，损伤由近及远扩展，在冲击角度为90°时，中间板产生完全受拉损伤的混凝土体积约67%，而处于完全受压损伤状态的混凝土仅有不到20%。在离碰撞中心较远的两边板，拉压损伤基本对称，且损伤程度远低于中间板。图22为混凝土桩的损伤情况，总体上桩受到的拉压损伤要轻于板，但在90°的冲击角度下，中间两根桩仍有56%的混凝土处于完全受拉损伤状态，有60%的混凝土处于重度受拉损伤状态。图23～图24给出了在90°冲击角度下不同冲击速度桩、板的损伤情况。在25m·s-1，90°冲击工况下中间板有75%混凝土处于完全受拉损伤状态，42%混凝土处于重度受压损伤状态。对于4根混凝土桩而言，中间两桩在20m·s-1 及25m·s-1 的冲击速度下损伤程度较严重，完全受拉损伤混凝土体积百分比最大为72%，同时50%的混凝土由于压应力而导致重度损伤。

3 结 论

结合经典弹塑性理论，采用ABAQUS动力有限元，在考虑钢筋与混凝土率相关性的前提下，模拟多种工况下滚石冲击桩板拦石墙的动力响应，得到以下结论：

(1)各种计算工况下，系统伪应变能所占比例均在5%以下。冲击过程中弹-塑性压入、弹性回跳的冲击过程与弹塑性理论完全吻合，以上均表明本文计算是合理的。

(2)碰撞过程中，法向冲击力与切向冲击力并非同时达到最大；法向冲击力随冲击角度的增加而增大，而切向冲击力随冲击角度的增加而减小；相对于法向冲击力，随冲击角度的增加，切向冲击力占合冲击力的比例愈来愈小。

(3)冲击力与速度基本线性相关，与几种经典的冲击力计算模型对比发现，相同情况下，本文计算结果与其他模型均处于相同数量级，当速度大于15m·s-1 时，本文FEM解与Kawahara模型解趋于一致。

(4)根据碰撞过程各构件产生的塑性应变能可知，缓冲层中土体部分是最主要的耗能部件。相对于土体缓冲层，格栅整体吸收的塑性应变能较小，最大仅占总塑性应变能的6%左右。相对永久塑性变形而言，弹性恢复的位移很小。

(5)根据混凝土损伤因子的定义，提出了一种量化结构损伤程度的方法，根据本文分析，在1072KJ～1675KJ冲击能量下，结构损伤较严重，存在完全丧失承载力的风险。考虑到天然条件下达到以上冲击能量为极端情况，在603KJ的中等冲击能量下结构整体损伤程度一般，本文研究的桩板拦石墙是可靠的。

后续将进一步开展模型试验研究，系统对冲击力、冲击位移等关键特征参量进行对比验证，以及对本文提出的混凝土损伤等级分类进行试验标定。