三峽庫區分汊河段細沙起動試驗研究及數值模擬

2019-09-12 08:09:58李文杰楊勝發

重慶交通大學學報(自然科學版) 2019年9期

肖毅，李文杰，楊勝發

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

0 引言

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

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

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

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

1 試驗

1.1 數據采集

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

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

1.2 細沙起動水槽試驗

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

圖3 細沙起動試驗水槽布置Fig. 3 Layout of flume for fine sediment incipient experiment

1.3 三峽庫區細沙起動流速公式分析

1.3.1 細沙起動流速實測值與理論計算值對比

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

(1)

(2)

(3)

(4)

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

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

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

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

由表1可見：

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

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

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

1.3.2 改進的細沙起動流速公式

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

(5)

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

(6)

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

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

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

2 平面二維水沙模型

2.1 現有平面二維水沙模型

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

(7)

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

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

2.2 改進細沙模塊

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

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

(8)

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

3 模型驗證

3.1 網格劃分

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

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

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

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

3.2 水流模塊驗證與結果分析

3.2.1 水位

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

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

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

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

3.2.2 流速

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

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

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

由圖11可知，江心洲與二次流的共同作用使得分汊河道的流速分布變為不對稱，汊道進口處的右汊流速為左汊流速的3倍。從斷面S202與S203流速的實測值與計算值分布的對比可以看出，該平面二維模型的水流模塊能模擬彎曲分汊型河道水動力過程。

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

3.3 細沙模塊控制方程驗證分析

3.3.1 細沙沖淤量驗證

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

由表2可以看出：

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

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

3.3.2 河道地形驗證

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

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