大壩下游河床潛流帶溫度場的影響因素研究

任 杰，劉豪杰

(西安理工大學(xué) 水利水電學(xué)院， 陜西 西安 710048)

河流中大壩建成后，壩前后的河流水位及流速發(fā)生改變，壩前水位抬升引起水溫呈現(xiàn)垂直分層現(xiàn)象，分層水溫通常為上層水溫高，下層水溫低[1]。水庫下泄水溫與天然河流水溫存在顯著差異，對河床潛流帶內(nèi)水熱交換有一定影響。水溫也是影響河流及其鄰近地區(qū)(如河岸帶、洪泛區(qū)和潛流帶)的水生生態(tài)系統(tǒng)重要的因素之一[2]。因此，越來越多的研究學(xué)者認識到地表水-地下水相互作用的重要性。陳孝兵等[3]構(gòu)建了循環(huán)式水槽裝置，基于擴散模型，研究了不同河床地形的潛流交換與地表水動力及河床滲透特性之間的關(guān)系。劉東升等[4]以新安江大壩下游為研究對象，對比分析了該水庫下游河岸帶冬夏季潛流交換特征及溫度場分布規(guī)律。陳斌等[5]揭示了滲流儀的不同裝置結(jié)構(gòu)對潛流交換的影響。李英玉等[6]研究了新安江大壩下游河岸帶溫度場時空分布特征，并基于一維瞬態(tài)解析模型對地下水流速進行了求解。國外專家，Norman等[2]構(gòu)建了循環(huán)式水槽裝置，研究了3種河床地形、3種水流流量和2種滲透特性情況下河床潛流帶熱流傳輸規(guī)律。Mamer等[7]通過構(gòu)建一個循環(huán)式水槽裝置，并利用一維瞬態(tài)解析模型對地下水和地表水的相互作用進行了量化，結(jié)果表明一維瞬態(tài)解析模型方法在精細尺度上能夠準確定位和量化潛流交換通量。Sawyer等[8]使用了實驗室水槽試驗和數(shù)值模擬兩種方法對潛流和熱交換進行了量化。上述大都是針對單一因素影響下河床潛流帶的規(guī)律研究，并未針對多種因素共同影響下河床潛流帶的溫度場影響規(guī)律變化問題進行研究分析。

本文構(gòu)建一個概化水槽模型，并以水溫作為一種示蹤劑，基于雷諾平均方程(RANS)和k-ω湍流模型，利用商業(yè)有限元軟件CFD-Fluent和COMSOL Multiphysics進行求解，并對模型參數(shù)進行Morris全局靈敏度分析，研究大壩下游河床潛流帶內(nèi)溫度場在水流流速、河道水位、沙坡高度和河床滲透系數(shù)等單因素和多因素影響下的變化規(guī)律。

1 數(shù)學(xué)模型

1.1 流體控制方程

Patel等[9]和Yoon等[10]通過二維水槽試驗與數(shù)值模擬對比了相同邊界條件下沙波上方流速分布、漩渦大小、阻力系數(shù)等特征因素，結(jié)果顯示k-ω湍流模型對刻畫沙波上方紊流結(jié)構(gòu)具有較好的適用性。

對于不可壓縮流體，雷諾平均方程為[11-12]：

(1)

(2)

(3)

雷諾應(yīng)力與湍流動能k比率耗散系數(shù)ω有關(guān)，其計算式為：

(4)

vt=k/ω

(5)

式中：式中δij為克羅里克符號。

湍流動能k方程為：

(6)

耗散系數(shù)ω方程為：

(7)

1.2 多孔介質(zhì)滲流控制方程

河床飽和多孔介質(zhì)中滲流連續(xù)控制方程為[13]：

(8)

式中：ρf為水體密度，kg/m3；v為滲流速度矢量，m/s。

1.3 多孔介質(zhì)溫度場控制方程

河床多孔介質(zhì)中的溫度傳播采用對流傳熱方程，其計算式為[14]：

(9)

式中：εp為多孔介質(zhì)孔隙率；cw為水體導(dǎo)熱系數(shù)，J/(kg·℃)；ρs為多孔介質(zhì)密度，kg/m3；cs為多孔介質(zhì)導(dǎo)熱系數(shù)，J/(kg·℃)；Tw為水溫，℃；λij為多孔介質(zhì)導(dǎo)熱系數(shù)，W/(m·℃)；QT為熱(匯)源項。

2 計算模型

參考Norman等[2]的室內(nèi)水槽物理模型，河床為5.0 m(長)×0.7 m(高)，床面有9個沙坡，沙坡有效范圍為0～4.5 m。單個沙坡0.5 m(長)×0.05 m(高)。地表水水位為0.1 m。地表水-地下水數(shù)值模型的概念模型如圖1所示。以2.0 m～3.5 m之間的區(qū)域作為研究對象，研究區(qū)域內(nèi)河床監(jiān)測點布設(shè)5行×7列，如圖2所示。地表水模型采用CFD-FLUENT建模并求解，共劃分1 058 125個網(wǎng)格單元；地下水模型采用COMSOL Multiphysics建模并求解，共劃分97 158個網(wǎng)格單元。

2.1 邊界條件

水力邊界：河流水位ac邊為定水頭邊界，bd邊為自由出流，ab為對稱邊界，cd邊界為不可移動的墻邊界；cd邊界為壓力邊界，ce、ef、df邊界均為不透水邊界。

溫度邊界：多孔介質(zhì)模型中，ce、ef、df邊界均為絕熱邊界，cd邊界溫度與水體溫度相同。

初始條件：初始流速為0.14 m/s，水深為0.1 m。水流水溫為10℃，河床底質(zhì)滲透系數(shù)k和初始溫度分為0.033 m/s和20℃。

2.2 模型參數(shù)

模型所需的水力及熱力學(xué)特性參數(shù)從文獻[2，14-15]中進行選取，模型參數(shù)見表1。

表1 模型參數(shù)

2.3 全局靈敏度分析

考慮到所建立的模型參數(shù)靈敏度對輸出結(jié)果的影響，本文選取Morris方法[15]對河床潛流帶溫度場的影響因素進行全局靈敏度分析，具體計算步驟參考文獻[15]。本文選取坡高h，流速v，沉積物滲透系數(shù)k和水深H等4個因素，進行靈敏度分析。數(shù)值計算時，各因素按照20%增加進行求解，模擬時間為24 h。

3 結(jié)果分析

3.1 全局靈敏度分析

根據(jù)各因素排列方式得到10種因素組合，并進行數(shù)值求解，結(jié)果如表2所示。由表2可知：

(1) 針對單一因素分析，流速v和水深H是影響河床潛流帶溫度場的主要因素，但從流速v和水深H對河床溫度場作用的幅度來看，流速v對河床潛流帶溫度場的影響是負相關(guān)的，而水深H對河床溫度場的影響是正相關(guān)的。

(2) 針對多因素分析，組別7、8、10河床潛流帶溫度場的影響較大，均呈現(xiàn)負相關(guān)。組別9影響最小。

(3) 考慮這10組因素組合，單個因素變化或因素組合變化均會影響河床潛流帶的溫度場。

3.2 溫度場分析

由圖3可知，研究區(qū)域內(nèi)的溫度場等值線在模擬結(jié)束時均會在沙坡下面會出現(xiàn)一個半圓型的低溫區(qū)域。結(jié)合圖4和圖5可知，迎水面壓力大于背水面，導(dǎo)致地下水從壓力大處向壓力小處運動，使沙坡下出現(xiàn)一個半圓的低溫區(qū)域。由圖3(a)、圖3(b)可知，當流速v增大時，研究區(qū)域內(nèi)溫度場等值線均有所向下降的趨勢。圖3(c)是流速v、水深H、河床底質(zhì)滲透系數(shù)k和沙坡高度h等4個因素均發(fā)生變化情況下的河床研究區(qū)域內(nèi)的溫度場變化情況，相比額定溫度場(見圖3(a))和流速v發(fā)生變化(見圖3(b))情況而言，多因素共同變化對河床影響的更為明顯，說明考慮多因素共同影響更能反映河床真實狀態(tài)下的溫度場變化規(guī)律。

表2 全局靈敏度分析因素組合及溫度變幅

4 結(jié) 論

本文基于雷諾平均方程(N-S方程)、k-ω湍流模型，并利用CFD-Fluent及COMSOL Multiphysics軟件，構(gòu)建地下水-地表水耦合模型；通過Morris全局靈敏度分析方法探討單因素影響及多因素影響下河床潛流帶內(nèi)溫度場變化規(guī)律，結(jié)論如下：流速v、水深H、河床底質(zhì)滲透系數(shù)k和沙坡高度h等。

(1) 針對單一因素，河床潛流帶溫度場受水流流速v、河道水位H的影響最為明顯，其次為沙坡高度h和河床滲透系數(shù)k；針對組合因素，坡高h和滲透系數(shù)k組合、流速v、水位H和坡高h組合以及流速v、水位H、坡高h和滲透系數(shù)k組合等3個組合對河床潛流帶溫度場影響較大。

(2) 河床潛流帶的溫度場受河床表面壓力分布影響，且在沙坡下會出現(xiàn)一個半圓的低溫區(qū)域；在多因素共同影響下，河床潛流帶溫度變化較為明顯。