锚杆拉拔模拟中平衡状态对模拟结果的影响研究

杜学领

(贵州理工学院 矿业工程学院，贵州 贵阳 550003)

0 引 言

锚杆是支护工程中最常用的材料之一，其历史可追溯到1913年Frohlich和Klupfel申请的岩石锚杆专利。但是，锚杆的研究工作却起步较晚，一般认为，1949—1974年澳大利亚的Snowy Mountains Scheme项目执行期间开展的一系列研究，使人们对锚杆有了更加科学的认识，并由此确立了锚杆在工程中的永久支护地位[1]。在我国井工开采的煤矿中，目前也广泛采用锚杆作为巷道支护材料，并研发出注浆锚杆、可伸长锚杆、桁架锚杆等不同支护锚杆，极大地满足了不同地质条件对工程质量的需求[2]。学者们提出了悬吊、组合梁、组合拱、松动圈、中性点、围岩加固、突破点等锚杆支护理论[3]，并通过物理实验、数值模拟、工程实践等多种手段研究了锚杆及锚固体系的性能，其中拉拔试验是评价锚杆性能的重要手段。韦四江等[4]通过拉拔试验证明，拉拔速率有提高锚杆拉拔峰值载荷和残余强度的作用；胡江春等[5]研究认为，多项式时效方程在描述锚固段的蠕变过程方面要优于伯格斯模型；宋义敏等[6]通过试验发现，锚固剂与围岩界面处是拉拔过程中变形的主要部位。与此同时，大量数值模拟研究被应用于锚杆及锚固系统的支护效果评价中，如马双文等[7]使用ANSYS分析了锚杆横肋间距对锚固效果的影响；张辉等[8]采用ABAQUS研究了钻孔扩孔后的响应，结果表明，增大围压后扩孔直径萎缩，导致扩孔效果变差；马振乾等[9]研究认为，锚索间距是影响泥质顶板巷道稳定的最主要因素；贾后省等[10]采用FLAC3D研究二次采动影响下不同锚杆锚索支护方案对围岩的控制效果，并认为超前工作面5 m的范围在二次采动影响下，提高支护强度对于改善顶板塑性破坏及顶板下沉的控制效果有限。在利用FLAC进行的数值模拟中，较多研究人员采用Cable结构单元模拟锚杆，如在隧道、边坡、巷道等支护研究中均有此类案例[11-13]。在数值模拟中，时步控制对模拟结果具有重要影响，上述研究却未能披露如何控制锚杆的安装时机、锚杆安装后与围岩的相互作用时步、锚杆与钻孔的关系等细节。对于实际工程而言，锚杆安装需要经历钻孔、锚杆锚固、预紧等环节，对于拉拔试验而言，需要锚杆锚固稳定后才可进行拉拔，涉及的局部或系统稳定问题与数值模拟中达到的平衡状态是近似的，不同的平衡状态会直接影响到模拟结果，并进一步影响到对锚固效果的科学评价。因此，有必要对数值模拟中平衡状态对模拟结果的影响展开研究。本文利用FLAC3D软件建立锚杆的拉拔试验模型，结合实际工程及模拟中存在的几种不同平衡状态，探究锚杆拉拔数值模拟中平衡状态对模拟结果的影响。

1 基本模型的建立

结合室内锚杆拉拔试验，建立数值模拟模型，几何尺寸如图1所示，其中，圆柱体用以模拟围岩，圆柱体轴向沿Y轴，高2.2 m，直径0.4 m。采用文献[14]中的岩石属性，围岩参数为：密度2 600 kg/m3，杨氏模量52.39 GPa，泊松比0.29，内聚力50 MPa，摩擦角30°，抗拉强度10 MPa。利用Cable结构单元建立全长锚固的锚杆，锚杆位于圆柱体的中轴线，参数为：密度7 850 kg/m3，弹性模量200 GPa，直径22 mm，压缩和拉伸屈服强度均为260 kN。锚固剂参数为：单位长度内聚力2×106N/m，刚度4.5 GPa，摩擦角30°，锚固周长9.42×10-2m。锚杆长2.2 m，其中外露段长0.1 m，深入围岩2.1 m，沿轴向共划分为22个Cable结构单元，并由拉拔端开始编号1～22，坐标原点位于锚杆的第2结点。

图1 数值模拟模型几何尺寸示意图

模型底部固定，沿Y轴正向施加重力加速度10 m/s2。锚杆采用全长锚固，固定锚杆外露段端点，并在端点处施加恒定拉拔速度10-6m/步，方向沿Y轴负向。为了提高运算速度，拉拔试验开始前的所有计算中设定围岩为弹性模型，拉拔试验开始后的围岩本构关系为莫尔-库伦模型，同时在拉拔试验开始前清空模型的位移、速度、状态量。以锚杆的分段标号记录锚杆的轴向力演化，每隔10步记录1次数据。

2 初始平衡状态对锚固效果的影响

在数值模拟建模中，特别是工程尺度下建模时，由于模型尺寸较大，相对于工程体而言，锚杆是小物件。因此，在一些支护工程的锚杆建模中，开挖巷道后立即或运行指定计算时步安装锚杆，锚杆安装之后，往往立即或运算指定时步后继续执行开挖和支护循环，以此评价锚杆与围岩之间的相互作用。在室内试验及现场拉拔测试中，锚杆安装后往往需要等待几小时甚至数天以后才能进行拉拔试验，这是为了锚杆与锚固剂、围岩进一步相互作用，并达到新的锚固系统平衡。针对以上两种情况，该文设置两种模拟方案：方案1，围岩初始平衡后立即安装锚杆，安装锚杆后马上进行拉拔试验；方案2，围岩初始平衡后安装锚杆，安装锚杆后继续利用弹性模型计算,获得锚固系统新的平衡，在新平衡的基础上进行拉拔试验。其中，围岩的初始平衡指在重力作用下围岩的初始应力条件，初始平衡状态指上述两方案中锚杆被拉拔之前所对应的不同平衡状态，该文以FLAC3D默认的平衡收敛条件作为达到相应平衡状态的判断条件。

两种方案中锚杆轴向力演化如图2所示，以拉拔试验开始作为运算时步的起始位置。由两种方案的结果可知，从轴向力的演化趋势而言，两种方案的趋势基本相同。锚杆轴向力随着远离拉拔端而逐渐降低，在拉拔过程中，靠近拉拔端的位置轴向力保持较高水平，而远离拉拔端的杆体内轴向力则相对较小。轴向力达到屈服极限后(轴向力达到最大值并保持不变，即图2中曲线的水平段)，在锚杆-围岩锚固系统中开始出现锚固失效。距离拉拔端越远，轴向力越小，此部分锚杆杆体内的轴向力尚未达到杆体破坏的屈服强度，即使锚固系统出现失效后，远离拉拔端的锚杆杆体依然保持完整，但因锚杆已经出现破坏，该部分完整的杆体也不能继续承受拉拔载荷。对于此问题，有学者提出通过Fish函数实现锚杆破断后受力归零[15-16]。由于此次研究仅针对锚固平衡状态展开，故无需采用该方法。

图2 方案1与方案2的锚杆轴向力演化图

以拉拔试验开始时作为两种方案的起始位置对相同监测点的轴向力做差，图3为方案1与方案2轴向力差值的演化图，两种方案中轴向力均为正值，故差值为正数，正数越大，说明方案1比方案2中的轴向力大；差值结果为负数时，负数的绝对值越大，说明方案1比方案2中的轴向力小；差值接近于0时，说明两方案中的轴向力差别不明显。

图3 方案1与方案2的轴向力差值演化图

由图3可知，方案1与方案2中锚杆轴向力的差别主要体现在靠近拉拔端的锚杆结构单元上，锚杆破断前，方案1中拉拔端的轴向力高于方案2的，而深入围岩杆体中的轴向力差别则不明显。距离拉拔端越远，轴向力越小，造成两种方案远离拉拔端的差值可以忽略不计。方案1的拉拔开始于19 310时步，方案2的拉拔开始于37 920时步，锚杆从拉拔到破坏仅需要不足1 000个时步。在一些研究中，试图通过在锚杆安装后运行一定时步以模拟锚杆与围岩的相互作用，但从锚杆与围岩重新达到平衡所需时步来看，这种效果并不显著，甚至在包含围压的条件下，有可能锚杆安装后就受到类似的拉拔作用，人为指定一定运算时步，还有可能造成锚杆提前进入拉拔状态，锚杆安装后受围岩的应力状态影响更为明显，因为围岩对锚杆施加了类似的拉拔作用。对于既要考虑开挖又要考虑支护的模拟，锚杆安装后与围岩的平衡问题很难克服。对于单根锚杆的室内拉拔试验，有必要考虑锚杆与围岩的二次平衡，以便更真实地反映整个锚固系统各物理量的演变特性。

3 钻孔内锚固平衡对模拟结果的影响

3.1 模型建立方式概述

另一个值得注意的问题是FLAC3D中采用Cable结构单元建立锚杆时的钻孔平衡问题。在实际工程或试验中，一般要经历围岩开挖自稳平衡、钻孔后的平衡、安装锚杆后的平衡等几个过程。但在FLAC3D的研究中，较少见到对相关细节的呈现。为此，设计3个模型进行对比研究：(1)模型1，与上述方案2相同，考虑围岩的初始平衡和安装锚杆后的锚杆-围岩平衡，不考虑模型中的钻孔问题；(2)模型2，围岩初始平衡后，在圆柱体的实体网格中开挖直径30 mm的钻孔，并计算钻孔平衡，然后安装锚杆计算锚杆与围岩的二次平衡，最后进行拉拔试验；(3)模型3，利用模型2空心钻孔后的平衡模型，将原本空心钻孔的位置建立锚固剂的实体网格，参照水泥砂浆锚固剂特性，赋予该部分参数，即密度2 100 kg/m3，弹性模量25 GPa，泊松比0.25，内聚力5 MPa，摩擦角30°，抗拉强度5 MPa。安装锚杆并计算新的平衡，再进行拉拔试验。

3.2 结果分析

3.2.1 锚杆轴向力的演化差异

图4为模型2与模型3的锚杆轴向力演化图，模型1的轴向力演化过程可参考图2(b)。综合图2和图4的结果可知，模型3与模型1模拟效果类似，但模型3提高了远离拉拔端锚杆结构单元的受力效果。与模型1中锚杆的受力部分主要以拉拔端为主相比，虽然模型3中依然是距拉拔端越近，轴向力越大，但轴向力并没有急速下降，从而可以避免模型1中的局部应力集中和应力差较大的问题。对于采用水泥砂浆类锚固剂，可考虑钻孔开挖后安装锚固剂、锚杆与围岩的新平衡状态可改善锚杆杆体的受力。模型2中，由于钻孔位置被挖空，锚杆杆体处于临空状态，因此,各监测点的轴向力变化几乎同步。即拉拔过程中，模型2中锚杆杆体是独立于围岩的，不能反映出围岩与锚杆的相互作用。因此，这种方案虽然考虑了实际工程中的钻孔，但是却并不能模拟出锚杆-围岩锚固系统的特性。原因在于，FLAC3D中锚杆特性由杆体、锚固剂、node、link共同构建，特别是node和link两元素，它们通过继承杆体、锚固剂的特性与围岩的网格单元发生作用，并将这种作用反馈给锚杆杆体[17]。模型2中的钻孔为空心，造成锚杆安装后无法建立有效的link，也就无法与围岩发生相互作用。换言之，FLAC3D中Cable结构单元的锚固剂并不是一种实体，而是需要依托node和link发挥作用，即便对锚杆结构单元中的锚固剂赋予了参数，但如果锚固剂范围内没有实体网格单元，也不能发挥出锚杆的锚固作用。因此，某种程度上可以认为模型2是错误的。另一方面，该文采用了全长锚固，工程中很多锚杆为非全长锚固。模型2也说明，对于非全长锚固，如果锚杆周围没有实体单元，则锚杆自由段内的轴向力演化基本一致，这与以往研究中只将自由段锚固参数设置到较低水平是有所差别的。因此，要想准确模拟出自由段轴力演化基本一致的效果(即模拟出理论中假设自由段中锚杆受力处处相等的效果)，可将锚杆自由段部分的实体单元设置为null，此时，相比于锚固段，自由段整体处于较高的拉拔力作用下，当破坏位置接近自由段与锚固段的交界面时，该部位的自由段杆体易发生应力集中并导致破坏。因此，工程中进行锚杆支护设计时，应避免锚固段与自由段的交界处于节理、裂隙等接触面之上。

图4 模型2与模型3中的轴向力演化图

3.2.2 沿锚杆轴向围岩的位移演化

图5为模型1和模型3相同位置围岩中沿锚杆轴向位移截面图(截取x=0平面)。由图5可知：模型1和模型3在围岩的轴向位移演化方面存在明显差异。模型1中的变形已经扩展到围岩当中，且距拉拔端越近，位移量越大。而在模型3中，位移主要发生在钻孔内，变形并未向围岩扩展，越靠近锚杆的轴心，位移量越大，超出锚固剂的范围后，围岩内的轴向位移量几乎为0。造成这一现象的可能原因为模型3中采用了实体网格作为锚固剂重新填充了钻孔，实体锚固剂的参数与围岩和锚杆相比，比较软弱，锚杆受到拉拔作用后，与软弱的锚固剂相互作用，锚固剂与围岩之间没有形成良好的锚固体系，此时锚杆的锚固作用主要发生在锚固剂之中，而对围岩的控制效果有限。另一方面也说明，如果锚固剂与围岩相比，锚固作用足够强，那么锚固系统的整体性将更好，锚杆也将更好地发挥协调围岩变形的作用，而不会直接从锚固剂中滑脱。

图5 围岩中沿锚杆轴向位移演化结果

为验证上述猜想，设置3种强化方案。强化方案1，将模型3中钻孔内锚固剂实体单元的部分参数进行强化，强化参数为：弹性模量80 GPa，内聚力80 MPa，抗拉强度20 MPa；强化方案2,将模型3中Cable结构单元的锚固剂参数强化，具体为：单位长度的内聚力20×106N/m，刚度80 GPa；强化方案3，对钻孔内锚固剂实体单元和Cable的锚固剂参数均强化，数值与前两种方案相同。除上述参数外，其他参数不变。

对模拟结果进行分析，发现强化方案1与强化方案3获得了与图5(a)类似的效果，变形均已扩展到围岩中，特别是强化方案3，变形向围岩中扩展的层次性更明显，其位移的演化结果如图6所示。而强化方案2的模拟结果与图5(b)中呈现的结果类似，变形主要集中在锚固剂中，并没有向围岩扩展。上述3种方案均对锚固剂进行了不同方式的强化，但从模拟结果而言，使用Cable结构单元时，如果不改变围岩属性，仅改变Cable的锚固参数，在拉拔条件下获得的锚固效果是类似的。对钻孔范围内的锚固剂实体单元进行强化，能更好地实现锚杆与围岩的协调运动。由上述结果对比可知，当锚固剂足够强时，锚杆受到拉拔作用后会通过锚固剂将拉拔作用传递给围岩，从而带动围岩协同运动。当更多的围岩能够参与到变形破坏过程时，整个锚固系统抵御破坏失效的整体能力就得到了提升，从而提高了锚固强度。但Cable锚固效果的好坏还有赖于与其临近的围岩，Cable周围围岩越坚硬，越有利于锚杆与围岩发挥出系统性的锚固作用。

图6 强化方案3中沿锚杆轴向位移演化结果

3.2.3 钻孔内锚固平衡状态对模拟结果的影响

钻孔内锚固平衡状态对模拟结果的影响主要体现在模拟时步上，以拉拔起始时步为例说明，模型1～3、强化模型1～3的拉拔起始时步分别为37 920，33 740，49 070，43 830，49 390，43 940，锚杆进入屈服状态所需的运算时步分别为910，7 200，3 900，700，4 170，540，结合前述分析可知：根据真实情况建立空心钻孔时，运算时步最长且获得的锚固效果最差(模型2)；采用实体单元模拟锚固剂充填钻孔后，实体单元属性要比Cable本身锚固剂参数的影响更明显，实体单元的属性越强，越有利于变形在拉拔端部向围岩层次性扩展(图6)，但同时也会造成位移加载条件下锚杆失效加速，如模型1、强化模型1、强化模型3中锚固剂的属性逐渐增强，相应的模型达到拉拔失效所需的时步逐渐减少；采用实体单元模拟锚固剂充填钻孔后，若仅改变Cable锚固剂的参数，则锚固剂属性越强，达到锚杆失效所需的时步越多，如强化模型2的时步大于模型3的；从相对时步而言，钻孔内锚固平衡状态对模拟结果的影响较小，但从绝对时步而言，可考虑钻孔内锚固平衡后运算时步会显著增加。因此，若考虑围岩应力释放和锚杆滞后发挥支护作用，则不能忽略钻孔内锚固平衡状态，若仅评价锚杆参数改变对支护效果的最终影响，则无需考虑钻孔内锚固平衡问题。

4 模拟结果的讨论

4.1 初始平衡问题

从前述研究可以看出，在FLAC3D中采用Cable进行锚杆模拟时，是否运算锚杆与围岩的平衡会影响到锚杆破坏前后的力学状态，从而影响分析的准确性。但对于工程建模而言，比如煤矿井下的巷道掘进和支护，锚杆只是该工程中的一小部分，特别是计算模型往往还要受到地应力影响，实际上不可能每次安装锚杆都在FLAC3D中进行锚杆与围岩之间的平衡运算，这将导致锚固系统失效分析并不精确。当安装锚杆后，企图通过人为指定一定运算时步模拟这种平衡，实际上效果很有限。最主要的原因在于指定的运算时步远达不到平衡所需的时步，即便接近了平衡状态，也很难在一个工程模型中将模型的位移、速度等平衡后的状态变量归零。诚然，可以通过物理参数的相对变化了解一些参数前后的变化，但由于FLAC3D中很多参数是矢量，网格之间存在交互影响，这种相对变化也未必足够精确。因此，对于复杂工程中的锚杆支护问题而言，在静力学分析状态下，FLAC3D更像是一种趋势性分析，与现实工程一一对应的精确定量化分析目前还非常困难，而在拉拔试验模拟中，则有必要考虑平衡状态对锚杆性能的影响。

4.2 工程围压问题

该文研究中，除重力外，围岩周围并没有设定围压，因此获得的位移量相对较小。但对于地下工程，地下岩体受到地应力的影响，地应力由FLAC3D中实体网格的运动提供动力，有可能改变应力、位移的数量级。已有研究表明[18]，一定围压范围内，围压增大后，锚杆的拉拔力随之增大。实际上，对于大部分室内拉拔试验而言，围岩被固定在箱体中，箱体施加于围岩的围压很小，从而导致数值模拟和物理拉拔试验中水平位移比较小。现场地下工程原位的拉拔试验情况更为复杂，即便是瓦斯抽采，也会造成锚杆拉拔特性的改变[19]。常规的围压研究一般以假三轴围压模式为主，未来还有必要结合地应力分布，进一步深入研究非等压、非均匀分布围岩等对锚杆性能的影响。

4.3 锚固剂问题

该文研究表明，FLAC3D中锚杆支护性能的发挥与锚固剂参数密不可分，实际工程中，锚固剂也是非常重要的参数。原因在于锚杆与围岩组成的锚固系统中，锚固剂是连接锚杆与围岩的媒介，而在锚固剂与围岩、锚杆与锚固剂之间又形成接触面，接触面的强度与锚固剂的性能也有重要关系。康红普等[20]调研表明，树脂锚杆的破坏主要发生在各种交界面。因此，锚固系统的性能不仅要关心锚杆杆体的性能，而且更要注重锚固剂的性能和施工质量，以使锚固剂能充分发挥出能效，避免发生锚杆从锚固剂中脱锚、锚固剂破碎、锚固剂与围岩接触不实的现象。

另一方面，Itasca公司提出可利用十分之一因素法确定数值模拟中锚固剂的刚度kg[17]。对于煤矿而言，锚杆直径通常为16～25 mm，一般钻孔直径与锚杆直径之差控制在6～10 mm[21]。对于直径为16～25 mm的锚杆，使用十分之一因素法后，kg为剪切模量G的1.3～2.9倍。该数值与不使用十分之一因素法相比，存在10～27.5倍的误差。尤其是G，其数量级一般是以109计算，如此大的数量级，也会产生十分之一因素法误差的放大效应。因此，FLAC3D中是否采用十分之一因素法确定锚固剂的刚度，应结合研究的具体问题确定，应该多进行一些对比研究，以观察不同参数组合在FLAC3D中的响应。

4.4 锚固平衡问题

近年来，对锚杆性能的评价主要集中于锚固系统平衡后的承载性能评价，但对于锚固平衡状态以及锚固平衡状态带来的后续影响研究则较少，以至于对于非全长锚固的锚杆，锚固后围岩中的应力场分布形式、分布范围、分布时效等方面的研究还不够深入,特别是工程中分区破裂、巷道不同深度应力非均匀分布、不同围岩属性等情况，目前的研究仍相对简化。加之FLAC3D中无法再现锚固剂的成型、固化、黏结等过程，对于锚固平衡问题，未来有必要在物理试验精确监测的基础上，在数值模拟中进一步深入研究。

5 结 论

(1)从轴向应力的演化趋势而言，是否计算锚杆安装后锚杆与围岩的二次平衡并不影响锚杆拉拔试验中轴向力的演化趋势，但会影响到轴向力演化的数量级和锚固系统失效的时间。对于既要考虑开挖又要考虑支护的模拟而言，试图通过锚杆安装后运行一定时步，以模拟锚杆与围岩的新的平衡效果并不显著，锚杆安装后受围岩应力状态的影响更为明显。

(2)FLAC3D中Cable结构单元的锚固剂并不是实体，需要依托node和link发挥作用，因此，锚杆锚固段不能建立在null的空间中。想要准确模拟出锚杆自由段轴力处处相等的效果，可将与锚杆自由段部分相接触的实体单元网格设置为null。在进行锚杆支护设计时，应避免锚固段与自由段的交界处于节理、裂隙等接触面上。

(3)对于相对时步，钻孔内锚固平衡对模拟结果的影响较小，但对于绝对时步，考虑钻孔内锚固平衡后运算时步会显著增加，因此，若考虑围岩应力释放和锚杆滞后发挥支护作用，则不能忽略钻孔内锚固平衡状态，若仅评价锚杆参数改变对支护效果的最终影响，则无需考虑钻孔内锚固平衡问题。与锚杆接触的实体单元属性要比Cable本身锚固剂参数对模拟结果的影响更为显著。