本發(fā)明涉及水利工程技術領域，具體的說是多沙河流水庫渾水明流與異重流耦合模擬方法。

背景技術：

水庫異重流是多沙河流水庫泥沙管理中引起普遍關注的問題。在水庫設計或運用階段對異重流問題考慮不足有可能影響其正常運轉甚至減損水庫壽命。目前對于水庫異重流在動力學機制上仍有認識不足，模型發(fā)展現(xiàn)狀還不能滿足實際工程需要。此外在水庫調水調沙過程中，渾水明流段的水沙運動對于異重流的形成和運動有重要影響，在此種情況下對整個庫區(qū)的水沙運動進行數值模擬存在很大的挑戰(zhàn)。

現(xiàn)有異重流模型通常分為解析異重流模型、一維和平面二維異重流模型，以及立面二維和三維異重流模型。

解析異重流模型只能對于無交界面摻混和無床面泥沙交換的理想情況給出近似解。目前解析模型的研究成果只能用于規(guī)則斷面形態(tài)下，具有特定上游邊界條件下(即總體積固定和恒定入流兩種情況)的異重流傳播計算。立面二維和三維異重流模型可以模擬出流速、含沙量沿垂向的不均勻分布。這類模型構建時大部分涉及到紊流封閉，也有一部分使用直接數值模擬。由于目前對于紊動現(xiàn)象的理解不足，不同的學者提出了許多紊流模式理論，相應的模型對于同一異重流過程給出的預測結果差異很大。此外立面二維和三維異重流模型計算成本很高，大部分應用仍局限于對水槽實驗過程的模擬。

一維和平面二維異重流模型是在對Navier-Stokes方程深度平均的基礎上得來的，雖然不能求解水沙要素的垂向分布，但是計算速度快且能夠確定工程應用中關心的傳播速度、沖淤量等問題。傳統(tǒng)的一維和平面二維異重流模型由于不考慮明流段與異重流段水沙運動的耦合關系，只能用于模擬潛入點下游的異重流傳播過程而很難模擬異重流的產生。另一方便對于庫區(qū)存在支流的多沙河流水庫，現(xiàn)有的模型沒有考慮干支流倒回灌的影響。

技術實現(xiàn)要素：

本發(fā)明的目的在于提供一種多沙河流水庫渾水明流與異重流耦合模擬方法，能預測水庫異重流的形成時間和潛入位置，以及異重流發(fā)生期間明流段和異重流段的水沙輸移過程，干支流倒回灌過程，水庫排沙比等，為可持續(xù)的水庫泥沙管理提供科學依據和技術支撐。

本發(fā)明的原理為：多沙河流水庫中的河床沖淤演變受到渾水明流段和異重流段水沙運動的共同作用，同時受到干支流倒灌的影響，有必要在模擬時全面考慮這些因素的影響。入庫泥沙在明流段的輸移或明流段的沖刷為異重流的形成提供了泥沙，洪水在明流段的演進影響了異重流中期潛入位置的變動，支流口門異重流的厚度和含沙量決定了異重流倒灌的強度，而支流向干流的回灌又影響了明流段和異重流段的水位。因此，本發(fā)明結合控制方程交替求解模式、異重流潛入判別條件、零維水庫法和異重流倒灌流量公式，提出多沙河流水庫渾水明流與異重流耦合模擬方法。

本發(fā)明多沙河流水庫渾水明流與異重流耦合模擬方法的數值計算方法，包括如下步驟：

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

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

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

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

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

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

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

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

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

步驟(4.2)，計算得到下一時刻的和N為模擬區(qū)域內的斷面?zhè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表示某一斷面上劃分的子斷面?zhèn)€數；h為水深；h_c為過流斷面形心高度；q_l為干支流倒回灌凈單寬流量，規(guī)定從干流到支流為正，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₀為交匯區(qū)干流異重流厚度；ξ為阻力損失系數；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分配到異重流段的斷面結點上，就得到了

本發(fā)明多沙河流水庫渾水明流與異重流耦合模擬方法的有益效果在于：

(1)采用了水沙耦合的渾水明流與異重流控制方程，能夠反映高含沙量或急劇床面沖淤對洪水演進的影響；

(2)可預測異重流的潛入位置和形成時間、異重流傳播過程，以及異重流到達壩前后渾水層抬高和排出的過程；

(3)可根據明流段和異重流段的模擬結果，計算泥沙沖淤總量，以及沖淤量在干流不同庫段和干支流間的分布，可以計算水庫排沙比。

附圖說明

圖1為本發(fā)明實施例計算方法的流程圖；

圖2為本發(fā)明實施例橫斷面灘槽節(jié)點劃分示意圖；

圖3為本發(fā)明實施例渾水明流段的橫斷面示意圖；

圖4為本發(fā)明實施例零維水庫法示意圖；

圖5為本發(fā)明實施例異重流倒灌的示意圖；

圖6為本發(fā)明實施例異重流段的橫斷面示意圖；

圖7為本發(fā)明實施例模擬的不同時刻庫區(qū)水位和交界面高程分布,(a)t＝317.57h，(b)t＝378h。

具體實施方式

下面結合附圖和實施例，對本發(fā)明做進一步說明。如圖1-圖7所示，多沙河流水庫渾水明流與異重流耦合模擬方法，具體步驟如下：

步驟1：根據水庫干支流實測斷面地形劃分干流斷面灘槽節(jié)點，計算支流底坡J和水位庫容關系V(z_s)；將模擬區(qū)域劃分為一組控制體；

由各個實測斷面的起點距高程數據繪制斷面形態(tài)圖，如圖2所示，圖中橫向底坡的轉折點定義為灘槽分界點，兩個灘槽分界點之間的節(jié)點定義為主槽節(jié)點，灘槽分界點至斷面左端點或右端點之間的節(jié)點定義為灘地節(jié)點，每兩個斷面節(jié)點之間為一個子斷面。如果某子斷面兩端均為主槽節(jié)點，則令該子斷面的糙率n_j等于主槽糙率n_MC；如果某子斷面兩端含有灘地節(jié)點，則令該子斷面的糙率n_j等于灘地糙率n_FP，一般情況下n_FP略大于n_MC。

對于庫區(qū)內每條支流，根據實測斷面數據確定支流底坡J和深泓線最低點，在最低點高程與水庫校核洪水位之間劃分若干水位，計算各水位下相應的庫容，得到各個支流的水位庫容關系V(z_s)，z_s表示水位。

模擬區(qū)域即為需要研究分析的庫區(qū)范圍，其下游邊界到達水庫壩址，上游邊界到達水庫回水范圍以上可獲得入庫水沙測量數據的測站位置?？刂企w即為數值模擬時對空間上連續(xù)的地形以及水位，流量等物理量進行離散化表示的基本單元。劃分控制體時，每個實測斷面對應一個控制體，每兩個相鄰實測斷面的中點位置定義為與這兩個斷面對應的控制體的交界面。

步驟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。

本發(fā)明采用考慮干支流倒回灌的水庫渾水明流與異重流耦合模型，模型包括渾水明流與異重流兩組控制方程，渾水明流控制方程如下：

式中：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表示某一斷面上劃分的子斷面?zhèn)€數；h為水深；h_c為過流斷面形心高度；q_l為干支流倒回灌凈單寬流量，規(guī)定從干流到支流為正，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₀為交匯區(qū)干流異重流厚度；ξ為阻力損失系數；J為步驟1中得到的支流底坡；L為潛入段長度；η＝(γ-γ₀)/γ，γ₀和γ分別為清渾水容重。部分變量定義如圖5所示。

步驟2.3：完成渾水明流控制方程組中式(1)第二三行計算得到U^m+1，完成渾水明流控制方程組中式(4)的計算得到河床沖淤量，將沖淤量分配到斷面各個節(jié)點從而修改(更新)河底高程z_b，修改后的高程暫記為

當床面發(fā)生淤積時，將淤積量平均分配到斷面各個節(jié)點；當床面沖刷時，將沖刷量只分配到主槽節(jié)點上。

步驟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和斷面形態(tài)計算出來。

步驟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為模擬區(qū)域內的斷面?zhèn)€數。

步驟5：計算干支流沖淤量及水庫排沙比，包括步驟5.1-5.2。

步驟5.1：計算庫區(qū)沖淤量和水庫排沙比。

模擬時段內入庫泥沙量W₁由實測入庫水沙資料計算。出庫泥沙量W_N根據模型輸出的結果計算，即：

式中：Δt為時間步長，M為模擬結束時的時間步數，m是時間步長數，取值1,2,3,…M，分別是干流出口斷面的流量和含沙量。

庫區(qū)沖淤總量為W₁-W_N，水庫排沙比為W_N/W₁。

步驟5.2：計算支流淤積量。

對于水庫內任意一條支流，按下式計算支流淤積量W_tl：

式中：Δt為時間步長，M為模擬結束時的時間步數，m是時間步長數，取值1,2,3,…M，b是支流口門寬度，即為式(7)計算的異重流倒灌流量，為倒灌異重流的含沙量，一般可取交匯處的干流異重流含沙量。

圖7表示的是使用本發(fā)明方法所模擬的小浪底水庫2004年調水調沙過程中產生的異重流。圖7(a)顯示的是異重流形成4小時后整個庫區(qū)的縱剖面，模型預測的異重流潛入點距壩約46km，與當年實測結果一致。潛入點下游的清渾水交界面將異重流段的水體分為上下兩層。在圖7(a)所示時刻，入庫流量Q_in為1014m³/s，壩前水位為233.79m。在圖7(b)所示時刻，異重流已經到達壩前，壩前水位下降到229.08m，但是由于入庫流量Q_in增加到4434m³/s，渾水明流段水深有所增加，同時異重流潛入點向下游移動了約2km。

應當理解的是，本說明書未詳細闡述的部分均屬于現(xiàn)有技術。

應當理解的是，上述針對較佳實施例的描述較為詳細，并不能因此而認為是對本發(fā)明專利保護范圍的限制，本領域的普通技術人員在本發(fā)明的啟示下，在不脫離本發(fā)明權利要求所保護的范圍情況下，還可以做出替換或變形，均落入本發(fā)明的保護范圍之內，本發(fā)明的請求保護范圍應以所附權利要求為準。