水库枢纽工程是保障流域水安全的重要战略措施, 在国民经济建设中起着重要作用。水库蓄水运行后, 泥沙落淤, 其淤积形态与水库自身特点、入库水沙条件以及水库运用方式等因素有关。而淤积形态亦会影响库容分布、水库输沙流态等, 进而影响水库设计功能的有效发挥。因此, 深入掌握水库淤积形态的形成机理及其模拟方法, 具有重要意义, 可为水库的科学调度提供技术支撑。

泥沙淤积是水库建设和运行需解决的一项重要关键问题^[1-4], 尤其是针对多沙河流的黄河小浪底水库枢纽工程等开展过大量研究。在模拟试验方面, 文献[5]采用准二维数学模型, 对小浪底水库运用初期库区淤积过程进行了数值模拟研究；李涛等^[6]采用小浪底水库实体模型, 开展了多个组次的水库降水库区发生溯源冲刷的试验, 并分析了水量、侵蚀基准面等因素与库区溯源冲刷量的关系；李书霞等^[7]改进了描述异重流运动的动量方程, 提出了新的异重流潜入条件判别式, 并用多组室内及野外实测资料对该判别条件进行率定与验证。随着水库的实际运用, 积累了丰富的原型资料, 为泥沙研究提供了最直接的素材, 张帅等^[8]分析了小浪底水库汛期排沙比与入库平均含沙量、进出库平均流量、平均库容等影响因素之间的相关程度, 确定了影响汛期平均排沙比的主要因素, 采用回归分析方法建立了拟合效果较好的排沙比公式；王婷等^[9]对小浪底水库蓄水运用以来(2000—2016年)库区泥沙淤积特征进行了统计分析；孙东坡和吴默溪^[10]基于原型观测资料, 分析了小浪底水库入库水沙、水库运用以及相应库区泥沙淤积的变化特征。张俊华等^[11]从异重流潜入与运动规律、干支流倒灌机制、异重流高效排沙试验、水库高效输沙调度原则等方面, 分析了小浪底水库异重流高效输沙理论与调控关键技术, 也指出了水库细颗粒浮泥层(淤积物)动力学机制的复杂性等需要进一步研究的相关问题^[12]。小浪底水库近坝段淤积泥沙粒径极细, 以往研究较少关注淤积物自身的流动特性, 而对于细颗粒泥沙淤积形态而言, 传统模拟结果与实际情况存在较明显的差异^[13]；针对细颗粒淤积物的流动特性问题, 假冬冬等^[14]认为三峡水库近坝区淤积形态的成因既与坝前水流流态的三维特性有关, 又与水库蓄水初期坝前淤积物颗粒特性有关, 水库蓄水初期沉积在近坝区的泥沙颗粒较细, 使得淤积物表现出一定的流动特性, 并在此基础上给出了概念性的临界坡降计算方法^[13]以及借鉴浅水流动方程来描述失稳后细颗粒淤积物重力驱动流动过程, 但并未对细颗粒淤积物的流变特性及其运动本构方程进行分析与模拟。

本文以黄河小浪底水库为例, 针对细颗粒泥沙淤积特点, 揭示水库细颗粒淤积物的流变特性与流型特征, 通过引入水、淤积物、床面之间的界面受力分析, 构建细颗粒淤积物失稳流动描述模式, 并与水沙输移模型相耦合, 建立考虑细颗粒淤积物流动特征的水库淤积形态模拟方法, 在此基础上, 对小浪底水库淤积形态进行模拟分析。研究成果可为水库淤积形态形成机理及其对水沙调控的响应研究提供技术支撑。

1   水库细颗粒淤积物流变特性与流型特征

小浪底水库近坝段淤积泥沙粒径极细, 中值粒径一般为0.008~0.015 mm^[10], 此类淤积物具有较强的流动性、淤积面平坦, 亦是影响水库淤积形态的一项重要因素。

1.1   流变特性

通过采集小浪底水库近坝段典型淤积物样品(级配如图 1所示, 中值粒径约0.01 mm), 采用数显流变仪进行不同密度条件下的细颗粒淤积物流变特性试验, 在此基础上分析流变特性参数屈服应力(τ_B)和黏度系数(μ)与密度(ρ_m)之间的关系。不同密度淤积物是采用一定量的细颗粒泥沙样品加入适当水体配备并搅拌均匀而成, 在1.08~1.35 g/cm³之间, 总共进行7种不同密度的试验, 随后点绘流变特性参数与密度的关系, 如图 2、图 3所示。由实测结果可见:水库细颗粒淤积物的τ_B和μ随淤泥密度的不同而发生变化, 当密度较小时, 流变参数随密度的变化较缓, 当密度较大时, 流变参数随密度的变化较快, 中间存在从缓变到急变的转化过程。从实测结果可看出, 当密度大于1.25 g/cm³后, 流变参数随密度的变化速率都快速增大。因此, 可认为当淤积物密度大于1.25 g/cm³后, 淤积物的流动性已快速减弱, 将不易流动。

图  1  小浪底水库典型细颗粒淤积物级配

Fig.  1  Typical sediment gradation in the Xiaolangdi Reservoir

下载: 全尺寸图片 幻灯片

图  2  淤积物屈服应力与密度之间关系

Fig.  2  Relationship between τ_B and ρ_m

下载: 全尺寸图片 幻灯片

图  3  黏度系数与密度之间关系

Fig.  3  Relationship between μ and ρ_m

下载: 全尺寸图片 幻灯片

根据van Kessel和Kranenburg^[15]的研究, 表征细颗粒淤积物流变特性的参数τ_B和μ可表示为:

式中: ${C_v} = \frac{{{\rho _m} - {\rho _f}}}{{{\rho _s} - {\rho _f}}}$为细颗粒淤积物的体积含沙量, ρ_m和ρ_f分别为细颗粒淤积物和水体的密度, ρ_s为泥沙容重, 一般可取值2.65 g/cm³；μ_w为清水黏度；c₁、c₂、c₃、c₄为参数, 根据本文试验数据率定, 分别取660、3.7、220和1.6, 具体对比结果如图 2、图 3所示, 其中屈服应力计算式与实测值吻合较好(R²=0.98), 黏度系数计算式与实测值吻合程度略差(R²=0.85)。

1.2   流型特征

为分析小浪底水库细颗粒淤积物的流型特征, 分别对不同密度的淤积物样品(1 082 kg/m³、1 112 kg/m³、1 150 kg/m³)在不同剪切速度下的μ进行了试验分析。研究结果如图 4所示, 随着密度的增加, μ逐渐增加；μ随剪切速度的变化较小, 随着剪切速度的增加, μ略有减小, 但在6~60 r/min的转速范围内, 其变化值一般在5%以内, 可基本认为在同一密度条件下μ接近一定值。密度较小的水库细颗粒淤积物在克服自身屈服应力后, 即可发生流动, 且μ随剪切速度的变化(较小)基本不变, 由此可认为未密实的低密度水库细颗粒淤积物为典型宾汉型流体。

图  4  黏度系数随剪切速度变化

Fig.  4  Relationship betweenμand shear rate

下载: 全尺寸图片 幻灯片

2   考虑细颗粒淤积物流动的水库淤积形态数值模拟方法

尚未密实的低密度水库细颗粒淤积物为典型宾汉型流体, 具有较明显的流动特性。因此, 对于由细颗粒泥沙特征决定的水库淤积形态而言, 除了考虑传统水沙输移过程外, 还需对已落淤的细颗粒淤积物自身流动过程进行模拟分析。关于水库泥沙输移的三维动力学模型, 可参见相关文献[13]。本文将结合细颗粒淤积物的流变特征, 对考虑细颗粒淤积物流动的水库淤积形态数值模拟方法进行探讨。

2.1   细颗粒淤积物失稳流动模式

前述分析表明，沉积历时较短的水库细颗粒淤积物属于典型宾汉体。其运动之前需克服自身屈服应力τ_B, 本文采用临界坡度作为细颗粒淤积物失稳流动的判别标准^[14]；通过引入水、淤积物、床面之间的界面受力分析, 建立细颗粒淤积物流动模式, 描述其失稳后的流动及其泥沙重新分配过程。

2.1.1   失稳流动判别

细颗粒淤积物具有一定临界坡度J_c, 则在J>J_c时(J_c=τ_B/γ_sΔz, τ_B为屈服应力, γ_s、Δz分别为细颗粒淤积物的容重和厚度), 淤积物将以其自重就可以沿着坡度方向运动^[16]。对于不同容重的细颗粒淤积物而言, 具有不同的屈服应力, 本文根据前述试验分析建立的细颗粒淤积物屈服应力的计算公式(式(1))进行失稳判别。即首先计算地形坡降J, 再根据时段内落淤细颗粒淤积物的厚度Δz以及其屈服应力τ_B(式(1))和容重γ_s来计算其临界坡降J_c, 地形坡降大于临界坡降(J>J_c)时, 即可判断细颗粒淤积物发生失稳流动, 失稳后的流动及其泥沙重新分配过程由下文流动方程控制。

2.1.2   流动方程

细颗粒淤积物失稳流动后, 由于细颗粒淤积物厚度相对较小, 流动过程具有近似平面水流的性质, 流场特性在深度方向的变化量远小于平面方向的变化量^[14], 通过沿深度平均简化可得细颗粒淤积物流动方程组。

连续方程:

通过引入水、淤积物、床面之间的界面受力分析, 且忽略扩散项时, 下部细颗粒淤积运动方程可表示为:

式中: Z为细颗粒淤积物表面高程；M=uh, N=vh, h为细颗粒淤积物厚度, u、v分别为细颗粒淤积物在x、y方向的流动速度；q为单位面积上细颗粒淤积物的源汇强度, 由各单元格失稳判别后计算得到, 即q=(Δz－Δz_c)/dt, Δz_c为可停留的临界厚度(Δz_c=τ_B/J_cγ_s)；g为有效重力加速度；τ_{m, bx}和τ_{m, by}为细颗粒淤积物和下部床面间切应力在x、y方向的分量；τ_{m, sx}和τ_{m, sy}为细颗粒淤积物和上部水体间切应力在x、y方向的分量。

上述细颗粒淤积物与下部床面以及上部水体之间的切应力由下式计算^[17] :

式中: u_f和v_f为淤积物上部水体的流速在x、y方向的分量；f_s和f_m分别为细颗粒淤积物与上层水体和下部床面之间的摩阻系数。

2.2   数值求解方法

与三维水沙输移模型一致, 细颗粒淤积物流动方程基于有限体积法进行离散, 该方法具有质量守恒的优点, 采用SIMPLEC算法进行模型求解。考虑细颗粒淤积物流动特性水库淤积形态的求解思路是通过水沙输移计算, 模拟水库冲淤分布, 在此基础上对细颗粒淤积物自身流动特性加以模拟, 获得淤积物的重新分布, 在此基础上再次求解水沙输移, 如此循环直至求解结束。对于细颗粒淤积物自身流动过程的模拟与求解与一般浅水方程求解类似, 前一计算时段计算获得的床面作为当前计算时段的边界床面, 当前计算时段内淤积物厚度h的初始值由本时段内计算的冲淤量Δz给定, 其求解亦是获得本时段内淤积物的重新分布。

3   小浪底水库淤积形态数值模拟

水库泥沙淤积形态与水库自身特点(大小、地形等)、入库水沙条件以及水库运用方式等因素有关。以黄河小浪底水库为例, 采用本文数值模拟方法, 根据实测地形资料建立数学模型, 并依据实际出入库水沙条件以及水库运行方式作为模型计算边界, 在此基础上, 对水库细颗粒泥沙淤积形态开展数值模拟验证分析。

3.1   小浪底水库概况

小浪底水利枢纽是一座以防洪(防凌)、减淤为主, 兼顾供水、灌溉、发电, 除害兴利的综合利用枢纽工程, 在黄河治理开发的总体布局中具有重要的战略地位。小浪底水库总体处于峡谷地带, 平面形态狭长弯曲, 入汇支流较多, 大支流与干流交接处多为开阔地带, 如图 5所示。小浪底水库大坝上距三门峡大坝130 km, 控制流域面积69.4万km², 占黄河流域面积的92.3%, 水库275 m高程原始库容127.5亿m³, 长期有效库容51.00亿m³。1994年9月水库主体工程开工, 1997年10月截流, 1999年10月下闸蓄水, 2000年5月正式投入运用。小浪底水库建成以后, 库区淤积了大量泥沙, 2015年4月与1999年9月相比, 库区累计淤积泥沙30.49亿m³(断面法计算结果), 约占设计拦沙库容的42.1%^[8]。

图  5  小浪底水库概况

Fig.  5  River regime of the Xiaolangdi Reservoir

下载: 全尺寸图片 幻灯片

3.2   模型计算范围及网格划分

模型计算范围上起HH56断面、下至小浪底水库枢纽, 如图 5、图 6所示, 模拟黄河干流河道长约123 km。模型采用贴体曲线非正交网格与复杂岸线贴合良好, 计算网格为1 800×240, 平均网格尺度沿水流方向约为70 m、沿断面方向约为10 m, 垂向网格为13层, 平面计算网格示意情况见图 6。

图  6  模型计算网格

Fig.  6  Computational meshes of the Xiaolangdi Reservoir

下载: 全尺寸图片 幻灯片

3.3   淤积形态验证分析

模型依据小浪底水库2010年汛前(2010年4月)实测地形建立, 水库床面淤积物粒径采用2010年沿程实测资料, 采用2010年4—10月实际入库水沙条件为模型进口边界, 以水库运行方式(出库水沙过程)作为模型出口边界, 在此基础上, 开展了小浪底水库的冲淤验证工作。挟沙力公式采用多沙河流常用的计算公式^[18], 泥沙计算时间步长60 s, 细颗粒淤积物与上层水体之间的摩阻系数f_s以及细颗粒淤积物与下部床面之间的摩阻系数f_m分别取0.02和5, 泥沙落淤后, 淤积物密度取1 080 kg/m³。

小浪底水库干流剖面形态为典型三角洲淤积形态, 冲淤分布统计结果见表 1, 实测深泓线纵向淤积形态与模拟结果对比见图 7, HH56断面至坝前淤积总量的实测与模拟结果分别为1.050亿m³和0.970亿m³, 两者相对偏差在10%以内, 数学模型计算结果可较好地反映小浪底水库的冲淤特征。

表  1  2010年4—10月干流冲淤量分布对比 亿m³

Table  1  Comparisons between the simulated and observed volume of sediment deposition

区间	实测	计算
HH37断面以上	-0.110	-0.120
HH37断面—三角洲顶点	-0.570	-0.430
三角洲顶点—坝前	1.730	1.520
HH56—坝前(总量)	1.050	0.970

下载: 导出CSV 

| 显示表格

图  7  水库纵向淤积形态

Fig.  7  Comparisons between the simulated and observed sedimentation patterns

下载: 全尺寸图片 幻灯片

三角洲形态及顶点位置随着水库的运行而发生变化。2010年汛前(4月), 三角洲顶点在位于距坝约24 km附近处(图 7), 高程约为219.6 m。2010年汛期, 坝前水位降低运行(图 8), 三角洲洲面发生较明显冲刷, 水库排沙比38.9%, 前坡段与坝前段则出现明显泥沙淤积；2010年汛后(10月), 三角洲顶点向下游推进到距坝约20 km附近处(图 7), 三角洲顶点高程约为215.6 m, 较汛前高程降低约4.0 m。

图  8  2010年汛期坝前水位运行情况

Fig.  8  Variations of water stages in the Xiaolangdi Reservoir during the flood season of 2010

下载: 全尺寸图片 幻灯片

另外, 还选取水库干流典型断面(各断面位置见图 5), 对实测与计算的断面冲淤形态进行验证分析, 具体如图 9所示。结合图 7和图 9可见, 水库冲淤主要发生在高程相对较低的主槽部分, 水库上游段的断面(距坝里程在22 km以上)略有冲刷, 主要发生在主槽部位, 而近坝段断面(距坝里程在20 km以内, HH13断面以下)则淤积相对明显, 数学模型计算的断面冲淤分布与实测结果定性基本一致, 具体断面在定量上还存在一定差异, 但从整个库区的冲淤分布统计来看(表 1), 总体计算结果与实测结果基本吻合。

图  9  典型断面冲淤分布验证

Fig.  9  Simulated and measured results of deposition and erosion distribution at typical cross-sections

下载: 全尺寸图片 幻灯片

此外, 从小浪底水库淤积物粒径沿程分布来看, 除局部存在一定起伏变化外, 总体而言, 还是受泥沙分选作用, 从库尾至坝前沿程泥沙颗粒逐渐变细(图 10)。相对而言, 坝前5 km范围内的淤积物粒径较细, 中值粒径基本在0.01 mm以下；距坝5~20 km范围内, 淤积物中值粒径一般在0.015 mm以下；距坝20~57 km范围内, 淤积物中值粒径略为变粗, 一般在0.02 mm以下；距坝57 km以上, 淤积物粒径明显变粗, 中值粒径由0.02 mm逐渐变粗为0.25 mm左右, 有的还在0.25 mm以上。由此可见, 水库近坝段淤积在主槽内的泥沙粒径较细, 此类淤积物在密实前具有较强的流动性, 平衡坡降极小, 本次模拟过程对粒径小于0.02 mm的细颗粒淤积物自身失稳流动特性进行了考虑, 因此, 近坝段模拟的断面淤积形态均较为平坦(图 9中断面HH2、HH6、HH13、HH17), 这一结果与实际观测一致, 而未考虑淤积物自身流动特性时则断面坡度相对较大；另外, 从断面深泓高程来看, 考虑细颗粒淤积物流动特性后深泓高程要高于未考虑其流动特性的结果, 与实际情况亦相对吻合；由此说明, 对于细颗粒泥沙淤积形态的模拟, 除了考虑传统水沙输移外, 还需考虑淤积物自身的流动特性。相对而言, 50 km以上断面(如断面HH37、HH44)淤积物相对较粗, 模拟过程中不考虑其流动特性, 因此断面坡度略大(冲刷后淤积物可保持较明显的边坡), 实际情况亦与此基本一致。但同时也应指出, 水库细颗粒泥沙运动非常复杂, 本研究仅考虑水库冲淤过程以及落淤后细颗粒淤积物自身的流动过程, 即通过水沙输移计算, 模拟水库冲淤总量的分布相似, 在此基础上对细颗粒淤积物自身流动特性加以模拟, 还未涉及前期淤积物在后期水沙运动(异重流)过程中的起动(运动)特性, 这对水库输沙以及淤积形态也会产生进一步影响, 还有待进一步深入研究。此外, 小浪底水库支流众多, 干支流交汇, 本文仅考虑干流水库冲淤过程, 并未对支流冲淤过程进行模拟, 同样也会对计算结果造成一定影响, 亦有待进一步精细分析。

图  10  小浪底水库淤积物中值粒径(d₅₀)沿程变化

Fig.  10  Variations of sedimentation diameter (d₅₀) in the Xiaolangdi Reservoir

下载: 全尺寸图片 幻灯片

4   结论

(1) 低密度水库细颗粒淤积物为宾汉型流体, 其屈服应力τ_B和黏度系数μ随淤泥密度的不同而发生变化, 中间存在从缓变到急变的转化过程, 当密度大于1.25 g/cm³后, 流变参数随密度的变化速率快速增大, 因此, 可认为当淤积物密度大于1.25 g/cm³后, 淤积物的流动性已快速减弱, 将不易流动。

(2) 采用临界坡度作为细颗粒淤积物失稳流动的判别标准, 即J>J_c；通过引入水、淤积物、床面之间的界面受力分析, 建立了细颗粒淤积物流动模式, 描述其失稳后的运动及其泥沙重新分配过程；通过与传统水沙模型相耦合, 建立了考虑细颗粒淤积物流动特性的水库淤积形态数值模拟方法。

(3) 从小浪底水库淤积物粒径沿程分布来看, 受泥沙分选作用, 从库尾至坝前泥沙粒径逐渐细化, 尤其是近坝段淤积在主槽内的泥沙较细(中值粒径在0.01 mm左右), 此类淤积物在密实前具有较强的流动性, 平衡坡降极小, 本次模拟过程对细颗粒淤积物自身的失稳流动特性进行了考虑, 因此, 近坝段模拟的断面淤积形态均较为平坦, 这一结果与实际观测一致；相对而言, 对于淤积物相对较粗的水库中上段, 模拟过程中不考虑其流动特性, 因此断面坡度略大, 实际情况亦与此基本一致。