基于有限元的重力坝抗滑稳定静动力可靠度快速求解方法

刘银勇 林潮宁 刘晓青 杜效鹄 周兴波

(1.河海大学 水利水电学院, 南京 210098;2.水电水利规划设计总院, 北京 100120)

重力坝的坝基抗滑稳定是重力坝设计的重要内容,也是重力坝工程能够安全运行的关键[1-2].目前,重力坝坝基抗滑稳定的分析方法主要分为确定性分析方法和可靠度分析方法[3].

确定性分析方法虽然在工程实践中得到非常广泛的应用,但其假设影响重力坝坝基失稳的因素都是确定性的,这与工程实际情况并不相符,因此存在一定的局限性.相比于确定性的分析方法,以概率统计理论为基础的可靠度分析方法,可以对影响重力坝坝基抗滑稳定的不确定性因素加以量化分析,故更加符合工程的实际情况.

可靠度分析方法在重力坝坝基抗滑稳定分析中应用广泛.江胜华等[4]根据加权响应面法拟合隐式功能函数,计算了重力坝坝基各断层滑动破坏的可靠指标;CAO 等[5]研究了失效准则、结构参数的随机性与模糊性,结合蒙特卡罗法给出了重力坝抗滑稳定分析的计算方法和程序;王晓玲等[6]提出了重力坝抗滑稳定分析的PLS-ELM 动态响应面法,采用经过偏最小二乘法优化的极限学习机模型构建响应面,通过蒙特卡罗法对响应面进行动态更新并得到可靠指标;CHEN 等[7]提出了基于加权回归的改进响应面法,进行了偏差系数的修正,并将其用于重力坝的坝基抗滑稳定分析;朱晓斌等[8]提出了基于猫群算法的加权动态响应面法,采用猫群算法取代了传统响应面法中应用广泛的蒙特卡罗法,提高了数据利用率.

通常,重力坝坝基抗滑稳定的功能函数为隐式的,需要通过有限元等数值分析方法进行求解,进而结合蒙特卡罗法或响应面法计算可靠度指标.这种方法在实际工程的应用中仍存在一些不足之处:①蒙特卡罗法模拟抽样次数多,需要反复调用有限元计算,计算成本高,计算效率低下,花费机时较多;②对于复杂的可靠度问题,由于传统响应面法拟合的隐式功能函数是高维高非线性的,计算可靠指标时容易出现拟合精度不高、实用性不强等[9]问题;采用智能算法模型构建响应面时,需要合理确定模型所需的一些最优初始参数,计算复杂程度较高.

针对上述分析方法在求解重力坝坝基抗滑稳定可靠指标时存在的问题,本文提出一种基于有限元的改进一次二阶矩法进行重力坝抗滑稳定可靠度计算.将坝基岩体抗剪断强度参数作为随机变量,定义功能函数为重力坝抗滑稳定极限平衡状态下阻滑力与滑动力的差值,在极限平衡状态条件下推导了功能函数中随机变量的偏导数,并基于改进一次二阶矩法编写了求解重力坝坝基抗滑稳定可靠指标的计算程序.以某重力坝工程为例开展抗滑稳定的静动力可靠度分析,采用该方法与常见可靠度分析方法对抗滑稳定静力可靠指标进行对比分析,验证本文方法在求解效率方面的可行性;采用振型分解反应谱法计算地震作用,考虑地震荷载的随机性,采用基于全概率公式的数值拟合积分方法对基于本文方法求得的特定地震作用下的抗滑稳定可靠指标进行概率统计分析,计算重力坝在动力作用下的抗滑稳定可靠度.

1 重力坝抗滑稳定静动力可靠度计算方法

1.1 抗滑稳定静力可靠度计算方法

由于力具有矢量性[10],将重力坝抗滑稳定的功能函数定义为总阻滑力矢在滑动趋势方向上投影的代数和R与总滑动力矢在此方向上投影的代数和T的差值.

式中:c为材料凝聚力;σn为滑动面上任一点法向应力值;φ为材料内摩擦角为滑动面上任一点抗滑剪应力单位方向为滑动面上任一点法向应力矢量为滑动趋势方向为滑动面上任一点应力矢量.

当重力坝处于抗滑稳定极限平衡状态时,功能函数是极限阻滑力与滑动力代数和的差值.将重力坝滑动面上的滑动力达到阻滑力作为重力坝抗滑稳定极限平衡状态的标志,据此,建立滑动面的极限状态方程为:

式中:n为滑动面单元总数;Wi、Pi分别为第i个滑动面单元上的法向力和切向力;fi、ci分别为第i个滑动面单元上的抗剪断摩擦系数和抗剪断凝聚力;si为第i个滑动面单元上扣除受拉区后的面积.

对随机变量求偏导数:

采用基于改进一次二阶矩法编写的抗滑稳定可靠指标计算程序进行重力坝抗滑稳定静力可靠度分析.该程序可以实现有限元的调用计算、验算点以及功能函数的迭代更新,其中改进一次二阶矩法的基本原理[11]如下:

设结构的极限状态方程为

再设x*=(,,…,)T为极限状态面上的一点,即

在点x*处将式(7)按Taylor级数展开并取至一次项,有

在随机变量X空间,方程ZL=0为过点x*处的极限状态面的切平面.利用相互独立正态分布随机变量线性组合的性质,ZL的均值和标准差分别为:

根据可靠指标的定义,可得结构的可靠指标

将式(9)对应的极限状态方程ZL=0用Xi的标准化变量Yi=(Xi-μX)/σX改写,并用式(11)作为法化因子遍除,整理后得:

定义变量Xi的灵敏度系数如下:

设计验算点在标准化正态变量Y空间中的坐标为:

在原始X空间中的坐标为

将式(8)、(12)、(14)和(16)联立可迭代求解β和x*.重力坝抗滑稳定的静力可靠度计算流程如图1所示.

1.2 考虑地震荷载随机性的动力可靠度计算方法

动力可靠度计算方法与静力可靠度计算方法具有相似性.动力可靠度本质上是先采用动力分析理论进行结构动力响应计算,再基于响应分析的结果结合概率统计理论进行概率统计分析,从而计算出结构在动力作用下的可靠指标.在结构动力可靠度分析中,普遍将水平地震加速度系数视为随机变量,考虑其变异性对可靠度的影响.水平地震加速度系数A是水平向地震动峰值加速度ah与重力加速度g的比值,即:

采用振型分解反应谱法计算地震作用,在特定地震动峰值加速度作用下,基于本文方法求得的是不同水平地震加速度系数条件下重力坝抗滑稳定动力可靠度,与此相对应的失效概率是条件失效概率P(f|A).在已知水平地震加速度系数的概率分布函数为f(A)的前提下,可根据全概率公式求得重力坝抗滑稳定的失效概率Pf以及相应的动力可靠度指标β:

式中:Φ(·)为标准正态分布函数.考虑到难以用一个准确的数学解析式表达P(f|A),故采取数值拟合积分的方法对上式进行计算[12].重力坝抗滑稳定的动力可靠度计算流程如图2所示.

2 工程应用

2.1 工程概述

某水电站大坝为碾压混凝土重力坝,坝顶高程153 m,最大坝高112 m,坝顶宽度6 m,建基面高程41 m,自左岸向右岸分为10个坝段.上游正常蓄水位高程150 m,相应下游尾水位高程41 m,泥沙淤积高程67 m.坝区地震设计烈度为8度.

2.2 计算模型

选取5号坝段的挡水坝段进行抗滑稳定静动力可靠度分析计算,计算原点取在坝踵处,x轴向右岸为正,y轴向下游为正,z轴为高程方向,向上为正.坝基范围向坝基底部、上下游分别延伸2倍坝高,基础宽为21 m.模型共剖分7195个节点,5328个单元,三维有限元模型如图3所示.

图3 某重力坝三维有限元模型

选择对重力坝抗滑稳定影响较大的坝基岩体抗剪断强度参数作为随机变量,认为其服从独立对数正态分布.各材料密度ρ、弹性模量E、泊松比υ、抗剪断摩擦系数f、抗剪断凝聚力c的均值见表1.其中岩体的f、c变异系数分别为0.25、0.3.

表1 某重力坝材料参数均值

2.3 静力可靠度计算分析

采用本文提出的基于有限元的改进一次二阶矩法对坝体沿建基面的抗滑稳定静力可靠度进行计算分析,考虑的荷载包括上下游水压力、坝体自重、坝前淤沙压力和坝基面扬压力,扬压力折减系数取0.25.通过非线性有限元法计算获得滑动面上单元的法向力和切向力,构建重力坝沿建基面的抗滑稳定功能函数,根据改进一次二阶矩法的迭代步骤进行可靠度的分析计算.采用编写的可靠度计算程序,经过5轮的迭代求解就能达到精度要求,计算效率高,收敛速度快,得到的失效概率为2.802 1×10-6,相应的可靠指标为4.540 8.

在结构可靠度分析中,蒙特卡罗法模拟的次数多,计算精度高,可将其计算结果视为标准解[13].为了验证本文方法的可行性,采用蒙特卡罗法进行对比分析.通过1×106次的蒙特卡罗模拟计算,得到该重力坝的失效概率为3×10-6,相应的可靠指标为4.526 4,计算时长约为1.7×105s.比较两者的计算结果,误差为0.32%,本文方法仅需5轮迭代求解就能达到相同精度,说明本文方法具有更高的计算效率.

除此之外,采用其他可靠度分析方法对该重力坝进行可靠度分析,计算结果见表2.根据本文构建的功能函数,采用中心点法进行可靠度计算,得到的可靠指标为4.119 9,与蒙特卡罗法计算结果相比,误差为8.98%.采用不含交叉项的二次多项式响应面法进行可靠度计算,计算复杂度较高,拟合功能函数构建响应面时需要调用9次有限元计算,通过27次有限元计算,迭代更新3次响应面,计算得到的可靠指标为4.453 2,比较其计算结果与蒙特卡罗法计算结果,误差为1.62%.研究结果表明,本文方法计算过程简易,在保证计算精度情况下具有更高的求解效率和收敛速度,且基于有限元法进行的计算,能够考虑坝体和基岩复杂的几何构造、材料分区的不均匀性及滑动面上真实的应力、变形分布情况.

表2 可靠指标计算结果统计值

2.4 动力可靠度计算分析

考虑的荷载包括坝体自重、上下游水压力、坝基面扬压力、坝前淤沙压力及地震荷载.将地震作用进行离散化处理,在特定的水平地震加速度作用下,考虑水平向和竖向地震的共同作用,竖向地震加速度取水平向的2/3.采用振型分解反应谱法计算地震作用,取最大地震加速度反应谱值为2.0,阻尼比为0.1,反应谱特征周期为0.2 s,并根据本文方法,求得不同水平地震加速度作用下的建基面抗滑稳定可靠指标,计算结果见表3.

表3 不同地震水平加速度系数下的动力响应

表3数据表明,在地震作用下,建基面抗滑稳定可靠指标随水平地震加速度系数的增大而减小,这与工程实际情况是一致的.

对表3的重力坝抗滑稳定失效概率与水平地震加速度系数的关系进行多项式回归拟合,拟合结果如图4所示,计算得到的拟合函数为

图4 重力坝失效概率与水平地震加速度系数的拟合关系

依据文献[14]的研究成果,水平地震加速度系数的概率密度函数f(A)和概率分布函数F(A)服从下列分布:

式中:B、b为待定系数,需要依照场地有关地震的统计资料确定.参照文献[15-17],暂取B=0.049 9,b=19.948.水平地震加速度系数的概率密度函数f(A)和概率分布函数F(A)如图5～6所示.

图5 水平地震加速度系数的概率密度函数PDF

图6 水平地震加速度系数的概率分布函数CDF

将式(20)和(21)代入式(18)并进行水平地震加速度系数为0.5的右截尾计算,求得在地震作用下坝体沿建基面抗滑稳定的失效概率为1.772 8×10-4,相应的可靠指标为3.571 8.

该重力坝在地震作用下的坝基抗滑稳定可靠指标与无地震作用下的可靠指标相比,下降幅度约为21.34%.依据该工程规模大小和规范[18]规定,静力工况下,重力坝坝基抗滑稳定目标可靠度指标建议值为3.7,动力工况下的目标可靠度指标可低于静力工况的目标可靠度指标.在静力工况下可靠指标均能满足规范要求,表明该坝段满足抗滑稳定要求.

3 结 论

基于有限元法的可靠度分析,往往无法构建显式的功能函数,仅能采用蒙特卡罗法和响应面法求解,但采用蒙特卡罗法求解可靠指标时计算效率低、采用响应面法时求解精度不高,针对此问题,本文提出了一种基于有限元的改进一次二阶矩法,结论如下:

1)在保证计算精度情况下改进一次二阶矩法具有更高的求解效率和更快的收敛速度的优点.

2)基于本文方法,求得特定地震作用下的抗滑稳定可靠指标,根据全概率公式,采取数值拟合并积分的方式对某重力坝工程的抗滑稳定动力可靠度进行分析的组合方法具有明确清晰的数学含义,可以为相关工程动力可靠度的求解提供方法借鉴.

3)对于重力坝工程来说,影响其坝基抗滑稳定的随机因素有很多,除了坝基岩体抗剪断强度参数外,其它参数对其抗滑稳定的影响需要进一步研究.