应力腐蚀裂纹尖端的有限元数值分析与可视化
周晓敏 温庆阳
摘要： 针对在微观状态下结构力学行为会受尺度效应影响的问题，在偶应力理论中考虑微观结构的旋转梯度可以较好解释结构的尺度效应.建立基于一般偶应力理论的有限元法的基本方程，并在裂纹尖端引入奇异单元，计算受单向拉伸的中心斜裂纹板裂纹尖端场的应力强度因子（Stress Intensity Factor，SIF），分析特征长度变化对SIF的影响，对比偶应力理论下的结果与经典理论下的结果.结果表明：在裂纹尖端引入奇异单元可以提高计算精度和稳定性；偶应力使得裂纹尖端SIF比经典理论下的值小，并且SIF随着特征长度增大而减小.
关键词： 微观结构； 偶应力； 旋转梯度； 奇异单元； 应力强度因子； 有限元法； 尺度效应
中图分类号： TU452； O344.3文献标志码： A
Couple stress finite element method using
singular element at crack tip
ZHOU Xiaomin， WEN Qingyang
（School of Architectural and Surveying & Mapping Engineering， Jiangxi University of Science and Technology，
Ganzhou 361000， Jiangxi， China）
Abstract： As to the issue that the mechanical behavior of microstructure may be affected by the scale effect， the couple stress theory can explain the scale effect very well by considering the effect of rotating gradient. The basic equation of finite element method is built based on the general couple stress theory and the singular element is introduced into crack tip. The Stress Intensity Factor（SIF） of the crack tip field of a plate with incline center crack in uniaxial tension is calculated， and the effect of characteristic length on the SIF is analyzed； the results obtained by the couple stress theory and the classical theory are compared. The results show that， the calculation accuracy and stability can be improved by introducing singular element into crack tip， the SIF obtained by couple stress is less than that determined by the classical theory， and the SIF increases with the decrease of the value of characteristic length.
Key words： microstructure； couple stress； rotating gradient； singular element； stress intensity factor； finite element method； size effect
收稿日期： 2015[KG*9〗11[KG*9〗29修回日期： 2016[KG*9〗01[KG*9〗19
基金项目： 国家自然科学基金（41462009）；江西省自然科学基金（20151BAB206025）； 江西省教育厅科学技术研究项目（GJJ150628， GJJ150629）
作者简介： 周晓敏（1988—），女，江西吉水人，助教，硕士，研究方向为计算力学和结构强度分析，（Email） zhouxmim@126.com0引言
在微观状态下，结构的特征尺寸与材料的特征长度处于相近数量级，材料缺陷（微裂纹）处不同点之间的相互作用不可忽略，因此结构将出现传统连续介质理论无法解释的尺度效应现象.此时，建立基于连续介质框架、考虑尺度效应的本构模型就成为联系经典连续介质力学和微观力学之间的桥梁.偶应力理论就是在此背景下发展起来的应变梯度理论的一种.COSSERAT兄弟[1]提出一般偶应力理论，将材料的特征长度和第二切变模量引入结构的本构关系中，系统建立考虑偶应力影响时的本构模型，解释材料的尺度效应现象.基于偶应力理论的力学问题只有极少数问题有解析解，一般采用有限元法及扩展有限元法（Extended Finite Element Method， XFEM）研究.赵勇等[2]和齐磊等[3]采用有限元法研究考虑偶应力影响下的小孔应力集中问题.王胜军[4]基于有限元法给出一般偶应力理论与转动约束偶应力理论的异同.张敦福等[56]和吴延峰等[7]采用有限元法分析考虑偶应力理论情况下的岩石剪切带等问题.赵冰等[8]采用有限元法分析微梁弯曲的尺度效应.徐慧等[9]和侯日立等[10]采用有限元法研究单向拉伸裂纹板.茹忠亮等[11]给出扩展有限元法在断裂力学问题上的应用.吴圣川等[1213]对断裂力学问题及其CAE软件的现状和发展给出系统的综述，并且基于成熟的扩展有限元法和自主研发的虚节点法推出具有完全自主知识产权的三维疲劳裂纹扩展分析软件ALOF.霍金东等[14]采用Abaqus的XFEM功能对折弯片的断裂问题进行仿真.陈星文[15]运用XFEM分析三维高压管道中的裂纹扩展问题.
由于裂纹尖端应力场具有奇异性，计算结果受网格划分影响较大.BARSOUM[16]提出奇异单元，并应用于断裂力学中裂纹尖端应力场的求解，提高有限元法计算断裂问题的精度和稳定性，避免裂纹尖端区域网格划分困难.沈辉等[17]研究奇异单元有限元法计算精度受网格划分的影响情况.李尧臣等[18]提出引入位移协调奇异单元的有限元法.林广平等[19]构造三维奇异单元.段庆全等[20]采用奇异单元分析裂纹扩展问题.
本文研究裂纹尖端应力强度因子（Stress Intensity Factor，SIF）的尺度效应现象，将一般偶应力理论与有限元法相结合，并在裂纹尖端引入奇异单元，分别分析含边裂纹的三点受力梁和受单向拉伸的中心斜裂纹板的裂纹尖端SIF，给出偶应力和奇异单元对SIF的影响，并分析特征长度与裂纹长度比值变化对SIF的影响规律.
1分析方法
基于一般偶应力理论的平面应力单元（见图1），微元上作用的单位面积力和单位面积力偶的直角分量分别为（σx，σy，τxy，τyx）和（μx，μy）.
图 1考虑偶应力的微元体
Fig.1Microelement considering couple stress
分析考虑偶应力影响的平面弹性问题时，切应力不再具有互等性，需在每个平面微元中增加2个偶应力项，因此偶应力理论下的平衡方程需要在传统弹性力学的基础上进行相应修改，即为σxx+τyxy=0
τxyx+σyy=0
μxx+μyy+τxy-τyx=0 （1）剪切应变εyx和εxy示意见图2.在微元体中，由于偶应力的存在，微元体会有微转动出现，产生弯曲效应，并引发微转角ωz.
a） εyxb） εxy图 2剪切应变εyx和εxy
Fig.2Shear strains εyx and εxy
弯曲效应所产生的曲率使εxy≠εyx，剪应力互等定理不再成立.微元几何方程为εx
式中：E为弹性模量；v为泊松比；G为剪切模量；Gc为第二剪切模量；l为材料的特征长度，是为平衡一般偶应力理论中应变和应变梯度的量纲而引入的材料常数.l具有长度的量纲，是一种依赖于材料微观结构形态的特征常数，其大小取决于材料的微结构，不受外界约束或载荷影响.[5]
力边界条件为fx
fy
=[nm]σxτxyμx
τyxσyμy （4）式中：n和m为边界的外法线方向余弦；fx和fy为边界上已知的面力；为边界上已知的面力偶.
位移边界条件为u=
ωz=z （5）式中：和z分别为边界上已知的位移和转角.
为更好地描述裂纹尖端处应力场的奇异性，在裂纹尖端引入奇异单元.裂纹尖端划分的等参奇异三角形单元见图3，坐标原点设于裂纹尖端点，采用极坐标（r，θ）表示，将单元中与裂纹尖端相连边界上的节点移至边长的1/4处.a）奇异单元节点分布b）裂纹尖端坐标系及节点布置
三角形奇异单元的形函数可以由划线法得到.形函数为N1=（4ξ+4η-1）（ξ+η-1）
2数值算例
2.1含边裂纹的3点受力梁
含边裂纹3点受力梁模型见图4，其中：梁长L=8 m，梁高h=1 m，截面为正方形；在梁的中心沿高度方向有一长度为a=0.5 m的单边裂纹，在裂纹延长线与上表面的交点处作用有大小为F=10 kN的集中载荷.材料的弹性模量E=200 GPa，切变模量G=E/2（1+v），泊松比v=0.3，第二切变模量取Gc=G，特征长度与梁高比l/h=0.001.
图 4含边裂纹的3点受力梁模型
Fig.4Model of threepoints bending beam with edge crack
经典理论解、传统有限元解以及考虑奇异单元和偶应力下计算得到的含边裂纹三点受力梁模型SIF的结果比较见表1.由此可以看出：相对于未考虑奇异单元法的结果，采用奇异单元法得到的SIF更接近理论解；考虑一般偶应力理论后，裂纹尖端SIF要小于未考虑偶应力的值.
表 1不同方法得到应力强度因子
Tab.1Results of SIF obtained by different methods
kPa/m2求解方法SIF经典理论[21]213.0传统有限元191.8奇异单元211.8偶应力有限元177.5偶应力奇异单元201.3
2.2单向受拉中心斜裂纹板
两端受大小为σ=10 kPa均布力拉伸的中心斜裂纹板见图5，其中：板宽w=2 m，高h=2.5 m，裂纹长2a=0.2 m，与水平面夹角α=45°.图 5单向受拉中心斜裂纹板模型
Fig.5Model of plate with incline center crack in
uniaxial tension
裂纹尖端附近单元有限元网格划分见图6.材料的弹性模量E=200 GPa， 泊松比v=0.3，取特征长度与裂纹半长的比值l/a=1，第二切变模量与切变模量比值Gc/G=0.5[5]；KI=KI/K0，KII=KII/K0分别为无量纲化I型和II型SIF，其中K0=σπa.
图 6裂纹尖端附近单元有限元网格划分
Fig.6Finite element mesh around crack tip
在不考虑偶应力影响的条件下，采用解析法得到的经典理论解见表2.整体采用8节点四边形单元划分得到的8节点等参元解，整体采用8节点四边形单元划分后在裂纹尖端采用奇异单元划分得到奇异单元解.由表2可见：在相同网格划分密度下，采用奇异单元求得的复合型裂纹的I型和II型SIF相对于仅采用8节点等参元得到的结果具有更高的精度.
表 2I型和II型SIF
Tab.2SIF of type I and II求解方法KIKII经典理论[21]0.5000.5008节点等参元0.4140.416奇异单元0.4830.500
考虑偶应力影响的情况下不同特征长度对I型和II型SIF的影响见表3.由此可知：在不同的特征长度下，考虑偶应力影响计算得到的I型和II型SIF会随着特征长度与裂纹长度之比的增加而减小；当特征长度与结构尺寸处于相近数量级时，尺度效应明显；当特征长度远小于结构尺寸时，偶应力对计算结果的影响较小.
3结论
本文将奇异单元、偶应力理论和有限元法相结合，给出在裂纹尖端使用奇异单元的一般偶应力理论有限元法，得出如下主要结论.
（1）在同等网格划分密度下，采用奇异单元法得到的裂纹尖端SIF相对于仅采用8节点四边形单元得到的结果更接近理论解.
（2）考虑一般偶应力理论后会减小裂纹尖端SIF的值；随着特征长度的减小，SIF减小；在特征长度与裂纹半宽处于相近数量级时，SIF受特征长度影响明显；当特征长度远大于裂纹半宽时，SIF受特征长度改变影响较小.参考文献：
[1]COSSERAT E， COSSERAT F. Théorie des corps déformables[M]. Paris： Cornell University Library， 1909： 242246.
[2]赵勇， 张若京. 基于Cosserat理论的小孔应力集中问题的有限元分析[J]. 力学季刊， 2009， 30（3）： 410414.
ZHAO Y， ZHANG R J. Finite element analysis of stress concentration problem based on Cosserat theory[J]. Chinese Quarterly of Mechanics， 2009， 30（3）： 410414.
[3]齐磊， 张若京. 基于Cosserat理论的平面等参元[J]. 计算机辅助工程， 2010， 19（1）： 1216. DOI： 10.3969/j.issn.10060871.2010.01.005.
QI L， ZHANG R J. Plane isoparametric element based on Cosserat theory[J]. Computer Aided Engineering， 2010， 19（1）： 1216. DOI： 10.3969/j.issn.10060871.2010.01.005.
[4]王胜军. 关于弹性偶应力C0和C1理论及有限元[D]. 大连： 大连理工大学， 2008.
[5]张敦福， 王相玉， 李术才. 偶应力对含充填层状岩体边界层效应的影响[J]. 岩石力学与工程学报， 2011， 30（S1）： 32873294.
ZHANG D F， WANG X Y， LI S C. Interfaces boundary layer effect analysis of jointed rock mass with inclusion based on couple stress theory[J]. Chinese Journal of Rock Mechanics and Engineering， 2011， 30（S1）： 32873294.
[6]张敦福， 王相玉， 朱家明， 等. 偶应力对层状岩体结构面边界层效应的影响[J]. 岩土力学， 2012， 33（7）： 21812188. DOI： 10.3969/j.issn.10007598.2012.07.039.
ZHANG D F， WANG X Y， ZHU J M， et al. Influence of couple stress on interfaces boundary layer effect of layered rock mass[J]. Rock and Soil Mechanics， 2012， 33（7）： 21812188. DOI： 10.3969/j.issn.10007598.2012.07.039.
[7]吴延峰， 张敦福， 张波， 等. 拉力型锚杆锚固界面剪应力分布偶应力理论分析[J]. 岩土力学， 2013， 34（S1）： 187192.
WU Y F， ZHANG D F， ZHANG B， et al. Analysis of shear stress near anchorage interface of tensiletype anchor with couple stress theory[J]. Rock and Soil Mechanics， 2013， 34（S1）： 187192.
[8]赵冰， 郑颖人， 陈淳. Cosserat理论的有限元实现及微梁弯曲的尺寸效应模拟[J]. 长沙理工大学学报（自然科学版）， 2010， 7（4）： 3742. DOI： 10.3969/j.issn.16729331.2010.04.007.
ZHAO B， ZHENG Y R， CHEN C. Finite element method implementation of Cosserat theory and scale effect simulation of microbeam[J]. Journal of Changsha University of Science and Technology （Natural Science）， 2010， 7（4）： 3742. DOI： 10.3969/j.issn.16729331.2010.04.007.
[9]徐慧， 伍晓赞， 程仕平， 等. 复合裂纹的应力强度因子有限元分析[J]. 中南大学学报（自然科学版）， 2007， 38（1）： 7983. DOI： 10.3969/j.issn.16727207.2007.01.016.
XU H， WU X Z， CHENG S P， et al. Finite element analysis of stress intensity factor in composite mode crack[J]. Journal of Center South University （Science and Technology）， 2007， 38（1）： 7983. DOI： 10.3969/j.issn.16727207.2007.01.016.
[10]侯日立， 赵阳， 赵东， 等. 带裂纹板单向拉伸极限承载力有限元研究[J]. 固体力学学报， 2010， 31（S1）： 270274.
HOU R L， ZHAO Y， ZHAO D， et al. The finite element study of ultimate bearing capacity of centercracked plate[J]. Chinese Journal of Solid Mechanics， 2010， 31（S1）： 270274.
[11]茹忠亮， 朱传锐， 张友良， 等. 断裂问题的扩展有限元法研究[J]. 岩土力学， 2011， 32（7）： 21712176. DOI： 10.3969/j.issn.10007598.2011.07.042.
RU Z L， ZHU C R， ZHANG Y L， et al. Study of fracture problem with extend finite element method[J]. Rock and Soil Mechanics， 2011， 32（7）： 21712176. DOI： 10.3969/j.issn.10007598.2011.07.042.
[12]吴圣川， 吴玉程. 断裂分析及CAE软件的现状与发展[J]. 计算机辅助工程， 2011， 20（1）： 12. DOI： 10.3969/j.issn.10060871.2011.01.001.
WU S C， WU Y C. Present situation and development of fracture analysis and CAE software[J]. Computer Aided Engineering， 2011， 20（1）： 12. DOI： 10.3969/j.issn.10060871.2011.01.001.
[13]吴圣川， 吴玉程. ALOF： 新一代三维疲劳裂纹扩展分析软件[J]. 计算机辅助工程， 2011， 20（1）： 136140. DOI： 10.3969/j.issn.10060871.2011.01.031.
WU S C， WU Y C. ALOF： new 3D fatigue crack propagation analysis software[J]. Computer Aided Engineering， 2011， 20（1）： 136140. DOI： 10.3969/j.issn.10060871.2011.01.031.
[14]霍金东， 王涛， 柳占立， 等. 基于XFEM的折弯片断裂仿真[J]. 计算机辅助工程， 2013， 22（S2）： 1114. DOI： 10.3969/j.issn.10060871.2013.z2.003.
HUO J D， WANG T， LIU Z L， et al. Fracture simulation of bending sheet based on XFEM[J]. Computer Aided Engineering， 2013， 22（S2）： 1114. DOI： 10.3969/j.issn.10060871.2013.z2.003.
[15]陈星文. XFEM在高压管道裂纹扩展计算中的应用[J]. 计算机辅助工程， 2015， 24（4）： 5660. DOI： 10.13340/j.cae.2015.04.012.
CHEN X W. Application of XFEM in highpressure pipe crack propagation calculation[J]. Computer Aided Engineering， 2015， 24（4）： 5660. DOI： 10.13340/j.cae.2015.04.012.
[16]BARSOUM R S. Triangular quarterpoint elements as elastic and perfectlyplastic crack tip elements[J]. International Journal for Numerical Methods in Engineering， 1977， 11（1）： 8598. DOI： 10.1002/nme.1620110109.
[17]沈辉， 周储伟. 平面六节点等参应力奇异单元的精度研究[J]. 机械强度， 2012， 34（3）： 410414.
SHEN H. ZHOU C W. Study on the accuracy of plane sixnode isoparametric stress singular element[J]. Journal of Mechanical Strength， 2012， 34（3）： 410414.
[18]李尧臣， 王自强. 在裂纹尖端引入位移协调奇异单元的有限单元法[J]. 上海力学， 1983， 4（3）： 7999. DOI： 10.15959/j.cnki.0254 0053.1983.03.009.
LI Y C， WANG Z Q. The finite element method by using the singular element with concordant displacement at the crack tip[J]. Shanghai Mechanics， 1983， 4（3）： 7999. DOI： 10.15959/j.cnki.0254 0053.1983.03.009.
[19]林广平， 梅甫良. 12结点三维等参奇异单元的构造和应用[J]. 计算力学学报， 2009， 26（2）： 282286.
LING G P， MEI F L. Creation and application of a new type of threedimensional isoparametric singular element with 12 nodes[J]. Chinese Journal of Computational Mechanics， 2009， 26（2）： 282286.
[20]段庆全， 于国鹏， 张宏. 管道疲劳裂纹扩展数值模拟[J]. 石油化工设备， 2009， 38（2）： 2327. DOI： 10.3969/j.issn.10007466.2009.02.008.
DUAN Q Q， YU G P， ZHANG H. Numerical simulation of pipe fatigue crack propagation[J]. PetroChemical Equipment， 2009， 38（2）： 2327. DOI： 10.3969/j.issn.10007466.2009.02.008.
[21]中国航空研究院. 应力强度因子手册[M]. 北京： 科学出版社， 1981： 188189.
（编辑武晓英）
