功能梯度材料圆筒(以下简称 FGM 圆筒)是一种可设计性非均匀复合材料结构，因其消除了传统复合材料结构中存在的界面问题而被广泛应用于航空航天、汽车、船舶与海洋工程等领域[1-2]。由于FGM圆筒服役期间通常处于极端环境，极易出现裂纹损伤。裂纹的存在不仅会影响结构性能，还有可能导致整个工程结构失效[3]。因此，对含裂纹FGM 圆筒结构进行力学分析是保障工程结构安全性与先进性的先决条件。
裂纹尖端的应力强度因子是确定含裂纹构件安全性的重要参数[4-7]。近年来，已有大量研究论文涉及了这一领域[8-13]，基于权函数法[8-9]、能量释放率法[12]等方法研究含裂纹FGM圆筒的应力强度因子。已有研究结果显示，可以将针对含裂纹均匀材料的理论分析方法加以改进，应用于含裂纹FGM圆筒应力强度因子的研究中[8, 11]。但是，由于FGM性质的特殊性，FGM结构力学行为比均匀材料结构复杂很多。因此，尽管均匀材料结构应力强度因子已经可由经验公式得到[14]，FGM 结构应力强度因子计算过程却仍然难以摆脱复杂矩阵运算和数值积分[8-13, 15]。课题组前期研究发现，可以基于参考问题模型与未知问题模型之间的稳定应力关系简化奇异应力场强度的分析步骤，但是研究对象仅限于胶接结构，且界面端的材料组合以及几何尺寸需完全相同[16-19]。如果能够找到均匀材料结构与FGM结构应力强度因子之间存在的稳定关系，便有可能将复杂的FGM结构应力强度因子求解问题转化为简单的系数计算以及均匀材料结构应力强度因子经验公式求解问题。
为此，本文在课题组已有研究[16-19]的基础上，探寻FGM圆筒与均匀材料圆筒这两者裂纹尖端应力与应力强度因子之间存在的稳定关系，最终基于裂纹尖端应力比值，提出了一种适用于含环状裂纹FGM圆筒的应力强度因子高精度快速计算方法。为了适应工程需要、降低计算难度，本文研究时不考虑裂纹尖端非奇异应力的影响[20]。将本文计算结果与参考文献计算结果进行对比，验证了本文方法的可行性和优越性。
1 力学模型
采用有限元软件MSC. Marc Mentat 2012.1.0建立含环状裂纹FGM圆筒2D轴对称模型(见图1)，裂纹深度为a，圆筒长度为2L，厚度为h，内外半径分别为Ri和Ro，材料沿厚度方向连续性变化，圆筒两端施加拉力hσ。
图1 FGM圆筒结构示意图
Fig.1 Dimension for FGM hollow cylinder
2 理论推导
2.1 功能梯度材料
假设功能梯度材料性质沿厚度r′方向按照指数函数分布形式连续性变化[11]，则FGM圆筒内、外表面材料的体积分数之间存在以下关系：
式中：Vi(r′)为内表面材料体积分数；Vo(r′)为外表面材料体积分数；p为材料梯度指数(0 ≤ p ≤ ∞ )。当p=0时，FGM圆筒退化为均匀材料圆筒。
当FGM圆筒由两种材料组成时，圆筒的材料性质表现为从圆筒内表面材料到圆筒外表面材料性能的连续性变化。基于线性混合法则[21]，圆筒的弹性模量E和泊松比ν可表示为：
式中，Ei和Eo、iν和oν分别为FGM圆筒内外表面材料的弹性模量、泊松比。
2.2 应力强度因子
根据线弹性断裂力学，通常情况下的FGM结构裂纹尖端应力场表达式如下所示[22]：
式中：i, j = 1 ,2,3对应空间坐标系的三个方向；r和θ为裂纹尖端极坐标(见图 2)；
为无量纲角函数；K I、 KII、 KIII为应力强度因子，与材料梯度、载荷以及几何尺寸有关。
图2 圆筒裂纹尖端几何模型
Fig.2 Schematic of crack tip in hollow cylinder
对于图1所示的两端施加均匀轴向拉力的含环向裂纹FGM圆柱壳，其所对应的Ⅰ型裂纹尖端应力场可表示为：
材料梯度会影响应力强度因子数值，但是不会影响应力奇异性以及无量纲角函数[23]。因此，理论上可以利用“比例”形式消除奇异性以实现在保证精度的基础上降低计算难度、简化计算过程。
无量纲应力强度因子FI和应力强度因子KI满足关系式：
当FGM圆筒和均匀材料圆筒裂纹受力情况、裂纹类型、裂纹尖端附近网格划分情况等均相同时，这两者的裂纹尖端附近应力比值与相应应力强度因子比值之间存在如下关系：
由式(5)和式(6)可知，当FGM圆筒与均匀材料圆筒两端受力hσ相同时，
和裂纹尖端应力σy ( r)之间存在如下比例关系：
即FGM圆筒的
值可由式(8)计算得到：
式中： a F GM 和a均匀材料可根据实际情况自由选择，两者数值可不相同；
为远距均匀轴向拉伸均匀材料圆筒的无量纲应力强度因子(见式(9))[14]。
3 算例与讨论
建立图3所示1/4圆筒模型即可满足计算需求，圆筒两端施加拉应力hσ=1 MPa，h=1 mm，使用四节点四边形单元，裂纹尖端附近网格加密，em in=1/38 mm、 R i /Ro=0.8时单元总数为23 289，其中裂纹尖端附近单元数约为20 000。本文中，设定FGM圆筒外表面材料弹性模量Eo数值固定，内表面材料弹性模量Ei可变。假设泊松比为常数，ν=0.3。当
o/ i
E E=1时，FGM退化为均匀材料。
图3 模型网格划分示意图Fig.3 Mesh pattern for model
3.1 裂纹尖端应力分布
FGM圆筒外表面材料体积分数Vo沿厚度方向随体积函数幂指数p变化(见图4)，其可设计性可以通过控制内外两种材料组分比例实现[24]，Vo=1-Vi。当p≤0.2或p≥10时，大部分Vo值接近1或者0。以p≥10为例，虽然材料性能仍然由内表面材料沿厚度方向连续性变化至外表面材料，但Vo初始变化幅度非常小，圆筒内侧材料对材料性能起主导作用；靠近筒外壁时，Vo变化幅度接近最大。
图4 FGM圆筒外表面材料沿厚度方向体积分数Vo Fig.4 Variation of volume fraction of outside material Vo of FGM hollow cylinder in the thickness direction
由于裂纹尖端奇异性的存在，使用 Marc得到的应力值yσ不仅受圆筒材料性质影响，还会随裂纹尖端附近最小网格尺寸emin以及裂纹尖端距离r的改变而发生变化(见图 5)，不能直接使用。当h/ L ≥ 5 时，h/L对应力值的影响可以忽略不计。圆筒几何参数为：L/h＞5，R i /Ro =0.8，a/h=0.1。
图5 材料性质和emin对裂纹尖端应力分布的影响Fig. 5 The influences of material property and emin on the stress distribution at crack tip
当FGM圆筒与均匀材料圆筒厚度、载荷工况以及裂纹尖端附近网格划分等完全相同、emin又小到一定程度时，两者裂纹尖端应力比值稳定，可以通过比值形式消除应力奇异性。
由表 1中的 a)～b)可见，当 h / L ≥ 5 、em in=1/34 mm时，裂纹端部应力值随r/h增大而减小，应力比值
变化幅度较大。当em in=1/38 mm时，应力值仍然随r/h增大而减小，但是
已经基本恒定不变；表1中的c)，FGM 圆筒L/h=2，均匀材料圆筒仍然使用h/ L ≥ 5 。结果显示，应力比值变化规律与表1中的a)～b)中应力变化情况相同，仍然在emin=1/38 mm时达到稳定状态，并未受到筒体长度影响。因此，可以使用应力比值进行应力强度因子计算，并且仅使用裂纹尖端节点应力值 σ y0,FEM[16]即可满足计算需求，此结论与文献[16—19]相同。下文中，应力值σy均为σy 0, FEM 。
表1 不同网格尺寸时裂纹尖端应力分布
)
Table 1 Stress distributions for crack tip obtained by different mesh sizes (Ri /Ro = 0 .8, ( E o / E i ) F GM= 0 .4) of
a) / 5 L h≥ ,/ 0.1 a h= ofe = 1/38 mm min e = 1/34 mm min FGM ya/h0.1 σ = /σ均匀材料= /FGM y / 0.1 ya/h0.1σ σ FGM y / 0.1 r/h σ = /MPa均匀材料 r/h= =σ均匀材料= /FGM ya/h0.1ya/h0.1均匀材料MPa MPa y / 0.1 ah ah MPa σ σy / 0.1 ah ah= =→0 25.5710 41.719 0.613 →0 2.8120 4.5687 0.615 1/ 3812.7970 20.870 0.613 1/34 1.5169 2.3967 0.633 2/ 38 9.3436 15.233 0.613 2/34 1.1916 1.8388 0.648 3/ 38 7.4826 12.194 0.614 3/34 1.0298 1.5511 0.634 4/ 38 6.4610 10.525 0.614 4/34 0.9538 1.4049 0.679 b) / 5 L h≥ , / 0.8 ah= (FGM), / 0.1 ah= (均匀材料)e = 1/38 mm min e = 1/34 mm min σ σ σ = /FGM ya/h0.8σ均匀材料= /ya/h0.1σ σ FGM y a/h 0.8 ah= =σ = /MPa ofr/h ofr/h FGM ya/h0.8σ均匀材料= /ya/h0.1 FGM y / 0.8均匀材料 均匀材料MPa ofMPa ofy / 0.1 ofMPa ofy / 0.1 ah ah= =→0 251.070 41.719 6.01 →0 27.8580 4.5687 6.10 1/ 38 125.520 20.870 6.01 1/3413.8180 2.3967 5.77 2/ 38 91.548 15.233 6.01 2/34 9.9266 1.8389 5.40 3/ 38 73.230 12.194 6.01 3/34 7.7770 1.5511 5.01 4/ 38 63.160 10.525 6.01 4/34 6.5358 1.4049 4.65 c) / 2(FGM), / 5( )Lh Lh=＞ 均匀材料 ofe = 1/38 mm min e = 1/34 mm min σ = /σ均匀材料= /ya/h0.1 FGM ya/h0.1σ = / ofMPa ofσ σ FGM y / 0.1均匀材料 ofr/h ofr/h= =σ均匀材料= /FGM ya/h0.1ya/h0.1均匀材料MPa FGM y / 0.1 MPa σ σy / 0.1 ah ah MPa y / 0.1 ah ah= =→0 33.9930 41.719 0.815 →0 3.7346 4.5687 0.817 of1/ 38 17.0100 20.870 0.815 1/34 2.0012 2.3967 0.835 of2/ 38 12.4190 15.233 0.815 2/34 1.5608 1.8389 0.849 of3/ 38 9.9442 12.194 0.815 3/34 1.3385 1.5511 0.863 of4/ 38 8.5856 10.525 0.816 4/34 1.2305 1.4049 0.876 of
本文中，emin=1/38 mm已经可以满足计算需求。若进一步细化网格，采用更小的emin，单元数量明显增多，运算时间随之延长，但对计算结果基本没有影响。因此，不宜过度细化网格。
3.2 应力强度因子计算
求解时，由式(8)可得含裂纹均匀材料圆筒无量纲应力强度因子
，再将由 Marc得到的裂纹尖端应力值代入式(9)即可得到含裂纹FGM圆筒无量纲应力强度因子
，具体计算结果见表2。圆筒几何参数为： L / h ≥ 5 ， R i / Ro =0.8。
表2 无量纲应力强度因子
计算Table 2 Calculation of dimensionless stress intensity factor F ⅠFGM
含裂纹均匀材料圆筒 of 含裂纹FGM圆筒(o/i E E=0.4, p=1) ofa / h ofyσ均匀材料/ ofF均匀材料FGM yσ/ ofI F ofa/h ofMPa ofFGM I F MPa ofF均匀材料 FGM I I 0.1 of 25.572 of 0.6130 of 0.6959 0.2 of 41.556 of 0.7044 of 0.7996 0.3 of 59.597 of 0.8248 of 0.9363 0.1 of 41.719 of 1.1352 of0.4 of 81.484 of 0.9766 of 1.1087 0.5 of 108.800 of 1.1663 of 1.3240 0.6 of 143.400 of 1.4033 of 1.5931 0.7 of 188.400 of 1.7069 of 1.9378 0.8 of 251.070 of 2.1277 of 2.4155
本方法中，仅需一组均匀材料圆筒计算结果
及
)即可满足各种情况下的 F FGM计算需求。均匀材料圆筒a/h可根据计算需求自由选择，表2中仅以a/h=0.1为例进行说明。为了方便对比，除非特别说明，下文中均选用表 2中的
MPa进行
计算。
图6中针对FGM圆筒长厚比L/h、筒体内外半径比 R i / Ro 、裂纹深度与圆筒壁厚比a/h、FGM性质(内外表面材料弹性模量比 E o /Ei ，材料梯度指数 p)等
影响因素，将本文方法计算结果与参考文献[11—12]结果进行对比了分析。从中可以看出，本文计算方法结果与文献计算结果吻合度较好，验证了本文方法的正确性和有效性。
文献[12]中，材料常数 γ = E o /Ei = e β ，弹性模量 E (r′ ) =
。
本文方法是基于裂纹尖端应力比值，而文献采用位移有限元法、柔度法、有限差分法等方法，计算方法的不同是造成误差的主要原因。由图3可知，FGM圆筒材料性质沿厚度方向随幂指数p变化。当p=10时，靠近筒外壁处(a/h=0.1)Vo值变化幅度很大，对计算结果影响较为明显。FGM圆筒内、外表面材料性质差异较大时，p影响更为显著。
图6 无量纲应力强度因子
对比分析
Fig.6 Comparison of dimensionless stress intensity factor 
由图7可见，近90%的本文计算结果与参考文献计算结果之间的误差小于 4%，仅有一处数据误差值超过7%。当h/L＜5、 a / h ≥ 0 .2时，误差接近0%(见图7)。当γ减小或 E o /Ei 增大时，误差有增大趋势，但最大误差可控制在6.5%以内。选取较大误差位置所对应参数进行校核( E o /Ei =10、p=10、a/h=0.6)，计算结果
=1.9067，参考文献计算结果为1.9725，误差为-3.3367%，仍保持在合理范围内。壳体参数类型改变时，计算误差也并未出现明显的规律性变化，计算结果稳定。
图7 计算误差对比分析
Fig.7 Comparison of calculation errors
本文使用的
经验公式是由 Neuber的近似应力集中系数法导出，该经验公式对于实心轴情形(Ri→0)时与别的近似解存在百分之几以内误差[7]，因此不建议进行实心轴情形时的
计算。
计算时，采用不同的
可能会造成同一L / h ≥ 5 以及L/h＜5这两种情况，以图 7中误差最大位置处数据为例，给出了
不同时同一FGM圆筒的
计算结果。结果显示，当L/h＜5时，
1.31±0.069 (
为
的平均值±标准偏差)，该结果与参考文献计算结果[11]吻合较好；当 L / h ≥ 5 时， 
=1.73±0.091，该计算结果较参考文献计算结果[12]总体偏小，采用a/h较小时的
可有效减小本方法计算结果与参考文献结果之间的误差。因此，计算时可考虑采用a/h较小时的
进行
计算。
FGM圆筒的
计算结果不同。表3中分别针对
表3 
对
计算结果的影响
Table 3 The effect of 
on the results of 
含裂纹F GM圆筒 of/=2 E E ,p=10,i/o=0.8含裂纹均匀材料圆筒 ofL h≥ ，o/i=10/ 5 R R ,L h ,γ=2,i/o=0.5 R R , ofFGM FⅠ =1.973 [12] ofFGM FⅠ =1.352[11] ofa/h yσ均匀材料F均匀材料 a/h FGM yσ/FGMyσ/I F ofa/h ofFGM I FFGMI/MPa ofMPa ofMPa of0.1 41.72 1.135 1.907 1.440 0.2 63.65 1.171 1.823 1.377 0.3 85.55 1.228 1.742 1.315 0.4 109.40 1.310 1.677 1.267 0.5 136.40 1.426 1.638 1.237 0.6 167.80 1.596 1.632 1.233 0.7 206.00 1.855 1.670 1.261 0.8 258.20 2.303 0.6 171.60.1 52.92 1.7681.335
由于本文计算方法基于裂纹尖端应力比值，裂纹尖端附近网格划分需满足要求，需细化至裂纹尖端附近应力比值基本稳定不变时再进行计算，且FGM 圆筒与均匀材料圆筒裂纹尖端附近网格划分必须一致。
与文献中使用的计算方法相比，本文计算方法避免了复杂的矩阵运算以及数值积分，计算精度高，计算量小，适用范围广。计算时，均匀材料圆筒
可由经验公式获得，过程简便，适合工程应用。本文方法不仅适用于含外表面环状裂纹FGM圆筒，还可用于其他裂纹形式的FGM结构的应力强度因子计算，但是现有研究结果还很有限，本项目后续工作将进行深入研究。
4 结论
本文基于FGM圆筒与均匀材料圆筒裂纹尖端应力值与应力强度因子之间的稳定比例关系，提出了一种适用于含外表面环状裂纹的FGM圆筒应力强度因子计算方法。该方法在保证精度的基础上成功避免了复杂的矩阵运算和数值积分，仅通过使用含环状裂纹均匀材料圆筒应力强度因子经验公式、FGM圆筒与相应均匀材料圆筒裂纹尖端应力比值，即可得到FGM圆筒应力强度因子，计算过程简单。计算结果显示，本文方法适用于不同FGM类型、圆筒尺寸以及裂纹深度时的应力强度因子计算，计算时所用的
可任意选取，且仅需计算一次
即可满足不同情况下的
计算需求，计算精度高、结果稳定，便于工程应用。
[1] Li X, Du C C, Li Y H. Parametric instability of a functionally graded cylindrical thin shell subjected to both axial disturbance and thermal environment [J].Thin-Walled Structures, 2018, 123: 25-35.
[2] Nasirmanesh A, Mohammadi S. Eigenvalue buckling analysis of cracked functionally graded cylindrical shells in the framework of the extended finite element method[J]. Composite Structures, 2017, 159: 548-566.
[3] 彭亦鹏, 陈爱军. 热冲击下带裂纹功能梯度厚壁圆筒的分析[J]. 振动与冲击, 2017, 36(15): 64-76.Peng Yipeng, Chen Aijun. Numerical simulation for a functionally graded thick-walled cylinder with cracks under thermal shock [J]. Journal of Vibration and Shock,2017, 36(15): 64-76. (in Chinese)
[4] Farahani B V, Tavares P J, Moreira P M G P, et al. Stress intensity factor calculation through thermoelastic stress analysis, finite element and RPIM meshless method [J].Engineering Fracture Mechanics, 2017, 183: 66-78.
[5] Ramezani M K, Purbolaksono J, Andriyana A, et al.Empirical solutions for stress intensity factors of a surface crack in a solid cylinder under pure torsion [J].Engineering Fracture Mechanics, 2018, 193(15): 122-136.
[6] 文龙飞, 王理想, 田荣. 动载下裂纹应力强度因子计算的改进型扩展有限元法[J]. 力学学报, 2018, 50(3):599-610.Wen Longfei, Wang Lixiang, Tian Rong. Accurate computation on dynamic SIFs using improved XFEM [J].Chinese Journal of Theoretical and Applied Mechanics,2018, 50(3): 599-610. (in Chinese)
[7] 王珊. 基于Kirchhoff板中裂纹尖端辛本征解的有限元应力恢复方法[J]. 工程力学, 2018, 35(5): 10-16.Wang Shan. A stress recovery method for cracks in Kirchhoff plate based on the symplectic eigen-solutions near the crack tips [J]. Engineering Mechanics, 2018,35(5): 10-16. (in Chinese)
[8] Seifi R. Stress intensity factors for internal surface cracks in autofrettaged functionally graded thick cylinders using weight function method [J]. Theoretical and Applied Fracture Mechanics, 2015, 75: 113-123.
[9] Eshraghi I, Soltani N. Stress intensity factor calculation for internal circumferential cracks in functionally graded cylinders using the weight function approach [J].Engineering Fracture Mechanics, 2015, 134: 1-19.
[10] Ooi E T, Natarajan S, Song C, et al. Crack propagation modeling in functionally graded materials using scaled boundary polygons [J]. International Journal of Fracture,2015, 192(1): 87-105.
[11] Li C, Zou Z. Stress intensity factors for a functionally graded material cylinder with an external circumferential crack [J]. Fatigue & Fracture of Engineering Materials &Structures, 1998, 21: 1447-1457.
[12] Chen Y Z. Stress intensity factors in a finite cracked cylinder made of functionally graded materials [J].International Journal of Pressure Vessels and Piping,2004, 81(12): 941-947.
[13] Mahbadi H. Stress intensity factor of radial cracks in isotropic functionally graded solid cylinders [J].Engineering Fracture Mechanics, 2017, 180: 115-131.
[14] 中国航空研究院. 应力强度因子手册[M]. 北京: 科学出版社, 1981.Chinese Aeronautical Establishment. Handbook of stress intensity factor [M]. Beijing: China Science Publishing &Media Ltd, 1981. (in Chinese)
[15] Li J, Liu J Z, Korakianitis T, et al. Finite block method in fracture analysis with functionally graded materials [J].Engineering Analysis with Boundary Elements, 2017, 82:57-67.
[16] Noda N A, Miyazaki T, Li R, et al. Debonding strength evaluation in terms of the intensity of singular stress at the interface corner with and without fictitious crack [J].International Journal of Adhesion and Adhesives, 2015,61: 46-64.
[17] Li R, Noda N A, Takaki R, et al. Most suitable evaluation method for adhesive strength to minimize bend effect in lap joints in terms of the intensity of singular stress field[J]. International Journal of Adhesion and Adhesives,2018, 86: 45-58.
[18] 李戎, 梁斌, NODA Nao-Aki. 基于临界应力强度因子的胶接对接接头抗拉强度简便预测方法[J]. 工程力学,2017, 34(11): 218-224.Li Rong, Liang Bin, Noda Nao-Aki. A convenient adhesive strength prediction method for adhesive butt joint in terms of the critical stress intensity factor [J].Engineering Mechanics, 2017, 34(11): 218-224. (in Chinese)
[19] Noda N A, Li R, Miyazaki T, et al. Convenient adhesive strength evaluation method in terms of the intensity of singular stress field [J]. International Journal of Computational Methods, 2018, 15(1): 1850085.
[20] Weißgraeber P, Leguillon D, Becker W. A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress raisers [J]. Archive of Applied Mechanics, 2016, 86: 375-401.
[21] Markworth A J, Ramesh K S, Parks W P. Modeling studies applied to functionally graded materials [J].Journal of Materials Science, 1995, 30: 2183-2193.
[22] Jin Z H, Noda N. Crack tip singularity fields in nonhomogeneous materials [J]. Transactions of the ASME/Journal of Applied Mechanics, 1994, 61(3):738-740.
[23] Rousseau C E, Tippur H V. Evaluation of crack tip fields and stress intensity factors in functionally graded elastic materials: Cracks parallel to elastic gradient [J].International Journal of Fracture, 2002, 114: 87-111.
[24] Li R, Liu L X, Liang B, et al. An efficient method to improve the stability of submerged functionally graded cylindrical shell [J]. Journal of Mechanical Science and Technology, 2019, 33(6): 2527-2536.