四节点等参单元xfem程序设计

[ ]
b j1i b 2i
j
，j = 1 ～ 4
（ 2）
k11 au k11 kbju 11 uu k41 au k41 au k41
uu
ua k11 aa k11 b ja k11
ub j k11 ab j k11 b jb j k11
… … … …
XFEM Programming for Four Nodes Isoparametric Element and Its Application in Crack Problems
Su Yi，Wang Shengnan，Yan Xiaozhong
（ School of Aeronautics，Northwestern Polytechnical University，Xi'an 710072 ）
2012 年 第 31 卷
12 月 第 12 期
机械科学与技术 Mechanical Science and Technology for Aerospace Engineering
December Vol． 31
2012 No． 12
四结点等参元 XFEM 程序设计及 在裂纹问题中的应用
苏
苏 毅
毅， 王生楠， 闫晓中
x =
Ni xi ， y ∑ i =1
=
Ni yi ∑ i =1
（ 1）
XFEM 中四结点等参单元的位 关于裂纹问题， 移模式为
[u ] = ∑ N （ x） [ v ] + ∑ α N （ x） H [ a v
4
ui
i
4
a1i
2i
i
i
i
i =1 4
i =1
]
+
∑ β i N i （ x） φj
i =1
T
为结点与 Heaviside 函数有关的加强自由度 ； 第三 项只有在点 x 处在裂尖单元 、 单元结点需要改进时 当第三项有意义时 β i = 才有意义 ； β i 为指示函数 ， 1， ； （ x ） x 否则为零 φ j 为点 的改进函数 ｛ φ（ r ， r sin θ ， r cos θ ， r sin θ sin θ ， θ） ｝ = ｛ 槡 槡 槡 2 2 2 r cos θ sin θ ｝ ( b j1 i 槡 2
1 2 3 4 Ω T 48 ×3
［ D］ B1 ， B2 ， B3 ， B4］ 3 ×3［ 3 ×48 tdxdy =
uu k14 au k14 ua k14 aa k14 b ja k14 ub j k14 ab j k14 b jb j k14 j =1 ～4 ， ub j k44 ab j k44 b jb j k44
苏
毅等： 四结点等参元 XFEM 程序设计及在裂纹问题中的应用
1951
单元荷载列阵 f u a 1 2 3 fi e = ｛ fi ， fi ， f ib， f ib， f ib， f ib4 ｝ T
扩展有限元以常规有限元和单位分解法为基 础， 保留了 传 统 有 限 元 的 优 点 。 应 用 扩 展 有 限 元 分析断裂问题时 ， 先不考虑裂纹面的位置 ， 直接划 分网格 ， 然后应用单位分解法思想 ， 在位移表达式 在裂纹贯穿区 中增加了附 加 函 数 提 高 了 自 由 度 ， ［10］ 域的单元节点用广义 Heaviside 函数 加强 ， 来反 映裂纹面的位移不连续性 ， 对包含裂尖的单元节点 用裂尖渐进位移场函数来加强， 来反映裂尖区域的 局部特性。这样就可以用附加函数来间接反映不连 从而在划分网格时候， 裂纹这种强不连续 续的存在， 的存在与有限元网格相互独立， 因此克服了在裂纹 尖端等高应力区和变形集中区进行高密度网格剖分 带来的困难。 1 ） 位移场表达式 对于四结点等参单元， 其坐标变换式同标准有 限元的位移模式， 即
收稿日期： 2011 － 11 － 04 作者简介： 苏 毅（ 1983 － ） ， 博士研究生， 研究方向为飞机结构疲劳 suxiaoyi12@ 126． com； 王生楠（ 联 断裂可靠性及损伤容限， wangshna@ nwpu． edu． cn 系人） ， 教授， 博士生导师，
［1 ］ 不连 续 问 题 。 XFEM 是 基 于 单 位 分 解 的 方 法 （ PUM） 对单元的形函数加以改进， 从而考虑所研究
问题的不连续、 奇异性等特性。 XFEM 所使用的网 格与结构内部几何或物理界面无关， 从而克服了裂 纹尖端等高应力和变形集中区内网格划分的困难 ， 使得在模拟裂纹生长过程中无需对网格进行重新划
1950
机械科学与技术
第 31 卷
分。自 XFEM 问世以来， 在国际上得到了很快的发 ［2 ～ 6 ］ 。 文献［ 7］ 展和广泛的应用 综述了 XFEM 在静 态 和 扩 展 裂 纹 问 题 中 的 应 用， 并与广义有限元 （ GFEM） 进行了比较， 8］ 结果符合的很好。文献［ 采 用 XFEM 研究了双材料界面裂纹问题的应力强度 9］综述了 XFEM 的基本思想、 因子的计算。 文献［ 实施步骤及其应用， 初步展望了该领域发展需要研 究的课题。 尽管 XFEM 在处理裂纹等这种强不连续性问 题时理论上是成功的， 但在方法的实现和应用上， 尚 存在许多技术如网格密度、 结点加强区域范围等问 题值得研究。为了验证扩展有限元法在计算裂纹应 力强度因子上的有效性， 在已有 XFEM 理论的基础 上， 推导了四结点等参元相应的计算公式 ， 编写了完 对 XFEM 整的 Matlab 程序。 针对典型含裂纹平板， 的应用进行了研究， 并进一步研究了其计算结果对 网格密度的敏感性。 1 XFEM 中四结点等参元的计算 SIF
b ju … … k14
k
ua 41 aa k41 b ja k41
k
ub j 41 ab j k41 b jb j k41
… … … …
k
uu 44 au k44

ua k44 aa k44 b ja k44
b ju … … k44
（ 5）
第 12 期
( )
( )
( )
( )
b j2 i ) 为改进结点两个方向的
T
只有第 附加自由度。当结点为常规单元的结点时， 一项； 当结点为被裂纹贯穿单元的结点时 ， 有第一项 和第二项； 当结点为裂纹尖端所在单元的结点时 ， 有 第一项和第三项； 结点同时处于裂纹贯穿单元和裂 纹尖端所在单元， 应优先属于裂纹尖端所在单元， 加 强方式选择裂纹尖端所在单元结点加强方式 。 2 ） 虚功方程和支配方程
4 4
∫ σ·δεdΩ = ∫ t·δudΓ + ∫ b·δudΩ
Ω Γt Ω
（ 3）
b 分别为分布体力、 式中： t， 分布面力。 位移模式构 造后， 和常规有限元方法一样， 将式 （ 2 ） 代入式 （ 3 ） 由虚功原理推导扩展有限元的支配方程 Kδ = F （ 4）
式中： δ 是结点未知向量； K 和 F 分别为总体刚度矩 K 和 F 先逐个单元计算， 阵和总体荷载列阵， 再按常 单元层次上的 k 和 f 分别为 K I 单元 规步骤进行组装。 劲度矩阵 k ［ k ij］ = B， B， B］ ［B ，
v i 为点 x 所 式中 ： N i （ x ） 为标准有限元形函数 ； u i 、 4 （ u 、 v 在单元 个结点位移分量 i i 并非最终的结点 单 位移值 ） ； 第二项只有在点 x 处在裂纹贯穿单元 、 元结点需要改进时才有意义 ； α i 为指示函数 ， 当第 二项有意义时 α i = 1 ， 否则为零 ； H 为跳跃函Байду номын сангаас ， 在 局部坐标系下定义为 H = +1 {－ 1 y ＞ 0 y ＜ 0 ； （ a1i ， a2i ）
Abstract： It was successful that extended finite element method （ XFEM ） was theoretically applied to the strong discontinuity problems such as cracks so that the development and application of the method are being increased rapidly． However，in practical applications，there are still many technical issues such as the mesh density to be studied． In order to verify and improve the effectiveness of XFEM on the calculation of crack stress intensity factor， the corresponding formula of XFEM with fournode isoparametric element was derived and a complete Matlab code was also edited aiming at calculating the crack tip stress intensity factor of a plate with crack． Based upon the program，the crack tip stress intensity factors for a typical plate with crack were calculated and the results by XFEM were compared with those by the traditional finite element method （ FEM ） ． The effect of the mesh parameters on the XFEM results was further studied． The study shows that XFEM has a very good accuracy in the calculation of crack tip stress intensity factor，but the results by XFEM are sensitive to the density of the mesh，this should be paid more attention to practical applications． Key words： XFEM； crack tip； stress intensity factor； four nodes isoparametric element； mesh density； discontinuity problem； cracks 1999 年， 以美国西北大学 Belytschko 教授为代 表的研究组首先提出用扩展有限元 （ XFEM ） 来解决
（ 西北工业大学 航空学院， 西安 710072 ）
要： 虽然扩展有限元法( XFEM) 在处理裂纹等这种强不连续性 问题 时 理 论 上 是 成 功 的， 得到很 快的发展和广泛的应用， 但在实际应用中， 尚存在许多技术问题如网格密度等值得研究。为了验证 摘 和提高 XFEM 在 计算 裂 纹 应力 强 度因子 上 的 有 效性， 针 对 四 结 点 等 参 元 推 导 了 XFEM 的 相 应 公 式， 编写了用于计算含裂 纹板 裂 纹 尖端 应力 强 度因子的 完 整 的 Matlab 程 序。 针 对 典 型 含裂 纹 平 板， 采用本文编写的程序计算了裂纹尖端应力强度因子， 与采用传统有限元法的结果进行了对比分 XFEM 在 计算 裂 纹 尖端 应力 强 度 析， 并进一步研究了网格参数， 对 XFEM 结果 的 影响。 结果 表明， 因子上有很好的计算精度， 但其计算结果对网格密度较为敏感， 在实际应用中应当引起重视。 关 键 词： 扩展有限元; 裂纹尖端; 应力强度因子; 四结点等参元; 网格密度 中图分类号： O346. 1 文献标识码： A 8728 （ 2012 ） 12194906 文章编号： 1003-