引 言
复合材料是由两种或多种不同性质的材料用物理或化学方法在宏观尺度上组成的具有新性能的材料．复合材料较单一材料强度高、刚度大、重量轻，并具有抗疲劳、减振、可设计等优点，数十年来在机械、建筑、航空、航天等领域日益得到广泛的应用[1]．复合材料的断裂破坏一般都始于材料界面的破坏，如钢筋混凝土是由钢筋和混凝土混合而成，钢筋和混凝土之间的界面为层间薄弱面，数值法和分析法均揭示裂纹易沿着界面发展[2-5]．裂纹结构全域的应力场是判断裂纹是否扩展和结构是否安全的重要参考依据．
目前已有大量文献研究双材料界面的应力强度因子[6-10]和裂纹尖端附近的应力场[11-14]．文献[11]采用插值矩阵法，求解得到了一般的塑性材料V形切口和裂纹前若干阶应力奇异指数和相应的特征函数；文献[12]通过构造新的应力函数，推出了Ⅰ型界面裂纹尖端的应力场、位移场的理论公式，其公式仅适用于求解裂纹尖端附近的应力场；文献[13]得出了仅考虑两相材料切口的两个奇异项来表征近应力场可能导致显著的误差，添加第一阶非奇异应力项可显著提高切口尖端近应力场的准确性；文献[14]讨论了两相材料黏合带应力强度因子的变化，但未获得裂纹结构的应力场．这些研究适用于裂纹尖端区域，但并未涉及远处乃至结构全域的应力场，对于不同弹性模量比值的双材料裂纹结构全域应力场准确性的分析也不多．如果一个方法能够给出两相材料裂纹结构位移和应力场，解的精度对于裂纹尖端附近和远处均准确适用，则其是圆满的解决路径．
本文基于裂纹尖端区域位移场的渐近级数展开式[15]，联合边界元法求解裂纹结构完整的位移和应力场，称之为扩展边界元法．其不仅可获得裂纹尖端区域的应力场，也可获得远处乃至结构全区域的应力场．本文将通过两个算例展现扩展边界元法求解两相材料裂纹结构全域应力场的准确性和适用性．
1 两相材料裂纹结构的应力奇性指数分析
图1 两相材料裂纹
Fig. 1 A bi-material crack
图1所示的两相材料平面裂纹，由具有不同材料特性的两个域Ω1和Ω2组成．Γ12(Γ21)是两相材料的结合界面，Ei和νi是域Ωi(i=1,2)的弹性模量和Poisson比．在裂纹尖端处同时定义一个xOy的直角坐标系和一个rOθ的极坐标系，坐标原点在裂纹尖端．
在裂纹尖端O点附近，以O为圆心挖去一半径为ρ的扇形域Ω″1∪Ω″2，如图2所示，扇形圆弧边界为Γ″1ρ∪Γ″2ρ，缺口处两径向边界为Γ″1和Γ″2，黏结交界为Γ″12(Γ″21)，如图2(b)所示．剩余结构区域为Ω′1∪Ω′2，域Ω′1的边界为Γ′1+Γ′1ρ+Γ′12+Γ3，域Ω′2的边界为Γ′2+Γ′2ρ+Γ′21+Γ4，如图2(a)所示．显然有Ωi=Ω′i∪Ω″i，Γj=Γ′j ∪Γ″j， j=1,2,12,21．
沿黏结面Γ″12(Γ″21)将图2(b)的扇形域剖分为两个域
Ω″1和Ω″2．根据线弹性理论分析，图2(b)所示扇形域Ω″1∪Ω″2在极坐标系rOθ(θ以逆时针方向为正)下的位移场可以表达成关于径向距离r的一系列级数渐近展开：
(1)
式中N为截取的级数项数，λk为应力奇性指数，
和
是沿r和θ方向的位移特征角函数，Ak为相应的位移幅值系数．
(a) 挖去扇形域后剩余结构(b) 裂纹尖端扇形域
(a) The remained structure without the sector(b) The sector around the crack tip
图2 黏结裂纹尖端挖去微小扇形
Fig. 2 A sector is removed from the bonded cracked structure
注意到式(1)的位移特征指数λk和位移特征角函数仅取决于裂纹的材料性质和楔形边界条件，不依赖载荷，因此应力特征分析不用考虑体力．取式(1)中的典型项代入弹性力学控制微分方程，通过一系列推导，可得裂纹尖端附近区域应力奇异性特征值问题的常微分方程组：
(2)
(3)
对于理想结合的黏结材料，在界面Γ″12(Γ″21)上满足位移和应力连续性条件：
(4)
设在Γ″1和Γ″2上面力自由，即有
(5)
如某Γ″i上视为固定约束，即有
(6)
因此，黏结材料裂纹尖端应力奇性指数λk的计算变成求解常微分方程组(2)和(3)；若两径向边界(Γ1和Γ3)自由，式(4)和(5)是相应的边界条件；若径向边界Γ1自由和径向边界Γ2固定，则式(4)～(6)是相应的边界条件．采用插值矩阵法求解方程得到裂纹的应力奇性指数λk及相应的特征角函数
及其导函数．值得指出，本文采用的插值矩阵法优点之一是求出的应力奇性指数λk、位移特征角函数和应力特征角函数
解的精度是同阶的．
2 扩展边界元法分析平面裂纹尖端附近应力场
将式(1)代入弹性力学几何方程和本构方程，可得尖端区域的应力分量：
(7)
式(7)中应力特征角函数为
(8)
由图2可知，边界Γ′1ρ∪Γ′2ρ和Γ″1ρ∪Γ″2ρ是分属于不同域的同一条边界，因此边界Γ′1ρ∪Γ′2ρ和Γ″1ρ∪Γ″2ρ上的位移连续，面力相等．若以uil,til(l=1,2)表示Γ′1ρ∪Γ′2ρ上的位移和面力分量，
表示Γ″1ρ∪Γ″2ρ上的位移和面力分量，则有
(9)
一般情况下
和
是复数，则图2(a)中边界Γ′1ρ∪Γ′2ρ上点的位移uil和面力til(l=1,2)在直角坐标系的表达式分别为
(10)
(11)
其中符号(·)R表示该复数的实部，(·)I表示其虚部．如果λk是单个实根，系数ζ=1；如果λk是共轭复根，系数ζ=2，因需计及复数共轭项对位移场和应力场的贡献．注意，式(10)、(11)中幅值系数AkR和AkI(k=1,2,…，N)待求．
图2(a)为无应力奇异性的普通弹性区域，沿其边界所有节点建立的边界积分方程如下：
y∈Γ3+Γ4,
(12)
式中
称为Cauchy主值积分；Cij(y)(i, j=1,2)称为位移奇性系数，参见文献[16]；bj为体力分量，y(y1,y2)是源点，x(x1,x2)是场点，域积分采用径向积分法[17]计算，积分核函数
是弹性力学Navier方程基本解，i, j=1,2．
弹性体区域Ω′i(i=1,2)内点的应力可分别由应力边界积分方程表达：
y ∈Ω′i,
(13)
式中积分核函数
和
为的线性组合，参见文献[18]．
扩展边界元法的分析步骤如下：
1) 将图1所示裂纹结构分成挖去扇形域的剩余结构(图2(a))和一个半径r=ρ的扇形域(图2(b))两部分．图2(a)区域没有应力奇异性，采用常规边界积分方程(12)，图2(b)区域为应力奇异性区域，用式(10)和(11)表示弧线边界Γ″1ρ∪Γ″2ρ上的位移和面力分量，两者联立求解得图2(a)边界上所有节点的未知面力、位移分量及式(1)中的幅值系数AkR,AkI(k=1,2,…,N)．扩展边界元法的优点之一是求得的AkR，AkI和对应的位移、应力函数的计算精度是同阶的，而常规有限元法和差分法求解的应力值精度比位移值低一个量级．
2) 将获得的边界节点位移和面力代入内点边界积分方程(13)，可得剩余结构区域Ω′1∪Ω′2内任意内点的应力．将幅值系数AkR，AkI代入渐近级数展开式(7)可求得尖端扇形附近区域的应力场．
至此，裂纹结构全域的应力场已完全获得．
3 算 例
目前对两相材料裂纹结构的研究局限于应力强度因子和裂纹尖端附近的应力场，而几乎没有对于其全域应力场的准确性研究．本文将通过算例展现扩展边界元法求解两相材料裂纹结构全域应力场的准确性和适用性．
含裂纹结构的单向拉伸，见图3．试件长H=190 mm，宽W=190 mm，裂纹深度为L=10 mm，拉伸面力P=100 MPa．线弹性平面应力问题，ν1=0.2,ν2=0.3,E1=206 GPa,E2/E1是变化的．
图3 复合型裂纹试件受拉伸和受剪切作用
Fig. 3 The bi-material crack subjected to
tension and shear load
本文用扩展边界元法对两相材料复合型裂纹结构E2/E1比值的极大情形(E2/E1=100)、较大情形(E2/E1=10，5)、相当情形(E2/E1=2)进行计算，展现扩展边界元法求解复合型裂纹结构弹性模量比值不同时全域应力场的准确性．
对于理想结合的黏结(连续性条件)材料，采用插值矩阵法计算常微分方程式(2)和(3)与边界条件式(4)和(5)，这里取N=12，获得裂纹尖端区域前12阶应力奇性指数λk，见表1．
当E2/E1=100时，域Ω2是域Ω1(参见图3)弹性模量的100倍，域Ω2相对于域Ω1是刚性的．因此对域Ω1单独计算应力奇异性，径向边界(Γ″1)自由，黏结边界(Γ″12)考虑为固支．计算常微分方程式(2)与边界条件式(4)～(6)，取N=5，获得域Ω1尖端前5阶应力奇性指数λk，见表2．对域Ω2单独计算应力奇异性，径向和黏结边界(Γ″2和Γ″21)自由．计算常微分方程式(3)和边界条件式(4)和(5)，取N=7，获得域Ω2裂纹尖端前7阶应力奇性指数λk，见表3．
表1 理想结合黏结材料的应力奇性指数λk
Table 1 Stress singularity indexes λk of the ideal bonded bi-material specimen
E2/E1λ1λ2λ3λ4λ5λ6100-1-10-0.499+i0.132-0.499-i0.1320E2/E1λ7λ8λ9λ10λ11λ121000.5+i0.1320.5-i0.1320.99911.5+i0.1321.5-i0.132
表2 域Ω1(一边自由和一边固结)裂纹尖端的应力奇性指数λk
Table 2 Stress singularity indexes λk in region Ω1 with free-fixed edges
λ1λ2λ3λ4λ5-0.499+i0.132-0.499-i0.13200.9991
表3 域Ω2(两边自由)裂纹尖端的应力奇性指数λk
Table 3 Stress singularity indexes λk in region Ω2 with free-free edges
λ1λ2λ3λ4λ5λ6λ7-1-100.5+i0.1320.5-i0.1321.5+i0.1321.5-i0.132
将表2、3和表1计算的应力奇性指数对比发现，两相材料黏结(连续性条件)计算的应力奇性指数λk(见表1)囊括了域Ω1(刚性小)和域Ω2(刚性大)计算的应力奇性指数λk(见表2、3)．
表2、3还显示了域Ω1和域Ω2计算的应力奇性指数各自占优，因此对域Ω1和域Ω2应力场的求解要选用不同的应力特征对
下面分为4种选择方案，分别采用扩展边界元法进行计算，然后对计算结果进行分析对比得出较优方案．
方案1 域Ω1(刚性小)和域Ω2(刚性大)均采用理想结合黏结材料的应力特征对；
方案2 域Ω1(刚性小)采用理想结合黏结材料的应力特征对，域Ω2(刚性大)抛去理想结合黏结材料的第4阶应力特征对，后续项全使用；
方案3 域Ω1(刚性小)采用域Ω1裂纹尖端的应力特征对，域Ω2(刚性大)采用理想结合黏结材料的应力特征对；
方案4 域Ω1(刚性小)采用理想结合黏结材料的应力特征对，域Ω2(刚性大)抛去理想结合黏结材料的第1～3阶应力特征对，后续项全使用．
采用扩展边界元法分别对图3裂纹结构按照4种方案计算全域的应力场，这里取图3尖端处小扇形半径ρ=0.018 mm，扩展边界元法对ρ值不敏感[19]，划分了336个边界节点和168个二次边界单元，获得式(1)中的位移幅值系数Ak，见表4．
观察4种方案的位移幅值系数计算值，方案3和4出现过大的幅值系数，且其过渡不平稳．方案1和2未出现过大的幅值系数，且其过渡较平稳．因此，方案3和4的计算结果不准确，其应力特征对选择不合理，但无法判断方案1和2的优劣，需对方案1和2的计算结果进行进一步的研究．
采用方案1和2对两相材料裂纹结构E2/E1=100,10,5,2进行计算，计算所得应力强度因子KⅠ值见表5．
表4 扩展边界元法分析两相材料裂纹的位移幅值系数Ak/mm-λk
Table 4 Amplitude coefficients Ak/mm-λk of the crack from the XBEM analysis
kλkRλkIscheme 1AkRAkIscheme 2AkRAkIscheme 3AkRAkIscheme 4AkRAkI1-10-8.6E-407.28E-40-4.49E+602.83E-302-102.71E-40-7.6E-40-2E+70-2.5E-303-0.4990.132-7.5E-4-1.1E-33.71E-4-3.92E-3-6.2E+5-7.03E+6-1.08E-2-7.87E-14002.06E-50-2.15E-20-1.27E+906.48E-20500-3.85E-40-6.13E-201.98E+906.53E-1060.50.1323.86E-5-4.56E-5-1.4E-1-2.24E-1-1.6E-2-1.6E-4-44.293-231.07370.99904.23E-30-2.68203.05E-30-3.64908103.05E-203.1910-1.83E-2047.6110
表5 扩展边界元法分析两相材料裂纹的应力强度因子KⅠ/(N · mm-2-λ3)
Table 5 Stress intensity factors KⅠ/(N · mm-2-λ3) of the crack from the XBEM analysis
E2/E1XBEMscheme 1scheme 2ref. [20]ref. [18]2650.9-1 641.6658.9658.95640.5-634.1655.1-10701.2635.9644.9643.8100740.6627.4623.7629.8
分析表5的数据可知，当E2/E1=2，5时，方案1计算所得的KⅠ和文献[18,20]较方案2接近；当E2/E1=10，100时，方案2计算所得的KⅠ和文献[18，20]较方案1接近．因此，初步断定，方案1较方案2更适合计算E2/E1=2，5时裂纹结构全域的应力场，方案2较方案1更适合计算E2/E1=10，100时裂纹结构全域的应力场．
为进一步证明上述结论的准确性，采用有限元法对裂纹尖端划分非常细的网格进行应力场计算，并将其计算结果和方案1、2进行分析对比．方案1、2和有限元法计算两相裂纹结构E2/E1=2，5，10，100时，尖端O附近区域沿y=0 mm的界面σy值如图4所示．
(a) E2/E1=2(b) E2/E1=5
(c) E2/E1=10(d) E2/E1=100
图4 沿y=0 mm的界面σy应力
Fig. 4 Variation of interface stress σy along y=0 mm
扩展边界元法划分了168个二次边界单元，其解可以很好地描述尖端根部的应力场．有限元法采用四边形4节点等参元，划分了400 000个单元，在尖端O附近单元很密，其使用的四边形4节点单元数是扩展边界元法二次单元数的2 380倍，极大地增加了计算机计算量．
观察图4可知：
1) 在y=0 mm，x<14.0 mm的尖端区域，由于尖端O附近出现强应力集中，有限元解已失真．随着距尖端距离的增大，有限元法与扩展边界元法计算相同点的应力值差值减小．
2) 在y=0 mm，x<0.02 mm的范围内，方案1和2计算的界面σy值相差较大，因为方案2对域Ω2(刚性大)不考虑第1阶应力特征对．在y=0 mm，x>0.05 mm的渐远区域方案1和2的计算值趋同，说明该范围区域，方案1和2的解均是准确的．
3) 图4(a)、4(b)显示在y=0 mm，x=0.005 mm时，方案2计算的界面σy值有突变，而方案1计算的界面σy值连续，进一步证明了方案1较方案2更适合计算双相材料E2/E1=2，5时Ⅰ型裂纹结构全域的应力场．
4) 图4(c)、4(d)显示在y=0 mm，x=0.02 mm时，方案1计算的界面σy值有突变，而方案2计算的界面σy值连续，进一步证明了方案2较方案1更适合计算双相材料E2/E1=10，100时Ⅰ型裂纹结构全域的应力场．
基于以上结论，对图3结构施加拉伸面力P=100 MPa，剪切面力τ=100 MPa，域Ω1采用硬铝合金材料，E1=7.0×104MPa，ν1=0.3，域Ω2采用合金钢材料，E2=2.06×105MPa，ν2=0.3，两种材料的E2/E1=2.94，为相当情形．采用方案1和有限元法计算其完整的应力场，在尖端O附近区域沿y=0 mm的σy应力和τxy应力见表6．
表6 钢铝结合材料裂纹尖端附近内点应力σy和τxy
Table 6 Stresses σy and τxy near the crack tip with the bonded steel and aluminum
inner point coordinates(x,y)/mmσy/MPaXBEM scheme 1FEMτxy/MPaXBEM scheme 1FEM(0.002,0)4 530.43 597.3(0.01,0)1 913.71 880.9(0.234,0)546.1493.6(0.99,0)259.6481.7268.4552.4(3.19,0)153295.6170.3333.4(14,0)104.8106118.4123.9(61.9,0)98.298.2103.5105(130,0)95.596.896.698.4(-4.8,5)9.112.876.3164.9(4.5,5)140.9156106.3112.5(33.7,5)100.2102.5108.7109.4(103.2,5)98.1100.599.1101.3(-4.8,-5)50.4100.3117224.7(4.5,-5)145.5154.48198.6(33.7,-5)98.1100.1105.2107.2(103.2,-5)95.498.69899
采用扩展边界元法对图3裂纹结构按照方案1计算全域的应力场，这里取图3尖端O处小扇形半径ρ=0.018 mm．观察表6中的第2～第4列中的数值可知：沿y=0 mm，扩展边界元法可以获得离尖端O极近区域(如x≤0.005 mm)的应力场；而有限元法无法获得尖端根部(x<0.99 mm)区域的应力场；只有当x≥14 mm时，有限元解方为有效，体现了扩展边界元法求解尖端应力场的优越性．
域Ω2(合金钢)内点(-4.8 mm，-5 mm)附近的σy值约为域Ω1(合金铝)内点(-4.8 mm，5 mm)的σy值的5.54倍；域Ω2内点(-4.8 mm，-5 mm)附近的τxy值约为域Ω1内点(-4.8 mm，5 mm)的τxy值的1.53倍，符合力学规律．因此方案1的扩展边界元法可以准确获得两相材料复合型裂纹结构E2/E1=2.94比值相当情形的全域应力场．本文对两相材料裂纹结构全域的应力场结果可评估裂纹体破坏行为．
4 结 论
本文采用扩展边界元法分析了两相材料裂纹结构全区域应力场，特别是尖端区域奇异应力场，主要结论如下：
1) 求解裂纹全域的应力场时，扩展边界元法可使用较少的单元数准确获得全域的应力场，而有限元法不能获得尖端根部区域的应力场．利用尖端区域位移渐近展开式求解两相材料裂纹结构尖端区域奇异应力场时，根据两者弹性模量比值的不同，应对两材料区域分别取用适配的位移渐近展开式项数．
2) 针对E2/E1比值的极大情形(如E2/E1=100),域Ω1(刚性小)采用理想结合黏结材料的应力特征对，域Ω2(刚性大)舍去理想结合黏结材料的第4阶应力奇异项，后续若干高阶项全使用(即方案2)．然后用扩展边界元法求解，可以得到全域的应力场．
3) 针对E2/E1比值的较大情形(如E2/E1=10，5),域Ω1(刚性小)采用理想结合黏结材料的特征对，当E2/E1=10比值靠近极大情形时，域Ω2(刚性大)抛去理想结合黏结材料的第4阶应力奇异项，后续若干高阶项全使用(即方案2)；当E2/E1=5比值靠近相当情形时，域Ω2和域Ω1均采用理想结合黏结材料的应力特征对(即方案1)．然后用扩展边界元法求解，可以得到全域的应力场．
4) 针对E2/E1比值的相当情形(如E2/E1=2),域Ω1和域Ω2均采用理想结合黏结材料的应力特征对(即方案2)．然后用扩展边界元法求解，可以准确得到全域的应力场．
5) 对E2/E1任何比值，域Ω1和域Ω2采用方案1和2求解的离尖端较远区域应力场均准确．而有限元法在尖端附近采用非常密的网格求出的应力值仅在距离尖端为裂纹长度1～2倍的远场方为有效．
[1] 沈观林, 胡更开. 复合材料力学[M]. 北京: 清华大学出版社, 2006: 3-5.(SHEN Guanlin, HU Gengkai. Mechanics of Composites[M]. Beijing: Tsinghua University Press, 2006: 3-5.(in Chinese))
[2] 葛仁余, 牛忠荣, 程长征, 等. 边界元法分析二维线弹性裂纹扩展[J]. 计算物理, 2015, 32(3): 310-320.(GE Renyu, NIU Zhongrong, CHENG Changzheng, et al. Propagation analysis of two-dimensional linear elastic crack with boundary element method[J]. Chinese Journal of Computational Physics, 2015, 32(3): 310-320.(in Chinese))
[3] 王振, 余天堂. 模拟三维裂纹问题的自适应多尺度扩展有限元法[J]. 工程力学, 2016, 33(1): 32-38.(WANG Zhen, YU Tiantang. Adaptive multiscale extended finite element method for modeling three-dimensional crack problems[J]. Engineering Mechanics, 2016, 33(1): 32-38.(in Chinese))
[4] 秦洪远, 黄丹, 刘一鸣, 等. 基于改进型近场动力学方法的多裂纹扩展分析[J]. 工程力学, 2017, 34(12): 31-38.(QIN Hongyuan, HUANG Dan, LIU Yiming, et al. An extended peridynamic approach for analysis of multiple crack growth[J]. Engineering Mechanics, 2017, 34(12): 31-38.(in Chinese))
[5] 江守燕, 杜成斌, 顾冲时, 等. 求解双材料界面裂纹应力强度因子的扩展有限元法[J]. 工程力学, 2015, 32(3): 22-27, 40.(JIANG Shouyan, DU Chengbin, GU Chongshi, et al. Computation of stress intensity factors for interface cracks between two dissimilar materials using extended finite element methods[J]. Engineering Mechanics, 2015, 32(3): 22-27, 40.(in Chinese))
[6] 吕君, 柴国钟. 双材料裂纹问题的积分方程方法[J]. 浙江工业大学学报, 2017, 45(2): 130-136, 158.(LÜ Jun, CHAI Guozhong. Integral equation methods for the problem of bi-material crack[J]. Journal of Zhejiang University of Technology, 2017, 45(2): 130-136, 158.(in Chinese))
[7] 牛忠荣, 葛仁余, RECHO N, 等. 平面V形切口塑性应力奇异性分析[J]. 中国科学: 物理学 力学 天文学, 2014, 44(1): 79-90.(NIU Zhongrong, GE Renyu, RECHO N, et al. Evaluation of plastic stress singularities of plane V-notches and cracks in hardening materials[J]. Scientia Sinica: Physica, Mechanica & Astronomica, 2014, 44(1): 79-90.(in Chinese))
[8] WILLIAMS M L. Stress singularities resulting from various boundary conditions in angular corners of plates in extension[J]. Journal of Applied Mechanics, 1952, 19(1): 287-298.
[9] 李有堂, 王勇. 功能梯度材料V型缺口根部裂纹场强特性分析[J]. 兰州理工大学学报, 2018, 44(4): 162-166.(LI Youtang, WANG Yong. Characteristic analysis of crack stress field intensity at V-shaped notch root of functionally gradient materials[J]. Journal of Lanzhou University of Technology, 2018, 44(4): 162-166.(in Chinese))
[10] MIRSAYAR MM, ALIHA M R M, SAMAEI A T. On fracture initiation angle near bi-material notches-effects of first non-singular stress term[J]. Engineering Fracture Mechanics, 2014, 119: 124-131.
[11] 李俊林, 张少琴, 杨维阳, 等. 正交异性双材料界面裂纹尖端应力场[J]. 应用数学和力学, 2008, 29(8): 947-953.(LI Junlin, ZHANG Shaoqin, YANG Weiyang, et al. Study of stress fields near interface crack tip of double dissimilar orthotropic composite materials[J]. Applied Mathematics and Mechanics, 2008, 29(8): 947-953.(in Chinese))
[12] AYATOLLAHI M R, MIRSAYAR M M, NEJATI M. Evaluation of first non-singular stress term in bi-material notches[J]. Computational Materials Science, 2010, 50(2): 752-760.
[13] 范海军, 肖盛燮. Ⅱ型平面应力裂纹线场弹塑性极坐标精确解[J]. 应用力学学报, 2015, 32(4): 543-548, 701.(FAN Haijun, XIAO Shengxie. Elastic-plastic exact solution of polar coordinates on mode Ⅱ plane stress cracking in stress field[J]. Chinese Journal of Applied Mechanics, 2015, 32(4): 543-548, 701.(in Chinese))
[14] LAN X, NODA N A, MITHINAKA K, et al. The effect of material combinations and relative crack size to the stress intensity factors at the crack tip of a bi-material bonded strip[J]. Engineering Fracture Mechanics, 2011, 78(14): 2572-2584.
[15] 牛忠荣, 程长征, 胡宗军, 等. V形切口应力强度因子的一种边界元分析方法[J]. 力学学报, 2008, 40(6): 849-857.(NIU Zhongrong, CHENG Changzheng, HU Zongjun, et al. Boundary element analysis of the stress intensity factors for V-notched Structures[J]. Chinese Journal of Theoretical and applied Mechanics, 2008, 40(6): 849-857.(in Chinese))
[16] 王有成. 工程中的边界元方法[M]. 北京: 中国水利水电出版社, 1996.(WANG Youcheng. Boundary Element Method in Engineering[M]. Beijing: China Water & Power Press, 1996.(in Chinese))
[17] 高效伟, 郑保敬, 刘健. 功能梯度材料动态断裂力学的径向积分边界元法[J]. 力学学报, 2015, 47(5): 868-873.(GAO Xiaowei, ZHENG Baojing, LIU Jian. Dynamic fracture analysis of functionally graded materials by radial integration BEM[J]. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(5): 868-873.(in Chinese))
[18] YUUKI R, CHAO S B. Efficient boundary element analysis of stress intensity factors for interface cracks in dissimilar materials[J]. Engineering Fracture Mechanics, 1989, 34(1): 179-188.
[19] 李聪, 牛忠荣, 胡斌, 等. 扩展边界元法分析切口和裂纹结构应力场的准确性[J]. 固体力学学报, 2018, 39(5): 539-551.(LI Cong, NIU Zhongrong, HU Bin, et al. Accuracy of the extended boundary element method analyzing the stress fields of V-notched/cracked structures[J]. Chinese Journal of Solid Mechanics, 2018, 39(5): 539-551.(in Chinese))
[20] 张明, 姚振汉, 杜庆华. 双材料界面裂纹应力强度因子的边界元分析[J]. 应用力学学报, 1999, 16(1): 21-26.(ZHANG Ming, YAO Zhenhan, DU Qinghua. Boundary element analysis of stress intensity factors of bimaterial interface crack[J]. Chinese Journal of Applied Mechanics, 1999, 16(1): 21-26.(in Chinese))