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

祁勇峰，苏海东，龚亚琦

(长江科学院 材料与结构研究所；水利部水工程安全和病害防治工程技术研究中心，武汉 430010)

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

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

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

1 裂纹尖端位移场的Williams级数解

如图1所示的裂纹尖端区域，Williams位移级数解

图1 裂纹尖端坐标系

(1)

(2)

2 裂纹尖端位移场的解析解覆盖

在六面体网格内，将式(1)写成矩阵形式。

对于i=0:

(3)

对于i≥1:

(4)

将二维数值覆盖强制约束方法[13]推至三维，形成独立的三维解析覆盖，权函数wj均为结点1的权函数w1。其中，m为覆盖函数项数，下同。

对于第i项(i≥1)，其应变子矩阵为

(5)

(6)

将各项偏导数展开

(7)

式中：

另外，

刚度子矩阵

(8)

式中：D为弹性矩阵；V为流形元体积。

3 裂纹周边网格的解析解与数值解覆盖

在裂纹周边的各网格内，位移统一表示为式(9)，其中i≥1。

(9)

式中：dnkl、enkl、fnkl为待求的多项式系数；n、k、l为多项式阶次(p为多项式的最高阶次)；wj1、wj2为相应于各结点的权函数；J和L的个数根据网格而定，但保持J+L=8。显然，取L=0则退化到解析解覆盖式(4)。

如图2所示(图中窄条为部分重叠区域，长方体区域为独立覆盖，阴影为裂纹)，采用强制约束的方法，将网格结点28、40、39、33、34、45、46约束到结点27，即令

图2 六面体数学网格中的解析解覆盖和数值解覆盖

对于周边的其他网格，比如：网格1-2-14-13-7-8-20-19中的结点均为数值解的独立覆盖；解析覆盖的相邻网格39-40-52-51-45-46-58-57中的结点39、40、45、46为解析结点，其他4个结点为数值结点；而裂纹自由端所在网格25-26-38-37-31-32-44-43在竖直方向采用了两个独立覆盖，实现裂纹两边的独立运动。

对于解析覆盖结点j处的第i项，其应变矩阵为

(10)

(11)

令f(x,y,z)=xn-k-lylzk，则式(11)可写为

(12)

解析覆盖与数值覆盖相关的刚度子矩阵为

(13)

由式(10)和式(12)可见，刚度矩阵积分中具有非多项式的函数，基于多形式函数的单纯形积分公式无法应用，必须采用数值积分。因此，采用四面体区域的Hammer积分[14]，将流形元边界上的三角片与其形心相连形成四面体，然后再进行积分计算。

4 算例分析

以含边界裂纹的无限长柱体为例，考虑两种不同类型荷载。

4.1 两端受均布拉力

图3 两端受均布拉力的无限长柱体内的边界裂缝

整体的流形元网格如图4所示，共划分6个独立覆盖和9个部分重叠覆盖(窄条区域)，每个独立覆盖的大小基本相同。裂纹所在的独立覆盖区域大小为0.65 m×0.2 m×1 m(长×宽×厚)，柱体底面施加法向约束。

图4 部分重叠覆盖的流形元网格

应力强度因子的计算结果见表1，考虑了裂纹尖端所在的独立覆盖取Williams级数的不同项数，独立覆盖周边的数值解覆盖可取多项式的不同阶数。

表1 应力强度因子

随着覆盖函数阶数的增加，三维应力强度因子计算值基本接近理论值。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 应力强度因子(n=2)

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

4.2 两端受剪切荷载作用

图5 两端受剪力的无限长柱体内的边界裂缝

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

图6 部分重叠覆盖的流形元网格

表3 应力强度因子

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

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

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

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

5 结论

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

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

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

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

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