全断面隧道掘进机刀盘裂纹尖端应力强度因子收敛性分析

凌静秀 孙 伟 杨晓静 童 昕

1.福建工程学院机械与汽车工程学院,福州，350118 2.大连理工大学机械工程学院，大连，116024

全断面隧道掘进机刀盘裂纹尖端应力强度因子收敛性分析

凌静秀1孙 伟2杨晓静1童 昕1

1.福建工程学院机械与汽车工程学院,福州，350118 2.大连理工大学机械工程学院，大连，116024

针对全断面隧道掘进机刀盘裂纹损伤及寿命预测等工程问题，提出了基于子模型技术的应力强度因子求解方法，并用含裂纹的矩形钢板对该方法进行了验证，分析了裂纹网格参数对刀盘裂纹尖端应力强度因子的收敛性影响。结果表明，钢板的应力强度因子数值和理论计算结果最大相对误差为3.6%。同时得到了保证刀盘应力强度因子求解精度和效率的裂纹单元网格参数，为结构的裂纹扩展寿命预测提供了参考。

全断面隧道掘进机刀盘；裂纹；应力强度因子；子模型技术；有限元法

0 引言

全断面隧道掘进机(tunnel boring machine，TBM) 是利用回转刀具开挖及破碎洞内围岩来开凿岩石隧道，形成完整隧道断面的一种大型高端隧道工程装备。刀盘作为TBM关键部件，与围岩直接接触，是支撑滚刀切削岩石的箱体结构。由于地质环境及掘进参数的多变性，整机在掘进过程中振动极其剧烈，导致刀盘盘体开裂、支撑结构大变形及主轴承失效等问题[1-2]。刀盘结构在加工、制造及装配过程中，不可避免地存在初始缺陷。这些初始缺陷会演化为疲劳破坏的裂纹源，在破岩载荷作用下裂纹不断扩展，直至临界尺寸，导致结构破坏，失去服役能力，因此，确定初始裂纹扩展到失效尺寸的时间，即预测刀盘的裂纹扩展寿命、合理评估含表面裂纹刀盘的抗断裂能力显得极为重要。

在进行结构断裂分析时，应力强度因子是其中一个重要指标，是判断裂纹扩展规律和驱动裂纹扩展的关键参数。应力强度因子求解方法可归结为理论计算和试验两类方法，而理论方法包括解析法和数值法。解析法只能计算相对简单的模型；数值法包含有限元法(FEM)、有限差分法(FDM)及边界配置法(BGNP)等。解析法受到很多限制，数值法则在断裂力学中得到广泛的应用[3]。在用有限元法计算应力强度因子时，裂纹网格参数会对结果有一定的影响，需要合理选择网格参数，进而得到稳定、准确的数值解。

在刀盘载荷及系统设计方面，相关学者已开展了大量的研究工作，且取得了一定的研究成果。谢启江等[4]通过建立刀盘、机械系统和撑靴接触刚度耦合力传递模型，研究了刀盘载荷波动与岩石接触界面刚度的耦合关系。韩美东等[5]、夏毅敏等[6]采用ABAQUS 有限元软件模拟刀盘掘进破碎岩石的过程，分析了刀盘载荷变化规律及不同参数对其统计值的影响，并用工程实际值进行了验证。霍军周等[7]建立了带复杂性能约束的刀具布置优化模型，并用工程案例进行了验证。另外，文献[8-10]采用集中参数质量法、多体动力学仿真法及现场实测等手段对TBM刀盘系统的振动特性展开研究，分析了掘进参数、结构参数等不同参数对刀盘振动的影响，并与实测振动数据进行对比，来验证理论动力学模型的有效性。

关于刀盘结构损伤的研究国内外基本处于空白，笔者基于系统动力学及有限元法分析了不同参数对刀盘裂纹应力强度因子的影响，提出了适用长厚板在复杂应力状态下的新型失效判据准则，并预测了刀盘的裂纹扩展寿命[2,11]，但没有分析裂纹网格参数对刀盘应力强度因子的影响，仅采用默认参数进行研究。本文以辽西北引水工程TBM刀盘为研究对象，采用奇异单元对刀盘分体结构裂纹尖端附近区域进行精细网格划分，应用基于子模型技术的有限元方法直接计算刀盘裂纹尖端的应力强度因子，分析裂纹网格参数对应力强度因子收敛性的影响，进而确定较为稳定可靠的裂纹奇异单元网格参数。

1 应力强度因子求解方法及有效性验证

1.1 应力强度因子求解方法

数值方法能够求解复杂结构的三维应力强度因子，而有限元法中的奇异单元法是数值方法中求解应力强度因子适用面最广的一种方法。构件中任意三维裂纹及尖端坐标系可抽象成图1所示坐标系。

图1 裂纹尖端局部坐标系Fig.1 Local coordinate system of crack tip

根据线弹性断裂力学的解析解，由任意外载作用所产生的裂纹尖端附近区域的位移场和应力强度因子的关系可表示如下[12]：

(1)

其中，r、θ为局部柱坐标系中的两个坐标分量；u、v、w分别为裂纹尖端任一点的径向位移、法向位移和切向位移；KⅠ、KⅡ、KⅢ分别为Ⅰ型、Ⅱ型及Ⅲ型应力强度因子；G为剪切模量；k为与材料泊松比μ有关的常数，对于平面应变问题，k=3-4μ，而对于平面应力问题，k=(3-λμ)/(1+μ)。

对于平面应变状态，Ⅰ型裂纹的应力强度因子KⅠ可用裂纹面的法向位移表示，即

(2)

式中，E为弹性模量。

有限元软件ANSYS中提供了一种求解三维裂纹问题的奇异单元，如图2所示。裂纹尖端大部分区域处于平面应变状态，在求得应力应变解后，取裂纹尖端1/4节点的裂纹张开位移v(1/4)代替尖端位移，代入式(2)中即可求得应力强度因子[12]：

(3)

式中，r(1/4)为1/4节点到裂纹尖端的距离。

图2 ANSYS提供的三维奇异单元Fig.2 Three dimensional singular element provided by ANSYS

1.2 本方法的有效性验证

采用上述有限元子模型及计算方法，可方便快捷地求解出刀盘分体表面裂纹的应力强度因子，但计算结果是否正确，需要进行验证和分析。然而，对于这类复杂结构的三维表面裂纹问题，目前没有切实可行的验证方法，因此，本文以中心部位含有半椭圆型三维表面裂纹的矩形钢板为例，两端承受均匀拉应力σt=100 MPa，如图3所示。图中，钢板长度L=45 mm，宽W=70 mm，厚T=8 mm；裂纹长半轴c=4 mm，短半轴a=2 mm，θ为裂纹离心角。用ANSYS软件求解该类裂纹问题的应力强度因子，同时采用工程界广泛认可的Newman-Raju经验公式计算[13]，二者对比，以检验本文奇异单元法的可靠性。

(a)含表面裂纹的钢板

(b)裂纹离心角图3 含表面裂纹钢板的力学模型Fig.3 Mechanical model of steel plate with a surface crack

采用四面体单元划分钢板整体网格，裂纹区域引入奇异单元，模拟裂纹尖端应力场的奇异性，设置合适的裂纹网格参数，建立有限元模型，如图4所示。在钢板一端施加拉应力，一端固定，等效两端受拉情况，载荷及边界条件见图3。裂纹尖端椭圆曲线共划分为60份，含60组单元和61个节点，每个节点对应一个离心角，可得到61个应力强度因子。裂纹尖端径向划分为6层单元，即裂纹尖端各节点可得到6个应力强度因子。一般来说，越靠近内层的计算结果波动越大，越靠近外层的值则趋于稳定。为保证求解精度，程序默认每个裂纹尖端节点从最内层开始，会依次输出6个应力强度因子值，本文选取最外层的值作为数值计算结果，提取结果如图5所示。

(a)网格划分(b)载荷与边界条件图4 有限元分析模型Fig.4 Finite element analysis model

图5 ANSYS应力强度因子计算结果Fig.5 Stress intensity factor results calculated by ANSYS

由计算结果可知，应力强度因子基本呈现左右对称的趋势，即裂纹最深处应力强度因子最大，两个表面点的值相对较小，这和实际情况相符。因此，采用Newman-Raju公式计算时，仅计算裂纹离心角θ在0～90°的应力强度因子值K1，并与数值计算结果进行比较，结果如图6所示。

图6 有限元与Newman-Raju公式计算结果对比Fig.6 Comparison of calculation results between FEM and Newman-Raju formula

由以上对比结果可知，总体而言，有限元结果与Newman-Raju公式的计算结果基本相符，最大相对误差为3.6%，在工程允许的计算误差范围内，结果较符合实际。有限元计算结果与经验公式结果在裂纹表面点处基本重合，随着离心角的增大，有限元结果逐渐大于经验公式结果，然后二者又慢慢接近，靠近裂纹尖端最深处区域时，经验公式结果又稍微大些，这些计算结果的误差变化很小，基本可以忽略。这说明按照本文的有限元方法，应用ANSYS软件对刀盘分体三维表面裂纹的应力强度因子进行分析,理论上是可行的。

2 刀盘分体有限元模型

基于上述有限元子模型技术，合理等效刀盘分体结构，删除滚刀、焊缝等细小结构，分割出裂纹子模型结构，利用ANSYS软件对其进行网格划分，采用四面体单元划分裂纹子模型，奇异单元表征裂纹尖端的奇异性，分体其余实体结构采用六面体单元划分。同时经过严格控制整体网格大小和裂纹尖端的网格参数，精细划分刀盘分体结构，提高求解效率和精度。建立含半椭圆型三维表面裂纹的刀盘分体有限元模型，如图7所示，同时加载刀盘载荷及位移边界约束[2]。

图7 含裂纹特征的刀盘分体有限元模型[2]Fig.7 FEM model of cutterhead piece with crack[2]

3 裂纹网格参数对应力强度因子的影响

在上述有限元模型的基础上，即可求出不同载荷作用下的应力强度因子值，但裂纹尖端附近的网格参数变化会对结果有所影响，已有研究针对不同问题提出了相应的裂纹网格参数选择要求[13-14]。由此，需要结合刀盘有限元模型及实际载荷边界，通过收敛性检验确定出合适的裂纹网格参数。ANSYS中裂纹尖端区域的网格参数如图8所示，本文主要分析裂纹尖端周向网格份数nC、裂纹尖端径向网格份数nM、最大轮廓半径Ltip及裂纹尖端曲线网格份数nCF这四个裂纹网格控制参数的变化对刀盘裂纹尖端应力强度因子的影响。

图8 裂纹尖端区域的网格参数Fig.8 Mesh parameters in the crack tip region

3.1 周向网格份数

在图7所示的有限元模型中，取裂纹短半轴a=30 mm，长半轴c=150 mm。在裂纹网格参数nM=6，Ltip=10,nCF=100的情况下，分析裂纹尖端周向网格份数分别为8、16、24、40时的应力强度因子分布，图9所示为不同nC值对应的裂纹尖端区域的横截面，计算得到不同nC值的应力强度因子值如图10所示。

(a)nC=8 (b)nC=16

(c)nC=24 (d)nC=40图9 不同nC值的裂纹尖端区域横截面Fig.9 Cross sections of crack tip region with different nC values

图10 不同nC值时的应力强度因子Fig.10 Stress intensity factors with different nC values

由图10可知，裂纹尖端周向网格份数的变化对应力强度因子几乎没有影响，随着nC值的变化，计算结果最大幅度变化相对误差在1%以内。为了提高求解效率，后续建立刀盘裂纹有限元模型时采用nC=8即可满足精度要求。

3.2 裂纹尖端径向网格份数

保持上述有限元模型及边界条件不变，取nC=8，Ltip=10，nCF=100，分析裂纹尖端径向网格数nM分别为6、8、10、20时的应力强度因子变化，不同nM值对应的裂纹尖端区域的横截面如图11所示。在第1谱级载荷幅值作用下，不同nM值的有限元模型计算得到的应力强度因子值如图12所示。

(a)nM=6 (b)nM=8

(c)nM=10 (d)nM=20图11 不同nM值对应的裂纹尖端区域的横截面Fig.11 Cross sections of crack tip region with different nM values

图12 不同nM值对应的应力强度因子Fig.12 Stress intensity factors with different nM values

由图12可知，裂纹尖端径向网格份数的变化对应力强度因子的影响基本可以忽略，随着nM值的变化，应力强度因子最大值与最小值之间相对误差仅1.4%，在后续有限元建模时nM可取6～20之间的任意整数，本文取nM=8。

3.3 最大轮廓半径

在ANSYS软件中，裂纹尖端区域大小与裂纹的最大轮廓半径成一定比例关系，程序会根据人为设定的最大轮廓半径Ltip自动控制裂纹尖端区域的大小。同样在以上模型的基础上，取nC=8，nM=8，nCF=100，分析最大轮廓半径Ltip分别为2、6、10、14、18时的应力强度因子变化，结果如图13所示。

图13 不同Ltip值对应的应力强度因子Fig.13 Stress intensity factors with different Ltip values

由图13可知，裂纹的最大轮廓半径变化对应力强度因子有一定的影响，而当Ltip值大于6时，计算结果几乎保持不变，说明随着Ltip的增大，应力强度因子逐渐趋于稳定。综合考虑计算效率和求解精度，在后续有限元求解时可取Ltip=10。

3.4 曲线网格份数

同样保持上述有限元模型不变，取nC=8，nM=8，Ltip=10，分析裂纹尖端曲线网格份数nCF分别为40、60、70、100时的应力强度因子变化，得到结果如图14所示。

图14 不同nCF值对应的应力强度因子Fig.14 Stress intensity factors with different nCF values

由图14可知，裂纹尖端曲线网格份数的变化对应力强度因子影响很小，随着nCF值的增大，应力强度因子最大值与最小值之间相对误差仅1.2%，在后续有限元建模时nCF可取40～100之间的任意整数，本文取nCF=40即可满足精度要求。

通过分析对比有限元模型中裂纹网格参数的变化对应力强度因子的影响可知，除裂纹的最大轮廓半径外，其余裂纹网格参数对计算结果几乎没有影响。为保证后续分析结果的稳定可靠及提高计算效率，最终确定适合本文刀盘裂纹模型的网格参数如下：nC=8，nM=8，Ltip=10，nCF=40。

4 结论

(1)针对TBM刀盘的载荷特殊性及结构的复杂性，提出基于有限元子模型技术的裂纹应力强度因子求解方法，并用含裂纹的矩形钢板进行验证，得到数值和理论计算结果最大相对误差为3.6%，在工程允许的误差范围内。

(2) 将提出的应力强度因子求解方法应用于实际引水工程的TBM刀盘结构裂纹，同时分析不同的裂纹网格参数对应力强度因子的影响，除裂纹的最大轮廓半径外，其余参数对计算结果影响不大。

(3)通过分析得到了适用于求解实例刀盘应力强度因子的裂纹单元网格参数，为刀盘的裂纹扩展寿命预测提供了技术支撑。

[1] 黄文鹏.节理密集带地质硬岩TBM刀盘损坏形式及对策[J].隧道建设,2012,32(4):587-593.HANGWenpeng.StudyonDamageofHardRockTBMInducedbyTunnelingthroughJoint-densely-developedHardRockSectionandCountermeasures[J].TunnelConstruction, 2012, 32(4):587-593.

[2]LINGJingxiu,SUNWei,HUOJunzhou,etal.StudyofTBMCutterheadFatigueCrackPropagationLifeBasedonMulti-degreeofFreedomCouplingSystemDynamics[J].Computers&IndustrialEngineering, 2015, 83:1-14.

[3] 衣振华.疲劳裂纹扩展研究及在装载机横梁寿命估算中的应用[D].济南:山东大学，2011.YIZhenhua.ResearchonFatigueCrackPropagationandItsApplicationsinLoaderBeamLife[D].Jinan:ShandongUniversity, 2011.

[4] 谢启江, 余海东. 硬岩掘进机刀盘载荷与撑靴接触界面刚度的耦合关系[J]. 上海交通大学学报, 2015, 49(9):1269-1275.XIEQijiang,YUHaidong.CouplingRelationshipbetweenLoadsonCutterheadofTunnelBoringMachineandContactStiffnessofGripperShoesandRocks[J].JournalofShanghaiJiaotongUniversity, 2015, 49(9):1269-1275.

[5] 韩美东, 曲传咏, 蔡宗熙,等. 刀盘掘进过程动态仿真[J]. 哈尔滨工程大学学报, 2015(8):1098-1102.HANMeidong，QUChuanyong，CAIZongxi,etal.DynamicNumericalSimulationofTunnelingbytheTBMCutterHead[J].JournalofHarbinEngineeringUniversity, 2015(8):1098-1102.

[6] 夏毅敏, 朱湘衡, 林赉贶,等.TBM刀盘掘进载荷影响因素研究[J]. 现代制造工程, 2015(9):1-6.XIAYimin，ZHUXiangheng，LINLaikuang,etal.AnalysisofInfluenceFactorsonLoadofTBMCutterhead[J].ModernManufacturingEngineering, 2015(9):1-6.

[7] 霍军周, 史彦军, 滕弘飞,等. 全断面岩石掘进机刀具布置设计方法[J]. 中国机械工程, 2008, 19(15):1832-1836.HUOJunzhou,SHIYanjun,TENGHongfei,etal.CutterLayoutDesignofFull-faceRockTunnelBoringMachine(TBM) [J].ChinaMechanicalEngineering, 2008, 19(15):1832-1836.

[8]SUNWei,LINGJingxiu,HUOJunzhou,etal.DynamicCharacteristicsStudywithMultidegree-of-freedomCouplinginTBMCutterheadSystemBasedonComplexFactors[J].MathematicalProblemsinEngineering, 2013(3):657-675.

[9]HUOJunzhou,WUHangyang,YANGJing,etal.Multi-directionalCouplingDynamicCharacteristicsAnalysisofTBMCutterheadSystemBasedonTunnellingFieldTest[J].JournalofMechanicalScience&Technology, 2015, 29(8):3043-3058.

[10] 霍军周, 欧阳湘宇, 王亚杰,等. 重载冲击激励下TBM刀盘振动特性的影响因素分析[J]. 哈尔滨工程大学学报, 2015(4):555-559.HUOJunzhou,OUYANGXiangyu,WANGYajie,etal.AnalysisofInfluencingFactorsofVibrationBehaviorsofTBMCutterheadunderHeavyImpactLoads[J].JournalofHarbinEngineeringUniversity, 2015(4):555-559.

[11]SUNWei,LINGJing,HUOJunzhou,etal.StudyofTBMCutterheadFatigueDamageMechanismsBasedonaSegmentedComprehensiveFailureCriterion[J].EngineeringFailureAnalysis, 2015, 58:64-82.

[12]LINXB,SMITHRA.FiniteElementModelingofFatigueCrackGrowthofSurfaceCrackedPlates—PartI:TheNumericalTechnique[J].EngineeringFractureMechanics, 1999, 63(5):503-522.

[13]NAMIMR.StressIntensityFactorsinaRotatingImpellerContainingSemi-ellipticalSurfaceCrack[J].MechanicsBasedDesignofStructuresandMachines, 2012, 40(1):1-18.

[14]CHIEWSP,LIEST,LEECK,etal.StressIntensityFactorsforaSurfaceCrackinaTubularT-joint[J].InternationalJournalofPressureVesselsandPiping, 2001, 78(10):677-685.

(编辑 王旻玥)

Convergence Analysis of Stress Intensity Factor at Crack Tips for TBM Cutterheads

LING Jingxiu1SUN Wei2YANG Xiaojing1TONG Xin1

1.School of Mechanical and Automotive Engineering,Fujian University of Technology,Fuzhou,350118 2.School of Mechanical Engineering,Dalian University of Technology,Dalian,Liaoning,116024

Aiming at the engineering problems such as crack damage and life prediction of TBM cutterheads, a solution method of stress intensity factors was proposed based on the sub-model technique, and it was validated by a rectangular steel plate with a crack. Then the convergence effects of crack mesh parameters on stress intensity factors at the crack tips of cutterheads were analyzed. The analysis results show that the maximum relative errors of numerical and theoretical results are as 3.6%. Besides, the crack element mesh parameters were obtained that ensured both of the accuracy and efficiency for cutterhead stress intensity factors solutions, which may provide a reference for the crack growth life prediction of structures.

tunnel boring machine(TBM) cutterhead; crack; stress intensity factor; sub-model technique; finite element method(FEM)

2016-06-28

国家重点基础研究发展计划(973计划)资助项目(2013CB035402);福建省科技创新平台建设项目(2014H2002)；福建省自然科学基金资助项目(2016J01722,2017J01675)；福建工程学院科研启动基金资助项目(GY-Z160048)

TP391.9

10.3969/j.issn.1004-132X.2017.10.002

凌静秀，男，1985年生。福建工程学院机械与汽车工程学院讲师、博士。主要研究方向为复杂机械装备动力学分析及疲劳寿命预测。E-mail：ljxyxj@fjut.edu.cn。孙 伟，男，1967年生。大连理工大学机械工程学院教授、博士研究生导师。杨晓静，女，1985年生。福建工程学院机械与汽车工程学院助理实验师。童 昕，男，1964年生。福建工程学院机械与汽车工程学院教授、博士研究生导师。