Bui Manh Tuan 陈云飞

(1东南大学机械工程学院，南京210081)

(2濉和工业学院机械工程学院，越南濉和56000)

在进行断裂安全分析时，应力强度因子(stress intensity factor，SIF)K是判断裂纹体在载荷作用下是否产生破坏的主要参数.目前，确定应力强度因子的方法较多，典型的有解析法、有限元法(FEM)、虚拟裂纹扩展技术［1］(virtual crack extension techniques)、修改裂纹闭合积分技术［2］(modified crack closure integral techniques)、J-集成技术［3-4］(J-integration techniques)、节点位移技术［5］等.由于有限元法计算结果相对精确，因此对于结构或裂纹形状复杂和受复杂载荷作用的结构件，一般采用有限单元法.常用的有限元软件包括ANSYS，ABAQUS，NASTRAN，MSC，Hyperworks等.与其他方法相比，有限元法具有单元的布局灵活、节点的配置方式比较任意，对裂纹形状、位置无特殊限制等优点.因此，对几何形状和载荷比较复杂的裂纹体，有限元法是确定应力强度因子最为有效的方法.本文基于有限元法和节点位移方法，研究了在中间裂纹与侧裂纹处裂纹尖端的应力强度因子和应力场.使用8节点四边形等参元和1/4节点退化单元，对网格密度和裂纹长度对计算精度的影响进行了研究.

1 裂纹尖端附近的应力和位移

在断裂力学中，按裂纹体受力(即裂纹面相对位移的方向)可把裂纹分为Ⅰ，Ⅱ，Ⅲ型，如图1所示.

图1 三种基本裂纹模式

设 u，v，w 是裂纹沿 x，y，z方向的位移分量，将裂纹面看成位移的间断面，即空间坐标相同裂纹的2个表面上点的位移是不连续的.

1)Ⅰ型裂纹(张开形).受垂直于裂纹面的拉应力作用，裂纹表面在x-y面内法向张开.

2)Ⅱ型裂纹(滑开形).受平行于裂纹面且垂直于裂纹前缘的剪应力作用，裂纹面在x-z面内沿x轴方向相对滑动.

3)Ⅲ型裂纹(撕开形).受平行于裂纹面和裂纹前缘的剪应力作用，裂纹面在x-z面内沿z轴方向相对滑动.

对任意裂纹前沿为曲线的表面裂纹，在裂纹尖端取坐标原点，局部直角坐标系的x-y平面为裂纹前缘任一点的法平面，y-z平面为其切平面;局部柱坐标系的r→0平面为裂纹前缘任一点的法平面，如图2所示.根据Westergaard复变函数方法可求得线弹性断裂力学中Ⅰ，Ⅱ，Ⅲ型裂纹尖端(r→0)的应力应变场为［6-7］:

1)Ⅰ型

2)Ⅱ型

3)Ⅲ型

图2 1/4节点单元

式中，r，θ为局部柱坐标系中的2个坐标分量;G为剪切弹性模量，G=E/2(1+μ)，E为弹性模量，μ为材料泊松比;σx，σy为法向应力;τxy，τxz，τyz为剪切应力;KⅠ，KⅡ，KⅢ为Ⅰ，Ⅱ，Ⅲ型应力强度因子;k 为与材料泊松比μ有关的常数，k=(3－μ)/(1+μ).

实际结构中的裂纹并不是受单一荷载作用的.因此，由任意外载荷作用所产生的裂纹尖端附近区域的应力场和位移场，可分别由上述的3种独立变形型式的应力场和位移场线形叠加给出，即

在常见的裂纹尖端场表达式中，一般高阶量O(r)，O(r0)都忽略，因为当r与面内尺寸相比很小时，这些高阶量与主导项相比，可以略去不计.因此，通常采用主导项来表达线弹性裂纹尖端区的应力位移场.

由裂纹尖端区域的应力场式(17)和位移场式(18)可知:如果给定KⅠ，KⅡ，KⅢ三个参数，则裂纹尖端区域的应力场和位移场也就可以完全确定.KⅠ，KⅡ，KⅢ与坐标(r，θ)无关，它们表征了靠近裂纹尖端区的应力场与位移场奇异性的强度，其值由裂纹体的几何和所受载荷决定，因此可以用来度量相应的应力应变场强度，从而判断裂纹是否会扩展以及扩展的方向［6，8］.

2 应力强度因子计算方法

2.1 1/4节点退化单元

奇异单元(quarter-point element)是绕裂纹尖端附近构造的一种特殊单元，裂纹尖端区域除采用奇异单元之外，其余采用常规单元.奇异单元的特点是位移插值函数中直接包含所需要的奇异项，因而能够反映裂纹尖端附近的奇异性.目前，已有许多不同的奇异单元用于求解裂纹尖端应力应变场，其中应用较为广泛的是8节点等参数奇异单元，其节点退化如图2所示［6-8］.

2.2 有限元位移法求解应力强度因子

采用有限元法和节点位移技术来计算应力强度因子K的方法，也被称为直接方法.

当 θ= ± π，r=ra－b时，应力强度因子可以用下式计算:

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

式中，ua和ub为y方向在点a和点b上的位移;rab为裂纹尖端到点 b的距离(即点 a到点 b的距离).

有限元法计算应力强度因子的准确度依赖于网格的细度和单元类型，如果采用1/4节点单元，该方法准确度将得到改善.位移沿裂纹边对1/4节点退化单元的值为

式中，vupper，vlower分别为1/4节点沿abc处和ade处的裂纹位移，如图2所示.

根据有限元法计算结果，裂纹张开位移(crack open displacements，COD)为

应力强度因子 K 为［5，7］

KⅠ，KⅡ，KⅢ的计算公式相同，只要代入不同方向的位移量即可计算Ⅰ，Ⅱ，Ⅲ型的应力强度因子.

有限元法计算的流程如图3所示.

图3 有限元分析流程图

3 算例

3.1 计算模型1(中间裂纹)［9-10］

计算模型1如图4所示，弹性模量 E=20 MN/cm2;泊松比μ=0.3;法向应力P=30 N/cm.由于模型对称，只建一半模型.模型网格划分采用8节点四边形等参元和1/4节点退化单元.在裂纹尖端附近网格的细度划分如图5所示.

图4 计算模型1

图5 网格划分与在裂纹尖端附近网格的细度

对于中间裂纹，应力强度因子理论计算公式［9］为

式中，aL为裂纹长度;bL为板跨度;Ktheory为裂纹理论计算强度应子;KFEM为有限元计算强度应子;Kerror为有限元计算与理论计算强度应子的误差.

应力强度因子计算精度与网格划分及1/4节点单元退化的尺寸有关.采用有限元法和理论计算法计算的Ⅰ型应力强度因子如图6所示.

假设裂纹长度 aL=5，10，…，50 cm，则应力强度因子KⅠ与裂纹长度aL之间的关系如图7所示.

假设荷载 q=10，20，…，50 N/cm，应力强度因子KⅠ与不同荷载之间的关系如图8所示.

图6 不同网格数2种计算方法得到的KⅠ比较

图7 应力强度因子KⅠ与裂纹长度之间的关系

图8 不同载荷下的应力强度因子KⅠ

由图6～图8可以看出，精确度(收敛程度)与网格密度、1/4节点单元退化的尺寸和裂纹长度有关.在网格数为54(网格密度低)时，有限元法计算结果与理论值的误差为6.861 9%，提高划分网格密度值，当网格数大于868时，误差小于1.777 9%.当裂纹长度为15 cm时，有限元位移法计算结果与理论结果相比，误差为最大(1.508 5%);当裂纹长度为35 cm时，计算结果与理论结果相比，误差最小(0.004 4%).不同载荷下，有限元位移法计算结果与理论结果计算一致.

3.2 计算模型 2(侧裂纹)［9，11-12］

计算模型2如图9所示，弹性模量 E=20 MN/cm2;泊松比 μ=0.3，切向应力Q=30 N/cm.

图9 计算模型2

对于侧裂纹理论计算公式为［7］

应力强度因子计算精度与网格划分及1/4节点退化单元的尺寸有关.采用有限元位移法和理论计算法计算的Ⅱ型应力强度因子如图10所示.

假设裂纹长度aL为10，20，…，80 cm，则应力强度因子KⅡ与裂纹长度aL之间的关系如图11所示.

假设荷载 Q=10，20，…，50 N/cm，则应力强度因子KⅡ与不同载荷分布力之间的关系如图12所示.

由图8～图12可以看出，KⅡ有限元法计算结果与理论值的误差大部分小于0.773 2%.在网格数为454时，误差小于0.266 0%.不同裂纹长度应力强度因子KⅡ的有限元计算结果与理论结果相比，误差很小.当裂纹长度为80 cm时，有限元位移法计算结果与理论结果相比误差为最大(1.049 1%).当裂纹长度为25 cm时，有限元位移法计算结果与理论结果相比误差为最小(0.061 8%).不同载荷下，有限元位移法计算结果与理论结果计算一致.

图10 不同网格数2种计算方法应力强度因子KⅡ比较

图11 应力强度因子KⅡ与裂纹长度之间的关系

图12 不同载荷下的应力强度因子KⅡ

4 裂纹尖端附近的位移场与应力场计算

除了在裂纹尖端附近网格细度加密外，网格的划分采用8个节点单元类型和1/4节点单元退化.

位移场与应力场的中间裂纹有限元模拟结果如图13所示.位移场与应力场的侧裂纹有限元模拟结果如图14所示.

图13 位移场和应力场的中间裂纹

由位移场和应力场计算结果可见，在裂纹尖端附近存在严重退化，且应力最大.

5 结语

本文基于有限元法和节点位移方法采用8节点四边形等参元和1/4节点退化单元，对中间裂纹和侧裂纹尖端的强度应子进行了计算.采用整体较稀疏的网格划分，但在裂纹尖端附近区域增加网格和节点的密度，能够大大减少网格数量，便于方程组的求解，节省计算时间，同时确保了问题的准确性.

通过有限元模型计算结果与理论计算结果相比表明:有限元模型计算精度与网格密度、1/4节点单元退化的尺寸和裂纹尖端的裂纹长度有关，其中网格划分是最主要的影响因数，提高划分密度可以有效提高计算精度.位移场和应力场计算结果表明，中间裂纹和侧裂纹的裂纹尖端附近存在严重退化，且应力最大.

图14 位移场和应力场的侧裂纹

References)

［1］Liu P F，Hou S J，Chu J K，et al.Finite element analysis of postbuckling and delamination of composite laminates using virtual crack closure technique［J］.Composite Structures，2011，93(6):1549-1560.

［2］Muthu N，Falzon B G，Maiti S K，et al.Modified crack closure integral technique for extraction of SIFs in meshfree methods［J］.Finite Elements in Analysis and Design，2014，78(5):25-39.

［3］Okada Hiroshi，Ohata Shogo.Three-dimensional J-integral evaluation for cracks with arbitrary curvatures and kinks based on domain integral method for quadratic tetrahedral finite element［J］.Engineering Fracture Mechanics，2013，109(9):58-77.

［4］刘明尧，柯梦龙，周祖德，等.裂纹尖端应力强度因子的有限元计算方法分析［J］.武汉理工大学学报，2011，33(6):116-121.Liu Mingyao，Ke Menglong，Zhou Zude，et al.Analysis of finite element calculation methods for crack-tip stress intensity factor［J］.Journal of Wuhan University of Technology，2011，33(6):116-121.(in Chinese)

［5］Fehl Barry D，Truman Kevin Z.An evaluation of fracture mechanics quarter-point displacement techniques used for computing stress intensity factors［J］.Engineering Structures，1999，21(5):406-415.

［6］张俊清.高速列车空心车轴表面裂纹应力强度因子研究［D］.北京:北京交通大学交通运输学院，2011.

［7］Tracy Dennis M.Finite elements for determination of crack tip elastic stress intensity factors［J］.Engineering Fracture Mechanics，1971，3(3):255-266.

［8］Wu Zhixue.The shape of a surface crack in a plate based on a given stress intensity factor distribution［J］.International Journal of Pressure Vessels and Piping，2006，83(3):168-180.

［9］Tada H，Paris P，Irwin G.The stress analysis of cracks handbook［M］.Washington DC，USA:Washington U-niversity，1957:41-52.

［10］Brennan F P，Teh L S.Determination of crack tip stress intensity factors in complex geometries by composition of weight function solutions［J］.Fatigue ＆Fracture of Engineering Materials＆Structures，2004，27(1):1-7.

［11］Muthu N，Maiti S K，Falzon B G，et al.A comparison of stress intensity factors obtained through crack closure integral and other approaches using extended element-free Galerkin method［J］.Computational Mechanics，2013，52(3):587-605.

［12］Salari R H，Rahimi D M，Rahimi P S，et al.Meshless EFG simulation of linear elastic fracture propagation under various loadings［J］.Arabian Journal for Science and Engineering，2011，36(7):1381-1392.