高能河流二維輸沙模擬在渭河上的應用

董禮瑋

（西安水文水資源勘測中心,陜西 西安 710000 ）

1 引言

地下水是農(nóng)村和城市地區(qū)最重要的水資源之一。由于淡水資源與河流在數(shù)量和質量上的強烈交換,淡水資源極其脆弱。隨著與氣候變化相關的水資源的減少和需求的增長,大多數(shù)國家將地下水和地表水視為單一資源[1-2]。事實上,這兩個實體是相互關聯(lián)的,根據(jù)其數(shù)量或質量,它們之間關系密切。河流形態(tài)動力學是對控制和維持河流環(huán)境的物理過程的重要研究。因此,了解這些復雜的形態(tài)變化是了解河流和含水層之間交換區(qū)的關鍵問題。河流匯流對于河流系統(tǒng)中的水和沉積物的路由以及形態(tài)演變非常重要[3-4]。作為河網(wǎng)連接的河道中的分流現(xiàn)象通過這些交叉點相互影響,比單個河道中的分流現(xiàn)象復雜得多。對于沙質河流而言,在高含沙匯合洪水期間,其形態(tài)動力學在時間和空間域上都會發(fā)生顯著變化,這意味著水流運動和泥沙輸移之間存在著積極的相互作用,這給河流模型的開發(fā)帶來了更大的挑戰(zhàn)。因此,利用水動力模型模擬高濃度匯流洪水仍然是一項困難的任務。

2 研究區(qū)概況

2.1 渭河下游

渭河下游是黃河最大的支流。研究段長134.4 km,寬6 km~10 km,河道縱坡0.3‰~0.6‰,北洛河在潼關上游12 km~20.5 km處注入渭河。三門峽水庫蓄水后10 年間,渭河下游河道曲流頻繁,北羅河入口位置發(fā)生突變。沖積含水層是當?shù)鼐用竦闹饕嬘盟Y源之一。事實上,最近幾年的統(tǒng)計結果提到,有5000 萬m3的水是從這種非承壓含水層中抽取的[5]。幾十年以來,人類活動（沉積物提取、城市化等。）影響著地下水位的變化。

2.2 渭河下游形態(tài)動力學的關鍵問題

渭河下游的農(nóng)業(yè)發(fā)展、城市化和沉積物提取導致了河流形態(tài)變化。事實上,為改善渭河下游農(nóng)業(yè)發(fā)展等,修建了許多水利工程。修筑堤壩是為了保護人口和農(nóng)業(yè)用地免受洪水侵襲；為降低河床,從河床中提取了1500萬m3（相當于150年的自然沉積物供應）沉積物。此外,還建造了攔水堰,通過地表水下滲來增加地下水位。河床從150 m延伸到280 m,這導致了水流速度的增加。因此,在下游沿河許多地方可以觀察到河岸侵蝕不斷發(fā)生,這些物理過程影響了河流泥沙和水流之間的運移情況。

除了人為活動引起的河道形態(tài)變化外,河床的物理動力受洪水期影響,導致沉積物重組。為了解形態(tài)動態(tài)及其對河流含水層交換區(qū)的影響,必須建立數(shù)值模型。本文介紹了實現(xiàn)數(shù)值模型的方法和用于模擬渭河下游河道形態(tài)變化。

3 數(shù)據(jù)和方法

河流演變的建模需要計算模塊來模擬流體動力學和床層演化。如果以適當?shù)姆绞焦芾黹_發(fā)過程,2D模型可以將演變過程進行可視化。

3.1 數(shù)據(jù)

2009年~2013年期間渭河下游河勢變化頻繁,不僅河流形態(tài)發(fā)生了變化,洪水頻發(fā),導致渭河流域徑流量呈枯豐交替。利用2011年和2013年的兩個數(shù)字高程模型對渭河下游的形態(tài)動態(tài)進行評價,這種高分辨率和全維度的河道地形模型是開發(fā)表征河流特征和理解水流挾沙運動機制的新方法,為河流輸沙數(shù)值模擬提供了強有力的信息。2011年洪水事件引起的形態(tài)變化可以通過比較兩個地面模型來識別（圖1）。雖然這種方法在幾何角度方面有很多優(yōu)點,但是受數(shù)據(jù)的準確性影響,這種方法也有一定的局限性。

圖1 2011年洪水事件中識別沉降和侵蝕的比較

3.2 流體動力學模型

MIKE 21 FM的HD模型,通過求解質量和動量的二維淺水方程,進行了流體動力學模擬。為了使用有限體積數(shù)值方法求解這些方程,使用非結構化網(wǎng)格來定義研究區(qū)域的地形。模擬在單層流體中產(chǎn)生不穩(wěn)定的二維流動（垂直均勻）。流量和水位變化由Saint-Venant方程（垂直方向上積分的質量和動量守恒）描述,如下所示：

式中：h是水深,m；w=R-I是凈流量,m/s；z是表面高度,m；（p,q）=（hu,hv）分別是x和y方向的通量密度,m2/s；g是重力加速度,m/s2；pa是大氣壓,kg/（m·s2）；pw是水密度,kg/m3；τxx,τyy,τxy是有效剪切應力的組成部分。

3.3 泥沙輸送模型

使用MIKE 21 FM中的輸沙模塊（ST）進行了床層演變的模擬,該模塊計算沉積物輸送能力,河床層水平輸沙的初始速率以及由于水流運動引起的非粘性沉積物的形態(tài)變化。沉積物輸送計算基于流體動力學條件和沉積物性質。ST模塊考慮了河床層負荷和懸浮負荷。河床層負荷由剪切應力或每單位流量控制,并與流動瞬間反應。懸浮負荷是與流量相比,運輸中的相位滯后。河床荷載不需要對流色散模擬,但必須考慮兩個重要因素：河床剪切應力方向的偏差和傾斜河床的影響[5]；另一方面,流體中非粘性懸浮沉積物的建模可以通過體積沉積物濃度的運輸方程來描述。

式中：Sbl為推移質,m3/s；Ssl是懸移質,m3/s；T是取決于臨界摩擦速度和有效摩擦速度的無量綱傳輸階段；D*是無量綱粒子參數(shù)；S是沉積物的相對密度；Ca體積床濃度,m；f基于床層荷載的懸浮荷載修正系數(shù)；V是速度,m/s。

確定水平變化的關鍵參數(shù)是元素單元中心的床層水平變化率。河床的這種演化是用沉積物連續(xù)性方程Exner方程得到的：

式中：n是河床的孔隙度；z是河床的高度,m3；Sx是河床的推移質或者總的推移質沿著x方向,m2/s ；Sy是河床的推移質或者總的推移質沿著y方向,m2/s；ΔS是泥沙的匯或者源。

根據(jù)估計的床面變化率,通過形態(tài)學模擬對床面進行連續(xù)更新。泥沙輸送部分的目的是比較建立一個模型的合適方法,以MIKE 21 FM和ST模塊能夠代表渭河下游河道的形態(tài)變化,以及用于求解河床形態(tài)變化的方程。

4 結果和討論

4.1 流體動力學結果

首先,比較不同幾何形狀對流體動力學模擬的影響,根據(jù)網(wǎng)格分辨率和作業(yè)管理的需要,討論了5 種分辨率的網(wǎng)格模擬。為了得到一個高效的MIKE 21 FM 模型,對同一個洪水事件進行了多個網(wǎng)格的模擬,研究區(qū)域采用不同分辨率（5 m、10 m、15 m、20 m、25 m）的三角形離散化表示。然后將三角形網(wǎng)格與四邊形網(wǎng)格進行了比較。最后,利用 villemonte 堰流公式實現(xiàn)了水工建筑物的施工計算,這是一個經(jīng)驗公式[6]。本文選擇2015年10月3日~2015年10月6日一系列的洪水事件進行模型的測試。模型與觀測資料的比較（圖2）表明10 m的分辨率是最準確的。采用三角形離散化方法模擬溢流堰的流動更為有效（圖3）。模型中溢流堰結構化的實現(xiàn)需要許多參數(shù),并不能真正代表堰效應（圖4）。最后,將三角形網(wǎng)格,結構的三角形網(wǎng)格和四邊形網(wǎng)格的模擬速度進行比較,得出三角形網(wǎng)格以及四邊形網(wǎng)格模擬效果要優(yōu)于結構三角形網(wǎng)格模擬。

圖2 不同分辨率的模型與觀測資料的水深比較

圖3 模型中采用三角形網(wǎng)格與四邊形網(wǎng)格計算比較

圖4 模型中采用三角形網(wǎng)格與結構三角形網(wǎng)格計算比較

圖5 各網(wǎng)格的速度比較

4.2 河床形態(tài)演變

渭河存在的主要問題之一是下游淤積。泥沙輸送模型考慮了推移質和懸移質的影響。因此,2011年洪水事件的模擬再現(xiàn)了包括堰在內(nèi)的河床水位。流體動力學研究結果表明,結構的實施不同程度地影響流場。針對堰體對地貌的影響,通過對三種模型的比較,找到解決方案: 一是模型以均勻的顆粒尺寸定義,二是模型以結構化網(wǎng)格劃分實現(xiàn),三是模型以非均勻的顆粒尺寸,可定義為阻止破壞堰體的推移質。三種模型的河床變化比較（圖6）表明,對于渭河河床模擬的這一部分來說,結構的實施是模擬河床物理過程的最合適的方法。河床演變的計算是基于水動力學和泥沙特性。在這里,只有沉積物的特性不足以描述河床形態(tài)變化。

圖6 各模型模擬的床層變化比較

5 結語

渭河由于頻繁的高含沙洪水,河道演變十分活躍。基于水動力和形態(tài)過程耦合在高濃度合流洪水模擬中的重要性,本文提出2D水力模型對河流內(nèi)泥沙輸移過程進行模擬。利用水動力學模型的離散化方法,建立了泥沙輸送模型。模擬結果表明,三角形網(wǎng)格（10 m分辨率）雖然是描述堰上流場的最佳方法,但不適合模擬堰下的河床演變。因此,它需要在實施模擬的結構的基礎上提高模擬性能,同時結合相關河流輸沙經(jīng)驗公式法,應該著重分析該模型使用不同泥沙輸送公式的靈敏度,并將這種方法應用河于道的蜿蜒部分。