馬 強，徐立榮，董曉知，徐 晶

(濟南大學 水利與環(huán)境學院，山東 濟南 250022)

黃河是世界上輸沙量最大和含沙量最大的河流，因此引黃灌區(qū)具有“引黃必引沙”的重要特點[1]。引黃灌區(qū)自開始引水以來，泥沙問題一直伴隨著灌區(qū)發(fā)展[2]。雖然自2002年小浪底水庫進行調水調沙試驗，并轉入正式生產運行后，黃河下游特別是山東境內的引黃條件發(fā)生了很大變化，如引黃能力下降，輸沙量變小等，但是引黃灌渠排沙問題仍未減輕，如何排沙依舊是黃河下游治理的重中之重，需要重新提出科學的水沙調度和合理的工程措施，從而延長沉沙池使用壽命，減少骨干渠道淤積，更好地發(fā)揮灌區(qū)的灌溉效益。

渠道是灌區(qū)灌排系統(tǒng)組成的重要部分，而渠道輸水輸沙能力的有效實現是整個灌區(qū)正常運行的關鍵[3-4]。為了解決灌區(qū)的泥沙問題，提高渠道的輸水輸沙能力，減少泥沙淤積，保證灌區(qū)的正常運行，人們通常以控制渠道泥沙淤積為目標，建立相應水沙數學模型。一般對河道泥沙淤積的數學模型研究多采用一維水沙模型，模擬斷面各水沙要素的總平均值[5]，相比之下，平面二維水沙模型能更清楚地表達河床的平面變形，因此被廣泛應用于河道水沙運動研究[6-7]。

目前，對于引黃灌渠輸沙渠道二維水沙的泥沙淤積研究成果較少。本文中以山東省聊城市位山灌區(qū)一支改造后的輸沙干渠——東輸沙渠作為研究對象，以控制輸沙渠渠道泥沙淤積為目標，利用水力學模型MIKE 21具有的簡單、快速、精確、可視性好等優(yōu)點[8]，構建輸沙渠二維水沙模型，在經過模型率定與驗證后，設計典型工況對水流特性及泥沙淤積分布規(guī)律進行模擬預測，篩選保水輸沙的最佳方案，為黃河下游引黃灌區(qū)干渠泥沙治理和利用提供參考。

1 研究區(qū)概況

山東省聊城市位山灌區(qū)地處黃河下游，是全國第五大灌區(qū)、山東省最大灌區(qū)[9-10]。灌區(qū)共包括四區(qū)一市五縣，土地總面積為5 380 km2，總耕地面積為3 793 km2，設計灌溉面積為3 387 km2[11]。灌區(qū)南靠黃河，北臨衛(wèi)運河，東與山東省德州市相接，馬頰河和徒駭河從灌區(qū)自西南向東北貫穿而過。灌區(qū)屬黃河沖積形成的平原，地勢平坦開闊，走向略微向東北方向傾斜[9，11]。

東輸沙渠于2019年改造完成，全長14.5 km，擔負著一干渠供水供沙的任務，具體位置如圖1所示。東輸沙渠下接進入東沉沙區(qū)，有沉沙池3個，已還耕1個，主要使用1#、3#池輪用，2個沉沙池長度分別為3.0、3.2 km。渠道斷面被改造成梯形復式斷面，為了使水流順暢、流態(tài)穩(wěn)定，渠道中心線保持不變。渠道比降仍采用原比降1/6 000，邊坡坡比為1∶2.0，復式斷面底寬為16.0 m，渠道開挖深度為1.5 m，淺灘寬度為2 m。

位山灌區(qū)地理位置地圖從國家測繪地理信息局標準地圖服務網站下載，地圖審批編號為GS(2019)3333(http://bzdt.ch.mnr.gov.cn/browse.html?picId=%224o28b0625501ad13015501ad2bfc0210%22)，經過ArcGIS10.3軟件數字化處理得到。圖1 位山灌區(qū)輸沙干渠圖

2 二維水沙模型建立及驗證

2.1 水動力模型理論

水力學模型MIKE 21中的水動力(HD)模塊是必需選擇的核心模塊，是所有模擬的基礎。該模塊建立在基于二維數值求解方法的淺水方程基礎上，在深度上集成了三向不可壓縮雷諾(Reynolds)平均Navier-Stokes方程，并服從Boussinesq假設和靜水壓假設[12-15]。

二維非恒定淺水方程組包括連續(xù)性方程，X、Y方向動量方程，即

連續(xù)性方程

(1)

X方向動量方程

(2)

Y方向動量方程

(3)

該模型是由動量守恒方程和水流連續(xù)方程組成，在水平面上可以使用笛卡爾坐標。

2.2 泥沙模型理論

水力學模型MIKE 21中的輸泥模塊(MT)結合了多粒徑級和底床分層，描述了黏聚性泥沙(淤泥或黏土)在波浪和水流作用下的沖刷、傳輸和沉積[16]。

1)基本方程。二維輸泥模塊的對流擴散控制方程為

(4)

2)泥沙沉積過程。泥沙在水體中通常處于沉積和懸浮狀態(tài)，泥沙淤積是泥沙顆粒從水體沉降至床面的輸移過程。當床面切應力小于淤積臨界切應力時，產生泥沙淤積；反之，床底發(fā)生沖刷。細顆粒泥沙的沉降速度取決于顆粒大小、溫度、泥沙濃度等。黏性泥沙的淤積的控制方程為

D=ωsρbPD，

(5)

式中：D為沉積率；ωs為泥沙沉降速度，m/s；ρb為泥沙的近底質量濃度，kg/m3；PD為泥沙沉降在河床的概率函數，

(6)

其中τb、τcd分別為河床剪切應力、河床淤積臨界剪切應力，N/m2。

3)河床的沖刷過程。河床剪切應力τb大于臨界沖刷剪切應力τce時，就會發(fā)生河床層的沖刷。床底沖刷程度表達式為

(7)

式中：SE為侵蝕率；E為床底沖刷程度，kg/(s·m2)；n為侵蝕程度。

2.3 模型的建立

2.3.1 模型的范圍與網格

模型范圍是上起位山閘口(東經116.122 905°，北緯36.141 67°)，下至東沉砂池入口(東經116.142 618°，北緯36.264 66°)。

河段沖淤演變計算所采用的網格劃分為普遍用于各種河道的非結構三角網格，渠道邊界條件簡單，因此不對網格進行局部加密，河道形態(tài)存在彎道河段，模擬計算時間設置為30 d。經過反復試算，網格劃分參數設置如下：三角網格邊長為6 m，網格總個數為35 012，總節(jié)點個數為19 919，最終網格文件由網格生成器生成，如圖2所示。

2.3.2 邊界條件及參數設置

模型計算采用上游位山引黃閘處流量及含沙量作為進口邊界，東沉沙池入口處水位作為下游出口邊界，模擬總時間步長為3 600 s，設置時間步數為700，空間和時間的積分都采用低階、快速運算(low order，fast algorithm)的求解格式來求解淺水方程。模型建立完畢后，進行水動力模塊和泥沙運輸模塊的參數設置，見表1。

表1 水動力模型邊界條件與模型參數設置

2.4 模型驗證

模型驗證的目的是檢驗所建模型與實際工況的相似度，對模型中相關參數的率定結果進行檢驗和修正。

2.4.1 水動力模塊驗證

位山灌區(qū)多年沒有發(fā)生洪水，只存在關閘不引水而斷流的情況，河道經歷了開挖、淤積、清淤、改造，處于不斷變化的過程，缺乏長系列配套的水文地形資料。對于2019年下半年改造后的東輸沙渠，根據現有的水文資料情況，主要驗證內容為沿程水面線。使用改造實施方案中的設計規(guī)劃圖所規(guī)定的渠底高程、設計水位及典型斷面的構造等數據為驗證參數，最終得出驗證水面線和設計水面線的數值，見圖3。由圖可以看出，經過多次的率定驗證后，二維水動力模型輸出的模擬計算水位與河道設計水位相比，水位值存在較小正負偏差，水位的絕對誤差為0～0.07 m，均方根誤差(RMSE)為0.054，水面線的驗證模擬計算結果與河道設計水位過程線基本吻合，整體誤差較小，模型擬合度較好，表明模型的構建及參數設定合理。

在各典型斷面的站點附近上、下游200 m處為起點、終點截取水位圖，如圖4所示，在關山東站(距閘口1 100 m)水位為38.2～38.4 m，張廣站(距閘口5 846 m)水位為37.4～37.6 m，獅子宋站(距閘口8 000 m)水位為37.2～37.4 m，后張站(距閘口10 000 m)、固堆王站(距閘口10 900 m)水位為37.0～37.2 m，王小樓站(距閘口13 290 m)水位為36.6～36.8 m，模擬計算結果合理，符合水流運動規(guī)律，該模型可以用于下一階段對泥沙數值模擬及預測。

2.4.2 泥沙模塊驗證

截取渠道的6個典型橫斷面模擬淤積1 a后的泥沙淤積厚度進行對比分析(見圖5)。由圖可以看出，最大泥沙淤積厚度為2.0 m，紅色區(qū)域淤積情況最為嚴重，主要淤積分布在自閘口至距閘口500 m處渠道中線處，渠道中段的淤積厚度在1.2～1.5 m之間，渠尾的平均淤積厚度為1.0 m。同一斷面的淤積先在地形高程較低的位置形成，在兩側原渠底的部分淤積厚度為0.52～0.86 m，中間新建挖深后子槽渠底的淤積厚度為1.12～2.17 m。從各斷面渠底高程、實測高程與計算高程的對比可以看出，淤積量計算高程線基本擬合于實際測定高程線，演變態(tài)勢大致相同，表明耦合泥沙模塊后的數學模型計算的淤積量與橫截面上的實際淤積量基本吻合。

所建立的泥沙數學模型基本上可以反映東輸沙渠的沖淤演變規(guī)律，整體趨勢符合較好，沖淤量計算結果與實測結果相近，可應用于研究不同流量下東輸沙渠的沖淤過程中水流沿程變化及橫向、縱向泥沙淤積分布及淤積量的變化。

3 不同方案沖淤數值模擬對比

根據東輸沙渠水位及流量關系，以水流條件作為自變量，分別選取設計流量為10、20、40、60 m3/s對應的水位，由實測資料及合理率定得出，具體計算工況條件見表2。

表2 模擬工況數據設置

以率定驗證后的水動力泥沙耦合模型作為工具，根據4種方案的流量水位作為起始條件運行模型，得到4種預測情景下東輸沙渠流速、水位和含沙量等模擬數據，觀察運行過程中水流條件沿程變化、河道主流變化情況、淤積體平面淤積形態(tài)及其變化過程以及橫斷面淤積特性。

3.1 水流特征變化

模擬東輸沙渠不同方案時的泥沙淤積過程，提取整條渠道流速、水位沿程變化曲線，結果如圖6所示。4個方案的渠道中段距閘口6 000～10 000 m處的平均流速依次為0.527、0.648、0.868、1.103 m/s，流速隨流量增大而增加。由于流速受過水斷面的影響，因此在同一流量時整體的渠道流速有所減小，距閘口5 000 m處的流速有小幅增加，原因是中游區(qū)域斷面的過水面積有所減小，但實際過流面積變化率不大，流速增幅僅為0.10～0.15 m/s。方案1的流速沿程變化幅度不大，方案2的流速有沿程逐漸減小的趨勢但波動不明顯，方案3、4的流速變化趨勢相似，均是上段流速相對較大，下段流速相對較小，流速變化最大的位置大致相同。

圖6 不同方案的淤積后流速沿程變化

4個方案的水體流向、水體流速變化趨勢基本保持一致。引水口附近無漩渦，整個河道上也沒有出現漩渦和回流，流態(tài)比較平穩(wěn)。圖7是距閘口5 000 m處的局部水力特性圖，清晰地展示在不同方案下渠道中段區(qū)域的水流速大小及流向。由于渠道較為平直單一，彎道處的曲率半徑較大，因此整體流速變化幅度不大，流向沒有變化。

不同方案的淤積后水面線沿程變化如圖8所示。由圖可見，除方案1外，渠道入口周圍的水流受河岸邊界的影響產生壅水，絕大部分水流集中到渠道中，隨著渠道底坡的沿程降低變化，河道水面高程逐漸降低，渠道內部水面高程均勻下降，水位高程降低速率大致相同，距閘口4 500 m處為轉折點，出現較大水位變化，在該段出現水位下降的原因主要是河道主流出現側向收縮，水位的下降速率增大。由于渠道引水流量的增加，河道水面線隨著流量的增加不斷抬升，渠道距閘口5 000～8 000 m處的渠道的中間段水位抬升幅度較大，其中方案4最為明顯，最大水位值出現在距閘口6 000 m處，隨著流量的增加，方案3、4的水位抬升幅度分別為0.81、0.68 m，渠道下游距出口2 000 m抬升幅度逐漸減小。由于方案1的水位受到泥沙淤積厚度的影響，因此變化規(guī)律與其他3個方案不同，水面線降低幅度較大，前后水位相差3.6 m。

圖8 不同方案的淤積后水面線沿程變化

3.2 泥沙淤積特征變化

3.2.1 泥沙淤積縱向特征變化

泥沙淤積縱向形態(tài)的發(fā)展直接體現了渠道沿程被泥沙侵占的過程及縱向分布。不同方案提取的渠道內等距的縱向深泓點沿程變化如圖9所示。從圖中可以看出，在入口起至閘前500 m處是淤積最為嚴重的地區(qū)，方案1的深泓點的淤積高度最高達2 m，方案2的淤積高度為0.5～1 m，方案3、4淤積不明顯。距閘口1 000～12 200 m是渠道的穩(wěn)定沖淤范圍，方案1的泥沙的形態(tài)最為明顯，淤積高度從1.6 m到1.0 m逐層向后遞減，方案2的泥沙沿程淤積較為均勻，淤積厚度為0.1～0.2 m，在8 000～10 000 m有沖刷的現象。方案3、4的沿程淤積較少，渠道多存在沖刷現象，渠道內泥沙淤積厚度最大為0.1 m左右。流量增大后，渠道內的泥沙顯著減少。通過對比分析整個縱斷面的河床淤積變化，可以看出方案3、4的淤積較為均勻、穩(wěn)定。

圖9 不同方案的渠道縱向深泓點高程變化

3.2.2 泥沙淤積總量特征變化

泥沙隨水流從閘口進入渠道后，一部分隨水流向前流動，當淤積剪切應力小于臨界淤積剪切應力后開始出現沉積，一部分泥沙在渠底淤積，閘后至500 m處渠道前段的泥沙淤積量較大，之后渠道中后段含沙量逐漸減少，泥沙懸浮量較少，水流攜帶細顆粒泥沙沿河槽向前輸移，在渠道處呈點狀或片狀擴散的趨勢，待流速穩(wěn)定過后落淤至渠底，泥沙的淤積量也相應減少。

流量的加大使得這一過程的水流攜帶的泥沙大部分處于懸浮狀態(tài)，流量為60 m3/s時還可能對渠底形成沖刷，渠道內泥沙淤積較少，水流將攜帶混合大量泥沙，最后將泥沙送至沉砂池發(fā)生沉降。

4個方案的泥沙淤積量計算結果見表3。從表中數據可看出，最大、最小淤積高程都隨著與閘距離的增加而減小，與渠道渠底高程有關，隨著流量的增大相同河段的淤積量不斷減少，方案3、4之間的變化不明顯。淤積厚度中，除方案1渠前段最大，大量泥沙在渠首充分落淤，其余3個方案均是后段大于前段和中段。方案1的平均淤積厚度大于1 m，泥沙淤積嚴重，方案3、4的淤積厚度達到較為理想狀態(tài)。方案2的泥沙淤積量比方案1的大253 723 m3，方案3的泥沙淤積量比方案2的大35 167 m3，方案4的泥沙淤積量比方案3的大12 173 m3，方案3、4的累積淤積變化量最小，大流量下的泥沙淤積在渠道內的部分較少。

表3 不同方案的泥沙淤積計算結果

4 結論

1)本文中基于水力學模型MIKE 21構建了二維水沙模型，首先對水動力模塊進行率定和驗證，模擬計算水位與河道設計水位值的絕對誤差為0～0.07 m，均方根誤差為0.054，模擬計算結果與河道設計水位過程線基本吻合，模型的構建及參數設定合理；其次對泥沙模塊進行率定和驗證，淤積量計算高程線基本擬合于實際測定高程線，演變態(tài)勢大致相同，耦合泥沙模塊后的數學模型計算的淤積量與橫截面上的實際淤積量基本吻合，表明水流與泥沙模型具有較高的精度，適用于研究輸沙渠水沙運移特性。

2)根據水流特征變化可知，本文中4個方案的渠道內的水位、流速隨上游河道流量的增加而增加，隨渠底高程的降低而沿程下降。其中，渠道水流形態(tài)和流速方向基本一致，不同流量時流速變化最大的位置大致相同，流量為40、60 m3/s的變化趨勢最為相似。除方案1受到泥沙淤積量的影響，水位降幅較大，其他3種方案的水位高程降低速率大致相同。

3)從泥沙淤積的縱向形態(tài)發(fā)展變化來看，方案1的泥沙淤積厚度從1.6 m至1.0 m逐層向后遞減，方案2的泥沙沿程淤積較為均勻，淤積厚度為0.1～0.2 m，方案3、4的泥沙沿程淤積較少，渠道多存在沖刷現象，渠道內泥沙淤積厚度最大在0.1 m左右，表明流量增大后渠道內的泥沙顯著減少。

4)通過對泥沙淤積總量的匯總及分析對比各方案的結果，方案3、4的累積淤積量變化不大，可以認為基本接近沖淤平衡狀態(tài)，可以達到沖沙的目的，但方案4存在明顯的沖刷渠底的現象，因此可選用流量為40 m3/s作為渠道沖沙的優(yōu)選方案。