任 杰，劉豪杰

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

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

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

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

1.1 流體控制方程

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

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

(1)

(2)

(3)

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

(4)

vt=k/ω

(5)

式中：式中δij為克羅里克符號(hào)。

湍流動(dòng)能k方程為：

(6)

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

(7)

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

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

(8)

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

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

河床多孔介質(zhì)中的溫度傳播采用對(duì)流傳熱方程，其計(jì)算式為[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為熱(匯)源項(xiàng)。

2 計(jì)算模型

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

2.1 邊界條件

水力邊界：河流水位ac邊為定水頭邊界，bd邊為自由出流，ab為對(duì)稱邊界，cd邊界為不可移動(dòng)的墻邊界；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ù)從文獻(xiàn)[2，14-15]中進(jìn)行選取，模型參數(shù)見表1。

表1 模型參數(shù)

2.3 全局靈敏度分析

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

3 結(jié)果分析

3.1 全局靈敏度分析

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

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

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

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

3.2 溫度場(chǎng)分析

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

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

4 結(jié) 論

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

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

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