功能梯度材料^[1](functionally graded material，FGM)是一种材料属性呈连续变化的新型功能材料，因其具有缓解热应力、释放残余应力及提高黏结的能力，被广泛应用于航空航天、能源、建筑与生物医学等领域.受生产过程及工作环境等因素影响，FGM与均匀材料一样不可避免地产生裂纹，且其材料属性的非均匀性导致FGM结构中裂纹的断裂行为更为复杂，通常呈现为混合形式.在外荷载作用下，这些裂纹易进一步扩展，严重时会导致FGM结构的整体破坏.因此，为提高带裂纹FGM结构的安全性，需深入研究裂尖断裂力学参数的变化规律，而应力强度因子(stress intensity factors，SIFs)作为结构断裂分析的重要参数之一，表征了裂尖应力-应变场的强弱程度，对深入研究带裂纹FGM结构的破坏机理及提高结构的安全性具有重要意义.

近年来，众多学者围绕带裂纹FGM结构中裂尖SIFs问题进行了大量研究.Guo等^[2-3]建立了分段指数模型(PE模型)，通过求解奇异积分方程组，研究了功能梯度分层结构中穿越界面任意方向的裂纹，发现其夹角对混合型裂纹尖端SIFs的影响较大.Pan等^[4]在Guo等^[2-3]建立的PE模型基础上，同样采用求解奇异积分方程组研究了含共线裂纹的任意热机械属性功能梯度条I型热应力，分析了非均匀常数和裂纹几何参数对裂尖热应力强度因子(TSIF)的影响.Wang等^[5]基于开尔文基本解建立了FGM结构断裂分析的边界-域积分方程，采用径向积分边界元法分析含圆盘状裂纹的FGM三维体中裂尖SIFs，研究表明材料梯度平行于裂纹面方向变化对I型和II型SIFs均有显著影响.Zhang等^[6]应用径向积分边界元法对三维连续非均匀各向同性线弹性FGM结构进行了断裂分析，研究了材料梯度对裂纹张开位移和SIFs的影响.前述文献采用的奇异积分方程组和径向积分边界元法均是半解析半数值方法，需经大量繁琐的理论推导过程，难以适用于非规则边界和复杂荷载条件的结构.随着计算机运算能力的提升，数值方法得到了迅速发展，Abdollahifar等^[7]首次将无网格局部彼得洛夫伽辽金法用于分析带边界裂纹功能梯度条的I型裂纹裂尖SIFs；随后，又采用此方法对带裂纹FGM板的I型和II型裂纹SIFs、最大能量释放率、裂纹扩展角和非均匀载荷作用下功能梯度板的动态SIFs进行了研究.结果表明，此方法具有较高的精度，相比传统伽辽金法可节省大量计算时间^[8-9]，但计算效率仍不高.滕子浩等^[10]对含有间断的非均匀材料的断裂问题，将虚节点多边形单元的形函数引入到扩展有限元(extended finite element method，XFEM)中，提出了一种基于四叉树结构的动态网格细化方法，并对非均质材料中的裂纹扩展问题进行研究，结果表明该方法相比传统XFEM具有更好的精度、收敛性以及计算效率，但需对裂纹附近的网格进行细化，增加了工作量.由此可见，FGM结构的断裂行为备受关注，而现行的分析方法还存在诸多问题，半解析法应用范围有限，数值法虽有广泛的适用性和较好的精度及计算效率，但尚有改进空间.

本课题组研究建立的广义参数Williams单元^[11](简记为W单元)可直接求解裂尖SIFs，无需网格加密和后处理，精度很高，且已开发了相应的快速解法程序块^[12].目前，W单元的研究成果主要集中于均匀材料的断裂分析，鉴于当前FGM结构的断裂行为研究尚需深入，本文建立了带裂纹FGM薄板I-II复合型裂纹尖端SIFs分析的W单元新格式，并分析了弹性模量E的变化形式、裂纹长度及裂纹倾角对裂尖SIFs的影响，为带裂纹的FGM薄板裂尖SIFs求解提供了新思路.

1 FGM薄板的广义参数有限元模型

含中心斜裂纹均匀拉伸FGM薄板如图 1所示，板宽为W，高为H，厚为t，裂纹长度为2a，倾斜角为γ，薄板上下边界承受均布拉应力σ₀.以裂纹中心作为整体坐标原点o，水平方向为x轴，逆时针旋转90°为y轴，建立直角坐标系xoy.以任意裂尖为坐标原点o_i，沿裂纹扩展方向为x_i轴，逆时针旋转90°为y_i轴，建立裂尖局部坐标系x_io_iy_i(图中i取1或2).文献[2]对含裂纹FGM薄板研究表明:泊松比μ对裂尖SIFs的影响较小，可忽略，因此本文研究取μ为常数，仅考虑弹性模量E随坐标轴分布，并假设其为整体坐标的函数.

因裂尖应力场具有奇异性，将整个计算模型划分为裂尖奇异区域和外围常规区域.FGM矩形薄板常规区网格可利用ANSYS有限元软件中四边形8结点等参单元自动离散，如图 2所示，因此类单元中不同位置E值不同，称之为非均匀单元.裂尖奇异区网格离散如图 3所示，其离散过程如下:选取以裂尖o_i为中心点，边长为l的正方形微小区域为奇异区，将o_i与各边中点和角点连接形成8个相同三角形条元，利用一组环绕裂尖o_i且相互平行的折线将每个条元划分为n个梯形微单元和1个裂尖三角形微单元(红色区域)，任意相邻折线Γ_n，Γ_n-1到裂尖o_i的距离之比为常数α(α < 1)，α称为奇异区径向离散比例因子，n称为奇异区径向离散单元层数.选取条元o_iL₁L₂作为典型单元对其组成进行详细说明:最外层梯形微单元L₁L₂L₃L₄的边L₁L₂为常规区与奇异区的公共边界，其上3个黑色实心圆点表示常规单元结点，另外的5个空心圆点位于奇异区内，为虚结点，故称微单元L₁L₂L₃L₄为过渡单元；最内层三角形微单元因面积过小，可忽略其刚度贡献，因此，可用条元L₁L₂L₅L₆的刚度代替o_iL₁L₂整个条元的刚度，因条元内各微单元的位移场受Williams级数控制，故称之为W单元.

2 FGM薄板裂尖奇异区位移场与SIFs

裂纹尖端奇异区位移场采用Williams级数表示，并截取前m+1项，即

(1)

式中:u_i，v_i分别表示i裂尖奇异区内任一点的整体位移分量；G(x, y)=E(x, y)/[2(1+μ)]为剪切模量，由整体坐标和局部坐标关系可知x=x(r_i, θ_i)，y=y(r_i, θ_i)，(r_i, θ_i)表示i裂尖局部极坐标系下的坐标；f_h1^x，f_h2^x，f_h1^y，f_h2^y取值见文献[11]；a_ih，b_ih为待定系数，因无特定物理意义，称为广义参数.由文献[11]知，应力强度因子K_{i, I}，K_{i, II}与a_i1，b_i1有直接关系，即

(2)

3 建立FGM薄板整体控制方程 3.1 FGM薄板常规区单元控制方程

在FGM薄板中，8结点等参单元各结点处的弹性模量值不相等，根据有限单元法可将单元位移函数表示成结点位移函数，将FGM薄板中常规区非均匀单元内部任意位置处的弹性模量表示为结点弹性模量的函数，即

(3)

式中:E_p为结点p所对应的弹性模量；为单元插值函数，与位移插值函数一致.

根据有限单元法，建立FGM薄板中常规区非均匀单元控制方程:

(4)

式中: 为单元刚度矩阵；为单元结点位移向量；为单元结点荷载向量.其中，

(5)

其中: 为单元应变矩阵；为单元插值函数矩阵；f为作用在二维体边界上的面积力; Ω^e为积分区域；D(x, y)为弹性矩阵，因单元材料属性的非均匀性，其值与整体坐标(x, y)相关.

常规区中部分结点位于常规区与奇异区公共边界，即公共结点；其余结点位于外围常规区，则FGM薄板常规区单元整体控制方程可分块表示:

(6)

式中: 为刚度矩阵；为结点位移列阵；为常规区结点荷载列阵；下标r和b分别表示非公共结点和公共结点.

3.2 FGM薄板裂尖奇异区W单元控制方程

以图 3中的条元L₁L₂L₅L₆为例，说明i裂尖1个W单元控制方程的建立过程.

将四边形8结点等参单元位移向量表示为

(7)

令W单元中所有子单元的8个结点按相同顺序排列，则相邻两层梯形微单元内层与外层相同编号结点到裂尖的距离之比为α，由各层子单元的相似性可知其极坐标关系为.

(8)

根据广义参数有限元法，结点的局部位移场服从整体位移场，由式(1)和式(8)可知W单元中第k层四边形8结点等参单元结点位移向量为

(9)

式中：k(k=1, 2, 3, …, n)表示W单元中子单元所在层数；T_i^(k)为第k层子单元的转换矩阵;

为简化计算，避免重复计算W单元各层子单元转换矩阵中的T_{q, i}^(k)，利用式(9)，可得第k层子单元转换矩阵与第1层子单元转换矩阵中的T_{q, i}⁽¹⁾关系:

(10)

式中：

将式(10)带入式(9)可得

(11)

将式(11)代入非均匀单元控制方程式(4)可得

(12)

方程两边同乘可得W单元中第k层梯形微单元广义控制方程:

(13)

式中:

由式(14)可知，所有梯形微单元广义控制方程均含有相同的广义参数列阵φ_i，但式中K_i^(k)含有与剪切模量相关的量ω_i^(k)，与均匀材料奇异区控制方程直接等比叠加的集成方式不同，需将第2至第n层梯形微单元广义控制方程直接叠加，可得梯形条元L₃L₄L₆L₅的控制方程:

(14)

式中:

因过渡单元部分结点位于公共边界，有结点号，而内部5个结点位于奇异区，为虚结点，其位移场需进行广义变换.根据结点所处区域不同，将过渡单元控制方程进行分块:

(15)

式中： = ；下标s表示奇异区内虚结点.

将式(14)与式(15)进行集成，即得1个W单元控制方程，其中因奇异区内无荷载，有，即

(16)

式中：.

3.3 整体控制方程

同一裂尖，周围所有W单元均可按照式(16)形成控制方程，则i裂尖8个W单元控制方程可集成为

(17)

式中：为i裂尖所有过渡单元公共结点刚度矩阵；和为i裂尖所有过渡单元公共结点与虚结点关联刚度矩阵；为i裂尖奇异区内广义刚度之和；裂尖所有公共结点位移向量.

对于含有λ个裂尖的FGM薄板，根据常规区控制方程式(6)与单个奇异区所有W单元控制方程式(17)，将所有局部坐标下的裂尖奇异区控制方程转换为整体坐标后与常规区控制方程进行集成，可得模型整体控制方程:

(18)

式中： + , ，.

引入应力与位移边界条件，求解式(18)，即可求得i裂尖广义参数a_i1，b_i1，进而根据式(2)可以直接求得各裂尖应力强度因子K_{i, Ⅰ}，K_{i, Ⅱ}.为便于与文献解进行对比分析，对SIFs进行无量纲化，无量纲化SIFs可定义为 / .

4 算例分析 4.1 中心斜裂纹

取一含中心斜裂纹的FGM矩形薄板，如图 1所示.板宽W=0.4 m，高H=0.6 m，厚度t取单位厚度，裂纹长度2a=0.2 m，倾斜角为γ，薄板上下边界承受均布拉应力σ₀，以裂纹中心点为坐标原点o，水平方向为x轴，逆时针旋转90°为y轴，建立整体坐标系xoy，泊松比μ取0.3.分别取弹性模量E呈指数型和线性型分布进行分析:

(19)

式中:β₁和β₂为材料非均匀参数；E₀表示坐标原点处弹性模量.另外，板左、右边界处的弹性模量分别用E_L，E_R表示，上、下边界处弹性模量分别用E_T，E_D表示.

将FGM板进行有限元网格离散，奇异区尺寸取0.02 m，每个奇异区离散为8个扇形条元，外围常规区采用大型通用有限元分析软件ANSYS离散为684个四边形8结点等参单元，整个计算模型共有2 172个结点.计算结果如下:

1) 弹性模量E沿x轴呈指数型分布:取板的左边界弹性模量E_L=380 GPa，右边界弹性模量E_R=116 GPa，式(19)中β₂=0，经计算可得材料非均匀参数β₁=-0.029，选取不同的裂纹倾角γ，研究其与k_{i, g}的变化规律，此模型与文献[13]一致，并比较其结果，如图 4所示.

文献[13]通过大型通用有限元软件ANSYS中的1/4奇异单元对模型进行求解，由图 4可知:本文解与文献[13]解吻合很好，证明了本文方法用于求解FGM裂纹问题的正确性.当弹性模量E沿x轴呈指数型分布时，左右裂尖的k_{i, Ⅰ}及k_{i, Ⅱ}变化趋势一致；随着裂纹倾角γ的增大，左右裂尖k_{i, Ⅰ}均逐渐减小，且弹性模量大的一侧k_{i, Ⅰ}越大，k_{i, Ⅱ}则均呈现先增大后减小的趋势，左右裂尖值较接近.

2) 弹性模量E呈线性型分布:取板的左边界弹性模量E_L=380 GPa，整体坐标原点处的弹性模量E₀=E(0, 0)=210 GPa，研究弹性模量E沿x轴呈线性型分布时对k_{i, g}的影响规律，计算结果如图 5所示.

对比图 5与图 4的计算结果发现:二者只在值的大小上有细微差别，变化趋势一致.由此可知:弹性模量E沿x轴的任意单调变化形式均不改变k_{i, g}的变化趋势，且当弹性模量E沿x轴呈指数型或线性型分布时，仅对k_{i, I}有较为显著的影响，且弹性模量较大一侧k_{i, I}较大，弹性模量较小一侧k_{i, I}较小.

3) 弹性模量E沿x，y轴同时呈指数型分布:取E(-W/2, -H/2)=380 GPa，E(W/2, H/2)=116 GPa，E₀=210 GPa，经计算可得材料非均匀参数β₁=β₂=-0.029.研究弹性模量E满足式(19)，沿x，y轴同时呈指数型分布时，k_{i, g}的变化规律，计算结果如图 6所示.

将图 4与图 6进行比较可知:二者k_{i, Ⅰ}的值相同，而k_{i, Ⅱ}的值出现了较大差异，说明沿y轴方向分布的弹性模量E改变时，主要对两裂尖的k_{i, Ⅱ}值产生影响.进一步分析可知:当弹性模量E分布形式呈单调变化且梯度与荷载方向平行或垂直时，分别使中心斜裂纹两裂尖的k_{i, Ⅰ}或k_{i, Ⅱ}值产生差异.

4.2 边界裂纹

取一含边界裂纹的FGM矩形薄板，如图 7所示.板宽W=1 m，高H=2 m，厚度t为单位厚度，裂纹长度为a，倾斜角度γ=0°，薄板上下边界承受均布拉应力σ₀，以裂纹与左边界的交点为坐标原点o，水平方向为x轴，逆时针旋转90°为y轴，建立整体坐标系xoy.泊松比μ=0.23，分别取E为指数型和线性型进行分析，其表达式为

(20)

将FGM板进行有限元网格离散，奇异区尺寸为0.025 m，奇异区离散为8个扇形条元，外围常规区采用大型通用有限元软件ANSYS离散为817个四边形8结点等参单元，整个计算模型共有2 588个结点.计算结果如图 7所示.

1) 令E₁=E_R=dE(0, 0)，分别取d=0.1，0.2，0.5，2.0，5.0，10.0，通过设置不同的裂纹长度a，研究不同裂纹长度下弹性模量梯度对裂尖SIFs的影响，此模型与文献[14-15]一致，并将计算结果分别与文献[14-15]解进行比较，如图 8及图 9所示.

文献[14-15]分别采用奇异积分方程组及1/4奇异单元求解含边界裂纹FGM薄板裂尖SIFs，两种方法的计算精度均相对较高.由图 8、图 9可知:本文解与文献[14-15]解吻合得很好，证明了本文方法用于求解边界裂纹FGM薄板结果的正确性.当弹性模量E沿x轴呈指数型或线性型分布时，k_I的变化趋势均一致；当裂纹长度a为0.4，0.6 m时，k_I随着E₁/E₀的增大而逐渐减小；而当裂纹长度a=0.2 m时，k_I随E₁/E₀的增大先增大后减小，此时是由于裂纹尺寸过小，裂尖过于靠近边界，且E₁/E₀ < 0.5时梯度过大共同作用造成此现象.

2) 令E₁=E(W, H/2)=dE(0, 0)，分别取d=0.1，0.2，0.5，2.0，5.0，10.0，研究当弹性模量E沿x，y轴同时呈指数型分布时，由于弹性模量梯度与裂纹面不平行，导致k_I，k_II同时存在，本文计算结果如图 10所示.当弹性模量E沿x，y轴同时呈指数型分布时，k_I的变化趋势与弹性模量E沿x轴呈指数型或线性型分布时一致，k_I随着弹性模量比E₁/E₀的增大而逐渐减小，但当裂纹长度a=0.2 m时，k_I在E₁/E₀接近于0时，会随着E₁/E₀的增大而增大；而k_II则随着弹性模量比E₁/E₀的增大而逐渐增大.

5 结论

1) 结合材料属性变化的W单元能有效地解决FGM薄板平面断裂问题，以裂尖SIFs作为基本未知量在刚度方程中直接求解，得到高精度结果.

2) FGM结构中裂尖位置处弹性模量越大的一侧裂尖SIFs越大.分布形式呈单调变化的弹性模量E，当其梯度与荷载方向平行或垂直时，将使中心斜裂纹两个裂尖的I型或II型SIFs值产生差异.

3) 弹性模量E沿水平方向呈线性型分布、指数型分布或水平及竖直双向呈指数分布时，SIFs随裂纹倾角γ的变化趋势相同，即在含裂纹的FGM薄板断裂分析中，当弹性模量E分布形式呈单调变化时，裂尖SIFs随裂纹倾角γ的变化规律不受弹性模量E分布形式的影响.