功能梯度材料板内部I型裂纹数值分析

功能梯度材料层合板内跨界面I 型裂纹数值分析2007.3.9材料属性随着空间坐标连续变化的功能梯度材料具有很好的耐热、耐磨、强韧等性能，克服了传统单一或组合材料属性不连续的缺点，已经在航空航天等领域引起了广泛关注。

材料在使用过程中失效主要来源于断裂破坏，因此功能梯度材料的断裂性能倍受关注。

上世纪80年代以来，许多学者对功能梯度材料断裂性能进行了理论研究。

1988年，Delale 和Erdogan [1]首先研究了无限大非均匀板中的内部裂纹问题。

发现功能梯度材料中裂纹尖端的应力场仍具有逆平方根奇异性，并且泊松比对应力强度因子的影响很小。

1994年，Jin 和Noda [2]对比了均匀材料和功能梯度材料裂纹前沿附近的物理场，并指出：对于一般的非均匀材料，只要材料属性连续且分段可微，则非均匀材料裂纹尖端场的奇异性与均匀材料的一致。

Erdogan 和Wu [3]研究了不同载荷下含有垂直于边界的裂纹的各向同性功能梯度条，给出了应力强度因子随着裂纹长度等参数变化的规律。

Wang [4]等人则考虑了具有垂直于边界的裂纹的功能梯度条，他们通过把材料沿着梯度方向分成了一些均匀材料层进行分析。

Kim 和Paulino 领导的科研组对功能梯度材料数值计算进行了系统的研究，Jin 和Paulino [5]提出了将单元节点处的材料属性使用与位移相同的形函数进行插值，得到单元内部任意点处的材料信息。

本文采用基于非均匀单元的有限元方法[6]研究了含有垂直于边界的边裂纹的三维各向同性功能梯度材料在均匀载荷作用下断裂问题，利用位移法求解了功能梯度材料层合板内跨界面I 型裂纹前沿的应力强度因子，给出了裂纹尖端到达层板界面处时，应力强度因子的变化规律。

基于非均匀单元的有限元方法有限元方法中形成单元刚度矩阵时，需要如下积分：(,,)e e T VK B D x y z BdV =⎰ (1)这里B 为应变转换矩阵，D 为材料矩阵，e V 为单元体积。

本文采用20节点等参单元，如图1-a 所示。

同时为了方便数值计算，将积分离散，本文采用三点Gauss 积分，公式(1)中的单元刚度矩阵e K 经过坐标变换后可以写成：333,,111((,,))i j keT j k i i j k K W W W B D B Jξξηηζζξηζ=======∑∑∑ (2)其中，ξ、η、ζ为等参坐标系O-ξηζ中的坐标，J 为坐标系O-xyz 变换到坐标系O-ξηζ的Jacobian 矩阵，(,,)i j k ξηζ为Gauss 积分点，i W 、j W 和k W 分别是三个坐标方向的权重值。

传统有限元整个单元取同一材料属性。

这里，对于梯度材料属性，令公式(2)中材料矩阵(,,)D ξηζ直接取Gauss 积分点(,,)i j k ξηζ处对应直角坐标系中(,,)i j k x y z 处的实际材料属性[6]。

每个单元内部含有27个高斯积分点，同一单元内部可以模拟材料的梯度变化，如图1-b 所示。

文[6]取一个单层功能梯度板，杨氏模量为指数形式，材料参数21ln(/)h E E δ=，将本方法计算出的应力强度因子值与Kadioglu [7]和Guo [8]的理论结果进行了对比，结果吻和很好。

图1 20节点梯度单元和单元内部的高斯积分点图(a)(b)图2 梯度层合板模型及几何参数x无量纲应力强度因子k 1/k 0a 0图3. 含有内部裂纹各向同性功能梯度板无量纲应力图4 具有形单元梯度层合板模型及数值求解含有垂直于纸面方向（z 方向）贯穿裂纹的各向同性功能梯度板，裂纹平行于梯度方向。

如图2所示（仅显示x-y 平面），功能梯度材料板长为L ，厚度分别为H 1和H 2。

厚度方向板界面处的材料杨氏模量分别为1E 、2E 和3E ，由于泊松比对功能梯度材料裂纹尖端应力强度因子的影响很小[1]，因此泊松比取为常数0.3ν=。

裂纹尖端的x 坐标为分别为a 和b 。

定义材料参数仅为x 的函数，如公式(3)所示。

2132ln(/)11ln(/)()212(0)()()E E x HE E x H HE ex H E x E e H x H -⎧≤<⎪=⎨⎪≤≤⎩ (3) 本文所研究的层板厚度21H H H ==的无限长板，取平面应变情况。

取长度与厚度的比值122/()8L H H +=足以模拟无限长板。

定义无量纲裂纹中心和裂纹长度分别为0()/2c b a H =+和0()/2a b a H =-，及无量纲应力强度因子0/In I K K K =，其中0K σ=0σ为裂纹面上均匀压载荷。

由于本程序为三维程序，因此在垂直于纸面方向（z 方向）取一个单元，限制所有节点的z 方向自由度。

裂纹前沿布置应力具有1/2r -奇异性的楔形单元，如图4所示。

使用位移法计算裂纹前沿的应力强度因子。

取00.2a =，针对公式(3)描述的材料属性，梯度参数21/E E 分别取0.1、1和10，32/E E 取为定值1.0。

裂纹面上施加均匀压载荷0σ。

裂纹尖端的应力强度因子随着裂纹中心位置0c 的变化曲线如图5所示。

无量纲应力强度因子K I n由图5可以发现，对于32/E E 取为定值1.0，当21/1E E <时，裂纹a 、b 两端的应力强度因子a K 和b K 都随着0c 增加而单调减小，在两层界面处出现尖锐的极致，过了界面后随着0c 增加而单调增加。

与之相反，当21/1E E >时，裂纹a 、b 两端的应力强度因子a K 和b K 都随着0c 增加而单调增大，在两层界面处出现尖锐的极致，过了界面后随着0c 增加而单调减小。

同时发现，随着裂纹中心0c 的增加，裂纹两端应力强度因子差值逐渐减小。

而当21/1E E =，即两梯度层模量按相同的梯度参数变化，不存在上述规律，应力强度因子变化比较平缓，仅与裂纹中心位置有关。

为了进一步考察材料参数对应力强度因子的影响，仍取00.2a =，针对公式(3)描述的材料属性，分别取21/E E 为0.2和5.0两种情况。

仍然在裂纹面上施加均匀载荷0σ，裂纹尖端应力强度因子随32/E E 的变化关系如图6。

由图6可以发现，对于21/E E 为0.2和0.5，裂纹a 端的应力强度因子a K 都随着32ln(/)E E 的增加而单调减小，b 端的应力强度因子随着32ln(/)E E 的增加或单调增加或单调减小。

同时发现，两层模量都为指数形式的材料也具有图5的性质，即在裂纹尖端到达两层界面时，应力强度因子出现极值。

而且发现两层材料模量以相同函数形式且梯度参数32/E E 与21/E E 相同时，裂纹尖端的应力强度因子值随裂纹位置改变的变化量最小。

图6 00.2a =不同位置的裂纹无量纲应力强度因子In K 随着材料参数32/E E 的变化关系：(a) 21/0.2E E =；(b) 21/ 5.0E E =ln(E 3/E 2)K I n(a) K I nln(E 3/E 2)(b)接下来验证两层材料模量为不同函数形式，应力强度因子是否在层板界面处出现极值。

仍取00.2a =，针对如公式(4)描述的材料属性，分别取21/E E 为0.2和5.0两种情况，梯度参数32/E E 分别取0.2、1和5.0。

仍然在裂纹面上施加均匀载荷0σ，裂纹尖端的应力强度因子随着裂纹中心位置0c 的变化曲线如图7。

21ln(/)1123212(0)()()(1)]()E E x HE e x H E x xE E E H x H H⎧≤<⎪=⎨⎪+--≤≤⎩ (4)由图7可以发现，对于两梯度层内模量按不同函数形式变化的情况，无论32/E E 是否等于21/E E 时，裂纹两端的应力强度因子a K 和b K 都会在各自尖端到达界面时出现极值。

而且从图7(b)可以发现，应力强度因子在未到界面时并不像图5那样具有随裂纹位置改变的单调性了。

经过进一步的计算发现：仅在梯度板两层模量为相同函数形式时，应力强度因子（在裂纹尖端未达层间界面时）才具有随裂纹位置变化的单调性；在梯度板层间界面处模量本身与其导数都连续的情况下，裂纹尖端到达界面时才不出现尖锐的极值，并且对于同一裂纹长度，裂纹位置对应力强度因子的影响较小。

结论本文采用基于非均匀单元的有限元方法研究了层合梯度板内跨界面I 裂纹随裂纹位置的变化关系。

结果显示：1）当裂纹尖端位于层间界面时，应K I nc 0K I nc0(a)(b)图7 两层梯度板内材料模量按不同函数形式变化的无量纲应力强度因子InK随裂纹位置变化关系：(a) 21/0.2E E =；(b) 21/ 5.0E E =力强度因子出现极值现象；仅在层间界面处模量本身与其导数都连续的情况下，才不出现尖锐的极值，并且对于同一裂纹长度，裂纹位置对应力强度因子的影响较小。

2）梯度板两层模量为相同函数形式且模量在层间连续时，应力强度因子（在裂纹尖端未达界面时）才具有随裂纹位置变化的单调性。

3）同一裂纹位置，应力强度因子随着梯度参数改变而单调变化。

参考文献[1] F. Delale, F. Erdogan. The crack problem for a nonhomogeneous plane. Journal of Applied Mechanics.1983, 50: 609-614[2] Z. H. Jin and N. Noda. Crack-Tip Singular Fields in Nonhomogeneous Materials. Journal of AppliedMechanics. 1994, 61:738~740[3] F. Erdogan and B.H. Wu. The Surface Crack Problem for a Plate with Functionally Graded Properties.Journal of Applied Mechanics. 1997, 64: 449-456[4] B.L. Wang, Y.W. Mai and N. Noda. Fracture mechanics analysis model for functionally gradedmaterials with arbitrarily distributed properties. International Journal of Fracture. 2002, 116: 161-177 [5] Z. H. Jin, G. H. Paulino and R. H. Dodds Jr. Finite element investigation of quasi-static crack growth infunctionally graded materials using a novel cohesive zone fracture model. Journal of AppliedMechanics. 2002, 69:370-379[6] 于红军, 吴林志, 果立成基于非均匀单元有限元方法的功能梯度材料三维断裂行为的数值模拟. 第十四届复合材料学术会论文集(下册) 2006, 714-717 中国宇航出版社.[7] S. Kadioglu, S. Dag, S. Yahsi. Crack Problem for a Functionally Graded Layer on an ElasticFoundation. International Journal of Fracture. 1998, 94:63-77[8] L.C. Guo, L.Z. Wu, T. Zeng,, L. Ma. Mode I crack problem for a functionally graded orthotropic strip.European Journal of Mechanics A/Solids. 2004, 23:219-234。