基于COMSOL Multiphysics的某水库均质坝渗流与稳定性模拟分析研究

李宏伟

(乌鲁木齐市水利勘测设计院(有限责任公司)，乌鲁木齐 830001)

1 概 述

我国水资源分布面积较广[1-3]，如何有效利用水资源是水利工作者夜以继日思考的问题。为此，水利工程技术人员考虑修建水库大坝等水工建筑物，用以水资源调度[4-5]。这些水工建筑物给人类的生存发展提供重要推动力，但其安全性一直是业内研究考虑的重要方面。已有设计人员与学者基于室内水工模型试验，设计安全稳定性较高的水工建筑，提出水库大坝类水工结构安全设计施工的要点[6-9]。还有一些学者借助数值有限元软件，模拟计算水工建筑物在长期运营荷载下安全稳定性，准确预判坝体失稳[10-12]。不论采用何种手段，其目的均是为了保证水工建筑物的安全状态。针对具体的工程案例，采用数值模拟软件，综合分析获得渗流与结构稳定性，为工程施工提供重要的理论参考。

2 渗流与稳定性分析理论

2.1 渗流理论

渗流基本理论方程为达西公式，其表达式如下：

(1)

式中流速v与水力坡降J成正比关系。

基于达西基本理论，联系岩土材料各向异性特点，获得各方向上达西渗流微分表达式：

(2)

针对不可压缩流体，获得容纳黏滞参数的应力方程为：

(3)

而在多孔介质中渗流运动与孔隙度有关，运动方程可表述为：

(4)

图1 单元体中渗流示意图

基于质量守恒定律，单元体水质量变化曲线为：

(5)

式中：n为介质孔隙率；ρ为流体密度；V为单元体体积。

联系渗流量公式，并考虑流体不可压缩性，可得到：

(6)

式(6)即为刚度无限大的介质中流体运动微分方程，任意时刻任意点的流量值为零。基于二维流体方程，推广至三维维度内，稳定渗流微分方程为：

(7)

当介质为各向同性时，微分方程基于Laplace变换，可得：

(8)

若流体运动处于非稳定状态，则渗流运动势必需考虑介质状态与含水层，运动微分方程可表述为：

(9)

当介质材料为各项同性时，式(9)可简化得到：

(10)

不论是稳定渗流或是非稳定，运动微分方程的求解主要基于初始条件与边界条件，获得符合条件的收敛值，有限元计算收敛值即是依靠计算机强大计算能力，获得精度解。在前述分析基础上，渗流运动方程基本表达式为：

(11)

式中：h=h(x，y，z)为水头函数关系式；Ss为蓄水量。

基于变分与等参量变换，获得单元渗流矩阵表达式：

(12)

联系渗流基本微分方程与单元体水质量本构，等参变换有：

(13)

利用高斯积分与简化的雅可比矩阵，积分解为：

(14)

式中：φ(ξ，η，ζ)为积分函数在高斯积分点处的最高精度函数值。

通过微分方程解与渗流边界条件联系，即可获得渗流场中渗流速度、水头压力等特征参数值。

2.2 稳定性理论

COMSOL Multiphysics有限元软件计算平台计算水库坝坡稳定性可采用强度折减法，以边坡失稳和有限元计算不收敛同时出现，作为解精度要求的依据。强度折减法本质上是以M-C强度准则为内在依据，公式如下：

(15)

其中：

(16)

基于有限元软件计算坝坡直至不收敛时，岩土体抗剪参数临界值及安全稳定系数作为分析重要基础，为评价坝坡安全稳定性提供参量。

3 工程概况

某水库为当地蓄水调度重要枢纽，满足当地工农业基本发展需求，并在旱季提供水资源输送渠道，灌溉农田面积超过2.667×104hm2，水库库容为500×104m3，有一条苗河提供水源，该河流流域面积超过10 km2，河流淤泥沉积悬浮沙过多时，影响到下游水力发电后，该水库可作为反调节水利工程，引水调节河流水质。蓄水工程包括泄洪道、涵闸室输水孔与消力池等水工建筑，水库正常水位设计为920.3 m，洪水位为921.3 m，泄洪道等水利设施均为Ⅳ级建筑，调度水平均流量为1.2 m3/s。坝体主轴线设计修建为与河道正交，设计高程为922 m，坝顶宽度为2.5 m，上下游坝体坡度为0.5和0.4，坝身属均质土。

现场地质踏勘表明，坝址场地内地质构造属平静状态，无较活跃断层带与断裂带，但可见部分出露地表岩石节理发育较丰富，夹有风化碎屑岩层。依据钻孔资料得知，场地内覆盖层包括人工填土与粉质黏土等，人工填土以种植土为主，韧性较差，夹有第四系河流搬运作用形成的砾石，粒径为2～10 mm，磨圆度较高，由于经过压实作用，含水量较低，仅为15%。区域内地形地貌以河漫滩冲积平原为特点，两侧岸坡依附在山体中，山体坡度为25°～41°，岩层以半风化花岗岩为主，中粗粒结构。基岩为砂砾岩，下伏另有二长花岗岩与之不整合接触，标准承载力超过300 kPa。为保证坝体安全稳定，水库沿线岸坡均已完成清淤清障，提升坝基承载力。

4 均质坝渗流与稳定模拟计算

4.1 模型建立与研究工况

均质坝体断面图见图2(a)。本文将基于断面形态与岩土材料建立几何模型，后以sat格式形式导入至COMSOL计算平台中，几何模型见图2(b)。

图2 断面图与几何模型

以单元网格划分几何模型，获得三角形单元网格为主的数值模型，共划分出单元网格2 512个，节点数2 165个，单元质量最小为0.45，见图3(a)。图3(b)为模型中单元质量直方图，直方值较高区域集中于后半部分，即坝体渗流与稳定性可在求解之时获得较高精度解。

图3 数值模型与直方图

所有边界荷载参数均参照实际工程勘察资料报告，边界约束荷载与工况有关，渗流场以正常蓄水位、水位回落期开展；稳定性研究工况与坝体修建以及水库运营有关，故以稳定渗流期、水位回落期及坝体施工期开展计算分析。

4.2 渗流分析

水位回落曲线与时间关系曲线见图4。在COMSOL中将定义荷载与时间相关函数，分析水位回落曲线特征可知，水位陡降位于第24 d，速率为8 m/d。笔者认为水位回落后半阶段属骤降工况，应考虑对上游坝坡渗流与稳定性影响。

图4 水位回落曲线与时间关系曲线

基于上述分析分别获得正常蓄水位与水位回落期两工况下渗流场特征参数解，见图5、图6。从图5、图6中可看出，正常蓄水位稳定渗流状态下浸润面并未蔓延至坝顶区域，迎水侧坝身大部分区域渗透压力均为正值；从渗透水压力分布云图亦可看出，下游侧背水侧及坝踵等区域渗透压力为负值以下，最小为-0.2 MPa，最大渗透水压力为1.4 MPa，位于上游迎水侧坝基上覆盖层等区域，坝坡与下游侧渗透水压力差符合渗流特征，初步判断上游坝坡安全。结合坝体渗透水压力矢量分布图可知，渗透方向在坝体内部由上游侧指向下游侧坝基，另基于COMSOL后处理分析获得坝体渗流水头等势线，可看出上游侧逐渐递减至下游侧，坝顶中心线等势值为50%，水头递减规律显著。正常蓄水位坝体最大水力坡降值为2.89，位于底部沟槽处，浸润线上部起始点处水力坡降值为0.08。综合正常蓄水位各特征参数值可认为该均值坝体在正常蓄水位下渗流处于安全状态，坝体不会出现渗透破坏等现象。

图5 浸润线(正常蓄水位)

图6 渗流特征参数分布图(正常蓄水位)

为研究水位回落期渗流场特征，给出坝体各个时间段内水力坡降在坝坡高程上分布，见图7，图7中每一条曲线即为一个时间段。分析水力坡降在水位回落过程中表现可知，非稳定渗流场中水力坡降值随水位回落，逐渐增大，在第12 d～18 d每下降1 m高程，水力坡降增大0.009；当第18 d～25 d，每下降1 m高程，水力坡降增大0.034，速率增大近1个量级，即水力坡降增大速率与降落速度有关。

图7 各时间段水力坡降在坝体上分布曲线

在t=27 d后，水力坡降值超过0.53，即越过安全临界水力坡降值，坝坡渗流安全性受到挑战，为此利用COMSOL提取出t=27 d后渗流特征参数分布，见图8、图9。

图8 浸润线(水位回落期)

图9 渗流特征参数分布图(水位回落期)

从图8、图9中可分析看出，浸润线最高点高程为917.5 m，逼近正常蓄水位高程；底部沟槽处渗透水压力接近于零，最大渗透水压力位于上游侧坝基处，但坝体内部渗透水压力分布较大，下游侧坝基处渗透水压力由增大趋势。由渗流矢量方向可看出，渗流方向指向分为两类：上游侧坝坡与坝体内部，即上游侧坝坡渗透压力会在一定程度促进坝坡渗流运动，导致坝坡渗漏量增加或渗透破坏。最大水力坡降值为5.61，处于上游侧坝坡内，相比正常蓄水位下增大94.1%，水位回落过程中非稳定渗流场各特征参数种种表现表明，该工况下坝坡内渗流活动处于危险区间，应考虑针对上游坝坡增加防渗措施结构，在坝身表面铺设土工防渗膜，合理调控水库水位回落，保证坝体渗流安全。

4.3 稳定性分析

稳定性分析按照前述理论基础，借助COMSOL有限元软件开展正常蓄水位、水位回落(正常蓄水位回落至死水期)、施工期3个工况下坝坡稳定性分析，获得3种工况下坝坡失稳滑弧云图，见图10。其中，水位回落工况共采用有效应力强度折减法与最小应力组合两种计算方案。正常蓄水期与施工期坝坡模拟计算出的滑弧均出现在下游侧，而水位回落期滑弧出现在上游侧，表明水位回落期坝坡失稳首先出现在迎水侧。经计算获得3种工况下安全稳定系数，见表1，3种工况下安全系数分别为1.54、1.33和1.73。对比规范要求临界最小安全系数，3种工况下均处于安全状态。

图10 坝坡失稳滑弧云图

表1 3种工况下安全稳定系数

5 结 论

引入三维渗流分析理论与强度折减方法计算稳定性，基于某水库均质坝体工程资料分析，借助COMSOL Multiphysics数值软件建立分析模型，研究不同工况下坝体渗流与稳定性，结论如下：

1) 获得了正常蓄水位下坝坡的上下游侧渗透压力差分布、浸润线、水头等势线均由上游侧递减至下游侧特征，最大水力坡降值为2.89，位于底部沟槽处，渗流处于安全，坝体不会出现渗透破坏。

2) 水位回落期中水力坡降随水位回落逐渐增大，且增大速率与降落水位有关；浸润线最高点高程为917.5 m，渗流方向指向上游侧坝坡与坝体内部；最大水力坡降值为5.61，处于上游侧坝坡内，坝体渗流处于危险状态，应考虑布设防渗结构。

3) 分析了正常蓄水位、水位回落期、施工期3个工况下坝坡危险滑弧，得到3种工况下安全系数分别为1.54、1.33和1.73，均高于临界最小安全系数，坝体稳定性处于安全状态。