技術(shù)特征：

1.一種多沙河流水庫(kù)渾水明流與異重流耦合模擬方法，其特征在于，包括以下步驟：

步驟(1)，根據(jù)水庫(kù)干支流實(shí)測(cè)斷面地形劃分干流斷面灘槽節(jié)點(diǎn)，計(jì)算支流底坡J和水位庫(kù)容關(guān)系V(z_s)；將模擬區(qū)域劃分為一組控制體；

步驟(2)，求解渾水明流控制方程組，使用零維水庫(kù)法和考慮底坡的異重流倒灌流量公式計(jì)算干支流倒回灌流量，具體步驟如下：

步驟(2.1)，根據(jù)下一時(shí)刻的邊界條件以及當(dāng)前時(shí)刻的地形和守恒變量U^m，計(jì)算各控制體交界面上的數(shù)值通量F；

步驟(2.2)，對(duì)于和支流相連的控制單元，使用零維水庫(kù)法計(jì)算下一時(shí)刻水位及干支流倒回灌凈單寬流量q_l，如果該控制單元位置還存在異重流，則用考慮底坡的異重流倒灌流量公式計(jì)算倒灌流量q_tl；對(duì)于不和支流相連的控制體，直接按顯式離散格式計(jì)算得到

步驟(2.3)，完成渾水明流控制方程中第二分量流量、第三分量含沙量的計(jì)算得到U^m+1，更新后的河底高程z_b暫記為

步驟(3)，使用潛入判別條件判斷異重流潛入位置；

從整個(gè)模擬區(qū)域的上游第一個(gè)斷面位置開始，逐個(gè)計(jì)算每一控制體的U^m+1是否滿足潛入判別條件，若滿足則記錄下斷面號(hào)p，將明流段河底高程更新為進(jìn)入步驟(4)；若所有斷面位置的U^m+1均不滿足潛入判別條件，則令整個(gè)計(jì)算域內(nèi)m＝m+1，返回步驟(2)進(jìn)行下一時(shí)段的計(jì)算；

步驟(4)，求解異重流控制方程組；

步驟(4.1)，將作為異重流段上游邊界條件，計(jì)算異重流控制方程中的數(shù)值通量和源項(xiàng)；

步驟(4.2)，計(jì)算得到下一時(shí)刻的和i＝p+1,p+2,...,N，N為模擬區(qū)域內(nèi)的斷面?zhèn)€數(shù)，如果已經(jīng)達(dá)到模擬的最終時(shí)刻，則進(jìn)入步驟(5)，否則令m＝m+1，返回步驟(2)進(jìn)行下一時(shí)段的計(jì)算；

步驟(5)，計(jì)算干支流沖淤量及水庫(kù)排沙比。

2.根據(jù)權(quán)利要求1所述的多沙河流水庫(kù)渾水明流與異重流耦合模擬方法，其特征在于：

上述步驟(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表示時(shí)間；x表示沿程距離；U是守恒變量；F是數(shù)值通量；S是源項(xiàng)；A是斷面過(guò)流面積；Q是斷面流量；C為體積比含沙量；B為水面寬度；E為床沙上揚(yáng)通量；D為懸沙沉降通量；U是斷面平均流速；ρ_w，ρ分別為清混水的密度；ρ_b為床沙飽和濕密度；ρ_s為泥沙顆粒密度；p為床沙的孔隙率；g為重力加速度；S_b為底坡；S_f為阻力坡度，其表達(dá)式為其中n_j為步驟(1)中指定的子斷面糙率，A_j和R_j為子斷面的面積和水力半徑，NP表示某一斷面上劃分的子斷面?zhèn)€數(shù)；h為水深；h_c為過(guò)流斷面形心高度；q_l為干支流倒回灌凈單寬流量，規(guī)定從干流到支流為正，C_l為與這部分流量對(duì)應(yīng)的體積比含沙量；A_b表示基準(zhǔn)高程以上的河床橫斷面面積；為沖淤導(dǎo)致的斷面面積的變化速率。

3.根據(jù)權(quán)利要求2所述的多沙河流水庫(kù)渾水明流與異重流耦合模擬方法，其特征在于：

上述步驟(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是控制體編號(hào)，m是時(shí)間步長(zhǎng)數(shù)，Δt為時(shí)間步長(zhǎng)，Δx為空間步長(zhǎng)，和分別表示第i個(gè)控制體左右兩側(cè)的數(shù)值通量。

4.根據(jù)權(quán)利要求2所述的多沙河流水庫(kù)渾水明流與異重流耦合模擬方法，其特征在于：

上述步驟(2)中，零維水庫(kù)法如下：

在更新渾水明流方程中U的第一個(gè)分量即A時(shí)，對(duì)于與支流相連的控制體，用下標(biāo)為r的控制體表示，將其與支流作為一個(gè)整體，記A_*為這一擴(kuò)大的控制體的斷面面積；A_*的定義為：

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

式中：z_s,r為下標(biāo)為r的控制體中的水位，V(z_s,r)表示水位為z_s,r時(shí)與下標(biāo)為r的控制體相連的支流的庫(kù)容，A_r(z_s,r)表示水位為z_s,r時(shí)下標(biāo)為r的控制體的斷面面積，A_*(z_s,r)表示水位為z_s,r時(shí)擴(kuò)大的控制體的斷面面積，Δx為空間步長(zhǎng)；對(duì)于上述擴(kuò)大的控制體，在計(jì)算其與兩側(cè)普通的控制體之間的數(shù)值通量時(shí)，不考慮支流的存在，然后按下式將A_*演進(jìn)到下一時(shí)間步：

$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為時(shí)間步長(zhǎng)，Δx為空間步長(zhǎng)，表示當(dāng)前時(shí)刻下標(biāo)為r的控制體右側(cè)的數(shù)值通量的第一個(gè)分量，表示當(dāng)前時(shí)刻下標(biāo)為r的控制體左側(cè)的數(shù)值通量的第一個(gè)分量，表示當(dāng)前時(shí)刻下標(biāo)為r的控制體的源項(xiàng)的第一個(gè)分量，表示當(dāng)前時(shí)刻擴(kuò)大的控制體的斷面面積，表示下一時(shí)刻擴(kuò)大的控制體的斷面面積；這里的中不包括干支流倒回灌凈單寬流量q_l，因?yàn)榈够毓嗟乃魇菙U(kuò)大的控制體內(nèi)部的流動(dòng)；通過(guò)式(7)得到之后，干支流交匯處的水位即下標(biāo)為r的控制體中的水位z_s,r就可以由式(6)得到，那么原來(lái)干流中第r個(gè)控制體內(nèi)的斷面面積即為

5.根據(jù)權(quán)利要求2所述的多沙河流水庫(kù)渾水明流與異重流耦合模擬方法，其特征在于：

上述步驟(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₀為交匯區(qū)干流異重流厚度；ξ為阻力損失系數(shù)；J為支流底坡；L為潛入段長(zhǎng)度；η＝(γ-γ₀)/γ，γ₀和γ分別為清渾水容重；g為重力加速度。

6.根據(jù)權(quán)利要求1所述的多沙河流水庫(kù)渾水明流與異重流耦合模擬方法，其特征在于：

上述步驟(3)中，潛入判別條件為：

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

式中：Fr_p為潛入點(diǎn)弗汝德數(shù)，S_v為渾水水流的體積比含沙量；如果由某控制體的U^m+1計(jì)算的弗汝德數(shù)Fr小于潛入點(diǎn)弗汝德數(shù)Fr_p，即認(rèn)為滿足潛入判別條件。

7.根據(jù)權(quán)利要求1至6中任一項(xiàng)所述的多沙河流水庫(kù)渾水明流與異重流耦合模擬方法，其特征在于：

上述步驟(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表示時(shí)間；x表示沿程距離；T是守恒變量；G是數(shù)值通量；R是源項(xiàng)；A_t，Q_t，C_t是異重流層的面積，流量和含沙量；B_t是清渾水層交界面的寬度；ρ_w為清水的密度；ρ_t，U_t是異重流層的密度與斷面平均流速；h_t是異重流層厚度；E_t和D_t是異重流與床面的泥沙交換通量；e_w是清水摻入系數(shù)；p為床沙的孔隙率；g為重力加速度；g′＝(ρ_t-ρ_w)g/ρ_t為有效重力加速度；S_b為底坡；S′_f為考慮交界面阻力后的綜合阻力坡度；h_ct為異重流斷面形心高度；A_b表示基準(zhǔn)高程以上的河床橫斷面面積；為沖淤導(dǎo)致的斷面面積的變化速率；q_tl表示考慮底坡的異重流倒灌流量；

計(jì)算源項(xiàng)R時(shí)，其第二個(gè)分量M_t中的水面梯度用步驟(2.2)中得到的計(jì)算。

8.根據(jù)權(quán)利要求7所述的多沙河流水庫(kù)渾水明流與異重流耦合模擬方法，其特征在于：

上述步驟(4)中，由顯式離散格式計(jì)算得到：

$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為時(shí)間步長(zhǎng)，Δx為空間步長(zhǎng)，和表示當(dāng)前時(shí)刻第i個(gè)控制體右側(cè)和左側(cè)的數(shù)值通量，表示當(dāng)前時(shí)刻第i個(gè)控制體的源項(xiàng)，和表示當(dāng)前和下一時(shí)刻第i個(gè)控制體的守恒變量。

9.根據(jù)權(quán)利要求7所述的多沙河流水庫(kù)渾水明流與異重流耦合模擬方法，其特征在于：

上述步驟(4)中，的計(jì)算方法如下，對(duì)式(12)直接按照時(shí)間前差進(jìn)行離散：

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

將上式得到的沖淤量ΔA_b分配到異重流段的斷面結(jié)點(diǎn)上，就得到了