MacCormack法應用于輸沙計算的研究

高波 黃金柏

摘 要：本研究探討將MacCormack方法應用于河流的推移質計算，為流域輸沙計算提供科學的數值方法。研究基于MacCormack法的差分方案，將水流計算分為預測計算和修正階段，而后將推移質輸沙計算公式與MacCormack計算方案耦合，并應用于設定條件的事例計算。計算結果形象地再現了河流在改變流量條件下輸沙量隨時間的變化，該方法適用于計算常流和變流量條件下的河道輸沙量、也可用于河床演變，河道疏浚的計算。

關鍵詞：MacCormack法，有限差分，推移質，輸沙

中圖分類號：TV148文獻標識碼： A 文章編號：

1前言

流域的水沙過程包括降雨產流、匯流的全過程, 以及泥沙侵蝕、輸移、沉積的不同運動形態. 對這一物理過程的全面描述有其必要性和現實意義[1]。輸水的同時輸沙, 對河流起著泄洪排沙、維持河道正常演變的作用[2]。泥沙運動既應滿足一定的力學規律, 同時也是一種隨機現象[3]。輸沙計算，一般先求解水流條件, 得到水深、流速、流量等水力因子, 將水力因子代入輸沙方程中, 得到含沙量或輸沙率過程[4]。

在現代河流泥沙物理模型及數學模型中,經常遇見輸沙計算問題, 但由于河流動力學所面臨的自然現象, 是邊界復雜多變的天然河流及水流泥沙運動的多相體, 許多問題的運動規律尚無定論, 特別是水流挾沙力的計算, 人們還不能全面地把握它的力學機制, 因此, 多少年來, 泥沙界的學者們對水力輸沙特性在不懈地探索[5]。崔俠等（1987）從沖淤過程的機理出發, 推導出了二度恒定非均勻流中泥沙擴敬方程, 求出了含沙量沿程變化規律, 著重討論了水流挾沙力與含沙量的關系[6]。張洪武等（1992）從水流能量消耗和泥沙懸浮力之間的關系出發, 考慮了泥沙存在對卡門常數和泥沙沉速等的影響, 給出了半經驗半理論的水流挾沙力公式并在黃河的輸沙計算中應用[7]。胡海明等（1996）較為系統地研究了非均勻沙運動機理及輸沙率的計算方法, 在建立泥沙運動交換模式的基礎上, 指出現有非均勻泥沙起動流速公式存在的不足并對此作了修正[8]。黃永健等（1997）結合灌渠的結構特點及水流泥沙特性, 對一維非耦合、非均勻、不平衡輸沙基本方程組進行了合理的簡化, 根據實體模型試驗及實測資料, 確定了自流與提灌時干、支渠分流的不同分沙分水比并應用于黃河下游的輸沙計算[9]。迄今為止的輸沙計算方法普遍存在兩個不足，一是經驗型公式（模型）較多，其應用范圍常被限制在研究開展的流域或河段；二是計算精度普遍較低，對輸沙過程的模擬精度不高。本研究探討將MacCormack方法用于推移質的輸沙計算，為河流輸沙，河床演變以及河道疏浚等提供科學的計算方法。

2 MacCormack方法

MacCormack方法由Robert W. MacCormack在1969年提出，是計算流體力學中用來求解雙曲形偏微分方程數值解的一個普遍離散方法，為二階有限差分方法。該計算方法非常簡潔，易于理解和實現編成計算，應用MacCormack方法計算時包括兩個階段：預測階段和修正階段，該方法非常適合線性方程的求解。由于采用二階差分對預測階段和修正階段進行計算，數值計算結果的精度較高。對于非線性方程，應用MacCormack方法也可得到較好的結果。

3 算法建立

3.1 基礎方程式

用于徑流輸沙過程計算的基礎方程式包括地表徑流連續方程式，運動方程式以及輸沙量計算公式[10]。

·地表徑流連續方程式（一維穩定流）

（1）

·運動方程式

（2）

·輸沙量計算公式（Ashida / Michiwue 經驗公式）

（3）

式中，Δt為計算的時間步長，1s；Δx為河流縱向計算的單位步長，50 m；b為河道寬度，m；Q為流量；g為重力加速度，m/s2；h為水深，m；ib為河道平均坡度；ie為水力坡度；qbi為單位時間內單寬河道輸沙量，m3/s；s為沙在水中的比重，kg/m3；d為流沙的平均粒徑，0.5 mm，為推移質；τ*c為界限推移力，kg；τ*為推移力，kg；u*c為界限摩擦速度，m/s；u*為摩擦速度，m/s。

3.2 有限差分

計算時，需要依據地表徑流連續方程式及運動方程式耦合輸沙計算在時間上離散化，即有限差分。因為計算依賴于一定的初始條件和邊界條件，所以預測階段采用后退差分法，修正階段采用中間差分法[11,12]。

·預測階段差分式

（4）

·修正階段差分式

（5）

式中，i為計算的柵格編號，n為計算的時間編號，A為過流面積，m2; kv為人工粘性系數，0.00001；其它因子同上。

3.3 計算程序編譯

實現數值計算的程序采用計算機高級語言Fortran編譯，此部分內容在此略去。

4 計算事例

4.1 計算條件

圖1 計算河道長度及形狀

計算河道的形式如圖1所示，河口前長度為2 km，河床坡度ib為0.001。自河口以下計算長度為1k m，水流自河口開始以與河道中心線成20o角向兩側擴散，河口以下坡度為0.01， 河口的標高為0。河道年均流量為 50 m3/s，計算開始5年后發生了100年一遇的洪水，其設計流量為2500 m3/s，洪峰持續時間為5小時。求解問題為5年后發生洪水前河床的形態，洪水發生后每隔1小時河床的演變情況。

4.2 結果及討論

利用上述的數值計算方法和題設條件進行計算，根據求解問題的要求，當年平均流量為50m3/s時 5年后的河床形狀、水面線以及初始狀態的河床形式如圖2所示。

圖2 河床的初始形態與5年后形態的對比

洪水發生后，每隔1小時河床的形態與水面線的數值計算結果如圖3所示。

圖3洪水發生后每1小時間隔的河床及水面線形狀

由圖2可知，因輸沙的影響，河床的初始形態與5年后的形態明顯不同，特別是河口部發生了明顯的堆沙現象。

圖3描述了洪水發生后，水面線及河床變化的過程，河床形態以及水面線在河口部以前呈現相對穩定的狀態，沒有發生明顯的變化，即河床演變在河口部上游區域表現并不明顯。在河口部，洪水發生后導致推移質隨時間逐漸前移，并在河口部堆積，其泥沙的堆積區域主要在距河口700 m范圍之內，而后泥沙堆積的程度逐漸變小，在距河口1000m之處幾乎沒有發生泥沙堆積現象。

利用數值計算計算的結果，可以推求輸入河口的泥沙量，即河口的堆沙量。其計算方法是比較河口以下發生洪水后的河床高程與計算初始時河床各段的高程（Δx=50 m），由二者的差值結合河道的寬度可近似計算出輸沙量。從而為河口部泥沙疏浚提供基礎數據。

5 結論

研究實現了MacCormack方法與推移質輸沙計算耦合，為河流推移質計算提供了一種實用的方法。

計算結果可以很好地模擬河道推移質移動的過程，利用MacCormack方法進行河流輸沙計算可以得到準確度較高的結果。

研究可為河流輸沙量的準確推求，河床演變以及河道疏浚等提供數值計算方法。