冲沙河道演变过程的数值研究

申 俊

(塔里木河干流管理局，新疆 库尔勒 841000)

本文采用三维数值模型，模拟了新疆地区塔里木河流域希尼尔水库2020年冲洗作业过程中河床随表面速度场的演变趋势，并将模拟结果与现场测量结果进行比较。此外，还定量评估了冲出泥沙量在不同冲洗作业阶段的变化情况。

1 材料与方法

1.1 案例说明

本文所研究的大坝水库总库容和有效库容分别为9.01×106m3和1.66×106m3。图1显示了2020年6月进行的水库冲洗作业之前测量的河床水位，以及用于进一步定量评估河床演变模式的横截面位置(A-F)。将近2 km长的研究案例分为三个区域，即区域Ⅰ、Ⅱ和Ⅲ，根据每个区域存在的河床物质类型分析研究。Ⅰ区河床物质类型为粗粒型，Ⅲ区和Ⅱ区后半部分河床物质类型则变为细粒型。研究区侵蚀沉积物的总体积为408.7×103m3。在该水库中，中间区域的左岸(即区域Ⅱ)沉积的细泥沙(图1中矩形框所示)无法有效清除。表1显示了冲洗作业前不同横截面的平均泥沙粒径分布。选取0.37 ～ 316 mm的7种沉积物尺寸作为代表性粒径。图2显示了2020年冲洗作业期间的进入流量和水位变化，其中，主要的泥沙流入量是冲洗负荷。开始作业后，水位初步下降时间为8～24 h，自由流动状态为24～38 h。

图1 水库冲水作业前的实测河床地形以及用于评估河床演变的横截面A-F位置

图2 冲洗作业期间的水位和进水流量变化

1.2 数值模型

数值模型求解质量守恒方程和雷诺平均N-S方程来计算湍流的水运动，如式(1)、式(2)：

(1)

(2)

本文采用有限体积法进行离散化[1-5]，采用二阶迎风格式求解雷诺平均N-S方程中的对流项。湍流由标准k-ε模型建模，网格是非正交的、非结构化的和自适应的，并随着河床和自由水面的变化而垂直移动。因此，仅对水体进行建模。网格在每次时间步长后重新生成，采用的润湿/干燥算法使模型能够根据水和床位的变化具有不同数量的网格单元(例如在垂直方向上)。润湿/干燥算法消除了水深小于定义边界的单元。

进水采用第一类边界条件(对数速度分布)，出水和泥沙采用零梯度边界条件。对于河床和墙壁的边界条件，在没有水通量的情况下，采用壁面定律。河床形式，即沙丘效应和波纹效应，也通过经验公式考虑。模拟形态变化的泥沙输移计算分为悬移质输移和推移质输移。悬移质通过求解瞬态对流-扩散方程公式计算得出，推移质可通过修正的MPM公式或Van Rijn公式来模拟。MPM公式更适用于陡峭的河流，这些河流主要输送靠近河床的粗泥沙，如式(3)所示。

(3)

式中：qb，i为单位宽度推移质第i部分的泥沙输送速率，kg/s；g为重力加速度，m/s2；ρs为泥沙密度，kg/m3；ρw为水的密度，kg/m3；d50为特征泥沙粒径，mm；r为水力半径，mm；I为能量线的斜率。

Van Rijn公式如式(4)为：

(4)

式中：di为第i部分的直径，m；τ为剪切应力，N/m2；τc，i为di的临界剪切应力，N/m2；v为运动黏度，Pa·s。模型中引入了在单个时间步长内可被侵蚀的上部活性层厚度以及为活性层提供足够沉积物的下部非活性层厚度。

该模型利用初始运动临界剪应力的折减函数和沙滑算法来考虑边坡效应。当河床坡度超过河床物质的休止角时，沙滑算法将修正河床坡度。如果没有内聚力影响，泥沙颗粒可以在边坡上相互独立移动，并可以适当模拟河岸侵蚀。

2 结果与讨论

2.1 模型校准

根据图1所示冲洗前测量的河床水位，制作计算网格。水流方向和横向的网格单元尺寸分别为10～20 m和5～10 m。假定河床物质为2650 kg/m3。图2显示了模拟中用作水动力边界条件的水位和进水流量波动。此外，还利用7种具有代表性的沉积物尺寸将水库内的非均匀河床物质尺寸分布引入模型。根据可用河床物质样品的尺寸插值获得库区内的河床物质尺寸分布。由于冲洗负荷是假定在水库中无沉积的情况下被输送，因此，在计算中忽略了冲洗负荷对河床演变的影响。

为了揭示以经验为主参数的影响趋势(如活性层厚度、河床物质的含水量以及沙滑算法的临界休止角)，首先，建立了一个参考案例，假设上述参数的一般值。 然后，改变假设值，计算冲出泥沙的总体积(以下称为TVFS)并与测量值(即4008.7×103m3)进行比较，以校准模型。表2显示了与参考案例相比，TVFS对选定的主要经验参数的敏感性分析。

表2 参考案例对所选经验参数的敏感性分析 m3

可以看出，TVFS随着活性层厚度和含水量的增加而增强，随着临界休止角的增加而降低。当大量的河床物质(即更厚的活性层)可以在一个时间步长内被侵蚀时，预计泥沙的侵蚀量会更大。假设泥沙沉积中较高的含水量(如50%)会降低河床物质的淹没密度，从而导致较高的夹带泥沙。在每个时间步长后，较高的临界休止角会使冲沙河道的侧岸更陡，导致水道进一步加深，这并不是增加TVFS的有效方法。相反，较低的临界休止角可以实现冲沙河道的横向发展，有利于增加TVFS。除了TVFS，最终将模拟的河床地形模式与测量结果进行了定性比较，以提供Ⅲ区下部周围侵蚀的信息。之后，对参考案例进行了修改并使用更新后的值进行验证(例如活性层厚度、含水量和临界休止角分别设置为1 m、40%和33°)。

2.2 冲沙河道演变模拟

2020年6月冲洗作业后数值模拟的河床水位与冲洗后的实测水位对比结果：在水位下降过程中，数值模型在自由流动条件下观察到独特的冲沙河道，其特征与原型相似。

图3显示了采用MPM公式、Van Rijn公式(分别为图3(a)和3(b))和实测结果(图3(c))的冲洗后的模拟河床水位。可以看出，应用Van Rijn公式可以降低Ⅲ区的侵蚀，而Ⅱ区的侵蚀被高估了。使用MPM公式计算的TVFS为313.14×103m3，而使用Van Rijn公式计算的TVFS为333.91 ×103m3。

图3 模拟河床水位

当采用MPM公式时，Ⅲ区被侵蚀的粗粒型物质可能会再次沉积到Ⅰ区。Ⅰ区水位较深导致流速降低。此外，在Ⅰ区，无论所用推移质输沙公式类型如何，模拟结果都显示出与观测结果相比更窄的冲沙河道。这是由于在靠近大坝的大弯曲区域是由二次流和反向流形成的复杂流场。该区域形成的河床形式(即沙丘和波纹)将增加水流的复杂性，进而增加泥沙输移模式。此外，该区域的输沙能力达到了最高水平，但仍低于原型中的原始值。因此，河床水位与实测水位不同。

由于水库不同河段的推移质类型和水位不同，每个推移质输运公式(即MPM和 Van Rijn)都可以对特定段的河床演变模式提供更合理的预测。因此，将MPM和Van Rijn公式模拟的最终河床高程与图4中横截面A-A和E-E的测量值进行了比较。选择MPM公式对数值结果进行进一步的定量评估。图5还给出了使用MPM公式在横截面A-A、B-B、D-D和F-F中数值模拟的河床演变模式。

图4 最终河床高程

图5 各个横截面的河床演变模式

在河床物质较粗的Ⅲ区，泥沙侵蚀被低估，模拟的冲沙河道宽度和深度小于测量值(图5(a)和图5(b))。由于粗糙河床物质的复杂流速剖面和有效水动力作用，特别是在低水头的自由流动状态下，模型无法很好地模拟该区域的侵蚀模式。另一个原因是在Ⅲ区下部左岸的计算域内消除了具有浅水头(即干涸)的单元，其中已经测量到明显的河床退化(图5(b))。在Ⅱ区，即水库的宽中段，通过应用MPM公式得到的冲沙河道形状和位置与观测结果基本一致。然而，在数值结果中，沿右路堤的河床退化被高估了(图5(c))。Ⅱ区的冲沙河道特征对保持水库蓄水能力和提高冲洗效率有很大贡献。在自由流动冲洗作业期间，狭窄的冲沙河道位置偏向Ⅱ区的右岸区域，因此，左岸区域中大量沉积物的侵蚀可能性较低。还应考虑到研究水库是在雨季进行冲沙，并且在冲洗前后不进行河床水位测量。因此，测量的河床水位可能存在不确定性。此外，沿非活性层也没有高空间分辨率的泥沙粒径分布测量，但这会显著影响数值模型结果的质量。

图6显示了在自由流动条件下具有二次流矢量的水流速度。在冲沙河道中，流速升至4.5 m/s，多个区域出现超临界流(图6(a)和图6(b))。三维模型还捕捉了河道弯道中河床发育的特征(即弯道外侧的侵蚀和内侧的沉积)，并再现了宽度上的非对称速度剖面，以及弯道顶点的侧向水面倾斜(图6(c))。

图6 自由流动条件下具有二次流矢量的水流速度

2.3 冲沙量变化

根据冲洗过程中每隔一小时冲洗掉的侵蚀河床物质的计算通量可知，平均流速和湍流波动的增加，即初流阶段和自由流动阶段的最后部分，往往会增强泥沙的流动性和随后的泥沙夹带量。 较小的沉积物尺寸往往比较大的沉积物尺寸更早被侵蚀和冲走，如直径为3.7 mm的沉积物在t=5 h的洪峰后立即出现高浓度。直径较小的沉积物对水位和流量波动都很敏感。另一方面，直径较大的沉积物对水位变化比对流量变化更敏感。

3 结 论

(1)三维数值模型能够较好地模拟真实边界条件(如非稳态水动力条件)下的冲沙河道演变模式。Van Rijn和MPM推移质输沙公式都可以在一定范围内代表冲洗过程中冲出的泥沙量。

(2)MPM公式能较好地预测由较粗物质覆盖的上游区域的河床水位变化。Van Rijn公式在靠近大坝的区域表现出更好的性能，其中，河床主要被细粒物质覆盖，而且水库也更深。

(3)三维数值模型可以识别弯道中的螺旋流以及由此产生的侵蚀和沉积模式，比二维建模更具优势。

(4)就冲洗过程中冲出的泥沙量来看，较粗的河床物质主要在初步水位下降结束时和自由流动状态下冲出。因此，在自由流动条件下引入额外的流量和小幅波动的水位降低可能会增加泥沙侵蚀的可能性。