复合材料界面层的损伤破坏是导致其结构整体失效的重要因素之一. 在实际工程中，界面层的破坏可由多种荷载情况造成，如切向荷载、法向荷载以及复合型荷载. 复合材料界面层损伤破坏的力学机理复杂，破坏模式多样，因此建立能够准确描述界面破坏机理的理论模型是学术界和工程界共同关注的焦点^[1-2].

通过试验准确分析并研究界面层的非线性损伤断裂机理有较大困难，因此理论模型分析和数值模拟成为研究界面层力学性能的重要手段. Wong等^[3-5]采用内聚力模型(cohesive zone model， CZM)模拟界面层的非线性损伤破坏现象. 在假定损伤仅在界面内产生及扩展的前提下，内聚力模型将界面中的牵引力与界面上下表面的相对位移进行关联，以模拟界面层的黏结滑移关系^[6]. 有关CZM的研究通常假定牵引力与相对位移间的损伤关系是弹性的^[7]，即牵引力在初始状态下呈弹性增长，直至达到起始损伤强度，损伤演化才开始，牵引力降低直至为零. 内聚力理论模型大多基于“牵引力准则”定义^[8-9]. 在一维单一方向荷载下(法向或切向)，CZM中牵引力与相对位移的本构关系相对简单直接，且过往研究中采用的单向本构较为类似^[10]. 在复合加载下，CZM本构根据不同的损伤表达方式可分为独立(decoupled)损伤模型和耦合(coupled)损伤模型^[11]. 独立损伤模型在各加载方向上的损伤标量相互独立^[12]，但此模型存在2个明显的缺陷: 1)某一方向上的损伤标量仅与该方向的加载历史有关，与其他方向上的加载历史无关2)界面的整体破坏通过特定的破坏准则控制，在数值模型实现过程中可能存在不连续性问题. 对于耦合损伤模型，一种常用的方法是对法向和切向的相对位移进行耦合，形成有效相对位移，进而计算单一损伤标量^[7]；另一种方法是耦合各向的能量，推导定义单一损伤标量. 例如Needleman等^[8-9]通过对耦合的势能表达式进行关于各向位移差的偏导，得到切向和法向的牵引力，进而得到单一损伤标量. 此方法仅可用于相对位移单调增长且混合模态比恒定的状况下损伤的定义，在卸载及循环加载时并不适用. 为了解决这一问题，Nairn等^[13]提出从“强度理论”角度重新定义破坏模型. 在复合加载下，基于强度理论的内聚力模型(strength theory cohesive zone model， ST-CZM)不需要定义特定的破坏准则，其混合模态破坏面由各方向上相互独立的强度模型决定，混合模态比在加载和卸载过程中随加载历史改变. ST-CZM的提出很好地解决了上述各种CZM模型中存在的问题，且ST-CZM在准脆性材料断裂模拟复合材料界面层损伤模拟中具有一定的应用前景.

本研究拓展二维强度理论，得到基于强度理论的三维内聚力模型. 在此基础上，使用Abaqus平台的用户单元子程序(user-defined element， UEL)^[14]进行有限元实现. 建立3个验证模型，模拟多种破坏模式，验证模型的有效性和准确性，将STM-CZM应用于FRP-混凝土试件的单剪试验和四点弯曲梁试验数值模型中，模拟现实试验的复杂情况.

在一维单一方向荷载下(法向或切向)，CZM中的内聚力法则 $ F(u) $是牵引力关于黏结面相对位移 $ u $的函数. 在强度模型中，采用描述牵引力强度变化的函数 $ S(\delta ) $替换 $ F(u) $. 其中 $ \delta $为描述黏聚力损伤状态的变量， $\delta = \max \;(u)$. 双线性内聚力函数样例如图1所示. 图中，纵坐标表示黏结单元应力，横坐标表示位移. 当采用强度模型时，在相对位移 $ u $达到损伤起始位移 $ u_{\rm{e}} $之前，损伤尚未开始，其强度 $ S(\delta ) $保持为峰值应力 $ {\sigma _{\rm{c}}} $不变. 随着损伤的演变， $ S(\delta ) $的表达式为

(1) $ S(\delta ) = {\sigma _{\rm{c}}}\frac{{u_{\rm{c}} - \delta }}{{u_{\rm{c}} - u_{\rm{e}}}} . $

牵引力 $ T $关于相对位移的表达式为

(2) $ T(u) = \frac{{S(\delta )}}{\delta }u . $

定义2个损伤相关变量：

(3) $ D = \frac{{{u_{{\rm{dam}}}}}}{u} , $

(4) $ \zeta = \delta - {u_{{\rm{ela}}}} = \delta - \frac{{S(\delta )}}{k} = D\delta . $

式中： $ D $为损伤参数， $ {u_{{\rm{dam}}}} $为损伤相关的位移， $ {u_{{\rm{ela}}}} $为黏结区弹性位移， $ k $为初始刚度.

一维单一方向荷载下 $ F(u) $与 $ S(\delta ) $仅存在语义表达上的差异，但在混合模态推导过程中，二者差异显著，且区分F(u)与 $ S(\delta ) $的定义非常重要^[13]. 在复合加载状态下，此推导过程须分别定义代表法向和切向的独立强度^[15]. 牵引力矩阵 $ {\boldsymbol{T}} $和应变能矩阵 $ {\boldsymbol{U}} $关于相对位移矩阵、损伤矩阵的表达式为

(5) $ {\boldsymbol{T}} = {\boldsymbol{K}}({\boldsymbol{I}} - {\boldsymbol{D}}){\boldsymbol{u}} , $

(6) $ {\boldsymbol{U}}({\boldsymbol{u}},{\boldsymbol{D}}) = 0.5{\boldsymbol{K}}({\boldsymbol{I}} - {\boldsymbol{D}}){\boldsymbol{u}} {\boldsymbol{u}}^{\rm{T}} ， $

(7) $ {\boldsymbol{u}} = {\left[ {{u_n},{u_{t1}},{u_{t2}}} \right]^{\rm{T}}} , $

(8) $ {\boldsymbol{K}}={\rm{diag}}\;({{K}}_{{n}},{{K}}_{t}{}_{1},{{K}}_{t}{}_{2}) , $

(9) $ {\boldsymbol{D}} = {\rm{diag}}\;({{{D}}_{{n}}},{{{D}}_t}_1,{{{D}}_t}_2) . $

式中： $ {\boldsymbol{u}} $为相对位移矩阵，包含相对位移法向分量 $ {u_n} $，位移切向分量 $ {u_{t1}} $、 $ {u_{t2}} $； $ {\boldsymbol{K}} $为初始刚度矩阵，包含法向刚度 $ {{{K}}_{{n}}} $，切向刚度 $ {{{K}}_{{t}}}_{\text{1}} $、 $ {{{K}}_{{t}}}_{\text{2}} $； $ {\boldsymbol{D}} $为损伤矩阵，包含法向损伤标量 $ {{{D}}_{{n}}} $，切向损伤标量 $ {{{D}}_{{t}}}_{{1}} $、 $ {{{D}}_{{t}}}_{{2}} $. 损伤矩阵适用于独立损伤模型. 在耦合损伤模型中，采用单一损伤标量D定义模型的整体损伤.

二维相关矩阵由三维矩阵减少1个平面外切向分量得到. 假设二维破坏面如图2所示，三维破坏面基于二维拓展，如图3所示，表达式分别为

(10) $ {\left( {\frac{{{{\left\langle {{T_n}} \right\rangle }_+}}}{{{S_n}\left( {{\delta _n}} \right)}}} \right)^2}+{\left( {\frac{{{T_t}}}{{{S_t}\left( {{\delta _t}} \right)}}} \right)^2} = 1 , $

(11) $ {\left( {\frac{{{{\left\langle {{T_n}} \right\rangle }_+}}}{{{S_n}\left( {{\delta _n}} \right)}}} \right)^2}+{\left( {\frac{{{T_{t1}}}}{{{S_{t1}}\left( {{\delta _{t1}}} \right)}}} \right)^2}+{\left( {\frac{{{T_{t2}}}}{{{S_{t2}}\left( {{\delta _{t2}}} \right)}}} \right)^2} = 1 . $

式中： $ {T_n} $、 $ {T_t} $、 $ {T_{t1}} $、 $ {T_{t2}} $分别为各个方向的牵引力；在压缩条件下， $ {\left\langle {{T_n}} \right\rangle _+} $=0；在拉伸条件下， $ {\left\langle {{T_n}} \right\rangle _+} = {T_n} $.

混合模态须定义模型的混合模态比 $ \theta $，在本研究中，混合模态的计算仅与各相的相对位移比相关. 在二维状态下，

(12) $ \tan \theta = {{{u_n}}}/{{{{{u}}_t}}} \Rightarrow {{{u}}_n} = \left\| {\boldsymbol{u}} \right\|\sin \theta ，{u_t} = \left\| {\boldsymbol{u}} \right\|\cos \theta . $

在三维状态下，由标量 $ \varphi $、 $ \theta $表示，

(13) $ \begin{gathered} \tan \varphi = {{{u_n}}}/{{({u_{t1}^2+u_{t2}^2})^{1/2} }}，\tan \theta = {{{u_{t1}}}}/{{{u_{t2}}}} \Rightarrow \\ {u_n} = \left\| {\boldsymbol{u}} \right\|\sin \varphi ，{u_{t1}} = \left\| {\boldsymbol{u}} \right\|\cos \varphi \sin \theta ，{u_{t2}} = \left\| {\boldsymbol{u}} \right\|\cos \varphi \cos \theta . \end{gathered} $

引入损伤起始状态标量 ${\delta _{\rm{d}}}$，当 $\left\| {\boldsymbol{u }}\right\| \to {\delta _{\rm{d}}}$时，表示牵引力达到破坏面，

(14) $ {\left( {\frac{{{T_n}}}{{{S_n}\left( {{\delta _n}} \right)}}} \right)^2}+{\left( {\frac{{{T_t}}}{{{S_t}\left( {{\delta _t}} \right)}}} \right)^2} = \delta _{\rm{d}}^2\left( {\frac{{{{\sin }^2}\theta }}{{\delta _n^2}}+\frac{{{{\cos }^2}\theta }}{{\delta _t^2}}} \right) = 1 . $

当牵引力达到三维破坏面时，表示为

(15) $ \begin{split} {\left( {\frac{{{T_n}}}{{{S_n}\left( {{\delta _n}} \right)}}} \right)^2}+{\left( {\frac{{{T_{t1}}}}{{{S_{t1}}\left( {{\delta _{t1}}} \right)}}} \right)^2}+{\left( {\frac{{{T_{t2}}}}{{{S_{t2}}\left( {{\delta _{t2}}} \right)}}} \right)^2} = \\ \delta _{\rm{d}}^2\left( {\frac{{{{\sin }^2}\varphi }}{{\delta _n^2}}+\frac{{{{\cos }^2}\varphi {{\sin }^2}\theta }}{{\delta _{t1}^2}}+\frac{{{{\cos }^2}\varphi {{\cos }^2}\theta }}{{\delta _{t2}^2}}} \right) = 1, \end{split} $

(16) $ {\delta _{\rm{d}}} = \left(\dfrac{{{{\sin }^2}\varphi }}{{\delta _n^2}}+\dfrac{{{{\cos }^2}\varphi {{\sin }^2}\theta }}{{\delta _{t1}^2}}+\dfrac{{{{\cos }^2}\varphi {{\cos }^2}\theta }}{{\delta _{t2}^2}}\right)^{-1/2} . $

基于上述推导，得到在二维复合加载下强度S关于 $ {\delta _n} $、 $ {\delta _t} $和 $ \theta $的表达式为

(17) $ S\left( {{\delta _n},{\delta _t},\theta } \right) = \left({\dfrac{{\dfrac{{{S_n}{{\left( {{\delta _n}} \right)}^2}}}{{\delta _n^2}}{{\sin }^2}\theta +\dfrac{{{S_t}{{\left( {{\delta _t}} \right)}^2}}}{{\delta _t^2}}{{\cos }^2}\theta }}{{\dfrac{{{{\sin }^2}\theta }}{{\delta _n^2}}+\dfrac{{{{\cos }^2}\theta }}{{\delta _t^2}}}}}\right)^{1/2} . $

三维复合加载下强度S关于 $ {\delta _n} $、 $ {\delta _{t1}} $、 $ {\delta _{t2}} $、 $ \theta $和 $ \varphi $的表达式为

(18) $ \begin{split} {U_1}{\text{ = }}& \frac{{{S_n}{{\left( {{\delta _n}} \right)}^2}}}{{\delta _n^2}}{\sin ^2}\varphi +\frac{{{S_{t1}}{{\left( {{\delta _{t1}}} \right)}^2}}}{{\delta _{t1}^2}}{\cos ^2}\varphi {\sin ^2}\theta + \\ & \frac{{{S_{t2}}{{\left( {{\delta _{t2}}} \right)}^2}}}{{\delta _{t2}^2}}{\cos ^2}\varphi {\cos ^2}\theta , \end{split} $

(19) $ {U_2}{\text{ = }}\frac{{{{\sin }^2}\varphi }}{{\delta _n^2}}+\frac{{{{\cos }^2}\varphi {{\sin }^2}\theta }}{{\delta _{t1}^2}}+\frac{{{{\cos }^2}\varphi {{\cos }^2}\theta }}{{\delta _{t2}^2}} , $

(20) $ S\left( {{\delta _n},{\delta _{t1}},{\delta _{t2}},\varphi ,\theta } \right) = ({{{{U_1}}}/{{{U_2}}}})^{1/2} . $

采用牵引力回归的方法完成对损伤矩阵的更新，即让超出破坏面的试牵引力回归更新的破坏面，以求解当前状态下的牵引力值. 该方法二维状态下的求解式为

(21) $ \begin{split} {\left( {\frac{{{T_n}\left( {{u_n}+{\rm{d}}{u_n}} \right)}}{{{S_n}\left( {{\delta _n}+{\rm{d}}{\delta _n}} \right)}}} \right)^2}+{\left( {\frac{{{T_t}\left( {{u_t}+{\rm{d}}{u_t}} \right)}}{{{S_t}\left( {{\delta _t}+{\rm{d}}{\delta _t}} \right)}}} \right)^2} = \\\quad {\left( {\frac{{{u_n}+{\rm{d}}{u_n}}}{{{\delta _n}+{\rm{d}}{\delta _n}}}} \right)^2}+{\left( {\frac{{{u_t}+{\rm{d}}{u_t}}}{{{\delta _t}+{\rm{d}}{\delta _t}}}} \right)^2} = 1. \end{split} $

三维状态下的求解式为

(22) $ \begin{split} {\left( {\frac{{{T_n}\left( {{u_n}+{\rm{d}}{u_n}} \right)}}{{{S_n}\left( {{\delta _n}+{\rm{d}}{\delta _n}} \right)}}} \right)^2}+{\left( {\frac{{{T_{t1}}\left( {{u_{t1}}+{\rm{d}}{u_{t1}}} \right)}}{{{S_{t1}}\left( {{\delta _{t1}}+{\rm{d}}{\delta _{t1}}} \right)}}} \right)^2}+ \\\quad {\left( {\frac{{{T_{t2}}\left( {{u_{t2}}+{\rm{d}}{u_{t2}}} \right)}}{{{S_{t2}}\left( {{\delta _{t2}}+{\rm{d}}{\delta _{t2}}} \right)}}} \right)^2} = {\left( {\frac{{{u_n}+{\rm{d}}{u_n}}}{{{\delta _n}+{\rm{d}}{\delta _n}}}} \right)^2}+ \\ {\left( {\frac{{{u_{t1}}+{\rm{d}}{u_{t1}}}}{{{\delta _{t1}}+{\rm{d}}{\delta _{t1}}}}} \right)^2}+{\left( {\frac{{{u_{t2}}+{\rm{d}}{u_{t2}}}}{{{\delta _{{t2}}}+{\rm{d}}{\delta _{{t2}}}}}} \right)^2} = 1. \end{split} $

将式(21)、(22)分别以泰勒级数展开，在增量中仅保留一阶项^[13]，分别得到

(23) $ \;\; \frac{{u_n^2}}{{\delta _n^3}}{\rm{d}}{\delta _n}+\frac{{u_t^2}}{{\delta _t^3}}{\rm{d}}{\delta _t} = \frac{{{u_n}}}{{\delta _n^2}}{\rm{d}}{u_n}+\frac{{{u_t}}}{{\delta _t^2}}{\rm{d}}{u_t} , $

(24) $ \begin{split} & \frac{{u_n^2}}{{\delta _n^3}}{\rm{d}}{\delta _n}+\frac{{u_{t1}^2}}{{\delta _{t1}^3}}{\rm{d}}{\delta _{{t1}}}+\frac{{u_{t2}^2}}{{\delta _{t2}^3}}{\rm{d}}{\delta _{{t2}}} = \\ & \qquad \frac{{{u_n}}}{{\delta _n^2}}{\rm{d}}{u_n}+\frac{{{u_{{t1}}}}}{{\delta _{{t1}}^2}}{\rm{d}}{u_{{t1}}}+\frac{{{u_{{t2}}}}}{{\delta _{{t2}}^2}}{\rm{d}}{u_{{t2}}} . \end{split} $

式(23)、(24)的结果并不能完成对多个损伤标量的更新，还需要另外的求解条件. 因此在耦合损伤模型的假定下，假设各向的损伤标量相同且可以用统一的损伤标量D定义，得到 $ {D_n} = {D_{t1}} = {D_{t2}} = D $.定义函数 $ R(\delta ) $以关联损伤标量的增量 $ {\rm{d}}D $和损伤状态变量的增量 $ {\rm{d}}\delta $.

(25) $ R(\delta ) = {{{\rm{d}}D}}/{{{\rm{d}}\delta }} = {{\varphi (\delta )}}/{{k{\delta ^2}}} , $

(26) $ \varphi (\delta ) = S(\delta ) - \delta S'(\delta ) . $

3个方向的 $ R(\delta ) $函数分别为 ${R_n}({\delta _n}) = {\rm{d}}{D_n}/{\rm{d}}{\delta _n}$、 ${R_{{t1}}}({\delta _{{t1}}}) = {{\rm{d}}}{D_{{t1}}}/{{\rm{d}}}{\delta _{{t1}}}、{R_{{t2}}}({\delta _{{t2}}}) = {{\rm{d}}}{D_{{t2}}}/{{\rm{d}}}{\delta _{{t2}}}$. 根据以上定义得到损伤标量最终的增量表达式为

(27) $ {\rm{d}}D = \dfrac{{\dfrac{{{u_n}}}{{\delta _n^2}}{\rm{d}}{u_n}+\dfrac{{{u_t}}}{{\delta _t^2}}{\rm{d}}{u_t}}}{{\dfrac{1}{{{R_n}\left( {{\delta _n}} \right)}}\dfrac{{u_n^2}}{{\delta _n^3}}+\dfrac{1}{{{R_t}\left( {{\delta _t}} \right)}}\dfrac{{u_t^2}}{{\delta _t^3}}}} . $

三维损伤标量的增量表达式为

(28) $ {{\rm{d}}}D = \dfrac{{\dfrac{{{u_n}}}{{\delta _n^2}}{\rm{d}}{u_n}+\dfrac{{{u_{{t1}}}}}{{\delta _{{t1}}^2}}{\rm{d}}{u_{{t1}}}+\dfrac{{{u_{{t2}}}}}{{\delta _{{t2}}^2}}{\rm{d}}{u_{{t2}}}}}{{\dfrac{1}{{{R_n}\left( {{\delta _n}} \right)}}\dfrac{{u_n^2}}{{\delta _n^3}}+\dfrac{1}{{{R_{{t1}}}\left( {{\delta _{{t1}}}} \right)}}\dfrac{{u_{{t1}}^2}}{{\delta _{{t1}}^3}}+\dfrac{1}{{{R_{{t2}}}\left( {{\delta _{{t2}}}} \right)}}\dfrac{{u_{{t2}}^2}}{{\delta _{{t2}}^3}}}} . $

通过式(27)、(28)分别更新二维混合模态和三维混合模态下的损伤标量.

为了充分利用Abaqus中前后处理系统和非线性求解器，采用Abaqus用户单元子程序UEL实现基于强度理论的内聚力模型. UEL实现过程中最重要的是定义ST-CZM单元的雅可比刚度矩阵及右手系向量数组，其中决定雅可比刚度矩阵的关键是得到联系内聚力增量和界面相对位移增量的本构关系，即内聚力切向刚度矩阵. 在损伤演化过程中，切向刚度矩阵与损伤标量D同步更新.

(29) $ {\rm{d}}{\boldsymbol{T}} = (1 - D){\boldsymbol{K}}{\rm{d}}{\boldsymbol{u}} - {\boldsymbol{K}}{\boldsymbol{u}}{\rm{d}}D , $

将 $ R(\delta ) $与 $ {\rm{d}}D $代入式(29)，得到二维切向刚度矩阵、三维切向刚度矩阵分别为

(30) $ \begin{aligned} &\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{T_n}} \\ {{\rm{d}}{T_t}} \end{array}} \right] = \left({{\dfrac{{{k_n}u_n^2}}{{{\varphi _n}{\delta _n}}}+\dfrac{{{k_t}u_t^2}}{{{\varphi _t}{\delta _t}}}}}\right)^{-1}\times\\& \left[\begin{array}{cc} \dfrac{S_n' k_n u_n^2}{\varphi_n \delta_n}+\dfrac{S_n k_t u_t^2}{\varphi_t \delta_t \delta_n} & -\dfrac{k_n u_n u_t}{\delta_t^2} \\ -\dfrac{k_t u_t u_n}{\delta_n^2} & \dfrac{S_t' k_t u_t^2}{\varphi_t \delta_t}+\dfrac{S_t k_n u_n^2}{\varphi_n \delta_n \delta_t} \end{array}\right]\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{u_n}} \\ {{\rm{d}}{u_t}} \end{array}} \right] . \end{aligned} $

(31) $ A = \left({{\dfrac{{{k_n}u_n^2}}{{{\varphi _n}{\delta _n}}}+\dfrac{{{k_{{t1}}}u_{{t1}}^2}}{{{\varphi _{{t1}}}{\delta _{{t1}}}}}+\dfrac{{{k_{{t2}}}u_{{t2}}^2}}{{{\varphi _{{t2}}}{\delta _{{t2}}}}}}} \right)^{-1}. $

(32) $ \left[ {\begin{array}{*{20}{c}} {{\rm{d}}{T_n}} \\ {\rm{d}}{T_{{t1}}} \\ {\rm{d}}T_{t2} \end{array}} \right] = A\left[ {\begin{array}{*{20}{c}} {\dfrac{{S_n^\prime {k_n}u_n^2}}{{{\varphi _n}{\delta _n}}}+\dfrac{{{S_n}{k_{{t1}}}u_{t1}^2}}{{{\varphi _{{t1}}}{\delta _{{t1}}}{\delta _n}}}+\dfrac{{{S_n}{k_{{t2}}}u_{{t2}}^2}}{{{\varphi _{{t2}}}{\delta _{{t2}}}{\delta _n}}}}&{ - \dfrac{{{k_n}{u_n}{u_{{t1}}}}}{{\delta _{{t1}}^2}}}&{ - \dfrac{{{k_n}{u_n}{u_{{t2}}}}}{{\delta _{{t2}}^2}}} \\ { - \dfrac{{{k_{{t1}}}{u_{{t1}}}{u_n}}}{{\delta _n^2}}}&{\dfrac{{S_{{t1}}^\prime {k_{{t1}}}u_{{t1}}^2}}{{{\varphi _{{t1}}}{\delta _{{t1}}}}}+\dfrac{{{S_{{t1}}}{k_n}u_n^2}}{{{\varphi _n}{\delta _n}{\delta _{{t1}}}}}+\dfrac{{{S_{{t1}}}{k_{{t2}}}u_{{t2}}^2}}{{{\varphi _{{t2}}}{\delta _{{t1}}}{\delta _{{t2}}}}}}&{ - \dfrac{{{k_{{t1}}}{u_{{t1}}}{u_{{t2}}}}}{{\delta _{{t2}}^2}}} \\ { - \dfrac{{{k_{{t2}}}{u_{{t2}}}{u_n}}}{{\delta _n^2}}}&{ - \dfrac{{{k_{{t2}}}{u_{{t2}}}{u_{{t1}}}}}{{\delta _{{t1}}^2}}}&{\dfrac{{S_{{t2}}^\prime {k_{{t2}}}u_{{t2}}^2}}{{{\varphi _{{t2}}}{\delta _{{t2}}}}}+\dfrac{{{S_{{t2}}}{k_n}u_n^2}}{{{\varphi _n}{\delta _n}{\delta _{{t2}}}}}+\dfrac{{{S_{{t2}}}{k_{{t1}}}u_{{t1}}^2}}{{{\varphi _{{t1}}}{\delta _{{t1}}}{\delta _{{t2}}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{u_n}} \\ {\rm{d}}{u_{{t1}}}\\ {\rm{d}}{u_{{t2}}} \end{array}} \right] . $

根据式(5)和切向刚度矩阵得到内聚力矩阵.

经典验证模型分别为Garg等^[16]的双悬臂梁模型(double cantilever beam， DCB)、固定比率混合模型(fixed ratio mixed mode， FRMM)以及Park等^[17]的混合模态弯曲模型(mixed mode bending， MMB). 将ST-CZM与传统双线性内聚力模型及试验数据进行对比. 内聚力单元的界面参数如表1^[18]所示. 表中， ${\sigma _{\rm{I}}}$为法向峰值应力， ${u_{\rm{eI}}}$为法向损伤起始位移， ${K_{\rm{I}}}$为法向初始刚度， ${G_{\rm{I}}}$为法向断裂能， ${\sigma _{\rm{II}}}$为切向峰值应力， ${u_{\rm{eII}}}$为切向损伤起始位移， ${K_{\rm{II}}}$为切向初始刚度， ${G_{\rm{II}}}$为切向断裂能. 验证算例模型的各项材料参数如表2所示^[16]. 表中， ${E_{{\text{11}}}}$、 ${E_{{\text{22}}}}$、 ${E_{{\text{33}}}}$分别为3个方向的弹性模量， ${G_{{\text{12}}}}$、 ${G_{{\text{13}}}}$、 ${G_{{\text{23}}}}$分别为3个方向的切变模量， ${v_{{\text{12}}}}$、 ${v_{{\text{13}}}}$、 ${v_{{\text{23}}}}$分别为3个方向的泊松比.

DCB数值模型试件尺寸如图4所示，数值模型边界条件按试验实际取用^[16]. 如图5所示为DCB采用ST-CZM计算的数值结果与Abaqus自带双线性内聚力单元结果、试验数据的对比图. 图中， ${P}$为荷载， ${s}$为张开位移，DCB模型计算结果使用自带单元模拟，所得峰值相较于实验结果误差为4.44%，ST-CZM相较于实验结果误差为0.41%. 明显ST-CZM对峰值的预测更接近试验的实际情况，其对损伤阶段的预测符合试验的破坏模式.

FRMM的几何尺寸与DCB的几何尺寸相同，边界条件不同. FRMM上端受到荷载作用，下端不受约束. FRMM的边界条件按试验实际取用，如图6所示^[16]. 如图7所示为FRMM采用ST-CZM计算的数值结果与自带双线性内聚力单元以及LS-DYNA^[19]模拟结果的对比图. 可以看出，3种结果破坏模式相近，三者峰值强度相差不超过2 N，但LS-DYNA计算得到的结果曲线与其余2条相差较大. 原因是LS-DYNA求解采用显示动力学算法模拟静态加载时得到的结果，因此受加载速率影响较大^[19].

MMB几何尺寸及边界条件如图8所示^[17]. 图中， $\varDelta $为开口位移；L为试件的长度；L=102 mm；单梁高1.56 mm；d_s为试件厚度，d_s=25.4 mm. 如图9所示为采用ST-CZM计算的数值结果与基于势能的内聚单元模型(potential-based cohesive zone model， PPR)单元^[17]的结果、解析解的对比图，其中解析解基于线性梁理论以及线弹性断裂力学得到. 可以看出，MMB计算结果中使用PPR单元所得峰值相较于解析解误差为4.07%，自带单元相较于解析解误差为11.35%，ST-CZM相较于解析解误差为3.1%. ST-CZM对峰值的预测与Park等^[17]的结果相近，且破坏模式相同，损伤阶段结果与解析解符合较好. 当改变MMB的界面单元参数时(如修改 ${\sigma _{\rm{I}}}$=80 MPa)，使用Abaqus自带CZM单元无法收敛. 原因在于加载过程中MMB的混合模态比是变化的，自带单元很难灵活应对. 自带单元对损伤的定义式^[20]为

(33) $ d = \frac{{\delta _{\rm{m}}^{\rm{f}}\left( {\delta _{\rm{m}}^{\max } - \delta _{\rm{m}}^{\rm{o}}} \right)}}{{\delta _{\rm{m}}^{\max }\left( {\delta _{\rm{m}}^{\rm{f}} - \delta _{\rm{m}}^{\rm{o}}} \right)}} . $

式中： $ \delta _{\rm{m}}^{\rm{f}} $为混合模态相对位移^[20]， $ \delta _{\rm{m}}^{\max } $为混合模态破坏位移， $ \delta _{\rm{m}}^{\rm{o}} $为混合模态损伤起始位移. 由式可知，只有当 $ \delta _{\rm{m}}^{\max }/\delta _{\rm{m}}^{\rm{o}} $不变时，混合模态比的变化才不会影响损伤的定义. 否则当损伤产生时，刚度矩阵出现突变^[13]. ST-CZM没有这个缺陷，表明ST-CZM相较自带单元拥有更灵活的耦合方式.

以上3个经典算例的验证表明了ST-CZM的准确性，且ST-CZM模拟精度比传统内聚力模型有较大的提高. 为了进一步研究ST-CZM在实际工程中的应用，使用实用算例验证ST-CZM的应用效果.

纤维增强复合材料(fiber reinforced polymer， FRP)加固钢筋混凝土结构的技术应用范围广、研究历史长、破坏模式复杂，与黏结界面损伤力学模型的研究联系紧密^[21-22]. 选取几组典型的FRP-混凝土试验作为实用算例，验证ST-CZM的可靠性. 所建单剪试验数值模型及四点弯曲梁试验数值模型分别代表素混凝土和钢筋混凝土试件试验，参照Zhang等^[23-24]的试验进行结果对比. 采用的界面参数如表3所示.

剪切拉拔是常见的研究FRP-混凝土界面的试验方法^[25-26]. 选取单剪试验结果^[23]与数值模型对照，以确定ST-CZM在单一模态下的可靠性. 素混凝土数值模型长300 mm，宽100 mm，高50 mm. 标准测试得到的试件抗压强度为29.3 MPa. FRP板数值模型厚度为2 mm， FRP试件测试得到的拉伸强度为4.1 GPa，弹性模量为256 GPa，黏结长度为180 mm，黏结宽度为20 mm. 文献[27]~[29]中的试验也采用类似试件. 试件几何尺寸如图10所示，数值模型一切边界条件按照试验实际取用^[23].

如图11所示为应用ST-CZM的单剪试验数值结果与试验数据的对比图. 荷载取自FRP末端加载力，平均加载位移FRP与混凝土黏结末端的相对位移. 由图可以得到，ST-CZM模拟结果和试验结果接近，且破坏过程都分为损伤开始前和开始后2个阶段. 如图12所示为破坏过程中FRP-混凝土界面的4个具有代表性的损伤云图. 该云图只表示损伤的发展趋势，损伤具体数值受到试件实际情况的影响. 由图可以得到，模拟和试验的破坏模式相同且破坏过程相近，说明ST-CZM对剪切破坏的模拟有效且准确.

FRP-混凝土梁黏结界面的失效主要是跨中局部脱黏导致的^[30-32]. 采用的四点弯曲梁试验模型根据Fu等^[24]的试验确定，混凝土内嵌钢筋. 圆柱体抗压强度为47 MPa，四点弯曲梁混凝土试件的长4 000 mm，宽200 mm，高450 mm. 加强钢筋采用2根16 mm受拉钢筋和2根16 mm受压钢筋，采用8 mm箍筋，箍筋间距为100 mm. 底部粘贴的FRP板长3 800 mm，宽100 mm，厚0.999 mm，试件几何尺寸如图13所示. 图中，a为加载点与试件边缘的距离. 将a=1 250 mm的参考试件命名为L2D1250，a=1 750 mm的参考试件命名为L2D1750.

模拟L2D1250试验时，采用的自带单元是稳定且相对成熟的双线性内聚力单元. 采用本研究推导的ST-CZM实现UEL. L2D1250试验裂缝图、局部脱黏标记图以及模拟损伤云图分别如图14~16所示. 其中图15(b)、(c)为四点弯曲梁模型的底面损伤云图. 选取黏结单元首次出现集中局部损伤的计算步，根据厚度方向损伤平均值得到损伤云图. 观察图15可以得到，试验得到的局部脱黏位置与ST-CZM预测的黏结单元最早出现集中损伤的位置接近. 如图17所示为L2D1250采用ST-CZM、自带单元的模拟结果与试验结果的对比图. 图中， $w$为跨中挠度， $M$为跨中弯矩. 可以得到，试验曲线与ST-CZM的模拟结果吻合情况较好，说明了ST-CZM在实际工况模拟应用中的可靠性. FRP的加强作用并没有被自带单元准确估计，具体表现为脱黏发生得更早且界面的强度被低估.

过于靠近的加载端可能会影响试验以及模拟的结果^[32]，本研究建立加载跨度较大的L2D1750，直观表现具有多种破坏模式的复合模态下的界面破坏情况.

模型L2D1750采用的自带单元与L2D1250的相同. L2D1750试验的试验裂缝图、局部脱黏标记图以及模拟损伤云图分别如图18~20所示. 其中图19(b)、(c)分别为黏结单元损伤云图和基体损伤云图，得到云图的方法与L2D1250相同. 观察图19可以得到，试验得到的局部脱黏位置与ST-CZM预测的黏结单元最早出现集中损伤的位置接近. 如图21所示为L2D1750的试件分别采用自带单元、ST-CZM单元的模拟结果与试验结果的对比图. 分析可以得到1）试件破坏过程可以明显地划分为3阶段. 阶段1为斜率几乎不变的弹性阶段；阶段2开始损伤出现，曲线的斜率第1次减小，混凝土内部钢筋和FRP逐渐发挥加强作用；阶段3混凝土裂缝基本成形，斜率第2次减小. 观察模拟结果可以得出，自带单元对界面强度的估计值偏低，且计算结果更不稳定. 2）通过对比图18~20的结果可以看出，试验结果和数值结果的裂缝发展趋势接近. 分析原因，受到弯矩作用，在裂缝向上发展之前混凝土下半部分已经完全破坏，这是裂缝的相互作用导致的最终破坏结果^[31]. 3）数值模型中破坏集中的位置与试验得到的局部脱黏位置接近. 局部脱黏是导致 FRP 最终的脱落和失效的主要原因. FRP-混凝土黏结界面的破坏模式复杂，传统模型在这种模式下想要准确得到局部脱黏的位置有比较大的困难，而ST-CZM能够较准确地复现试验的破坏情况. 这也是本研究实现的ST-CZM的优越性之一.

ST-CZM从强度的角度重新定义了破坏模型，解决了部分模型的弊端、在复合荷载下，ST-CZM的破坏面由各方向上相互独立的强度模型决定，具有更灵活的耦合方式. 在实际应用中，与传统的黏结单元本构相比，ST-CZM能够更准确地模拟复合加载下FRP的局部脱黏位置. 本研究通过理论推导将基于强度理论的内聚力模型从二维拓展到三维，并在Abaqus平台利用UEL子程序开发了ST-CZM模型的用户自定义单元，实现了基于强度理论的内聚力模型在通用有限元平台的应用. 经过双悬臂梁模型、固定比率混合模型、混合模态弯曲模型的验证，证明了ST-CZM的准确性. 在此基础上，本研究对FRP加固混凝土结构的单剪试验和四点弯曲梁试验进行数值模拟，并与试验结果对比，证明了ST-CZM在模拟FRP与混凝土界面复杂破坏模式的有效性. 在未来的研究中，计划将ST-CZM拓展应用到不同类型材料界面的模拟破坏. 本研究为此模型的推广和应用提供了可靠的实现路径.