基于拓展有限元的应力强度因子方法研究
【摘要】本文应用拓展有限元方法，首先研究不同参数对于结果精确性的影响。并编写了基于拓展有限元的Matlab程序,改进了判定应力强度因子的方法,并与Abaqus得出的结果进行比较，结果发现，比用Abaqus判定得到的结果更为准确。
【关键词】拓展有限元，Matlab，Abaqus，应力强度因子，SIF 0 前言
在工程中的许多问题都能采用基本的物理方法来解决。这些基本的物理方法加上偏微分方程，或者是我们平时所说的控制方程。然而在大多数情况下使用控制方程并不能得到精确地解，主要还是因为处理现实生活中的问题面临着更复杂的情况和边界条件。普通的基于h-p方法[1]的有限元方法在处理连续场问题中可以表现的非常优异，然而对于不连续场问题为了获得精确地解，使用多项式近似法的有限元方法会需要花费更多的精力在网格的重划分上，这就是说每步的裂缝增长都需要依靠有限元软件的对于裂缝开展间断的重新划分。
为了突破有限元对于不连续场问题的限制，由Belytschko 和Black[2]共同提出了拓展有限元方法。这种方法克服了有限元方法在解决不连续场问题中重划分问题。拓展有限元方法能很好的应用在断裂力学中的不连续问题，如裂缝扩展，存在于复合材料截面裂缝问题。如今拓展有限元已经能应用在工程中的各个领域，然而为了能模拟更多的问题，拓展有限元还需要更多的进步，如多尺度的模拟。
本文主要目的是使用现在流行的商用有限元软件Abaqus的拓展有限元模块进行裂缝模拟，提出改进的裂缝扩展的计算公式，编写出Matlab程序然后进行比较。通过两者计算结果的比较得出程序的准确性。
1拓展有限元基本原理
一般来说拓展有限元大部分是基于有限元的方法，所以应用在有限元的方程还是可以使用的，然而对于不连续的部分，拓展有限元提出了新的概念-子区间。简单来说，子区间就是使用一组新的目标函数对位移函数来内插值。最终的位移函数如下式
(1)
其中N 是节点的个数，Nl是标准的形函数,H(x)是富集函数，也可称为跳远函数，它遵循着单位分解原则。 和 分别是标准状况下和富集状态下的节点自由度。而4代表了裂尖富集公式的个数。公式(1)的内插可以用在标准的伽辽金（Galerkin） 方法。
跳跃方程H(x),使用在不连续的部分，也就是在裂缝的处（而非裂尖）。跳跃方程是按照裂缝发展的方向建立坐标系，并且规定裂缝一侧为单位正值，另一侧为单位负值。如下列公式：
作为另一个加入的富集部分，用来描述裂缝在裂尖的部分，可以用公式（2）来表示。（2）
其中和r 是以裂尖为坐标原点而得到的，r 为单元端点到裂缝尖端的距离，为裂缝和坐标系的夹角，正如图1。
图1：基于XFEM 的裂缝模拟
2.Matlab程序的应用思路:
Matlab 应用的方法是基于Rice 提出的J-积分方法，这个方法起源于二维方法的线积分。在裂缝的两个面上有一个封闭的围线，对于这个围线进行积分，Rice发现这个为围线积分是恰好是裂缝开展的能量释放值。如公式（3）
（3）
其中Ｕ应变能密度，ｔ是拉力向量，其中定义为向外的拉力为正向量，ｄ是位移向量，ｄｓ是沿着积分的单元长度。其中作为线性分析下，J-积分能直接用来应力强度因子如公式（4）
（4）
我们注意到G为应变能释放率，我们可以通过文献（3）提出等效区域积分方法计算求得，因此就可以的到应力强度因子的值。
3.有限元模拟
本文算例采用Moes(1999)[4]中应用的例子，算例的尺寸如图2.其中长宽均为l=0.4M,半长裂缝与边长的比值不断改变（由0.1到0.9）。其中例子的计算参数如表格1。
表1.材料特性
图2 问题1的边界条件
这个例子中将会给出在不同的角度和裂缝长度下求得不同的应力强度因子。这个带有裂缝平板将受到一对拉应力。这个裂缝居于板中间有着角度。所有的模拟方法都是基于拓展有限元的方法。精确地应力强度因子由，Sih (et.al.1973)[5]提出的计算公式得到。公式如（5）（6）：
（5）
（6）
这个算例中，首先取不同的裂缝长度，然后选出误差率最小裂缝长度，然后在此长度下选取不同的角度（150,300,450,600），并且通过提取计算结果中的应力强度因子来比较两种方法的准确性。
表2 不同裂缝长度下的应力强度因子
图3 两种方法下的误差率
这表2和图3 展示了Abaqus和Matlab 下在准静态问题上的应力强度因子的结果。
我们从表2 和图3、4 可以得到如下结论：
通过改进的Matlab程序获得的应力强度因子能更好的吻合有Sic提出的理论解，有更高的精度。
就裂缝长度来看，在裂缝长度相对整个板长的比例较小的情况下，应力强度因子会更接近Sic提出的解析解。在模拟此类问题时，应让裂缝保持与板长的比例为0.1，这样得到的应力强度因子是最准确的。
随着裂缝角度的不断变化，可以看出在Matlab模拟下的应力强度因子相对于解析解依旧是十分准确。