部分冰蓋下水流水力特性數值模擬

陳 剛, 金 棟, 賴彥丞, 李映槿

(1.河海大學 水文水資源與水利工程科學國家重點實驗室, 江蘇 南京 210098; 2.云南省水利水電勘測設計研究院, 云南 昆明 650021)

1 研究背景

河冰是陸地冰凍圈的重要組成部分[1],在其結冰-封凍-解凍周期性地生消演變中影響著水流流動邊界[2],形成明流、完全冰封、部分冰封的不同流態[3-4]。封凍期水流表面被連續冰蓋完全覆蓋時,連續的附加無滑移邊界主要在垂向上改變水流的運動規律,形成以最大流速點為界的兩層結構[5-6];在結冰期和解凍期,水流表面一部分被冰覆蓋,另一部分為自由水面,不連續的附加無滑移邊界同時在橫向和垂向上改變水流的運動規律,此時水流兼具明流水流和完全冰封水流的水力特性[7-8]。岸冰是部分冰蓋的常見形式[9],按其形成時間和條件的不同,分為初生岸冰(結冰期)、固定岸冰(結冰期)、沖積岸冰(結冰期、融冰期)、再生岸冰(融冰期)、殘余岸冰(融冰期)[10]。然而,相比明流和完全冰封水流,關于部分冰封水流水力特性的研究較少。水流流態在水力、熱力的共同作用下急劇轉變是冰凌致災的主要原因,全面掌握包括部分冰封在內的所有流態的水力特性,才能更好地為寒區河流開發利用提供理論指導。

獲取河冰生消演變對水流特性影響最直接的方法是開展原型觀測,但在上游來水、河冰生消、河床演變的共同影響下,流動邊界條件千差萬別,原型觀測得到的水流運動規律經驗性較強,難以推廣應用[11]。水槽試驗能夠直觀地反映特定工況下水流流動特性,但是受限于斷面測點布置和儀器量程,往往難以揭示整個流場的全部細節。數值計算可全面地反映所研究流場的各種細節,能獲取許多原型觀測和模型試驗研究難以測量的數據,已被廣泛應用于冰蓋下水流結構的研究[12-13]。Lau等[14]通過建立k-ε湍流模型,研究了冰蓋形成對水流垂向混合特性的影響;茅澤育等[15]采用k-ε湍流模型模擬冰蓋下水流和物質輸移的垂向二維分布,模擬結果表明連續冰蓋下水流紊動黏滯系數并未出現零值點,水流縱向流速在流動核心區較為均勻,說明雙對數律的假定存在問題;王軍[16]采用標準k-ε湍流模型研究冰蓋形成對水流流速分布的影響,通過模擬試驗發現最大流速點的偏移范圍有限;王志興等[17]采用三維k-ε紊流模型模擬連續冰蓋下水流縱向流速垂向分布,模擬結果表明雙冪律待定參數的取值具有規律性;蘇磊等[18]建立基于k-ε紊流模型的冰蓋下水流垂向二維數值模型,探討了連續冰蓋下水流斷面流速分布的主要影響因素?？梢?現有研究多針對連續冰蓋,而采用數值模擬方法研究部分冰蓋下水流水力特性的成果還較少。

本研究采用Fluent軟件建立部分冰蓋下水流的數值模擬模型,并采用水槽試驗數據對模型進行校驗,從縱向流速分布、二次流結構、雷諾應力分布、紊動能分布等方面揭示部分冰蓋下水流的水力特性。

2 數學模型

2.1 控制性方程

對于笛卡爾坐標系中的不可壓縮流體的定常流動,紊流時均控制方程為[12]:

(1)

(2)

式中:ρ為流體的密度,kg/m3;u為雷諾平均流速,m/s;u′為脈動流速,m/s;p為壓強,Pa;σij為應力張量分量,Pa; 下標i,j分別表示相應變量在xi,xj方向上的分量。

根據牛頓流體的應力-應變率本構關系,有:

(3)

式中:δij為Kronecker符號;μ為分子黏性系數,Pa·s; 下標l表示相應變量在xl方向上的分量。

由于控制性方程不封閉,因而在雷諾應力模型中引入雷諾輸運方程、紊動能方程和紊動能耗散率方程使其封閉。紊動能和紊動能耗散率方程分別為:

(4)

(5)

式中:k為紊動能;ε為紊動能耗散率;Gk為紊動能生成項(剪切產生率);Cμ為無量綱常數(Cμ=0.09),Cε1和Cε2為源項中的經驗常數,取Cε1=1.44,Cε2=1.92。

雷諾應力輸運方程為[12]:

∏ij-Eij+Fij

(6)

(7)

(8)

(9)

(10)

(11)

(12)

式中:DT, ij為紊動擴散項;DL, ij為分子黏性引起的擴散項;Pij為應力產生項; ∏ij為壓力應變項;Eij為黏性耗散項;Fij為系統旋轉產生項。

2.2 邊界條件

(1)進口邊界條件:水流進口采用質量流進口邊界,數值模擬與水槽試驗相應工況的邊界條件相同。進口邊界的紊動能和紊動能耗散率分別采用公式(13)、(14)確定。

(13)

(14)

式中:I為紊動強度,結合試驗得到的紊動強度分布規律參照Nezu等[19]的公式綜合擬定;Cμ為常數,取Cμ=0.09;l為紊流特征尺度,l=0.07DH,其中DH為水力直徑。

DH的計算公式為:

DH=4A/χ

(15)

(2)出口邊界條件:采用壓力出口邊界。

(3)壁面邊界條件:床面、邊壁和冰蓋均采用無滑移邊界,近壁區采用標準壁面函數進行處理。床面和冰蓋的阻力系數采用粗糙高度表示,兩者的轉換如公式(16)所示[20]。

(16)

式中:n1、n2分別為床面和冰蓋的曼寧糙率系數,s·m-1/3;δ1、δ2分別為床面和冰蓋的粗糙高度,m。

(4)對稱邊界條件:在對稱邊界上,垂直邊界的速度取為0,其他物理量的值在對稱邊界內外相等。水面采用滑移壁面條件,即切應力為零。

2.3 計算工況

為便于采用水槽試驗實測數據校驗模型,數值模擬計算工況與水槽試驗相同。試驗在昆明理工大學水流與水工結構重點實驗室的玻璃水槽試驗系統中進行,玻璃水槽長度為15 m、寬度為0.8 m、高度為0.35 m,坡降固定為0.5‰,采用泡沫板制作模擬冰蓋。試驗時,采用矩形薄壁堰測量流量、水位測針測量斷面水深、聲學多普勒點式三維測速儀Vectrion測量流速。冰蓋模擬為對稱岸冰,以通過冰蓋邊緣的垂線為界,劃分為被岸冰覆蓋的冰蓋區和表面為自由水面的明流區。定義冰蓋覆蓋度F為冰蓋寬度與斷面寬度的比值,即F=2b1/B(b1為單側冰蓋寬度;B為斷面寬度),冰蓋寬度b1依次取為0.10、0.15、0.20和0.30 m,相應的冰蓋覆蓋度F分別為0.250、0.375、0.500和0.750。為進行對比,增加冰蓋覆蓋F分別為0、1的特殊工況,分別對應明流工況和完全冰封工況,相應的冰蓋寬度分別為0和0.40 m,數值計算工況見表1。測量時,每間隔5 cm布置1條垂線,依次記為Y05、Y10、Y15、Y20、Y25、Y30、Y35、Y40共8條垂線,每條垂線上布置10個測點。由于流場為對稱分布,為節省計算資源,僅計算半寬區域,上游計算長度為斷面寬度的4倍,下游計算長度為斷面寬度的6倍。

2.4 方程離散與網格劃分

采用計算流體軟件Fluent進行計算,其控制方程的離散方法為有限體積法,選取二階迎風格式作為空間離散格式,選取SIMPLEC(semi-implicit method for pressure linked equations consistent)算法處理壓力與速度的耦合關系。計算區域采用笛卡爾坐標系,坐標原點設置在測量斷面右側,采用非均勻結構化網格進行劃分,固壁(床面、邊壁、冰蓋)附近網格進行適當加密,計算網格為六面體結構化正交網格(圖1)。經網格無關性檢驗,選取網格總數為3.6×106(60×60×1000)的非均勻網格。各方程的收斂控制精度均為1×10-5。

圖1 模型計算網格劃分示意圖

2.5 模型校驗

為了檢驗數值模型的可靠性,采用縱向流速的數值模擬結果與實測值進行對比。為便于作圖,采用縱向時均流速垂向分布律對各垂線上的實測數據進行擬合。

無冰蓋的明流工況O00采用單冪律擬合[21],其數學表達式為:

(17)

式中:umax為最大流速,m/s;H為水深,m;z為測點距床面的距離,m; 1/m為指數,其值與雷諾數有關。umax、m值可通過水槽試驗實測資料擬合確定。

完全冰封工況采用雙冪律進行擬合[22],其數學表達式為:

(18)

式中:k0為與流量有關的常數,m/s;m1,m2分別為與床面、冰蓋粗糙高度有關的冪指數,可采用水槽試驗實測數據擬合確定。

采用平均相對誤差和擬合度指標評價理論公式性能。平均百分比誤差(mean absolute percentage error,MAPE)為相對誤差絕對值的算術平均值,其計算公式為:

(19)

式中:i為數據點數;N為數據量;us為模型模擬值;um為理論計算值。MAPE最小值為0,其值越小說明模型模擬值越接近參考數據。

擬合度指標是常用于判定非線性回歸方程擬合度的統計參數,其計算公式為:

(20)

式中:FOI為擬合度指標,FOI最大值為1,其值越接近1說明模型模擬值與實測值的吻合程度越好。

各工況縱向時均流速模擬值與實測值對比見圖2。通過計算得出各工況縱向時均流速的MAPE分別為4.16%、3.80%、2.57%、2.60%、4.00%和4.17%,FOI平均值分別為0.957 1、0.958 8、0.966 2、0.968 8、0.946 6和0.940 2,表明各工況的數值模型均得以校驗。

3 結果分析與討論

3.1 縱向時均流速分布

采用已校驗的模型模擬各計算工況的縱向時均流速斷面分布,如圖3所示。在明流工況O00中(圖3(a)),除近壁區流速較小外,在流核區流速相對均勻,模型驗證結果表明縱向時均流速在水深方向上滿足單冪律分布(圖2(a)),最大流速位于自由水面附近。對于完全冰封C40工況(圖3(f)),連續冰蓋作為附加無滑移邊界代替自由水面后,最大流速偏向斷面中部,縱向時均流速在水深方向呈不對稱“”型分布,可采用雙冪律描述(圖2(f))。對于部分冰封工況(圖3(b)～3(e)),岸冰使得流速在斷面上重新分配,總體上明流區流速增大,冰蓋區流速減小,兩者間存在明顯的橫向流速差。從圖2、3可以看出,部分冰封工況的冰蓋區流速垂向分布規律類似于完全冰封工況,可采用雙冪律描述;明流區部分垂線流速分布類似于明流工況,可采用單冪律描述;在靠近冰蓋邊緣的區域存在過渡區,最大流速點偏向斷面中部,縱向時均流速垂向分布規律不能采用單冪律描述,特別是冰蓋覆蓋度較大的工況,縱向時均流速在垂向上以不對稱“”型分布為主,具有明顯的二維特征。

圖3 不同工況冰蓋下水流縱向時均流速等值線斷面分布

Peters等[3, 7]對不同寬度岸冰下水流流速分布的試驗測量結果表明,明流區和冰蓋區間存在過渡區,其流速垂向分布既不同于明流區的單冪律分布,也不同于冰蓋區的雙冪律分布。圖3所示的模型模擬結果很好地展示了3種工況類型的水流縱向時均流速分布特性,特別是部分冰蓋水流同時存在明流區、冰蓋區和過渡區的二維特性。這表明采用的雷諾應力模型能夠很好地模擬3種不同冰蓋類型下水流的縱向時均流速分布特性。

3.2 二次流結構

圖4給出了不同工況冰蓋下二次流結構的模擬結果。圖4(a)所示的明流工況O00斷面上共有2個渦:底部渦A和表面渦B,以渦A為主,通過2個反向渦將流核區的高能流體帶到邊角區。隨著冰蓋覆蓋度的增加,渦體形狀、數量、大小和位置發生明顯變化,流速等值線發生彎曲。對于B10工況(圖4(b)),斷面上共形成3個渦體,其中渦A和渦B為橫跨冰蓋區和明流區的反向渦體,以冰蓋附近的渦B為主,將明流區的高能流量帶至冰蓋區,流速等值線凸向冰蓋區,渦C作用較弱。對于B15工況(圖4(c)),渦數量仍為3個,但冰蓋附近的渦B主要位于冰蓋區,以渦A為主,流速等值線凸向冰蓋區。對于B20工況(圖4(d)),冰蓋附近的渦B完全位于岸冰下,并偏向冰蓋邊角,渦C逐步增強,位置偏向冰蓋區和明流區交界的冰蓋邊緣附近。對于B30工況(圖4(e)),渦B強度減弱,位于冰蓋邊角,并且在冰蓋邊緣處形成新的渦體D,此時以床面邊角附近的渦A為主。對于C40工況(圖4(f)),僅存在兩個反向渦體,通過反向渦將流核區高能流體帶至邊壁附近。

圖4 不同工況冰蓋下水流二次流結構

部分冰封B10、B15、B20和B30工況渦的數量均多于明流O00工況和完全冰封C40工況的,且渦體的大小、位置各不相同。這表明冰蓋下水流二次流結構復雜,二次流結構隨冰蓋覆蓋度的增大而變化;由于冰蓋區流速小而明流區流速大,兩者間存在橫向流速梯度,通過二次流形成橫向動量交換,這也驗證了Peters等[3, 7]的試驗測量結果。

3.3 雷諾應力分布

雷諾應力是流體質點脈動導致紊動微團在相鄰流層間交換所產生的附加切應力。圖5給出了不同工況冰蓋下的雷諾應力的模擬值。

圖5 不同工況冰蓋下水流雷諾應力等值線分布

明流O00工況的雷諾應力最大值位于床面附近(圖5(a)),隨著距床面距離的增大而逐漸減小,總體符合線性分布。忽略寬淺斷面中橫向變化和二次流項,通過推導可得到明流雷諾應力的計算方法如公式(21)所示[23]。公式(21)表明紊流核心區雷諾應力在水深方向上呈線性分布,圖5(a)所示的模型模擬結果與理論分析所揭示的規律相符。

(21)

式中:τxz為x-z平面上的切應力,Pa;τ01為床面切應力,Pa;H為水深,m;z為距床面的距離,m。

對于C40工況(圖5(f)),雷諾應力在床面和冰蓋附近取得極值,床面附近為正值,冰蓋附近為負值,在水深方向上總體呈線性分布,零值點位于流核區。同理,推導得出的完全冰封水流雷諾應力計算方法如公式(22)所示[23]。公式(22)表明紊流核心區雷諾應力在水深方向上也呈線性分布,圖5(f)所示的模型模擬結果與理論分析所揭示的規律相符。

(22)

式中:τ02為冰蓋切應力,Pa。

對于部分冰封工況(圖5(b)～5(e)),冰蓋附近雷諾應力為負值,雷諾應力負值等值線凸向明流區,雷諾應力在斷面上分布不均勻,隨著冰蓋覆蓋度的增加,負值區范圍逐步擴大,雷諾應力垂向分布不符合線性分布。這是因為水流橫向動量交換較顯著,使二次流項的影響不可忽略,從側面反映了部分冰蓋工況水流結構的復雜性。

3.4 紊動能分布

紊動能是常用于直觀表征水流整體紊動狀況的物理量。不同工況冰蓋下水流紊動能的數值模擬結果見圖6。明流O00工況(圖6(a)),紊動能在固壁床面和邊壁附近取得極值,這表明固壁附近瞬時流速在其均值附近較分散,反映出水體流速脈動在固壁附近較強[24]。對于完全冰封C40工況(圖6(f)),紊動能在床面和冰蓋附近取得最大值,流核區紊動最弱,紊動能大小與固壁的粗糙度有關,固壁越粗糙,紊動能越大。部分冰封工況的紊動能在斷面上呈現明顯的差異(圖6(b)～6(e)),受床面、冰蓋和邊壁的共同影響,冰蓋區紊動能明顯大于明流區紊動能,紊動能等值線凸向冰蓋區,明流區紊動能在床面附近取得最大值,在自由水面附近紊動最弱。對比部分冰蓋B10、B15、B20和B30工況,岸冰替代部分自由水面成為附加的無滑移邊界,水流流動阻力增大,橫向流速差所產生的橫向動量交換也會消耗水流的部分動能,使得水流的動能減小,這類似于復式斷面明渠水流[4]。因此,在估算部分冰封斷面過流能力時,應考慮橫向動量交換的影響。

4 結 論

采用計算流體力學軟件Fluent中的雷諾應力模型模擬不同覆蓋度下恒定均勻流的水力特性,采用經水槽試驗測量數據校驗的數值模型,分析縱向時均流速分布、二次流結構、雷諾應力、紊動能的斷面分布規律,揭示出部分冰蓋下水流水力特性,主要結論如下:

(1)岸冰的形成會導致冰蓋下水流流速在斷面上重新分配,明流區水流縱向時均流速平均值增大,冰蓋區水流縱向時均流速平均值減小,明流區與冰蓋區之間存在明顯的流速差;相比明流和連續冰蓋下水流,部分冰蓋下的水流結構更加復雜,二次流渦體的形狀、數量、大小和位置隨著岸冰覆蓋度的增加而發生明顯變化。

(2)明流工況(F=0)和完全冰封工況(F=1)的水流雷諾應力沿水深方向呈線性分布。部分冰封工況(0

(3)水流紊動能在固壁附近取得極大值,相比明流工況(F=0)和完全冰封工況(F=1),由于部分冰封工況(0