陳海文， 李捷， 徐立， 陳迪林， 陶鋮， 孫強(qiáng)

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

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

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

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

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

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

1.1 連續(xù)性方程

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

(1)

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

αs+αl=1

(2)

1.2 動(dòng)量守恒方程

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

(3)

(4)

冰晶顆粒動(dòng)量守恒：

(5)

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

1.3 能量守恒方程

(6)

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

1.4 湍流模型

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

(7)

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

(8)

式中：

(9)

1.5 相間作用力

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

Fdrag.l=ksl(vs-vl)

(10)

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

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

(11)

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

(12)

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

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

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

(13)

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

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

(14)

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

(15)

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

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

(16)

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

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

(17)

(18)

(19)

1.7 壁面振動(dòng)模型

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

(20)

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

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

2.1 物理模型

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

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

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

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

2.2 邊界和初始條件

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

2.3 模型驗(yàn)證

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

圖3 海水冰漿流動(dòng)實(shí)驗(yàn)臺(tái)實(shí)物照片

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

表3 實(shí)驗(yàn)測(cè)量裝置參數(shù)

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

圖5 實(shí)驗(yàn)制取的海水冰漿

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

圖6 不同振幅條件下流動(dòng)壓降的實(shí)驗(yàn)結(jié)果與模型計(jì)算結(jié)果對(duì)比

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

3.1 振動(dòng)對(duì)海水-冰晶顆粒分布影響

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

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

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

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

3.2 振動(dòng)對(duì)速度分布影響

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

圖9 振動(dòng)直管中心處流速變化和湍流動(dòng)能的變化

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

3.3 流動(dòng)阻力變化

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

(21)

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

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

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

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

圖11 不同振動(dòng)工況下壓降的變化情況

4 結(jié)論

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

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

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