FLAC3D中锚杆剪切断裂的数值实现与应用

张志强，刘 银，赵梓彤，李化云，2

(1.西南交通大学 交通隧道工程教育部重点实验室, 四川 成都 610031；2.西华大学 土木建筑与环境学院，四川 成都 610039)

作为新奥法隧道支护体系中唯一能深入围岩内部的支护结构，锚杆在地下工程围岩稳定性控制方面有着广泛的应用[1-4]。对处于大埋深、高地应力等复杂地质条件下的地下工程(交通隧道、煤矿巷道、水工隧道等)而言，当围岩中存在软弱结构面时，结构面两侧岩体的相互错动会使穿过结构面的锚杆承受较大的剪切作用，此时，锚杆容易产生剪切断裂而失去对围岩形变的控制能力，甚至导致整个支护体系失效，引发围岩大变形灾害[5-7]。

国内外学者针对锚固体中锚杆的抗剪机制以及地下工程中锚杆的剪断机理进行了许多研究。张伟等[8]基于锚杆的切向和轴向形变特征，建立了节理岩体锚固体中锚杆的剪切变形分析模型；Li等[9]通过室内试验，对比分析了玻璃纤维锚杆与螺纹钢锚杆锚固体在剪切荷载作用下的力学行为；陈文强等[10]建立了加锚结构面抗剪强度计算模型，对剪切过程中锚杆的轴向和横向作用进行了分析；Forbes等[11]通过新的光纤测试技术，对锚固体剪切过程中锚杆的力学行为进行了分析；肖同强等[6]通过建立肩角锚杆力学模型，揭示了深部构造应力区煤巷肩角锚杆剪断机制；缪协兴等[12]分析了综放沿空巷道顶部全长黏结式锚杆承受剪切变形的机理；刘金海等[13]研究发现，深埋巷道肩角、肩窝处锚杆最容易在剪应力作用下产生断裂。同时，数值计算也常用于对锚杆的支护效果进行分析。但目前的数值计算中[14-16]，均未考虑锚杆在剪切荷载作用下的断裂效应。因此，导致计算结果与实际支护效果不符，夸大锚杆的作用。

FLAC3D软件在岩土工程大变形计算方面有着广泛的应用。根据特定研究目的，研究者采用FISH语言对软件自带的结构单元进行了针对性的二次开发。如李为腾等[17]将拉断准则引入FLAC3D中的CABLE单元，实现了锚杆的局部拉伸断裂；刘少虹[18]提出了爆破震动下锚杆单元破坏损伤判据，建立了基于FLAC3D平台的锚杆轴力动力损失修正模型；言志信[19]等将界面破坏时砂浆的应力软化效应引入锚杆结构单元，分析了地震作用下岩质边坡锚固界面的剪切失效过程。但这些二次开发模型均未考虑锚杆的剪断效应。

因此，本文基于FLAC3D软件平台，对PILE结构单元进行修正，使其能准确地反映锚固体中锚杆剪切形变的全过程；然后将修正PILE单元用于地下工程围岩稳定性分析中，对锚杆的剪断与拉断失效机制进行研究。

1 FLAC3D中PILE单元力学特性

1.1 PILE单元基本特性

FLAC3D中，单个PILE单元由2个节点以及节点间的结构单元构成，结构单元具有相同的截面形状及物理力学特性。当PILE单元被锚固剂包裹时，其与锚固剂之间通过法向与切向弹簧—滑块体系来实现形变的耦合与内力的传递，如图1所示。

图1 锚杆单元结构组成

1.2 PILE单元剪力学特性

剪切荷载作用下，PILE单元承受的剪力随着剪切位移的增加而线性增长。剪力与剪切应变之间的关系为

Fs=γGA

(1)

式中：Fs为锚杆杆体所承受的剪力；γ为剪切应变；G为杆体的剪切刚度；A为杆体的截面积。

1.3 PILE单元拉伸力学特性

PILE单元中，引入锚杆的拉断准则，即

εt≥εtmax

(2)

式中：εt为杆体的应变；εtmax为杆体的极限应变，其取值与材料的形变能力有关。

杆体的轴向本构关系如图2所示，其拉伸表达式为

(3)

式中：Ft为杆体的拉力；Et为杆体的弹性模量；Ftmax为杆体的极限抗拉承载力。

可见，锚杆杆体存在以下3种受拉力学状态。

状态1：杆体所受到的外力小于其极限抗拉能力，且杆体的应变小于其极限应变。此时杆体处于弹性阶段，杆体轴力随着形变值的增加而线性增长。

状态2：杆体承受的外力达到其极限承载能力，而杆体的应变小于其极限应变。此时杆体轴力保持不变，应变继续增加。

状态3：杆体的应变达到极限应变。此时杆体产生断裂，失去承载能力。

通过上述对PILE单元的分析可以发现，PILE单元可以较好地模拟锚杆的拉伸力学效应，但不能真实地反映锚杆的剪力学效应。

图2 拉伸力学特性

2 PILE单元剪断力学模型的实现

2.1 锚杆剪力学模型

Li等[7]通过对锚杆进行剪切试验，发现锚杆的剪力与剪切位移之间存在着图3(a)所示的关系。可以看出，在杆体剪力达到极限承载能力之前，剪力随着剪切位移的增加而近似线性增大(OA段)；当杆体剪力达到剪切极限力时，杆体的剪力会突降至其最大残余承载能力(AB段)；在BC段，随着剪切位移的增加，杆体承担的剪力出现波动，但整体上呈现出减小的趋势；当杆体的剪切位移超过其极限剪切位移时，锚杆产生断裂，失去承载能力。

根据文献[6]的试验结果，可得到简化的锚杆剪力学模型如图3(b)所示，其表达式为

(4)

式中：Fs,max为杆体所能承担的极限剪切荷载；γmax为杆体极限剪切应变；Fs1和Fs2为峰后特征点B和C的剪力；γs1为剪力峰值点对应的剪切应变。

图3 锚杆剪力学模型

2.2 PILE单元剪断力学模型的数值实现

为了使PILE单元能够真实地反映锚杆的剪力学效应，采用FISH语言，以式(4)为基础对PILE单元的剪力学行为进行修正，修正流程如图4所示。具体步骤如下。

(1)采用主程序计算得到每一个杆体单元的剪力与剪切应变。

(2)判断杆体单元受到的剪力是否超过其极限承载能力，若超过，则根据杆体单元的剪切应变通过插值计算得到此时杆体所能承担的最大剪力；同时，判断杆体的剪切应变是否超过其极限剪切应变，若超过，则将杆体的剪切极限力设为0；锚杆单元的轴力按式(3)进行计算。

(3)判断计算是否收敛，若收敛，则计算结束；若不收敛，则重复步骤(1)和步骤(2)。

图4 PILE单元剪断开发流程

3 锚杆剪断力学模型的验证

以文献[7]的锚杆剪切试验结果为基础，对本文锚杆剪断力学模型的有效性进行验证。建立的数值模型如图5所示，混凝土试样采用实体单元进行模拟，尺寸为0.4 m×0.4 m×0.4 m。混凝土单元的下端固定，上端为自由面，其余侧面施加相应方向的法向约束。锚杆直径为28 mm，所能承受的最大剪切荷载Fsmax=279 kN，峰后剪力特征值Fs1=150 kN，Fs2=80 kN，极限抗拉承载力Ftmax=558 kN，极限剪切应变γmax=83 mm。锚杆长度L=0.55 m，划分成15个PILE杆体单元，杆体的CID编号为1～15。杆体与锚固剂参数见表1和表2。

表1 锚杆杆体力学参数

表2 锚固剂力学参数

图5 剪切数值模型

计算过程中，在A点施加沿x方向的恒定速度，大小为1×10-6m·计算步-1，对编号为11的PILE单元，记录其不同剪切位移的剪力计算值，结果如图6所示；同时，将试验值也绘于图6中。

图6 剪力与剪切位移曲线室内试验与数值计算结果对比

从图6可以看出：修正模型可以较好地反映出锚杆剪切过程中的弹性阶段、塑性变形阶段以及破坏阶段；弹性阶段内，当杆体单元产生同样的剪切位移时，数值计算得到的剪力小于室内试验结果，且两者之间的差值随着剪切位移的增加而减小；塑性变形阶段，室内试验的波动幅度较数值计算大，剪力整体上均呈现出减小的特征；破坏阶段，锚杆的剪力变为0。可以看出，虽然本文的修正模型认为锚杆的残余抗剪强度与剪切位移之间为线性关系(图3(b)中BC段)，但由于锚杆在塑性剪切变形阶段内形变值会持续波动，导致数值计算结果得到的剪力在该阶段内也持续波动，但整体上体现出线性减小的趋势。

原始PILE单元、修正PILE单元最终形变与室内试验实物的结果对比如图7所示。可以看出：原始PILE单元不能反映锚杆在剪力作用下的非连续大变形特征，而修正PILE单元可以较好地反映这一形变特征；修正PILE单元得到的计算结果与室内试验结果吻合较好，可以较好地反映锚杆剪切形变全过程的力学特征。

图7 最终形变室内试验与数值计算结果对比

4 锚杆剪断力学模型的应用

单弱面层状岩体是隧道工程建设中常见的一类岩体，采用本文建立的修正PILE单元，分析该隧道中锚杆剪断及拉断效应对围岩稳定性的影响。

4.1 数值模型

建立的数值模型如图8所示。模型左右宽度均为40 m(约为4倍开挖跨度)，仰拱以下高度为40m，模型厚度为2 m。前、后、左、右边界以及上边界施加相应方向的法向应力边界，下边界约束竖向位移，地应力取值为σx=σy=σz=10 MPa。层理弱面与水平方向的夹角θ=0°，45°，90°，层理面的强度与刚度取值为围岩基质体的1/2。围岩基质体参数见表3。

图8 层状岩体隧道计算模型

表3 复合体及围岩基质体的参数

喷射混凝土采用实体单元模拟，钢拱架的作用按其弹性模量折算给喷射混凝土[20]，即

(5)

式中：E为喷射混凝土—钢拱架复合体弹性模量；E0为喷射混凝土弹性模量；Ag为钢拱架截面积；Eg为钢材弹性模量；Sc为喷射混凝土截面积。

等效后的复合体力学参数见表3-复合体初期支护的厚度为20 cm。锚杆长度为3.7 m，其布置方式如图8(b)所示，锚杆的环向间距为1 m，纵向间距为1.2 m。为了后续描述方便，对锚杆从左下角沿逆时针方向依次进行编号，编号为C1—C24。锚杆为全长黏结式锚杆，杆体和锚固剂的力学参数分别见表1和表2。采用PILE 结构单元模拟锚杆，为了提高模拟精度，每根锚杆划分为37 段，每段长0.1 m，杆体所能承担的极限剪切荷载Fsmax=140 kN，抗拉载荷为255 kN，拉伸延伸率为16%。

4.2 计算结果分析

4.2.1 锚杆断裂力学特性

以θ=45°为例，通过对锚杆剪断单元(锚杆编号为C3，CID编号为77)与锚杆拉断单元(锚杆编号为C2，CID编号为43)的拉力与剪力随计算步的变化情况进行记录，分析锚杆的断裂力学特性，结果如图9所示，可以看出，剪断锚杆单元与拉断锚杆单元的力学行为有着较大的差异。对于拉伸断裂单元，其计算过程中的拉力最大值为其抗拉能力极限值；对于剪切断裂单元，计算过程中的拉力最大值小于抗拉极限值。

1)剪断锚杆单元

单元的剪力变化过程(图9(a))为：初始变形阶段，剪力随着计算步的增加而近似线性增加；当计算至1 120步时，剪力为140 kN，达到剪切破坏临界值；此后，单元在下一计算步内断裂，剪力骤降至残余抗剪能力最大值(60 kN)；当计算至1 310步时，剪力降至40 kN，此时单元的应变达到其剪切极限应变；随着计算步的继续增加，锚杆出现剪切断裂，剪力从40 kN降至0 kN，之后随着计算步的增加保持不变。

单元的拉力变化过程(图9(b))为：初始形变阶段(计算步小于400步)，拉力近似线性增加；当计算步大于400步后，拉力出现剧烈的波动上升；当锚杆剪切断裂时，单元的拉力从190 kN降至0 kN，之后随着计算步的增加保持不变。

可见，当锚杆单元剪断时，单元的拉力并未达到抗拉能力极限值(225 kN)。因此，在数值计算中，对于受剪切荷载作用较大的锚杆而言，若不考虑锚杆剪断效应，则会同时夸大锚杆的抗拉与抗剪支护能力。

2)锚杆拉断单元

单元的拉力变化过程(图9(c))为：拉力先近似线性增加(计算步小于400步)，后波动增加(400步到890步)；当计算至890步时，单元所受到的拉力达到了抗拉能力最大值，即255 kN，此时，由于单元的拉伸应变未达到其极限拉伸应变，因此随着计算步的进一步增加，单元的拉力保持在255 kN 不变；当计算至1 390步时，单元的应变达到极限应变，此时锚杆拉伸断裂，拉力降至0 kN，后随着计算步的增加拉力保持不变。

单元的剪力变化过程(图9(d))为：在锚杆拉伸断裂之前，单元所受到的剪力均较小，最大值仅为4.8 kN，且随着计算步的增加波动变化。当锚杆拉伸断裂时，单元的剪力也降为0 kN，后随着计算步的增加剪力保持不变。可见，对于以抗拉作用为主的锚杆单元而言，当剪切荷载小于极限抗剪能力时，不考虑锚杆的剪断效应不会对锚杆的支护效果产生影响。

3)内力沿杆体的分布规律

剪切断裂锚杆(C3)与拉伸断裂锚杆(C2)各单元内力沿杆体的分布规律如图10所示。从图中可以看出：当锚杆发生破断后，锚杆的剪力沿长度方向基本全变为0 kN，说明此时锚杆在全长范围内失去了抗剪切的能力；而锚杆的拉力沿长度方向的分布规律为中间大两端小，说明此时锚杆还能对围岩的局部拉伸变形起到一定的抑制作用。

图9 锚杆单元剪断和拉断力学特性

图10 锚杆内力分布

4)锚杆支护体系断裂过程

计算至平衡状态时，所有产生断裂的锚杆其断裂位置处单元的拉力随时间变化曲线如图11所示。从图可以看出：锚杆的断裂顺序与断裂模式为C3→C22→C2→C4→C20→C5→C21→C23→C9→C16→C18→C8(粗体为拉伸断裂，下同)，即首先是左右边墙处的锚杆被剪断，之后剪断锚杆相邻位置处的锚杆被拉断，然后左右拱腰处锚杆被拉断，最后左右拱肩处锚杆被剪断。可见，该工况下，边墙和拱腰附近锚杆容易产生拉伸或剪切破坏，拱肩附近锚杆容易产生剪切破坏。

4.2.2 不同倾角条件下锚杆的断裂过程

θ=0°，90°时，锚杆的断裂过程如图12所示。从图12(a)可以看出：当θ=0°时，围岩未产生大的相对剪切滑动，锚杆只产生了拉伸断裂，断裂顺序为C3(C21)→C5(C18)→C8(C16)→C9(C15，C10，C14)，即边墙附近锚杆首先发生断裂，然后，锚杆的断裂部位集中分布在拱腰与拱顶之间。当θ=90°时，锚杆的断裂顺序为C10→C15→C11→C9→C16→C19(C18)→C6(C7)→C14，即拱肩与拱顶之间的锚杆首先产生剪切断裂，然后拱腰与拱肩之间的锚杆产生拉伸断裂。

4.2.3 与原始PILE单元结果的比较

1)锚杆受力

以θ=45°为例，采用原始PILE单元的锚杆受力如图13所示。由图13可知：由于未考虑锚杆剪断，C3锚杆编号为77的杆体单元所承受的剪力随着计算步的增加持续增加，承受的拉力先持续增加至抗拉能力极限值，之后保持拉力极限值不变；从锚杆断裂过程来看，锚杆的拉伸断裂顺序为C2→C4→C6(C22,C21)→C20→C17，集中在边墙与拱腰附近。

图11 锚杆破断过程

图12 不同倾角下锚杆破断过程

图13 原始模型锚杆受力特性(θ=45°)

将图13与采用改进PILE单元的图9和图10对比可知：锚杆断裂数量少于改进PILE单元计算结果，且断裂开始的时间要滞后于改进PILE单元的计算结果。

2)围岩特征位移

各种工况下围岩的特征位移见表4，可以看出，对于不同层理面倾角，修正PILE单元与原始PILE单元计算结果之间的差异具有不同的规律。当θ=0°时，由于锚杆未产生剪切断裂，修正模型与原始模型的计算结果相等。当θ=45°时，修正模型的拱顶沉降与洞周收敛比原始模型分别小17和44 mm，洞周收敛之间的差值更大；这是由于锚杆的剪切断裂主要发生在左右边墙附近，锚杆剪断会更多地减弱支护系统对围岩水平形变的约束能力。当θ=90°，修正模型的拱顶沉降与洞周收敛比原始模型分别小27和17 mm，拱顶沉降之间的差值更大；这是由于锚杆的剪切断裂主要发生在拱肩与拱顶附近，锚杆剪断会更多地减弱支护系统对围岩拱顶竖向形变的约束能力。

表4 围岩特征位移值

5 结 论

本文基于FLAC3D平台，将锚杆剪切断裂模型引入软件自带的PILE结构单元，对PILE单元的剪力学行为进行修正。采用修正的PILE单元对单弱面层状岩体中锚杆的剪断与拉断失效机制进行了研究，得到以下结论。

(1)提出了PILE单元在剪切荷载作用下的全过程修正力学模型。该模型增加了锚杆剪切变形的塑性阶段与破裂阶段，即当剪力达到抗剪能力极限值时，锚杆进入塑性残余变形阶段；当剪切应变达到极限剪切应变时，锚杆剪切破断，杆体的轴力和剪力突降至零并保持恒定。

(2)锚杆剪断力学模型的验证结果表明，修正PILE单元得到的计算结果与室内试验结果吻合较好，可以较好地反映锚杆剪切形变全过程的力学特征以及剪切荷载作用处锚杆的非连续大变形。

(3)单弱面层状岩体隧道锚杆支护应用的对比分析结果表明，对于受剪切荷载作用较大的锚杆而言，不考虑锚杆的剪断效应会同时夸大锚杆的抗拉与抗剪支护能力。与不考虑锚杆剪断相比，锚杆剪断后，支护系统对围岩形变约束能力会有所减弱。但由于锚杆剪断位置的差异，不同层理面倾角条件下，支护系统约束能力减弱的程度与层理面的倾角有关。