裂纹应力强度因子的有限元计算

第３０卷第５期　西安科技大　学　学报　Ｖｏ１．３０　Ｎｏ．５　２０１０年９月　ＪＯＵＲＮＡＬ　ＯＦ　ＸＩ　ＡＮ　ＵＮＩＶＥＲＳＦＩ３￣ＯＦ　ＳＣＩＥＮＣＥ　ＡＮＤ　ＴＥＣＨＮＯＬＯＧＹ　Ｓｅｐｔ．２０１０　文章编号：１６７２—９３１５（２０１０）０５—０６２９—０４　裂纹应力强度因子的有限元计算　乔宝明　（西安科技大学理学院，陕西西安７１００５４）　摘要：基于解析法和数值法，对有限长平板中存在的中心穿透裂纹，分析其附近的位移场、应力　场分布，得到计算裂纹尖端的应力强度因子的公式。通过对求得的应力强度因子值与经验值的　比较，表明本文给出的计算应力强度因子的方法具有较好的精度。　关键词：裂纹；应力强度因子；Ｆｏｕｒｉｅｒ变换；有限元；最小二乘拟合　中图分类号：Ｏ　３４６．１　文献标志码：Ａ　０　引　言　传统强度理论是在假设材料无缺陷，无裂纹的基础上建立起来的，在生产实践中经受了长期的考验，　但是在有些情况下，运用传统强度理论设计的结构，却仍然出现了断裂事故。断裂处的应力往往不高，甚　至远远低于材料的屈服极限，这说明传统的强度和韧度指标及计算结果虽然能满足设计要求，但并不能　保证结构的安全。在工程实际中，机械设备和金属结构的构件中往往存在因制造、使用、或材料本身缺陷　所致的宏观裂纹，这时，要确定构件能否继续安全使用，最重要的就是判断裂纹是否会失稳扩展从而导致　结构和设备的破坏。近代断裂力学正是Ｉｒｗｉｎ在对裂纹尖端的应力场进行分析后，引出了应力强度因子　的概念，并与Ｇｒｉｆｉｆｔｈ能量观点联系之后逐步发展起来的¨Ｊ。　应力强度因子ＳＩＦ（ｓｔｒｅｓｓ　ｉｎｔｅｎｓｉｔｙ　ｆａｃｔｏｒ）是断裂力学中表征裂纹尖端应力应变场强度的一个极为重　要的参数，与裂纹尖端能量释放率有着必然的联系。按断裂力学的观点：裂纹尖端的应力强度因子　若　小于材料的断裂韧性　，则构件是安全的，否则，构件是危险的　Ｊ。因此，如何准确、有效地求得构件裂纹　尖端的应力强度因子一直是工程技术人员关注的问题。在应力强度因子的计算中，利用解析法（例如　Ｆｏｕｒｉｅｒ方法）从理论上说可以得到精确解。但是，在工程实际中，由于结构和构件的形状及受力的复杂性　以及裂纹形态的多样性，解析法往往在数学上难以描述和求解。而有限元法　虽然是一种数值计算的近　似方法，但由于其具有强大的建模功能和能够充分利用计算机的计算能力，并且随着网格划分的不断细　化，可以不断接近精确解，从而满足工程实际的需要。　文中首先介绍Ｆｏｕｒｉｅｒ变换　求解应力强度因子解析表达式的方法，并将通过实例介绍基于有限元理　论的计算裂纹尖端应力强度因子的方法。　１　Ｆｏｕｒｉｅｒ变换法求解应力强度因子　在弹性力学理论中，常体力下弹性平面问题存在的应力函数　，称为Ａｉｒｙ应力函数，为双调和函　数　引。对于平面问题可以采用Ｆｏｕｒｉｅｒ变换或者Ｌａｐｌａｃｅ变换来计算应力强度因子。所需求解方程为　Ｖ　＝０．　（１）　记应力函数　关于变量　的Ｆｏｕｒｉｅｒ变换为　（　，ｙ）即　・收稿日期：２０１０—０４—２７　基金项目：陕西省教育厅专项基金项目（２００８ＪＫ３６６）；西安科技大学培育基金项目（２００６４５）；陕西省教育厅专项基金项目　（２００８ＪＫ３６５）　作者简介：乔宝明（１９６２一），男，河南省长葛人，副教授，主要从事数学教学及偏微分方程、小波理论应用方向的研究．　６３０　西安科技大　学　学报　２０１０年　（　，Ｙ）＝Ｊ　（　，Ｙ）ｅ－ｉｗｘｄｘ．　（２）　于是利用Ｆｏｕｒｉｅｒ变换微分性质得到　Ｊｆ　一∞　Ｖ　　（４　（　）ｘ，Ｙｙ）　＝ｆｅｉＷＸｄｘ＝（、　熹　）＾厂　，一　１　（　，）＝．，Ｙ）＝０．　（３）　求解此微分方程可以得到　（Ｗ，Ｙ）＝（ｃｌ＋ｃ２ｙ）ｅ一　＋（ｃ３＋ｃ４Ｙ）ｅ　，　（４）　因此应力分量的Ｆｏｕｒｉｅｒ变换为　＝　０－　（　ｄｙ　，ｙ），　＝一／．６，２　（　，，，），　（５）　蒡　．　利用Ｆｏｕｒｉｅｒ反变换即可得到　，　，　，　等。进一步还可计算得到位移场等。以椭圆形裂缝问题为例，　则可以得到应力强度因子　＝ｌｉ　２１Ｔｒ　（　，０）＝　￣／叮ｒ　．　（６）　在给定材料参数的前提下，则可进一步计算得到Ｋ　２　计算应力强度因子的数值方法　裂纹尖端应力强度因子的有限元计算方法有位移法，应力法和．，积分法［６］等，前面两种也称为直接　法，后者称为间接法。　２．１位移法　所谓位移法，就是先利用有限元方法求得裂纹尖端附近节点位移的数值解，再利用位移场的表达式　得到相应节点的应力强度因子，最后利用数值方法外推即可得到裂纹尖端的应力强度因子。下面给出位　移法的计算方法及过程。不妨假设已由有限元方法计算得到位移“和　，则根据断裂力学理论，裂纹尖端附　近位移场的表达式为　＝　・　【ｃ２Ｋ－１　ｓ　３—０］，　㈩　＝　・　［ｃ２Ｋ＋ｌ，ｓｉｎ　３—０］．　㈣　其中Ｅ为弹性模量，且平面应变时Ｋ＝３—４ｇ，平面应力时Ｋ＝孚　．　于是根据（８）式可得距裂纹尖端为ｒ，角度为　时的应力强度因子值为　２Ｅ　／２叮Ｔ　～／——　‘　（９）　其中　为沿任一方向　距离裂纹尖端为ｒ处的节点上Ｙ方向位移的数值解，也可以取作沿不同　方向距离　裂纹尖端为ｒ的节点上　的平均值。这里值得注意的是上述位移场表达式成立要求ｒ较小，而普通有限元单　元在ｒ较小时则无法反映裂纹尖端的奇异性。　事实上，从理论上来说，（９）式的值仅在卜斗０时才是真实的。即　２Ｅ　／２竹　一，一　Ｋ一～　（…ｌｉｍ　÷　５一ｇ）８　ｓｉｎ导ｓｉ一．“　３０　　（１０）由于节点位移的有限元解在ｒ＝０时误差很大，因此无法直接利用（１０）式来进行裂纹尖端（ｒ＝０处）　第５期　乔宝明：裂纹应力强度因子的有限元计算　６３１　的应力强度因子计算。因此，可以先求出裂纹尖端附近节点处的应力强度因子，然后再利用数值方法外推　得到裂纹尖端的应力强度因子。应力法的计算方法和过程同理可得。　２．２　Ｊ积分法　由Ｒｉｃｅ提出的．，积分法是计算围绕裂纹尖端的与路径无关的线积分（．，积分），其物理意义即为裂纹　尖端的能量释放率，再利用相应公式间接的求取应力强度因子，故称之为间接法。由此可以得到裂纹扩展　的应变能为　Ｊ＝ｈｍｆ（Ｗｄｙ一　Ｏ　ｕｄｓ），　（１１）　Ｉ　—ｒ　其中　，是积分路径；Ｗ为应变能密度；ｒＯ为积分路径边界上的应力　－　●　●　　ｌ　ｌｈ　矢量；ｄｓ为积分路径上的位移；　，），为裂纹的局部坐标；　为位移矢　　●　＿－－＂ｉ　２ａ；卜　量。　．．　．．．　．，．Ｊ－—．．．一　—　一　Ｉ　　Ｉ　●弹性情况下，ｌ，积分与ＳＩＦ有如下关系　Ｉ　ｆ　　●　Ｉ一２ｂ：一　ｈ　●　　Ｉ　＿一ｆ　．＿●●●　一　中心裂缝，受到均匀对称拉伸作用，平面应力状态下，　ｃｒ０　材料的弹性模量为Ｅ＝２０．７　ＧＰａ，泊松比为　：　０．３，　ｆｆｆｆｆｆｆｆｆ　０：３．８９０　Ｐａ．　显然此构件是对称的，只需对其１／４部分（不妨取　右上角部分，如图２（ａ）所示）作出分析即可，图２，（ｂ）　为利用数学软件ｍａｔｌａｂ对构件的１／４部分所作的三角　有限单元剖分。　由此进行有限元求解（这里采用应力法），得到裂　（ｂ１　纹尖端附近的　，然后代人　图２　有限长板条的１／４部分及其三角有限单元剖分　：　．　（１３）　Ｆｉｇ．２　１／４　ｐａｎｓ　ｏｆ　ｆｉｎｉｔｅ－ｌｅｎｇｔｈ　ｓｔｒｉｐｓ　ａｎｄ　进行计算。由于在ｒ较小的时候可以取ＫＩ　＝Ａ＋　，故　ｔｒｉａｎｇｕｌａｒ　ｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ｓｕｂｄｉｖｉｓｉｏｎ　裂纹尖端应力强度因子为　ｌ｝Ｉ　：！Ｆｉｍｇｌ４４Ｊ　＊（ｒ）＝Ａ．　（１４）　其中ＫＩ＊ｒ）可以用（１３）式近似计算得到，即Ｋｌ＊ｒ）：　．这相当于对所计算得到的一系列　（ｒ）进行最小二乘拟合。具体计算结果见表１和图３．　表１三角有限单元模型计算应力强度因子　Ｔａｂ．１　Ｔｒｉａｎｇｕｌａｒ　ｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ｍｏｄｅｌ　ｆｏｒ　ｃａｌｃｕｌａｔｉｏｎ　ｏｆ　ｓｔｒｅｓｓ　ｉｎｔｅｎｓｉｔｙ　ｆａｃｔｏｒ　图３应力强度因子数值解和经验值　Ｆｉｇ．３　Ｎｕｍｅｒｉｃａｌ　ｓｏｌｕｔｉｏｎ　ａｎｄ　ｅｘｐｅｒｉｍｅｎｔａｌ　ｖａｌｕｅ　ｏｆ　ｓｔｒｅｓｓ　ｉｎｔｅｎｓｉｔｙ　ｆａｃｔｏｒ　６３２　西安科］｛技心　大　学　学报　］Ｊ　１　２０１０生　］Ｊ　从表１结论可知，与经验方法得出的应力强度因子相比，通过三角有限单元法求得数值解并由最小　二乘拟合得到的裂纹尖端应力强度因子　存在误差，误差百分比与裂纹的宽径比ａ／ｂ有关，ａ／ｂ比值在　０．４附近时，采用本文数值方法求得的应力强度因子误差最小。当　解出后，根据需要可以继续对　进　行分析计算。　４结论　通过对以上计算过程和计算结果的分析比较，可知本文采用的裂纹尖端应力强度因子的三角有限单　元法计算方法具有良好的实际可行性和较高的精度。　如何针对其他更为复杂的构件和裂纹，计算应力强度因子，还有待于进一步的研究。　参考文献　Ｒｅｆｅｒｅｎｃｅｓ　褚武扬．断裂力学基础［Ｍ］．北京：科学出版社，１９７９．　ＣＨＵ　Ｗｕ—ｙａｎｇ．Ｔｈｅ　Ｂａｓｉｓ　ｏｆ　Ｆｒａｃｔｕｒｅ　Ｍｅｃｈａｎｉｃｓ［Ｍ］．Ｂｅｒｉｎｇ：Ｓｃｉｅｎｃｅ　Ｐｒｅｓｓ，１９７９．　曹毅中．工程断裂力学［Ｍ］．西安：西安交通大学出版社，１９９１．　ＣＡＯ　Ｙｉ—ｚｈｏｎｇ．Ｅｎｇｉｎｅｅｒｉｎｇ　Ｆｒａｃｔｕｒｅ　Ｍｅｃｈａｎｉｃｓ［Ｍ］．Ｘｉ　ａｎ：Ｘｉ　ａｎ　Ｊｉａｏｔｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ　Ｐｒｅｓｓ，１９９１．　尹奇志，肖金生，吕运冰，等．孔边应力集中和裂纹尖端应力强度因子的有限元分析［Ｊ］．武汉理工大学学报（交通科　学与工程版），２００２，２６（１）：４７—５０．　‘　ⅥＱｉ—ｚｈｉ，ＸＩＡＯ　Ｊｉｎ－ｓｈｅｎｇ，ＬＶ　Ｙｕｎ—ｂｉｎｇ，ｅｔ　ａ１．Ｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ａｎａｌｙｓｉｓ　ｏｆ　ｈｏｌｅ　ｅｄｇｅ　ｓｔｒｅｓｓ　ｃｏｎｃｅｎｔｒａｔｉｏｎ　ｎｄ　ａｃｒａｃｋ　ｔｉｐ　ｓｔｒｅｓｓ　ｉｎｔｅｎ－　ｓｉｔｙ　ｆａｃｔｏｒｓ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｗｕｈａｎ　Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ｔｅｃｈｎｏｌｏｇｙ（Ｔｒａｆｉｃ　Ｓｃｉｅｎｃｅ　ｎｄ　ａＥｎｇｉｎｅｅｒｉｎｇ　Ｅｄｉｉｔｏｎ），２Ｏ０２，２６（１）：４７—５０．　肖筱南，赵来军，丁先林，等．现代数值计算方法［Ｍ］．北京：北京大学出版社，２００２．　ＸＩＡＯ　Ｘｉａｏ—ｎａｎ，ＺＨＡＯ　Ｌａｉ－ｊｕｎ，ＤＩＮＧ　Ｘｉｎ—ａｌｉｎ，ｅｔ　ａ１．Ｍｏｄｅｍ　Ｎｕｍｅｒｉｃａｌ　Ｃａｌｃｕｌａｔｉｏｎ　Ｍｅｔｈｏｄ［Ｍ］．Ｂｅｒｉｎｇ：Ｂｅｒｉｎｇ　Ｕｎｉｖｅｒ－　ｓｉｔｙ　Ｐｒｅｓｓ，２００２．　刘怀恒．结构及弹性力学有限元法［Ｍ］．西安：西北工业大学出版社，２００７．　ＬＩＵ　Ｈｕａｉ—ｈｅｎｇ．Ｓｔｒｕｃｔｕｒｅ　ａｎｄ　Ｅｌａｓｔｉｃｉｙｔ　Ｆｉｉｎｔｅ　Ｅｌｅｍｅｎｔ　Ｍｅｔｈｏｄ［Ｍ］．Ｘｉ　ｎ：ａＮｏｒｔｈｗｅｓｔｅｍ　Ｐｏｌｙｔｅｃｈｎｉｃａｌ　Ｕｎｉｖｅｓｒｉｙｔ　Ｐｒｅｓｓ，２ＯＯ７．　尹双增．断裂・损伤理论及应用［Ｍ］．北京：清华大学出版社，１９９２．　ＹＩＮ　Ｓｈｕａｎｇ－ｚｅｎｇ．Ｆｒａｃｔｕｒｅ，Ｄａｍａｇｅ　Ｔｈｅｏｒｙ　ａｎｄ　Ｉｔｓ　Ａｐｐｌｉｃａｔｉｏｎ［Ｍ］．ＢｅＯｉｎｇ：Ｔｓｉｎｇｈｕａ　Ｕｎｉｖｅｒｓｉｔｙ　Ｐｒｅｓｓ，１９９２．　闫月梅，郭秉山．Ｙ型偏心支撑钢框架受力性能有限元分析［Ｊ］．西安科技大学学报，２００９，２９（２）：１５４—１５８．　ＹＡＮ　Ｙｕｅ—ｍｅｉ，ＧＵＯ　Ｂｉｎｇ—ｓｈａｎ．Ｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ａｎａｌｙｓｉｓ　ｏｆ　Ｙ—ｔｙｐｅ　ｅｃｃｅｎｔｒｉｃａｌｌｙ　ｂｒａｃｅｄ　ｓｔｅｅｌ　ｆｒａｍｅ　ｍｅｃｈａｎｉｃａｌ　ｐｒｏｐｅｒｔｉｅｓ［Ｊ］．　Ｊｏｕｎｍａｌ　ｏｆ　Ｘｉ　ａｎ　Ｕｎｉｖｅｓｉｒｔｙ　ｏｆ　Ｓｃｉｅｎｃｅ　ａｎｄ　Ｔｅｃｈｎｏｌｏｇｙ，２００９，２９（２）：１５４—１５８．　蔡来生，俞焕然．拉压模量不同弹性物质的本构［Ｊ］．西安科技大学学报，２００９，２９（１）：１７—２１．　ＣＡＩ　Ｌａｉ—ｓｈｅｎｇ．ＹＵ　Ｈｕａｎ—ｒａｎ．Ｃｏｎｓｔｉｔｕｔｉｖｅ　ｒｅｌａｔｉｏｎ　ｏｆ　ｅｌａｓｔｉｃ　ｍａｔｅｒｉａｌｓ　ｗｉｔ｝ｌ　ｄｉｆｅｒｅｎｔ　ｅｌａｓｔｉｃ　ｍｏｄｕｌｉ　ｉｎ　ｔｅｎｓｉｏｎ　ａｎｄ　ｃｏｍｐｒｅｓ－　ｓｉｏｎ［Ｊ］．Ｊｏｕｍａｌ　ｏｆＸｉ　ａｎ　Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ｓｃｉｅｎｃｅ　ａｎｄ　Ｔｅｃｈｎｏｌｏｇｙ，２００９，２９（１）：１７—２１．　Ｃｏｍｐｕｔａｔｉｏｎ　ｏｆ　ｃｒａｃｋ　ｓｔｒｅｓｓ　ｉｎｔｅｎｓｉｔｙ　ｆａｃｔｏｒ　ｂａｓｅｄ　ｏｎ　ＦＥＭ　ＱＩＡＯ　Ｂａｏ—ｍｉｎｇ　（Ｃｏｌｌｅｇｅ　ｏｆ　Ｓｃｉｅｎｃｅｓ，Ｘｉ　ａｎ　Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ｓｃｉｅｎｃｅ　ａｎｄ　Ｔｅｃｈｎｏｌｏｇｙ，Ｘｉ　ａｌｌ，７１００５４，Ｃｈｉｎａ）　Ａｂｓｔｒａｃｔ：Ｔｈｉｓ　ｐａｐｅｒ　ａｎａ１），ｚｅｓ　ｔｈｅ　ｓｈｉｔ，ｆｓｔｒｅｓｓ　ａｎｄ　ｓｔｒａｉｎ　ｆｉｅｌｄ　ｄｉｓｔｒｉｂｕｔｉｏｎ　ｏｆ　ｃｅｎｔｒａｌ—ｔｈｒｏｕｇｈ－ｃｒａｃｋ　ｔｉｐ　ｏｎ　ａ　ｉｆｎｉｔｅ　ｆｉａｔ　ｐｌａｔｅ　ｂａｓｅｄ　ｏｎ　ｔｈｅｏｒｅｔｉｃａｌ　ｍｅｔｈｏｄ　ａｎｄ　ｎｕｍｅｒｉｃａｌ　ｍｅｔｈｏｄ，ａｎｄ　ｃａｌｃｕｌａｔｅｓ　ｔｈｅ　ｓｔｒｅｓｓ　ｉｎｔｅｎｓｉｔｙ　ｆａｃｔｏｒ　ｏｆ　ｔｈｅ　ｃｒａｃｋ　ｔｉｐ．Ｂｙ　ｃｏｍｐａｒｉｎｇ　ｔｈｅ　ｎｕｍｅｒｉｃａｌ　ｒｅｓｌｔｕ　ｗｉｔｈ　ｔｈｅ　ｔｈｅｏｒｅｔｉｃａｌ　ｒｅｓｕｌｔ，ｉｔ　ｈａｓ　ｐｒｏｖｅｄ　ｔｈａｔ　ｔｉａｒｎｇｌｅ　ｅｌｅ－　ｍｅｎｔｓ　ａｎｄ　ＦＥＭ　ａｒｅ　ｃｏｎｓｉｄｅｒａｂｌｙ　ａｃｃｕｒａｔｅ　ａｎｄ　ｅａｓｉｅｒ　ｔｏ　ｕｓｅ　ｉｎ　ｔｈｅ　ｃａｌｃｕｌａｔｉｏｎ　ｏｆ　ｓｔｒｅｓｓ　ｉｎｔｅｎｓｉｔｙ　ｆａｃｔｏｒ．　Ｋｅｙ　ｗｏｒｄｓ：ｃｒａｃｋ；ｓｔｒｅｓｓ　ｉｎｔｅｎｓｉｔｙ　ｆａｃｔｏｒ；ｆｏｕｒｉｅｒ　ｔｒａｎｓｆｏｒｍ；ｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ｍｅｔｈｏｄ；ｌｅａｓｔ—ｓｑｕａｒｅ　ａｐ—　ｐｒｏｘｉｍａｔｉｏｎ　Ｃｏｒｒｅｓｐｏｎｄｉｎｇ　ａｕｔｈｏｒ：ＱＩＡＯ　Ｂｏａ—ｍｉｎｇ，Ａｓｓｏｃｉａｔｅ　Ｐｒｏｆｅｓｓｏｒ，Ｘｉ　ａｎ　７１（１０５５，Ｐ．Ｒ．Ｃｈｉｎａ，Ｔｅｌ：ＯＯ８６—１５９９１７０６０５８，Ｅ—ｍａｉｌ：ｑｉａｏｂｍ＠ｘｕｓｔ．ｅｄｕ．ｃｎ