本发明属于岩石力学与岩石工程领域，具体涉及有限元-离散元耦合(fdem)恒阻大变形锚杆隧道围岩加固数值模拟方法。
背景技术：
在深埋软岩地层中掘进隧道后将引起围岩的大变形，并表现出变形速率快、变形量大、变形持续时间长等特征。为了控制围岩的大变形，实际工程中经常采用锚网喷 衬砌联合支护措施。针对传统锚杆存在延伸率低、难以适用于软岩大变形加固控制等不足，何满潮等人发明了“恒阻大变形锚杆”(申请号：cn201010196197.2)，通过在自由端端头设置恒阻大变形装置的方法，使得该类锚杆具有延伸率大、阻力恒定等特性，实现了软岩隧道大变形的“让压”支护，因而在诸多深埋软岩隧(巷)道工程加固中得到了成功运用，如兰渝铁路木寨岭隧道、新疆沙吉海煤矿巷道和龙口北皂煤矿巷道等。
为了探究恒阻大变形锚杆力学特性并获取相关力学参数，可采用室内拉伸试验的方法；然而，为研究恒阻大变形锚杆对深埋软弱隧道围岩的加固机理并为得到合适的锚杆力学参数，往往需要进行现场测试研究。对现场测试而言，存在测试周期漫长、测试成本高昂、对隧道工程正常施工造成干扰等不足，且监测结果仅对本测试段有效，当围岩和赋存环境发生了较大变化时，又需要进行新的现场监测，难以得到普适性的、规律性的和机理性的研究成果。因此，可采用数值模拟方法进行恒阻大变形锚杆围岩加固的研究。
一般地，在岩石力学与岩石工程领域，数值模拟可分为连续性方法、非连续性方法和耦合方法。耦合方法中的fdem方法由于能够模拟岩体弹塑性连续变形至裂纹萌生-扩展-交汇贯通非连续变形全过程，再至断裂破碎块体间的接触挤压效应，使其在模拟隧道围岩破裂碎胀大变形失稳灾变过程及锚杆加固方面具有天然优势。在当前的fdem锚杆加固研究中，通过设定无厚度线性单元的方法，锚杆连续穿插三角形单元和四边形节理单元而不占据空间体积，锚杆与三角形单元块体绑定，三角形内的锚杆随三角形单元的变形一同发生弹性变形；锚杆的屈服与断裂通过四边形节理单元内的杆体单元来反映。
上述锚杆的表征方法只能模拟传统锚杆的全长锚固，而无法实现端部锚固及预应力锚固，也无法模拟新型恒阻大变形锚杆的锚固，且忽略了锚固剂的力学作用和锚固剂-锚杆和锚固剂-围岩间的粘结作用。因此，本发明基于实体建模的方法，提出设置特定的功能型四边形节理单元的方法，实现恒阻大变形锚杆在fdem中的表征；提出锚固剂力学本构方程、锚固剂-锚杆界面本构方程和锚固剂-岩石界面本构方程获取方法，为恒阻大变形锚杆围岩加固机理fdem数值模拟提供研究基础。
技术实现要素：
为了解决上述技术问题，本发明的目的是提出一种有限元-离散元耦合(fdem)恒阻大变形锚杆隧道围岩加固数值模拟方法，实现不同工况下恒阻大变形锚杆加固效果的有效预测，并为揭示恒阻大变形锚杆加固机理研究提供可靠的数值模拟手段。
本发明提供的方案如下：
一种有限元-离散元耦合恒阻大变形锚杆隧道围岩加固数值模拟方法，包括以下步骤：
(1)采用实体建模的方法建立锚杆、锚固剂和钢带数值模型；其中，锚杆、锚固剂分别采用结构化网格，而钢带采用非结构化网格；对锚杆而言，在自由端端头设置功能型四边形节理单元为恒阻大变形单元，并赋予相应的力学参数和本构方程；而除该单元外的其余三角形单元和四边形单元均采用线弹性本构方程，如图1所示，并根据锚杆在隧道中的纵向安装间距计算锚杆等效弹性模量和等效剪切模量；
(2)通过锚固剂力学试验和锚固剂-锚杆及锚固剂-岩石界面力学试验，获取相应的力学本构方程，实现各部件的模拟；
(3)建立锚杆、锚固剂和钢带的模型，通过破裂节理单元粘结修复算法，实现锚杆滞后安装模拟。
进一步，所述锚杆力学实验采用陶志刚等人提出的实验方法(《恒阻大变形锚杆/索力学特性数值模拟研究》.陶志刚，赵帅，张明旭，孟志刚，何满潮，赖豪.采矿与安全工程学报，2018年1月第1期)。
进一步，利用fdem中材料被划分为三角形单元 初始无厚度四边形节理单元的特性，采用结构化网格建立锚杆实体模型，在锚杆自由端端头设置功能型四边形节理单元实现锚杆恒阻大变形功效，如图1所示，并赋予室内试验得到的荷载-位移本构方程：
式中，k1、k2分别为锚杆刚度，u1、u2、umax分别为不同阶段位移拐点，pave为恒阻阶段荷载，如图2所示，δu为锚杆位移。在求得荷载后，恒阻大变形单元的轴向应力为：
σ0＝p·a(2)
式中，σ0为恒阻单元的轴力，a为锚杆横截面积。
进一步，对于恒阻大变形锚杆而言，通过室内拉伸试验可知(图2)，它经历了线弹性加载阶段、恒阻变形阶段和破坏阶段，且破坏发生于恒阻装置内，锚杆杆体全程仅发生弹性变形，因此，除了恒阻大变形单元，其余三角形单元和四边形节理单元均仅发生弹性变形，其本构方程为：
σ＝e′·ε(3)
式中，σ、e′、ε分别为除恒阻大变形单元外的锚杆轴力、锚杆等效弹性模量和轴向应变。在二维平面应变状态下，锚杆杆体厚度为1.0m，实际中的锚杆厚度为其直径，因此，等效弹性模量e′为：
式中，d为锚杆杆体直径，l为锚杆纵向安装间距，单位均为m；e为锚杆实际弹性模量，可取为210gpa。锚杆的剪切本构方程(包括恒阻力单元)采用如下形式：
τ＝g′·γ(5)
式中，τ、g′、γ分别为剪应力、等效剪切模量和剪切应变。等效剪切模量为：
式中，g为锚杆实际剪切模量。
进一步，锚固剂拉伸与剪切本构方程获取的方法如下：在锚固端，通过锚固剂将锚杆固定，锚固剂的拉伸与剪切本构方程获取方法如下：将锚固剂用水稀释，采用试验用砂浆搅拌机将稀释后的锚固剂充分搅拌，浇筑至50mm×100mm和50mm×50mm圆柱形磨具中，浇筑过程中不断振动，确保内部不再有气泡产生，保证浇筑质量。将浇筑完毕的试样在养护箱内养护28d，养护完毕后开展直接拉伸和变角剪切试验，将试验结果采用王中伟等人提出的公式进行拟合(黏结裂纹模型中的岩石软化曲线确定方法及装置，申请号：cn202010263316.5)，得到锚固剂拉伸与剪切本构方程，如下：
式中，σ、ε为任意时刻应力、应变；ft、εr为峰值应力、应变；c1，c2，c3，c4为拟合参数。
进一步，锚杆-锚固剂界面本构方程获取方法如下：
将锚杆-锚固剂-高强混凝土进行浇筑，混凝土采用高标号水泥，预留锚杆孔，待混凝土养护7d后贯入稀释后的锚固剂同时插入锚杆，轻微振动确保锚固剂与混凝土及锚杆充分接触粘结、且锚固剂内部无气泡；在养护箱内养护28d；
浇筑后的结构(图3所示)为锚杆通过锚固剂浇筑于高强混凝土中，高强混凝土通过固定螺杆连接到固定承压板上；
养护完成以后进行锚杆拉拔试验，得到锚杆-锚固剂界面剪切应力-位移本构方程。锚杆-锚固剂界面拉伸本构方程采用上述得到的锚固剂拉伸本构方程。
进一步，锚固剂-岩石界面处的拉伸及剪切本构方程的获得方法如下：将岩样分别加工成70mm×100mm的圆柱形试样和50mm×70mm×70mm的长方体试样，如图4-5所示。采用浇筑的方式制作相同尺度的锚固剂试样，在浇筑过程中，将制作好的岩样放入磨具底部，使得锚固剂能够与岩样得到充分接触粘结，在浇筑过程中轻微振动，保证锚固剂浇筑质量。待锚固剂-岩石样品达到终凝后拆除磨具，放入养护箱养护28d，在交界面处加工如图4-5所示楔形缺口，确保样品在加载过程中从交界面处起裂。分别进行室内直接拉伸和变角剪切试验，与锚固剂的方法相似，得到锚固剂-岩石界面处的拉伸及剪切本构方程。
进一步，采用gmsh软件建立如图6所示的模型，
其中，恒阻大变形锚杆垂直于隧道表面，深入围岩内部一侧的为锚杆锚固端，通过锚固剂固定；临近隧道表面的为锚杆自由端，并在靠近隧道表面的自由端指定一功能型节理单元为恒阻大变形单元；在锚杆与围岩间通过安装钢带分散锚杆自由端端头的集中应力；
锚杆、锚固剂和钢带采用实体建模的方式，钢带采用与锚杆相同的本构方程和力学参数，二者以接触效用发挥相互作用。锚杆采用结构化网格，使得恒阻大变形单元垂直于隧道轮廓线，锚杆上部为锚固端(通过锚固剂固定)、下部为自由端，在自由端设置功能型四边形节理单元为恒阻大变形单元，锚杆自由端和围岩间通过钢带分散锚杆端头的集中应力。锚固剂采用结构化网格，钢带采用非结构化网格。在锚杆激活前，锚杆、锚固剂与围岩属性相同，为围岩的一部分，一同发生变形与破裂；锚固剂与隧道核心材料相同，一同发生软化卸荷。当锚杆激活后，锚杆、锚固剂和钢带应力、位移和速度进行清零处理，将发生破坏的节理单元粘结修复重新发挥作用，如图7所示，赋予各自的本构方程和力学参数，开始发挥加固作用，节理单元在锚杆激活前发生的位移不计入在内，即在锚杆激活后，锚杆和锚固剂的拉伸/压缩与剪切位移为：
式中，o、o2、o1分别为锚杆激活后的拉伸/压缩位移、总位移和锚杆激活前的拉伸/压缩位移；s、s2和s1为对应的剪切位移。
本发明的有益效果：
与现有的fdem锚杆加固数值模拟方法相比，本发明基于实体建模的方法，设置功能型节理单元，实现恒阻大变形锚杆在fdem中的表征；对锚杆其余单元采用线弹性本构方程，提出了等效弹性模量和等效剪切模量计算方法；提出锚固剂、锚固剂-锚杆界面和锚固剂-岩石界面本构方程获取方法，顺利实现恒阻大变形锚杆加固fdem数值模拟；采用破裂节理单元(仅对锚杆和锚固剂而言)粘结修复算法，实现锚杆滞后安装模拟；本发明为恒阻大变形锚杆围岩加固效果预测及揭示恒阻大变形锚杆加固力学机理提供数值模拟手段。
附图说明
图1为恒阻大变形锚杆在fdem中的表征原理；
图2为实验和模拟得到的恒阻大变形锚杆荷载-位移曲线；
图3为锚杆拉拔试验装置示意图；
图4为锚固剂-岩石界面拉伸示意图；
图5为锚固剂-岩石界面剪切示意图；
图6为恒阻大变形锚杆隧道支护模拟数值模型图；
图7为破裂节理单元粘结修复示意图；
图8为恒阻大变形隧道加固fdem模拟效果。
具体实施方式
以下结合附图对本发明涉及有限元-离散元耦合(fdem)恒阻大变形锚杆隧道围岩加固数值模拟方法进行详细说明，本发明的内容完全不限于此。
实施例
恒阻大变形锚杆隧道围岩加固数值模拟方法如下：
步骤1.开展恒阻大变形锚杆拉伸试验，获取恒阻大变形锚杆力学参数及本构方程，或采用公开发表文献中已获取的本构方程，以《恒阻大变形锚杆/索力学特性数值模拟研究》(陶志刚，赵帅，张明旭，孟志刚，何满潮，赖豪.采矿与安全工程学报，2018年1月第1期)得到的拉伸-位移曲线为例，建立如图1所示的恒阻大变形锚杆实体模型；该拉伸本构方程如下：
式中，k1、k2分别为锚杆刚度，u1、u2、umax分别为不同阶段位移拐点，pave为恒阻阶段荷载，如图2所示，δu为锚杆位移。在图2中，k1＝12.8kn/mm，k2＝2.58kn/mm，pave＝128.9kn，u1＝10mm，u2＝590mm，umax＝640mm。
步骤2.本实施例中，锚杆纵向安装间距设置为l＝1.0m，锚杆直径d＝22mm，弹性模量e＝210gpa，剪切模量g＝84gpa，因此本实施例中锚杆等效弹性模量和等效剪切模量为：
式中，e′、g′分别为锚杆的等效弹性模量和等效剪切模量。根据上式得到等效弹性模量和等效剪切模量分别为e′＝4.62gpa、g′＝1.848gpa。锚杆采用如下线弹性本构方程：
σ＝e′·ε(4)
τ＝g′·γ(5)
式中，σ、ε分别为锚杆轴向应力和轴向应变；τ、γ分别为剪应力和剪切应变。此外，对于恒阻大变形单元而言，其剪切本构方程亦采用式(5)的形式，且g′＝1.848gpa。
步骤3.将锚固剂用水稀释，采用试验用砂浆搅拌机将稀释后的锚固剂充分搅拌，浇筑至50mm×100mm和50mm×50mm圆柱形磨具中，浇筑过程中不断振动，确保内部不再有气泡产生，保证浇筑质量。将浇筑完毕的试样在养护箱内养护28d，养护完毕后开展直接拉伸和变角剪切试验，将试验结果采用王中伟等人提出的公式进行拟合(黏结裂纹模型中的岩石软化曲线确定方法及装置，申请号：cn202010263316.5)，得到锚固剂拉伸与剪切本构方程，如下：
式中，σ、ε为任意时刻应力、应变；ft、εr为峰值应力、应变；c1，c2，c3，c4为拟合参数。
步骤4.将锚杆-锚固剂-高强混凝土浇筑成如图3所示试样模型，混凝土采用高标号水泥以确保混凝土不产生开裂，预留锚杆孔，待混凝土养护7d后贯入稀释后的锚固剂同时插入锚杆，轻微振动确保锚固剂与混凝土及锚杆充分接触粘结、且锚固剂内部无气泡。在养护箱内养护28d，进行锚杆拉拔试验，得到锚杆-锚固剂界面剪切应力-位移本构方程。锚杆-锚固剂界面拉伸本构方程采用步骤3得到的锚固剂拉伸本构方程。
步骤5.将岩样分别加工成70mm×100mm的圆柱形试样和50mm×70mm×70mm的长方体试样，如图4-5所示。采用浇筑的方式制作相同尺度的锚固剂试样，在浇筑过程中，将制作好的岩样放入磨具底部，使得锚固剂能够与岩样得到充分接触粘结，在浇筑过程中轻微振动，保证锚固剂浇筑质量。待锚固剂-岩石样品达到终凝后拆除磨具，放入养护箱养护28d，在交界面处加工如图4-5所示缺口，确保样品在加载过程中从交界面处起裂。分别进行室内直接拉伸和变角剪切试验，与步骤3所述方法相似，得到锚固剂-岩石界面处的拉伸及剪切本构方程。
步骤6.采用gmsh软件建立如图6所示的恒阻大变形锚杆隧道围岩加固数值模型，以矩形隧道为例，隧道宽6.0m、高3.0m，锚杆长4.0m，直径22mm，横向安装间距0.9m、纵向安装间距1.0m，锚固端长度1.0m、自由端长度3.0m，钢带厚0.122m，水平和垂直地应力均为36mpa。为了对比加固效果，仅对隧道顶板进行了恒阻大变形锚杆加固处理，围岩力学参数表1所示：
表1围岩参数
步骤7.在锚杆激活前，锚杆、锚固剂力学参数和本构方程与围岩相同，且与围岩一同发生变形及破裂；钢带与核心材料属性相同，与核心材料一起发生卸荷软化。当锚杆激活后，锚杆、锚固剂已破裂的节理单元进行粘结修复处理(重新赋予属性，进入节理单元力学计算)，前期发生的位移、应力和速度清零，锚杆激活后的位移为：
式中，o、o2、o1分别为锚杆激活后的拉伸/压缩位移、总位移和锚杆激活前的拉伸/压缩位移；s、s2和s1为对应的剪切位移，如图7所示。
通过上述步骤的模拟，得到了如图8所示的隧道围岩恒阻大变形锚杆加固模拟结果，表明恒阻大变形锚杆对隧道顶板加固效果明显。根据本发明提出的方法，实现了恒阻大变形锚杆在fdem中的表征，反映了恒阻大变形锚杆力学特性；建立了破裂单元粘结修复算法，实现了锚杆滞后安装模拟。
以上所述，仅为本发明较佳的具体实施方式，但本发明保护的范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内所做的任何修改，等同替换和改进等，均应包含在发明的保护范围之内。
1.一种有限元-离散元耦合(fdem)恒阻大变形锚杆隧道围岩加固数值模拟方法，其特征在于，包括以下步骤：
(1)采用实体建模的方法建立锚杆、锚固剂和钢带数值模型；其中，锚杆、锚固剂分别采用结构化网格，而钢带采用非结构化网格；对锚杆而言，在自由端端头设置功能型四边形节理单元为恒阻大变形单元，并赋予相应的力学参数和本构方程；而除该单元外的其余三角形单元和四边形单元均采用线弹性本构方程，并根据锚杆在隧道中的纵向安装间距计算锚杆等效弹性模量和等效剪切模量；
(2)通过锚固剂力学试验和锚固剂-锚杆及锚固剂-岩石界面力学试验，获取相应的力学本构方程，实现各部件的模拟；
(3)建立锚杆、锚固剂和钢带的模型，通过破裂节理单元粘结修复算法，实现锚杆滞后安装模拟。
2.根据权利要求1所述的方法，其特征在于：所述锚杆力学实验采用陶志刚等人提出的实验方法(《恒阻大变形锚杆/索力学特性数值模拟研究》.陶志刚，赵帅，张明旭，孟志刚，何满潮，赖豪.采矿与安全工程学报，2018年1月第1期)。
3.根据权利要求1所述的方法，其特征在于：通过设置功能型四边形节理单元的方法，实现恒阻大变形锚杆的数值表征，并将试验得到的本构方程和参数赋予该恒阻大变形单元，如下：
式中，k1、k2分别为锚杆刚度，u1、u2、umax分别为不同阶段位移拐点，pave为恒阻阶段荷载，δu为锚杆位移。
4.根据权利要求1所述的方法，其特征在于：所述锚杆除恒阻大变形单元外的其余单元采用线弹性本构方程，如下：
σ＝e′·ε
式中，σ、e′、ε分别为除恒阻大变形单元外的锚杆轴力、锚杆等效弹性模量和轴向应变。在二维平面应变状态下，锚杆杆体厚度为1.0m，实际中的锚杆厚度为其直径，因此，等效弹性模量e′为：
式中，d为锚杆杆体直径，l为锚杆纵向安装间距，单位均为m；e为锚杆实际弹性模量，可取为210gpa。锚杆的剪切本构方程(包括恒阻力单元)采用如下形式：
τ＝g′·γ
式中，τ、g′、γ分别为剪应力、等效剪切模量和剪切应变。等效剪切模量为：
式中，g为锚杆实际剪切模量。
5.根据权利要求1所述的方法，其特征在于：所述步骤(2)中，所述锚固剂力学试验的方法如下：将锚固剂用水稀释，采用试验用砂浆搅拌机将稀释后的锚固剂充分搅拌进行浇筑，浇筑成圆柱形试件，然后对试样进行养护，养护完毕后开展直接拉伸和变角剪切试验。
6.根据权利要求1所述的方法，其特征在于：所述步骤(2)中，所述锚固剂-锚杆力学试验方法为：
将锚杆-锚固剂-高强混凝土进行浇筑，混凝土采用高标号水泥，预留锚杆孔，待混凝土养护7d后贯入稀释后的锚固剂同时插入锚杆，轻微振动确保锚固剂与混凝土及锚杆充分接触粘结、且锚固剂内部无气泡；在养护箱内养护28d；
浇筑后的结构为锚杆通过锚固剂浇筑于高强混凝土中，高强混凝土通过固定螺杆连接到固定承压板上；
进行锚杆拉拔试验，得到锚杆-锚固剂界面剪切应力-位移本构方程。
7.根据权利要求1所述的方法，其特征在于：所述步骤(2)中，所述锚固剂-岩石界面学试验方法为：将岩样和锚固剂分别加工成圆柱形试样和长方体试样，在浇筑过程中，将制作好的岩样放入磨具底部，使得锚固剂能够与岩样得到充分接触粘结，终凝后得到锚固剂-岩石界面结构；然后放入养护箱养护28d；再将接触界面加工出楔形缺口以确保样品从界面处起裂；最后分别进行室内直接拉伸和变角剪切试验。
8.根据权利要求1所述的方法，其特征在于：所述步骤(2)中，锚固剂力学本构方程为：
式中，σ、ε为任意时刻应力、应变；ft、εr为峰值应力、应变；c1，c2，c3，c4为拟合参数。
9.根据权利要求1所述的方法，其特征在于，所述步骤(3)中，建模的模型结构如下：
恒阻大变形锚杆垂直于隧道表面，深入围岩内部一侧的为锚杆锚固端，通过锚固剂固定；临近隧道表面的为锚杆自由端，并在靠近隧道表面的自由端指定一功能型节理单元为恒阻大变形单元；在锚杆与围岩间通过安装钢带分散锚杆自由端端头的集中应力。
10.根据权利要求1所述的方法，其特征在于：所述步骤(3)中，在锚杆激活前，锚杆、锚固剂为围岩的一部分，与其他围岩发生变形与破裂；钢带为核心材料的一部分，随其他核心材料一同卸荷软化；
当锚杆激活后，将锚杆、锚固剂和钢带单元应力、位移和速度清零，赋予各自的力学本构方程及参数；破裂的锚杆和锚固剂节理单元重新修复，发挥节理单元的粘结作用，修复后的节理单元位移算法如下：
式中，o、o2、o1分别为锚杆激活后的拉伸/压缩位移、总位移和锚杆激活前的拉伸/压缩位移；s、s2和s1为对应的剪切位移。
技术总结