坝前滑坡涌浪压力分布对震损拱坝损伤影响的研究

李玉倩，陈健云，2，徐 强，李 静

(1．大连理工大学 工程抗震研究所，辽宁 大连 116024；2．大连理工大学 海岸与近海工程国家重点实验室，辽宁 大连 116024)

1 研究背景

大型水库的库岸滑坡涌浪可能会产生严重的灾害影响。例如，发生在1963年的Vajont滑坡涌浪漫坝灾难，约2.7亿m3的滑坡体滑入水库中产生了约270 m的巨大涌浪，造成了超过2000名的居民丧生[1]。Vajont大坝是在完好状态下承受到巨大涌浪的冲击，幸存的大坝在一定程度上降低了洪水和泥石流灾害。如果是大坝在强震后已经损伤情况下再次经历滑坡涌浪冲击而造成的溃坝，那后果将更加不堪设想。强震后受损边坡在相同的诱因下发生滑坡的可能性更高，比如2008年的汶川地震引发了数万次的山体滑坡[2]，其中紫坪铺水库的滑坡产生的涌浪高约25 m，造成了70多人死亡[3]。

强震后震损的山体，在余震或暴雨等诱发因素下比震前更容易滑坡[4-5]，从而对震后受损的大坝形成威胁。大坝在强震作用下会发生不同程度的破坏，再次遭受滑坡涌浪冲击将会进一步加剧损伤，一旦发生溃坝，将会对下游造成毁灭性的灾难。

在高拱坝的安全评价研究中大多仅考虑了地震的影响，如黄会宝等[6]基于大坝震后的监测原观数据，对大坝在地震前后的变形和横缝开度等进行了对比分析，为震后大坝安全评估提供了参考；秦礼君等[7]以白鹤滩拱坝为例，提出了基于抗滑安全系数及安全系数持时的高拱坝抗震安全评价指标；梁辉等[8]对高拱坝-地基系统整体进行安全评价，构建了概率地震风险分析模型；范书立等[9]考虑了地震动和材料的不确定性，建立了概率地震需求模型，并基于向量地震动强度指标进行拱坝的地震易损性分析；李静等[10]基于地震动IDA分析方法，建立高拱坝非线性数值分析模型，提出了大坝抗震性能评价指标，即坝体损伤体积比和坝面损伤面积比。

一些学者研究了滑坡涌浪对大坝的影响：黄锦林等[11]通过物理模型实验得到了坝前的涌浪荷载计算模型，并对乐昌峡大坝进行了安全评估；卢国帅[12]基于黄锦林的涌浪压力模型，探讨了地震滑坡涌浪作用下Koyna重力坝的损伤情况；李静等[13]采用SPH方法研究了坝面涌浪冲击压力模式以及涌浪对挡水重力坝的影响。

根据当前现状，滑坡涌浪对坝体的影响研究大多是针对重力坝的，对于震后滑坡涌浪冲击下的震损拱坝研究还处于空白。然而在我国西南地区一系列高拱坝已相继建成，由于西南地区属于高烈度地震区，且坝址地质条件复杂，库区存在许多大的潜在滑坡，地震后极有可能导致大规模的滑坡失稳入水激起巨大涌浪。因此，研究震后滑坡涌浪对震损拱坝的影响很有必要，对于能更好地保障大坝安全具有重要意义。

近年来，国内外学者开始关注滑坡涌浪传播至坝前的涌浪特性。邓成进等[14]通过三维物理模型实验得到了坝前不同测点的涌浪高，从浪高随时间的变化曲线可以看出，坝面不同位置的首浪高度是不同的，且涌浪首波到达坝面的时间也有差别；华璐等[15]探究了近坝库区滑坡涌浪的漫坝特性，发现坝前首浪振幅沿坝体两岸向中间减小；Xu等[16]基于SPH-DEM耦合方法再现了Vajont滑坡涌浪，得到了坝面不同位置的涌浪高度的变化，不同位置的涌浪漫顶高度是不同的。基于以上研究可见，滑坡涌浪传播到坝前的涌浪高和到达时间是不同的。然而，现有研究通常假定涌浪为同时到达坝前的等高荷载，尚未考虑坝前涌浪空间分布和时间特性对大坝的影响。

据此，本文以某高拱坝为例，首先，建立了动力分析的数值仿真模型，以模拟地震动对大坝的影响；其次，基于滑坡涌浪浪高估算的中国水科院经验公式[17]，分析了坝前滑坡涌浪高沿不同坝段的分布规律，并建立了坝前涌浪压力的计算模型；进而，从滑坡涌浪对震损拱坝造成的损伤范围和考虑损伤程度的损伤范围两个方面，建立了坝前滑坡涌浪对震损拱坝损伤影响的度量指标；最后，依据数值仿真的计算结果，详细分析了坝前滑坡涌浪压力的不同分布形式对震损拱坝损伤的影响。

2 动力分析数值仿真模型

某高拱坝为混凝土双曲拱坝，最大坝高289 m，坝顶高程为834 m，坝前正常蓄水位高程为825 m。本文基于商业有限元软件ABAQUS建立坝体-地基有限元模型，如图1(a)所示，拱坝设置了30条横缝，共包含31个坝段，编号如图1(b)所示。该模型共包括38 426个实体单元和56 417个节点，为了更精确地模拟坝体损伤，有限元网格最小尺寸为2 m左右，将坝体细化为25 776个实体单元和39 410个节点。

图1 拱坝-地基的三维有限元模型及混凝土应力-应变关系

选取混凝土塑性损伤本构模型(CDP模型)[18]模拟混凝土拱坝的材料非线性。塑性定义参数如表1所示。由于混凝土的抗压强度明显高于抗拉强度，因此本文只考虑混凝土的拉伸损伤。混凝土在单轴拉伸下的应力-应变关系如图1(c)所示。其中ε和σ分别为混凝土应变和应力；ft为抗拉强度；E0为初始弹性模量；dt为拉伸损伤因子，εp、εt和εf分别为等效塑性应变、最大弹性应变和极限拉伸应变，Gf为混凝土断裂能，lc为混凝土特征长度，通常定义为最大骨料粒径的三倍，本文取lc为0.45 m[19]。坝体和地基的材料参数如表2所示。采用三维黏弹性人工边界[20]模拟无限地基的辐射阻尼。

表1 CDP模型中的塑性定义参数

表2 拱坝-地基的材料属性

在荷载施加方面，静荷载包括大坝自重和静水压力，动荷载包括地震作用和震后的滑坡涌浪冲击压力。在动力分析中，动水压力采用Westergaard附加质量法[21]进行模拟。场地反应谱及地震波时程采用中国水利水电科学研究院给出的研究结果[22]，其中地震波是根据场地反应谱人工生成的。100年超越概率2%的设防水准基岩水平峰值加速度为0.414g，竖向峰值加速度为水平向的2/3。图2(a)为归一化的三向地震波加速度时程。本文采用基于黏弹性边界的地震波输入方法[23]模拟拱坝的地震反应。本文假设滑坡入水点与拱坝的相对位置如图2(b)所示。

图2 归一化的地震波加速度时程及滑坡的相对位置

3 坝前滑坡涌浪压力分布模型

3.1 坝前涌浪高的分布规律分析国内外关于滑坡涌浪的浪高研究方法主要有理论与经验估算法、物理模型试验法和数值模拟法。物理模型试验大多基于小尺度的试验，在实际工程中的适用性有限；数值模拟方法建模计算耗时久，效率低，成本高。理论与经验估算法计算简单、快捷、成本低，工程应用广泛。目前常用的经验估算方法有美国土木工程学会推荐方法[24]、Node法[25]、潘家铮方法[26]、中国水科院经验公式法[17]等。其中中国水科院经验公式法是根据大量的涌浪试验资料和原型观测成果总结得到的经验公式，易于计算且包含远场特征，被广泛采用[27-30]。因此，本文采用中国水科院经验公式法初步分析坝前的涌浪高分布规律。

中国水科院经验公式中认为滑坡涌浪的高度与距滑坡体的距离、滑坡速度和滑坡体的体积有关。距滑坡体不同距离的涌浪高度公式为：

(1)

(2)

式中：H为距滑坡体L处的滑坡涌浪高度，m；v为滑坡速度，m/s；U为滑坡体入水体积，m3；g为重力加速度，g=9.8 m/s2；n为计算系数，n=1.3～1.5；k1为与距离L有关的影响系数。

假设滑坡从左岸滑入水中，如图2(b)所示，入水点到31个坝段的距离分别记为L1，L2，…，L31，从而可以计算出31个坝段的涌浪高度。

鉴于本文主要目的是研究不同的坝前滑坡涌浪分布规律对震损拱坝坝体损伤的影响。为了具有可比性，经过前期的模拟试验，本文假定31个坝段的坝前最大涌浪高相同，均为100 m。相应的中国水科院经验公式中各个参数的计算范围设定为：L1=0.5～10 km，v=10～100 m/s，U=1×107～2×108m3。在计算范围内，满足坝前最大涌浪高为100 m的滑坡距离、滑坡速度和入水体积组合如图3(a)所示。

图3 所有滑坡情况下的不同坝段的涌浪高和到达相对时间

本文近似采用孤立波速度方程来计算滑坡涌浪的波速c[31]，公式如下：

(3)

式中h0为水深，本文h0=280 m。

根据图3(a)中所有的滑坡距离、滑坡速度和体积组合，可以计算得到不同坝段的涌浪高，如图3(b)所示；根据滑坡入水点到不同坝段的距离，可以计算出涌浪传播到不同坝段的相对时间，如图3(c)所示。可以看出，不同距离和规模的滑坡传播到31个坝段时的涌浪高的分布规律有所区别。当滑坡距拱坝的水平距离L1为0.5 km时，涌浪高的差值可达到30 m，滑坡离拱坝越远，涌浪高的差值越小，当L1大于2.5 km时，不同坝段的涌浪高几乎是相同的。对于相对到达时间来说，L1为0.5 km时，不同坝段的时间差最大，滑坡离拱坝越远，不同坝段的时间差越小。

将31个坝段涌浪高度差ΔH和时间差Δt的定义如下：

ΔH=Hmax-Hmin；Δt=tmax-tmin

(4)

式中：Hmax和Hmin分别为31个坝段涌浪高度的最大值和最小值；tmax和tmin分别为涌浪首波到达31个坝段所需的最长时间和最短时间。

分析其规律发现，31个坝段的涌浪高度差ΔH和时间差Δt与L1密切相关，即与滑坡入水点到坝前的水平距离有关，如图4所示。随着L1的增大，ΔH和Δt都呈现出先迅速减小再接近平稳的趋势。本文假定涌浪高度差ΔH≤5 m时，可以认为31个坝段的涌浪高度是近似相同的，均以H=100 m考虑；假定涌浪时间差Δt≤3 s时，不考虑其时间差，即认为此时涌浪首波同时到达31个坝段。因此，可以根据滑坡入水点距坝址的水平距离分为三种情况：

图4 涌浪高度差和时间差与水平距离L1的关系

LA：0.5 km≤L1<2.5 km，此时31个坝段同时存在涌浪高度差和时间差；

LC：L1≥6.5 km，此时31个坝段可以认为是同时作用均匀的涌浪高度。

基于上述分析，本文选出3种典型震后的滑坡涌浪工况，如表3所示。

表3 三种典型震后滑坡涌浪工况

3.2 坝前涌浪压力的计算模型震后的滑坡涌浪传播至坝前产生爬高，此时坝体承受的水压力由静水压力和滑坡涌浪产生的动水压力组成。在滑坡涌浪对结构的冲击压力分布研究方面，黄锦林[32]通过物理模型试验提出了一种坝前的涌浪压力计算模型。Chen等[33]进行了冰川入水产生冲击波的物理试验，提出了作用在大坝上的最大压力荷载经验公式。Tan等[34]通过物理模型实验测得了滑坡对岸的涌浪压力分布模型。Xu等[16]基于耦合的SPH-DEM算法，模拟了滑坡涌浪从产生到冲击大坝的全过程，并得到了大坝上游面的水压力。其中，黄锦林得出的压力模型与Xu等[16]模拟得到的动水压力分布基本一致，且较为合理。因此，本文采用此涌浪压力模型，坝前沿高程的涌浪压力分布如图5(a)所示，图中ϒ表示水的容重。在静水位以下，坝体上游面的涌浪压力沿高程是相同的，大小为1.9倍的涌浪高水头对应的浪压力；在静水位以上，1.9倍浪高范围内的涌浪压力采用静水压力分布形式。

图5 沿高程的涌浪压力分布模型及坝前涌浪高的作用形式

坝前涌浪高是随着时间不断变化的，依据以往的滑坡涌浪数值模拟和物理模型实验中得到的坝前涌浪变化曲线[14，15，32，35-37]可知，滑坡涌浪产生的首波对大坝的动力作用最大，当滑坡涌浪的首浪冲击到大坝时，坝前的涌浪高和压力都呈现先增加后减小的趋势，且曲线接近于半正弦的模式。因此本文仅分析涌浪首波对大坝的影响，近似采用半正弦波的形式模拟坝前涌浪高的变化过程。通过对比本文与以往研究中的滑坡规模及产生的涌浪高，本文拟取涌浪首波的作用时间为15 s，如图5(b)所示。

滑坡涌浪压力作用于地震动之后，为了保证地震动作用后拱坝能够达到稳定状态，在地震动作用结束后和涌浪冲击作用之间设置间隔20 s。根据表3中三种典型的震后滑坡涌浪工况计算坝前涌浪压力作用时程，如图6所示。工况1中不同坝段的最大涌浪压力和涌浪到达时间不同，工况2中只有涌浪到达时间不同，工况3中涌浪同时到达相同的涌浪压力。

图6 不同坝段的涌浪压力作用时程

4 坝前滑坡涌浪对震损拱坝损伤影响的度量指标

本文分别采用三种损伤影响指标，以度量滑坡涌浪对震损拱坝造成的损伤范围和考虑损伤程度的损伤范围。

(1)仅考虑损伤范围时，采用损伤体积比(RV)和损伤面积比(RA)作为指标，定义如下：

(5)

(2)分析损伤范围的同时，考虑每个单元的损伤程度，以基于损伤加权的损伤体积比(WRV)和损伤面积比(WRA)作为指标，定义如下：

(6)

第三，加强人民群众的法治观念和法律意识。人民是国家的主人，是依法治国的主体，人民群众法律水平的高低直接影响着依法治国的进程。广大人民群众必须不断提高法律意识和法治观念，自觉的遵守法律，依法维护自己的合法权益，运用法律手段同违法犯罪行为作斗争。

(3)为定量描述震后滑坡涌浪对震损拱坝的损伤增量，定义无量纲损伤比增量指标φRV，φWRV和φRA，φWRA：

(7)

5 计算结果分析与讨论

图7为拱坝在地震后和涌浪后的损伤分布图，图8为三种涌浪工况后相对震后新增的损伤分布图。由图7、图8可知，地震对坝体造成的损伤主要分布在坝体下游面的中上部以及建基面。震后涌浪冲击下坝体的损伤进一步发展，坝体上游面左侧坝段，以及各坝段的建基面损伤增加较大。对比三种工况的结果可以发现，涌浪压力分布中存在高度差和时间差(工况1)时，坝体增加的损伤最多，涌浪压力分布仅考虑时间差时(工况2)次之，均匀分布的涌浪压力(工况3)下坝体增加的损伤较少。

图7 震后和涌浪后的损伤分布图

图8 三种涌浪工况后的新增损伤分布图

通过水压力超载法[38]来评价拱坝在地震和涌浪后不同损伤状态下的坝体安全性。损伤坝段(坝段5)的坝顶顺河向位移随水压力超载系数K的变化如图9所示。由图可见，随着超载系数K的逐步增加，坝体顺河向位移也逐渐增大，K小于一定值时，曲线斜率保存一定，K超过一定值时，单位加载步长内的位移明显增加。拱坝在震后和不同涌浪工况后的损伤状态不同，对应的转折点的K也不一样。初步判断其K值处于2.0～3.0之间，在此范围内以0.1倍步长进行超载计算，得到不同损伤状态下的位移，明显增加点对应的K如表4所示。涌浪压力分布中存在高度差和时间差(工况1)时，最薄弱坝段增加的损伤最大，其K值最小，相比地震后其水压力超载能力降低了24%；均匀分布的涌浪压力(工况3)下坝体增加的损伤最小，相比地震后的水压力超载能力降低了14%。

表4 地震后和不同涌浪工况后的位移明显增加点对应的K

图9 坝顶节点顺河向位移-超载系数K的关系曲线

图10为不同坝段顶部节点的顺河向位移时程曲线。涌浪作用下新增的损伤主要在坝体上游面的左侧坝段，以及各坝段的建基面处，且考虑高度差和时间差时的涌浪压力分布(工况1)下的左侧坝段损伤增加最大，因此坝体的左侧坝段如坝段5和坝段10的顺河向位移，考虑高度差和时间差(工况1)时最大，仅考虑时间差(工况2)次之，均匀分布(工况3)最小；工况1时存在6.3 s的时间差，由于涌浪先传播至左侧坝段，坝面的力是不均匀的，坝段26在涌浪作用的初始阶段会发生向上游的位移。分析损伤坝段(坝段5)的坝顶节点的最大顺河向位移，考虑高度差和时间差(工况1)时最大顺河向位移为0.155 m，仅考虑时间差(工况2)时最大顺河向位移为0.137 m，均匀分布(工况3)时最大顺河向位移为0.137 m。与均匀工况相比，工况1增加了13%，工况2仅考虑了3.8 s的时间差，对最大顺河向位移几乎没有影响。

图10 不同坝段顶点的顺河向位移时程曲线

表5为地震作用和震后滑坡涌浪作用下的坝体损伤比和损伤比增量。由表可知：

表5 地震后和不同工况滑坡涌浪后的损伤比及损伤比增量 单位：%

(1)地震作用下和震后涌浪作用下坝体的下游面和建基面的损伤范围较大；震后涌浪增加的损伤范围主要在坝体上游面和建基面处，上游面损伤比增量较大，滑坡距坝体水平距离为0.5 km时能达到89.05%；震后涌浪对坝体下游面的影响很小。

(2)地震作用下和震后涌浪作用下，坝体建基面的损伤程度较大；震后涌浪对建基面处的损伤程度增大较多，尽管震后涌浪对坝体上游面的损伤范围增大较多，但是损伤值增加较小，上游面增加的损伤程度没有建基面处的大。

(3)对比震后三种滑坡涌浪工况，从损伤范围和考虑损伤程度的损伤范围两个角度分析，三种涌浪压力分布下的结果都呈现出相同的规律，即同时存在高度差和时间差的涌浪压力分布的结果(工况1)最大；仅存在时间差的涌浪压力分布的结果(工况2)次之；均匀分布的涌浪压力分布的结果(工况3)最小。究其原因：拱坝作为近似对称的三维空间壳体结构，在设计时所考虑的水荷载在同一高程是相同的，工况3中的均匀涌浪压力分布与其相似，拱效应发挥比较充分；而在工况1中，31个坝段的涌浪存在31 m的高度差、存在6.3 s的时间差，即对于同一作用时刻，同一高程不同坝面的受力相差较大，坝面上的涌浪压力是偏心作用的力，且根据图6的涌浪压力作用时程和图10的位移时程曲线可以看出，涌浪作用力是先作用于坝体左侧的，且坝体左侧的压力比右侧的压力要大，此时拱效应发挥作用相对于均匀分布(工况3)时要小，导致损伤坝段即左侧的变形比均匀情况下的大，上游面左侧坝段，以及建基面的损伤增加较多。

(4)现有研究中通常采用均匀分布的涌浪压力计算涌浪对坝体的影响，然而根据(3)的分析，均布荷载结果偏危险，低估了涌浪对坝体损伤的影响。因此，根据本文的计算内容，给出了当采用均匀的涌浪压力分布进行分析震后涌浪对坝体损伤时的修正系数，在考虑坝前涌浪压力分布中的高度差和时间差的影响时(工况1)，采用修正系数α修正均匀分布(工况3)下的损伤增量；在仅考虑时间差的影响时(工况2)，采用修正系数β修正均匀分布(工况3)下的损伤增量，公式如下：

(8)

表6为基于无量纲损伤比增量指标下的修正系数。可以看出，当滑坡距离坝体较近时即工况1，修正系数较大，且对于上游面损伤程度和建基面的损伤范围的修正系数较大。当涌浪压力分布中仅考虑3.8 s的时间差时，与均匀分布的情况相差不大，修正系数接近1。

表6 基于损伤比增量指标的修正系数

6 结论

本文以某高拱坝为例，基于中国水科院涌浪高经验公式，识别出了典型的坝前涌浪压力分布形式，借助数值仿真模拟的手段，详细研究了震后坝前滑坡涌浪压力分布对震损拱坝损伤的影响，主要结论如下：

(1)设定地震动后拱坝已发生比较明显的损伤，震后滑坡产生的涌浪最大浪高达100 m。在涌浪冲击下，坝体损伤会进一步发展，坝体上游面的损伤范围增加较大，在滑坡距拱坝水平距离为0.5 km时增加了89.05%；坝体建基面的损伤程度增加较多，在滑坡距拱坝水平距离为0.5 km时增加了30.02%，震后滑坡涌浪对震损坝体的下游面影响较小。

(2)当坝前滑坡涌浪的压力分布不同时，对震损拱坝损伤的影响是不同的；当坝前涌浪压力分布存在高度差和时间差时，震后涌浪对震损拱坝的影响比均匀分布的涌浪影响要大。

(3)均匀分布的滑坡涌浪压力分布模型，会低估滑坡涌浪对震损拱坝损伤的影响。建议研究震后滑坡涌浪冲击大坝时，考虑涌浪压力分布中存在的高度差和时间差的影响。

(4)基于本文提出的三种涌浪压力分布，讨论并给出了对均匀分布涌浪压力的结果修正系数，滑坡距离坝体越近，涌浪压力分布中的高度差和时间差越大，修正系数越大。

因此，建议考虑坝前涌浪的分布规律，以更加合理、安全地评估滑坡涌浪对震损拱坝损伤所带来的影响。同时，考虑到滑坡涌浪对拱坝影响的复杂性，基于本文现有的研究成果，仍然需要在滑坡涌浪的压力分布模型、涌浪作用时间、对大坝的影响度量指标以及影响规律、库区的潜在滑坡可能造成的涌浪对最大可信地震下大坝抗震安全性的影响等方面，展开更加深入的研究。