彭 刚， 洪海丰， 陈灯红， 杨紫辉， 赵艺园， 曹文昱

(1. 三峡大学 防灾减灾湖北省重点实验室，湖北 宜昌 443002；2. 三峡大学 土木与建筑学院，湖北 宜昌 443002；3. 三峡大学 水利与环境学院，湖北 宜昌 443002)

我国在西南地区建设了一批300 m级的大库容超高拱坝，如锦屏一级(305 m)、小湾(294.5 m)、白鹤滩(289 m)等重大水利水电工程[1]。西南地区地震活动频繁，且该地区拥有我国约80%的水能资源[2-3]，在该区修建超高拱坝其抗震安全问题日益突出。横缝工作性态是超高拱坝抗震安全评价的重要组成部分[4]，考虑横缝动接触对大坝的抗震安全[5]和抗震设计具有十分重要的意义。

在拱坝建设过程中通常用带键槽的横缝把坝体分割成柱状浇筑块，待坝体冷却到一定的温度后，进行接缝灌浆处理，使坝体形成整体，提高超载能力[6]。Niwa等[7]于1982年首次指出拱坝在强震作用下横缝可能会发生张开、闭合及滑动的现象，这种现象将对大坝的位移和应力分布产生影响。横缝的变形和受力过程看作是摩擦接触现象，使得坝体静动力响应呈现非线性特性。此后，众多学者根据横缝的几何、变形特性进行了更深入的研究。在实际工程分析中，拱坝横缝的模拟力学模型主要有联结单元模型[8]、裂缝涂抹式模型和动接触力模型[9]。陈厚群等[10]对Fenves提出的界面单元模型中进行了改进，对小湾拱坝进行了动力分析。Kuo[11]于1982年提出了裂缝涂抹模型，引进接缝无滑动假定，认为接触面有足够的摩擦以阻止界面滑动。牛智勇等[12-13]在B-可微方程组形式的横缝接触模型基础上，针对法向、切向约束条件进行了修正，提出了考虑横缝初始抗拉强度、键槽切向抗剪力与张开度相关的修正接触模型。Alembagheri等[14]通过非线性增量动力分析，研究了有横缝拱坝的抗震性能和安全性。田硕等[15-16]采用扩大基础等措施会降低坝体中部的横缝开度；采用考虑键槽咬合作用、横缝设置等的非线性动力分析模型，并对白鹤滩拱坝进行地震易损性分析。Liang等[17-18]采用拉格朗日乘子法建立了动态接触模型以考虑横缝接触的非线性。李海枫等[19]采用带开-闭迭代功能的接触单元模拟横缝的力学特性。上述学者基本是采用单一的横缝力学模型。联结单元模型可以模拟横缝的非线性关系，但不能保证接触体之间不发生相互嵌入，李静等[20]根据经验和判断选取单元的法向和切向刚度系数。动接触模型把界面影响看作是接触问题，通过在接触面上施加约束条件，来确定满足接触约束条件和平衡条件的接触力，常忽略键槽的咬合作用。因此，本文将探究动接触-联结混合模型模拟考虑键槽咬合作用的横缝接触非线性。

目前，工程实际常用的键槽形状有矩形键槽[21]、梯形键槽、弧形键槽、球形键槽。考虑梯形键槽咬合作用的横缝模型有精细键槽模型、键槽深度和坡度不变的简化键槽模型、键槽坡度不变的简化键槽模型[22-23]、平缝模型(附加切向弹簧刚度)[24]四类。本文在文献[25]的基础上，将从动接触-联结混合模型的力学模型出发，以精细键槽模型在位移荷载作用下参考点的相对位移为基准，对比分析了简化模型一、简化模型二和平缝模型的键槽咬合作用，并对平缝模型附加切向弹簧刚度系数的取值及其稳定性展开研究，以锦屏一级拱坝为研究对象，建立三维有限元分析模型，考虑键槽咬合作用的接触横缝模型研究其对拱坝的整体性作用和拱坝横缝的切向位移及横缝的法向开度影响。

1 动接触-联结混合模型

横缝接触状态的改变在力学上属于边界非线性。如图1所示的动接触边界，由主面和从面构成，从面上的结点S′在主面上有且仅有一个确定的结点M′(锚点)与之对应，锚点与该结点的连线方向即为接触的法线方向；从面结点与锚点之间沿法线和切线方向的相对距离，分别表征动接触的横缝开度与切向位移。通过从面结点与主面锚点之间的相对位移可以唯一确定接触的状态。

对于任意从面结点S′在主面区间(M-1,M)的对应点M′，参数βS可以通过式(1)来确定。

(1)

从面结点S′到主面区间(M-1,M)的法向和切向距离可由式(2)确定。当h>0时，主从面接触；当h<0时，主从面分离；当h=0时，主从面接触且相互之间没有力的作用。

(2)

由变分法推导出横缝开度h的变分为

δh=-n(δuM′-δxM-1-ξδv)

(3)

为改善迭代计算的收敛性，并且保证收敛条件满足精度要求，引入指数型应力-位移模型来描述接触面的法向接触本构，其数学表达式为

(4)

(5)

式中：c为初始间隙；h为接触面之间的相对位移(以嵌入为正)；p为接触点对上的接触压力；p0为特征接触力。式(5)反映了由法向开度的微小摄动引起的接触力变化，即法向的刚度贡献。若选取c为一合理的微量，并保证p0值足够大，可以将接触面间的嵌入深度控制为可接受的微量(近似满足接触面间法向不可嵌入)，能够满足动力计算缝两侧坝体在张开状态下无相互作用，在缝面接触时保持接触面之间的法向特性光滑过渡，从而使得采用隐式迭代算法时的收敛性相对较好。

应用变分原理，得到接触边界上法向势能的变分方程

δ∏=pδh+hδp=[p+(∂p/∂h)h]δh

(6)

结合式(4)、式(5)，代入式(6)中，并运用链式求导法则，可推导出

(7)

将式(6)中的节点位移变分向量的系数集成到总体刚度矩阵中，形成包含接触边界的总体刚度矩阵，并由此推导出基于牛顿迭代的有限元隐式求解算法。在平缝模型中，通过附加切向弹簧约束相邻坝段间的顺河向位移，来考虑横缝键槽的咬合作用。采用联结单元模型中的点对弹簧单元(如图2所示)，沿缝面布置附加切向弹簧刚度为KS的线性弹簧单元，通过调整KS可达到控制切向位移的目的。本文第3章将对KS的取值进行了讨论。

2 梯形键槽咬合作用

2.1 梯形键槽的传力特点

梯形键槽传递剪力主要是通过啮合来进行的，单向梯形键槽(如图3(a)所示)传递剪力时，键槽两边斜面相互接触。坝体横缝灌浆前，接触面法向存在初始开度d0(如图3(b)所示)，梯形键槽的倾斜角为β，则键槽平行于缝面的初始间隙d′可由式(8)求得。

m=d′=d0tanβ

(8)

梯形键槽的三种接触状态：张开、接触、闭合。横缝接触面切向发生相对错动时：当|Δu|-(Δwn+d0)/m<0时，键槽面之间不接触，缝面自由(如图4(a)所示)；当|Δu|-(Δwn+d0)/m≥0时，键槽面之间会发生接触(如图4(b)所示)；当Δu=0且Δwn+d0=0时，缝面闭合不产生切向咬合力(如图4(c)所示)。

如图3(b)所示，横缝面上的荷载由法向荷载和切向荷载。其中，法向荷载引起的应力利用第1章的指数型应力-位移模型来描述；将切向荷载引起的应力分为两部分：切向摩擦引起的应力τf和附加切向约束(即键槽咬合)引起的切向应力τtt；将τf和τtt在横缝面上均匀化后用切向应力表达，即

τt=τf+τtt
τf=μσn

(9)

式中：τf为横缝面切向摩擦力引起的应力，可用库仑摩擦来表征；σn为横缝面法向接触力引起的应力；μ为摩擦因数。

由式(9)可知，当横缝张开时，σn=0,τf=0，因此，横缝面切向应力有键槽承担，体现单键槽的的传剪作用的附加切向应力τtt为

(10)

式中：Ktt反应键槽咬合作用的附件切向刚度系数；Δu为横缝面切向相对位移；Δwn为横缝面法向相对位移；d0为横缝初始开度；m为键槽坡度。当Δu>1，sign(Δu)=1，当Δu<1，sign(Δu)=-1。

2.2 简化键槽模型与精细键槽模型咬合作用对比

为探究简化键槽模型与精细键槽模型咬合作用，计算了键槽深度和坡度不变的简化模型(以下简称：简化模型一，如图5(a)所示)、键槽坡度不变的简化模型(以下简称：简化模型二，如图5(b)所示)、反映键槽实际尺寸的精细模型(以下简称：精细模型，如图5(c)所示)的参考点相对位移。

(11)

(12)

键槽为单向梯形键槽，键槽两边混凝土块体的弹性模量均为25.5 GPa，泊松比为0.189。模型的右侧面水平和竖直方向为固定边界；左侧为自由边界，对左侧面施加如图6的竖向位移时程。

计算模型网格尺寸为0.122 m×0.200 m，间距d0=0.005 m，n=2时参考点BB’、CC’、DD’的相对位移如图7(a)所示。无论采用简化模型一还是简化模型二，不考虑附加切向刚度影响，就无法使简化模型达到与精细模型同等的咬合作用。当n=2时，简化模型一、简化模型二、精细模型竖向位移的负向峰值分别为-0.028 8 m，-0.025 0 m，-0.018 0 m，简化模型一、简化模型二与精细模型竖向位移的相对误差分别为60.0%，38.9%。说明简化模型二的咬合作用优于简化模型一，但两种简化模型的相对误差较大。参考点BB’，CC’，DD’的法向开度如图7(b)所示，简化模型一、简化模型二、精细模型法向开度的峰值分别为0.016 1 m，0.016 5 m，0.016 8 m，简化模型一、简化模型二与精细模型法向开度的相对误差分别为4.2%，1.8%，对比三个模型法向开度变化很小。

3 平缝模型附加切向弹簧刚度取值研究

通过含有两个键槽的精细模型与同尺寸带附加弹簧平缝模型的数值试验来确定附加切向弹簧刚度KS的取值。带键槽横缝接触面切向相互作用力由两部分组成：库伦摩擦力[27]和附加切向力(键槽咬合力)。当横缝张开时，切向摩擦力为0，只有附加切向力发生作用。因此，横缝张开情况下，根据两侧的变形来确定平缝模型的附加切向弹簧刚度KS的值。

3.1 附加切向弹簧刚度的取值

为确定合理的刚度系数，建立了精细键槽模型(示意图如图5(c))和平缝模型(示意图如图8(a)，计算模型网格如图8(b))。计算模型分为左右两部分，中间设横缝，横缝接触面的摩擦因数取为0.8，横缝间隙d0=0.005 m，混凝土弹性模量为25.5 GPa，泊松比为0.189。施加的边界条件和位移加载方式同2.2节算例，竖向位移加载时程如图6所示。精细键槽模型竖直面不发生接触，左右块体间的作用力通过键槽倾斜面来传递。以精细键槽模型中横缝两侧的竖向变形为依据，为平缝模型选取合适的附加切向弹簧刚度，使两模型的竖向变形趋于一致。

当n=2时，平缝模型和精细模型参考点的相对竖向位移曲线，如图9所示，沿缝面共布置17个点对弹簧单元，其中附加切向弹簧刚度KS=0.15 GPa/m，单位与文献[28-29]中横缝属性的单位一致。此时参考点A，A′竖向相对位移正负向峰值为0.026 5 m,-0.019 4 m，与精细模型参考点竖向相对位移正负向峰值+0.024 2 m，-0.018 0 m的相对误差分别为8.6%，7.8%。

相较于简化模型一、简化模型二60.0%，38.9%的相对误差，附加切向弹簧平缝模型的相对误差明显减小，即在10%以内。对平缝模型附加切向弹簧刚度的试算中，参考点竖向相对位移正负峰值的相对误差很难控制在5%以内。因此，本文在这里将正负峰值相对误差调整至10%。综上所述，考虑到拱坝有限元分析时建模过程中的便利性，本文采用便于建模的平缝模型来模拟横缝，采用附加切向弹簧的方式来模拟键槽的咬合作用。

3.2 附加切向弹簧刚度取值稳定性的验证

为探究平缝模型附加切向弹簧刚度与单元特征尺寸(接触面单元竖向尺寸)的关系，改变图8(b)所示计算模型网格尺寸，由0.122 m×0.200 m改为0.244 m×0.200 m，0.367 m×0.200 m，横缝面竖向尺寸分别为0.122 m，0.244 m，0.367 m。通过调整附加切向弹簧刚度值，使得修改单元尺寸模型参考点相对位移与精细模型参考点相对位移的相对误差在10%以内，得出的附加切向弹簧刚度值分别为0.15 GPa/m(沿横缝面布置17个弹簧)、0.32 GPa/m(布置8个弹簧)、0.50 GPa/m(布置5个弹簧)。三种网格尺寸模型横缝面布置附加弹簧总刚度值分别为2.55 GPa/m，2.56 GPa/m，2.50 GPa/m，取平均值2.54 GPa/m作为模拟两个横缝键槽平缝模型的附加切向弹簧刚度值，则单个键槽咬合作用可由切向刚度值为1.27 GPa/m的附加弹簧替代，不同网格尺寸平缝模型与精细模型竖向相对位移对比，如图10所示。

3.3 三维平缝模型附加切向弹簧刚度系数取值

附加切向弹簧刚度系数是对实际键槽切向咬合作用在横缝接触面上的均匀化反映，与连接在横缝接触面对应节点上的集中弹簧相等效，其变形方向沿键槽横向，垂直于键槽法向。集中弹簧的附加弹簧刚度与共用节点单元面积有关。

在考虑横缝接触面法向指数模型和切向库伦摩擦模型的基础上，鉴于实际应用中横缝两侧单元网格剖分较一致的特点，利用两侧节点构造8节点单元，如图11，节点1～4和节点5～8分别分布在横缝的两个接触面上。局部坐标下该单元刚度矩阵中的元素可由式(13)求得

Kij=KttiAiδi(j-4)(i,j=1,…,4)

(13)

式中：Ktti为单位面积键槽模型咬合作用的等效附加切向弹簧刚度系数Ktti=1.27 GPa/m；Ai为节点1～4所围成的面积对节点i的贡献面积值；δij为克罗内克函数。

4 工程应用

考虑到空库工况为高拱坝横缝张开-闭合的控制工况，以锦屏一级高拱坝为研究对象，建立三维有限元分析模型。其中，采用考虑键槽咬合作用的动接触-联结混合模型模拟地震作用下横缝的非线性接触，研究键槽咬合作用对于拱坝地震作用下横缝张开-闭合-滑移的影响。

4.1 工程概况及计算模型

锦屏一级为双曲拱坝，坝顶高程1 885 m，建基面高程1 580 m，最大坝高305 m。电站正常蓄水位1 880 m，死水位1 800 m，拱冠梁顶厚16 m，拱冠梁底厚63 m，最大中心角93.12°，顶拱中心线弧长552.23 m，厚高比0.207，弧高比1.811。三维有限元分析模型沿坝体高度方向分为17层，沿坝体厚度方向分位5层，共离散6 105个单元，其中坝体单元2 205个，地基单元3 900个，单元均为八节点六面体单元。为突出研究重点，假定该拱坝设置3条横缝，将坝体分给4个坝段，横缝间距为136～142 m，横缝采用平缝模型，通过施加附加切向弹簧模拟键槽的咬合作用。有限元模型如图12所示，分缝位置如图13所示。坐标系采用笛卡尔坐标系，其中x轴正方向为横河向左岸指向右岸(上游面)，y轴正方向为顺河向，z轴正方向为竖直向。

计算模型材料参数如下：坝体不考虑材料分区，混凝土密度为2 400 kg/m3、静弹性模量为28.8 GPa、泊松比为0.167；为突出研究重点，地基按经典、实用的无质量地基模型考虑，静弹性模量为30.0 GPa、泊松比为0.2。动力计算时，拱坝结构阻尼采用瑞利阻尼，阻尼系数α=2ξω1ω2/(ω1+ω2)，β=2ξ/(ω1+ω2)，其中ξ为阻尼比，ω1和ω2为选取的模态频率上下限。取ξ=5%，结合前两阶自振频率得到空库拱坝的阻尼系数α=0.778 8 s-1，β=0.002 9 s。空库工况是横缝开度的控制工况，本模型不考虑库水的影响，计算荷载为坝体自质量和地震荷载。地震荷载采用边界输入的方式，输入地震荷载加速度时程曲线如图14所示，地震动峰值加速度为0.20g，时间步长为0.02 s，地震持时20 s。

横缝为单向键槽设置，模型计算考虑两种情况：①不考虑键槽咬合作用时，取KS=0；②考虑键槽咬合作用时，附加切向弹簧刚度系数取1.27 GPa/m，结合有限元网格每层坝体厚度和节点对横缝面竖向尺寸的贡献来取值。横缝接触面的非线性相互作用，其中，法向采用指数接触模型；切向采用库伦摩擦模型及考虑键槽咬合作用的附加刚度弹簧的点对联结单元模型。通过试算确定指数模型的c=5×10-7m，p0=5 GPa，可较好模拟横缝法向张开-闭合的力学行为；切向动摩擦因数f取0.6[30]、附加切向弹簧刚度系数为1.27 GPa/m，可较好模拟键槽咬合作用。取两种计算模型左、中、右三条横缝坝顶部位的法向开度与切向位移进行对比分析。

4.2 键槽咬合作用对横缝开度和切向位移的影响

地震作用下，上述两种工况下三条横缝法向开度和切向位移峰值，如表1所示。以中缝为例，上述两种工况下其上游坝顶参考点的法向开度和切向位移时程，如图15所示。

由表1可见，在地震动荷载作用下，考虑键槽咬合作用能增强分缝拱坝的整体性，减小横缝的法向开度，且能效限制横缝的切向位移。考虑键槽咬合作用后，左缝的法向开度峰值由2.609 cm减小至2.184 cm，降低幅度约16.3%；右缝的法向开度峰值由3.868 cm减小至3.003 cm，降低幅度约22.4%。如图15(a)所示，考虑键槽咬合作用，即以动接触-联结混合模型模拟横缝的非线性接触，中缝法向开度峰值降低约20.1%。如图15(b)、图15(c)所示，考虑键槽咬合作用后，中缝的切向位移大幅降低，中缝的切向位移限制在毫米级，该结论与江守燕等、Jiang等、苑举卫研究中的结论相一致。同样，考虑键槽咬合作用，对左、右缝有相同的作用效果。综上，地震荷载作用下，考虑键槽的咬合作用对有缝拱坝的横缝开度和切向位移响应均有较大影响。

表1 地震动作用下横缝法向开度和切向位移极值

5 结 论

本文根据键槽的传剪特征，建立了带键槽横缝的小试件有限元精细模型。分析结果表明：简化模型一和简化模型二与精细模型相比竖向位移的相对误差较大，而通过调节平缝模型的附加切向弹簧刚度系数可以实现与精细化模型精度相当的咬合作用。因此，推荐采用带附加切向弹簧刚度的平缝模型来模拟横缝键槽的咬合作用。

在锦屏一级拱坝的地震响应分析中，研究了键槽咬合作用对横缝的法向开度和切向位移的影响。研究表明，联结单元模型的引入使得横缝模型具备了模拟键槽咬合作用的能力。同时，联结单元的引入改善了仅以动接触力模型模拟动力作用下横缝的非线性接触中出现的两侧坝体互相嵌入。采用考虑键槽咬合作用的非线性接触横缝模型能够较好地模拟横缝拱坝的力学性能，并有效地限制横缝的切向位移，使得横缝的法向开度小幅降低。