封凍河道冰下流速分布規律研究

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

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

1 引言

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

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

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

2 控制方程與離散

2.1 控制方程

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

圖1 坐標系示意

連續方程:

x方向動量方程:

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

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

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

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

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

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

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

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

2.2 離散求解策略

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

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

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

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

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

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

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

圖2 垂向網格示意

式中:θ為權重系數；k為垂向分層數。

代入公式(9),得到

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

其中

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)］結合床面邊界條件[式(14)］、自由水面邊界條件[式(16)］,最終形成如下三對角方程組,該方程組可用“追趕法”求解。

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

3 計算結果分析

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

圖5 計算值與文獻[7］流速分布比較

圖6 計算值與文獻[8］流速分布比較

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

圖8 對數坐標系下流速分布結果(文獻[8］中,ni/nb＝0.76)

圖9 對數坐標系下流速分布結果(文獻[8］中,ni/nb＝1)

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

4 結論

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