李民康，冀鴻蘭，羅紅春，郜國明，張寶森

(1.內蒙古農(nóng)業(yè)大學水利與土木建筑工程學院，內蒙古 呼和浩特 010018； 2.黃河水利科學研究院國際合作與科技局，河南 鄭州 450003；3.黃河水利科學研究院防汛搶險技術研究所，河南 鄭州 450003)

封凍期河冰演變過程包括成冰、流凌、初封、穩(wěn)封、融冰、開河過程，冰花的遷移貫穿在整個封凍期，是河冰變化過程中不可或缺的一部分[1]，冰在水體中隨水流流動稱為流凌，可分為封凍前結冰流凌和開河期解凍流凌，前者以冰花為主，后者以碎冰為主。流凌的輸移演變改變了河道邊界條件和水力特征，河道堆冰容易形成冰塞冰壩，在彎道險工處尤為明顯，冰塞冰壩引起過流斷面減小，壅高上游水位，易形成凌汛險情。在數(shù)值模擬分析領域，目前對流凌和水內冰的研究頗具代表性。茅澤育等[2-3]采用河道冰水耦合模型及垂向二維紊流模型，對天然河道冰塞、水內冰演變、河冰過程等進行了模擬；樊霖等[4]根據(jù)冰/水力學、熱力學建立數(shù)值模型模擬伊丹河封凍期冰情演變；王軍等[5]基于多相流理論，對流體相和顆粒相建立數(shù)學模型，得到冰顆粒的遷移規(guī)律；李超等[6]基于二維DynaRICE河冰模型對三湖河口彎道水力特性及卡冰過程進行數(shù)值模擬，并分析了彎道卡冰機理；Shen[7]利用SPH法模擬河冰二維運動，在RICE和DynaRICE模型的基礎上改進了河冰模擬系統(tǒng)CRISSP，可模擬冰塞演變及冰顆粒遷移等河冰過程；倪寶玉等[8]利用FLUENT模擬了波浪與冰的耦合作用，獲得了多浮冰耦合效應；金占軍[9]利用RNGk-ε湍流模型和流體體積法(Volume of fluid，VOF)方法，對冰水兩相流浮冰受力、冰水兩相流中浮冰的翻轉下潛進行了研究；冀鴻蘭等[10]基于可變模糊聚類循環(huán)迭代模型對流凌的輸移發(fā)展進行分析，預報了巴彥高勒站流凌日期；張瑩[11]利用PIV實驗研究與SPH數(shù)模結合的方法對單/雙顆粒流場特性進行了研究；劉寒秋等[12]通過歐拉-拉格朗日雙向耦合法，對冰晶影響下的管道沖蝕磨損進行了研究；鄧義斌等[13]研究了冰水兩相流在不同因素影響下的管道阻力特性，對冰水兩相流動的阻力特性進行數(shù)值模擬研究；王曉玲等[14]建立了三維非穩(wěn)態(tài)歐拉兩相流模型，模擬分析了不同轉彎半徑對彎道排冰的影響；姜坤[15]利用彎道水槽，基于FLUENT-EDEM耦合的三維兩相流數(shù)值模型，對冰顆粒的輸移特性進行了研究。

以上研究主要是利用各類方法對冰水兩相流的內在機制進行探討，或在特定約束條件(彎道)下進行冰水耦合運動的分析。然而，彎道幾何特性如曲率半徑、彎曲度、寬深比等，對水流結構、泥沙輸移及多相流耦合過程有著重要影響，其中，吳新宇等[16]以若爾蓋盆地實測資料為基礎，建立水動力數(shù)學模型，分析了彎曲河流頸口裁彎過程中的水流運動特性。然而在流凌條件下，彎道幾何特性的研究還需進一步深入。本文利用VOF技術，采用RNGk-ε湍流模型，對彎道幾何特性中的彎曲度進行建模分析，利用CAD構建4種典型對稱彎道，類型為大曲率比寬淺明流，底坡為平底坡。將流凌簡化為球體顆粒，主要研究封凍前期結冰流凌的遷移，以及冰水耦合作用下彎道水力特性，探究冰水耦合作用下彎曲度不同彎道的流凌分布、水位、流速、湍動能(TKE)等指標的基本變化規(guī)律。

1 數(shù)值模型

1.1 水流控制方程

基于笛卡爾坐標的水流連續(xù)性方程和動量方程為

(1)

(2)

式中：ui為第i個方向(i=1，2，3)的速度分量；xi為第i個方向的空間坐標；p為壓強；ρ為流體密度；Gi為第i個方向的重力加速度；τij為黏性應力；τbi為壁面剪應力，其在靠近實體邊界處激活?？墒褂糜邢摅w積法結合流體體積法(VOF)來求解水流控制方程，以跟蹤液體表面，并在液體表面應用適當?shù)膭恿吔鐥l件。還可利用廣義最小殘差求解器(GMRES)求解壓力和速度。

1.2 湍流模型

RNGk-ε湍流模型可考慮旋流及反映平均應變率的主流[17]，適應于高剪切和分離區(qū)域，尤其適用于顆粒流情況下復雜水流的模擬，因此采用了RNGk-ε湍流模型模擬冰水兩相流。

湍流模型由以下脈動能k和耗散率ε的輸運方程組成：

(3)

(4)

式中：vt為流體渦黏度系數(shù)；μ為流體動力黏滯系數(shù)；模型參數(shù)cε1=1.42，cε2=1.68，c3=0.012,η0=4.38，cμ=0.085，αk=αε=1.39。模型中R的存在使湍流對平均應變速率更敏感。

1.3 粒子的動力學控制方程

(5)

拉格朗日粒子的擴散根據(jù)Monte Carlo方法進行計算，粒子之間不會相互影響。采用基于經(jīng)驗公式導出的阻力模型，該模型將阻力系數(shù)和球體與其周圍流體的相對速度關聯(lián)以計算阻力[18]。

(6)

粒子-流體耦合模型將離散質量粒子的動量與連續(xù)流體隱式耦合，兩者的相互作用主要是動量交換，最終理論上顆粒的運動主要受重力、黏滯力、兩相拖曳力以及根據(jù)壓力梯度計算出的浮力等影響。描述每個相的方程式必須同時求解，兩相動量耦合表示為阻力系數(shù)與相間相對速度的乘積。兩相耦合基本方法：①根據(jù)未知的流體速度計算新的粒子速度；②將粒子速度代入流體動量方程中，并求解新的流體速度；③計算兩相阻力及粒子速度[19]。

而顆粒與邊界(壁邊界)亦會產(chǎn)生碰撞，存在能量損失，所以模擬考慮恢復系數(shù)。該系數(shù)將顆粒與壁面碰撞前后的速度分為切向和法向速度，反射后與反射前粒子的法向速度分量之比即為恢復系數(shù)r，當r=1.0時，則在反射過程中不會發(fā)生粒子動能損失。

1.4 FAVOR技術

自由表面邊界和流體界面使用流體體積法VOF處理：對任意網(wǎng)格點，若被流體占據(jù)則體積分數(shù)F=1，否則F=0。F是一個以0和1為邊界的連續(xù)函數(shù)，由對流傳輸方程控制：

(7)

模擬采用FAVOR技術(fractional area /volume obstacle representation)來定義計算域，可在歐拉網(wǎng)格內定義實體邊界，通過確定部分單元中開放面積和體積的分數(shù)模擬復雜形狀。定義邊界的過程獨立于網(wǎng)格進行，消除和順滑不平整區(qū)域，且可在特定區(qū)域做加密處理。

2 模型驗證

2.1 粒子輸移驗證

通過對流凌擬顆?；M行簡化，引入具有離散相的冰顆粒，在此參照文獻[20]，對引入顆粒的擴散性進行分析驗證。試驗內容為2 000個標記粒子初始條件下向x方向擴散，擴散系數(shù)為0.1 cm2/s，圖1中展示了顆粒濃度分布，以總粒子數(shù)作無量綱化處理，以理論高斯曲線分布作為對比?？梢钥闯?，理論曲線和模擬曲線吻合良好，拉格朗日粒子模型可有效描述粒子擴散。

圖1 不同時刻無量綱化點源粒子分布與理論高斯分布的比較

2.2 彎道水動力學模型驗證

以Shukry彎道水流試驗對CFD的水動力學模塊進行驗證[21]。Shukry矩形水槽彎道上游采取流量邊界，給定初始水位0.3 m，同時給定恒定流量0.072 m3/s，沿程使用壁邊界，下游使用壓力出口，控制出口水位0.28 m。該數(shù)值模擬激活重力、湍流模型及密度變化模型，采用RNGk-ε湍流模型，對數(shù)值運算步長及壓力迭代機制等進行參數(shù)控制，滿足收斂性條件，對彎道速度、水位的實測值與模擬值進行對比驗證，如圖2所示。由圖2可見，該模型平面流速與試驗等值線圖基本一致，最大流速在進口近凸岸處，在彎頂自凹岸向凸岸遞增，在出口近凹岸側。雖彎道速度最大值與試驗值略有出入，但能夠模擬出水動力軸線的遷移，較好地反映了實際水流水力特性。

圖2 Shukry彎道水位和流速對比

模擬與實測水位變化范圍均在0.20～0.33 cm之間，水位在彎道處呈現(xiàn)凹岸向凸岸遞減的趨勢，能夠模擬水位凹岸高、凸岸低的水位特征。因此，采用的水動力學模型可用于模擬彎道三維水流運動特性。

3 模型建立及邊界條件

參考周建銀[22]被廣為引用的系列水槽試驗，建立單一彎道模型，其曲率半徑與水槽模型一致，對入口、出口直道段及水槽高度進行適當延長加高，以滿足拉格朗日粒子的充分輸移。分別建立彎曲度不同的單彎模型，其中水槽彎曲段中心線半徑長為 8.53 m，彎曲度分別為60°、90°、150°和180°，中心線曲率比為3.645，彎道底坡為平底坡。兩個連接直道段長15 m，寬2.34 m，高1.2 m，初始水深0.8 m，模型尺寸如圖3所示。

上游均采用流量進口，設置流量q=1.872 m3/s，進口水位0.8 m，水流弗勞德數(shù)Fr=0.357，平均流速vav=1 m/s，Re=470 422.7。下游采用壓力出口，水位控制在0.65 m，流體經(jīng)過區(qū)域采用壁邊界，頂邊界為大氣壓。采用結構化矩形網(wǎng)格，最小網(wǎng)格尺寸為0.12 m×0.12 m，計算時間200 s。參考野外實測數(shù)據(jù)，流凌前期冰花尺寸在5 cm左右，拉格朗日粒子模型設置為相同尺寸的mass solid particles，其直徑為0.05 m，密度為917 kg/m3，與彎道邊界接觸后的恢復系數(shù)為0.5，阻力系數(shù)乘數(shù)設為1(意味計算域采用球體粒子充填)。水槽糙率設為0.015，靜摩擦系數(shù)為0.5。粒子設為隨水流運動，經(jīng)過驗算確定粒子源生成粒子速率為200個/秒。粒子和流體相互作用采用雙向動量耦合，并采用控制變量法將無冰顆粒的清水流作為對照試驗。

4 結果與分析

冰的密度小于水，主要分布于水體表層，對水流表層的水力特性影響較大，故以下主要從流凌分布、水面速度、橫向環(huán)流、水位及湍動能變化幾個方面探究流凌對彎道水流的影響。

4.1 彎道流凌分布

流凌在不同彎道輸移特點如圖4所示。由圖4可見，彎曲度不同的各彎道，在達到穩(wěn)定狀態(tài)后，流凌在進口段及進入彎道初期(進口斷面)的分布情況與直道分布相似。過彎道時，流凌前緣沿水流主軸線不斷向前發(fā)展，隨著粒子的集聚，逐漸向凹岸聚集，呈現(xiàn)凹岸多，凸岸少的特點，并沿凹岸輸出彎道。出彎后流凌沿水流中心線左右擺動，隨下游控制水位的影響，逐漸向兩岸集聚。

圖4 不同彎曲度彎道流凌輸移及分布

另外，各彎道流凌沿凹岸逐漸增多擴散，流凌分布多呈現(xiàn)楔形分布，這與李淑祎等[23]在試驗中的結論基本一致，因模擬無法對粒間相互作用力進行描述，該分布主要是流凌與槽間作用及彎道環(huán)流的影響所導致。各彎道速度較大流凌近凸岸，且隨著彎曲度的增大，其范圍逐漸縮小，高流速流凌數(shù)量沿彎道斷面依次遞減。流凌在凹岸的聚集，也解釋了冰塞現(xiàn)象易集中停滯在彎道處，尤其在凹岸及回水區(qū)低流速區(qū)。

4.2 彎道流速及橫向環(huán)流

4.2.1彎道流速

不同彎道冰水流和清水流表面流速等值線如圖5所示。

圖5 各彎道表面流速分布等值線

對于清水流，各彎道均在入口凸岸產(chǎn)生高流速區(qū)，出口凸岸產(chǎn)生低速區(qū)，這與Shukry彎道水槽試驗規(guī)律基本一致。隨著彎曲度的增大，主流越早偏向凸岸，彎道出口斷面主流區(qū)向凹岸偏移程度隨著彎曲度的增大不斷加大[24]。

冰水流表面流速基本規(guī)律與清水流基本一致，主要不同在凹岸流態(tài)及主流區(qū)的改變。流凌在凹岸形成長條狀低流速區(qū)，且沿水流該區(qū)域逐漸變寬。由上述冰凌運移規(guī)律可知，此處為流凌聚集處，聚集的流凌對凹岸流態(tài)產(chǎn)生了很大影響。隨著彎曲度的加大，低流速區(qū)發(fā)展越廣泛，冰花堆積越多。冰水流的主流亦因凹岸冰花的堵塞，向凸岸偏移程度變大，縮小主流過流面積。隨上游流凌的不斷堆積發(fā)展，此處極易形成冰塞體，并逐漸向上游遞進，為河段險工部位。綜上可見，由于流凌與水槽作用及彎道環(huán)流的存在，冰水流流速發(fā)展不充分，導致冰水流與清水流流速分布出現(xiàn)差異。

對各彎道冰水流及清水流流速進行定量分析，依據(jù)彎曲度的不同，對進口、彎頂、出口3個典型斷面沿凹岸到凸岸水面流速分布進行分析。

如圖6所示，60°彎道冰水流及清水流表面流速最大值均為v頂>v進>v出，出口斷面冰水流速最小值小于清水流。90°彎道冰水流及清水流表面流速最大值均為v進>v頂>v出，進口斷面兩類流動相似，出口斷面冰水流速最小值小于清水流，兩類流動凹岸向凸岸流速變化劇烈。150°彎道冰水流及清水流表面流速最大值v進>v頂>v出，進口及出口斷面兩類流動速度相似，彎頂斷面凹岸速度差異大，冰水流的速度更小。180°彎道進口及出口斷面，兩類流動表面流速分布基本一致，流速最大值v進>v出>v頂，凹岸向凸岸方向，彎頂斷面冰水流速最小值小于清水流。綜上可知，兩類流動進口斷面流速基本相等，彎頂及出口斷面流速變化強烈。各彎道最大進口流速均大于出口流速，主要是凹岸流速發(fā)生變化。60°和90°彎道冰水流和清水流進口與彎頂斷面水面流速基本一致，主要體現(xiàn)在出口斷面的不同，150°和180°彎道冰水流與清水流進口與出口斷面流速基本一致，主要體現(xiàn)在彎頂斷面流速的不同。這是由于隨彎曲度增大，凹岸低流速區(qū)離凸岸高流速區(qū)越遠，低流速集中區(qū)延后。隨著彎曲度增大，在彎頂斷面，冰水流沿凹岸到水流中心流速差呈增大趨勢，即冰水流在急彎河段彎頂部位，其低流速區(qū)范圍更廣泛。

圖6 各彎道典型斷面水面流速沿槽寬分布

如表1所示，除180°彎道，隨著彎曲度的增大，冰水流與清水流彎頂斷面凹岸流速最小值的差值呈增大趨勢，而出口斷面凹岸流速最小值的差值呈減少趨勢，說明彎曲度對典型斷面流速的影響逐漸由出口轉向彎頂斷面。

表1 不同彎道彎頂斷面及出口斷面凹岸流速變化 m/s

4.2.2 橫向環(huán)流分布

為分析流凌條件下水流垂向分布特征，對流凌發(fā)展充分且流速變化大的180°彎道的典型斷面進行分析，并與相同條件下清水流對比。不考慮縱向流速矢量，180°冰水流和清水流速度矢量圖如圖7所示?？梢园l(fā)現(xiàn)，在進口及出口斷面，兩類流動橫向環(huán)流速度矢量基本相似，其中進口斷面橫向環(huán)流呈現(xiàn)自凹岸向凸岸變化的特征，而出口斷面產(chǎn)生相似的對稱環(huán)流，冰水流凹岸下有小環(huán)流。對彎頂斷面分析可知，有流凌覆蓋的條件下，凹岸形成愈加強烈的冰下環(huán)流，與靠近凸岸的大環(huán)流是對稱關系，環(huán)流強度也小于后者。冰水流位于凸岸底部的大環(huán)流明顯小于清水流，可能是由于積聚的流凌影響水流過流能力，在一定程度上影響凸岸的大環(huán)流的發(fā)展演化。

圖7 不同彎道冰水流與清水流各斷面橫向環(huán)流分布

4.3 水位

為研究流凌對彎道水位的影響，對彎頂斷面橫向水位進行分析，由于離心力的存在而使自由水面的平衡狀態(tài)遭到破壞，進入彎段即有從凸岸到凹岸傾斜的橫比降Jr[25]，水深橫向變化如圖8和圖9所示。各彎道凹岸到凸岸水位都有降低趨勢，水面超高依次為0.03 m(60°彎道)、0.04 m(90°彎道)、0.06 m(150°彎道)和0.03 m(180°彎道)，說明彎曲度對寬淺明渠水深和水面超高有一定影響，但不明顯。彎頂斷面冰水流水位比清水流高，但隨著彎曲度的增大兩者差值變小。

圖8 不同彎道彎頂斷面冰水流與清水流水深變化對比

圖9 各彎道彎頂斷面的水深變化

4.4 彎道湍動能

從圖10中看出彎曲度為60°、90°、150°、180°的彎道，清水流和冰水流湍動強烈部位主要從彎道入口處發(fā)展，集中于近岸兩側。凹岸高湍動能區(qū)冰水流比清水流寬，湍動較強，凸岸側相似，凹岸向凸岸湍動能先減小后增大，呈U形分布。說明聚集的流凌冰在一定程度上增加了流體的湍動能，隨彎曲度增大，湍動能呈增大趨勢。

圖10 各彎道冰水流和清水流湍動能變化

5 結 論

a.過彎道時，流凌前緣沿水流主軸線不斷向前發(fā)展，隨著流凌的集聚，流凌逐漸向凹岸聚集，呈現(xiàn)凸岸少、凹岸多的特點，并逐漸向凹岸集聚輸出，彎道處流凌分布多呈現(xiàn)楔形分布。各彎道流速較大流凌顆粒近凸岸側，且隨著彎曲度的增大，其范圍逐漸縮小，高流速流凌數(shù)量分布呈現(xiàn)沿彎道斷面依次遞減。

b.各彎道最大進口流速均大于出口流速，隨著彎曲度的增大，彎頂斷面離凸岸高流速區(qū)越遠，低流速集中區(qū)較為延后，冰水流沿凹岸到水流中心流速差呈增大趨勢，即急彎河段彎頂部位滯冰區(qū)越典型、凹岸低流速區(qū)域越大。流凌的存在增大凹岸表層水流橫向環(huán)流，在一定程度上影響凸岸底部大環(huán)流的發(fā)展。

c.彎曲度對寬淺明渠水深和水面超高有一定影響，但不明顯。彎頂斷面水位冰水流比清水流高，但隨著彎曲度的增大兩者差值變小。冰水流與清水流相比，主要體現(xiàn)在凹岸湍動能變化，隨彎曲度的增大凹岸湍動能呈增大趨勢，說明聚集的流凌冰在一定程度上增加了流體的湍動能。