埋藏裂纹应力强度因子的有限元分析

岳应娟，李 宁，陈 飞，孙 钢

(第二炮兵工程大学理学院，陕西 西安，710025)

埋藏裂纹应力强度因子的有限元分析

岳应娟，李 宁，陈 飞，孙 钢

(第二炮兵工程大学理学院，陕西 西安，710025)

采用三维奇异单元模拟裂纹前缘应力应变场的奇异性，建立了计算球形压力容器中埋藏裂纹应力强度因子的有限元模型，有限元分析结果和经验公式计算结果吻合较好，证明所建有限元模型具有较高的准确性和可靠性。通过仿真计算，分析不同因素对应力强度因子的影响，得出了裂纹前缘应力强度因子随裂纹形状、容器壁厚以及裂纹中心与壁厚中心的偏移量与壁厚比值的变化规律。

埋藏裂纹；应力强度因子；球形压力容器；有限元

由于生产工艺和工作环境的影响，球形压力容器焊缝区域不可避免地存在不同程度的裂纹缺陷[1]，而且在长期的贮存使用过程中，压力容器还会因载荷、介质、充放气操作等各种因素的影响而萌生出新的裂纹。实践证明，并非所有的裂纹缺陷都会导致压力容器失效，为此世界各国相继制定了基于“合乎使用原则”的安全评定标准。应力强度因子K是压力容器安全评定中的重要参数，对其影响因素和变化规律进行研究具有重要意义。

目前，一些计算裂纹应力强度因子的经验公式大多是在含穿透裂纹的无限大平板计算模型的基础之上进行适当修正得到的[2]，有时与实际情况差距较大。而有限元法不受裂纹体几何形状及载荷形式的限制，因而在断裂力学研究中得到非常广泛的应用。本文在ANSYS软件的基础上，对裂纹应力强度因子的有限元求解、影响因素及变化规律进行研究。

1 基本理论

ANSYS中采用节点位移法计算应力强度因子，即先计算出位移，然后通过后处理命令完成应力强度因子的计算。根据线弹性断裂理论，三维裂纹前缘区域的位移场可表示为[3]：

(1)

式中：u、v、ω为裂纹前缘局部坐标的位移；KⅠ、KⅡ、KⅢ分别为Ⅰ型、Ⅱ型、Ⅲ型裂纹的应力强度因子；G为材料的切变模量；r、θ为裂纹前缘的极坐标；k是与泊松比μ有关的系数，其关系式为：

(2)

将θ=π代入式(1)，得到裂纹前缘附近各点的位移公式为：

(3)

从而得到：

(4)

式(4)中的u、v、ω可以通过有限元法计算得到。

由于裂纹前缘的应力和应变具有奇异性[4]，因此在进行有限单元建模或单元网格划分时，必须先在裂纹前缘定义奇异点，而且围绕奇异点的有限单元是二项式的奇异单元。本文选取20节点SOLID95六面体单元[5]，单元中心节点外推至1/4位置处，六面体单元退化成三棱柱体，以模拟裂纹前缘奇异区域及外围区域，如图1所示。

2 基于ANSYS的仿真计算

本文选取的研究对象是在球形压力容器焊缝中存在较多的埋藏裂纹。对埋藏裂纹进行研究时，一般将其规则化为椭圆裂纹，椭圆裂纹长轴和短轴的尺寸由埋藏裂纹外接矩形的长和高来确定[6]，见图2。图3为球形压力容器内部椭圆埋藏裂纹示意图，其中，B为球形压力容器壁厚，2a为裂纹高度，2c为裂纹长度，e为埋藏裂纹中心与壁厚中心的偏移量；A1代表靠近容器内壁的椭圆裂纹短轴端点，A2代表靠近容器外壁的椭圆裂纹短轴端点，C1、C2代表椭圆裂纹长轴端点。

球形压力容器主要承受内部高压气体产生的内压，其它外力可以忽略不计。由于裂纹和球形容器内部具有对称性，其受力也具有对称性，故仅需建立1/4模型，如图4所示。图5为裂纹区域模型。球形压力容器内径Di=1400 mm，壁厚B=75 mm，取椭圆裂纹尺寸c=15 mm，a=7.5 mm，偏移量e=0。材料为15MnMoVN，本构模型为线性、各向同性。15MnMoVN在室温条件下的弹性模量为196 GPa，泊松比为0.3。

Fig.3 Schematic diagram of embedding crack in spherical pressure vessel

对裂尖所在表面进行二维网格划分[7-8](裂纹尖端附近利用ANSYS前处理命令KSCON生成二维奇异单元)，采用扫略分网VSWEEP将所生成的面网格沿裂纹前缘进行扫略，生成三维奇异单元，如图6所示。对裂纹区以外部分进行自由网格划分，生成带裂纹球形压力容器有限元模型，如图7所示。对1/4模型内表面施加30MPa均布载荷，分别对面1、面2施加关于x-y平面、x-z平面的对称约束。

对模型求解后进行后处理，得到裂纹前缘各节点处的应力强度因子KⅠ最大值为528.65 MPa·mm1/2，KⅡ、KⅢ均为0，这是因为在只受内压的情况下，裂纹为Ⅰ型张开型，不存在由Ⅱ型滑开型和Ⅲ型撕开型引起的KⅡ、KⅢ值。根据GB/T 19624—2004[6]附录中含椭圆埋藏裂纹的板壳应力强度因子经验公式计算得到的KⅠ值为568.74 MPa·mm1/2。两者误差为7.04%，满足计算精度及安全储备要求。有限元计算结果小于板壳经验公式计算结果，主要是由于球壳曲率的影响，裂纹实际轴线会偏离椭圆形的轴线，使得所建有限元模型的椭圆率和相对深度都发生了改变，另外用ANSYS计算应力强度因子时，划分的网格尺寸大小也会对结果产生一定的影响。通过对网格进一步细化，所得结果与上述结果相差很小，表明文中网格划分具有足够的精度，能够保证分析结果的准确性。

3 应力强度因子的影响因素分析

3.1 a/c值

在固定椭圆裂纹短轴尺寸a值不变的条件下，改变裂纹圆度(a/c值)，保持其它参数不变，建立有限元模型，计算应力强度因子。

根据仿真计算结果绘制裂纹前缘不同位置处(图3中A1、A2、C1)的应力强度因子随a/c值的变化趋势图，见图8。

从图8可以看出：椭圆埋藏裂纹短轴端点的KⅠ值大于长轴端点的KⅠ值；靠近容器内壁短轴端点的KⅠ值大于靠近容器外壁短轴端点的KⅠ值，这是因为裂纹面受到的拉应力是非均匀的，且内壁处拉应力大于外壁处拉应力；随着a/c值的增大，即裂纹由椭圆形逐渐向圆形过渡，短轴端点处KⅠ值呈减小趋势，而长轴端点处KⅠ值呈增大趋势，这表明对于相同长度的裂纹，条状裂纹较圆形裂纹危害更大。

Fig.8 Stress intensity factor of cracks with differenta/cvalues

3.2 B值

在75～105 mm范围内改变球形压力容器壁厚B值，保持其它参数不变，得到裂纹前缘不同位置处的应力强度因子随B值的变化趋势，如图9所示。

Fig.9 Stress intensity factor of cracks with differentBvalues

从图9可以看出：随着B值的增大，埋藏裂纹短轴端点处和长轴端点处的KⅠ值均呈减小趋势，且减小趋势逐渐变缓。也就是说，在存在相同裂纹的情况下，壁厚越大，压力容器越安全。

3.3 e/B值

固定容器壁厚B值不变，改变埋藏裂纹中心与壁厚中心的偏移量e值(裂纹中心靠近内壁时，e为正；裂纹中心靠近外壁时，e为负)，得到裂纹前缘不同位置处的应力强度因子随e/B值的变化趋势，如图10所示。

从图10可以看出：总的来说，e/B值越大，即裂纹越靠近容器内壁，短轴端点处和长轴端点处的KⅠ值就越大，这表明越靠近容器内壁的裂纹对压力容器的安全所产生的危害就越大。

Fig.10 Stress intensity factor of cracks with differente/Bvalues

4 结语

本文运用ANSYS进行了埋藏裂纹应力强度因子的仿真计算，从裂纹圆度、球形压力容器壁厚以及裂纹中心与壁厚中心的偏移量与壁厚的比值这三个方面分析了其对应力强度因子的影响，结果表明采用有限元法可以很好地模拟计算不同条件下的应力强度因子，得到的应力强度因子变化规律与文献[6]的公式推导结果基本一致。

[1] 周嘉炜，江楠.球罐安全评定方法的比较[J].化工机械，2009，36(5)：521-524.

[2] 陈家权，沈炜良，徐家园，等.应力强度因子的有限元计算[J].装备制造技术，2003(4)：6-9.

[3] Suresh S.材料的疲劳[M]. 王中光，等译.北京：国防工业出版社，1993.

[4] Akin J E. The generation of elements with singularities[J]. International Journal for Numerical Methods in Engineering,1976, 10(6): 1249-1259.

[5] 祝效华，余志祥.ANSYS高级工程有限元分析范例精选[M].北京：电子工业出版社,2004.

[6] 全国锅炉压力容器标准化技术委员会.GB/T 19624—2004 在用含缺陷压力容器安全评定[S].北京：中国标准出版社，2005.

[7] 瞿伟廉，鲁丽君，李明.带三维穿透裂纹结构的有限元实体建模方法[J].武汉理工大学学报，2008，30(1)：87-90.

[8] 康国政.大型有限元程序的原理、结构与使用[M].成都：西南交通大学出版社，2004.

[责任编辑 尚 晶]

Finite element analysis of stress intensity factorof embedding cracks

YueYingjuan,LiNing,ChenFei,SunGang

(College of Science,the Second Artillery Engineering University, Xi’an 710025, China)

3D singular elements are employed to describe the singularity of stress-strain field of crack front, and the finite element model for calculating stress intensity factor of embedding cracks in spherical pressure vessel is established. The FEM results agree well with the data calculated from empirical formula, which proves the built finite element model is accurate and reliable. The influences of different factors on the stress intensity factor are analyzed by simulation, and the varying patterns of stress intensity factor at the crack front with crack shape, wall thickness of vessel and the ratio of the offset of crack center and wall center to wall thickness are obtained.

embedding crack；stress intensity factor；spherical pressure vessel；finite element

2014-12-04

岳应娟(1972-)，女，第二炮兵工程大学教授，博士.E-mail: yingjuanyue@163.com

O346.1

A

1674-3644(2015)02-0125-04