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

技术特征：

1.一种多沙河流水库浑水明流与异重流耦合模拟方法，其特征在于，包括以下步骤：

步骤(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均不满足潜入判别条件，则令整个计算域内m＝m+1，返回步骤(2)进行下一时段的计算；

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

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

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

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

2.根据权利要求1所述的多沙河流水库浑水明流与异重流耦合模拟方法，其特征在于：

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

$\frac{\∂U}{\∂t}+\frac{\∂F}{\∂x}=S---\left(1\right)$

$U=\left[\begin{array}{c}A\\ Q\\ AC\end{array}\right],F=\left[\begin{array}{c}Q\\ Q^2/A\\ QC\end{array}\right],S=\left[\begin{array}{c}B(E-D)/(1-p)-q_l\\ M_w\\ B(E-D)-C_l{q}_l\end{array}\right]---\left(2\right)$

$\begin{array}{c}{M}_w=-\frac{{\mathrm{\ρ}}_b-\mathrm{\ρ}}{\mathrm{\ρ}\left(1-\mathrm{\ρ}\right)}BU\left(E-D\right)-\frac{{\mathrm{\ρ}}_s-{\mathrm{\ρ}}_w}{\mathrm{\ρ}}\left({\mathrm{Cq}}_l-C_l{q}_l\right)U+\\ gA(S_b-S_f)-gA\frac{\∂h}{\∂x}-{\mathrm{gh}}_c\frac{A}{\mathrm{\ρ}}\frac{\∂\mathrm{\ρ}}{\∂x}\end{array}---\left(3\right)$

$\frac{\∂A_b}{\∂t}=\frac{B(D-E)}{1-p}---\left(4\right)$

式中：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.根据权利要求2所述的多沙河流水库浑水明流与异重流耦合模拟方法，其特征在于：

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

$U_i^{m+1}=U_i^m-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(F_{i+1/2}^m-F_{i-1/2}^m)+{\mathrm{\ΔtS}}_i^m---\left(5\right)$

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

4.根据权利要求2所述的多沙河流水库浑水明流与异重流耦合模拟方法，其特征在于：

上述步骤(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_*演进到下一时间步：

$A_{*}^{m+1}=A_{*}^m-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(F_{1,r+1/2}^m-F_{1,r-1/2}^m)+{\mathrm{\ΔtS}}_{1,r}^m---\left(7\right)$

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

5.根据权利要求2所述的多沙河流水库浑水明流与异重流耦合模拟方法，其特征在于：

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

$q_{tl}=h_0^{\frac{3}{2}}\sqrt{\frac{(2J\frac{L}{h_0}-4)\mathrm{\ξ}-6J\frac{L}{h_0}+4}{{(3-\mathrm{\ξ})}^3}(1+\mathrm{\ξ})\mathrm{\η}g}---\left(8\right)$

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

6.根据权利要求1所述的多沙河流水库浑水明流与异重流耦合模拟方法，其特征在于：

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

${\mathrm{Fr}}_p^2=0.191{\left(S_v\right)}^{0.75}---\left(9\right)$

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

7.根据权利要求1至6中任一项所述的多沙河流水库浑水明流与异重流耦合模拟方法，其特征在于：

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

$\frac{\∂T}{\∂t}+\frac{\∂G}{\∂x}=R---\left(10\right)$

$T=\left[\begin{array}{c}{A}_t\\ Q_t\\ A_t{C}_t\end{array}\right],G=\left[\begin{array}{c}{Q}_t\\ Q_t^2/A_t\\ Q_t{C}_t\end{array}\right],R=\left[\begin{array}{c}{B}_t\left(E_t-D_t\right)/\left(1-p\right)+B_t{e}_w{U}_t-q_{tl}\\ M_t\\ B_t\left(E_t-D_t\right)-C_t{q}_{tl}\end{array}\right]---\left(11\right)$

$\begin{array}{c}{M}_t=-\frac{{\mathrm{\ρ}}_b-{\mathrm{\ρ}}_t}{{\mathrm{\ρ}}_t\left(1-p\right)}{B}_t{U}_t\left(E_t-D_t\right)+\frac{{\mathrm{\ρ}}_t-{\mathrm{\ρ}}_w}{{\mathrm{\ρ}}_t}{B}_t{e}_w{U}_t^2+g^{\′}{A}_t{S}_b-\\ {\mathrm{gA}}_t{S}_f^{\′}-{\mathrm{gA}}_t\frac{{\mathrm{\ρ}}_w}{{\mathrm{\ρ}}_t}\frac{\∂z_s}{\∂x}-g^{\′}{A}_t\frac{\∂h_t}{\∂x}-{\mathrm{gh}}_{ct}\frac{A_t}{{\mathrm{\ρ}}_t}\frac{\∂{\mathrm{\ρ}}_t}{\∂x}\end{array}---\left(12\right)$

$\frac{\∂A_b}{\∂t}=\frac{B_t(D_t-E_t)}{1-p}---\left(13\right)$

式中：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)中得到的计算。

8.根据权利要求7所述的多沙河流水库浑水明流与异重流耦合模拟方法，其特征在于：

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

$T_i^{m+1}=T_i^m-\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(G_{i+1/2}^m-G_{i-1/2}^m)+{\mathrm{\ΔtR}}_i^m---\left(14\right)$

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

9.根据权利要求7所述的多沙河流水库浑水明流与异重流耦合模拟方法，其特征在于：

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

$\frac{{\mathrm{\ΔA}}_b}{\mathrm{\Δ}t}={\[\frac{B_t(D_t-E_t)}{1-p}\]}^m---\left(15\right)$

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