冰蓋下水流紊動(dòng)特性試驗(yàn)研究

陳 剛，董增川，王海軍，顧世祥，楊紅宣

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；2.云南省水利水電勘測(cè)設(shè)計(jì)研究院，云南 昆明 650021；3.昆明理工大學(xué)電力工程學(xué)院，云南 昆明 650500)

河冰是陸地冰凍圈的重要組成要素[1-2]，盡管覆蓋范圍相對(duì)較小，但對(duì)人類社會(huì)的影響卻最為顯著[3]。河冰生消演變會(huì)改變河流中物質(zhì)、能量的輸移規(guī)律[4-5]，在特定的水文氣象、河勢(shì)和水力條件下甚至造成冰凌災(zāi)害，危及沿岸人民的生命財(cái)產(chǎn)安全[6]。中國陸地冰凍圈占陸地國土面積的78%[1]，冰情影響范圍廣，冰凌災(zāi)害頻發(fā)[7]。冰蓋是河冰的重要存在形式，結(jié)冰蓋輸水是中國長距離調(diào)水、灌溉渠系工程實(shí)現(xiàn)冬季安全輸水的常用運(yùn)行方式[8]。冰蓋形成附加的無滑移邊界對(duì)河渠水流紊動(dòng)特性的影響是河渠冰凌災(zāi)害防治、輸水渠系冰期運(yùn)行、調(diào)水工程設(shè)計(jì)與運(yùn)行調(diào)度的重要理論基礎(chǔ)。因此，研究冰蓋下水流紊動(dòng)特性對(duì)自然災(zāi)害防治和供水安全保障具有重要意義。

根據(jù)河冰生消演變對(duì)水流表面的影響，水流流態(tài)可分為明流、完全冰封和部分冰封[4，9]。部分冰封是指水流表面部分被冰覆蓋(稱為冰蓋區(qū))，部分為自由水面(明流區(qū))，是前2種流態(tài)相互轉(zhuǎn)換的過渡流態(tài)，岸冰是其主要存在形式。在水流紊動(dòng)特性方面，相比明渠流已取得系統(tǒng)性研究成果，冰蓋下水流的研究相對(duì)較少，尤其是部分冰封。早期研究側(cè)重于揭示河冰生消的機(jī)理[10]，僅有少數(shù)研究者通過原型觀測(cè)分析冰蓋形成對(duì)流速分布和輸水能力的影響[11]。隨著先進(jìn)實(shí)驗(yàn)儀器的問世和室內(nèi)試驗(yàn)條件的改善，國內(nèi)外學(xué)者開始通過水槽試驗(yàn)研究冰蓋下水流的特性。Parthasarathy等[12]探究完全冰封水流縱向時(shí)均流速、雷諾應(yīng)力和紊動(dòng)能的分布規(guī)律；陳建國等[13]測(cè)量不同冰蓋組合下的水流結(jié)構(gòu)，探討冰蓋流的輸水能力；拾兵等[14]研究不同流速、水深、冰蓋糙率條件下完全冰封水流的縱向紊動(dòng)強(qiáng)度，發(fā)現(xiàn)冰蓋下水流縱向紊動(dòng)強(qiáng)度比明渠流大；Peters等[15]測(cè)量岸冰下水流冰蓋區(qū)和明流區(qū)的縱向流速分布，發(fā)現(xiàn)岸冰形成使流動(dòng)核心區(qū)向明流區(qū)收縮，導(dǎo)致明流區(qū)流速增加，冰蓋區(qū)流速減??；Nyantekyi-Kwakye等[16]采用PIV技術(shù)測(cè)量不同寬度岸冰下水流瞬時(shí)流場(chǎng)，用于檢驗(yàn)4種邊界切應(yīng)力計(jì)算方法的適用性；Wang等[17]測(cè)量岸冰下水流流速分布，用于校驗(yàn)冰蓋下水流橫向流速分布準(zhǔn)二維模型?？梢?，現(xiàn)有的研究多針對(duì)單一流態(tài)，對(duì)比研究3種流態(tài)紊動(dòng)特性的成果還較少。

本文在國內(nèi)外已有研究成果的基礎(chǔ)上，通過水槽試驗(yàn)測(cè)量3種流態(tài)的恒定均勻流場(chǎng)，從縱向時(shí)均流速、雷諾應(yīng)力和紊動(dòng)強(qiáng)度等方面分析冰蓋對(duì)水流紊動(dòng)特性的影響。

1 試驗(yàn)裝置和研究方法

1.1 試驗(yàn)裝置

試驗(yàn)在昆明理工大學(xué)水流與水工結(jié)構(gòu)重點(diǎn)實(shí)驗(yàn)室的自循環(huán)玻璃水槽試驗(yàn)系統(tǒng)中進(jìn)行，該系統(tǒng)由地下水池、泵站、量水堰、穩(wěn)流池、水槽和回水池組成，見圖1。水槽采用有機(jī)玻璃制作，長15 m，寬0.80 m，深0.35 m，縱坡坡降固定為0.5‰。水流通過泵站從地下水池提升至進(jìn)水渠，經(jīng)量水堰進(jìn)入穩(wěn)流池，溢流后進(jìn)入玻璃水槽，從水槽尾門出流后經(jīng)回水池進(jìn)入地下水池，形成水流循環(huán)。為保證水槽的水流流態(tài)穩(wěn)定，水槽進(jìn)口端為深1.4 m、長4.5 m、寬6.5 m的穩(wěn)流池，通過調(diào)節(jié)水槽末端的尾門來控制水槽中的水位，使水槽保持恒定均勻流。量水堰為無側(cè)收縮矩形薄壁堰，用于測(cè)量流量，堰上水頭由精度為0.1 mm的水位測(cè)針測(cè)量得到。流速采用挪威諾泰克公司(Nortek AS)生產(chǎn)的聲學(xué)多普勒點(diǎn)式三維測(cè)速儀Vectrion測(cè)量得到。

圖1 試驗(yàn)裝置示意Fig.1 Sketch of experimental set-up

1.2 試驗(yàn)材料

試驗(yàn)探究冰蓋下水流的紊動(dòng)特性，典型過流斷面見圖2(a)，其中，B為斷面寬度，b1為單側(cè)冰蓋的寬度，H為水深。以冰蓋邊緣為界劃分為冰蓋區(qū)和明流區(qū)，則斷面冰蓋覆蓋度定義為F=2b1/B，當(dāng)b1=0時(shí)，F(xiàn)=0；當(dāng)b1=B/2時(shí)，F(xiàn)=1，即明流、完全冰封可視為冰蓋覆蓋度分別為0、1的特例。試驗(yàn)以F為變量，共設(shè)計(jì)6種工況，F(xiàn)值分別為0、0.250、0.375、0.500、0.750、1.000，詳見表1。

在普通實(shí)驗(yàn)室中開展冰蓋下水流水力特性研究，通常采用模擬冰蓋替代真冰進(jìn)行水槽試驗(yàn)，利于形成穩(wěn)定流場(chǎng)。采用購置成本低、制作簡便的成型單塊泡沫板制作模擬冰蓋，使用火把T2型電熱絲泡沫切割機(jī)將泡沫板切割成長2.00 m，厚0.05 m，寬分別為0.10、0.15、0.20、0.30和0.40 m的條板。為確保水槽進(jìn)口水流流態(tài)平順穩(wěn)定，距進(jìn)口0.80 m段為槽內(nèi)穩(wěn)流段(圖1)，同時(shí)為避免尾門調(diào)節(jié)影響模擬冰蓋，距尾門0.20 m段不布置冰蓋，布置模擬冰蓋的渠段長14.0 m，每種工況各準(zhǔn)備14塊泡沫條板。

圖2 冰蓋下典型過流斷面和試驗(yàn)測(cè)量斷面布置示意Fig.2 Sketch of a typical under-ice section and measuring cross section

表1 試驗(yàn)工況

1.3 試驗(yàn)方案

試驗(yàn)時(shí)貼近水槽邊壁布置切割好的泡沫條板，并利用水槽橫梁上的旋轉(zhuǎn)螺絲桿進(jìn)行固定，使其貼近水槽邊壁，不能隨流波動(dòng)，但能上下浮動(dòng)。同時(shí)，由于真冰和泡沫板密度差異較大，通過旋轉(zhuǎn)螺絲桿施加垂向壓力，使泡沫板浸沒在水中的深度與同體積真冰的浸沒深度相等。

為確保試驗(yàn)在準(zhǔn)恒定均勻流條件下進(jìn)行，需要調(diào)整尾門使水槽內(nèi)水面近似平行于床面，為此布置5個(gè)水深測(cè)量斷面，分別距水槽進(jìn)口5.0、6.5、7.5、8.5和10.5 m，見圖2(b)。各斷面都安裝有水位測(cè)針，其零刻度線與水槽底板齊平，精度為±0.1 mm。在冰蓋寬度和流量固定時(shí)，通過調(diào)整尾門高度，使各斷面的水深相等。為確保水流充分發(fā)展和避免尾門的影響，流速測(cè)量斷面為水槽中段的C-C斷面，結(jié)合斷面的對(duì)稱性，僅測(cè)量半寬的斷面。自水槽右邊壁起，每間隔5 cm布置1條測(cè)線，共布置8條測(cè)線，依次記為Y05、Y10、Y15、Y20、Y25、Y30、Y35和Y40，每條垂線上均布置10個(gè)測(cè)點(diǎn)，見圖3。對(duì)于部分冰封工況，均布置1條測(cè)線與冰蓋邊緣垂線重合。采用標(biāo)準(zhǔn)Vectrion測(cè)量各測(cè)點(diǎn)的瞬時(shí)流速，采樣頻率設(shè)為25 Hz，采樣時(shí)間不少于60 s，即每個(gè)測(cè)點(diǎn)至少采樣1 500次。根據(jù)輸出數(shù)據(jù)的振幅數(shù)、信噪比和相關(guān)數(shù)，對(duì)采樣的瞬時(shí)流速數(shù)據(jù)點(diǎn)進(jìn)行質(zhì)量控制[18]，符合要求的用于計(jì)算時(shí)均流速與脈動(dòng)流速，進(jìn)而計(jì)算得到雷諾應(yīng)力和紊動(dòng)強(qiáng)度，用于揭示冰蓋對(duì)水流紊動(dòng)特性的影響。

圖3 C-C測(cè)量斷面上測(cè)點(diǎn)布置示意Fig.3 Sketch of measuring points on the C-C measuring section

2 結(jié)果與分析

2.1 縱向時(shí)均流速分布

時(shí)均流速計(jì)算公式如下：

(1)

式(1)表明時(shí)均流速受采樣時(shí)間的影響。以B20工況Y20測(cè)線上靠近冰蓋的測(cè)點(diǎn)為例，其縱向時(shí)均流速與采樣時(shí)間的關(guān)系見圖4，其中u′、v′、w′表示x、y、z方向的脈動(dòng)流速，為瞬時(shí)流速與時(shí)均流速之差。

圖4 采樣時(shí)間對(duì)時(shí)均流速的影響及其脈動(dòng)流速分布Fig.4 Influence of sampling time on time average velocity and its fluctuating velocity distribution

由圖4可見，當(dāng)采樣時(shí)間較短時(shí)，縱向時(shí)均流速波動(dòng)較大；當(dāng)采樣時(shí)間為30 s時(shí)，縱向時(shí)均流速趨于穩(wěn)定值，因此，試驗(yàn)確定采樣時(shí)間為60 s是合理的。同時(shí)，橫向、垂向時(shí)均流速均不等于0，這表明水流除縱向的主流動(dòng)外，存在疊加于主流之上的二次流；橫向時(shí)均流速小于0，表明該點(diǎn)的流量總體上是從明流區(qū)流向冰蓋區(qū)；垂向時(shí)均流速大于0，表明該點(diǎn)流量總體上是流向無滑移邊界的。在脈動(dòng)流速分布上，3個(gè)方向的瞬時(shí)流速在時(shí)均流速附近波動(dòng)，脈動(dòng)流速均滿足正態(tài)分布。

常用于描述冰蓋下縱向時(shí)均流速的分布律有雙對(duì)數(shù)律和雙冪律。雙對(duì)數(shù)律存在最大流速處不連續(xù)、當(dāng)量粗糙高度不易獲取、紊流核心區(qū)垂向流速分布不符合雙對(duì)數(shù)律[19]等問題，而雙冪律能描述冰蓋下水流垂向流速光滑連續(xù)的特性，與實(shí)測(cè)資料也較吻合[12，20]，故采用雙冪律描述連續(xù)冰蓋下水流縱向時(shí)均流速的垂向分布：

(2)

式中：k0為常數(shù)；mb、mi分別為與床面、冰蓋粗糙度有關(guān)的冪指數(shù)。特別地，當(dāng)mi→∞時(shí)，k0→umax，umax為最大流速。式(2)簡化為常用于描述明流縱向時(shí)均流速垂向分布的單冪律。采用平均百分比誤差和擬合度指標(biāo)評(píng)價(jià)雙冪律的模擬性能，其中平均百分比誤差(EMAP)為相對(duì)誤差絕對(duì)值的算術(shù)平均值，計(jì)算公式為

(3)

式中：i表示數(shù)據(jù)點(diǎn)數(shù)；n為數(shù)據(jù)量；Y為模型預(yù)報(bào)值；X為實(shí)測(cè)值。EMAP最小值為0，其值越小說明模型模擬值越接近參考數(shù)據(jù)。擬合度指標(biāo)是常用于判定非線性回歸方程擬合度的統(tǒng)計(jì)參數(shù)，有別于線性回歸的確定系數(shù)，故記為IFO，其計(jì)算公式為

(4)

式中：IFO最大值為1，其值越接近1，說明模型預(yù)報(bào)值與參考值的吻合程度越好。對(duì)于明流O00工況，各垂線的EMAP(u)均小于5.0%，IFO(u)均大于0.90，這表明明流工況縱向流速的垂向分布符合指數(shù)律，見圖5(a)。該工況的寬深比為4.17，受二次流影響最大流速點(diǎn)偏離自由水面[21]，特別是靠近邊壁的Y05測(cè)線。

圖5 O00工況和C40工況縱向時(shí)均流速垂向分布Fig.5 Vertical profiles of streamwise velocity for Run O00 and Run C40

完全冰封C40工況的實(shí)測(cè)流速與式(2)計(jì)算值進(jìn)行對(duì)比見圖5(b)。當(dāng)mb=9.69、mi=7.99時(shí)，各條垂線的EMAP(u)均小于5.0%，IFO(u)均大于0.90，這表明雙冪律可用于描述連續(xù)冰蓋下水流縱向時(shí)均流速的垂向分布。受冰蓋附加無滑移邊界的影響，各垂線最大流速點(diǎn)位置不同于明流工況，位于垂線中部區(qū)域。理論上，最大流速點(diǎn)的速度梯度為0，故將雙冪律對(duì)z求導(dǎo)，令導(dǎo)數(shù)為0得到最大流速點(diǎn)的相對(duì)水深為

(5)

部分冰蓋形成導(dǎo)致冰蓋區(qū)上邊界為無滑移邊界，明流區(qū)上邊界為自由水面，水流結(jié)構(gòu)較明流工況和完全冰封工況都復(fù)雜[15]。目前，針對(duì)部分冰蓋下水流流速垂向分布的理論分析較少。圖6給出B10、B15、B20和B30工況各垂線上的縱向時(shí)均流速的垂向分布。由圖6可見，部分冰封工況縱向時(shí)均流速垂向分布兼具明流和完全冰封工況的特征：① 在冰蓋區(qū)，所有垂線均類似于完全冰封工況，采用C40工況mi和mb的值，由式(2)得到冰蓋區(qū)縱向時(shí)均流速的計(jì)算值，并同實(shí)測(cè)值對(duì)比，結(jié)果表明除B20工況Y15垂線外，各垂線的EMAP(u)均小于5.0%，IFO(u)均大于0.90，表明部分冰蓋下水流冰蓋區(qū)縱向流速分布可采用雙冪律描述；② 明流區(qū)部分垂線類似于明流工況，采用明流O00工況冪指數(shù)mb的值，由式(2)得到明流區(qū)部分垂線縱向流速垂向分布，各垂線EMAP(u)均小于5.0%，IFO(u)均大于0.90，表明部分冰蓋下水流明流區(qū)靠近紊流核心區(qū)的流速垂向分布符合單冪律；③ 在兩者間存在過渡區(qū)，縱向時(shí)均流速垂向分布既不同于明流工況，也不同于完全冰封工況。根據(jù)實(shí)測(cè)值按雙冪律進(jìn)行擬合，得到的mb均小于C40工況的值，但相反地mi均大于C40工況的值。

圖6 部分冰封工況縱向時(shí)均流速垂向分布Fig.6 Vertical profiles of streamwise velocity in open channels with partial ice cover

2.2 雷諾應(yīng)力垂向分布

聯(lián)立連續(xù)性方程和動(dòng)量方程，得到恒定均勻流的控制方程為[23]

(6)

(7)

(8)

(9)

(10)

(11)

式中：u*為斷面摩阻流速，這表明紊流核心區(qū)雷諾應(yīng)力在水深方向上呈線性分布。圖7(a)對(duì)比了明流O00工況Y40垂線上雷諾應(yīng)力實(shí)測(cè)值和計(jì)算值，結(jié)果表明，明流工況紊流核心區(qū)的雷諾應(yīng)力在水深方向上呈線性分布，雷諾應(yīng)力的最大值出現(xiàn)在床面附近。但是，對(duì)應(yīng)于最大流速點(diǎn)因二次流影響偏離自由水面，雷諾應(yīng)力零值點(diǎn)也偏離自由水面。

(12)

采用斷面摩阻流速對(duì)切應(yīng)力進(jìn)行量綱一化。采用試驗(yàn)測(cè)量的斷面尺寸和比降，通過下列公式計(jì)算得到u*：

(13)

(14)

式中：hb為床面層的水深；c1、c2為回歸系數(shù)。顯然，切應(yīng)力零值點(diǎn)相對(duì)水深不同于最大流速點(diǎn)，兩者相對(duì)位置與床面、冰蓋的相對(duì)粗糙程度有關(guān)。以C40工況為例，由式(14)計(jì)算的切應(yīng)力零值點(diǎn)相對(duì)水深為0.399，位于最大流速點(diǎn)和相對(duì)光滑的固壁之間。

圖7 O00工況和C40工況雷諾應(yīng)力垂向分布Fig.7 Vertical profiles of shear stress for Run O00 and Run C40

對(duì)于部分冰封工況，由于冰蓋區(qū)、明流區(qū)的垂向流速分布差異較大，存在橫向動(dòng)量交換，由式(7)推求雷諾應(yīng)力時(shí)不能忽略二次流的影響。目前尚未建立雷諾應(yīng)力的數(shù)學(xué)模型，僅能從測(cè)量數(shù)據(jù)探究其垂向分布規(guī)律。以部分冰封B20、B30工況為例，冰蓋邊緣垂線上雷諾應(yīng)力分布見圖8。受橫向動(dòng)量交換影響，雷諾應(yīng)力在水深方向上散亂分布，不滿足線性分布的規(guī)律。

圖8 部分冰封工況雷諾應(yīng)力垂向分布Fig.8 Vertical profiles of shear stress in open channels with partial ice cover

2.3 紊動(dòng)強(qiáng)度垂向分布

紊動(dòng)強(qiáng)度是流場(chǎng)內(nèi)任一空間點(diǎn)脈動(dòng)流速的二階中心矩，表示瞬時(shí)流速在其均值附近的分散程度，反映水體流速脈動(dòng)的強(qiáng)弱程度。對(duì)于明流工況，Nezu等[25]采用指數(shù)律描述其縱向紊動(dòng)強(qiáng)度垂向分布：

(15)

式中：urms為縱向紊動(dòng)強(qiáng)度；Du和ξu為經(jīng)驗(yàn)系數(shù)。根據(jù)明流O00工況實(shí)測(cè)的瞬時(shí)流速，計(jì)算得到各垂線的縱向紊動(dòng)強(qiáng)度，見圖9(a)。從紊動(dòng)強(qiáng)度垂向分布來看，在床面附近存在最大值，隨水深增加逐漸減小。當(dāng)相對(duì)水深z/H>0.5后，紊動(dòng)強(qiáng)度總體上趨于穩(wěn)定，在自由水面附近達(dá)到最小。由圖9(a)可見，明流O00工況，各垂線紊動(dòng)強(qiáng)度垂向分布符合指數(shù)律，EMAP(urms)均值為9.43%，IFO(urms)均值為0.895 6。從橫向分布來看，總體上越靠近邊壁，紊動(dòng)強(qiáng)度越大，越靠近紊流核心區(qū)，紊動(dòng)強(qiáng)度越小。結(jié)合縱向流速的垂向分布，紊動(dòng)強(qiáng)度在最大流速附近達(dá)到最小值，說明最大流速處流速脈動(dòng)較小，相鄰流層間的動(dòng)量交換最弱。

圖9 O00工況和C40工況紊動(dòng)強(qiáng)度垂向分布Fig.9 Vertical profiles of turbulence intensity for Run O00 and Run C40

對(duì)于完全冰封工況，按照雙層假定將過流斷面在垂向上劃分為床面層和冰蓋層，在各流層中紊動(dòng)強(qiáng)度垂向分布均符合指數(shù)律[26]，則

(16)

式中：Du，b、Du，i、ξu，b和ξu，i為經(jīng)驗(yàn)系數(shù)。根據(jù)C40工況實(shí)測(cè)的瞬時(shí)流速，計(jì)算得到各垂線的紊動(dòng)強(qiáng)度，見圖9(b)。在垂向上，紊動(dòng)強(qiáng)度在床面和冰蓋附近存在最大值，隨著到無滑移邊界距離的增加而逐漸減小，在最大流速附近達(dá)到最小值。經(jīng)回歸分析，C40工況各垂線上床面層和冰蓋層的縱向紊動(dòng)強(qiáng)度均滿足指數(shù)分布的規(guī)律，EMAP(urms)=6.94%，IFO(urms)=0.920 3。

對(duì)于部分冰封工況，目前尚未提出描述其紊動(dòng)強(qiáng)度的分布律，僅能通過分析試驗(yàn)測(cè)量數(shù)據(jù)探究其垂向分布規(guī)律。圖10給出了B10、B15、B20和B30工況各測(cè)線上紊動(dòng)強(qiáng)度垂向分布，總體上無滑移邊界(床面、冰蓋)附近的紊動(dòng)強(qiáng)度相對(duì)較大，且邊界相對(duì)粗糙的冰蓋區(qū)紊動(dòng)強(qiáng)度大于邊界相對(duì)光滑的床面區(qū)。部分冰封4種工況冰蓋區(qū)urms/u*均值分別為1.68、1.87、1.90和1.79，明流區(qū)urms/u*均值分別為1.41、1.66、1.78和1.62，冰蓋區(qū)紊動(dòng)強(qiáng)度大于明流區(qū)紊動(dòng)強(qiáng)度。當(dāng)F≤0.5時(shí)，冰蓋區(qū)、明流區(qū)urms/u*均隨F的增加而增大；當(dāng)F>0.5時(shí)，相比B20工況，B30工況冰蓋區(qū)、明流區(qū)urms/u*均減小，紊動(dòng)強(qiáng)度減弱。

3 結(jié) 論

在普通實(shí)驗(yàn)室循環(huán)水槽系統(tǒng)中開展模型試驗(yàn)，采用三維聲學(xué)多普勒測(cè)速儀測(cè)量不同冰蓋覆蓋度下恒定均勻流場(chǎng)，從縱向流速、雷諾應(yīng)力和紊動(dòng)強(qiáng)度等方面分析冰蓋下水流紊動(dòng)特性，主要得到以下結(jié)論：

(2) 明流工況雷諾應(yīng)力在水深方向上符合線性分布，最大值出現(xiàn)在床面附近；完全冰封工況雷諾應(yīng)力在水深方向也符合線性分布，切應(yīng)力零值點(diǎn)位于最大流速點(diǎn)和相對(duì)光滑邊界之間；部分冰封工況受橫向動(dòng)量交換的影響，雷諾應(yīng)力在水深方向上不符合線性分布。

(3) 明流工況紊動(dòng)強(qiáng)度垂向分布符合指數(shù)律，完全冰封工況基于雙層假定劃分的冰蓋層和床面層中紊動(dòng)強(qiáng)度可分別采用指數(shù)律描述，部分冰封工況紊動(dòng)強(qiáng)度受冰蓋影響在橫向上存在明顯的差異。