钟　红，宋平平

(大连理工大学 建设工程学部，110624 辽宁 大连)



任意裂纹面荷载作用下界面断裂分析

钟红，宋平平

(大连理工大学 建设工程学部，110624 辽宁 大连)

摘要:为研究裂纹面上作用的荷载对裂纹稳定性的影响，本文基于比例边界有限元方法提出裂纹面作用有任意方向、任意大小面荷载的界面应力强度因子求解模型.界面裂纹具有复数形式的应力奇异性指数，在任意裂纹面荷载作用下其奇异应力场更为复杂.应用本模型，径向的位移和应力可解析求解，无需网格细分即可自动反映裂尖的应力奇异性.裂纹面上的任意荷载首先可分解成平行于裂纹面以及垂直于裂纹面的分量，并进一步分解成有限项幂函数的和.对每个幂函数荷载解析求解，基于线性叠加原理获得结构在全部荷载作用下的解.该模型对各向同性材料和各向异性材料均适用.文中通过板承受裂缝面荷载时的应力强度因子求解的多个算例对该模型进行了验证和应用，对板的几何尺寸和双材料参数进行了敏感性分析，并应用于重力坝坝踵界面裂缝在水压力作用下的应力强度因子求解.

关键词:比例边界有限元；界面断裂；任意裂纹面荷载；应力强度因子

很多实际工程都涉及界面，例如岩基上的混凝土坝，坝基和混凝土的交界面通常比较薄弱，易于出现裂纹，尤其是当水进入缝内后，水压力将对裂缝的进一步扩展起促进作用，从而劣化大坝的稳定性.在日常生活中，也常常遇到界面断裂问题，诸如焊接、粘接等结合材料，通常在结合处或者其附近首先开裂.这是因为结合材料界面附近不仅容易存在缺陷，导致结合强度的低下，而且会因界面的存在而引发应力集中并产生残余应力，使界面附近的材料处于较高的应力水平.随着复合材料应用范围的扩大，界面问题变得越来越重要，传统的强度分析和评价方法局限性也日益明显.

不同于均质材料断裂，界面断裂有一些特殊性.Williams[1]分析了界面裂纹尖端的奇异场，利用应力函数的分离变量形式，求得奇异性指数和奇异应力场，但是该奇异性指数(0.5±iε)不是实数而是复数，导致了裂纹尖端应力场的振荡奇异性和裂纹面的互相嵌入.振荡引起I型断裂和II型断裂耦合，对称结构内的裂缝即使处于对称荷载作用下，其断裂也是复合型的.常见的断裂力学求解方法，包括有限元法、边界元法、边界配置法[2]和扩展有限元[3]等，所采用的标准插值函数都是光滑的，与奇异应力场相差甚远[4].有限元法求解断裂问题时，为了得到更精确的应力解，划分有限元网格时需要在裂尖局部加密或引入奇异单元(如四分之一节点单元[5-6])进行求解.然而对于界面断裂问题，奇异应力场的近似解是非常复杂的，对单元进行改进的复杂程度远远大于求解断裂问题本身.Miyazaki等[7]提出M1积分方法求解双材料界面断裂问题，分别计算了含单边裂纹和中心斜裂纹双材料板的应力强度因子;Munz等[8]基于有限元方法描述了双材料界面裂纹处的应力分布特征;陈瑛等[9]综合评述和分析了多种断裂力学模型和实验方法，同时介绍了双材料界面断裂力学在FRP-混凝土复合结构中的应用.

界面裂缝的缝面荷载对裂缝的稳定性有至关重要的影响.在这种情况下，裂尖的奇异应力场和应力强度因子都将产生显著变化，从而对数值方法和数值模型提出了新的挑战.对于裂纹面上承受任意荷载的复杂情况研究较少，其中胡小飞[10]采用基于辛体系的解析奇异单元分析含裂纹的结构；刘钧玉[11]基于比例边界有限元法计算了一类面荷载作用下的裂缝奇异应力场；涂传林[2]利用边界元法研究了裂纹面上受均匀法向外荷载的断裂问题.以上研究所考虑的荷载形式和作用方向均较简单.

比例边界有限元法(scaled boundary finite element method，SBFEM)是一种新型的半解析数值方法，可计算多种材料交界面处[4]的奇异应力场，以及温度荷载[12]、动荷载[13]等作用下的奇异应力场，并已推广至非线性断裂模拟[14].本文采用比例边界有限元法，基于裂纹面荷载的幂级数展开和线性叠加原理，提出了求解任意裂纹面荷载作用下的界面断裂求解模型.将该模型应用于各向同性和各向异性双材料板的界面问题，通过与文献结果对比进行了验证.在此基础上开展了一定的参数敏感性分析.

1比例边界有限元方法的基本原理

整体坐标系下一点的坐标用比例边界有限元坐标表示为

(1)

图1　比例边界有限元模型和ξ-η坐标

(2)

应力为

(3)

式中[D]是材料的弹性矩阵，B1(η)和B2(η)是应变位移矩阵，参见文献[15].用位移表达的比例边界有限元方法的控制方程为

(4)

(5)

(6)

式中μ是个很小的数(如0.000 1)，径向的内部节点力为[17]

(7)

方程(4)可写成一阶常微分方程

(8)

式中[Z]是Hamiltonian系数矩阵[17]，特征值为λi和-λi，方程(8)可以通过[Z]阵特征向量进行解耦，运算中容易出现数值不稳定，产生对数奇异.

本文采用块对角Schur分解[17]

(9)

裂纹面荷载的节点位移模态为

(10)

相应的等效节点力为

(11)

则位移解为

(12)

对于给定的积分常数，边界上的节点位移为

(13)

则对应的等效边界节点力为

(14)

由式(13)可得积分常数用边界位移表达为

(15)

将方程(15)代入方程(14)得

(16)

式中[K]为刚度矩阵.通过边界条件，由式(16)解出边界节点位移{ub}，代入式(15)求得积分常数{c}，位移场由式(12)求出.求得的位移场代入式(3)，最后求出应力场：

(17)

上式可整理写成

(18)

其中

(19)

(20)

2应力强度因子

(21)

其中μi是剪切模量

对于各向同性双材料板，标准应力强度因子定义为

(23)

式中L为特征长度.方程(23)可表示为矩阵形式:

(24)

其中

(25)

对于各向异性双材料板，应力强度因子可定义为

(26)

其中W1，W2可由各向异性材料的弹性常数计算得出[18].本文采用的广义应力强度因子通过推导可表示为[15]

(27)

3数值算例

给出4个带裂缝平板的应力强度因子，考虑了各向同性和各向异性材料，裂纹面荷载考虑了法向和切向荷载.执行计算工作的计算机配置为：处理器Intel(R)Core(TM)i5-2300 CPU @ 2.80 GHz，4个内核，4个逻辑处理器，物理内存8.00 GB.

3.1各向同性单边裂纹单材料板承受法向裂纹面荷载

图2　单边裂纹单材料板

表1同时给出了半无限大板单边裂纹承受任意荷载的I型应力强度因子解析解[19]，作为本文的参考解.从结果可以看出，随着荷载指数n的增大，应力强度因子减小，随着板的尺寸W/a增大，应力强度因子减小，板的尺寸W/a足够大时，可以近似用来模拟单边裂纹半无限大板.当W/a=30，网格划分N=10时的计算结果与解析解很接近，表格最后一行给出了W/a=30，N=10计算结果与解析解之间的误差，误差范围均小于3%.图3给出了板尺寸W/a=30时，不同网格下的计算结果与解析解的对比，可以看出误差很小，粗细不同的3种网格计算结果相差较小，由此说明本方法的计算精度对网格粗细划分不敏感，较少的网格就可以达到计算精度.N=2时1.4 s即可完成整个计算过程，N=10时60 s完成计算过程.

表1各向同性单边裂纹单材料板承受法向裂纹面荷载的无量纲应力强度因子

W/a网格疏密n=0n=1n=2n=3N=21.3678820.5875600.3881170.2917255N=61.3678830.5875600.3881000.291690N=101.3678820.5875580.3881050.291812N=21.1892650.4792680.3102180.23250410N=61.1893290.4792920.3102720.231388N=101.1893280.4792890.3102250.230874N=21.1390520.4493770.2889690.21667820N=61.1398970.4496310.2897310.216867N=101.1398980.4496120.2880790.214919N=21.1285140.4438040.2905670.23177930N=61.1299020.4436630.2871290.223066N=101.1299470.4436840.2850070.203025解析解[19]1.1214710.4385480.2811460.208314误差/%0.75 1.17 1.37 2.54

3.2各向同性单边裂纹双材料板承受法向裂纹面荷载

图3　不同网格计算结果与解析解对比

图4　单边裂纹双材料板

表2各向同性单边裂纹双材料板承受法向裂纹面荷载的无量纲应力强度因子

η应力强度因子n=1n=3K*IK*IIK*IK*II1本文0.4434760 0.1950100 文献[10]0.4362640 0.2078400 2本文0.443854-0.0064820.217686-0.007837文献[10]0.438484-0.0049010.209213-0.0082825本文0.444680-0.0157390.206461-0.022349文献[10]0.438604-0.0127020.209820-0.01861310本文0.444920-0.0219410.213263-0.025749文献[10]0.438746-0.0182920.209939-0.024628

3.3正交各向异性单边裂纹双材料板承受法向裂纹面荷载

本文求得的应力强度因子见表3，无解析解可与之对比.可以看出虽然几何图形是对称的，裂纹面只承受对称法向荷载的作用，但是由于材料1和材料2的差异性会产生II型应力强度因子，并且随着指数n的增大，相应的I型和II型应力强度因子减小.当φ1=0°,90°时，材料1为正交各向异性材料，Ⅰ型和Ⅱ型应力强度因子绝对值小于φ1=30°、60°时应力强度因子的绝对值.原因是当材料为正交各向异性材料时式(26)中W2为0.

图5　承受法向裂纹面荷载的单边裂纹板

3.4正交各向异性单边裂纹双材料板承受切向和法向面荷载

对于裂纹表面既承受法向裂纹面荷载σ，又承受剪切荷载τ的问题，I型和II型应力强度因子不仅与法向荷载σ有关，还与剪切荷载τ有关.同时I型、II型应力强度因子也与两种材料的弹性常数有关，两种材料之间有干涉作用，4组材料中材料1的性能不变，材料2的y方向弹性模量在变化，随着y方向弹性模量的减小，I型应力强度因子增大，II型应力强度因子减小.应用本文界面断裂求解模型整个计算过程不超过5s.

表3正交各向异性单边裂纹双材料板承受法向裂纹面荷载的无量纲应力强度因子

φ1/(°)应力强度因子n=0n=1n=2n=30K*I1.3767080.6231680.4149620.312437K*II-0.219304-0.044082-0.021666-0.01360230K*I1.5299580.6641310.4380710.330323K*II-0.417664-0.147609-0.095482-0.07151760K*I1.6013020.6912450.4556050.342079K*II-0.654783-0.244728-0.160870-0.12118190K*I1.4892040.6510890.4299860.322961K*II-0.723042-0.240713-0.152680-0.112912

图6　承受法向和切向裂纹面荷载的单边裂纹板

3.5重力坝算例

表5给出了坝体和地基的不同模量比、不同缝面水压力分布时的应力强度因子.可看出对于不同水压力分布形式，随着坝体和坝基弹模比值的增大，KI均明显减小；对于给定弹模，λ=0时KI最小，随着n的减小，施加的缝内水压增大，KI增大.当坝体地基模量比较小时，水压力的差异对KI的影响更为重要，随着模量比增大，界面断裂的耦合效应影响加大.KII的大小主要取决于上游面水压力，因此受缝内水压分布影响不大，但当坝体地基模量比增大时，界面断裂耦合效应使得KII有所增大.由于此时KI减小，KII/KI呈增大趋势，裂尖剪切分量增大.

图7　重力坝尺寸(m)　图8　子域划分与网格剖分

表4正交各向异性单边裂纹双材料板承受切向和法向面荷载的无量纲应力强度因子

(E22/E11)II应力强度因子n=0n=1n=2n=31K*I1.1259920.5209490.3476320.262732K*II1.4748270.5172800.3214450.2344610.5K*I1.1961750.5394620.3587740.270020K*II1.4090670.5049300.3163230.2314790.3K*I1.2544330.5560230.3691320.278416K*II1.3453490.4927320.3108750.2285020.1K*I1.3864110.5988640.3959580.297984K*II1.1487650.4546000.2924090.216945

表5坝基界面裂缝的应力强度因子106N·m-3/2

E1/E2应力强度因子λ=0n=0n=1n=21KI1.2963811.6594831.3830281.335719KII1.4076561.4614531.4393591.4076562KI1.0101231.5044261.1356531.071478KII1.4223281.3930041.4260531.4272305KI0.6112961.2535520.7959350.712791KII1.5104721.3493691.4732531.48526510KI0.3471891.0668080.5697720.476663KII1.6038411.3603831.5400061.559327

4结语

基于比例边界有限元方法提出了裂纹面作用有任意方向、任意大小面荷载的界面断裂求解模型.首先给出了比例边界有限元方法的基本方程，针对任意裂纹面荷载问题，将荷载分解成平行于裂纹面以及垂直于裂纹面的分量，并各自分解成有限项幂函数的和，对每个幂函数荷载解析求解，基于线性叠加原理获得结构在全部荷载作用下的解.第一个算例单材料板的计算结果与解析解进行对比，验证了本模型有较高的计算精度和计算效率，网格剖分简单.接着3个算例双材料界面断裂问题，研究了几何尺寸和材料参数的变化对KI和KII的影响，本文计算模型可用于求解各向同性和各向异性双材料界面断裂问题.最后将本模型应用于重力坝坝踵裂缝承受水压力时的应力强度因子求解，发现缝内水压分布形式对KI影响较大；随着坝体和地基模量比的增大，KI明显减小，KII有所降低，裂尖的剪切分量比重增大，断裂模态复合的程度加剧.

参考文献

[1] WILLIAMS M L. The stresses around a fault or crack in dissimilar media[J]. Bulletin of the Seismological Society of America, 1959, 49(2): 199-204.

[2] 涂传林.裂缝面上受外荷载作用下的边界配置法及其应用[J].水利学报,1983,7: 9-16.

[3] BECHET E, MINNEBO H, MOES N, et al. Improved implementation and robustness study of the X-FEM for stress analysis around cracks[J]. International Journal for Numerical Methods in Engineering, 2005, 64(8): 1033-1056.

[4] SONG C, WOLF J P. Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multi-materials with the scaled boundary finite-element method[J]. Computers & Structures, 2002, 80(2): 183-197.

[5] BARSOUM R S. Application of quadratic isoparametric finite elements in linear fracture mechanics[J]. International Journal of Fracture, 1974, 10(4): 603-605.

[6] HENSHELL R D, SHAW K G. Crack tip finite elements are unnecessary[J]. International Journal for Numerical Methods in Engineering, 1975, 9(3): 495-507.

[7] MIYAZAKI N, IKEDA T, SODA T, et al. Stress intensity factor analysis of interface crack using boundary element method-application of contour-integral method[J]. Engineering Fracture Mechanics, 1993, 45(5): 599-610.

[8] MUNZ D, YANG Y Y. Stresses near the edge of bonded dissimilar materials described by two stress intensity factors[J]. International Journal of Fracture, 1993, 60(2): 169-177.

[9] 陈瑛,乔丕忠,姜弘道,等. 双材料界面断裂力学模型与实验方法[J].力学进展,2008,38(1):53-61.

[10]胡小飞.基于辛空间的解析奇异单元及其在断裂力学中的应用[D].大连:大连理工大学,2012.

[11]刘钧玉.裂纹内水压对重力坝断裂特性影响的研究[D]. 大连:大连理工大学，2008.

[12]SONG C. Analysis of singular stress fields at multi-material corners under thermal loading[J]. International journal for numerical methods in engineering, 2006, 65(5): 620-652.

[13]SONG C. A super-element for crack analysis in the time domain[J]. International journal for numerical methods in engineering, 2004, 61(8): 1332-1357.

[14]YANG Z J, DEEKS A J. Fully-automatic modelling of cohesive crack growth using a finite element-scaled boundary finite element coupled method[J]. Engineering Fracture Mechanics, 2007, 74(16): 2547-2573.

[15]WOLF J P, SONG C. The scaled boundary finite-element method-a primer: derivations[J]. Computers & Structures, 2000, 78(1): 191-210.

[16]SONG C, WOLF J P. The scaled boundary finite-element method-a primer: solution procedures[J]. Computers & Structures, 2000, 78(1): 211-225.

[17]SONG C. A matrix function solution for the scaled boundary finite-element equation in statics[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(23): 2325-2356.

[18]SONG C, TIN-LOI F, GAO W. A definition and evaluation procedure of generalized stress intensity factors at cracks and multi-material wedges[J]. Engineering Fracture Mechanics, 2010, 77(12): 2316-2336.

[19]丁遂栋.断裂力学[M].北京：机械工业出版社.1997.

(编辑赵丽莹)

Analysis of interface crack with arbitrary crack tractions

ZHONG Hong, SONG Pingping

(Faculty of Infrastructure Engineering, Dalian University of Technology, 110624 Dalian, Liaoning，China)

Abstract:This paper presents a model for solving interface crack with arbitrary crack traction based on the Scaled boundary finite element method to study the significant influence of the traction acting on the crack faces on the stability of a crack. The order of stress singularity is complex for an interface crack. With the existence of crack traction, the stress singularity is more complicated. Base on the proposed model, stress and displacement are solved analytically in the radial direction, and the stress singularity at crack tip is obtained with high precision without refined mesh. The arbitrary crack traction is firstly decomposed to one component parallel to the crack and the other one perpendicular to the crack, then both the two components are expressed as the sum of a limited number of power functions respectively. The effect of each power function is solved analytically. According to the Linear superposition principle, the solution of a structure with arbitrary crack traction can be obtained. The proposed model is effective for both anisotropic and isotropic materials. The model is verified by several plates with crack tractions, in which stress intensity factors are calculated. Sensitivity analysis is also performed concerning the plate geometry and material properties. Finally the model is applied to solve the stress intensity factors of an interface crack of a gravity dam filled with water.

Keywords:the scaled boundary finite element method; interface crack; arbitrary crack traction; stress intensity factor

中图分类号:TU311.1

文献标志码:A

文章编号:0367-6234(2016)02-0152-06

通信作者:钟红，hzhong@dlut.edu.cn.

作者简介:钟红(1981—)，女，副教授，硕士生导师.

基金项目:国家自然科学基金(51009019,51579033)；中央高校基本科研业务费专项资金(DUT14LK40)；中国博士后基金特别资助项目(2013T60283).

收稿日期:2014-11-08.

doi：10.11918/j.issn.0367-6234.2016.02.026