中国是水资源大国，总量丰富但时空分布极不均衡，其中西部地区水资源较为丰富。混凝土坝大多建设在多发地震的西北及西南地区，地质条件十分复杂。多次大地震中，大坝的震害使其抗震安全成为研究者关注的重点之一。混凝土坝的地震响应十分复杂，迄今为止，特别是混凝土高坝，遭受强震的震害实例很少，对混凝土坝的抗震性态仍在不断探索当中^[1]。坝体在强地震荷载作用下，易产生裂缝并极易扩展，造成坝体结构失效破坏^[2]。因此，分析混凝土坝在强地震荷载作用下的动力响应、损伤破坏过程及失效概率十分重要。

混凝土重力坝在强地震动下的安全性和可靠度分析至今依旧是一项难题。近年来，国内外学者在混凝土坝地震动输入和损伤破坏分析^[2]方面进行了全面研究，取得了阶段性的突破。目前，水电工程水工建筑物抗震设计规范》（NB35047—2015）^[3]指出，在地震作用下常用的响应计算方法有拟静力法、振型分解反应谱法和时程分析法。由于坝体在不同地震动强度下的非线性时程分析需要反复调幅计算，计算量较大且容易离散，由伊朗学者Estekanchi首次提出的耐震时程法（endurance time analysis，ETA）可以解决这一问题。Estekanchi、Valamanesh等^[4–6]详细介绍了耐震时程加速度曲线（ETA时程）的优化合成过程，并研究坝体损伤情况及ETA在线性、非线性分析计算中的可行性。Hariri-Ardebili、Riahi^[7–8]等以钢框架为模型，提出评价指标，研究并比较ETA法与增量动力分析法（incremental dynamic analysis，IDA）的抗震性能。白久林等^[9]合成了基于中国抗震反应谱的耐震时程加速度曲线，对钢框架结构的抗震性能进行分析和评估，并与IDA法的结论对比，得出ETA法能很好地预测结构的非线性动力响应及破坏概率。李杰和陈建兵等^[10–14]首次提出概率密度演化法（probability density evolution，PDE），考虑参数的随机性，建立广义概率密度演化方程^[11]，并提出数值求解方法，得到结构响应的概率密度曲线及概率密度演化过程，有效地获得了响应全过程的概率信息^[13]，为多自由度体系非线性结构的响应提供了新途径，在结构动力响应与失效概率方面获得了较大的成功。

基于以上方法，作者采用Koyna坝的场地谱合成了20组耐震时程加速度曲线，并作用于Koyna重力坝基底；对Koyna混凝土重力坝进行非线性时程分析，得出了坝体的动力响应及损伤破坏程度。通过动力损伤评价指标的选取，根据概率守恒原理，建立广义概率密度演化模型，得到损伤评价指标的概率密度演化过程，并依据概率密度曲线计算得出坝体的失效概率以评估坝体的安全性。通过对混凝土重力坝进行抗震安全评价及失效概率分析，为结构在强地震动作用下的设计提供了可靠的依据。

1 耐震时程法

耐震时程法（ETA）是一种简单的动力推覆过程，当结构受到预先施加的增强的动力激励时，可预测结构全过程的地震响应和抗震性能^[5]。耐震时程法的本质是生成随时间增加、强度逐渐增加的耐震时程加速度曲线，即ETA时程；且不同周期下的反应谱与目标反应谱存在成倍数拟合关系。ETA法给出耐震时程和目标时间：在目标时间下，反应谱与预先定义的目标反应谱一致；在其他时间下，与耐震时程成倍数关系。ETA法作为一种新的加速度时程输入方式，可模拟不同地震动强度下的响应，具有一定的适用性。

ETA时程的生成过程如下：先生成平稳地震动时程，再将时程乘以优化系数，使各时刻的加速度反应谱尽可能地与目标谱进行拟合，任意时段内的目标加速度反应谱与该时段的持续时间t呈线性关系^[9]：

${S_{\!\!{\rm{aT}}}}(T,t) = \frac{t}{{{t_{{\rm{target}}}}}}{S_{\!\!{\rm{aC}}}}(T)$

(1)

式中， ${t_{{\rm{target}}}}$ 为目标时间， ${S_{\!\!{\rm{aT}}}}(T,t)$ 为0～ $t$ 时刻的目标加速度反应谱， ${S_{\!\!{\rm{aC}}}}(T)$ 为0～ ${t_{{\rm{target}}}}$ 时刻的目标反应谱。

根据位移谱与加速度谱的关系，耐震时程目标位移反应谱表示为：

${S_{\!\!{\rm{uT}}}}(T,t) = \frac{t}{{{t_{{\rm{target}}}}}}{S_{\!\!{\rm{aC}}}}(T) \times \frac{{{T^2}}}{{4{{\text{π}}^2}}}$

(2)

式中， ${S_{\!\!{\rm{uT}}}}(T,t)$ 为0～ $t$ 时刻的目标位移反应谱。

以上过程只是生成不同周期下与目标反应谱大致拟合的ETA时程，还需进一步优化加速度曲线上的每一点。采用式（3）并利用MATLAB软件编制程序对其进一步进行无约束优化，使不同时刻的加速度反应谱值尽可能地在目标反应谱周围振荡，且拟合效果好^[15]。

$\begin{aligned}[b] {\rm{Min }}\; F({a_{{g}}}) = \int\nolimits_0^{{T_{{\rm{max}}}}} {\int\nolimits_0^{{t_{{\rm{max}}}}} {\left\{ {{{\Big[ {{S_{\!\!{\rm{a}}}}(T,t) - {S_{\!\!{\rm{aT}}}}(T,t)} \Big]}^2} + } \right.} } \\ \qquad \qquad\;\;\;\;\; \left. {{{\left. {\alpha [{S_{\!\!{\rm{u}}}}(T,t) - {S_{\!\!{\rm{uT}}}}(T,t)} \right]}^2}} \right\}{\rm{d}}t{\rm{d}}T \end{aligned}$

(3)

式中， ${a_{{g}}}$ 为初始生成的地震动， ${T_{\max }}$ 为反应谱时间最大值， ${t_{\max }}$ 为地震动时间最大值， $\alpha $ 为权重系数， ${S_{\!\!{\rm{a}}}}(T,t)$ 为周期T下0～t时刻的加速度反应谱， ${S_{\!\!{\rm{u}}}}(T,t)$ 为周期T下0～t时刻的位移反应谱。权重系数在文中取0，只考虑加速度反应谱的影响。

基于Koyna坝的场地谱合成了20组20 s的ETA时程，由于对任意时间下的反应谱拟合计算量较大，故选取0～5、0～10、0～15、0～20 s共4个时间段拟合反应谱，说明生成的ETA时程符合要求，如图1所示。

采用的加速度时程为ETA时程，并不是真实的地震动，仅采用ETA方法表征不同峰值加速度，表征坝体在不同地震动强度下的动力响应与损伤信息。在对Koyna坝施加的ETA时程中，0～5 s加速度反应谱表征0.15g的峰值加速度，0～10 s加速度反应谱表征0.30g的峰值加速度，0～15 s加速度反应谱表征0.45g的峰值加速度，0～20 s加速度反应谱表征0.60g的峰值加速度。通过对一条ETA时程进行分析可知，表示了不同峰值加速度下的坝体响应在不同时刻（对应不同峰值加速度）对应不同的损伤分布与响应信息。

耐震时程法具有两方面的特性：其一是峰值加速度随时间持续增加，其二是不同时程下的反应谱与标准反应谱成倍数关系。也即：耐震时程法的峰值加速度与ETA的时刻有关，一条ETA时程包含了许多条反应谱特性；而传统时程分析法，需要反复调幅生成地震动进行非线性时程分析，计算量大。由图1（b）可知IDA法与ETA法不同，ETA法的优点在于不同时刻对应不同峰值加速度，可得到不同峰值加速度下的动力响应，且计算量小，便于分析。

2 概率密度演化理论

多自由度结构在地震作用下的动力方程可表示为：

${{M}} \ddot {{X}} {\rm{ + }}{{C}} \dot {{X}} + {{F}}({{X}}){\rm{ = - }}{{M}}\mathop {{{I}}{\ddot {{X}}_{{g}}} } $

(4)

式中， ${{M}}$ 为质量矩阵， ${{C}}$ 为阻尼矩阵， ${{F}}({{M}})$ 为非线性恢复力， $\mathop { \ddot {{{X}} _{{g}}} } $ 为地震动加速度， $ \ddot {{X}} $ 、 $ \dot {{X}} $ 、 ${{X}}$ 分别为加速度、速度、位移向量， ${{I}}{\rm{ = (1,1,}}\cdots{\rm{,1}}{{\rm{)}}^{\rm T}}$ 为单位列向量。

s维随机向量 ${{\varTheta}} $ 代表结构的基本参数， ${{\varTheta}} {\rm{ = (}}{{{\varTheta}} _1}{\rm{,}}$ ${{{\varTheta}} _2}{\rm{,}} \cdot \cdot \cdot {\rm{,}}{{{\varTheta}} _{{s}}}{{{)}}^{\rm T}}$ ，即结构加速度、速度、位移值等基本物理量。除这些基本物理量之外，还需要求解一些感兴趣的物理量，用 ${{Z}}(t)$ 表示， ${{Z}}(t){\rm{ = (}}{{\textit{Z}}_{\rm{1}}}{\rm{,}}{{\textit{Z}}_{\rm{2}}}{\rm{,}} \cdot \cdot \cdot {\rm{,}}{Z_m}{{\rm{)}}^{\rm T}}$ ，m为感兴趣物理量的数目，感兴趣的物理量都可用基本物理量表示^[11–12]：

${{Z}}(t) = {{H}}({{\varTheta}} ,t)$

(5)

${{Z}}(t)$ 对时间的导数 ${\dot {{Z}}(t)} $ 为广义速度向量，也可表示为 ${{\varTheta}} $ 的函数：

${\dot{{Z}}(t)} = {{h}}({{\varTheta}} ,t)$

(6)

根据概率守恒原理得到^[10]：

$\frac{{\partial {p_{{{\textit{z}} \rm{\theta }}}}({\textit{z}},\theta ,t)}}{{\partial t}} + \sum\limits_{i = 1}^m {\dot {{{{Z}}_i}} (} \theta ,t)\frac{{\partial {p_{{{\textit{z}}\rm{\theta }}}}({\textit{z}},\theta ,t)}}{{\partial {{\textit{z}}_i}}} = 0$

(7)

特别地，当m=1时，1维广义概率密度演化方程为：

$\frac{{\partial {p_{{{\textit{z}}\rm{\theta }}}}({\textit{z}},\theta ,t)}}{{\partial t}} + \dot {{Z}} (\theta ,t)\frac{\partial {p_{{{\textit{z}}\rm{\theta }}}}({\textit{z}},\theta ,t)}{{\partial {\textit{z}}}} = 0$

(8)

式中， ${p_{{{\textit{z}}\rm{\theta }}}}({\textit{z}},\theta ,t)$ 为 ${{Z}}(t)$ 与 ${{\varTheta}} $ 的联合概率密度函数。

基于给定的初始条件和边界条件，求解偏微分方程，可得到 ${{Z}}(t)$ 的概率密度函数^[11–12]：

${p_{\textit{z}}}({\textit{z}},t) = \int {{p_{{{\textit{z}}{{{\varTheta}} }}}}({\textit{z}},\theta ,t){\rm{d}}\theta } $

(9)

目前，对于一般的动力系统，联合初始条件求解偏微分方程，得到其显示解答很难。文献[15]中采用有限差分方法求解该偏微分方程，可获得其数值解，作者采用MOL算法进行求解^[16]。取每组ETA时程作为一个样本点，取多组ETA时程（20组）建立样本集合，且假定每组ETA时程的发生概率相同。

3 数值模拟

Koyna混凝土重力坝是坝体在地震动作用下遭到破坏的实例之一，可在一定程度上反映坝体的损伤扩展状况，许多学者对Koyna混凝土重力坝在地震作用下的动力响应进行了深入研究^[17–19]。故作者选该混凝土重力坝，采用本文方法进行概率密度演化、失效概率分析。

Koyna混凝土重力坝高103 m，坝顶宽14.8 m，坝底宽70 m，坝高66.5 m处下游坝面出现折坡。为使计算结果具有较高的精度，采用在下游折坡处和坝踵处加密的有限元模型，如图2所示。本构关系采用混凝土塑性损伤本构模型，混凝土材料弹性模量为31 GPa，泊松比为0.2，密度为2 643 kg/m³，膨胀角为36.31°，初始抗压强度为13 MPa，极限抗压强度为24.1 MPa，抗拉强度为2.9 MPa，瑞利阻尼系数为0、0.003 23，断裂能为200 N/m。

外荷载主要考虑重力、静水压力（静水压力水位为91.75 m）、动水压力及地震作用，动水压力采用Westergaard附加质量法计算，地震作用采用ETA法合成的水平向峰值加速度0.6g和竖直向峰值加速度0.4g的ETA时程。

3.1 基本步骤

基于概率密度演化理论，给出工程实例在耐震时程加速度曲线下的概率密度演化及失效概率分析过程，具体实施过程为：

1）建立坝体有限元模型，施加荷载，根据混凝土重力坝在重力、静水压力静力荷载下的响应特点选取损伤评价指标作为随机变量；

2）分别计算Koyna混凝土重力坝在多组ETA时程下（文中取20组）的动力响应，输出坝体在ETA时程下的响应结果，得到坝体各损伤评价指标及其时间导数；

3）结合各损伤评价指标的边界条件，通过步骤2）得到的各损伤评价指标的时间导数，基于概率密度演化理论，采用MATLAB编制程序求解其概率密度演化方程，得到各损伤评价指标的概率密度曲线及演化过程；

4）对各损伤指标下的多组概率密度演化过程结果取平均（假定每组ETA时程的发生概率相同），得到概率密度演化结果；

5）设定各损伤评价指标的不同标准，根据其概率密度曲线及演化过程求解对应损伤评价指标的失效概率，对平均结果均匀化处理后得出失效概率的时程曲线（对应不同的峰值加速度）。

3.2 耐震时程加速度曲线输入

采用水平向和竖直向Koyna地震动场地谱放大倍数曲线，合成20组（最大水平向峰值加速度为0.6g、竖直向峰值加速度为0.4g）ETA时程，如图3所示。

采用MATLAB编制无约束优化程序，利用式（3）实现此过程，无约束优化需要对时程上的每个点进行。合成20 s的ETA时程加速度曲线，即根据不同周期下的目标反应谱，需优化2 000个点，反复迭代生成ETA时程。ETA时程0～5、0～10、0～15、0～20 s时间段下的反应谱拟合曲线如图4所示。由图4可知，在0～5、0～10、0～15、0～20 s时间上的反应谱均能与目标反应谱较好地拟合，能反映峰值加速度变化的情况下加速度反应谱的匀速变化，此方法合成的耐震时程加速度具有较高的精度，可反映出反应谱相似状态下峰值加速度变化对结构响应的影响。

3.3 混凝土塑性损伤本构

混凝土材料在初期浇筑时便存在一些细小的微裂缝，在长期荷载作用下，裂纹逐渐扩展至失稳。采用混凝土塑性损伤模型模拟混凝土材料在强地震作用下的裂缝扩展情况，由于混凝土的抗压强度明显高于抗拉强度，因此只考虑混凝土的拉损伤。单轴循环加载条件下的混凝土软化过程曲线见图5^[20]。

图5中， $\sigma $ 、 $\varepsilon $ 分别为混凝土应力、应变， ${E_0}$ 为初始弹性模量， ${{{d}}_{\rm{t}}}$ 为拉损伤因子， ${G_{\rm{f}}}$ 为断裂能， ${f_{\rm{t}}}$ 为抗拉强度， ${\varepsilon _{\rm{t}}}$ 、 ${\varepsilon _{\rm{f}}}$ 分别为最大弹性、极限拉应变， ${\varepsilon _{\rm{p}}}$ 为等效塑性应变， ${l_{\rm{c}}}$ 为混凝土单元的特征长度。

3.4 响应分析

通过ABAQUS软件非线性时程分析得到坝体在不同地震动强度下的响应，图6为典型时刻的损伤响应。

由图6可知：蓝色为损伤因子小于0.3的区域，全过程有许多条不同程度的宽裂缝出现。6 s（对应水平峰值加速度0.18g）时，下游折坡处和坝踵处出现轻微损伤；11.35 s（对应水平峰值加速度0.340 5g）时，下游折坡处损伤逐渐扩展，裂缝完全贯穿，此时坝体即将进入失效状态；20 s（对应水平峰值加速度0.6g）时，损伤体积逐渐增大，坝头裂缝逐渐扩展并新增加许多贯穿型裂缝，坝踵处损伤逐渐增加，裂缝达到约1/2倍坝底长度。在ETA时程的作用下，随着峰值加速度随时间按一定比例增大，坝体的损伤区域逐渐扩展。

在地震动作用下，采用概率密度演化方法可获得ETA时程响应全过程的概率信息，采用范书立等^[21]提出的能量指标（塑性耗散能、损伤耗散能）和损伤指标（损伤体积）作为动力响应指标。其中：损伤体积采用总损伤体积和基准损伤体积两种标准；总损伤体积定义为损伤因子大于0时的损伤体积；基准损伤体积选取0.7作为临界限值，定义为坝体损伤值大于0.7的损伤体积。通过分析得出4种指标的时程曲线，如图7所示，用于概率密度演化的分析和求解；图8为归一化的指标时程曲线。

由图7、8可知：20组塑性耗散能、损伤耗散能、损伤体积和基准损伤体积均随时间逐渐增加；且对于相同指标，离散性较小，具有相似的统计规律，归一化的4种指标时程在ETA时程分析中具有相似的变化规律。

3.5 概率密度演化分析

采用概率密度演化方法分析4种指标的概率密度演化过程。图9（a）、（b）、（c）、（d）分别为Koyna混凝土重力坝在ETA时程下的塑性耗散能、损伤耗散能、损伤体积及基准损伤体积的概率密度演化过程。

理想状态下，随着样本点数量的增加，概率密度演化过程结果会趋近真实的概率密度分布，为说明概率密度演化过程的正确性，分析并比较了4种响应指标下的20个样本点得到的概率密度演化过程及正态分布概率密度函数的结果。图10～13分别为各指标在2种方法下的比较结果。

由图10～13可知：对于每个指标，随着时间增加，概率密度曲线峰值逐渐减小，且逐渐变宽。选取12、16、20 s典型时刻，其正态分布结果与概率密度演化结果相近。相对于正态分布只能反映概率密度函数峰度的变化而言，概率密度演化方法可同时反映概率密度函数的峰度与偏度演化过程，优于假定的正态分布结果。演化理论可宏观地看出Koyna混凝土重力坝在ETA时程下动力响应全过程的概率密度演化信息。

上述对生成的20组ETA时程的分析中，样本数量存在一定的局限性，随着样本数量增加，概率分布更能反映真实的概率分布形态。

4 失效概率

结构的失效概率是判定结构损伤程度的有效指标，对于结构的失效概率，应给出结构响应的界限值要求，小于这一限值为可靠部分，大于响应界限值为失效部分，直接进行积分即可得到失效概率^[22]。若假设 ${X_i}$ 为不同时间下的响应值，X_L为响应的界限值（文中特指破坏的界限值，即为破坏标准），则有：

${p_{\rm{f}}} = P({X_i} \ge {X_{\rm{L}}}) =\displaystyle\int_{{X_{\rm{L}}}}^{ + \infty } {{p_x}({X_i},t)} {\rm d}x$

(10)

值得注意， ${p_x}({X_i},t)$ 为概率密度函数。式（10）建立了概率密度与失效概率的关系，有效地将概率密度演化方法与失效概率结合，得到了较为精确的计算结果。采取式（10）计算失效概率，需要选取破坏标准X_L，从ABAQUS损伤云图中得到损伤因子。中等破坏定义为损伤因子大于0时，上下游损伤贯穿；在结构的动力响应中找到不同指标结果，求得不同损伤指标下的平均界限值，采用2种方法分别计算不同指标下的失效概率，得到了4种指标在中等破坏程度下的失效概率时程曲线，如图14所示。

由图14可知，在不同峰值加速度响应下，4种指标的失效概率曲线具有较小的离散程度，塑性耗散能的失效概率值大于其他三者，这说明损伤耗散能、损伤体积与基准损伤体积作为大坝结构的中等破坏的失效概率评价指标差异性不大，具有一定的等效关系。

当损伤因子大于0.7时，一般认为已经发生了宏观裂缝，此时坝体损伤十分严重，故定义为严重破坏，上下游损伤发生贯穿。通过概率密度演化与正态分布方法的结果，分别得到不同指标的平均损伤界限值，根据概率密度曲线与失效概率的关系，得到超过这一界限值的失效概率，图15为严重破坏下的失效概率曲线。

由图15可知：失效概率随时间逐渐增加；随着时间增大，塑性耗散能与损伤耗散能的失效概率增长趋势大于损伤体积与基准损伤体积的失效概率增长趋势。这说明严重破坏条件下，各损伤评价指标有一定的离散型，塑性耗散与损伤耗散的失效概率可作为表征结构失效概率的上限，有一定的安全裕度。

通过比较图14、15中正态分布方法与概率密度演化方法可知，2种方法所得结果相似，验证了概率密度演化方法的正确性。同时，结果表明基于正态分布的计算结果偏于保守，高估了失效概率。

5 结　论

基于ETA法理论，通过MATLAB优化合成耐震时程加速度曲线，建立时间与峰值加速度的关系，满足不同周期下预先设定目标反应谱（场地谱）的要求。采用PDE理论，建立了广义概率密度演化方程，得到不同动力响应下的概率密度演化全过程及概率密度曲线，建立了时间与概率的关系，实现了坝体结构非线性动力响应的概率密度演化方法。

通过概率密度曲线计算得到失效概率，建立了时间与失效概率的关系。通过以上耐震时程法与概率密度演化法间接地建立了不同时刻峰值加速度与概率密度、失效概率的关系，有效地说明了本方法的适用性。以Koyna混凝土重力坝为模型，研究其结构动力响应与失效概率。研究结果表明：不同指标下，概率密度演化方法与正态分布方法结果近似，该方法具有较高的精度；在中等破坏下，塑性耗散能失效概率曲线增长趋势较为明显，损伤耗散能、损伤体积与基准损伤体积失效概率无明显差异；在严重破坏下，指标有一定的离散型，塑性耗散能与损伤耗散能可作为表征结构失效概率的上限，有一定的安全裕度，所以用能量指标更能保守估计宏观裂缝的扩展历程；概率密度演化理论具有较高的精度且能得出结构非线性响应的概率密度函数及其演化过程，失效概率计算结果比较稳定。

仅对生成的20组ETA时程进行分析，样本数量存在一定的局限性，使用PDE法的样本数量和真实概率分布的误差分析有待进一步研究。