基于刚体动力学-离散元耦合数值方法的复合地层盾构滚刀响应

徐汪豪， 姚志刚， 方 勇， 王飞阳， 周 健

(1.浙江杭海城际铁路有限公司， 海宁 314499； 2.西南交通大学交通隧道工程教育部重点实验室， 成都 610031； 3.中铁第四勘察设计院集团有限公司， 武汉 430063)

近年来，中国隧道及地下工程开始大规模建设，各领域的隧道总数及长度均快速增加[1]。在这些隧道中，较多的工程使用了全断面机械化掘进设备，如盾构机、TBM设备。但是，随着中国盾构机应用领域及应用数量的不断扩大，机械设备施工面临的困难与日俱增[2]。同时，随着盾构机面临愈加复杂的地层环境，刀盘多同时采用滚刀和刮刀[3]。

在复合地层中，采用滚刀刀具的盾构机将面临偏磨的挑战[4]，偏磨的产生与隧道地质情况、滚刀响应机理有关[5-11]。目前室内实验能较好得分析滚刀对岩土体的破坏作用。邹飞等[12]基于实验分析了滚刀作用下岩体位移与应变场。张桂菊等[13]采用数值分析和实验相结合的形式对滚刀切削力学进行了分析。由于滚刀运动状态在实际掘进过程中的不可见性，研究者多采用数值模拟的手段进行研究。van Wyk等[14]以Paarl花岗岩为宏观对象，通过离散元数值分析得到了刀具的切削力和切削深度以及刀具磨损后的刀具切削面大小之间的联系。孙伟等[15]使用三维离散元方法，发现刀刃宽的增大将导致滚刀受力和碎岩体积的增加，而刀刃角对这两者的影响较小；张魁等[16]基于二维离散元数值方法，认为滚刀破岩过程中存在4种基本的破岩模式，这4种破岩模式和围岩以及刀间距有关。Cho等[17]使用三维有限元方法，将数值分析中滚刀的受力行为和室内滚刀线性切割试验进行了对比，最终利用三维有限元方法优化TBM(tunnel boring machine)滚刀间距。

目前以数值模拟为手段对滚刀的研究中，依然存在一定不足。前人研究中，针对硬岩破岩效率的研究较多，针对复合地层滚刀运动规律的较少；在三维数值分析中，多将滚刀转动速度设定为恒定值[14-20]，认为滚刀的角速度与线速度存在理想匹配，无法实现滚刀转动的动态平衡，与实际工况存在差异性。为此，使用刚体动力学(rigid body dynamics，RBD)-离散元(discrete element method,DEM)耦合数值模拟方法实现滚刀的动态滚动行为，并以此进行分析。

1 滚刀偏磨与滚刀响应机理分析

随着盾构机制造技术的不断进步，大直径盾构机的广泛应用，刀盘内外侧滚刀的运动速度存在显著差异，如图1所示。

滚刀的自转状态是由地层摩擦受力的综合体现，如图2所示。在图2中，滚刀迎土面受到的摩擦力f、破岩力F等，综合体现为滚刀轴承处的法向力Fn和切向力Ft，以及维持滚刀转动的扭矩Tt和阻碍滚刀转动的阻碍扭矩Tci，在大多数的理想岩体中，Tt一般大于Tci或维持动态平衡，滚刀能够形成和线速度相匹配的角速度，导致滚刀接触面和岩面的相对速度较低。

但是在复合地层中滚刀面对离散性较大的软质岩土体，滚刀转动扭矩Tt和阻碍扭矩Tci并不能维持较好的动态平衡状态。如滚刀启动时无法提供足够的Tt，将导致滚刀角速度和线速度动态匹配差异性较大，使得滚刀在岩土体的界面处存在较大的相对速度，如图3所示。

滚刀与岩土体界面处的相对速度大小将决定滚刀刀圈的磨损速率，同时进一步影响滚刀的运动状态即角速度和线速度的动态匹配效果。在单一匀质岩体中，滚刀刀圈接触面处的实际线速度和理论线速度较为接近，可以处于较为理想的工作状态；在复合地层中，滚刀不能达到较好的动态匹配效果，接触面处相对速度较高，从而导致滚刀处于非正常工作状态。现有工程研究也已表明，当滚刀无法保证正常的自转运动，将导致刀具的异常磨损[5-11]，如图4所示。

图1 不同安装位置下的滚刀线速度方向及大小Fig.1 The direction and size of the hob linear velocity under different installation positions

∝表示角速度与线速度的匹配图2 滚刀的运动状态Fig.2 Movement state of disc cutter

图3 复合地层中滚刀界面速度的差异性Fig.3 Difference of interface velocity in composite stratum

图4 复合地层中的滚刀偏磨Fig.4 Eccentric wear of disc cutter in composite stratum

2 RBD-DEM耦合数值模拟方法

常见DEM数值模拟中，几何体对象使用较为固定的运动设定，而实际工况中的几何体存在较为复杂的动力学运动RBD，如被动转动、被动移动等。因此，使用RBD-DEM耦合数值模拟的方法[21]，将DEM中刚体受到的复杂受力解传递至RBD中进行计算以决定DEM中下一时刻刚体的运动状态，最终实现滚刀转动被动平衡的数值模拟。

图5 刚体动力学计算时的物体坐标变化Fig.5 Object coordinate change in RBD numerical simulation

刚体动力学的计算运算时的物体坐标变化如图5所示。图5中，t0为初始计算时刻，此时物体的方向坐标系(local coordinate system，LCS)为初始坐标系，方向坐标系原点的全局坐标(global coordinate system，GCS)为物体初始化时质心位置。经过一段时间的运动后，物体沿一定运动轨迹运动至新的位置，此时t1时刻的方向坐标系原点的全局坐标、方向角发生一定改变，此时用向量和方向角方式表示为

(1)

式(1)中：dt为计算时间步，dt=t1-t0；L为运动向量；Lt为下一计算循环的运动向量；V为运动速度向量；O_AP为方向坐标系原点的全局坐标；APt为下一计算循环方向坐标系原点的全局坐标。

通过每一计算循环对运动向量的修改，完成物体的移动计算，物体的转动计算公式为

Ot+Δt=OtMt

(2)

Mt(Aaxis,θ)=

(3)

(4)

式中：Ot+Δt为下一计算循环步t+Δt时刻的物体方向坐标系的矩阵；Ot为当前计算循环时物体方向坐标系的矩阵描述；Mt为当前计算循环时的旋转矩阵；旋转矩阵可通过旋转轴Aaxis和改变角θ表示；ωt为当前循环的角速度向量，可表示为

(5)

式(5)中：α为角加速度；T为离散元计算得到的滚刀转矩矩阵；ωt+Δt为上一计算循环步t-Δt时刻的物体交速度向量；I为惯性矩阵。

最终的计算流程图如图6所示，其中RBD计算使用基于C++的EDEM二次开发实现，DEM则由EDEM本身进行。

3 数值模拟分析

3.1 数值模拟模型

图8中为100、150、200、300 kPa围压下的实验室与离散元三轴试验对比，σ1、σ3为大主应力与小主应力。

在计算中使用17寸标准的滚刀作为滚刀计算对象，在具体计算中将滚刀详细尺寸进行简化，尺寸如图9所示。

离散元模型如图10所示，模型颗粒总计256 131个，模型尺寸为1 m(长)×0.2 m(宽)×0.1 m(高)，计算时间步为1.274 24×10-5s。在计算中，模型除顶部无约束外，其余模型面均使用平面几何体进行模型的约束固定。墙体、滚刀与颗粒的接触采用Hertz-Mindlin理论，墙体弹性模拟取60 MPa、泊松比取0.3，滚刀弹性模量取206 GPa、泊松比为

图6 RBD-DEM耦合计算流程Fig.6 Calculation process of RBD-DEM coupling numerical simulation

表1 离散元计算参数

图7 离散元三轴试验模拟Fig.7 Simulation of triaxial test with DEM

σ1为大主应力；σ3为小主应力图8 离散元模拟三轴和室内试验数据对比Fig.8 Comparison of DEM and practical test

图9 17寸滚刀简化尺寸图Fig.9 Simplified dimension drawing of 17 inch disc cutter

图10 离散元计算模型Fig.10 computational model in DEM

0.25，墙体-颗粒的摩擦系数取0.5、滚动摩擦系数取0.5，滚刀-颗粒的摩擦系数取0.2、滚动摩擦系数取0.01。

滚刀从模型A侧向B侧线性移动，移动速度按工况设定，滚刀转动由RBD耦合计算程序控制。

3.2 不同贯入度下数值模拟结果分析

研究表明，在华南地区的复合地层盾构施工过程中，滚刀贯入度在1～16 mm，随着硬岩比例的增加而降低贯入度[22]。为探究不同贯入度下滚刀的动态运动方式，分别设置2、4、7、8、9、10、12、14 mm贯入度工况进行计算，得到的滚刀转速如图11所示。可以看出，随着贯入度的不断提高，滚刀响应速度逐渐边块，滚刀能更快地达到理论转速附近，处于较好的动态运动状态。但是从图11中可以发现，在低贯入度时(2、4 mm)，滚刀一直处于较低的转速，没有进入理想的动态平衡，滚刀和岩土体间将具有较大的界面相对速度，符合图3的猜想。但当滚刀的贯入度达到10 mm及以上时，能较好地达到动态平衡。计算最终各工况模拟得到的滚刀滑移率，如图12所示，其中滑移率S计算方式为

(6)

式(6)中：S为滑移率；ω′为实际转动角速度；ωo为理论转动角速度。

通过滑移率S可以表征滚刀在运动中滚动运动的占比，即界面相对速度的占比。滑移率S越高，说明界面相对运动速度的占比越高，当滑移率S达到1时说明物体处于完全滑动的运动状态。

图11 滚刀不同贯入度下动态响应的变化规律Fig.11 The changing law of disc cutter’s dynamic response with different penetration

图12 各贯入度下最终得到的滚刀滑移率Fig.12 Slip ratio in different simulated conditions

从图12中可以发现，随着贯入度的提高滑移率逐渐下降，当贯入度达到10 mm及以上时滑移率下降趋势不再明显，而是趋于缓和，此时滚刀的运动状态中滑动运动占比较少，可以认为滚刀处于动态平衡当中。但是如果继续加大贯入度，滚刀的滑移率下降不明显，说明贯入度对于滚刀进入合适动态状态的贡献降低。因此，可以认为贯入度对于滚刀进入合适动态平衡的贡献只存在于一定区间内，当小于这个区间时滚刀存在较大的界面滑移速度；当贯入度大于这个区间时并不能进一步提高动态响应效果。

3.3 不同速度下数值模拟结果分析

目前实际工程中，不同盾构机型号及地层差异导致盾构掘进时刀盘转速存在较大不同。研究表明，在复合地层中盾构机的刀盘转速区间为1.2～2.0 r/min[23-26]，因此使用1.5 r/min作为数值分析时的刀盘转速。考虑目前盾构机制造技术发展的日新月异、大直径盾构的普遍应用，以6 m作为模拟计算半径的上限，即在1.5 r/min时线速度约为1 m/s。结合计算机性能的约束，最终设定0.3、0.5、0.8、1.0 m/s为模拟工况，并采用上述计算得到的10 mm贯入度作为贯入度设定以保证滚刀出于合适的动态平衡。

由于滚刀均能进入合适的动态平衡状态，因此控制滚刀运动的主要因素为滚刀受到的扭矩，提取得到各工况中滚刀扭矩时程分布如图13所示。

图13 滚刀不同速度时的扭矩时程分布Fig.13 Torque time history distribution of disc cutter at different speeds

从图13中可以看出，由于使用离散元作为岩土对象的模拟手段，因此扭矩分布存在一定离散性，但是仍具有一定规律。滚刀受到的扭矩先随滚刀运动而升高，到达某一峰值点附近后开始下降，随后在零值附近波动。这规律与图11中进入动态平衡工况显示的角速度曲线相吻合，即滚刀的扭矩先增大后减小，在滚动阶段则在零值附近波动以保持转动的稳定性。图3显示，滚刀扭矩的提升速度和线性运动速度有关，随着线性运动速度的提高，滚刀达到扭矩峰值点附近的时间越短；同样的，线性运动速度越高滚刀的峰值扭矩越大。

图13中，各线性速度下的扭矩时程分布曲线规律与Cho等[17]进行岩石线性切割实验测得的滚动力分布类似，与其采用理想滚动的数值模拟结果有较大不同。因此认为在数值分析中，使用RBD耦合控制滚刀运动优化了数值分析边界条件，更为接近现实情况。对图13的扭矩的时程分布进行拟合，对比拟合函数峰值点和峰值距离的差异性如图14所示。

图14 各工况峰值扭矩和半峰宽值差异性对比Fig.14 Comparison of peak torque and half peak width under different working conditions

图14显示峰值扭矩和运动速度间基本成线性规律，对于线性速度越低的滚刀，其峰值扭矩更低，运动速度为1 m/s的滚刀峰值扭矩是运动速度为0.3 m/s时的1.86倍。半峰宽度的差异性显示，速度越快的滚刀达到峰值扭矩的时间更低，运动速度为0.3 m/s的滚刀达到峰值扭矩时间是运动速度为1 m/s时的1.72倍。

结合图13和图14的结果，可以认为滚刀线性运动速度对于滚刀启动快慢有一定影响。滚刀启动峰值扭矩和运动速度正相关，运动速度较低的滚刀峰值扭矩更低；滚刀达到峰值扭矩的时间和运动速度负相关，运动速度较低的滚刀达到峰值扭矩的时间更长。因此在复合地层中，近隧道轴线中心处的滚刀对安装精度、渣土和易性更为敏感，由于理想状态下启动峰值扭矩更低、响应时间更长，外界的少许扰动将使滚刀不正常滚动，最终造成异常磨损。

4 结论

使用RBD-DEM耦合数值模拟技术，实现了滚刀在数值模拟中的动态转动，进一步使针对滚刀的数值模拟符合工程实际。通过不同贯入度、不同运动速度，研究了滚刀进入动态转动的影响原因及响应快慢，得出如下结论。

(1)存在最优贯入度区间，低于该区间时滚刀不能进入动态转动状态，大于该区间时改善效果不明显且增大滚刀受力。

(2)滚刀启动时受到的扭矩先随滚刀运动而升高，到达某一峰值点附近后开始下降，随后在零值附近波动。

(3)滚刀启动峰值扭矩和运动速度正相关，运动速度较低的滚刀峰值扭矩更低。

(4)滚刀达到峰值扭矩的时间和运动速度负相关，运动速度较低的滚刀达到峰值扭矩的时间更长。

因此在复合地层中，滚刀存在一个最优贯入度区间使得滚刀能够出于优秀的动态响应状态且磨耗较低。同时近隧道轴线中心处的滚刀由于运动速度低，因此对安装精度、渣土和易性更为敏感，因为理想状态下启动峰值扭矩更低、响应时间更长，外界的少许扰动将使滚刀不正常滚动，最终造成异常磨损。通过数值模拟验证了RBD-DEM的合理性，为后续相关滚刀动态研究提供了参考；初步研究成果为盾构刀盘刀具设计及盾构机的操作提供了参考。