基于三维有限元解的紧凑拉伸试样应力强度因子计算公式

贾 旭，胡绪腾，宋迎东

(南京航空航天大学能源与动力学院，江苏省航空动力系统重点实验室，

机械结构力学及控制国家重点实验室， 南京 210016)



基于三维有限元解的紧凑拉伸试样应力强度因子计算公式

贾 旭，胡绪腾，宋迎东

(南京航空航天大学能源与动力学院，江苏省航空动力系统重点实验室，

机械结构力学及控制国家重点实验室， 南京 210016)

摘要：对标准紧凑拉伸(CT)试样的二维和三维应力强度因子有限元解进行了对比分析，基于三维有限元虚拟裂纹闭合法拟合了一种新的标准CT试样的应力强度因子计算公式，并采用ANSYS软件内嵌位移外推法进行了验证。结果表明：基于二维分析所得的应力强度因子计算公式计算结果与实际三维CT试样的具有较大的差异；在加载孔等效分布力一定的条件下，CT试样裂纹前沿大部分区域的应力强度因子与中心点的相近，且与厚度无关；拟合得到的三维CT试样裂纹前沿中心点的应力强度因子计算公式具有很高的精度，其计算结果与ANSYS软件内嵌位移外推法的相对误差在0.5%以内。

关键词：三维分析；CT试样；应力强度因子；有限元法

0引言

为提高航空发动机的安全性，20世纪80年代以来损伤容限设计已逐渐成为国内外航空发动机轮盘等关键件的重要设计准则之一[1-2]。损伤容限设计的关键是对以断裂为主要破坏形式的关键件进行裂纹扩展分析，而裂纹扩展分析的基础是应力强度因子的计算。

工程中普遍采用标准紧凑拉伸(CT)试样进行试验来建立材料裂纹扩展速率模型、确定材料断裂韧性。1974年Newman[3]基于二维平面假设(平面应变假设)下的改进边界配置法拟合得到了标准CT试样在不同几何配置下的应力强度因子计算公式。1976年Srawley[4]对此公式进行了改进，拓宽了标准CT试样应力强度因子计算公式的适用范围。目前ASTM E399-12e3《平面应变断裂韧性测试标准》和ASTM E647-13ae1《疲劳裂纹扩展速率测试标准》仍在继续使用Srawley改进的标准CT试样应力强度因子计算公式。非标准试样的裂纹扩展速率研究[5-6]也仍在使用应力强度因子手册中的二维应力强度因子表达式。然而，二维平面假设毕竟存在一定的局限性，在实际的三维裂纹体中，即使是标准CT试样中形状规则、初始裂纹前缘呈直线、承受单一销钉载荷的张开型穿透裂纹，其裂纹前沿应力场也呈现出不同的应力状态。特别是裂纹与自由表面相交的区域，既不是平面应变状态也不是平面应力状态[7-8]。二维平面应力强度因子计算公式将整个裂纹前沿生硬地假设为平面应变或平面应力两种极端情况，势必对应力强度因子的计算造成不可避免的误差。目前少有学者对三维标准CT试样的应力强度因子计算公式进行研究，Towers[9]计算对比了裂纹前缘中心点、1/4厚度点、表面点的应力强度因子的不同；Neman[10]在平面应力状态假设下研究了CT试样的权函数。为了得到标准CT试样三维条件下的应力强度因子计算公式，作者建立了二维平面应变状态的CT试样的有限元模型，通过不同方法计算了应力强度因子并与目前普遍采用的应力强度因子计算公式进行了对比；然后建立了CT试样的三维有限元模型，得到了CT试样三维应力强度因子分布；最后基于CT试样的三维有限元解拟合得到了一种新的CT试样应力强度因子表达式，并与ANSYS软件内嵌的位移外推法进行了对比分析。

1CT试样应力强度因子的二维有限元分析

目前在ASTM标准中广泛采用的CT试样应力强度因子计算公式为：

(1)

式中：K为应力强度因子；a为裂纹长度；W为CT试样的宽度；B为CT试样的厚度(平面状态假设下B=1)；P为施加在加载销钉上的力。

以TC4合金为试样及销钉材料(弹性模量E为109 GPa,泊松比ν为0.34)，建立了宽度为0.04 m的CT试样的两种二维平面应变有限元1/2模型。第一种(①)在裂纹尖端采用六节点奇异单元plane82网格划分，单元个数2 837；第二种(②)在裂纹尖端采用八节点高阶单元plane82分网，单元个数3 007，如图1所示。图1中包围裂纹尖端的16个奇异单元的特征尺寸l为0.01W。为了说明所建立的有限元模型、加载方式、应力强度因子计算方法的准确性，分别采用了ANSYS软件内嵌位移外推法[11](Displacement Extrapolation Method，简称DEXM)和虚拟裂纹闭合法[12](Virtual Crack Closure Technique，简称VCCT)计算了两种不同加载方式下的应力强度因子，并分别与式(1)结果进行对比，见表1。耦合加载方法将加载孔上表面内所有节点与加载孔中心的质点进行刚性耦合后对质点进行集中力加载；接触加载方法是通过建立销钉模型，进行接触求解。相对于耦合加载，接触加载方式为非线性求解，虽然计算效率较低但最切合实际。表1中δ1，δ2，δ3分别为奇异单元划分/耦合加载时DEXM、奇异单元划分/接触加载时DEXM、非奇异单元划分/接触加载时VCCT计算结果与式(1)的相对误差。

图1　CT试样二维有限元模型Fig.1　2D finite element model about the CT specimen：(a) general graph and (b) local graph at crack front

从表1可以看出，奇异单元划分/耦合加载时，DEXM计算结果与式(1)的相对误差在±2%以内。奇异单元划分/接触加载下DEXM与非奇异元划分/接触加载VCCT计算结果相对式(1)的误差较一致，均在±1%以内。这一结果不仅说明了接触加载的方式更为合理并验证了有限元模型和方法的正确性，也说明了以平面假设为基础的式(1)在计算CT试样二维平面K时具有足够高的精度。

表1　CT试样二维平面应变假设下的应力强度因子计算结果及相对误差Tab.1　 Calculated stress intensity factors of CT specimens and relative errors under two-dimensional plane strain assumption

2CT试样应力强度因子的三维有限元分析

三维模型CT试样与销钉的材料及参数与二维分析时相同。CT试样的宽度W为0.04 m，宽厚比W/B分别为4，6，8，10，12，14，16，20。分三个区域对CT试样的三维模型进行了划分：靠近裂纹尖端采用五面体奇异单元或20节点等参元划分，远离裂纹尖端采用20节点等参元划分，过渡区采用四面体等参元划分，单元为solid 95，如图2所示。第一种(①)网格类型的单元数目为12 210，第二种(②)网格类型的单元数目为16 194。为了说明CT试样二维平面应力强度因子与三维模型的不同，首先以a/W为0.3绘制了不同宽厚比下裂纹前缘假设为平面应变与平面应力状态下的应力强度因子分布，如图3所示。其中平面应变状态和平面应力状态并非指三维裂纹尖端的应力状态，而是指应力强度因子数值计算方法的前提条件。然后以a/W为0.3，W/B为10为例对比了三维应力强度因子与式(1)的相对误差，见表2。其中坐标轴Z平行于裂纹线，Z/(B/2)=1表示CT试样的中心位置,δ4为平面应变状态下的结果与式(1)的相对误差，δ5为平面应力状态下的结果与式(1)的相对误差。从图3可以看出，以二维平面假设为基础的经验公式显然与三维平面应变或平面应力状态下的应力强度因子相差较大；同一裂纹尺寸、不同宽厚比下CT试样的应力强度因子在裂纹体表面时各有差异；在靠近裂纹体中心时，逐渐汇聚到一个数值。

图2　CT试样的三维有限元模型Fig.2　The 3D finite element mode of CT specimen:(a)general graph and local graph and (b) mesh at crack front

表2　CT试样三维有限元解与式(1)结果的对比Tab.2　Comparison of 3D finite element solutions and formula (1) results

图3　a/W为0.3时不同厚度下的K分布Fig.3　Stress intensity factor distribution at differentthicknesses with a/W equal to 0.3

从表2可以发现，虽然三维CT试样中心区域接近为平面应变状态，但中心区域内的应力强度因子与式(1)的结果也有很大的差别，最大相对误差达到10.36%。虽然三维CT试样表面点为平面应力状态，但与式(1)的结果也有很大的差别，最大相对误差达到-16.12%。

3基于三维有限元解的CT试样应力强度因子计算公式

由式(1)可以得到：

(2)

由式(2)中可以看出,P/(0.125W×2×B)为P等效在加载孔上的分布力，当加载孔上的等效分布力一定时，CT试样的应力强度因子只与W和a有关，与试样的厚度无关。为考证这一特性，并定性地描述不同厚度下三维CT试样裂纹前沿的应力强度因子分布特性，选用了计算效率更高的耦合加载方式以及ANSYS软件内嵌的DEXM法，在同一加载孔上等效分布力下，分别采用平面应变状态计算了三维CT试样W/B分别为4，6，8，10，12，14，16，18，20，a/W分别为0.2，0.3，0.4，0.5，0.6，0.7，0.8下的应力强度因子分布，如图4所示。

图4　三维CT试样裂纹前沿的应力强度因子分布Fig.4　Stress intensity factor distribution at the crack frontof the 3D CT specimen

由图4发现，对于同一裂纹尺寸、不同宽厚比的CT试样，位于试样中心区域的应力强度因子基本分布在同一直线上，直线区域约占整个厚度的2/3。由此可见，在加载孔等效分布力一定的条件下，CT试样裂纹前沿大部分区域的应力强度因子与中心点的应力强度因子接近，并且与厚度无关。

基于以上结论，作者在同一加载孔等效分布力的基础上，采用了更精确的接触加载方式，并采用了对网格密度依赖性较低的非奇异单元VCCT法计算了CT试样中心点的应力强度因子。以a/W为0.5，W/B为10为例，采用裂纹尖端不同单元尺寸(l)进行网格划分验证了应力强度因子计算的网格无关性，见表3。其中δ′表示表中后一项K值相对于前一项的误差。

由表3可知，l/a为0.01的K值相对于0.02时的误差仅有0.158%，0.005的相对于0.01时的误差仅有0.079%，相对误差较小并且随l/a的减小而减小。说明当l/a为0.01时，有限元法已经达到预期收敛效果，其计算结果与网格尺寸无关，为了方便划分网格和提高计算效率，选用l/a为0.01。

表3　网格无关性验证Tab.3　Mesh-independent verification

随后，采用既定划分网格方式计算了不同裂纹尺寸下的应力强度因子，并与式(1)的进行了对比，如图5所示。从图5中可以看出三维CT试样中心部位的应力强度因子随裂纹尺寸的变化趋势与式(1)的基本一致。

图5　VCCT法计算结果与式(1)结果的对比Fig.5　Comparison of VCCT calculation results andformula (1) results

鉴于目前采用的应力强度因子计算公式不适用于三维CT试样，为了建立一种准确的适用于三维CT试样的应力强度因子计算公式，在式(1)的基础上对VCCT法的计算结果进行了应力强度因子的拟合，结果如下：

(3)

为了验证式(3)的精度，对比了VCCT法(拟合数据点)的结果和式(1)与式(3)的相对误差，并采用接触加载下的奇异单元DEXM法对式(3)进行了验证，结果见表4。表中δ6为式(1)的结果相对式(3)的误差，δ7为VCCT法(拟合数据点)的结果相对式(3)的误差，δ8为DEXM数据相对式(3)的误差。

从表4可以看出，式(1)结果与式(3)的相对误差在-7.79%～-11.75%，式(1)结果比式(3)的小；式(3)结果与VCCT法的误差在±0.04%以内，说明式(3)具有很高的拟合精度；DEXM法的结果与式(3)的相对误差在±0.5%以内，说明新的CT试样应力强度因子表达式(3)在计算三维CT试样应力强度因子时具有足够高的精度。

表4　三维CT试样裂纹中心部位的应力强度因子的计算结果与相对误差Tab.4　Calculated stress intensity factors and relative errors at the crack center of the 3D CT specimen

4结论

(1) 目前广泛采用的应力强度因子计算公式是通过二维假设分析所得，其计算结果与实际三维CT试样的有很大差异，最大误差达到10%以上。

(2) 在加载孔等效分布力一定的条件下，CT试样裂纹前沿大部分区域的应力强度因子与中心点的相近，并且与厚度无关。

(3) 拟合得到的三维CT试样裂纹前沿中心点的应力强度因子计算公式具有很高的精度，相对于其他方法的误差在0.5%以内。

参考文献：

[1]原航空工业部第三零一研究所.航空发动机结构完整性指南:GJB/Z 101-97[S].北京:[出版者不详],1997：36-44.

[2]吕文林，陈俊粤，田德义. 航空涡喷、涡扇发动机结构设计准则(研究报告)第二册：轮盘[R]. 北京: 中国航空工业总公司发动机系统工程局，1997.

[3]NEWMAN J C. Stress analysis of compact specimens including the effects of pin loading[J].ASTM STP,1974，560:105-105.

[4]SRAWLEY J E. Wide range stress intensity factor expressions for ASTM method E399 standard fracture toughness specimens[J]. International Journal of Fracture, 1976，12(3): 475-476.

[5]李旭东，穆志韬，贾明明. 加载频率对航空铝合金腐蚀疲劳裂纹扩展速率的影响[J].机械工程材料,2014，38(7)：50-52.

[6]余圣甫, 王铁琦, 杨其良, 等. ZG20MnSi钢断裂韧度和疲劳裂纹扩展速率试验研究[J].机械工程材料,2005,29(6):17-19.

[7]GARCIA-MANRIQUE J, CAMAS D, LOPEZ-CRESPO P, et al. Stress intensity factor analysis of through thickness effects[J]. International Journal of Fatigue, 2013, 46: 58-66.

[8]BAZANT Z P, ESTENSSORO L F. Surface singularity and crack propagation[J]. International Journal of Solids and Structures, 1979, 15(5): 405-426.

[9]TOWERS O L, SMITH A P. Stress intensity factors for curved crack fronts in compact tension specimens[J]. International Journal of Fracture, 1984, 25(2): 43-48

[10]NEWMAN J C, YAMADA Y, JAMES M A. Stress-intensity-factor equations for compact specimen subjected to concentrated forces[J]. Engineering Fracture Mechanics, 2010, 77(6): 1025-1029.

[11]SLAVIK D C, MCCLAIN R D, LEWIS K. Stress intensity predictions with ANSYS for use in aircraft engine component life prediction[J]. Fatigue and Fracture Mechanics, 2000, 31: 371-390.

[12]解德，钱勤，李长安. 断裂力学中的数值计算方法及工程应用[M]. 北京: 科学出版社, 2009.

Calculation Formula for Stress Intensity Factors of CT Specimens

based on Three Dimensional Finite Element Solutions

JIA Xu, HU Xu-teng, SONG Ying-dong

(Jiangsu Province Key Laboratory of Aerospace Power Systems,

State Key Laboratory of Mechanics and Control of Mechanical Structures,

College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

Abstract:Two-dimensional and three-dimensional finite element calculations of stress intensity factors for standard compact tension(CT) specimens were compared and analyzed. And a new calculation formula of stress intensity factors for standard CT specimens was established based on the three-dimensional finite element virtual crack closure method and verified using displacement extrapolation method of ANSYS. The results show that great differences existed between the stress intensity factors from calculation formula based on the two-dimensional analysis and those of actual three-dimensional CT specimens. Under the certain equivalent distributed force on the load-hole, the stress intensity factors at the most crack front regions were similar to those at the center of the CT specimen and were independent of the thickness. The fitting calculation formula of stress intensity factors at the crack front center of three-dimensional CT specimens was of high precision. And the relative errors between the fitting calculation values and the results of the displacement extrapolation method of ANSYS were less than 0.5%.

Key words:three dimension analysis; CT specimen; stress intensity factor; finite element method

通讯作者：李落星教授

作者简介：王冠(1985-)，男，河南郑州人，讲师，博士。

基金项目：国家自然科学基金面上资助项目(51475156)；宁夏大学自然科学研究基金资助项目(ZR1403)；宁夏大学人才引进科研启动基金资助项目(BQD2014018)

收稿日期:2015-06-02；

修订日期:2015-10-23

DOI:10.11973/jxgccl201512009

中图分类号：O346.1

文献标志码：A

文章编号：1000-3738(2015)12-0030-05