鄭成龍， 陳胖胖， 王 軍

（合肥工業(yè)大學(xué) 土 木與水利工程學(xué)院，安徽 合 肥 230009）

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

對(duì)河冰過(guò)程研究的方法主要有原型觀(guān)測(cè)、試驗(yàn)研究和數(shù)值模擬等［2－5］。數(shù)值模擬近年來(lái)在國(guó)內(nèi)外發(fā)展很快。文獻(xiàn)［6］對(duì)冰蓋下冰塞堆積進(jìn)行了模擬；文獻(xiàn)［7］對(duì)河道的冰塞進(jìn)行了模擬；文獻(xiàn)［8］對(duì)江河的河冰演變進(jìn)行了模擬；文獻(xiàn)［9］對(duì)Lawrence河灣冰厚分布進(jìn)行了模擬。

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

本文基于多相流和傳熱學(xué)理論，把河流當(dāng)作大氣－冰體－水體的一個(gè)耦合系統(tǒng)來(lái)處理，考慮系統(tǒng)中各界面之間的熱量交換，并耦合冰花在冰蓋下堆積對(duì)冰蓋熱量交換的影響，通過(guò)數(shù)值模擬，分析了入口冰花濃度、入口流速、冰蓋上表面溫度等因素對(duì)冰蓋厚度增長(zhǎng)的影響。

1 控制方程

連續(xù)性方程：

其中，t為時(shí)間；αk為第k相的體積分?jǐn)?shù)；ρk為第k相的密度；uk為第k相的平均速度矢量；i為冰相，l為水相；水的密度ρl為997kg／m3，冰的密度ρi為917kg／m3。

動(dòng)量方程：

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

溫度方程：

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

邊界條件如下。

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

湍流動(dòng)能率方程為：

其中，Cμ為經(jīng)驗(yàn)常數(shù)，取Cμ＝0.09；Lm為特征長(zhǎng)度。

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

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

初始條件如下：

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

2 模型建立及網(wǎng)格劃分

建立長(zhǎng)1 000m，水深2m的矩形河流二維模型，采用結(jié)構(gòu)網(wǎng)格對(duì)模型進(jìn)行網(wǎng)格劃分，將網(wǎng)格劃分為0.1m×0.1m的正方形網(wǎng)格，共10 000×20＝20×104個(gè)網(wǎng)格。用有限體積法離散連續(xù)性方程、動(dòng)量方程、標(biāo)準(zhǔn)k－ε方程等，在近壁面處采用壁面定律以考慮邊壁糙率的影響，采用Simple算法。計(jì)算模型及模型局部網(wǎng)格劃分，如圖1、圖2所示。

圖1 計(jì)算模型

圖2 計(jì)算模型局部網(wǎng)格劃分

計(jì)算過(guò)程中紊流的計(jì)算采用標(biāo)準(zhǔn)的k－ε方程，模型中各通用常數(shù)為：Cμ＝0.09，C1＝1.44，C2＝1.92，σk＝1，σε＝1.33。

模擬時(shí)能量的松弛系數(shù)取0.8，體積分?jǐn)?shù)的松弛系數(shù)取0.2，其他保持缺省的松弛系數(shù)值。

3 不同條件下數(shù)值模擬分析

3.1 水溫分布模擬

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

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

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

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

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

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

3.2 沿程冰蓋厚度模擬

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

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

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

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

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

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

3.3 冰蓋下速度場(chǎng)的分布模擬

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

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

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

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

6 結(jié) 論

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

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

［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］ 王 軍.平衡冰塞輸冰的試驗(yàn)研究［J］.水力發(fā)電學(xué)報(bào)，2002，76（1）：61－67.

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

［7］ 楊開(kāi)林，劉之平，李桂芬，等.河道冰塞的模擬［J］.水利水電技術(shù)，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］ 王曉玲，周正印，蔣志勇，等.考慮氣溫變化影響的引水渠道水內(nèi)冰演變數(shù)值模擬［J］.天津大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，43（6）：516－522.