胶凝砂砾石坝极端暴雨响应研究

2022-05-24张金凤郭利霞申伟平郭宇航

中国农村水利水电 2022年5期

郭磊，张金凤，郭利霞，申伟平，郭宇航

（1.华北水利水电大学水利学院，郑州 450046；2.河南水谷研究院，郑州 450046；3.河南省水环境模拟与治理重点实验室，郑州 450002）

0 引言

极端暴雨使河道水量增多，河水流速变大，会在短时间内形成较大的洪峰流量，使库水位急剧上升。目前，国内外关于水库水位变化对混凝土坝的影响分析已有很多，BAYRAK T［1，2］采用一种新的动态变形分析方法，研究了亚穆拉大坝的位移和库水位之间的关系，结果表明库水位变化对坝体变形影响显著。PYTHAROULI S I 等［3］利用谱分析技术对拉东大坝30 余年的大地检测记录进行分析，发现坝顶中部的垂直位移和水平位移最大，水位荷载与坝体变形为非线性关系。SUN K M 等［4］采用强度折减方法对混凝土坝进行了稳定分析，结果表明随着水位荷载和降雨强度的增大，坝体内应力应变变化也越大，越容易发生失稳破坏。GUO L等［5］通过数值仿真手段研究了混凝土坝对极端暴雨洪水的响应，研究表明，水温变化对坝体应力影响较小，而水位变化对坝体应力影响显著。以上研究均是对静动荷载下的混凝土大坝受力特性进行分析。本文以胶凝砂砾石坝为研究对象，对极端暴雨作用下水库水位的影响进行分析，并在有限元分析中引入了胶凝砂砾石材料的弹塑性损伤模型，分析了暴雨历时过程中坝体的应力、位移和损伤破坏规律，可为应对气候变化的结构设计提供理论依据。

1 计算理论

1.1 流量-水位关系

降雨径流是一个复杂的过程，降雨通过产流计算扣除损失量之后剩下净雨，再通过汇流作用经河网汇集到流域河道形成断面流量（流域净流量），水库入库流量通常是通过测量库前河道断面流速经计算得到的［6-8］。通过河道的全断面流量Q为：

式中：A为河道断面面积，m2；dA为A内的单元面积（其宽为db，高为dh），m2；v(h，b)为垂直于dA的流速，m/s；B为水面宽度，m；h为水深，m；hb为水边到水面宽为b处的水深，m。

由上可知，水位-流量存在非线性关系，但受河道的断面形态、弯曲状态、底坡以及河床表面糙率等因素的影响，很难精确获取其数学模型［9，10］。可利用有限元流体软件模拟水流运动过程［11，12］，其一般微分方程为：

连续方程：

动量方程：

式中：ρ为水的密度；u为沿水流方向的水流流速；w为垂直方向上的水流流速；μ为动力黏度系数；p为压强；Su、Sw为广义源项；为雷诺应力项，通用形式为：

式中：μt为湍动系数；k为湍动能；ui为时均流速；δij克罗内克符号（当i=j时，δij= 1；当i≠j，δij= 0）。

1.2 考虑损伤的有限元结构分析

对于空间结构的有限元分析计算，是将建立的实体模型划分为有限个具有块状结构离散的实体单元，相邻单元之间通过节点以铰接的方式连接，利用相互作用的连续单元运用数学方法去模拟物体真实的物理系统，实现利用有限个未知量求得无限个未知量的物理系统的近似解。在空间问题中，弹性力学物理方程可表示为：

式中：εx、εy、εz分别为X、Y、Z方向上的垂直应变；σx、σy、σz分别为X、Y、Z方向上的垂直应力；γxy、γyz、γzx分别为X-Y、Y-Z、Z-X平面上的剪切应变；τxy、τyz、τzx分别为X-Y、Y-Z、Z-X平面上的剪切应力；E为弹性模量，可将E换成考虑损伤后的弹性模量=(1 -d)E，其中d为损伤参数，可根据试验得到。

2 计算模型及参数

2.1 有限元模型

选取土耳其的Oyuk 坝为研究对象，该坝位于在峡谷之间，两岸坝肩陡峭，坝基以片麻岩为主。坝体高度为100 m，坝顶宽度为7.5 m，上、下游坡比为1∶0.7。在建立有限元模型时，为了提高仿真计算的可靠度，需适当扩大计算区域，坝踵处的基岩向上游延伸100 m，坝趾处的基岩向下游延伸100 m，坝底基岩竖直向下延伸100 m。采用SOLID185 单元划分网格，模型共有56 896 个节点，计算单元共51 060 个，有限元模型如图1 所示。应力计算时在地基底部施加全约束，在地基四周施加法向约束。

图1 坝体有限元模型Fig.1 Finite element model of dam body

2.2 计算参数

依据工程资料对坝体胶凝砂砾石材料进行了室内试验，坝基为岩体，参数选取参考文献［13］，胶凝砂砾石坝仿真模型计算参数见表1所示。

表1 计算参数Tab.1 Calculation parameters

胶凝砂砾石是一种弹塑性材料［14］，选取《混凝土结构设计规范》［15］的分段曲线模型作为胶凝砂砾石材料本构模型，通过单轴抗拉和单轴抗压试验得到其损伤演化公式为：

受压损伤公式为：

式中：dc为胶凝砂砾石单轴受压损伤演化参数；x为胶凝砂砾石压应变与其单轴抗压强度代表值相应的峰值压应变的比值。

受拉损伤公式为：

式中：dt为胶凝砂砾石单轴受拉损伤演化参数；x为胶凝砂砾石拉应变与其单轴抗拉强度代表值相应的峰值拉应变的比值。

2.3 水位边界条件的获取

由于极端强降雨天气的发生，导致水库水位上升较快，水库蓄满后泄水建筑物无法宣泄多余洪水，使得水位持续升高出现漫顶甚至溃坝事故。本文选取当地某次典型极端暴雨洪水过程为研究依据，因缺乏该库区实测资料，故假设库区相对规则，入库洪水过程线如图2所示。

图2 入库净流速过程线Fig.2 Process line of net inflow velocity

为了获取极端暴雨最不利工况下的水位变化特征，利用Fluent 软件采用VOF（Volume of Fluid，VOF）法和标准k-ε 模型对坝前水流进行模拟。模型使用非结构化四边形网格划分，共44 583 个网格。流体域从坝体分别向上游、下游延伸150 m，竖直向上延伸50 m。边界条件设置：模型底部为默认的固体壁面边界，上部设置为压力型入口；模型上游分别设置速度型入口和压力型入口，坝前库水位为97 m，断面平均流速设置为0.68 m/s，气压为一个工程大气压，下游设置为压力型出口。计算模型如图3 所示，对流场进行初始化，初始速度均为0 m/s，计算时长为100 s。

图3 流体域网格剖分示意图Fig.3 Schematic Diagram of Grid Subdivision in Fluid Domain

根据仿真计算结果绘制出时间-水位变化曲线，随着汇流作用水位呈现先急速后缓慢上升趋势，最高水位达到105 m，高出坝体5 m。

由图4拟合出暴雨工况下时间-水位关系式为：

图4 暴雨工况下库水位变化曲线图Fig.4 Changes of water level under heavy rain conditions

3 计算结果分析

应力计算采用超载法，将漫顶水位荷载完全等效为作用在坝体上的水荷载。极端降雨工况下坝体应力计算结果如图5～10所示。

图5 胶凝砂砾石坝不同历时应力云图（t=0 s，单位：Pa）Fig.5 Stress nephogram of CSG dam with different durations（t=0 s）

图6 胶凝砂砾石坝不同历时应力云图（t=20 s，单位：Pa）Fig.6 Stress nephogram of CSG dam with different durations（t=20 s）

图7 胶凝砂砾石坝不同历时应力云图（t=40 s，单位：Pa）Fig.7 Stress nephogram of CSG dam with different durations（t=40 s）

图8 胶凝砂砾石坝不同历时应力云图（t=60 s，单位：Pa）Fig.8 Stress nephogram of CSG dam with different durations（t=60 s）

图9 胶凝砂砾石坝不同历时应力云图（t=80 s，单位：Pa）Fig.9 Stress nephogram of CSG dam with different durations（t=80 s）

图10 胶凝砂砾石坝不同历时应力云图（t=100 s，单位：Pa）Fig.10 Stress nephogram of CSG dam with different durations（t=100 s）

（1）应力响应分析。由以上计算结果可知，在上游水位荷载的作用下，坝踵处所受拉应力最大，坝趾所受压应力最大。由于坝体与坝基材料不同，因变形不协调致使二者接触面较为薄弱，存在较大的不利应力。随着水位的上升，坝趾处压应力有所增大，而坝踵处拉应力明显降低了很多，这是因为坝踵区域处的材料发生了损伤，并且随着损伤的累计，现有破坏部位便会在应力集中的作用下，沿着顺水流方向和坝体纵深方向进行进一步的扩展，最后形成较大的宏观裂缝致使建筑物失效。故在极端暴雨工况下，胶凝砂砾石坝可能会在坝踵处发生损伤破坏，且为拉裂破坏；坝体与坝基接触面处存在较大不利应力，但影响不大。

（2）位移响应分析。胶凝砂砾石坝顺水流方向上最大变形面的位移量随暴雨历时变化如图11所示。

图11 暴雨工况下胶凝砂砾石坝顺水流方向位移图Fig.11 Displacement diagram of downstream flow direction on top of CSG dam under rainstorm condition

由图11 可知，坝体位移量随着上游水位的上升而逐渐增大，其增长速率由大变小，增长趋势与坝体上游水位荷载增长趋势相同，为正相关性。从大坝整体来看，坝体中部位移量较大，坝顶次之，坝底位移量最小。

（3）损伤单元分析。为了更好表述坝体的受力状况和损伤演化规律，提取坝体若干损伤特征单元进行分析对比，特征单元的具体位置详见图1，其损伤变化如图12所示。

图12 暴雨工况下胶凝砂砾石坝特征点的弹性模量损伤图Fig.12 Elastic modulus damage diagram of characteristic points of CSG dam under rainstorm conditions

由图12 分析可知，随着上游水位不断的升高，作用在上游坝坡的荷载逐渐增大，内部的39 961 单元无损伤，其他特征单元均产生损伤，且内部39 751 单元的损伤值小于表层39 301 单元，可见坝体开始由无损状态向损伤状态过渡，损伤范围逐渐扩大。由于损伤不断累积，特征单元的弹性模量有减小的趋势，其趋势由大逐渐变小直至弹性模量为零，损伤量达到极值，即特征点处结构完全失效，坝体材料被破坏。单元损伤值随时间的变化同时也验证了坝体应力的变化。总之，随着水位的升高坝踵部位的损伤量会优先达到极值，坝踵区域处的浅层地基会优先发生塑性屈服，塑性屈服区会分别沿着竖直方向向下扩展和建基面方向向下游扩展。