蘇 磊，藺建銘，張寶森，謝志剛，田治宗，鄧 宇，于國卿

(1.黃河水利委員會 黃河水利科學(xué)研究院，河南 鄭州 450003;2.鄂爾多斯市水利事業(yè)發(fā)展中心，內(nèi)蒙古 鄂爾多斯 017000）

1 引言

高寒地區(qū)冬季氣溫低,河渠易結(jié)冰而形成冰蓋。冰蓋顯著改變河渠斷面的水流結(jié)構(gòu),使封凍河道各水力要素變得十分復(fù)雜[1］,河道邊界條件、水動力條件變化導(dǎo)致冰下水位、流量、流速變化,將影響冬季河渠輸水,甚至誘發(fā)凌汛災(zāi)害,因此研究封凍河道水流特性具有重要的理論及實際意義[2］。

國內(nèi)外學(xué)者對冰下流速分布問題進(jìn)行了大量的研究,冰凌數(shù)學(xué)模型方面,國內(nèi)外開發(fā)了一系列一維動/靜態(tài)模型(RICE、ICEJAM、RIVER1D、HEC-RAS等)、二維動/靜態(tài)模型(CRISSP2-D、DYNARICEDE[3］等)。王軍等[4］、楊開林等[5］、茅澤育等[6］結(jié)合水力學(xué)、熱力學(xué)原理改進(jìn)了河道冰凌數(shù)學(xué)模型,但針對冰蓋下水流垂向流速分布特性的研究相對較少。

本文對N-S方程進(jìn)行簡化、離散,采用k-ε紊流模型建立冰蓋下水流垂向二維數(shù)學(xué)模型,模擬封凍河道冰蓋下流速分布情況,并與實測資料進(jìn)行比較,分析、探討其變化規(guī)律。

2 控制方程與離散

2.1 控制方程

如圖1所示,x表示流動方向坐標(biāo),z表示垂直于河床方向坐標(biāo)。對于冰蓋下的明渠紊流,考慮靜壓假定的不可壓縮流體的控制方程無量綱化后,可以表示如下。

圖1 坐標(biāo)系示意

連續(xù)方程:

x方向動量方程:

在z方向上引入σ坐標(biāo),某一坐標(biāo)處的垂向質(zhì)點流速w可用垂向場流速ω表示為

式中:u、w分別為x、z方向的時均流速；t為時間；ζ為自由水面高程(ζ＝zb＋H),zb為床面高程；R0為羅斯比數(shù)；為科氏系數(shù),g為重力加速度；ρ為密度；p為壓強；AH為水平向紊動渦黏系數(shù)；AV為垂向紊動渦黏系數(shù)；H為水深。

此外,式中定義如下無量綱組合:①Rossby數(shù),R0＝Ur/fXr；Froude數(shù),Fr＝Ur/(gZr)1/2。②密度Froude數(shù),。③水平向Ekman數(shù),EH＝AHr/；垂向Ekman數(shù),EV＝AVr/。其中下標(biāo)r表示對應(yīng)變量的參考值。

垂向紊動渦黏系數(shù)AV采用垂向k-ε模型進(jìn)行計算,該模型可以較好地反映風(fēng)剪切應(yīng)力、床面剪切應(yīng)力、流速梯度產(chǎn)生項、耗散、擴散和分層等的影響。k-ε模型的基本思想是將垂向紊動渦黏系數(shù)AV與紊流動能k和紊流耗散率ε聯(lián)系起來,即

紊流動能k和紊流耗散率ε的輸運方程表示如下:

式中:cν、c1、c2、σk、σε為常數(shù),取值cν＝0.09、c1＝1.44、c2＝1.92、σk＝1.0、σε＝1.3；源項P為流速垂向梯度引起的紊流產(chǎn)生項；G為浮力產(chǎn)生或耗散項。

河床阻力的確定,以x方向為例,河床切應(yīng)力在x方向表示為,其中:U、V分別為x、y方向的平均流速；cf為(無量綱)床面摩阻系數(shù),cf＝gn2/H1/3＝g/C2,n為Manning糙率系數(shù),C為Chézy系數(shù)；mb反映了河床底坡的影響。

冰蓋阻力的確定,采用H.A.Einstein提出的阻力劃分方法可將流體分為兩層等效明流層(冰面層和床面層),假定上下兩層互不影響,分別只受冰蓋底部和床面粗糙度的影響,與表示河床阻力相同的方式,冰蓋底部的切應(yīng)力在x方向表示為。

2.2 離散求解策略

將x方向流速u表示為如下形式:

式中:u′為x方向流速相較于垂向平均流速的差值,即攝動流速。

u′所滿足的方程可通過將三維控制方程減去二維控制方程得到,這樣的處理移除了原始方程中的自由水面重力波項,因此求解u′后再通過式(7)間接計算u,能使模型獲得更高的穩(wěn)定性,保證計算質(zhì)量。

對于攝動動量方程,除時間項外,只有垂直擴散項被視為隱式的,其余各項(包括科氏力項、對流項、水平擴散項和斜壓項)均被視為顯式的,也就是說,攝動動量方程可以簡寫為如下形式:

式中:Rx為x方向的水平源項,包括科氏力項、對流項、水平擴散項和斜壓項4部分。

將式(8)采用有限體積法進(jìn)行離散,將其在一個時間步長Δt和一個控制體積Δσ上進(jìn)行積分可得

如圖2所示,假設(shè)變量在中心點的值代表該變量在整個控制體積Δσ上的積分平均值,并認(rèn)為隱式項在n、n＋1兩個時刻的加權(quán)平均值代表其在整個時間步長Δt上的積分平均值,顯式項在n時刻的值代表其在整個時間步長Δt上的積分平均值,則上式可表示為

圖2 垂向網(wǎng)格示意

式中:θ為權(quán)重系數(shù)；k為垂向分層數(shù)。

代入公式(9),得到

上式即式(8)的最終離散形式,并對k＝2,3,…,m-1成立,其中m為垂向分層總數(shù),并可寫成通用形式:

其中

2.3 定解條件及方程組的求解

2.3.1 邊界條件

(1)床面邊界條件(k＝1),如圖3所示。

圖3 床面邊界條件示意

式(8)在床面附近所需滿足的邊界條件為

代入式(9),得到

式(13)寫成通用形式為

其中

(2)冰蓋邊界條件(k＝n),如圖4所示。

圖4 冰蓋邊界條件示意

式(8)在冰蓋附近所需滿足的邊界條件為

代入式(9),得到

式(9)寫成通用形式為

其中

2.3.2 初始條件

2.3.3 方程組求解

離散方程[式(12)］結(jié)合床面邊界條件[式(14)］、自由水面邊界條件[式(16)］,最終形成如下三對角方程組,該方程組可用“追趕法”求解。

求解上述三對角方程組獲得攝動流速u′后,可由式(7)計算得到每層流速u。整個計算的迭代收斂條件為,其中ε取10-4。

3 計算結(jié)果分析

應(yīng)用實測資料對上述模型計算結(jié)果進(jìn)行了驗證,如圖5、圖6所示,圖中h為水深、y為相對水深、為斷面平均流速。其中,圖5為計算值與文獻(xiàn)[7］中現(xiàn)場實測資料(River Moskva上長為400 m的封凍順直河段第4號量測斷面垂線流速分布實測成果)的比較情況,圖6(a)為計算值與文獻(xiàn)[8］中室內(nèi)水槽實測資料(槽寬80 cm,水深8.92 cm)的比較情況,圖6(b)為計算值與文獻(xiàn)[8］中第12號水流條件(水深200 cm,冰蓋糙率與床面糙率相同)的比較情況。結(jié)果顯示,模型的計算值與實測值均吻合較好。同時,流速分布對于冰蓋糙率ni和床面糙率nb的影響較敏感,最大流速Umax的位置偏向糙率較小的一側(cè),且在相同的水流條件下,ni/nb的值越小,最大流速點的位置越靠近糙率較小的一側(cè)；當(dāng)冰蓋糙率和床面糙率相等時,流速分布以最大流速Umax為軸呈對稱分布,且在最大流速Umax處流速梯度幾乎為零[9］。

圖5 計算值與文獻(xiàn)[7］流速分布比較

圖6 計算值與文獻(xiàn)[8］流速分布比較

將冰蓋下水流以最大流速為界分為床面層和冰面層兩層,分層示意如圖7所示,hb為床面層水深,yb為床面層相對水深,hi為冰面層水深,yi為冰面層相對水深。以文獻(xiàn)[8］中兩組工況為例子,工況1為床面糙率和冰蓋糙率不相等時(ni/nb＝0.76),將其計算結(jié)果繪于圖8所示的對數(shù)坐標(biāo)系下；工況2為床面糙率和冰蓋糙率相等時(ni/nb＝1),將其計算結(jié)果繪于圖9所示的對數(shù)坐標(biāo)系下。由圖8可知,工況1當(dāng)ln(yb/hb)＜-0.4時(即yb/hb＜0.67時),床面層流速分布符合對數(shù)分布,當(dāng)ln(yb/hb)＞-0.4時(即yb/hb＞0.67時),床面層流速分布將不再符合對數(shù)分布,且流速值略小于符合對數(shù)分布的情況；在冰面層有同樣的規(guī)律但影響區(qū)域略有不同,當(dāng)ln(yi/hi)＜-0.32時(即yi/hi＜0.72時),冰面層流速分布符合對數(shù)分布；ln(yi/hi)＞-0.32時(即yi/hi＞0.72時),冰面層流速分布將不再符合對數(shù)分布,且流速值略小于符合對數(shù)分布的情況。同時,由圖9可知,工況2當(dāng)ln(yb/hb)＜-0.41時(即yb/hb＜0.66時),床面層和冰面層流速分布符合對數(shù)分布；ln(yb/hb)＞-0.41時(即yb/hb＞0.66時),床面層和冰面層流速分布將不再符合對數(shù)分布。

圖8 對數(shù)坐標(biāo)系下流速分布結(jié)果(文獻(xiàn)[8］中,ni/nb＝0.76)

圖9 對數(shù)坐標(biāo)系下流速分布結(jié)果(文獻(xiàn)[8］中,ni/nb＝1)

計算結(jié)果表明,封凍河道的斷面流速分布主要與床面糙率、冰蓋糙率及河道的水流條件等影響因素相關(guān)。床面糙率與冰蓋糙率相等和床面糙率與冰蓋糙率不相等兩種工況,冰面層和床面層均有70%左右的流速符合對數(shù)分布,其余靠近最大流速點約30%范圍內(nèi)的流速不再符合對數(shù)分布且流速值略小于對數(shù)分布的情況。工況1當(dāng)床面糙率與冰蓋糙率相等時,冰面阻力和床面阻力干擾影響的區(qū)域相等,最大流速Umax在水深中點的位置。工況2當(dāng)床面糙率與冰蓋糙率不相等時,冰面阻力和床面阻力干擾影響的區(qū)域不相等,且糙率越大阻力越大,阻力干擾影響的區(qū)域越大,從而導(dǎo)致最大流速Umax位置偏向相對光滑的一側(cè)；Umax的位置與ni/nb比值相關(guān),且ni/nb的值越小Umax位置越靠近相對光滑的一側(cè)。床面糙率與冰蓋糙率相等和不相等兩種情況,計算結(jié)果均表現(xiàn)出相同的規(guī)律,這與前人的試驗成果吻合。

4 結(jié)論

本文基于淺水方程,采用k-ε紊流模型建立冰蓋下水流二維數(shù)值模型,并將其應(yīng)用于封凍河道冰下流速分布的計算。結(jié)果表明,冰蓋下水流流速分布對于冰蓋糙率ni和床面糙率nb的影響較敏感,當(dāng)床面糙率與冰蓋糙率不相等時,冰面阻力和床面阻力干擾影響的區(qū)域不相等,從而導(dǎo)致最大流速Umax位置偏向相對光滑的一側(cè),Umax位置的偏離程度與ni/nb值直接相關(guān)；當(dāng)冰蓋糙率和床面糙率相等時,流速分布以最大流速Umax為軸呈對稱分布。床面糙率與冰蓋糙率相等和不相等兩種情況,計算結(jié)果均顯示在最大流速點附近30%的區(qū)域內(nèi)流速分布并不遵循對數(shù)分布規(guī)律,且流速值略小于對數(shù)分布的情況。