赵飞涛

（1.中国地质大学（武汉），湖北 武汉 430074；2.襄阳市建筑设计院，湖北 襄阳 441000）

0 引言

锚杆具有可施加预应力、施工简便等优点，在电力工程、水利、建筑、边坡、隧道等领域得到广泛的应用，持续受到工程人员和研究者关注［1-10］。随着时代的发展，锚杆材料从原始的钢筋、钢绞线逐渐引入复合纤维等新型材料，并发展出负泊松比锚杆、扩大头锚杆等多种类型［10-21］。

FLAC3D是常用的数值分析软件，基于有限差分法的开发，内置网格单元和beam、Cable、Pile 等结构单元，可分析各种小变形、大变形、几何非线性问题，广泛应用于锚杆承载变形特性的分析。根据目前的研究，FLAC3D中建立锚杆模型一般采用FLAC 中自带的Cable、Pile等结构单元，也有研究者采用实体网格配合接触面（Interface）单元进行模拟，另有研究者采用网格单元对钢筋横肋进行精细化模拟来研究钢筋与锚固体间的关系［5-15］。

Pile单元的轴向行为特性采用一维模型描述，粘结摩擦特征采用法向和剪切向耦合弹簧模拟，但一维线性模型很难对锚杆截面形态带来的影响进行充分考虑。为研究锚杆本身承载变形性能，构建一种网格单元与Pile 单元组合形成的锚杆，用Pile 单元模拟钢筋，用网格单元模拟锚固体，网格单元强度参数参考Pile单元耦合弹簧的强度参数取值，其抗拉试验模拟结果与单独采用Pile单元相比，最大抗拔力等相差较大。

1 一般粘结锚杆抗拉数值模拟

1.1 锚杆结构

锚杆一般由筋材和锚固体组成，一般分为自由段和锚固段，锚固段为注浆段与筋体紧密粘结段，锚固体一般由注浆体等凝固形成，筋材为钢绞线等柔性材料时也称为锚索。根据注浆体和筋体的关系又分为拉力型、压力型和两者的组合型。对于一般粘结拉力型锚杆，其主要结构如图1所示。

图1 一般粘结锚杆结构示意图Fig.1 Structural diagram of general bonded anchor rod

1.2 FLAC3D软件内置锚杆单元

FLAC3D软件模拟锚杆时，一般采用内置的Cable或Pile 结构单元。FLAC 软件手册将内置的Beam、Cable、Pile 等结构单元描述为由结构元素和节点组成的特定物理网格，与网格单元一样采用显式的拉格朗日算法求解，可与网格单元进行耦合。其中，锚杆、桩均由若干结构单元（SELs）和节点（node）组成，各个分段的质量都会集中在节点处。

锚杆结构单元每个节点有一个轴向的平动自由度，桩结构单元每个节点有6 个自由度，桩结构单元（pileSELs）除了提供梁（beamSELs）结构单元的结构性能外，桩与网格之间还存在法向（垂直于桩轴）和剪切（平行于桩轴）的摩擦作用［12］。

FLAC3D软件锚单元的轴向行为特性采用一维模型描述，粘结摩擦特征采用法向和剪切向耦合弹簧锚杆模拟（见图1），认为沿桩-网格界面的剪切力是界面粘结强度和界面（与压力相关）摩擦力的函数，并用式（1）对单位长度桩的最大剪切力进行了表述，可通过该公式对相关参数进行估算。另外，FLAC3D软件可通过自定义表格给出参数修正有效侧限应力，也可通过自定义cs_sctable 和cs_sftable 表来体现剪切向耦合-弹簧黏聚力和摩擦角的软化。

图2 FLAC3D中Pile结构单元的作用机理[4]Fig.2 Mechanism of pileSELs unit in FLAC3D[4]

Cable 结构单元可赋予的属性参数有12 个，每个桩结构单元包含6 个材料参数、4 个截面几何特性参数、10 个耦合特性参数。当使用锚杆逻辑时，还需要考虑另外的7 个参数，在模拟锚杆单元时一般需要输入的参数见表1。

表1 桩结构单元模拟锚杆时常用参数表Table 1 Common parameters of pileSELs unit in simulating anchor rod

FLAC3D软件手册指出，桩土相互作用的具体参数通常可以通过标准现场或实验室测试获取，其相互作用的剪切力能从岩土的参数中估算，与锚杆的截面积、弹性模量和屈服点等参数相比，获得有关水泥浆锚固体的参数显然更为困难。通常情况下，水泥浆锚固体剪切刚度Kbond都可通过拉拔试验测得，假设水泥浆环孔材料为理想弹塑性体，根据St.John，Van Dillen1983推导的水泥浆锚固体与岩石界面的剪应力方程也可近似获得理论公式（1）。

剪切刚度、钢筋表面到钻孔表面的相对剪切变形和剪切力之间的关系，经适当调整见公式（2）。

式（2）中：Kbond为水泥浆锚固体剪切刚度（cs_sk）；G为水泥浆锚固体的剪切模量；cssfric为剪切耦合弹簧摩擦角（cs_sfric）；D为筋材直径；t为水泥浆锚固体的厚度。

γ为考虑了发生在寄宿域网格点和钻孔表面之间产生的相对剪切位移影响时取0.1，不考虑时取1。

式（3）中：τmax为最大剪切强度。

FLAC3D软件手册指出，式（3）可在锚固系统的破坏发生在水泥浆和岩石的交界面时使用，在筋材拔出锚固体破坏时t按0考虑，在采用该公式进行计算预估最大剪切强度，粘结抗剪强度假定为浆液单轴抗压强度的一半［4］。

从这个角度看，软件中的锚固剂（锚固体）本质上也可看做是锚杆的剪切破裂面（带），在土体中计算锚固体的相关参数时需要根据剪切破裂面（带）设置，而非利用实际的锚固体强度参数进行计算。

1.3 锚杆模型试验方案

为研究锚杆抗拉的性能，首先根据软件内置的Pile 单元进行抗拔试验，相关参数根据软件手册要求设定，同时构建一种网格单元与Pile 结构单元组合形成的锚杆，用Pile结构单元模拟钢筋，用网格单元模拟锚固体，网格单元参数参考Pile单元锚固剂（剪切耦合弹簧）参数。

模型一为内置Pile 结构单元模型，将岩土体设置为弹塑性模型，锚固体设置为莫尔库伦模型，锚杆采用内置的Pile 结构单元模拟，模型岩土体的总尺寸为4 m×4 m×9.1 m，岩土体模型底部及周侧约束，模型顶部设置为不受约束的自由面，岩土体相关试验参数见表2；锚杆为一般粘结锚杆，布置在模型的中部，自由端长1.3 m，锚固体长度6.5 m，锚固体直径0.16 m，钢筋材直径1.524 cm，将锚固体周长设置为锚固体实际周长0.502 4 m，即软件默认设置方式，模拟所需的相关参数参考软件手册及案例设定，见表1、表2。

表2 岩土体物理力学参数表Table 2 Physical and mechanical parameters of rock-earth mass

模型二为网格单元与Pile结构单元组合形成的锚杆，用Pile 单元模拟钢筋，用网格单元模拟注浆体，网格采用莫尔库伦模型，将锚固体周长设置为0.068 m，即筋材横截面周长，其他设置与模型一保持一致，锚固体网格单元强度参数参考Pile 单元锚固剂（剪切耦合弹簧）强度参数。

首先，通过弹性模型求取初始应力场，清空速度位移场后模拟锚杆抗拔过程，在锚杆顶部节点设置一个Z（拉伸）方向的常速度边界（-1e-5m/step），并限制锚杆在X、Z（侧向）向的移动，通过锚杆顶部的拉伸，可监测Pile单元的轴力、剪应力变化，进而获取锚杆的最大抗拔力，计算中锚杆的顶部的位移等于计算开始后的步数（step）乘以常速度。

2 数值模拟结果分析

通过数值模拟试验可知（图3-图8），锚杆计算最大抗拔力均小于筋材的抗拉强度，锚杆均非拉断破坏。不同类型锚杆的计算抗拔力与时步关系曲线见图3，不同类型锚杆沿锚固段的轴力分布见图4，其最大抗拔力时的位移、轴力和剪应力云图见图5-图8。

图3 不同类型锚杆的计算抗拔力与时步关系曲线Fig.3 Relationship curve between calculated pullout force and time step of different types of anchor rods

图4 不同类型锚杆沿锚固段的轴力分布Fig.4 Axial force distribution of different-type anchor rods along the anchorage section

图5 Pile单元锚杆最大抗拔力时的位移云图Fig.5 Displacement nephogram of pileSELs unit anchor rod at maximum pullout force

图6 Pile单元和网格单元组合锚杆最大抗拔力时的位移云图Fig.6 Displacement nephogram of pileSELs unit and grid cell combined anchor rod at maximum pullout force

图7 Pile单元锚杆最大抗拔力时的剪应力云图Fig.7 Nephogram of shear stress at the maximum pullout force of pileSELs unit anchor rod

图8 Pile单元和网格单元组合锚杆最大抗拔力时的剪应力云图Fig.8 Nephogram of shear stress at the maximum pullout force of combined anchor rod of pileSELs unit and grid cell

Pile 单元锚杆在4 200 时步时开始拉伸，在5 960时步时获得的最大抗拔力是164 kN，锚杆顶部最大位移1.76 cm 时最大抗拔力是164 kN。网格单元与Pile单元组合形成的锚杆在4 200 时步时开始拉伸，在8 360时步时获得的最大抗拔力是416 kN，锚杆顶部最大位移4.16 cm时最大抗拔力是416 kN。

计算可知，组合结构锚杆与Pile单元锚杆，在采用相同参数计算时最大抗拔力相差比较大。

3 结论差异原因分析

根据FLAC 软件内置手册说明和相关研究可知，Cable、Pile单元的轴向行为特性采用一维模型描述，通过调整法向和剪切向耦合弹簧刚度、粘聚力、摩擦角来模拟锚杆结构单元和岩土体网格单元间的摩阻力。另外，FLAC 软件在考虑Cable、Pile 单元与网格单元的耦合时，假定锚杆的每个节点都位于某个FLAC 网格域内，该域被称为寄宿域，锚索节点到寄宿域各网格点的距离作为权重因子进行插值，权重因子的计算以满足弯矩平衡为原则，FLAC采用同样的插值方法将水泥浆与岩石交界面产生的力反馈到寄宿域网格点上［10］。

由此可见，锚杆与岩土体间都是通过Pile 单元节点（node）与网格单元（zone）进行相互作用的，节点上被赋予了桩的结构属性和耦合弹簧属性，所以锚杆单元节点与寄宿域网格点之间的相互联结方式、寄宿域网格节点特性等都会直接影响到计算结果。

4 结语

根据数值模拟结果，Pile 单元锚杆在锚杆顶部最大位移1.76 cm时得到最大抗拔力164 kN，网格单元与Pile 单元组合形成的锚杆在顶部最大位移4.16 cm 时最大抗拔力为416 kN，网格单元的强度参数按Pile 单元耦合弹簧的强度参数考虑。

FLAC 软件在考虑Pile 单元与岩土体网格单元的耦合时，假定Pile 单元的每个节点都位于某个网格域内，该域被称为寄宿域，Pile单元节点到寄宿域各网格点的距离作为权重因子进行插值。Pile单元节点上被赋予了桩的结构属性和耦合弹簧属性，所以Pile 单元节点与寄宿域网格点之间的相互联结方式、寄宿域网格节点特性等直接会影响到计算结果。

如果将Pile等结构单元按一维线性模型与岩土体进行耦合计算，锚杆结构物的尺寸效应将会被严重忽视，进而影响到计算结果，因此，在进行锚杆细部分析时需谨慎使用内置线性结构单元。