位山灌區(qū)輸沙渠泥沙運移特性研究

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

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

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

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

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

1 研究區(qū)概況

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

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

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

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

2.1 水動力模型理論

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

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

連續(xù)性方程

(1)

X方向動量方程

(2)

Y方向動量方程

(3)

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

2.2 泥沙模型理論

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

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

(4)

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

D=ωsρbPD，

(5)

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

(6)

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

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

(7)

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

2.3 模型的建立

2.3.1 模型的范圍與網(wǎng)格

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

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

2.3.2 邊界條件及參數(shù)設(shè)置

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

表1 水動力模型邊界條件與模型參數(shù)設(shè)置

2.4 模型驗證

模型驗證的目的是檢驗所建模型與實際工況的相似度，對模型中相關(guān)參數(shù)的率定結(jié)果進行檢驗和修正。

2.4.1 水動力模塊驗證

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

在各典型斷面的站點附近上、下游200 m處為起點、終點截取水位圖，如圖4所示，在關(guān)山東站(距閘口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，模擬計算結(jié)果合理，符合水流運動規(guī)律，該模型可以用于下一階段對泥沙數(shù)值模擬及預(yù)測。

2.4.2 泥沙模塊驗證

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

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

3 不同方案沖淤數(shù)值模擬對比

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

表2 模擬工況數(shù)據(jù)設(shè)置

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

3.1 水流特征變化

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

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

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

不同方案的淤積后水面線沿程變化如圖8所示。由圖可見，除方案1外，渠道入口周圍的水流受河岸邊界的影響產(chǎn)生壅水，絕大部分水流集中到渠道中，隨著渠道底坡的沿程降低變化，河道水面高程逐漸降低，渠道內(nèi)部水面高程均勻下降，水位高程降低速率大致相同，距閘口4 500 m處為轉(zhuǎn)折點，出現(xiàn)較大水位變化，在該段出現(xiàn)水位下降的原因主要是河道主流出現(xiàn)側(cè)向收縮，水位的下降速率增大。由于渠道引水流量的增加，河道水面線隨著流量的增加不斷抬升，渠道距閘口5 000～8 000 m處的渠道的中間段水位抬升幅度較大，其中方案4最為明顯，最大水位值出現(xiàn)在距閘口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ā)展直接體現(xiàn)了渠道沿程被泥沙侵占的過程及縱向分布。不同方案提取的渠道內(nèi)等距的縱向深泓點沿程變化如圖9所示。從圖中可以看出，在入口起至閘前500 m處是淤積最為嚴(yán)重的地區(qū)，方案1的深泓點的淤積高度最高達(dá)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有沖刷的現(xiàn)象。方案3、4的沿程淤積較少，渠道多存在沖刷現(xiàn)象，渠道內(nèi)泥沙淤積厚度最大為0.1 m左右。流量增大后，渠道內(nèi)的泥沙顯著減少。通過對比分析整個縱斷面的河床淤積變化，可以看出方案3、4的淤積較為均勻、穩(wěn)定。

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

3.2.2 泥沙淤積總量特征變化

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

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

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

表3 不同方案的泥沙淤積計算結(jié)果

4 結(jié)論

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

2)根據(jù)水流特征變化可知，本文中4個方案的渠道內(nèi)的水位、流速隨上游河道流量的增加而增加，隨渠底高程的降低而沿程下降。其中，渠道水流形態(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的泥沙沿程淤積較少，渠道多存在沖刷現(xiàn)象，渠道內(nèi)泥沙淤積厚度最大在0.1 m左右，表明流量增大后渠道內(nèi)的泥沙顯著減少。

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