基于边界有限元计算方法的多相场重力式大坝动力特性分析研究

张宏美

(浙江湖州宏盈建设工程有限公司，浙江 湖州 313100 )

1 概 述

水资源是人类生存发展的基础，作为水资源的重要调控设施，水库的稳定性是保障水资源高效利用的重要基础[1-3]。水库的重要组成部分即为大坝，水库的安全运营与大坝安全稳定性息息相关。已有很多学者针对水库大坝的静力稳定性分析作过相当多的研究，极大提高了水库大坝的建设与设计水平[4-6]。大坝作为重要水利枢纽工程，抗震能力是不可避免需要考虑的重要方面。不同于静力分析，大坝动力响应是在多相条件下，包括坝基与流水状态，既有固体场亦有流体场。国内外已有一些学者通过振动台模型试验等分析大坝动力反应特征[7-10]；当然在试验条件不足情况下，上世纪已有一些学者通过解析方法，求解坝体地震反应特征[11-12]。上述两种研究手段，一定程度上均需要消耗大量的人力物力成本，而数值分析作为高效求解手段，针对于动力反应分析具有较高的匹配性。因而，基于多相场的边界有限元分析方法[13-15]，求解大坝动力反应特征，为准确评估大坝动力响应稳定性提供重要的理论参考。

2 二维边界有限元理论

2.1 基本方程

在动力分析理论中，常常会引入边界元分析理论，并结合有限元离散计算方法，提升求解速率与精确度，为解决工程解析问题提供重要手段。在二维边界有限元理论中，重定义边界比例，中心点已不是传统的零点，其几何示意图见图1。

图1 比例边界坐标体系

在比例边界坐标体系中，传统直角坐标体系中坐标可表述为[16]：

(1)

根据斜率求导原则，可获得偏微分方程为：

(2)

而边界有限元坐标体系中，边界点的坐标可表示为：

(3)

引入雅可比矩阵J，将上式以矩阵形式表述为：

(4)

上式即为边界有限元中运动矩阵方程，根据上式可演变出静力问题下的应力与变形解，弹性变形下的应力与变形方程如下：

(5)

(6)

而变形区域内几何方程为：

(7)

式中：[L]为微分算子式，其代表式如下：

(8)

根据有限元插值理论可知，比例边界坐标系中边界与中心点距离可用径向插值位移表述，其插值服从下式：

u(ξ,η)=u(ξ)N(η)

(9)

联立上述应力应变以及几何方程，可得坐标系中几何统一方程为：

(10)

式中:B1、B2代表系数矩阵，其具体代表式如下：

(11)

物理运动方程可统一表述为：

(12)

2.2 动力控制方程

引入能量做功理论来表述动力运动特征，能量参数的泛函式如下：

(13)

根据材料虚功原理可知，能量泛函式左边为0，即可得到：

(14)

各分项做功分别为：

(15)

经简化后可得到：

(16)

式中：系数矩阵均为

(17)

针对二维状态下的动力平衡方程中，引入惯性力来表述频域，其控制方程为：

E0ξ2u(ξ),ξξ+(E0+(E1)T-E1)ξu(ξ),ξ-E2u(ξ)+ξ2ω2M0u(ξ)+F(ξ)=0

(18)

式中有：

(19)

动力控制方程的求解需要引入对偶变量q：

q(ξ)=E0ξu(ξ),ξ+(E1)Tu(ξ)

(20)

通过矩阵转换以及特征值换算，获得控制方程解为：

(21)

式中：c为常数矩阵。

由于本文分析的是多相场作用下动力反应特征，因而划分单元网格时考虑多场单元体模型，其形函数可表述为：

(22)

根据形函数表述形式，获得动力特征中应力应变转换矩阵表达式为：

(23)

因而，可在上述基础上，获得比例边界坐标体系中刚度矩阵K与质量矩阵M分别为：

(24)

3 工程概况

某水库建设有混凝土重力坝，坝顶最高处达80 m，坝顶宽度为3.5 m，已修建有通行公路，上下游坡度均为0.75，水库最大库容可达1 200×104m3，承担着需水灌溉及工农业生活用水，在枯水期作为空间上水资源调配的重要利器。大坝包括有液压控制闸门，启闭度以及高度均由电脑系统精确控制，输送水资源量可精确到0.01 m3/s，闸门上下游分别设置有陡坎与消力池，减弱大水流对坝体冲蚀作用。坝体整体防渗性较高，通过在上游坝基覆盖层以上部分设置有厚60 cm的防渗墙，坝身上安置有止水面板以及土工格栅，降低水流对坝体渗漏破坏性。水库正常蓄水位设计为72.8 m，所有水位设计均按照大坝允许承受荷载计算。由于混凝土施工关乎着重力式大坝安全可靠性，因而采用BIM技术对大坝混凝土施工均进行过多次模拟，获得混凝土最佳振捣次数为2次，减少由于混凝土裂缝带来的坝体失稳与渗漏威胁。

根据地址钻孔资料得知，地区内主要为河漫滩平原，地势较为平坦，最高点为两侧的丘陵山坡，坝体两翼岩体处于较稳定状态，滑坡及坍塌等地质灾害并不显著，基岩层为白云质灰岩，在局部坝体区段处还含有夹砾岩，粒径为15～60 mm，完整性较好，地表出露主要在右侧坡体上。表面覆盖土层包括人工填土、粉质壤土与砂土，其中粉质壤土夹有粉、黏粒，含水量达37%，部分钻孔取样显示有结核，具有韧性表现，分布范围基本涵盖坝体所有区段，厚度约为2.2～6.9 m；人工填土主要为种植土，较薄，最厚处仅为2.5 m，夹有细砂结构，析水性较强，即渗透系数较低，分布并不均匀，密实性相差较大，部分地段密实度较高；砂土与卵石相交织，其中取样表明卵石含量约占该深度取样总重的51.0%～88.4%，砂土属全新统，分析表明卵石可能是由于水流搬运作用，从上游补充水源河流冲击搬运至该处，并持续沉积。根据进一步调查发现，由于区域内所处板块间，因而考虑地震动响应特征很有必要。

4 多相作用下动力分析

4.1 模型建立

以坝体平面状态作为建模基础，坝基作为刚性变形体，以SOLID65作为其基本单元体模型，坝体采用PLANE42模型，利用ANSYS建立数值模型。图2为大坝桩号1+150～1+160区段的数值模型图，该区段划分单元体网格，共获得18 368个，节点数16 538个，动荷载通过引入地震荷载开展分析。

4.2 计算结果分析

为分析方便，本文给出地震动响应特征解析解与边界有限元计算结果对比，为评价边界有限元方法的科学合理性提供重要依据。以坝体-水压力、坝体-坝基的两相场以及坝体-水压力-坝基三相场分别进行计算分析，根据相场特点，与图2建模类似，分别取其中坝体段开展计算分析。

图2 数值模型图

4.2.1 坝体-流水两相场

通过施加地震动荷载后，获得各阶荷载下自振频率。本文以一阶自振频率下不同反应系数下动水压力的变化曲线作为对比分析参数，其计算结果见图3，其中一阶自振频率为32.31 rad/s。

图3 不同反应系数下动水压力的变化曲线

从图3可看出，反应系数为0.25时，曲线相关系数达到0.99，即边界有限元计算结果与解析解基本一致，动水压力整体上呈从坝踵至坝顶，动水压力逐渐降低，在坝踵处，基于边界有限元方法计算出的动水压力最大可达到0.026 5 MPa，与高程0.8H处时相比来说下降52.8%；当反应系数为0.5时，坝踵处动水压力为0.037 MPa，相比反应系数0.25时增大39.6%，反应系数0.5时坝体高程上动水压力分布基本与小反应系数时基本一致；当反应系数增大至0.75，甚至更大的0.925时，动水压力在高程上分布趋势并未发生较大改变，且坝踵处动水压力分布呈现随反应系数增大逐渐增大，其中反应系数0.925时坝踵动水压力相比反应系数0.25时增大326.4%，达到0.113 MPa。

图4为一组激励频率下坝体动水压力分布曲线。

图4 一组激励频率下坝体动水压力分布曲线

从图4中可看出，最大动水压力曲线出现在激励频率与自振频率相等时，即坝体出现自振现象，动水压力在坝体上分布处于最高，且坝踵处为动水压力中分布最大值。由于坝踵处动水压力分布最大，给出不同激励频率下坝踵处动水压力变化曲线，见图5。

图5 不同激励频率下坝踵处动水压力变化曲线

从图5中可看出，随着激励频率变化，动水压力整体变化呈先增大后逐渐降低，其中最大动水压力出现在一阶自振频率时，即频率5 Hz下动水压力最大，达0.058 MPa。在该频率后，坝踵动水压力逐渐降低，但当降低至二阶自振频率时，会出现陡降现象，即频率15 Hz时，下降幅度超过70%，且在二阶自振频率后动水压力持续处于较低水平。综上分析表明，两相场条件下，坝体动水压力分布受地震动荷载频率与动力反应系数影响，坝踵处动水压力乃坝体上最大值。

4.2.2 坝体-水压力-坝基三相场

前文计算结果为两相场条件下坝体动力反应特征，但坝体在一定程度上能否抗震，很大程度上取决于坝基的稳定性，故而将坝基纳入分析动力反应特性中很有必要。即构建起三相场条件下动力反应体系，坝体区段取图2所示的数值模型，相关参数按照土工试验结果给定。

为分析方便，本文给出大坝两个特征点处X、Y相位移时程曲线，两个点分别是坝顶与下游坡度变化点，见图6和图7。

图6 坝顶位移时程曲线

图7 下游变坡点位移时程曲线

从图6和图7中可看出，坝顶处X向位移值整体高于Y向位移，在时间为第二个峰值点时，即时间1.2 s,X向位移达到0.066 m，Y向位移仅为前者的30.3%；同理，在谷值点时亦是如此，在第二个谷值点处，Y向位移达到负方向的0.017 m，而X方向相比增长2.2倍，达到0.054 m。从X、Y向位移整体来看，各振荡次数下的峰值点位移、谷值点位移之间的周期基本相近，且振荡表现逐渐趋于稳定，即边界有限元所计算出来的位移时程结果整体是稳定的。相比坝顶处位移结果，下游变坡点位移整体表现与之一致，不论是X向亦或是Y向，均处于较为稳定状态；从位移量值来看，下游变坡点位移值亦是X向位移高于Y向位移，在第二个峰值位移点时，X向位移高出Y向位移4.5倍，且谷值点处亦是相差近5倍，即变坡点处X向与Y向位移间差距相比坝顶处更大。

5 结 论

针对重力式大坝动力响应特性问题，引入边界有限元动力分析方法，计算获得两相场与三相场条件下大坝动力特征，结论如下：

1) 基于边界有限元动力理论计算的结果与解析解基本一致，不论动力反应系数如何变化，两者结果均一致，相关系数均超过0.99，边界有限元理论对多相场重力坝求解具有较高的自适应性。

2) 一阶自振频率下动水压力从坝踵至坝顶逐渐降低，反应系数为0.25时坝踵动水压力最大，达0.026 5 MPa，而高程0.8 H时降低52.8%；当动力反应系数增大时，坝踵处最大动水压力逐渐增大，反应系数0.925时坝踵动水压力相比反应系数0.25时增大326.4%。

3) 当激励频率与自振频率一致时，坝体动水压力最高；随激励频率增大，坝踵处动水压力呈先增大后逐渐降低，其中最高点位于一阶自振频率时，达0.058 MPa，二阶自振频率时动水压力会出现陡降，降低幅度超过70%。

4) 三相场条件下X、Y向位移周期与相位均稳定，且坝顶与下游变坡点均呈X向位移高于Y向位移，第二个峰值点时，坝顶处X向位移达到0.066 m，Y向位移仅为前者的30.3%，但下游变坡点两者差距4.5倍，X、Y向位移间差距比坝顶处更大。