河流冰蓋熱力增長的數值模擬

鄭成龍， 陳胖胖， 王 軍

（合肥工業大學 土 木與水利工程學院，安徽 合 肥 230009）

河冰的變化過程是水力、熱力和力學平衡綜合作用的結果。當氣溫降至結冰點及其以下時，水體熱量由于紊動而散發，小的冰晶相互聚集形成冰花。冰花的形成是高緯度江河、湖泊、海洋中冰增長的重要組成［1］。水面冰向下游輸送時，會因水面流冰堵塞而停止，由于擁擠作用產生水面流冰堵塞，冰蓋開始形成。冰蓋厚度的變化取決于冰蓋上表面和底面上的熱量交換過程。冰花顆粒在冰蓋下面堆積可形成水內冰冰塞，即懸冰壩。冰塞是寒冷地區江河產生的嚴重冰情現象之一。

對河冰過程研究的方法主要有原型觀測、試驗研究和數值模擬等［2－5］。數值模擬近年來在國內外發展很快。文獻［6］對冰蓋下冰塞堆積進行了模擬；文獻［7］對河道的冰塞進行了模擬；文獻［8］對江河的河冰演變進行了模擬；文獻［9］對Lawrence河灣冰厚分布進行了模擬。

文獻［10］介紹了Athabasca River冬季的水內冰冰花顆粒尺寸在1～50mm。雖然冰花顆粒尺寸因水流條件、氣溫等條件而有所差異，但和河道宏觀尺寸相比，冰蓋下冰花顆粒運動符合一般固液兩相流的特征。

本文基于多相流和傳熱學理論，把河流當作大氣－冰體－水體的一個耦合系統來處理，考慮系統中各界面之間的熱量交換，并耦合冰花在冰蓋下堆積對冰蓋熱量交換的影響，通過數值模擬，分析了入口冰花濃度、入口流速、冰蓋上表面溫度等因素對冰蓋厚度增長的影響。

1 控制方程

連續性方程：

其中，t為時間；αk為第k相的體積分數；ρk為第k相的密度；uk為第k相的平均速度矢量；i為冰相，l為水相；水的密度ρl為997kg／m3，冰的密度ρi為917kg／m3。

動量方程：

其中，p為流體微元體上的壓力；g為重力加速度，9.8m／s2；τ為湍流應力；Mk為單位體積相間動量傳遞；SD為動量源項。

溫度方程：

其中，uj為縱向流速在方向上的分量；T為水溫；ve為有效黏度；σT為紊流模型系數；Sf為原項，水體與冰蓋之間的熱量傳遞通過原項Sf控制實現。

邊界條件如下。

（1）入口邊界條件。假設河流入口流速、冰花濃度等為均勻分布，給定入口水流平均速度v，則入口邊界的湍動能方程為：

湍流動能率方程為：

其中，Cμ為經驗常數，取Cμ＝0.09；Lm為特征長度。

（2）出口邊界條件。k、ε按照坐標局部單向化方式處理，速度v按照局部質量守衡，為使模擬計算更容易收斂，采用壓力出口邊界。

（3）壁面邊界。冰蓋下表面和河床均為無滑移壁面邊界，流體計算域邊界采用考慮粗糙度影響的壁面函數法。

初始條件如下：

文獻［11］對渠道內冰花的體積分數分布進行了數值模擬，渠道出口斷面冰花體積分數為0.005～0.1。假定冰花顆粒在進口均勻分布，冰花體積分數（下文中稱為“冰花濃度”）為0.01～0.07；考慮單一冰花粒徑，取1cm。河流入口處流速為0.5～1.2m／s，且進口冰花平均流速與水流平均流速相同；冰蓋上表面溫度設定為－3、－8、－13、－23℃；入口處冰花溫度設為0℃，入口水溫為0.02℃。

2 模型建立及網格劃分

建立長1 000m，水深2m的矩形河流二維模型，采用結構網格對模型進行網格劃分，將網格劃分為0.1m×0.1m的正方形網格，共10 000×20＝20×104個網格。用有限體積法離散連續性方程、動量方程、標準k－ε方程等，在近壁面處采用壁面定律以考慮邊壁糙率的影響，采用Simple算法。計算模型及模型局部網格劃分，如圖1、圖2所示。

圖1 計算模型

圖2 計算模型局部網格劃分

計算過程中紊流的計算采用標準的k－ε方程，模型中各通用常數為：Cμ＝0.09，C1＝1.44，C2＝1.92，σk＝1，σε＝1.33。

模擬時能量的松弛系數取0.8，體積分數的松弛系數取0.2，其他保持缺省的松弛系數值。

3 不同條件下數值模擬分析

3.1 水溫分布模擬

（1）不同入口冰花濃度模擬結果與分析。圖3所示為冰蓋上表面溫度為－23℃，入口流速0.8m／s，水深2m，不同入口冰花濃度條件下，距起點500m處斷面的水溫分布。入口冰花濃度不同則水流與冰蓋下面之間傳遞的熱量不同，因此流場的溫度存在差異，從圖3可以看出。入口冰花濃度越大，水流的溫度值越低。

圖3 不同入口冰花濃度斷面水溫分布

（2）不同冰蓋上表面溫度模擬結果與分析。圖4所示為入口冰花濃度0.05，入口流速0.8m／s，水深2m，不同冰蓋上表面溫度條件下，距起點500m處斷面的水溫分布。冰蓋上表面溫度不同則水流與壁面之間傳遞的熱量不同，因此水溫分布存在差異，從圖4可以看出，冰蓋上表面溫度越低，水溫越低。

圖4 不同冰蓋上表面溫度斷面水溫分布

（3）不同入口流速模擬結果與分析。圖5所示為入口冰花濃度0.05，冰蓋上表面溫度－23℃，水深2m，不同進口流速條件下，距起點500m處斷面的水溫分布，從圖5可以看出，不同流速時水溫分布存在差異，流速越小，水溫越低。

圖5 不同入口流速斷面的水溫分布

3.2 沿程冰蓋厚度模擬

（1）不同入口冰花濃度模擬結果與分析。圖6所示為冰蓋上表面溫度為－23℃，入口流速0.8m／s，水深2m，不同入口冰花濃度條件下沿程冰蓋厚度變化情況。入口冰花濃度影響傳遞熱量的大小，從圖6可以看出，入口冰花濃度越大，沿程冰蓋厚度越大，并且同一濃度條件下，沿程冰蓋厚度達到一定量值后變化不再明顯。

圖6 不同入口冰花濃度沿程冰厚

（2）不同冰蓋上表面溫度模擬結果與分析。圖7所示為入口冰花濃度0.05，入口流速0.8m／s，水深2m，不同冰蓋上表面溫度條件下沿程冰蓋厚度變化情況。冰蓋上表面溫度影響傳遞熱量的大小，從圖7可以看出，冰蓋上表面溫度越低，沿程冰蓋厚度越大。

圖7 不同冰蓋上表面溫度沿程冰厚

（3）不同入口流速模擬結果與分析。圖8所示為入口冰花濃度0.05，冰蓋上表面溫度－23℃，水深2m，不同進口流速條件下沿程的冰蓋厚度，在壁面導熱情況下，水流從入口處迅速凝結達到一定厚度。冰蓋的厚度沿程變化緩慢，從圖8可以看出，不同流速時沿程冰蓋厚度存在差異，流速越小，冰蓋厚度越大。

圖8 不同入口流速沿程冰厚

3.3 冰蓋下速度場的分布模擬

（1）不同入口冰花濃度模擬結果與分析。圖9a所示為冰蓋形成后距起點500m處斷面的流速分布。冰蓋形成后，過流斷面面積減小，斷面流速增大，從圖9a可以看出，入口冰花濃度越大，形成的冰蓋的厚度越大，斷面流速越大。

（2）不同冰蓋上表面溫度模擬結果與分析。圖9b所示為冰蓋形成后距起點500m處斷面流速分布。冰蓋形成后，過流斷面面積減小，斷面流速增大，從圖9b可以看出，冰蓋上表面溫度越低，形成的冰蓋的厚度越大，斷面流速越大。

（3）不同入口流速模擬結果與分析。圖9c所示為冰蓋形成后距起點500m處斷面的流速分布。冰蓋形成后，過流斷面面積減小，斷面流速增大，從圖9c可以看出，入口流速越小，形成的冰蓋的厚度越大，斷面流速變化越大。

圖11 冰蓋形成后斷面的流速分布

6 結 論

通過本文的研究得出了以下規律：冬季寒區河流冰蓋的厚度增長受入口冰花濃度、冰蓋上表面溫度、入口流速等因素的影響，入口冰花濃度越大、冰蓋上表面溫度越低、入口流速越小，則冰蓋下的水溫越低，形成的冰蓋厚度越大，斷面上的流速越大。

實際河流的結冰過程中，冰蓋與外界的熱交換還受到其他一些因素的影響，例如降水或降雪產生的熱交換，這些因素都可以在邊界條件的參數設置中做更深入的考慮。

［1］ Pegau W S，Page U S，Paulson C，Ronald J V.et al.Optical measurements of frazil concentration［J］.Cold Regions Science and Technology，1996，24：341－353.

［2］ Shen H T，Chiang L A .Simulation of growth and decay of river ice cover ［J］.Journal of Hydraulic Engineering，1984，110：958－971.

［3］ Sui Jueyi，Wang Jun，He Yun，et al.Velocity profiles and incipient motion of frazil particles under ice cover［J］.International Journal of Sediment Research，2010，25（1）：39－51.

［4］ Morse B，Richard M.A field study of suspended frazil ice particles［J］.Cold Regions Science and Technology，2009，55：86－102.

［5］ 王 軍.平衡冰塞輸冰的試驗研究［J］.水力發電學報，2002，76（1）：61－67.

［6］ 王 軍，陳胖胖，江 濤，等.冰蓋下冰塞堆積的數值模擬［J］.水利學報，2009，40（3）：348－354.

［7］ 楊開林，劉之平，李桂芬，等.河道冰塞的模擬［J］.水利水電技術，2002，33（10）：40－47.

［8］ Shen H T，Wang Desheng，Wasantha Lal A M.Numerical simulation of river ice processes ［J］.Journal of Cold Regions Engineering，1995，9（3）：107－118.

［9］ Kubat I，Sayed M，Savage S B，et al.Numerical simulations of ice thickness redistribution in the Gulf of St.Lawrence［J］.Cold Regions Science and Technology，2010，60：15－28.

［10］ Hicks F.An overview of river ice problems：CRIPE07 guest editorial［J］.Cold Regions Science and Technology，2009，55：175－185.

［11］ 王曉玲，周正印，蔣志勇，等.考慮氣溫變化影響的引水渠道水內冰演變數值模擬［J］.天津大學學報：自然科學版，2010，43（6）：516－522.