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

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

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

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

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

比例边界有限元控制方程的推导和求解见文献[15-17].本文主要介绍对任意裂纹面荷载的处理和奇异应力场的求解.图 1是具有V形裂纹口的比例边界有限元模型，O为比例中心，断裂问题比例中心通常选在裂尖处.定义ξ(0≤ξ≤1)为径向坐标，模型边界离散成一维线单元，ξ－η形成比例边界有限元坐标.

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

(1)

式中{N(η)}为形函数，{x}、{y}是节点坐标.{u(ξ)}，边界上的节点位移为{u}={u(ξ=1)}，则区域内任意一点位移函数可表示为

(2)

应力为

(3)

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

(4)

式中[E¹]、[E²]和[E³]为系数矩阵^[15]，各单元系数矩阵的计算和组装与有限元法类似.{F(ξ)}是外部荷载向量，包括裂纹面荷载{F_t(ξ)}.裂纹面上的任意荷载可以分解成有限项幂函数和的形式

(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)解出边界节点位移{u_b}，代入式(15)求得积分常数{c}，位移场由式(12)求出.求得的位移场代入式(3)，最后求出应力场：

(17)

上式可整理写成

(18)

其中

(19)

为应力模态，表示为

(20)

2 应力强度因子

应力强度因子通过极坐标下的奇异应力和来定义.对于双材料界面断裂问题，具有振荡奇异性，应力奇异性指数是一对共轭复数0.5±iε，ε为振荡指数，大小取决于两种结合材料参数的差异.

(21)

其中μ_i是剪切模量

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

(23)

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

(24)

其中

(25)

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

(26)

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

(27)

式中为奇异应力模态，是奇异应力项对应的积分常数.

3 数值算例

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

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

考虑如图 2所示的含单边裂纹平板，板的尺寸是W×2W，裂纹长度是a，比例中心为O.裂纹表面受到法向裂纹面荷载σ₀=σ₀(r)的作用，σ₀=λ(r/a)ⁿ，其中n=1, 2, …, n; n是任意常数，N表示W边划分的线单元个数，图 2给出N=6的网格图，线单元采用具有Gauss-Lobatto-Legendre形函数的11节点高阶单元.表 1给出不同板长，不同网格下(N=2, 6, 10)无量纲化应力强度因子的数值计算结果.

表 1同时给出了半无限大板单边裂纹承受任意荷载的Ⅰ型应力强度因子解析解^[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完成计算过程.

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

考虑如图 4所示的单边裂纹各向同性双材料板，界面处有一长度为a的裂纹，板的尺寸是W×2W，比例中心为O，裂纹表面作用有法向裂纹面荷载σ₀=λ(r/a)ⁿ，其中n=1, 2, …, n; n是任意常数.虽然几何对称、荷载对称，由于材料不对称，本题是个复合断裂问题.下面给出了无量纲化应力强度因子E₂/E₁，v=0.3, n取1和3，比例中心O选在裂尖，线单元采用具有Gauss-Lobatto-Legendre形函数的11节点高阶单元, W边划分4个线单元，网格剖分见图 4，计算结果见表 2.

胡小飞^[10]利用辛对偶体系构造高阶精度解析奇异单元，与常规有限元单元相结合求解应力强度因子.本文方法的计算结果与其计算结果很相近，随着η=E₂/E₁的增大，K_Ⅱ^*也相应增大，这是由于两种材料的不匹配加剧造成的.

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

如图 5所示单边裂纹双材料板，界面处有一长度为a的裂纹，板的尺寸是W×2W，W/a=5，比例中心为O，裂纹面上承受σ₀(r)=λ(r/a)ⁿ(n=1, 2, 3)的法向裂纹面荷载，材料1的属性为E₁₁=200 MPa，E₂₂=200 MPa, ν₁₂=0.4，G₁₂=29.41 MPa，材料2的属性E₁₁=10 MPa，E₂₂=100 MPa，ν₁₂=0.02，G₁₂=28.07 MPa，φ₁与φ₂是材料主轴与正方向的夹角.计算φ₂=0，不同φ₁下的应力强度因子.其中K_I^*=，W边划分4个线单元，网格划分见图 5.

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

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

如图 6所示的单边裂纹双材料板，位于界面上的裂纹a，板的尺寸是W×2W，W/a=5，比例中心为O，裂纹面上承受荷载σ₀(r)=λ(r/a)ⁿ，τ₀r=λ(r/a)ⁿ(n=1, 2, 3)，φ为材料弹性主方向与轴正向夹角，材料1的弹性参数E₁₁=100 GPa，E₂₂=10 GPa，ν₁₂=0.3, G₁₂=27.03 GPa，φ₁=0，材料2的弹性参数E₁₁=100 GPa，E₂₂/E₁₁=1、0.5、0.3、0.1，ν₁₂=0.3，G₁₂=27.03 GPa，φ₂=0，图 6给出了网格划分，W边划分4个线单元.表 4给出了无量纲化应力强度因子的计算结果.

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

3.5 重力坝算例

以Koyna重力坝(图 7)为例，计算了坝踵裂缝在任意分布水压力作用下的应力强度因子.假定裂缝长1.93 m，位于坝体与地基结合面.大坝承受满库水压力及自重，考虑平面应力状态.坝体混凝土与坝基均为各向同性材料，两者的泊松比均为0.25，混凝土密度为2 450 kg/m³，坝基自重不考虑，两者的弹性模量分别记为E₁和E₂，通过调整E₁/E₂研究结合材料差异对应力强度因子的影响.地基模拟范围为从坝体向上下游及向下延伸各2倍坝高.采用具有Gauss-Lobatto-Legendre形函数的11节点高阶线单元对坝和地基进行离散，整个模型共剖分7个子域，含108个线单元(图 8).以裂尖为原点，假定裂缝承受水压力σ₀(r)=λ(r/a)ⁿ，(n=0, 1, 2)，r为缝面上的节点与裂尖的距离，λ为裂缝口的静水压力.λ=0时对应缝内无水压，n=0, 1, 2分别代表水压力均匀分布、线性分布和按二次函数分布的情况.

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

4 结语

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