陳海文， 李捷， 徐立， 陳迪林， 陶鋮， 孫強

(1.武漢理工大學(xué) 高性能船舶技術(shù)教育部重點實驗室， 湖北 武漢， 430063；2.武漢理工大學(xué) 能源與動力工程學(xué)院， 湖北 武漢， 430063))

船舶在極地航行時，換熱管在冷卻系統(tǒng)運行中具有重要作用，必須使其能夠高效穩(wěn)定的運行。由于海水中存在冰晶介質(zhì)，一定條件下產(chǎn)生堵塞現(xiàn)象[1]，使換熱器無法正常工作；船舶上機械設(shè)備的運行產(chǎn)生固有頻率在0～100 Hz振動[2]，使管道內(nèi)海水-冰晶兩相流的流動狀態(tài)更為復(fù)雜。因此，有必要綜合研究振動對直管中海水-冰晶兩相流流動特性的影響。

近年來，隨著振動強化傳動技術(shù)的發(fā)展，Sun等[3]通過CFD(計算流體動力學(xué))中的離散相模型、沖蝕磨損模型和動網(wǎng)格技術(shù)，研究不同振動工況下海水-冰晶兩相流對管道沖蝕磨損特性的影響。Xu[4]建立了CFD模型和數(shù)群平衡模型(PBM)相耦合方法。研究水平管道中冰漿流動，表明冰粒大小對冰漿流動特性有明顯影響。劉圣春等[5]通過以顆粒相動力學(xué)為基礎(chǔ)的歐拉-歐拉模型，研究不同管道模型內(nèi)冰漿的壓降特性，表明直管、彎管、T形管壓力隨流速增大而增加，且T形管的壓降最大。王繼紅等[6]通過CFD研究液固兩相流在水平直管中的流動特性，顆粒的密度和湍流強度對液固兩相流流動狀態(tài)起到關(guān)鍵作用。Bordet等[7]通過實驗研究了入口平均軸向速度和冰體積分數(shù)對冰漿流動模式和壓降的影響，在層流和湍流狀態(tài)下，體積冰濃度高達18.4%。徐立等[8]基于熱焓多孔介質(zhì)模型建立數(shù)學(xué)模型，研究流速對管道換熱產(chǎn)生的作用。Man等[9]對管換熱器中插入紐帶實驗研究，在低雷諾數(shù)下提高熱交換器的傳熱效率，但可能造成壓降增大。但目前對于海水-冰晶兩相流的流動特性研究較少，對添加振動的海水-冰晶兩相流還未涉及，振動對于管道的摩擦損失也不容忽略。

基于船舶航行的實際工況中振動復(fù)雜及不規(guī)則性，實際中實驗仿真不易，相關(guān)研究較少。本文使用Euler雙流體模型和UDF自定義動網(wǎng)格技術(shù)，研究不同振動下水平直管中海水-冰晶兩相流的流動特性，探討振動對冰晶顆粒分布及流動阻力的影響。

1 數(shù)學(xué)模型

海水-冰晶兩相流流動具有一定的復(fù)雜性，本研究通過假設(shè)冰晶為球形、光滑非彈性的顆粒簡化數(shù)學(xué)模型，不考慮流動過程中的粘性耗散[10]。冰晶顆粒在管道中的流動分布通過歐拉-歐拉雙流體模型描述，將海水和冰晶看作可以相互貫穿的連續(xù)性介質(zhì)，分別建立相應(yīng)的Navier-Stokes方程。

1.1 連續(xù)性方程

質(zhì)量守恒是流體力學(xué)中每一個流動介質(zhì)微元體必須遵循的基本規(guī)律，表達式為：

(1)

局部固相與液相的體積分數(shù)相關(guān)性為：

αs+αl=1

(2)

1.2 動量守恒方程

由于固液相的作用力不同，用相間動量傳遞表達動量守恒。海水動量守恒為：

(3)

(4)

冰晶顆粒動量守恒：

(5)

式中：g表示重力加速度；P表示靜壓；Ps表示冰晶顆粒相間的正應(yīng)力；μ表示剪切黏度；τ表示剪切應(yīng)力；λ表示粘滯系數(shù)；M表示海水與冰晶顆粒兩相間的動量交換；I表示單位張量。

1.3 能量守恒方程

(6)

式中：ui、vi、wi分別表示各相在x、y、z方向上的速度，m/s。ai為各項的擴散率，m2/s。

1.4 湍流模型

模型以海水-冰晶混合液為研究對象，考慮液體相和顆粒相湍流動能作用對流動過程的影響。海水-冰晶兩相流在管道中的流動特性用混合相k-ε湍流模型描述，并通過壁面函數(shù)法解決壁面流動問題。k-ε湍流模型[11]中的湍流動能方程為：

(7)

式中σk=1。湍流動能的耗散率方程可表示為：

(8)

式中：

(9)

1.5 相間作用力

溶液流動時，固相和液相之間的動量交換是由相間力引起，其相間作用力主要考慮湍流擴散力和拖曳力。其中湍流擴散力Ftd,l通過Burns等[12]的模型獲得；拖曳力表述為：

Fdrag.l=ksl(vs-vl)

(10)

式中ksl表示動量交換系數(shù):

ksl=3CDρIαsαl/4dS|vs-vl|

(11)

拖曳力系數(shù)CD，表示為：

(12)

1.6 顆粒流動力學(xué)理論

固體顆粒的性質(zhì)通過顆粒流動力學(xué)理論描述，其基本思想為把固體顆粒處理為致密氣體分子。在國內(nèi)外兩相流研究中，Gidaspow[13]提出的基于氣體分子運動論的顆粒動力學(xué)理論模型被廣泛使用，表達式為：

(-PsI+τs):vs+·(kθsθs)+δls-γθs

(13)

式中：(-PsI+τs):vs是粒子應(yīng)力張量；kθs是擴散系數(shù)；kθsθs為擴散能；δls為相間轉(zhuǎn)換能量；γθs是顆粒碰撞所消耗的能量。

式(13)中，顆粒溫度θs用于描述顆粒波動能量，其表達式為：

(14)

式(13)中，顆粒壓力Ps的方程表示為:

(15)

式中：ess==0.9為粒子碰撞恢復(fù)系數(shù)；go是顆粒徑向分布函數(shù)；δls=-3kslθs。

體積粘度ξs表征顆粒抵抗變形的能力，表示為：

(16)

顆粒剪切粘度由粒間碰撞引起的顆粒粘性μs.c及動力粘度μs.k組成:

μs=μs.c+μs.k

(17)

(18)

(19)

1.7 壁面振動模型

如圖1為水平管內(nèi)冰晶流動示意圖，沿x軸為方向流動，在管壁外壁面上添加振動。振動具有簡諧運動規(guī)律，將振動相位角設(shè)為0，則圓管振動過程位移、速度及加速度分別表示為[14]:

(20)

圖1 水平直管內(nèi)冰晶流動示意

2 模型設(shè)置及驗證

2.1 物理模型

物理模型如圖2所示，選取長度為1 000 mm，外徑為25 mm，內(nèi)24 mm的水平直管。數(shù)值模擬水平直管相關(guān)參數(shù)如表1所示，材料為碳鋼管，密度為7 850 kg/m3,采用ANSYS ICEM對直管進行網(wǎng)格劃分，管道兩端采用O-Grid網(wǎng)格，將直管劃分至 75 411個六面體網(wǎng)格，數(shù)值計算結(jié)果符合網(wǎng)格無關(guān)性要求。通過C++語言編譯振動條件UDF自定義動網(wǎng)格程序，將動網(wǎng)格用于壁面進行求解。海水-冰晶熱物性參數(shù)與海水鹽度有關(guān)，本研究海水鹽度取為15‰。根據(jù)在北極實地航行的科學(xué)考察船“永盛輪”提供的資料和殼管式換熱器設(shè)計規(guī)范[15-16]，海水-冰晶熱物性參數(shù)如表2所示。

圖2 直管網(wǎng)格劃分

表1 水平直管相關(guān)參數(shù)

表2 海水和冰晶的熱物性參數(shù)

2.2 邊界和初始條件

為了研究海水-冰晶兩相流流動特性，本研究忽略直徑較小的冰晶顆粒，顆粒直徑統(tǒng)一為0.5 mm。在管道邊界條件方面，計算域入口流體采用速度入口，由于振動形成非穩(wěn)定流場，造成出口處壓強和速度改變，因此采用自由流出條件。管壁面為無滑移速度，近壁面處采用標準壁面函數(shù)，忽略顆粒在流動過程中的損耗，管壁壁面熱流量為0，絕熱情況。采用相耦合SIMPLE算法及二階迎風(fēng)精度求解離散方程，將時間步長設(shè)為0.02 s，計算殘差小于1.0×10-4時為收斂狀態(tài)。

2.3 模型驗證

為了驗證所建立的CFD數(shù)值模型的有效性，搭建了海水冰漿流動實驗平臺，如圖3所示為其實物圖，相對應(yīng)的海水冰漿制取流動實驗測試系統(tǒng)原理如圖4。實驗系統(tǒng)包括：冰漿的制取裝置、冰漿的儲存裝置、冰漿輸送系統(tǒng)、流動實驗測試段。本實驗通過對整個管道系統(tǒng)外壁包裹2 cm厚度的聚氨酯材料進行保溫絕熱，直管內(nèi)海水冰漿流動在絕熱條件下測試。測試段為長度1 m，管徑為24 mm的直管段，通過無紙記錄儀及差壓變送器讀取壓降值，實驗測量裝置參數(shù)如表3。

圖3 海水冰漿流動實驗臺實物照片

注：1.海水制冰機， 2.儲冰保溫槽， 3.冰漿泵， 4.泄水閥， 5.盲板閥， 6.電磁流量計， 7.溫度傳感器(a：測量入口溫度 b：測量出口溫度)， 8.振動調(diào)節(jié)器， 9.測試部分， 10.差壓變送器， 11.可視部分， 12.球閥(a、b、c：流量調(diào)節(jié)閥；d：取樣閥)， 13.排氣閥。

表3 實驗測量裝置參數(shù)

本實驗保證冰漿運行溫度在-1.7 ℃左右，冰漿樣品如圖5所示。冰粒直徑通過制冰機制取0.5 mm范圍內(nèi)，含冰率IPF控制約為15%，流速控制在1.5 m/s。

圖5 實驗制取的海水冰漿

實驗主要測試振動頻率為20 Hz條件下，通過振動調(diào)節(jié)器改變振幅大小，測量冰漿流動過程中壓降的結(jié)果，并與CFD雙歐拉模型模擬結(jié)果相比較，如圖6所示。圖6可以看出，實驗結(jié)果與仿真計算對于壓降影響的變化趨于一致，且兩者的最大相對誤差在20%范圍內(nèi)，這由于數(shù)值模擬忽略粘性耗散與摩擦效應(yīng)，未考慮實驗中存在的相變等因素，在實驗中隨著振幅的增大，顆粒與海水的能量交換越多。由此，數(shù)值模擬方法適用于研究海水-冰晶兩相流流動特性。

圖6 不同振幅條件下流動壓降的實驗結(jié)果與模型計算結(jié)果對比

3 數(shù)值模擬結(jié)果與討論

3.1 振動對海水-冰晶顆粒分布影響

直管出口處不同振動工況下冰晶顆粒體積分布如圖7所示，冰晶體積分數(shù)分布范圍0.03～0.21，沿管道圓形截面垂直方向上，冰晶體積分數(shù)呈對稱性分布。冰晶體積分布出現(xiàn)明顯分層，冰晶顆粒多聚集在上壁面，這主要因為冰晶顆粒密度小于海水密度，在浮升力的作用下顆粒群出現(xiàn)上浮現(xiàn)象。圖7表示振幅為0.5 mm下，上壁面冰晶顆粒體積分布明顯減少，由于振動頻率增加，導(dǎo)致冰晶顆粒之間碰撞頻繁，致使冰晶顆粒受到的浮升力小于冰晶顆粒碰撞產(chǎn)生的力。

圖7 出口處冰晶顆粒體積分布(A=0.5 mm)

圖8在頻率為20 Hz下，由于振幅的增大導(dǎo)致冰晶顆粒振動勢能增加，顆粒與壁面碰撞增強，冰晶體積分布逐漸減少。

圖8 出口處冰晶顆粒體積分布(ω=20 Hz)

3.2 振動對速度分布影響

如圖9(a)研究振動對軸線上冰晶顆粒速度分布，流速分布先增加再減小，最終趨于平緩。主要因為進口段擾動較大，冰粒與流體間相互作用增強，能量轉(zhuǎn)換增多使得壓降增加，造成流速迅速增加。隨著冰粒與壁面的磨損作用及冰粒間的碰撞損耗，冰粒速度逐漸降低，在0.5 m處趨于平緩。與靜態(tài)相比，隨著振幅與頻率的添加，流速有微小的減弱。主要因為振動對中心區(qū)域流速影響較小，流速基本趨于穩(wěn)定。如圖9(b)，沿直管軸線截面湍流動能先逐漸增加，0.5 m處湍流動能達到較大值，主要由于入口段壓降較大，顆粒與海水相互作用增強，導(dǎo)致軸向拖曳力增加。由于管內(nèi)粘性摩擦與機械摩擦作用，流動損耗增加，湍流動能逐漸降低至穩(wěn)定值。

圖9 振動直管中心處流速變化和湍流動能的變化

圖10為不同振動工況下振動對出口截面y軸顆粒流速分布系數(shù)的影響，從近壁面到主流區(qū)流速逐漸增加且主流區(qū)流速分布系數(shù)基本不變。這由于靠近壁面粘性摩擦與機械摩擦作用，同時冰粒在上壁面聚集，導(dǎo)致顆粒間的相互作用力阻礙顆粒的流動。隨著振動頻率及振幅的增加，近壁面流速分布系數(shù)逐漸增大。這由于頻率增加導(dǎo)致顆粒與壁面之間碰撞頻繁，振幅的增加促使顆粒與壁面碰撞作用增強，導(dǎo)致機械振動能量轉(zhuǎn)換為動能的量增大，顆粒的速度增加。

3.3 流動阻力變化

海水-冰晶流動過程中，由于固液流體之間的相互作用力及進口處擾動，流體流動阻力不容忽略。流體的流動阻力采用摩擦系數(shù)f表示，摩擦系數(shù)由流體在管內(nèi)流動過程中壓力降表示為[17-18]：

(21)

式中：D和L分別是直管直徑和長度；ρc是冰晶顆粒密度；u為流體流速。

流體流動阻力f與流道壓降P成正比關(guān)系，流動阻力的變化可通過分析壓力降得出。

圖10 不同振動工況下振動對出口處流速分布系數(shù)的影響

圖11為不同振動工況下壓降變化趨勢。隨著振動頻率及振幅的增加，直管進出口的壓降逐漸增大。由于振幅和頻率的增大，顆粒的機械振動能量增加，造成顆粒與壁面及顆粒間的相互作用增大，致使更多的振動勢能轉(zhuǎn)換為顆粒動能，通過式(21)壓降與速度平方成正比，壓降增加。

圖11 不同振動工況下壓降的變化情況

4 結(jié)論

1)冰晶兩相流在管道中流動由于浮升力作用產(chǎn)生明顯分層現(xiàn)象，上壁面附近冰晶顆粒積聚。隨著頻率從0增到50 Hz、振幅從0增加到1.2 mm，管壁附近區(qū)域流體的擾動增強，造成上壁面附近的冰晶顆粒體積分布逐漸減少。

2)由于顆粒與壁面間粘性摩擦與機械摩擦作用，壁面附近速度小于主流區(qū)速度。主流區(qū)流速及湍流動能受振動影響較小，在近壁面振動增強引起振動能量增加，轉(zhuǎn)換成動能量增大，流速與湍流動能增大。

3)通過數(shù)值模擬與實驗結(jié)果先比較，驗證了數(shù)值模擬的有效性。在頻率為20 Hz條件下，振幅從0增加到1.2 mm，海水冰晶兩相流的流速增大，造成管道進出口壓降增大。