肖 毅，李文杰，楊勝發(fā)

(重慶交通大學(xué) 國家內(nèi)河航道整治工程中心，重慶 400074)

0 引 言

三峽水庫蓄水以來，由于原有水動力與來沙補(bǔ)給條件發(fā)生了改變，致使庫區(qū)部分分汊河段產(chǎn)生大量細(xì)沙累積性淤積，一方面，引起主支汊易位，從而導(dǎo)致航道季節(jié)性改道[1]；另一方面，淤沙性淺灘增多引起航道尺度縮減，部分河段航道尺度甚至減小為運(yùn)行初期的30%[2]。而所關(guān)注庫區(qū)河段細(xì)顆粒淤積物粒徑中值僅為0.008 mm，具有明顯流動特性[3]，通常，此類細(xì)顆粒泥沙淤積多出現(xiàn)于水庫壩前及河口海岸，山區(qū)河道型水庫庫區(qū)河段較為少見。目前，三峽庫區(qū)的細(xì)沙淤積已打破部分分汊河段原有穩(wěn)定性，開始對三峽庫區(qū)局部河段消落期、汛期的航運(yùn)能力及長江“黃金水道”的全面建設(shè)產(chǎn)生不利影響。隨著未來三峽水庫上游水利樞紐的開發(fā)及水土保持工作的深入開展，入庫泥沙勢必進(jìn)一步減少且變細(xì)，庫區(qū)淤沙性分汊河段細(xì)沙淤積、分布等問題將更加嚴(yán)重地影響庫區(qū)航運(yùn)的暢通性。

三峽庫區(qū)內(nèi)，常年回水區(qū)的忠縣河段距三峽大壩上游約350 km，平面形態(tài)為彎曲分汊型(圖1)，其整體長度為24 km，河寬750 ～ 1 000 m。三峽水庫蓄水以來，常年回水區(qū)細(xì)沙已淤積了1.4億t，遠(yuǎn)遠(yuǎn)超過了論證階段的研究結(jié)果；其中，忠縣河段的左汊淤積較為嚴(yán)重，最大淤積厚度達(dá)40 m，汛期上行船只不得不改走右汊航道。

圖1 三峽庫區(qū)忠縣河段實(shí)際淤積分布Fig. 1 Distribution of actual sedimentation in Zhongxian reaches of Three Gorges Reservoir

筆者采用三峽庫區(qū)忠縣黃花城典型分汊河段的現(xiàn)場細(xì)沙，開展細(xì)沙樣本起動特性水槽試驗(yàn)，以研究細(xì)沙起動特性；提出了適用于三峽庫區(qū)的細(xì)沙起動流速公式，并將公式應(yīng)用于已構(gòu)建的平面二維水沙數(shù)值模型，使得改進(jìn)模型能準(zhǔn)確模擬三峽庫區(qū)淤沙性分汊河段細(xì)沙輸移。同時，筆者采用改進(jìn)模型對黃花城分汊河段細(xì)沙沖淤過程進(jìn)行了模擬計算，并與原模型計算值及實(shí)測值進(jìn)行了對比分析。結(jié)果表明：改進(jìn)模型能較為真實(shí)地反映三峽庫區(qū)細(xì)沙的運(yùn)動過程。

1 試 驗(yàn)

1.1 數(shù)據(jù)采集

2017年7月18日—8月2日，筆者所在研究團(tuán)隊(duì)在三峽庫區(qū)忠縣河段選取了8個斷面(圖2)，采集細(xì)沙樣本[4]，同時布置27條垂線，共計193個測點(diǎn)；用船上GPS定位，采用ADCP測量各斷面流速U以及水深h。

圖2 忠縣河段現(xiàn)場測量布置Fig. 2 Layout of field measurement of Zhongxian reaches

1.2 細(xì)沙起動水槽試驗(yàn)

筆者針對忠縣黃花城河段細(xì)沙樣本，開展了不同干容重γs下臨界起動條件的水槽試驗(yàn)研究。試驗(yàn)水槽的長×寬×高=6 m × 0.25 m × 0.5 m，可調(diào)節(jié)比降J=0.005～0.050；為便于觀測細(xì)沙運(yùn)動狀態(tài)，水槽兩側(cè)鑲嵌玻璃板。將黃花城河段細(xì)沙樣本(中值粒徑D50=0.008 mm)鋪設(shè)于水槽中，試驗(yàn)段長約1 m，鋪設(shè)后，試驗(yàn)段與前方水流調(diào)整段和尾門控制端的水槽床面保持齊平；在水槽1#、2#、3# 斷面處分別布置電磁流量計、超聲水位儀及流速計各1個(圖3)，以測量斷面流量Q、水位H及流速U。筆者共進(jìn)行了6組試驗(yàn)，表1為3組的試驗(yàn)結(jié)果。

圖3 細(xì)沙起動試驗(yàn)水槽布置Fig. 3 Layout of flume for fine sediment incipient experiment

1.3 三峽庫區(qū)細(xì)沙起動流速公式分析

1.3.1 細(xì)沙起動流速實(shí)測值與理論計算值對比

目前，常用的細(xì)沙起動流速Uc理論計算公式包括唐存本公式(1)[5]、竇國仁公式(2)[6]、張瑞瑾公式(3)[7]及竇國仁公式(4)[8]：

(1)

(2)

(3)

(4)

式中：h為水深，m；d為細(xì)沙代表粒徑，m；Ha為水大氣壓力的水柱高，m，通常Ha=10 m；γs、γ′s、γ′*s分別為細(xì)沙的干容重、實(shí)際干容重與穩(wěn)定后干容重，N/m3；γw為水的容重，N/m3；ρw為水密度，kg/m3；ks為粗糙高度，m；m、C、δ、D*、ε0分別為細(xì)沙起動參數(shù)〔天然河流下，m=6；水槽試驗(yàn)下，m=4.7×(h/D50)0.06〕。

根據(jù)水槽試驗(yàn)測得的參數(shù)，按照式(1)～式(4)得出的細(xì)沙起動流速理論計算值Uc及水槽試驗(yàn)值Uec，見表1。計算值與試驗(yàn)值的對比見圖4。

圖4 細(xì)沙啟動流速試驗(yàn)值與理論公式計算值的對比Fig. 4 Comparison between experimental and calculated values of incipient velocity of fine sediment

表1 細(xì)沙起動流速試驗(yàn)值與理論計算值對比Table 1 Comparison between experimental and calculated values of incipient velocity of fine sediment

由表1可見：

1)式(2)、式(3)得到的理論計算值Uc為試驗(yàn)值Umc的1.5～2.5倍。

2)由于考慮了水壓力與分子力，式(1)及式(4)得到的理論計算值均比試驗(yàn)值小，相對誤差在-40%～40%范圍內(nèi)。

由圖4可知，式(1)的理論計算值與6組樣本細(xì)沙試驗(yàn)起動流速趨勢基本一致。

1.3.2 改進(jìn)的細(xì)沙起動流速公式

由于三峽水庫常年回水區(qū)河段最大水深可達(dá)100 m，所以水壓力對細(xì)沙起動的影響極為顯著。因此，筆者針對與三峽庫區(qū)細(xì)沙起動試驗(yàn)相近的式(1)，提出了考慮大水深水壓力影響的細(xì)沙黏性作用力N的計算式(5)：

(5)

從而得到改進(jìn)細(xì)沙起動流速公式(6)：

(6)

式中：C、C1、k1、k2分別為細(xì)沙起動參數(shù)，由表1水槽試驗(yàn)測量數(shù)據(jù)得到：C=4.324×10-4g/cm、C1=2.252×10-9m2/s2、k1=k2=5；其他符號同前。

圖5為根據(jù)筆者提出的改進(jìn)細(xì)沙起動流速公式(6)、唐存本公式(1)[5]與竇國仁公式(4)[8]所計算的起動流速Uc，與2012年忠縣河段現(xiàn)場實(shí)測值的對比。由圖可見，改進(jìn)細(xì)沙起動流速公式(6)計算值的相關(guān)性R2提高至0.72，表明：改進(jìn)細(xì)沙起動流速公式能更好反映三峽庫區(qū)的細(xì)沙起動過程。

圖5 細(xì)沙起動流速實(shí)測值與理論公式計算值的對比Fig. 5 Comparison between field measured and calculated values of incipient velocity of fine sediment

2 平面二維水沙模型

2.1 現(xiàn)有平面二維水沙模型

現(xiàn)有的平面二維水沙模型包含了考慮彎道二次流的水流模塊與不平衡輸沙模塊[8]。輸沙模塊中的細(xì)沙模塊控制方程如式(7)：

(7)

式中：εξ、εη分別為計算方向的細(xì)沙擴(kuò)散系數(shù)；Sk為第k組細(xì)沙含沙量，kg/m3；Fk為床面細(xì)沙交換項(xiàng)，kgf/m3；Sok為河道細(xì)沙補(bǔ)給源項(xiàng)，由于研究河段并未細(xì)沙源項(xiàng)補(bǔ)給，因此不予考慮。

采用ADI有限差分法對該平面二維水沙模型的水流模塊控制方程進(jìn)行離散，除了連續(xù)方程中水位對時間偏導(dǎo)數(shù)項(xiàng)采用向前差分、以及動量方程中對流項(xiàng)采用一階迎風(fēng)和中心差分格式組合外(QUICK格式)，其余各項(xiàng)采用中心差分。細(xì)沙模塊中控制方程的離散采用類似水流方程的方法，只是在對流項(xiàng)中沒有采用QUICK格式，而是采用的一階迎風(fēng)[9]。

2.2 改進(jìn)細(xì)沙模塊

根據(jù)三峽水庫庫區(qū)忠縣河段的現(xiàn)場細(xì)沙運(yùn)動測量分析結(jié)果[2-3]，細(xì)沙運(yùn)動可分為沉積、沖刷以及輸運(yùn)3個階段。因此，細(xì)沙模塊控制方程(7)中床面細(xì)沙交換項(xiàng)Fk可寫為：

Fk=αkωk(β1S*-β2S)

(8)

式中：β1、β2分別為細(xì)沙運(yùn)動參數(shù)；Uc，k為計算第k組細(xì)沙起動流速，m/s，按式(6)計算；Uf為細(xì)沙臨界沉降速度，m/s，采用2012年忠縣河段現(xiàn)場測量值0.3 m/s[3]；αk第k組細(xì)沙的平衡非飽和系數(shù)；ωk為第k組計算細(xì)沙的沉速，m/s；S*為采用張瑞瑾輸沙公式第k組細(xì)沙的含沙量與輸沙率，kgf/m3。

3 模型驗(yàn)證

3.1 網(wǎng)格劃分

三峽庫區(qū)忠縣河段平面形態(tài)為彎曲型，計算初始地形采用2011年11月的實(shí)測資料，比尺為1∶2 000。計算的進(jìn)口條件采用距離忠縣段100 km的清溪場水文站流量水位觀測數(shù)據(jù)，出口條件采用忠縣河段下游60 km的萬縣水文站資料[7]。

為了確保邊界條件的準(zhǔn)確性，建模的計算區(qū)域長度為184 km，計算網(wǎng)格為1 600 × 21個，其中忠縣河段的網(wǎng)格為200 × 21 (圖6)，其上游網(wǎng)格為1 000 × 21，下游網(wǎng)格為400 × 21(圖7)。計算網(wǎng)格除了靠近河岸的少數(shù)網(wǎng)格外，在ξ與η方向正交性基本保持在+88 ～-92°之間，網(wǎng)格步長在ξ方向?yàn)?0～110 m、η方向?yàn)?5～40 m。

圖6 研究河段資料Fig. 6 Field data of reaches studied

圖7 研究河段上、下游計算網(wǎng)格示意Fig. 7 Grating scheme for upper and lower reaches studied

3.2 水流模塊驗(yàn)證與結(jié)果分析

3.2.1 水 位

選用2012年7月14日、24日及8月11日典型洪水工況的水流模塊作為驗(yàn)證資料。進(jìn)口邊界條件：清溪場的流量分別為34 000、63 000、232 000 m3/s；對應(yīng)出口邊界條件：萬縣站水位分別為157.78、159.81、157.28 m。模擬時間步長Δt=6 s。

圖8為63 000 m3/s流量下實(shí)測與計算的沿程水位分布。從圖中可知，計算值基本與實(shí)測值吻合，沿程4個水文站的計算水位誤差控制在0.08 m以內(nèi)；進(jìn)口附近的水位誤差相對較大，主要是由于進(jìn)口計算的不穩(wěn)定性引起的。

圖8 沿程水位H的實(shí)測值與理論計算值(63 000 m3/s)Fig. 8 Field measured and calculated water level H all along (63 000m3/s)

圖9為3個典型流量下沿程測站的水位計算誤差分布。從圖中可看出，隨著河道彎曲度的增加，二次流的影響不容忽視。由于平面二維水流模型在模擬二次流方面存在局限，因此，計算相對誤差與河道彎曲度成正比，最大相對誤差出現(xiàn)于63 000 m3/s下的忠縣河段。

3.2.2 流 速

圖10為23 200 m3/s流量下斷面S206的流速實(shí)測值與計算值分布。由圖可見，計算誤差小于0.05 m/s，最大誤差出現(xiàn)在彎道頂部；由于彎道二次流的影響，誤差隨著彎曲度的降低而減小。

圖9 沿程水位H的實(shí)測值與理論計算值誤差Fig. 9 Error between field measured and calculated water level H all along

圖10 斷面S206流速實(shí)測值與理論計算值分布Fig. 10 Distribution of field measured and calculated velocity at S206

由圖11可知，江心洲與二次流的共同作用使得分汊河道的流速分布變?yōu)椴粚ΨQ，汊道進(jìn)口處的右汊流速為左汊流速的3倍。從斷面S202與S203流速的實(shí)測值與計算值分布的對比可以看出，該平面二維模型的水流模塊能模擬彎曲分汊型河道水動力過程。

圖11 斷面S202、S203流速實(shí)測值與計算值分布Fig. 11 Distribution of field measured and calculated velocity at S202 & S203

3.3 細(xì)沙模塊控制方程驗(yàn)證分析

3.3.1 細(xì)沙沖淤量驗(yàn)證

驗(yàn)證資料采用2012年的水沙實(shí)測數(shù)據(jù)，見表2。進(jìn)出口邊界條件采用清溪場和萬縣站的日流量和水位觀測數(shù)據(jù)；考慮該河段的細(xì)沙中值粒徑D50=0.008 mm，模擬過程中采用均勻沙計算，進(jìn)口采用清溪場站的細(xì)沙含沙量[10]；模擬計算步長Δt=12 s。分別用改進(jìn)前后的細(xì)沙模塊對該河段進(jìn)行沖淤變化計算，計算結(jié)果見表2。

由表2可以看出：

1)由原細(xì)沙模塊計算的細(xì)沙淤積量僅為實(shí)測值的30%，而由改進(jìn)細(xì)沙模塊計算的細(xì)沙淤積量為實(shí)測值的84%，與實(shí)測值的誤差較改進(jìn)前降低了54%，說明改進(jìn)細(xì)沙模塊能較為真實(shí)地反映三峽庫區(qū)河段細(xì)沙輸移過程。

2)河道的平面形態(tài)與計算精度有一定的關(guān)系，彎曲段(S206—S203)較相對順直段(S203—S202、忠縣—斷面S206)的計算誤差有一定的增加。

3.3.2 河道地形驗(yàn)證

圖12為研究河段實(shí)測地形及改進(jìn)前后細(xì)沙模塊計算的河道平面地形。從圖12(b)可以看出：原細(xì)沙模塊模擬的沖淤分布與實(shí)測值相差較大，左汊的地形變化等高線明顯小于實(shí)測地形值，主要原因是原細(xì)沙模塊控制方程(7)未考慮三峽庫區(qū)細(xì)沙輸移特性；而改進(jìn)后細(xì)沙模塊模擬的地形變化〔圖12(c)〕則與實(shí)測值〔圖12(a)〕基本一致：淤積部位主要在忠縣黃花城的左汊，而最大厚度由原模塊的1.2 m提高至3.0 m，相當(dāng)于2012年實(shí)測淤積厚度的80%。

圖13為該河段斷面地形的實(shí)測值與計算值。由圖可見：改進(jìn)后細(xì)沙模塊計算值與實(shí)測值的變化趨勢基本一致，而原細(xì)沙模塊的計算結(jié)果遠(yuǎn)小于實(shí)測值的變化。在彎道的進(jìn)口處，由于左邊界的深潭引起的低流速帶，該處的細(xì)沙淤積較為顯著。汊道出口的左汊和右汊的斷面地形變化表明：淤積主要發(fā)生在左汊，其淤積厚度在S204增加最大，約2 m；而此時的右汊斷面S203基本沒有淤積。從沖淤變化整個過程分析，改進(jìn)細(xì)沙模塊能夠模擬三峽水庫淤沙性分汊河段的細(xì)沙淤積及河勢變化過程。

圖12 研究河段實(shí)測的、原細(xì)沙模塊模擬的及改進(jìn)細(xì)沙模塊模擬的平面地形Fig. 12 Field measured and simulated topographies with original and improved fine sediment models studied

圖13 研究河段S203、S204斷面實(shí)測的及理論計算的地形對比Fig. 13 Comparison of field measured and calculated bed level at S203 and S204

4 結(jié) 論

1)選取三峽庫區(qū)典型淤沙分汊型河段開展了現(xiàn)場實(shí)測和水槽試驗(yàn)；通過對觀測資料的分析，從內(nèi)在動力機(jī)制層面確定了適應(yīng)于三峽庫區(qū)細(xì)沙的起動流速公式與細(xì)沙輸移控制方程。

2)改進(jìn)了已有的平面二維水沙數(shù)值模型。

3)將改進(jìn)的平面二維水沙數(shù)值模型應(yīng)用于三峽庫區(qū)典型淤沙分汊河段(忠縣黃花城)的細(xì)沙輸移過程模擬。對該河段沖淤總量、平面沖淤分布以及斷面地形變化的分析表明，改進(jìn)后的模型能夠較為真實(shí)、準(zhǔn)確地反映該河段的細(xì)沙淤積過程；研究結(jié)論能為三峽庫區(qū)淤沙河段演變趨勢預(yù)測與相關(guān)航道整治工程提供技術(shù)支撐。