Ⅲ型应力强度因子分析的Williams单元

徐华¹，杨绿峰^1,2，佘振平¹

(1. 广西大学 土木建筑工程学院 工程防灾与结构安全教育部重点实验室，广西 南宁，530004；

2. 广西壮族自治区住房与城乡建设厅，广西 南宁，530028)

摘要：基于反平面断裂是工程结构的一类重要断裂类型，利用常规单元或奇异单元分析Ⅲ型裂尖应力强度因子(简记为SIF)时存在诸多困难，首先建立Ⅲ型裂纹裂尖奇异区单元总体位移场的Williams级数表达式，利用普通有限元形函数建立奇异域子单元的局部位移场，然后根据总体场控制局部场的原则，建立裂尖奇异域分析的Williams单元，以直接确定III型裂纹裂尖的SIF。在此基础上，研究Williams单元的径向离散因子、离散数和级数项这3个重要参数以及裂纹长度对裂尖SIF的影响，并给出这3个重要参数的建议值。研究结果表明：对于III型裂纹应力强度因子，Williams单元具有很高的计算精度和效率。

关键词：Ⅲ型裂纹；应力强度因子；Williams单元；广义参数有限元法

中图分类号：TU318          文献标志码：A         文章编号：1672-7207(2012)08-3237-07

Williams element for mode-III stress intensity factor

XU Hua¹, YANG Lu-feng^1,2, SHE Zhen-ping¹

(1. Key Laboratory of Disaster Prevention and Structural Safety of China, Ministry of Education,

School of Civil Engineering and Architecture, Guangxi University, Nanning 530004, China;

2. Department of Housing and Urban-Rural Development, Guangxi Zhuang Autonomous Region, Nanning 530028, China)

Abstract: Based on the fact that the anti-plane crack is one of the basic and crucial fracture type for engineering structures, and there are some limitations for traditional finite element in determination of mode-III stress intensity factor (SIF). Williams series for global displacement fields in singular region near the tip of mode-III crack were presented. The local displacement field of subelement were approximated by means of the shaped function of conventional finite element method.According to the principle that the global field determines the local one, the Williams element was presented for singular region near the crack tip, so that the mode-III SIF could be directly determined by the coefficients of the Williams element. Three important parameters, i.e the radial scale factor, the number of subelement and the series term, were discussed for Williams element. The influence of the relative length of crack on the SIF was investigated. The results show that the Williams element can achieve high accuracy and efficiency in evaluation of SIF.

Key words: mode-III crack; stress intensity factor; Williams element; method of finite element with generalized degrees of freedom

裂纹是结构构件中常见的缺陷形式，其基本类型有平面内的Ⅰ型(张开型)和Ⅱ型(滑开型)、平面外的Ⅲ型(撕开型)共3种。随着结构形式和受力特性的多样化，仅对平面内的断裂问题进行研究尚不能满足工程需求。目前，虽然国内外关于Ⅲ型断裂问题的研究主要集中在金属、复合材料和聚合物等领域，?但在地下工程、边坡工程、路面结构等建筑工程中，带裂隙岩体不仅承受Ⅰ型、Ⅱ型加载，而且往往会承受Ⅲ型加   载^[1]。在含纵向裂缝的道路结构中，车辆沿着连续的纵向裂缝行驶将引起剪切型移动，而在纵缝端部将导致边缘发生撕裂型移动^[2]。因此，开展含裂纹结构Ⅲ型断裂研究对工程结构的断裂分析和维护加固具有重要意义。在断裂力学的工程应用中，裂尖应力强度因子(简记为SIF)由于能够表征裂纹尖端应力场的强  度^[3-4]，成为结构失效判断、安全评估和寿命预测的重要力学参数。具有简单几何形状和加载方式的含裂纹体可以通过解析解或应力强度因子手册^[5]获得SIF，但对于一些复杂结构必须采用数值方法，其中有限元法因为具有很强的适应能力而得到广泛应用。由于裂纹尖端存在应力奇异性，应变场梯度大，常规有限元法需要采用很稠密的离散网格，导致计算速度慢、精度不高，解的收敛性也没保证。裂尖奇异单元^[6]尽管能够反映裂尖附近的应力奇异性，但推导复杂，且需要通过裂尖附近奇异域的SIF拟合并外推计算裂尖处的SIF，影响计算精度和计算效率。钱俊等^[7]应用分区混合有限元法，讨论了Ⅲ型复合材料切口应力应变的奇异性，并计算了切口的SIF。崔玉红等^[8]基于修正变分原理推导了混合Trefftz有限元法，计算反平面断裂问题的SIF。Williams^[9]提出了一种平面Ⅰ和Ⅱ型裂纹裂尖局部位移场和应力的级数表达式，能够很好地反映裂纹尖端应力场的奇异性，因而时常用于构造裂尖奇异区单元。Cheung等^[10]应用Williams级数建立了裂尖SIF分析的有限条模型；Leung等^[11]利用Williams级数建立了分形有限元法，开展断裂分析，计算裂尖SIF。Yang等^[12-13]研究建立了结构分析的广义参数有限元法，利用连续级数或光滑度较高的形函数建立结构分析的总体位移场，并通过总体场控制局部位移场，从而大大降低结构离散未知量，并应用于结构静力和动力分析^[14-15]。为此，本文作者推导Ⅲ型裂纹裂尖Williams级数解，分别利用Williams级数和普通有限元形函数建立裂尖奇异区单元的总体位移场和子单元的局部位移场，由总体场控制局部场的节点位移，建立广义参数子单元，进而利用有限项等比级数求和方法集合子单元刚度方程，形成裂尖奇异域Williams单元(简记为W单元)，以便高效、准确地计算Ⅲ型裂纹的裂尖SIF。

1  Ⅲ型裂纹裂尖奇异区位移场Williams级数解

以下分中心裂纹和边界裂纹2种情况，推导Ⅲ型裂纹裂尖奇异区位移场的Williams级数解。

1.1  中心裂纹裂尖奇异区位移场

图1所示为带中心裂纹板受面外剪切荷载示意图，含有中心裂纹的弹性薄板，其裂纹长度为2c，受面外剪切荷载作用处于反平面应变状态。以裂纹中心为平面坐标原点，建立Williams级数表示的裂尖奇异区位移场。

由于该弹性薄板沿x和y轴的位移满足u=0，v=0，沿z轴的出平面位移w≠0，因此，可得Laplace控制方程：

              (1)

式(1)为调和方程，其解可表示为：

         (2)

式中：z=x+iy为复变量；g(z)为复变解析函数^[16]，

       (3)

式中：A_n为待定复数，A_n=a_n+ib_n；a_n，b_n和λ_n均为实数。由此可得到应力分量的复变函数表达式：

，

     (4)

图1  带中心裂纹板受面外剪切荷载

Fig.1  Center cracked plate under anti-plane shear loading

将图1中的坐标原点从裂纹中心移到裂尖，则坐标变换：

式中：是极坐标，且有。此时，式(3)可表示为：

(5)

将式(5)代入式(2)和(4)，经整理得：

(6)

  (7)

  (8)

上述方程满足裂纹自由面边界条件：

        (9)

式中：上标“+”和“-”分别表示裂纹上表面和下表面。

引入边界条件，联立方程组解得：

     (10)

式中：c为裂纹半长；(r，θ)分别为裂尖奇异区任意点的极坐标；b_i为待定参数。式(10)为反平面III型中心裂纹断裂问题的Williams级数解，它满足裂纹尖端附近应力场的奇异性。

1.2  边界裂纹裂尖奇异区位移场

含有边界裂纹的弹性薄板，裂纹长度为c，受面外剪切荷载作用处于反平面应变状态，如图2所示。以裂尖为坐标原点建立坐标系，仿照中心裂纹裂尖的位移场推导过程，可求得边界裂纹的裂尖位移场：

         (11)

式中：b_i为待定参数；m为截取的Williams级数项数。

从式(10)和式(11)可以看出：与i=0相对应的项在裂纹尖端都具有奇异性。所以，系数b₀直接与SIF相关，也就是说，Ⅲ型裂纹的裂尖SIF可直接由系数b₀确定：

  (12)

图2  带边界裂纹板承受面外剪切荷载

Fig.2  Edge cracked plate under shear loading

2  裂尖奇异区广义参数Williams   单元

围绕裂尖O确定1个特定区域，它主要由应力奇异区构成(如图3所示)。并以裂尖O为中心用一组射线沿转角方向将奇异区离散为多个扇形条元，称为Williams单元(简记为W单元)。任取其中的1个W单元OAB作为典型单元，说明广义参数有限元法形成单元刚度矩阵的过程。利用1组环绕裂尖O且相互平行的折线(记为Γ₀，Γ₁，，Γ_n)将W单元离散为1组几何形状相似的梯形子单元和1个裂尖三角形OCD。任意2个相邻折线Γ_i-1和Γ_i(i=1，2，…，n)到裂尖的距离之比为常数α，称为径向离散比例因子。Γ_i-1和Γ_i两边的射线一起构成扇形条元内第i层子单元。

首先利用式(11)定义W单元OAB的总体位移场，并以矩阵形式表示：

             (13)

式中：

 (14)

其中：

         (15)

k_i为常数，对于中心裂纹和边界裂纹，分别等于(2c)^i+1/2和2。

图3  裂尖奇异域的离散

Fig.3  Discretization in singular region around crack tip

W单元OAB内包含n个子单元，每个子单元可根据常规有限元法建立8个节点的四边形单元位移场：

            (16)

式中：{w}^e为子单元上任一点位移列阵；为单元插值函数；{d}^e为单元节点位移列阵。根据广义参数有限元法思想，局部单元位移场必须与总体场保持位移协调，则子单元节点位移列阵{d}^e同单元总体位移场的待定参数之间满足如下关系：

              (17)

式中：[T]^e为单元转换矩阵。对于四边形，对于8节点四边形子单元，有：

 (18)

其中：(r_i，θ_i)为子单元节点i的极坐标。将式(17)代入(16)，可得：

       (19)

式中：；{b}中参数不局限于特定的物理意义，是广义参数。

根据式(19)，利用变分原理建立W单元OAB内第p层子单元的刚度方程：

          (20)

式中：下标W表示W单元；；；[K]^e和{f}^e分别表示常规8节点四边形等参元的刚度矩阵和荷载列阵。由于同一个W单元内各层子单元之间几何相似，所以，它们具有相同的[K]^e和{f}^e。

将式(14)代入式(18)，可得第p层子单元的[T]^(p)与第1层子单元的[T]⁽¹⁾之间存在如下关系：

[T]^(p)=[T]⁽¹⁾[diag(β_j)]；p=2，3，…，n     (21)

其中：[diag(β_j)]表示对角矩阵，其第j行对角元素为β_j，且有

β_j=α^(p-1)(j-1/2)；j=1，2，…，m+1      (22)

根据式(20)，将W单元内第2至第n层子单元的刚度方程进行叠加，可得：

[K′]_W{b}={f′}_W             (23)

式中：；。

由于W单元的最外层子单元ABFE部分位于奇异区，部分位于常规域，所以，称为过渡子单元，其单元节点编号如图4所示。因此，其节点位移列阵可以分块表示为：

(24)

式中：下标b代表位于W单元与常规域交界上的节点；下标s代表W单元的过渡子单元上位于奇异域的节点，且；

  (25)

图4  过渡子单元

Fig.4  Transitional subelement

根据式(24)，可以建立分块表示的过渡子单元刚度方程：

       (26)

式中：； ；，和均为常规等参元刚度矩阵分块表示时的子矩阵，且。

将式(23)与式(26)分块叠加，可形成W单元OAB的刚度方程：

 (27)

3  集成总体刚度方程

集成奇异域内所有W单元刚度方程，得：

      (28)

式中：{d}_b表示全体W单元与常规单元交界线上的节点位移参数列阵；和分别表示奇异区全部W单元集成的刚度和总体荷载列阵。

与W单元相邻的常规单元的刚度方程可以分块表示为：

       (29)

式中：下标r代表位于常规域的节点。

按照有限元法步骤，将全部W单元和常规单元刚度矩阵进行集成，得到分块表示的总体刚度方程：

   (30)

通过求解式(30)中待定参数{b}，并将其向量中b₀代入式(12)，即可求得应力强度因子K_III。

4  算例分析

   例1  含边界裂纹的半无限大弹性薄板如图5所示，裂纹长度c=16 cm，板宽l=160 cm，板高2h=600 cm，弹性模量E=2×10⁵ MPa，泊松比μ=0.167，沿z轴方向作用均布荷载τ=1 kN/cm。建立有限元计算模型时，要保证裂尖奇异区的离散网格与常规区的离散网格相协调，即W单元与相邻的常规单元在其分界线上共用相同的离散节点。为了减少有限元计算工作量，沿对称轴x轴取半板进行计算，此时，裂尖奇异区离散为4个扇形W单元即可保证计算结果的精度和收敛性。奇异区尺寸取0.5 cm，常规区离散为297个四边形8节点等参单元，整个计算模型共有984个节点，25个受约束节点位移。为验证W单元的计算精度，将计算结果与精确解^[5]

             (31)

相比较。同时，为了分析裂纹长度的影响，在固定板宽情况下，讨论裂纹相对长度c/l改变时的裂尖SIF，

其结果见表1。

图5  带边界裂纹的弹性板

Fig.5  An edge cracked elastic plate

表1  应力强度因子K_III与裂纹相对长度的关系

Table 1  Relationship between SIF and crack length

由表1可见：W单元的计算结果与精确解之间的误差很小，最大在0.63%以内，证明W单元具有很高的计算精度;当裂纹长度增大时，应力强度因子同步增大。

下面分析W单元中3个重要参数(径向离散因子α、离散数n以及级数项数m)对应力强度因子的影响。

K_III随着n和m变化如图6和图7所示。从图6和图7可见：随着奇异区径向离散比例因子α、径向离散数n逐渐增大，应力强度因子K_III的计算结果快速收敛，在α＞0.6和n≥100时K_III趋于稳定收敛。根据刚度矩阵形成过程可知：总体计算量不随α和n的增大而增加，所以，建议取α=0.9，n=300。从图8所示可以看出：随着Williams级数项数m的增加，计算结果逐渐收敛，当m≥4时K_III的结果趋于稳定。由于m增大时W单元的计算工作量随之增大，所以，建议取m=4。

例2  含中心裂纹的矩形弹性板如图9所示，裂纹长度2c=40 cm，板宽l=160 cm，板高2h=200 cm，弹性模量E=2×10⁵ MPa，泊松比μ=0.167。在       -c≤x≤c范围内沿z轴方向受均布荷载τ=1 kN/cm。沿着板的2个对称轴(x轴和y轴)取1/4板进行计算，其有限元离散模型的选取方法与原则与例1的相同。将裂尖奇异区离散为4个扇形W单元，奇异区长度取0.5 cm，常规区离散为28个四边形8节点等参单元，整个计算模型共有111个节点、11个受约束节点位移。为验证W单元的计算精度，将计算结果与解析解^[5]

             (32)

相比较。式中：F为矩形板修正系数^[5]，本算例中取F=0.565，采用W单元法和解析解所得应力强度因子K_III 分别为4.482 3 MPa·m^1/2和4.478 6 MPa·m^1/2，可见：本文方法计算结果与解析解较吻合，相对误差为0.83%，证明了W单元在不同裂纹位置、不同加载方式下所具有的良好适应性。这里研究W单元的3个重要参数α，n和m对计算结果的影响，见图6~8。从图6~8可以看出：对于中心裂纹，裂尖应力强度因子K_III的计算结果与α，n和m之间密切相关，并呈现出与例1相同的规律，验证了例1所建议的α，n和m取值在不同情况下的适应性。

图6  K_III与α的关系

Fig.6  Relationship between K_III and α

图7  K_III与n的关系

Fig.7  Relationship between K_III and n

图8  K_III与m的关系

Fig.8  Relationship between K_III and m

图9  带中心裂纹弹性板

Fig.9  A central cracked elastic plate

另外，作者对以上2个算例中的奇异区尺寸进行改变，通过分析发现：SIF计算结果对奇异区尺寸不太敏感；同时，在裂尖奇异区采用W单元代替普通单元，可以避免在该区域加密网格带来的离散未知量增大、计算效率降低的问题；而且W单元不像普通有限元那样需要通过插值拟合并外推计算确定裂尖SIF，而是根据单元的待定参数直接确定裂尖SIF，所以，计算简便，且计算精度高。

5  结论

(1) 建立了Ⅲ型断裂SIF分析的广义参数W单元，能直接计算应力强度因子，且具有较高的计算精度和计算效率；同时，奇异区尺寸的选择对计算结果无明显影响。

(2) W单元中的3个重要参数对SIF计算均有影响：径向离散比例因子α和离散数n增大时计算结果的精度随之提高，且不增加总的计算工作量；级数项数m增加时，计算精度和计算工作量同步增大。根据算例分析，建议取n=300，α=0.9，m=4。

参考文献：

[1] 谢海峰, 饶秋华, 王志. 反平面剪切(III型)加载下脆性岩石的断口分析[J]. 岩石力学与工程学报, 2007, 26(9): 1832-1839.
XIE Hai-feng, RAO Qiu-hua, WANG Zhi. Fracture morphology analysis of brittle rock under Anti-plane shear (Mode III) loading[J].Chinese Journal of Rock Mechanics and Engineering, 2007, 26(9): 1832-1839.

[2] 郑健龙, 周志刚, 张起森. 沥青路面抗裂设计理论与方法[M].北京: 人民交通出版社, 2003: 3-4.
ZHENG Jian-long, ZHOU Zhi-gang, ZHANG Qi-sen. Design theory and method for prevent crack of asphalt pavement[M]. BeiJing: China Communications Press , 2003: 3-4.

[3] Ghajar R, Moghaddam A S. Numerical investigation of the mode III stress intensity factors in FGMs considering the effect of graded Poisson’s ratio[J]. Engineering Fracture Mechanics, 2011, 78(7): 1478-1486.

[4] Cohen Y, Procaccia I. Stress intensity factor of mode-III cracks in thin sheets[J]. Physical Review E, 2011, 83(2): 026106(1-5).

[5] 中国航空研究院. 应力强度因子手册[M]. 增订版. 北京: 科学出版社, 1993: 494, 506-507.
China Aviation Institute. The stress intensity factor handbook[M]. Updated Version. Beijing: Science Press, 1993: 494, 506-507.

[6] Henshell R D, Shaw K G. Crack tip finite elements are unnecessary, International[J]. Journal for Numerical Methods in Engineering, 1975, 9(3): 495-507.

[7] 钱俊, 龙驭球. 分区混合有限元法计算反平面复合材料切口应力强度因子[J]. 计算结构力学及其应用, 1991, 8(3): 325-330.
QIAN Jun, LONG Yu-qiu. Sub region mixed FEM for calculating stress intensity factor of anti-plane notch in Bi-material[J]. Computational Structural Mechanics and Applications, 1991, 8(3): 325-330.

[8] 崔玉红, 秦庆华.混合Trefftz有限元法反平面断裂问题的探讨[J]. 固体力学学报, 2006, 27(2): 167-174.
CUI Yu-hong, QIN Qin-hua. Discussion on anti-plane crack problems by hybrid Trefftz finite element approach[J]. Acta Mechanica Solida Sinica, 2006, 27(2): 167-174.

[9] Williams M L. On the stress distribution at the base of a stationary crack[J]. Journal of Applied Mechanics, 1957, 24(1): 109-114.

[10] Cheung Y K, Jiang C P. Application of the finite strip method to plane fracture problems[J]. Engineering Fracture Mechanics, 1996, 53(1): 89-96.

[11] Leung A Y T, Tsang K L. Mode III two-Dimensional crack problem by two-level finite element method[J]. International Journal of Fracture, 2000, 102(3): 245-258.

[12] YANG Lu-feng, ZHAO Yan-Ling, LI Gui-qing. Finite elements with generalized coefficients for analysis of beams and plates[C]//Proceedings of 6th International Conference on Education and Practice of Computational Methods in Engineering and Science. Guangzhou, 1997: 560-565.

[13] LI Qiu-sheng, YANG Lu-feng, OU Xiao-duo, et al. The quintic finite element and finite strip with generalized degrees of freedom in structural analysis[J]. International Journal of Solids and Structures, 2001, 38(30/31): 5355-5372.

[14] LI Qiu-sheng, YANG Lu-feng, ZHAO Yan-Ling, et al. Dynamic analysis of non-uniform beams and plates by elements with generalized degrees of freedom[J]. International Journal of Mechanical Sciences, 2003, 45(5): 813-830.

[15] LI Qiu-sheng, YANG Lu-feng, LI Gui-qing. The quadratic finite element/strip with generalized degrees of freedom and their application[J]. Finite Elements in Analysis and Design, 2001, 37: 325-339.

[16] 张行, 陈宜歌. 反平面载荷作用下带界面裂纹的不同正交异性层板的应力分析[J]. 复合材料学报, 1989, 6(2): 1-8.
ZHANG Xing, CHEN Yi-ge. Stress intensity factor of dissimilar orthotropic laminate with an edge crack along the interface under anti-plane loading[J]. Acta Materiae Compositae Sinica, 1989, 6(2): 1-8.

(编辑 陈灿华)

收稿日期：2011-12-22；修回日期：2012-02-25

基金项目：国家自然科学基金资助项目 (51168003)；广西壮族自治区主席基金资助项目(2010GXNSFD169008)；工程防灾与结构安全教育部重点实验室主任基金资助项目(2009TMZR004)；广西大学科研基金资助项目(XGL090004)

通信作者：杨绿峰(1966-)，男，河南鲁山人，教授，博士生导师，从事结构工程的研究；电话：0771-3275070；E-mail：lfyang@gxu.edu.cn