基于新型流形法的三维应力强度因子求解

祁勇峰 苏海东 龚亚琦

摘 要：为保证在裂纹尖端具有最佳的逼近，裂纹尖端的解析解与其周边数值解联合应用的新型流形法可用来求解应力强度因子，但仅限于平面问题的Ⅰ型和Ⅱ型裂纹。沿用解析解与数值解联合应用的思路，以三维穿透直裂纹为研究对象，在裂纹尖端引入Williams解析解级数，应用解析覆盖与周边数值覆盖联合求解三维应力强度因子，相对于其他数值方法而言，计算精度高。推导相应的刚度矩阵和应变矩阵的表达式，通过典型算例验证了三维应力强度因子精确求解方法的有效性。

关键词：三维应力强度因子;数值流形方法;裂纹尖端Williams解析级数;解析覆盖;数值覆盖

Abstract： This paper proposes a new manifold method combined with the analytical solutions of the crack tip and its surrounding numerical solutions is used to solve the stress intensity factor to ensure the best approximation at the crack tip. But this method is limited to type I and type II cracks in the plane problem. This method follows the idea of combined application of analytical and numerical solutions. Taking the three-dimensional penetrating straight crack as the research object， the Williams analytical solution series is introduced at the crack tip， and the analytical covers and the surrounding numerical covers are combined to solve the 3D stress intensity factor. Compared with other numerical methods， this method has higher computational precision. Meanwhile， the corresponding formulas of the stiffness matrix and the strain matrix have been deduced，and typical numerical examples shows the validity of the method in solving the exact solution of the 3D stress intensity factor.

Keywords： 3D stress intensity factor; numerical manifold method; Williams series for crack tip; analytical covers; numerical covers

应力强度因子准则[1]是目前最常用的结构裂缝扩展准则之一，基于线弹性断裂力学的应力强度因子（Stress Intensity Factor，简称SIF）是用来表征裂纹尖端附近应力场和应变场强度，是控制裂纹尖端应力场和应变场强度的关键参数，在裂纹扩展分析中有着极其重要的地位。由于应力强度因子取决于外力的大小和分布、结构的几何条件以及裂缝的形状和位置，实际上只有少数问题存在解析解，对于复杂几何形状和加载条件的问题，只能通过数值方法来计算。目前裂缝扩展分析的主流数值计算方法有有限元法、扩展有限元法、数值流形法等。

有限元法[2-4]是目前分析裂缝等不连续问题的主要方法，为体现裂纹尖端（下简称裂尖）周围的应力集中和奇异性，往往需要在裂尖附近的复杂应力区布置高密度的单元网格，导致单元数目非常庞大;另外，在模拟裂缝扩展过程时，需不断重构网格，因此，有限元法对网格的要求和依赖性极大地影响了计算效率。扩展有限元法[5-7]通过在单元插值函数中引入能够反映裂缝面特性的不连续阶跃函数及反映裂尖局部特性的裂尖渐进位移场函数，裂缝可以穿过单元内部，裂缝扩展以后无需重新划分单元网格，采用同一网格就可以分析任意位置的裂縫问题，克服了常规有限元进行裂缝扩展分析的缺点，极大地简化了前处理工作。数值流形方法[8-11]引入不连续覆盖模拟裂缝，裂缝可以在网格内部穿过，巧妙地解决了常规有限元法裂缝面必须与单元边一致、裂缝扩展后需要重新划分网格的问题。相对于扩展有限元法设立不连续阶跃函数的方式而言，这种方式效果更好，在裂缝非常靠近单元边界时不会产生后者容易出现的数值误差。但无论是常规有限元法在裂尖布置细密网格的方式，还是扩展有限元法引入裂尖渐近位移场的方式，其主要目的都是为了提高裂尖附近的求解精度，从而提高应力强度因子的计算精度，这些方法都有改进的余地。即使是目前认为最适合于裂缝扩展分析的扩展有限元法，由于其只是使用了裂尖渐近位移场的部分特征函数，严格地说，还不能构成对裂尖位移场的最佳逼近。文献[12]提出的新型数值流形方法，实现了裂纹尖端的解析解与其周边数值解联合应用以求解应力强度因子，能够采用裂尖真实位移场的最佳逼近，并直接得到应力强度因子，计算精度高，但仅限于平面问题的Ⅰ型和Ⅱ型裂纹。

在上述研究的基础上，沿用解析解与数值解联合应用的思路，在裂纹尖端直接引入Williams解析解作为数学覆盖，应用裂纹尖端的解析解与周边数值解的三维流形覆盖联系技术，可直接得出裂纹尖端的三维应力强度因子，精度高且计算收敛快。

随着覆盖函数阶数的增加，三维应力强度因子计算值基本接近理论值。Williams级数的阶数以及周边数值覆盖阶数对计算值有较大影响。

1）Williams级数的阶数（m）影响最大。当m≤3，则KI与理论值相差较大;当m≥4时，KI接近理论值。

2）周边数值覆盖阶数升高有利于提高解的精度。当周边数值解取1阶时，KI的计算值普遍小于取2阶的情况，当m≥4时，KI与理论值接近，但仍有差距，仅当m取为7时才达到1.42，与理论值1.45最为接近。而当周边数值覆盖取2阶，m≥7时，计算值与理论值基本一致。

前期平面问题研究表明，在大的覆盖中单纯依靠提高覆盖函数阶次的方法往往会带来计算结果的振荡跳跃。反之，如果采用较小的覆盖而用相对简单的低阶多项式，则可以更好地逼近实际复杂的分布情况。表2的计算结果也说明，三维问题中，基于大覆盖，仅仅采用提高级数解及相邻覆盖函数的阶数的做法来提高计算精度，计算也表现出一定的不稳定性，要取得满意的计算精度，裂纹尖端及其周边的覆盖函数的阶数不小于7阶。

考虑到裂纹所在的独立覆盖区域较大，因此，将裂纹所在的独立覆盖进行局部加密，将裂纹尖端覆盖分别加密1倍及2倍，采用局部覆盖加密技术[16]，重新计算应力强度因子，如表2所示。

表2结果表明，当裂纹尖端独立覆盖加密1倍后、Williams级数的阶数≥4或当裂纹尖端独立覆盖加密2倍后，Williams级数的阶数≥3时，KI与理论解十分接近，且随着Williams级数阶数的提高，计算结果趋于稳定。

4.2 两端受剪切荷载作用

考虑三维荷载作用下的算例。如图5所示，无限长柱沿z方向施加均布荷载Q=100 kN/m2，裂纹长度a=0.02 m，柱体尺寸以及材料参数同上，其应力强度因子理论值为KIII=Qπc[15]，其中c=aw。KIII的理论值为2.51。

流形元网格如图6所示，共划分18个独立覆盖（方块区域）和25个部分重叠覆盖（窄条区域），每个独立覆盖的大小基本相同。裂纹所在的独立覆盖区域大小为0.2 m×0.2 m×1 m（长×宽×厚），柱体底面施加法向约束。应力强度因子的计算结果见表3。

表3计算结果表明：采用图6所示计算网格，当m≥7时，周边数值覆盖阶数取2、3的多项式阶数时，计算结果与理论值比较符合。当m≤7，计算结果与理论值有一定差别，局部数值出现跳跃，表明裂纹附近的网格还没有达到足够的密度。当网格加密一倍后，周边数值覆盖阶数均取2阶，当m≥7时，计算值与理论值十分接近，且随着阶数的提高，计算结果趋于稳定。

以上算例验证了三维裂缝计算公式和程序的正确性，表明裂纹尖端解析解覆盖和周边数值解覆盖联合应用求解三维线弹性断裂力学问题可行。与常规有限元方法相比，无需在裂纹尖端布置细密的网格，计算精度高，收敛相对较快。

裂纹尖端独立覆盖的密度、解析覆盖的级数以及相邻数值覆盖的阶数是影响应力强度因子计算精度的重要因素，但在保证独立覆盖有一定密度的情况下，提高与独立覆盖相邻数值覆盖的阶数可以得到应力强度因子的精确解。

裂纹尖端独立覆盖的合理布置对应力强度因子的计算精度及稳定性有一定的影响，因此，下一步要重点研究裂纹尖端附近的覆盖自动布置及密度问题，以保证方法的收敛性，便于开展三维裂缝扩展的动态模拟研究。

5 结论

将裂纹尖端解析解覆盖和周边数值解覆盖联合应用，分析三维线弹性断裂力学问题，得到以下主要结论：

1）在包含裂纹尖端的解析覆盖中，应用裂纹尖端附近的Williams位移解析解作为覆盖函数，并采用高阶多项式三维覆盖函数与解析覆盖的条形连接技术，实现了在解析覆盖中直接求得裂纹尖端的三维应力强度因子。

2）典型的张开型和撕开型的裂纹算例表明，应力强度因子的计算精度较高。鉴于三维裂缝扩展问题的复杂性，裂纹尖端周边数值覆盖阶数以及独立覆盖网格密度对应力强度因子计算精度的影响较二维问题更大。因此，协调独立覆盖密度、阶数与周边三维数值覆盖阶数的关系，来保证高精度求解收敛性的快速、稳定是下一步研究的重点。

考虑到解析级数是裂尖附近真实位移场的最佳逼近，相比其他方法而言，可以认为该方法在应力强度因子求解方面逼近效果更好、收敛更快，同时，由于网格布置根据不同区域的精度要求，只在裂尖附近进行覆盖加密，因此，相比采用均匀网格的扩展有限元而言，计算效率将有所提高，可以实现大规模计算。另外，应力强度因子SIF本身就是裂尖解析级数的未知数，在求解系统方程组时一并得到，而不需要像其它方法那样通过所谓的“直接”法或“间接”法来推求，不仅方便，而且不会引入额外误差，这也是该方法的优势所在。

该方法可以同时求解Ⅰ型、Ⅱ型、Ⅲ型（撕开型）裂纹的应力强度因子，应用复合型裂纹扩展准则就可以判斷其是否继续开裂，因此，该方法在三维裂缝扩展的动态模拟方面极具应用前景。

参考文献：

[1] PEREZ N. Linear-elastic fracture mechanics [M]∥ Fracture Mechanics. Cham： Springer International Publishing， 2016： 79-130.

[2] FAGEEHI Y A， ALSHOAIBI A M. Nonplanar crack growth simulation of multiple cracks using finite element method [J]. Advances in Materials Science and Engineering， 2020， 2020： 1-12.

[3] ALSHOAIBI M. Finite element simulation of crack growth path and stress intensity factors evaluation in linear elastic materials [J]. Jounal of Computational and Applied Research in Mechanical Engineering， 2019.

[4] HAJIAN M， MORADI M. An improved approach for computation of stress intensity factors using the finite element method [J]. Engineering Analysis with Boundary Elements， 2019（1）： 54-63.

[5] 李录贤， 王铁军. 扩展有限元法（XFEM）及其应用[J]. 力学进展， 2005， 35（1）： 5-20.

LI L X， WANG T J. The extended finite element method and its applications： A review [J]. Advances in Mechanics， 2005， 35（1）： 5-20. （in Chinese）

[6] GINER E， SUKUMAR N， DENIA F D， et al. Extended finite element method for fretting fatigue crack propagation [J]. International Journal of Solids and Structures， 2008， 45（22/23）： 5675-5687.

[7] ELGUEDJ T， DE SAINT MAURICE R P， COMBESCURE A， et al. Extended finite element modeling of 3D dynamic crack growth under impact loading [J]. Finite Elements in Analysis and Design， 2018， 151： 1-17.

[8] 石根华. 数值流形方法与非连续变形分析[M]. 裴觉民， 译. 北京： 清华大学出版社， 1997.

SHI G H.Numerical manifold method and discontinuous deformation analysis [M]. PEI J M， translated. Beijing： Tsinghua University Press， 1997.

[9] YANG S K， CAO M S， REN X H， et al. 3D crack propagation by the numerical manifold method [J]. Computers & Structures， 2018， 194： 116-129.

[10] ZHANG G X， LI X， LI H F. Simulation of hydraulic fracture utilizing numerical manifold method [J]. Science China Technological Sciences， 2015， 58（9）： 1542-1557.

[11] ZHENG H， LIU F， DU X L. Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method [J]. Computer Methods in Applied Mechanics and Engineering， 2015， 295： 150-171.

[12] 苏海东， 祁勇峰， 龚亚琦. 裂纹尖端解析解与周边数值解联合求解应力强度因子[J]. 长江科学院院报， 2013， 30（6）： 83-89.

SU H D， QI Y F， GONG Y Q. Compute stress intensity factors via combining analytical solutions around crack tips with surrounding numerical solutions [J]. Journal of Yangtze River Scientific Research Institute， 2013， 30（6）： 83-89.（in Chinese）

[13] 祁勇峰， 苏海东， 崔建华. 部分重叠覆盖的数值流形方法初步研究[J]. 长江科学院院报， 2013， 30（1）： 65-70.

QI Y F， SU H D， CUI J H. Preliminary study on numerical manifold method with partially overlapping covers [J]. Journal of Yangtze River Scientific Research Institute， 2013， 30（1）： 65-70. （in Chinese）

[14] 王勖成. 有限单元法[M]. 北京： 清华大学出版社， 2003.

WANG X C. Finite element method [M]. Beijing： Tsinghua University Press， 2003.（in Chinese）

[15] 中国航空研究院. 应力强度因子手册 [M]. 北京： 科学出版社， 1993.

Chinese Aeronautical Establishment. Stress intensity factors handbook [M]. Beijing： Science Press， 1993.（in Chinese）

[16] 苏海东， 龚亚琦， 颉志强， 等. 基于矩形独立覆盖初步实现结构静力分析的自动计算[J]. 长江科学院院报， 2016， 33（2）： 144-150.

SU H D， GONG Y Q， XIE Z Q， et al. Preliminary implementation of automatic computation for static analysis of structures using NMM based on independent rectangular covers [J]. Journal of Yangtze River Scientific Research Institute， 2016， 33（2）： 144-150. （in Chinese）

（編辑 章润红）