有限元是目前解决工程实际问题最有效的数值方法，但其存在某些固有的缺陷^[1-2]：①求解裂纹类强间断问题需细分网格；②模拟大变形问题时网格需不断地重构；③处理夹杂问题时需沿夹杂和基体的界面处划分网格；④刚度矩阵过刚，位移解偏小等.为克服前三点困难提出了扩展有限元，为改进解的精度提出了光滑有限元.

XFEM由Belytschko等^[3]提出，是目前求解含断裂问题最有效的数值方法.XFEM基于单位分解法，在位移场中引入扩展项，其计算网格独立于结构的任何内部细节点，具有计算精度高、网格划分简单等特点.Moës，Sukumar等^[4-5]将该方法推广到了三维，Asadpoure等^[6]利用该方法研究了正交材料中的静态裂纹问题.Esnaashari等^[7]提出了求解裂纹问题的BXFEM.Motamedi等^[8]在动态裂纹扩展方面进行了研究.方修君等^[9]将XFEM嵌套于ABAQUS软件中，对含裂纹混凝土结构进行了研究；余天堂^[10]将XFEM与线性互补法相结合，求解了裂纹面非线性接触问题.

SFEM(smoothed finite element method)是Liu等^[11]将光滑应变措施引入有限元法，改进有限元法刚度结构的一种方法，具有形函数简单、对网格要求低、计算效率高等优点，现已广泛应用于各个领域^[12-15].

本文基于CSFEM，结合XFEM，提出了CS-XFEM(cell-based smoothed element method)，对含中心裂纹、斜裂纹的正交各向材料板进行了模拟，并与FEM，XFEM和BXFEM计算结果进行了对比.

1 复合材料断裂力学

正交各向异性材料的正轴应变与应力关系为

(1)

式中：和分别为应力和应变列阵；S为材料的柔度矩阵，二维空间中为

(2)

式中E, ν和G分别为弹性模量、泊松比和剪切弹性模量.

如图 1所示，考虑一个等厚度、均匀的正交各向异性体含一条穿透裂纹的情况，满足力边界和位移边界 ，(x, y)为全局坐标，(x′, y′)为局部坐标，(r, θ)为极坐标，假定弹性主方向与参考坐标轴一致时，平面应力状态下应力函数F应满足的变形协调方程为

(3)

特征方程的根为λ₁, ₁, λ₂, ₂. ₁, ₂分别为λ₁和λ₂的共轭复数，则裂纹尖端应力场和位移场的渐近解^[6]如下：

I型：

(4)

(5)

(6)

(7)

(8)

II型:

(9)

(10)

(11)

(12)

(13)

式中：Re表示取实部；K_I和K_II分别为I型和II型裂纹的应力强度因子；

(14)

(15)

2 Cell-based光滑扩展有限元法

Cell-based光滑扩展有限元法的位移模式与扩展有限元表达形式一致，即

(16)

式中：I为节点(图 2中‘•’)，J为被裂纹完全贯穿单元的节点(图 2中‘□’)，K为裂尖单元的节点(图 2中‘○’)；N_I^u(x), N_J^a(x)和N_K^b(x)分别为相应节点的形函数，u_I, a_J和b_K分别为相应节点的位移；N^CS-FEM, N^CS-c, N^CS-f分别为节点I, J, K的集合.H(x)为Heaviside函数：

(17)

式中：x^*为裂纹面节点坐标；n为外法向向量.

F_l(x)为裂尖处扩展函数：

如图 2所示，将求解域Ω离散为N_e个四边形单元，节点个数为N_d, Ω=∪_i=1^N_eΩ_i^e, Ω_i^e∩Ω_j^e=, i≠j, 为空集，再将Ω^e_i划分为n_c个光滑区域，共N_s个光滑子域.

应变满足：

(19)

式中 _I^u(x_k), _J^a(x_k)和 _K^b(x_k)分别为相应I, J, K节点的光滑应变矩阵，可统一表示为

(20)

式中：

(21)

(22)

(23)

式中：h=x, y；l=1, 2, 3, 4；N_seg为边界Γ^s_k的个数；N_gau为每段边界高斯点的个数；w_{m, n}为高斯权函数；n_x和n_y为积分段外法向向量的分量；x_{m, n}为第m段边界处的第n个高斯点；A^s_k为第k光滑区域的面积：

(24)

离散方程为

(25)

式中：

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

式中：

(35)

式中：C为弹性矩阵；为体力；为面力；N为有限元形函数.

3 交互积分

考虑两种独立的平衡状态：状态1(σ_ij⁽¹⁾, ε_ij⁽¹⁾, u_i⁽¹⁾)为真实物理场状态，状态2(σ_ij⁽²⁾, ε_ij⁽²⁾, u_i⁽²⁾)为辅助物理场状态.叠加状态1和状态2可得到另一状态的J积分^[16]

(36)

式中：δ_1j为克罗内克函数；A为求解域；q为任一可微函数.

整理式(36)，得

(37)

式中：M⁽¹⁺²⁾为交互积分，

(38)

(39)

式(38)可化为

(40)

式中：K_I⁽¹⁾和K_II⁽¹⁾为真实场下的I型和II型应力强度; K_I⁽²⁾和K_II⁽²⁾为辅助场下的I型和II型应力强度;

(41)

(42)

(43)

式中Im表示取虚部.

取K_I⁽²⁾=1, K_II⁽²⁾=0，式(38)为

(44)

取K_I⁽²⁾=0, K_II⁽²⁾=1，式(38)为

(45)

4 数值算例 4.1 算例1

含中心裂纹的正交各向异性材料板受均布载荷作用，裂纹长度为2a，单位板厚、几何构型、加载方式，以及网格划分为4 900时的情况如图 3所示.材料参数：E₁₁=114.8 GPa，E₂₂=11.7GPa，G₁₂=9.66 GPa，ν₁₂=0.21.

表 1给出了FEM，XFEM，CS-XFEM和BXFEM求解含中心裂纹复合材料板应力强度因子K_I的结果，其中n_cell为单元数.从表中可以看出CS-XFEM具有较高的计算精度，与BXFEM和XFEM所得结果十分接近，远高于FEM求解精度；也可看出，积分区域c的选取对计算结果影响不大.CS-XFEM不仅具有XFEM的优点：单元与裂纹面相互独立, 裂尖不必为单元节点，裂尖处也不需要网格加密，还具有CSFEM形函数简单、对网格要求低的特点.

图 4给出了CS-XFEM得到的应力云图，很明显地表现出了应力场的不连续性和正交特性效应，从而也说明了CS-XFEM的正确性.

4.2 算例2

含斜裂纹的正交各向异性材料板受均布载荷作用，裂纹长度为a, φ=45°，单位板厚、几何构型、加载方式和单元划分如图 5所示.材料参数：E₁₁=0.81GPa，E₂₂=11.84GPa，G₁₂=0.63GPa，ν₁₂=0.38.

表 2给出了CS-XFEM和XFEM求解应力强度因子的结果，可见两者精度基本一致，证明了CS-XFEM的正确性与有效性.图 6给出了c/a=0.3，0.4，0.5，0.6，0.7，0.8时，采用CS-XFEM计算所得的应力强度因子K_I和K_II，可见CS-XFEM对c/a不敏感，具有较高的求解精度.

5 结论

1) CS-XFEM的计算精度同XFEM和BXFEM精度基本相同，远高于FEM求解精度.

2) CS-XFEM兼具CSFEM和XFEM的优点.

3) CS-XFEM对c/a的取值不敏感.