摘要：本文介紹了一維非恒定泥沙數學模型的原理和計算方法，并且對于泥沙數學模型中的一些關鍵技術等問題進行了較詳細的討論。利用已建立的一維非恒定泥沙數學模型，對黃河小北干流1981年洪水資料進行了模擬計算，計算結果表明：黃河干流、渭河、北洛河計算水位、流量過程線與實測值符合良好，可用于黃河的洪水演進計算。

關鍵詞：非恒定流 泥沙數學模型 河床變形

1 泥沙數學模型基本方程

明渠或天然河流常被考慮作為一維流動，根據洪水波運動的圣維南方程、泥沙連續方程和泥沙擴散方程，可以簡化推導出一維非恒定流泥沙數學模型的基本方程，其形式為水流連續方程

(1)

水流動量方程

(2)

泥沙連續方程

(3)

不平衡輸沙方程

(4)

挾沙力方程

S*=f(Q，A，B，ω，S)

(5)

其中 Q為流量，A為斷面面積，B為斷面寬度，Z為水位，K為流量模數，S為斷面平均含沙量，S*為斷面平均挾沙力，g為重力加速度，α為泥沙非平衡恢復飽和系數，ω為泥沙顆粒沉速，Ad為斷面沖淤面積，γ′s為泥沙干容重。

2 計算方法簡介

一維非恒定泥沙數學模型的計算采用非耦合方法，首先求解水流連續方程和動量方程，然后求解水流挾沙力、泥沙不平衡輸沙方程和泥沙連續方程，具體求解過程如下。

2.1 水流方程的求解

首先利用Preissmann隱式差分格式將水流連續方程和動量方程轉化為差分方程，然后對差分方程進行線性化，在線性化過程中，略去增量的乘積項，最后得到以下線性方程組

A1jΔQj+B1jΔZj+C1jΔQj+1+D1jΔZj+1=E1j

(6)

A2jΔQj+B2jΔZj+C2jΔQj+1+D2jΔZj+1=E2j

(7)

其中 Aij、Bij、Cij、Dij、Eij（i=1，2）為第j單元河段差分方程的系數（j=1，2，......，N-1，其中N為斷面個數）。

給定邊界條件

ΔQ1=Q1 n+1-Qn1=Q1(tn+1)-Q1n

(8)

ΔZN=ZNn+1-ZNn=ZN(tn+1)-ZNn

(9)

方程(6)、(7)式及邊界條件，共有2N個未知數，2N個方程，可以求解。由于差分方程中的系數包含有未知數，方程求解不能直接求出未知變量，因此方程求解時必須進行迭代處理。下面給出用追趕法求解的步驟，追趕方程為

ΔQj=FjΔZj+Gj

(10)

ΔZj=HjΔQj+1+IjΔZj+1+Jj

(11)

其中Hj、Ij、Jj、，Fj、Gj為追趕系數。

2.2 水流挾沙力、動床阻力、河寬變化、床沙級配調整采用文獻[2]的方法

2.3 不平衡輸沙方程求解

利用迎風格式，將(3)式離散為差分方程，整理后得

(12)

當Q≥0時，利用上邊界條件，自上而下計算各斷面含沙量；當Q<0時，利用下邊界條件由下至上計算各斷面含沙量。

2.4 河床變形及淤積量計算

由式(3)與式(4)相減可以得到河床變形方程為

(13)

將上式寫成差分方程，各斷面淤積面積為

(14)

第j河段的淤積量為

ΔWj=(ΔAj+ΔAj+1)Δxj/2

(15)

3 水流內邊界的處理

水流內邊界是指河道的幾何形狀的不連續或水力特性的不連續點。例如，河流的匯合點、河流分流、局部河段內生產堤決口等等。在這些內部邊界處，圣維南方程組和單一河道泥沙不平衡輸沙方程等都不再適用，必須根據其水力特性作特殊處理。內邊界條件通常包含兩個相容條件：即流量的連續條件和能量守恒條件（或動量守恒條件）。本模型主要考慮了以下幾個類型的內邊界處理。

3.1 水沙的匯入或匯出

如圖1所示，假設匯入或匯出點上下斷面滿足以下條件

3.2 支流從干流分流

如圖2所示，干流和支流上斷面之間應滿足連續方程和能量方程

圖1 水沙的匯入或匯出示意圖 Sketch of water/sediment inflow or outflow

圖2 支流從干流分流示意圖 Outflow from the main stem

圖3 支流匯入干流示意圖 Inflow into the main stem

支流從干流分流時，將干流分流斷面按干支流流量比分為兩部分，忽略時變項，將方程(4)式直接寫成差分形式，求得干流和支流下游斷面的含沙量為

其中φ為分流系數，由實測資料確定。

3.3 支流匯入干流

當支流匯入干流時，與支流從干流分流類似，干流和支流上斷面之間也滿足連續方程和能量方程

將匯流斷面按干支流流量比分為兩部分，忽略時變項，將方程(4)式直接寫成差分形式，求得匯流斷面的含沙量為

其中 φ為分流系數，由實測資料確定。

4 泥沙數學模型的應用

4.1 計算區域和時段

計算區域為龍門、華縣、河津、狀頭至潼關河段，黃河干流上有渭河和汾河匯入，在渭河上有北洛河匯入。本模型同時模擬黃河干流、渭河、北洛河三條河流的洪水演進過程，各匯入點作為內邊界處理，汾河僅作為已知水、沙過程線匯入黃河干流。上邊界條件為龍門、華縣、河津、NFDA1頭各站，水沙為進口控制的已知條件，下邊界條件為潼關站出口控制水位。

計算時段為1981年汛初第一場洪水，洪水時間為7月3號至7月14號，洪水持續時間為12天，龍門最大流量為6400m3/s、最大含沙量為298.0kg/m3；華縣最大流量為970m3/s、最大含沙量為117.0kg/m3。

原始大斷面資料采用1981年汛前實測大斷面資料。

4.2 計算結果與實測值的比較

圖4為潼關站計算流量與實測值的比較、計算含沙量與實測值的比較。由圖4可以看出，計算的潼關出口流量過程線與實測過程線比較符合，計算洪峰最大值和相應洪峰傳播時間與實測值比較接近;計算潼關出口含沙量過程線與實測值也是比較一致的。由圖4分析得出，本模型能夠比較好地模擬出三條河流的洪水傳播過程和泥沙沖淤調整過程。

圖4 潼關站流量、含沙量計算值與實測值的比較

Comparison of calculated and measured discharge and concentrations at Tongguan station

圖5 華陰站流量、水位計算值與實測值的比較

Comparison of calculated discharge and water levels with measured at Huayin station

圖5為華陰站計算水位與實測值的比較、計算流量與實測值的比較。由圖中可以看出，計算水位、流量過程線與相應實測過程線比較符合。同時也看出，在洪水初期，華陰站出現了倒灌現象，流量出現了負值，本模型也比較好地模擬出了黃河干流倒灌渭河的現象。

圖6為朝邑站計算水位與實測值的比較、計算流量與實測值的比較。可以看出，計算水位、流量過程線與相應實測過程線比較一致。由于計算北洛河各斷面的水位和流量，是通過一系列內邊界條件計算得到的，因此也可以說明本模型所采用的內邊界處理技術是比較合理的。

5 結語

本文較詳細地介紹了一維非恒定泥沙數學模型的原理和計算方法，并對泥沙數學模型中的一些關鍵技術等問題進行了討論。利用一維非恒定泥沙數學模型，對黃河小北干流1981年洪水資料進行了模擬計算，結果表明：模型比較好地模擬了黃河干流、渭河、北洛河洪水傳播特性，洪峰傳播時間和洪峰最大值與實測值比較一致；模型對于匯流區水流和泥沙計算所采用的處理技術也是比較合理的，對于干流倒灌支流或支流倒灌干流的現象，模型也是比較能夠適應的，因此，模型可用于多條河流的洪水演進計算。

圖6 朝邑站流量、水位計算值與實測值的比較

Comparison of measured and calculated discharge and water levels at Chaoyi station

參 考 文 獻

[1] Lyn，D.A. Goodwin，P. Stabling of a eneral Preissmann scheme. J. of Hydra. Eng. Vol. 113，No.1，1987.

[2] 梁國亭，張仁。黃河小北干流一維泥沙沖淤數學模型，人民黃河，1996年9月。