(1 天津商業(yè)大學(xué) 天津市制冷技術(shù)重點(diǎn)實(shí)驗(yàn)室 天津 300134； 2 北京理工大學(xué) 北京 100081)

冰漿是一種水、冰晶粒子以及凝固抑制劑(比如氯化鈉、乙醇等)組成絮狀混合物，其中冰粒直徑≤1×10-3m，一般冰粒子濃度大于5%時(shí)認(rèn)為是非牛頓流體。冰漿作蓄冷介質(zhì)具有良好的熱物性、傳輸性及相變特性[1]，因此在建筑制冷、礦井降溫、食品冷藏、醫(yī)療衛(wèi)生和消防運(yùn)輸?shù)阮I(lǐng)域都有應(yīng)用[2-6]。近年來國(guó)內(nèi)外許多研究者對(duì)冰漿的流動(dòng)特性進(jìn)行了研究，但主要是結(jié)合冰漿流動(dòng)實(shí)驗(yàn)現(xiàn)象與非牛頓流體流變特性建立冰漿阻力特性預(yù)測(cè)模型。常用的非牛頓流體流變模型[7]有Bingham、Herschel-Bulkley、Powder Law和Casson，一些研究者針對(duì)這4種流變模型研究冰漿水平管道流動(dòng)壓降模型，M. Grozdek等[8]實(shí)驗(yàn)分析了在水平直管段內(nèi)含冰率0%～30%條件下的壓降特性。F. Illán等[9]以9%的NaCl溶液作為載流流體，分析了直管內(nèi)的傳熱和壓降特性。V. Ayel等[10]分析了冰漿的流動(dòng)和傳熱特性，其中流動(dòng)特性主要研究了直管內(nèi)的壓降特性。A. Kitanovski等[11]同樣分析了冰漿的流動(dòng)特性。然而，冰漿流動(dòng)壓降與管徑、管型、流速、冰粒子濃度等有關(guān)，這些模型不能兼顧所有因素對(duì)流動(dòng)壓降的影響，而且由于模型采用的勻質(zhì)流假設(shè)，只能用于低速分層流動(dòng)，通用性和精度不高。

隨著計(jì)算機(jī)技術(shù)和計(jì)算流體力學(xué)的發(fā)展，固液兩相流的數(shù)值模擬也隨之有了很大發(fā)展。常用的固液兩相模型有以離散單元法(DEM)為基礎(chǔ)的Euler-Langrian模型[12]和以顆粒動(dòng)力學(xué)為基礎(chǔ)的Euler-Euler模型[13]。Euler-Euler模型將顆粒相作為連續(xù)相，相對(duì)于Euler-Langrian模型計(jì)算量大為減小。將冰漿看作固液兩相流體，利用Ansys Fluent研究冰漿的流動(dòng)壓降特性，但此方面研究尚不多見[14]。B. Niezgoda-Zelasko等[15]提出利用兩相雙流模型研究冰漿阻力特性，但冰粒子黏度模型仍將冰漿等效為均質(zhì)流體，未考慮冰粒子的非均勻分布對(duì)黏度的影響，因此數(shù)學(xué)描述與實(shí)際流動(dòng)存在一定的偏差。本文采用以顆粒相動(dòng)力學(xué)為基礎(chǔ)的Euler-Euler模型模擬冰漿在水平直管、90°彎管和T型管中的無相變流動(dòng)壓降，并用實(shí)驗(yàn)進(jìn)行模型驗(yàn)證。

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

描述冰漿流體等溫流動(dòng)的雙流體計(jì)算模型[13]，數(shù)學(xué)描述可表示為：

1)連續(xù)性方程

(1)

由質(zhì)量守恒方程，得:

(2)

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

將冰漿看做液固混合的顆粒相流體，用液相動(dòng)力參數(shù)描述動(dòng)量方程[16]：

αqρqg+αq(Fq+Flift,q+FVm,q)+

(3)

固體顆粒相動(dòng)量方程為：

αsρsg+αs(Fs+Flift,s+FVm,s)+

(4)

液體相與顆粒相之間的作用力可用以下Gidaspow模型來[13]確定。

當(dāng)αl>0.8時(shí)，液體-固體交換系數(shù)Ksl有如下形式：

(5)

其中：

(6)

當(dāng)αl≤0.8時(shí)，Ksl有如下形式：

(7)

式中：Re為相對(duì)雷諾數(shù)。

主相q和第二相p的相對(duì)雷諾數(shù)：

(8)

3)湍流方程

本文采用k-ε湍流模型，該模型以N-S方程組為基礎(chǔ)，同時(shí)考慮顆粒相和液相的相互作用[16]。

(9)

(10)

其中：

(11)

Gk,m=μt,m(um+

(12)

μt,m=ρmCμk2/ε

(13)

(14)

式中：C1ε=1.44；C2ε=1.92；σε=1.3；σk=1；Cμ=0.09;下標(biāo)m為混合相；k為湍流脈動(dòng)動(dòng)能；μt,m為混合相湍流動(dòng)力黏度， Pa·s；σk為湍流脈動(dòng)動(dòng)能的普朗特?cái)?shù)；ε為耗散項(xiàng)。

4)黏度模型

冰漿是一種固液兩相混合物，物理密度、動(dòng)力黏度以及結(jié)冰點(diǎn)不僅受添加物濃度影響，還與含冰率有關(guān)，D.G.Thomas[17]通過反復(fù)實(shí)驗(yàn)和理論推算，得出黏度關(guān)于含冰率的公式：

μ=μl(1+2.5IPF+10.05IPF2+

0.002 73(e16.6IPF-1))

(15)

式中：IPF為含冰率。

2 數(shù)值模擬結(jié)果分析

借助模擬軟件FLUENT進(jìn)行流體計(jì)算，設(shè)定含冰率為10%和20%，冰晶粒子直徑約為1×10-4m，管徑為0.02 m；雙精度并考慮重力影響，流速為0.2～10 m/s；邊界條件在入口處設(shè)為速度入口，出口設(shè)為出流，壁面設(shè)為無滑移；求解器采用一階迎風(fēng)格式的SIMPLE算法；不考慮相變模型但開啟能量源項(xiàng)，收斂精度設(shè)置為10-4，冰漿的基本參數(shù)見表1。

表1 冰漿的基本參數(shù)

1)水平直管

冰漿通過水平直管時(shí)，管道縱截面平均壓力如圖1所示，可以看出冰漿在水平直管中流動(dòng)時(shí)壓力一直下降，且入口的壓降明顯比出口大，造成這種現(xiàn)象的原因是入口處擾動(dòng)較大，能量損失嚴(yán)重。

圖1 水平直管縱截面平均壓力分布云圖Fig.1 The contour of pressure on the longitudinal section of horizontal tube

2)90°彎管

圖2所示為90°彎管縱截面平均壓力的分布云圖?？梢钥闯?，拐角處的壓力變化劇烈，壓力損失嚴(yán)重，這是由于在拐彎處存在摩擦和二次流現(xiàn)象[18]，引起能量的損失。由于入口處擾動(dòng)的影響，使入口壓降大于出口壓降。

圖2 90°彎管縱截面平均壓力分布云圖Fig.2 The contour of pressure on the longitudinal section of the elbow tube

圖3所示為彎管彎曲處冰粒子的速度矢量圖，由此看出在彎管的彎曲處存在明顯的二次流現(xiàn)象。當(dāng)流體運(yùn)動(dòng)時(shí)遇到變向阻礙，彎道外側(cè)近壁面邊界層速度減小，離心力減小，而中心層速度較高，提供向心力的壓力也較高，因此造成中心層向外層流動(dòng)的趨勢(shì)，形成附加流動(dòng)，與此同時(shí)，還形成外側(cè)向內(nèi)側(cè)的流動(dòng)以保持質(zhì)量守恒。二次流是主流與螺旋流動(dòng)的疊加，因此能量損失極大。

圖3 90°彎管的拐彎處管道橫截面速度矢量云圖Fig.3 The contour of velocity on the cross section of elbow tube

圖4所示為90°彎管的管道縱截面冰漿速度分布云圖，可以看出在流體轉(zhuǎn)彎之前，固體冰粒速度逐漸增大，且有明顯的邊界層。在彎管彎曲處，有一個(gè)明顯的高速漩渦出現(xiàn)，這是因?yàn)樵陔x心力作用下流體向管道外側(cè)運(yùn)動(dòng)，并逐漸外側(cè)堆積，進(jìn)而出現(xiàn)外側(cè)向內(nèi)側(cè)的流動(dòng)。流體離開轉(zhuǎn)彎后，彎管外側(cè)的速度較大，而內(nèi)側(cè)速度較小。

3)T型管

T型管形狀較上兩種都復(fù)雜，圖5所示為T型管縱截面壓力分布云圖，左側(cè)為冰漿進(jìn)口，流體在三叉處分流至兩個(gè)出口，其中上側(cè)為出口1，右側(cè)為出口2。

圖4 90°彎管的管道縱截面冰漿速度分布云圖Fig.4 The contour of velocity for the ice crystals on the longitudinal section of the tubes

圖5 T型管縱截面壓力分布云圖Fig.5 The contour of pressure on the longitudinal section of T-type tube

由圖5可以看出，沿著冰漿流動(dòng)方向壓力變化趨勢(shì)分明，從流動(dòng)進(jìn)口到分流之前，壓力逐漸下降。在拐角處分流至兩個(gè)方向，一個(gè)方向是垂直于來流方向，另一個(gè)方向?yàn)檠刂俣确较蛳蚯?。拐角分流處液體壓力劇烈變化，可以看出冰漿流經(jīng)此處時(shí)有劇烈的湍流現(xiàn)象和能量交換，經(jīng)過動(dòng)量交換，冰漿在沿著來流方向的分流方向上壓力劇增，然后又逐漸下降；而在垂直于來流方向的分流方向上，壓力基本保持不變。值得注意的是垂直分流方向上的兩個(gè)拐角出現(xiàn)一個(gè)壓力最高點(diǎn)和最低點(diǎn)，最低點(diǎn)的出現(xiàn)是由于邊界層脫落分離，造成壓力突然降低；最高點(diǎn)的出現(xiàn)與90°彎管外側(cè)的情況相似，都是由于能量交換使速度減小而壓力增大。

圖6所示為T型管縱截面冰漿平均速度分布云圖，冰漿的平均速度除了在入口處有明顯的變化外，在分流前基本保持不變，在分流處急劇下降至一定值，兩個(gè)分流管內(nèi)的流體平均速度基本相等。

圖6 T型管縱截面冰漿平均速度分布云圖Fig.6 The contour of average velocity for ice slurry

3 實(shí)驗(yàn)方法及過程

參考H. Kumano等[19-20]的實(shí)驗(yàn)方案設(shè)計(jì)如圖7所示，實(shí)驗(yàn)系統(tǒng)包括4部分：冰漿存儲(chǔ)裝置、冰漿輸送系統(tǒng)、流動(dòng)實(shí)驗(yàn)測(cè)試段、數(shù)據(jù)采集系統(tǒng)。其中，冰漿存儲(chǔ)裝置是兩個(gè)尺寸為0.395 m×0.395 m×0.025 m的不銹鋼箱，外部包有10 mm厚的保溫材料。冰漿輸送系統(tǒng)由冰漿輸送泵、過濾器、流量調(diào)節(jié)裝置和輸送管構(gòu)成，冰漿輸送泵采用可用于漿體的磁力驅(qū)動(dòng)循環(huán)泵，過濾裝置用于篩除粒徑大于1×10-3m的冰粒子。流動(dòng)實(shí)驗(yàn)段需要測(cè)量溫度、壓降和流量等物理量，實(shí)驗(yàn)測(cè)量裝置的參數(shù)見表2。數(shù)據(jù)采集主要由GP10便攜式無紙記錄儀記錄。測(cè)量段的3種管型采用的是管徑為2×10-2m的黃銅管，管總長(zhǎng)為1 m，且彎管和T型管均在中間位置有改變，彎管曲率半徑為0.1 m，測(cè)試管水平放置，外包裹厚度約為1×10-2m的聚氨酯材料進(jìn)行保溫。

圖7 流動(dòng)實(shí)驗(yàn)段原理Fig.7 Principle of the flow experiments

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

實(shí)驗(yàn)前，先對(duì)儀表進(jìn)行標(biāo)定，溫度校準(zhǔn)設(shè)備為OPTI-CAL溫濕度校驗(yàn)儀，精度±0.1 ℃；壓力變送器經(jīng)過Beamex多功能校驗(yàn)儀標(biāo)定，INT160壓力模塊的量程為0～16 MPa，分辨率為0.000 1，精度為0.005%FS+0.0125%RDG。

實(shí)驗(yàn)時(shí)，儲(chǔ)存箱1中的冰漿溫度在-2 ℃左右，冰粒直徑約為4×10-4～8×10-4m。IPF為10%～20%，流速在0.3～1.5 m/s選定5組。

冰漿的含冰率IPF用密度測(cè)量法計(jì)算為：

(16)

式中：ρslurry為溶液平均密度(質(zhì)量與體積的比值，體積用量筒測(cè)量)，kg/m3；ρice、ρwater分別為冰和水的密度，kg/m3。

4 實(shí)驗(yàn)結(jié)果分析

1)管型對(duì)壓降的影響

管道形狀分為直管、90°彎管和T型管，3者之間對(duì)比分析的意義在于研究阻力部位的形狀和阻力程度的關(guān)系。

圖8 不同含冰率下3種管型內(nèi)的壓降隨流速變化Fig.8 Variation trend of pressure drop in three piping shapes

由圖8可以看出，壓降受管型影響較大。含冰率為10%和20%，壓降值從大到小依次為：T型管>彎管>直管。流速較低時(shí)，3種管道的壓降值較接近，相差不大。流速越高，由于管型阻力不同，造成同流速下T型管與90°彎管、90°彎管與直管間壓降差值變大。直管中沒有阻力部件，因此壓降主要是沿程損失，且直管部分的曲線可看出速度越高，壓力損失越大。彎管中有一個(gè)彎曲部位，流動(dòng)方向改變時(shí)，發(fā)生二次流，冰漿會(huì)撞擊外側(cè)壁面造成動(dòng)量的損失，在豎直部分流動(dòng)時(shí)，外側(cè)流速大，質(zhì)量多，故摩擦劇烈，壓力損失較大。T型管分流后一支路沿原速度方向，另一支路轉(zhuǎn)向90°，冰漿流經(jīng)分流處遇到岔口，會(huì)形成漩渦，擾動(dòng)更大。因此，在冰漿輸送中，應(yīng)盡量使用直管和彎管，T型管壓降過大，盡量減少使用。

圖9 直管壓力損失系數(shù)與流速的關(guān)系Fig.9 Pressure lose coefficients in horizontal piping

流體流動(dòng)時(shí)會(huì)產(chǎn)生沿程阻力，為了評(píng)價(jià)速度對(duì)壓力損失的影響，參考流體力學(xué)中沿程損失和文獻(xiàn)[21]，研究壓力損失系數(shù)：

(17)

圖9所示為直管壓力損失系數(shù)與流速的關(guān)系，可以看出，水平直管中流速越大，壓力損失系數(shù)越小，即對(duì)壓降的影響越小。而流速與雷諾數(shù)成正比，且密度和黏度變化不大時(shí)，壓力損失系數(shù)與雷諾數(shù)的關(guān)系與圖9相同。

2)含冰率IPF對(duì)壓降的影響

如圖8所示，3種管型，IPF為20%的壓降比IPF為10%的大。這是由于IPF越大，冰漿的非牛頓流體特征越明顯，黏度也隨固體顆粒的增多而增大。3種管型在高流速時(shí)的壓降差高于低流速。這是因?yàn)樗俣刃?，冰水分層現(xiàn)象較為明顯；隨著速度增大，冰水混合逐漸均勻，顆粒效應(yīng)使得黏度變大，因此壓降差增大。當(dāng)流速為0.9 m/s時(shí)，IPF從10%增大到20%時(shí)，直管壓降從1.1 kPa/m增大到1.5 kPa/m； 90°彎管壓降從1.9 kPa/m增大到2.1 kPa/m； T型管壓降從2.27 kPa/m增大到2.73 kPa/m。由圖8還可以看出，當(dāng)速度增大到一定值時(shí)，彎管和T型管壓降差值趨于恒定，原因是阻力部件存在，使得冰水摻混劇烈，易達(dá)到充分混合狀態(tài)，此時(shí)，顆粒效應(yīng)穩(wěn)定，黏度差恒定，壓降差基本不變。

3)實(shí)驗(yàn)與模擬結(jié)果對(duì)比

圖10 壓降模擬值與實(shí)驗(yàn)值對(duì)比Fig.10 Comparison of simulation results and experimental data for the pressure drop

圖10所示為水平直管中實(shí)驗(yàn)數(shù)據(jù)與模擬數(shù)據(jù)的對(duì)比。由圖可知數(shù)值模擬和實(shí)驗(yàn)的壓降隨流速的變化趨勢(shì)一致，但在流速較大時(shí)，模擬值要大于實(shí)驗(yàn)值且誤差達(dá)到20%，這是由于流速越大，冰漿流動(dòng)的湍流程度就越大，冰顆粒與水的相互碰撞越劇烈，二者之間能量交換越多，相變影響越大，而數(shù)值模擬不考慮相變，所以模擬值與實(shí)驗(yàn)值偏差較大。

5 結(jié)論

本文從模擬和實(shí)驗(yàn)兩個(gè)角度對(duì)冰漿在帶有阻力部件的管道中流動(dòng)特性進(jìn)行了研究。管道的形狀有直管、90°彎管和T型管3種，其中直管不含局部阻力部件，彎管帶有一個(gè)彎曲處，T型管的連接部位是三通，帶有兩個(gè)彎曲處。主要研究了壓降受含冰率、流速、管型的影響，結(jié)論如下：

1)從數(shù)值模擬的結(jié)果看，不考慮相變的影響時(shí)，3種管型在入口處都有明顯的壓降。對(duì)于彎管，在拐彎處出現(xiàn)二次流，壓力變化較大；而T型管在三通處壓力變化明顯，甚至在垂直分流方向上的兩個(gè)拐角分別出現(xiàn)壓力最高點(diǎn)和最低點(diǎn)。

2)當(dāng)實(shí)驗(yàn)條件和含冰率相同時(shí)，直管的壓降最小，90°彎管次之，而T型管的壓降最為嚴(yán)重。因此，在實(shí)際應(yīng)用過程中應(yīng)盡量避免T型管的使用。

3)冰漿各管型中水平方向流動(dòng)時(shí)，壓降都隨著含冰率和流速的增大而增大。

4)水平直管中計(jì)算的沿程壓力損失系數(shù)隨速度的增大而降低。

5)直管的流動(dòng)模擬結(jié)果與實(shí)驗(yàn)值誤差在20%以內(nèi)，在低流速時(shí)吻合較好。

本文受天津市應(yīng)用基礎(chǔ)與前沿技術(shù)研究計(jì)劃項(xiàng)目(15JCYBJC21600)和天津市高等學(xué)校創(chuàng)新團(tuán)隊(duì)項(xiàng)目(TD12-5048)資助。(The project was supported by the Tianjin Research Program of Application Foundation and Advanced Technology (No.15JCYBJC21600) and the Science Research Innovation Team Project in Tianjin (No.TD12-5048).)

[1] EGOLF P W, KAUFFELD M. From physical properties of ice slurries to industrial ice slurry applications[J]. International Journal of Refrigeration, 2005, 28(1):4-12.

[2] DAVIES T W. Slurry ice as a heat transfer fluid with a large number of application domains[J]. International Journal of Refrigeration, 2005, 28(1):108-114.

[3] KAUFFELD M, WANG M J, GOLDSTEIN V, et al. Ice slurry applications[J]. International Journal of Refrigeration, 2010, 33(8):1491.

[4] BELLAS I, TASSOU S A. Present and future applications of ice slurries[J]. International Journal of Refrigeration, 2005, 28(1):115-121.

[5] 杜衛(wèi)新, 鄭光相. 冰漿制冷降溫系統(tǒng)在平煤六礦的應(yīng)用[J]. 煤炭工程, 2009, 365(4):63-65. (DU Weixin, ZHENG Guangxiang. The application of the ice slurry system for six coal mine in Pingdingshan Coal Group Co. Ltd.[J]. Coal Engineering, 2009, 365(4):63-65.)

[6] OTAKE H, SHITE J, PAREDES O L, et al. Catheter-based transcoronary myocardial hypothermia attenuates arrhythmia and myocardial necrosis in pigs with acute myocardial infarction[J]. Journal of the American College of Cardiology, 2007, 49(2):250-260.

[7] 何國(guó)庚, 王忠衡. 冰漿流體流動(dòng)與換熱研究綜述[J]. 制冷學(xué)報(bào), 2005, 26(4):1-5. (HE Guogeng, WANG Zhongheng. Review of study on flow and heat transfer of ice slurry[J]. Journal of Refrigeration, 2005, 26(4):1-5.)

[8] GROZDEK M, KHODABANDEH R, LUNDQVIST P. Experimental investigation of ice slurry flow pressure drop in horizontal tubes[J]. International Journal of Refrigeration, 2010, 33(2):357-370.

[10] AYEL V, LOTTIN O, PEERHOSSAINI H. Rheology, flow behaviour and heat transfer of ice slurries: a review of the state of the art[J]. International Journal of Refrigeration, 2003, 26(1):95-107.

[11] KITANOVSKI A, VUARNOZ D, ATA-CAESAR D, et al. The fluid dynamics of ice slurry[J]. International Journal of Refrigeration, 2005, 28(1):37-50.

[12] KUO H P, KNIGHT P C, PARKER D J, et al. The influence of DEM simulation parameters on the particle behaviour in a V-mixer[J]. Chemical Engineering Science, 2002, 57(17):3621-3638.

[13] GIDASPOW D. Multiphase flow and fluidization[J]. Multiphase Flow & Fluidization, 1994(95):1-29.

[14] MONTEIRO A C S, BANSAL P K. Pressure drop characteristics and rheological modeling of ice slurry flow in pipes[J]. International Journal of Refrigeration, 2010, 33(8):1523-1532.

[15] NIEZGODA-ZELASKO B, ZALEWSKI W. Momentum transfer of ice slurry flows in tubes, experimental investigations[J]. International Journal of Refrigeration, 2006, 29(3):418-428.

[16] 王繼紅, 王樹剛, 張騰飛, 等. 水平管道內(nèi)冰漿流體阻力特性CFD模擬[J]. 大連理工大學(xué)學(xué)報(bào), 2012, 52(6):845-849. (WANG Jihong, WANG Shugang, ZHANG Tengfei, et al. Simulation of pressure drop for ice slurry flow in horizontal pipes by CFD[J]. Journal of Dalian University of Technology, 2012, 52(6):845-849.)

[17] THOMAS D G. Transport characteristics of suspension: VIII. A note on the viscosity of Newtonian suspensions of uniform spherical particles[J]. Journal of Colloid Science, 1965, 20(3):267-277.

[18] 王繼紅, 王樹剛, 張騰飛, 等. 水平90°彎管內(nèi)冰漿流體流動(dòng)特性的數(shù)值模擬[J]. 高校化學(xué)工程學(xué)報(bào), 2012, 26(4):581-586. (WANG Jihong, WANG Shugang, ZHANG Tengfei, et al. Numerical simulation of ice slurry flow in a horizontal 90° elbow pipe[J]. Journal of Chemical Engineering of Chinese Universities, 2012, 26(4):581-586.)

[19] KUMANO H, HIRATA T, SHIRAKAWA M, et al. Flow characteristics of ice slurry in narrow tubes[J]. International Journal of Refrigeration, 2010, 33(8): 1513-1522.

[20] MIKA. Ice slurry flow in a poppet-type flow control valve[J]. Experimental Thermal & Fluid Science, 2013, 45:128-135.

[21] MIKA. Pressure loss coefficients of ice slurry in horizontally installed flow dividers[J]. Experimental Thermal and Fluid Science, 2013, 45(2): 249-258.