高波 黃金柏

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

關(guān)鍵詞：MacCormack法，有限差分，推移質(zhì)，輸沙

中圖分類號：TV148文獻(xiàn)標(biāo)識碼： A 文章編號：

1前言

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

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

2 MacCormack方法

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

3 算法建立

3.1 基礎(chǔ)方程式

用于徑流輸沙過程計(jì)算的基礎(chǔ)方程式包括地表徑流連續(xù)方程式，運(yùn)動方程式以及輸沙量計(jì)算公式[10]。

·地表徑流連續(xù)方程式（一維穩(wěn)定流）

（1）

·運(yùn)動方程式

（2）

·輸沙量計(jì)算公式（Ashida / Michiwue 經(jīng)驗(yàn)公式）

（3）

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

3.2 有限差分

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

·預(yù)測階段差分式

（4）

·修正階段差分式

（5）

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

3.3 計(jì)算程序編譯

實(shí)現(xiàn)數(shù)值計(jì)算的程序采用計(jì)算機(jī)高級語言Fortran編譯，此部分內(nèi)容在此略去。

4 計(jì)算事例

4.1 計(jì)算條件

圖1 計(jì)算河道長度及形狀

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

4.2 結(jié)果及討論

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

圖2 河床的初始形態(tài)與5年后形態(tài)的對比

洪水發(fā)生后，每隔1小時河床的形態(tài)與水面線的數(shù)值計(jì)算結(jié)果如圖3所示。

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

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

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

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

5 結(jié)論

研究實(shí)現(xiàn)了MacCormack方法與推移質(zhì)輸沙計(jì)算耦合，為河流推移質(zhì)計(jì)算提供了一種實(shí)用的方法。

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

研究可為河流輸沙量的準(zhǔn)確推求，河床演變以及河道疏浚等提供數(shù)值計(jì)算方法。