多沙河流水库浑水明流与异重流耦合模拟方法与流程

本发明涉及水利工程技术领域，具体的说是多沙河流水库浑水明流与异重流耦合模拟方法。

背景技术：

水库异重流是多沙河流水库泥沙管理中引起普遍关注的问题。在水库设计或运用阶段对异重流问题考虑不足有可能影响其正常运转甚至减损水库寿命。目前对于水库异重流在动力学机制上仍有认识不足，模型发展现状还不能满足实际工程需要。此外在水库调水调沙过程中，浑水明流段的水沙运动对于异重流的形成和运动有重要影响，在此种情况下对整个库区的水沙运动进行数值模拟存在很大的挑战。

现有异重流模型通常分为解析异重流模型、一维和平面二维异重流模型，以及立面二维和三维异重流模型。

解析异重流模型只能对于无交界面掺混和无床面泥沙交换的理想情况给出近似解。目前解析模型的研究成果只能用于规则断面形态下，具有特定上游边界条件下(即总体积固定和恒定入流两种情况)的异重流传播计算。立面二维和三维异重流模型可以模拟出流速、含沙量沿垂向的不均匀分布。这类模型构建时大部分涉及到紊流封闭，也有一部分使用直接数值模拟。由于目前对于紊动现象的理解不足，不同的学者提出了许多紊流模式理论，相应的模型对于同一异重流过程给出的预测结果差异很大。此外立面二维和三维异重流模型计算成本很高，大部分应用仍局限于对水槽实验过程的模拟。

一维和平面二维异重流模型是在对Navier-Stokes方程深度平均的基础上得来的，虽然不能求解水沙要素的垂向分布，但是计算速度快且能够确定工程应用中关心的传播速度、冲淤量等问题。传统的一维和平面二维异重流模型由于不考虑明流段与异重流段水沙运动的耦合关系，只能用于模拟潜入点下游的异重流传播过程而很难模拟异重流的产生。另一方便对于库区存在支流的多沙河流水库，现有的模型没有考虑干支流倒回灌的影响。

技术实现要素：

本发明的目的在于提供一种多沙河流水库浑水明流与异重流耦合模拟方法，能预测水库异重流的形成时间和潜入位置，以及异重流发生期间明流段和异重流段的水沙输移过程，干支流倒回灌过程，水库排沙比等，为可持续的水库泥沙管理提供科学依据和技术支撑。

本发明的原理为：多沙河流水库中的河床冲淤演变受到浑水明流段和异重流段水沙运动的共同作用，同时受到干支流倒灌的影响，有必要在模拟时全面考虑这些因素的影响。入库泥沙在明流段的输移或明流段的冲刷为异重流的形成提供了泥沙，洪水在明流段的演进影响了异重流中期潜入位置的变动，支流口门异重流的厚度和含沙量决定了异重流倒灌的强度，而支流向干流的回灌又影响了明流段和异重流段的水位。因此，本发明结合控制方程交替求解模式、异重流潜入判别条件、零维水库法和异重流倒灌流量公式，提出多沙河流水库浑水明流与异重流耦合模拟方法。

本发明多沙河流水库浑水明流与异重流耦合模拟方法的数值计算方法，包括如下步骤：

步骤(1)，根据水库干支流实测断面地形划分干流断面滩槽节点，计算支流底坡J和水位库容关系V(z_s)；将模拟区域划分为一组控制体；

步骤(2)，求解浑水明流控制方程组，使用零维水库法和考虑底坡的异重流倒灌流量公式计算干支流倒回灌流量，具体步骤如下：

步骤(2.1)，根据下一时刻的边界条件以及当前时刻的地形和守恒变量U^m，计算各控制体交界面上的数值通量F；

步骤(2.2)，对于和支流相连的控制单元，使用零维水库法计算下一时刻水位及干支流倒回灌净单宽流量q_l，如果该控制单元位置还存在异重流，则用考虑底坡的异重流倒灌流量公式计算倒灌流量q_tl；对于不和支流相连的控制体，直接按显式离散格式计算得到

步骤(2.3)，完成浑水明流控制方程中第二分量流量、第三分量含沙量的计算得到U^m+1，更新后的河底高程z_b暂记为

步骤(3)，使用潜入判别条件判断异重流潜入位置；

从整个模拟区域的上游第一个断面位置开始，逐个计算每一控制体的U^m+1是否满足潜入判别条件，若满足则记录下断面号p，将明流段河底高程更新为进入步骤(4)；若所有断面位置的U^m+1均不满足潜入判别条件，则令整个计算域内返回步骤(2)进行下一时段的计算。

步骤(4)，求解异重流控制方程组；

步骤(4.1)，将作为异重流段上游边界条件，计算异重流控制方程中的数值通量和源项；

步骤(4.2)，计算得到下一时刻的和N为模拟区域内的断面个数，如果已经达到模拟的最终时刻，则进入步骤(5)，否则令m＝m+1，返回步骤(2)进行下一时段的计算；

步骤(5)，计算干支流冲淤量及水库排沙比。

上述步骤(2)中，浑水明流控制方程组如下：

式中：t表示时间；x表示沿程距离；U是守恒变量；F是数值通量；S是源项；A是断面过流面积；Q是断面流量；C为体积比含沙量；B为水面宽度；E为床沙上扬通量；D为悬沙沉降通量；U是断面平均流速；ρ_w，ρ分别为清混水的密度；ρ_b为床沙饱和湿密度；ρ_s为泥沙颗粒密度；p为床沙的孔隙率；g为重力加速度；S_b为底坡；S_f为阻力坡度，其表达式为其中n_j为步骤(1)中指定的子断面糙率，A_j和R_j为子断面的面积和水力半径，NP表示某一断面上划分的子断面个数；h为水深；h_c为过流断面形心高度；q_l为干支流倒回灌净单宽流量，规定从干流到支流为正，C_l为与这部分流量对应的体积比含沙量；A_b表示基准高程以上的河床横断面面积；为冲淤导致的断面面积的变化速率。

上述步骤(2)中，采用显式离散格式求解控制方程，即将式(1)写为：

其中i是控制体编号，m是时间步长数，Δt为时间步长，Δx为空间步长，和分别表示第i个控制体左右两侧的数值通量。

上述步骤(2)中，零维水库法如下：

在更新浑水明流方程中U的第一个分量即A时，对于与支流相连的控制体，用下标为r的控制体表示，将其与支流作为一个整体，记A_*为这一扩大的控制体的断面面积；A_*的定义为：

A_*(z_s,r)＝A_r(z_s,r)+V(z_s,r)/Δx (6)

式中：z_s,r为下标为r的控制体中的水位，V(z_s,r)表示水位为z_s,r时与下标为r的控制体相连的支流的库容，A_r(z_s,r)表示水位为z_s,r时下标为r的控制体的断面面积，A_*(z_s,r)表示水位为z_s,r时扩大的控制体的断面面积，Δx为空间步长；对于上述扩大的控制体，在计算其与两侧普通的控制体之间的数值通量时，不考虑支流的存在，然后按下式将A_*演进到下一时间步：

其中Δt为时间步长，Δx为空间步长，表示当前时刻下标为r的控制体右侧的数值通量的第一个分量，表示当前时刻下标为r的控制体左侧的数值通量的第一个分量，表示当前时刻下标为r的控制体的源项的第一个分量，表示当前时刻扩大的控制体的断面面积，表示下一时刻扩大的控制体的断面面积；这里的中不包括干支流倒回灌净单宽流量q_l，因为倒回灌的水流是扩大的控制体内部的流动；通过式(7)得到之后，干支流交汇处的水位即下标为r的控制体中的水位z_s,r就可以由式(6)得到，那么原来干流中第r个控制体内的断面面积即为

上述步骤(2)中，考虑底坡的异重流倒灌流量公式为：

其中q_tl为倒灌流量；h₀为交汇区干流异重流厚度；ξ为阻力损失系数；J为支流底坡；L为潜入段长度；η＝(γ-γ₀)/γ，γ₀和γ分别为清浑水容重；g为重力加速度。

上述步骤(3)中，潜入判别条件为：

式中：Fr_p为潜入点弗汝德数，S_v为浑水水流的体积比含沙量；如果由某控制体的U^m+1计算的弗汝德数Fr小于潜入点弗汝德数Fr_p，即认为满足潜入判别条件。

上述步骤(4)中，异重流控制方程组如下：

式中：t表示时间；x表示沿程距离；T是守恒变量；G是数值通量；R是源项；A_t，Q_t，C_t是异重流层的面积，流量和含沙量；B_t是清浑水层交界面的宽度；ρ_w为清水的密度；ρ_t，U_t是异重流层的密度与断面平均流速；h_t是异重流层厚度；E_t和D_t是异重流与床面的泥沙交换通量；e_w是清水掺入系数；p为床沙的孔隙率；g为重力加速度；g′＝(ρ_t-ρ_w)g/ρ_t为有效重力加速度；S_b为底坡；S′_f为考虑交界面阻力后的综合阻力坡度；h_ct为异重流断面形心高度；A_b表示基准高程以上的河床横断面面积；为冲淤导致的断面面积的变化速率；q_tl表示考虑底坡的异重流倒灌流量；

计算源项R时，其第二个分量M_t中的水面梯度用步骤(2.2)中得到的计算。

上述步骤(4)中，由显式离散格式计算得到：

其中Δt为时间步长，Δx为空间步长，和表示当前时刻第i个控制体右侧和左侧的数值通量，表示当前时刻第i个控制体的源项，和表示当前和下一时刻第i个控制体的守恒变量。

上述步骤(4)中，的计算方法如下，对式(12)直接按照时间前差进行离散：

将上式得到的冲淤量ΔA_b分配到异重流段的断面结点上，就得到了

本发明多沙河流水库浑水明流与异重流耦合模拟方法的有益效果在于：

(1)采用了水沙耦合的浑水明流与异重流控制方程，能够反映高含沙量或急剧床面冲淤对洪水演进的影响；

(2)可预测异重流的潜入位置和形成时间、异重流传播过程，以及异重流到达坝前后浑水层抬高和排出的过程；

(3)可根据明流段和异重流段的模拟结果，计算泥沙冲淤总量，以及冲淤量在干流不同库段和干支流间的分布，可以计算水库排沙比。

附图说明

图1为本发明实施例计算方法的流程图；

图2为本发明实施例横断面滩槽节点划分示意图；

图3为本发明实施例浑水明流段的横断面示意图；

图4为本发明实施例零维水库法示意图；

图5为本发明实施例异重流倒灌的示意图；

图6为本发明实施例异重流段的横断面示意图；

图7为本发明实施例模拟的不同时刻库区水位和交界面高程分布,(a)t＝317.57h，(b)t＝378h。

具体实施方式

下面结合附图和实施例，对本发明做进一步说明。如图1-图7所示，多沙河流水库浑水明流与异重流耦合模拟方法，具体步骤如下：

步骤1：根据水库干支流实测断面地形划分干流断面滩槽节点，计算支流底坡J和水位库容关系V(z_s)；将模拟区域划分为一组控制体；

由各个实测断面的起点距高程数据绘制断面形态图，如图2所示，图中横向底坡的转折点定义为滩槽分界点，两个滩槽分界点之间的节点定义为主槽节点，滩槽分界点至断面左端点或右端点之间的节点定义为滩地节点，每两个断面节点之间为一个子断面。如果某子断面两端均为主槽节点，则令该子断面的糙率n_j等于主槽糙率n_MC；如果某子断面两端含有滩地节点，则令该子断面的糙率n_j等于滩地糙率n_FP，一般情况下n_FP略大于n_MC。

对于库区内每条支流，根据实测断面数据确定支流底坡J和深泓线最低点，在最低点高程与水库校核洪水位之间划分若干水位，计算各水位下相应的库容，得到各个支流的水位库容关系V(z_s)，z_s表示水位。

模拟区域即为需要研究分析的库区范围，其下游边界到达水库坝址，上游边界到达水库回水范围以上可获得入库水沙测量数据的测站位置。控制体即为数值模拟时对空间上连续的地形以及水位，流量等物理量进行离散化表示的基本单元。划分控制体时，每个实测断面对应一个控制体，每两个相邻实测断面的中点位置定义为与这两个断面对应的控制体的交界面。

步骤2：求解浑水明流控制方程组，使用零维水库法和考虑底坡的异重流倒灌流量公式计算干支流倒回灌流量，包括步骤2.1-2.3。(步骤2需要求解的是方程(1)中的A,Q,C和方程(4)中的A_b，步骤2.1计算F，步骤2.2完成式(1)第一行的求解，步骤2.3完成式(1)第二三行和式(4)的求解)

步骤2.1：根据下一时刻的边界条件以及当前时刻的地形和守恒变量U^m，计算各控制体交界面上的数值通量F。

边界条件包括上游边界条件和下游边界条件，上游边界条件指入库流量和含沙量，下游边界条件指水库坝前水位。地形指每个控制体所对应的实测断面，地形用来计算底坡Sb，以及由过流面积A计算水深h。

本发明采用考虑干支流倒回灌的水库浑水明流与异重流耦合模型，模型包括浑水明流与异重流两组控制方程，浑水明流控制方程如下：

式中：t表示时间；x表示沿程距离；U是守恒变量；F是数值通量；S是源项；A是断面过流面积；Q是断面流量；C为体积比含沙量；B为水面宽度；E为床沙上扬通量；D为悬沙沉降通量；U是断面平均流速；ρ_w，ρ分别为清混水的密度；ρ_b为床沙饱和湿密度；ρ_s为泥沙颗粒密度；p为床沙的孔隙率；g为重力加速度；S_b为底坡；S_f为阻力坡度，其表达式为其中n_j为步骤1中指定的子断面糙率，A_j和R_j为子断面的面积和水力半径，NP表示某一断面上划分的子断面个数；h为水深；h_c为过流断面形心高度；q_l为干支流倒回灌净单宽流量，规定从干流到支流为正，C_l为与这部分流量对应的体积比含沙量；A_b表示基准高程以上的河床横断面面积；为冲淤导致的断面面积的变化速率。部分变量的定义如图3所示。

本模拟方法采用显式离散格式求解控制方程，即将式(1)写为：

其中i是控制体编号，m是时间步长数，Δt为时间步长，Δx为空间步长，和分别表示第i个控制体左右两侧的数值通量。

步骤2.2：对于和支流相连的控制体，使用零维水库法计算下一时刻水位及干支流倒回灌净单宽流量q_l，如果该控制体位置还存在异重流，则用考虑底坡的异重流倒灌流量公式计算倒灌流量q_tl；对于不和支流相连的控制体，直接按显式离散格式计算得到

零维水库法的具体实施方法如下。

在更新浑水明流方程中U的第一个分量即A时，对于与支流相连的控制体(图4中下标为r的控制体)，将其与支流作为一个整体，记A_*为这一扩大的控制体的断面面积。A_*的定义为：

A_*(z_s,r)＝A_r(z_s,r)+V(z_s,r)/Δx (6)

式中：z_s,r为下标为r的控制体中的水位，V(z_s,r)表示水位为z_s,r时与下标为r的控制体相连的支流的库容，A_r(z_s,r)表示水位为z_s,r时下标为r的控制体的断面面积，A_*(z_s,r)表示水位为z_s,r时扩大的控制体的断面面积，Δx为空间步长。对于图4中粗实线内的控制体，在计算其与两侧普通的控制体之间的数值通量时，不考虑支流的存在，然后按下式将A_*演进到下一时间步：

其中Δt为时间步长，Δx为空间步长，表示当前时刻下标为r的控制体右侧的数值通量的第一个分量，表示当前时刻下标为r的控制体左侧的数值通量的第一个分量，表示当前时刻下标为r的控制体的源项的第一个分量，表示当前时刻扩大的控制体的断面面积，表示下一时刻扩大的控制体的断面面积。这里的中不包括干支流倒回灌净单宽流量q_l，因为倒回灌的水流是扩大的控制体内部的流动。通过式(7)得到之后，干支流交汇处的水位(即下标为r的控制体中的水位z_s,r)就可以由式(6)得到，那么原来干流中第r个控制体内的断面面积即为

考虑底坡的异重流倒灌流量公式为：

其中h₀为交汇区干流异重流厚度；ξ为阻力损失系数；J为步骤1中得到的支流底坡；L为潜入段长度；η＝(γ-γ₀)/γ，γ₀和γ分别为清浑水容重。部分变量定义如图5所示。

步骤2.3：完成浑水明流控制方程组中式(1)第二三行计算得到U^m+1，完成浑水明流控制方程组中式(4)的计算得到河床冲淤量，将冲淤量分配到断面各个节点从而修改(更新)河底高程z_b，修改后的高程暂记为

当床面发生淤积时，将淤积量平均分配到断面各个节点；当床面冲刷时，将冲刷量只分配到主槽节点上。

步骤3：使用潜入判别条件判断异重流潜入位置。

从整个计算域的上游第一个断面位置开始，逐个计算每一控制体的U^m+1是否满足潜入判别条件，若满足则记录下断面号p，将明流段河底高程更新为进入步骤4；若所有断面位置的U^m+1均不满足潜入判别条件，则令整个计算域内返回步骤2进行下一时段的计算。

这里用到的潜入判别条件为：

式中：Fr_p为潜入点弗汝德数，S_v为浑水水流的体积比含沙量。如果由某控制体的U^m+1计算的弗汝德数小于潜入点弗汝德数Fr_p，即认为满足潜入判别条件。某控制体的弗汝德数表达式为v＝Q/A,Q和A都包含在U^m+1，h可以根据A和断面形态计算出来。

步骤4：求解异重流控制方程组。包括步骤4.1-4.2。(步骤4用于求解方程(10)中的At,Qt,Ct和方程(13)中的Ab。)

步骤4.1：将作为异重流段上游边界条件，计算异重流控制方程中的数值通量和源项。

异重流控制方程如下：

式中：T是守恒变量；G是数值通量；R是源项；A_t，Q_t，C_t是异重流层的面积，流量和含沙量；B_t是清浑水层交界面的宽度；ρ_t，U_t是异重流层的密度与断面平均流速；h_t是异重流层厚度；E_t和D_t是异重流与床面的泥沙交换通量；e_w是清水掺入系数；g′＝(ρ_t-ρ_w)g/ρ_t为有效重力加速度；S′_f为考虑交界面阻力后的综合阻力坡度；h_ct为异重流断面形心高度。部分变量定义如图6所示。

计算源项R时，其第二个分量M_t中的水面梯度用步骤2.2中得到的计算。

步骤4.2：计算得到下一时刻的和如果已经达到模拟的最终时刻，则进入步骤5，否则令m＝m+1，返回步骤2进行下一时段的计算。

由显式离散格式计算得到，即：

其中Δt为时间步长，Δx为空间步长，和表示当前时刻第i个控制体右侧和左侧的数值通量，表示当前时刻第i个控制体的源项，和表示当前和下一时刻第i个控制体的守恒变量。

对式(12)直接按照时间前差进行离散：

将上式得到的冲淤量ΔA_b按照步骤2.3中的方法分配到异重流段的断面结点上，就得到了这里N为模拟区域内的断面个数。

步骤5：计算干支流冲淤量及水库排沙比，包括步骤5.1-5.2。

步骤5.1：计算库区冲淤量和水库排沙比。

模拟时段内入库泥沙量W₁由实测入库水沙资料计算。出库泥沙量W_N根据模型输出的结果计算，即：

式中：Δt为时间步长，M为模拟结束时的时间步数，m是时间步长数，取值1,2,3,…M，分别是干流出口断面的流量和含沙量。

库区冲淤总量为W₁-W_N，水库排沙比为W_N/W₁。

步骤5.2：计算支流淤积量。

对于水库内任意一条支流，按下式计算支流淤积量W_tl：

式中：Δt为时间步长，M为模拟结束时的时间步数，m是时间步长数，取值1,2,3,…M，b是支流口门宽度，即为式(7)计算的异重流倒灌流量，为倒灌异重流的含沙量，一般可取交汇处的干流异重流含沙量。

图7表示的是使用本发明方法所模拟的小浪底水库2004年调水调沙过程中产生的异重流。图7(a)显示的是异重流形成4小时后整个库区的纵剖面，模型预测的异重流潜入点距坝约46km，与当年实测结果一致。潜入点下游的清浑水交界面将异重流段的水体分为上下两层。在图7(a)所示时刻，入库流量Q_in为1014m³/s，坝前水位为233.79m。在图7(b)所示时刻，异重流已经到达坝前，坝前水位下降到229.08m，但是由于入库流量Q_in增加到4434m³/s，浑水明流段水深有所增加，同时异重流潜入点向下游移动了约2km。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。