劉獻(xiàn)飛, 王恒，王方, 李志強(qiáng), 朱彩霞, 張浩飛，王仕元
中原工學(xué)院能源與環(huán)境學(xué)院

       【摘  要】單螺桿膨脹機(jī)螺旋槽道內(nèi)存在多種流型，掌握槽道內(nèi)兩相流動(dòng)特性膜厚分布對(duì)減小氣體泄漏，保證膨脹機(jī)高效穩(wěn)定運(yùn)行至關(guān)重要。本文采用VOF兩相流數(shù)值模型，對(duì)垂直放置螺旋槽道內(nèi)兩相流動(dòng)特性及膜厚分布進(jìn)行了計(jì)算流體動(dòng)力學(xué)分析。根據(jù)前期實(shí)驗(yàn)獲得的螺旋槽內(nèi)兩相環(huán)狀流的實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證了數(shù)值模型的準(zhǔn)確性，揭示了不同液相流體入射角度螺旋槽道內(nèi)液相工質(zhì)演變過(guò)程，理論分析不同液相流體入射角度下螺旋槽道外側(cè)的液膜厚度分布規(guī)律，發(fā)現(xiàn)液相流體入射角度為0.9倍的螺旋升角時(shí)，液膜厚度分布最均勻。

       【關(guān)鍵詞】螺旋槽道; 兩相流; VOF模型; 流動(dòng)特征；液膜厚度

       【基金項(xiàng)目】國(guó)家自然科學(xué)基金項(xiàng)目（51906265）；河南省高校科技創(chuàng)新人才支持計(jì)劃（20HASTIT019）

引    言

       單螺桿膨脹機(jī)具有在低溫下效率高以及對(duì)于工質(zhì)的狀態(tài)不敏感的優(yōu)點(diǎn)，高壓氣體、飽和蒸汽、汽液兩相或熱液都可作為膨脹工質(zhì)，大大拓寬了低溫?zé)嵩蠢玫臏囟确秶４送猓瑔温輻U膨脹機(jī)所適用的功率范圍在1kW～1000kW 范圍內(nèi)，克服了傳統(tǒng)蒸汽輪機(jī)和燃?xì)廨啓C(jī)不能做小的缺陷。因此，單螺桿膨脹機(jī)在工業(yè)節(jié)能、可再生能源利用、建筑節(jié)能和熱泵系統(tǒng)等眾多領(lǐng)域具有非常廣闊的應(yīng)用前景^[1]。目前關(guān)于單螺桿膨脹機(jī)的研究開發(fā)得到了國(guó)內(nèi)外學(xué)者的高度重視，但由于單螺旋膨脹機(jī)涉及熱力學(xué)、傳熱學(xué)、流體力學(xué)、機(jī)械設(shè)計(jì)、材料學(xué)及控制等多學(xué)科，理論研究難度很大。從國(guó)內(nèi)外研究現(xiàn)狀看，目前仍然側(cè)重于膨脹機(jī)的加工和性能研究^[2,3]，對(duì)膨脹機(jī)螺旋槽道泄露間隙的液膜流動(dòng)及分布特性研究數(shù)據(jù)缺乏、理論不足，而這正是制約膨脹機(jī)效率提高的關(guān)鍵因素^[4]。

       螺旋管內(nèi)兩相流動(dòng)特性相對(duì)較小。Murai等人^[5]說(shuō)明了空氣-水兩相流通過(guò)螺旋管時(shí)，離心加速度對(duì)流動(dòng)結(jié)構(gòu)分布的影響。Zhu等人^[6]基于雙傳感器電導(dǎo)探針研究了離心力對(duì)相分布的影響。Zhang等人^[7] 研究了螺旋管內(nèi)油水分離機(jī)理，指出分離效率與旋轉(zhuǎn)半徑、壓降、相密度差和液滴直徑有關(guān)。Colombo等人^[8]采用歐拉模型對(duì)螺旋管內(nèi)氣液兩相混合流動(dòng)進(jìn)行了數(shù)值模擬，得到了離心力對(duì)流場(chǎng)和相分布的影響。Mota和Pagano^[9]研究了不同流動(dòng)參數(shù)條件和不同管道幾何結(jié)構(gòu)下螺旋管內(nèi)氣相分布。Jayakumar等人^[10]研究了螺旋管內(nèi)兩相流動(dòng)時(shí)，螺旋管直徑和螺距對(duì)熱工水力特性的影響。Vashith和Nigam^[11]研究了兩相流通過(guò)螺旋管時(shí)的速度分布、平均摩擦系數(shù)、相分布和入口長(zhǎng)度。

       我們前期對(duì)螺旋矩形通道內(nèi)環(huán)狀流的流型演變和無(wú)量綱液膜厚度進(jìn)行了實(shí)驗(yàn)研究^[12]，并且基于均相流模型^[13-15]，研究了不同結(jié)構(gòu)參數(shù)對(duì)兩相壓降和相分布的影響，給出了螺旋通道內(nèi)氣液兩相流體受離心力作用后相分布和液膜厚度。考慮到單螺桿膨脹機(jī)潤(rùn)滑油由油泵驅(qū)動(dòng)，單獨(dú)流入膨脹機(jī)間隙，用于潤(rùn)滑和密封^[1]。而目前對(duì)于潤(rùn)滑油單獨(dú)進(jìn)入間隙后，在螺旋矩形通道內(nèi)的相演化行為、界面現(xiàn)象和膜厚分布的研究還很缺乏，本文利用VOF模型，分別定義氣相和液相流動(dòng)參數(shù)條件，給出了膨脹機(jī)螺旋槽道內(nèi)液相演變規(guī)律及液膜厚度分布，這對(duì)如何調(diào)控膨脹機(jī)復(fù)雜截面螺旋槽道內(nèi)液相流動(dòng)、相界面特征及液膜厚度，避免高壓氣體向低壓區(qū)域的漏泄損失，保證螺桿膨脹機(jī)高效、穩(wěn)定運(yùn)行具有至關(guān)重要的理論指導(dǎo)意義。 

1. 數(shù)值方法

       1.1 螺旋槽道物理結(jié)構(gòu)示意圖

       圖1給出了螺旋槽道及其主要幾何參數(shù)的示意圖，矩形截面的寬度和長(zhǎng)度分別為a=18mm、b=30mm。螺旋直徑（從通道中心到軸測(cè)量）用D=195mm表示，兩個(gè)相鄰匝間的距離稱為節(jié)距H=306mm，圖中顯示半節(jié)距h。螺旋線圈一圈的投影與垂直于軸線的平面形成的角度稱為螺旋角β[^16]。流體入射角度（η）是指定義的流動(dòng)方向與法向邊界之間的夾角。螺旋槽道壁面最靠近線圈軸的一側(cè)稱為內(nèi)側(cè)，反之稱為外側(cè)。轉(zhuǎn)角（θ）螺旋通道任何位置的局部截面與入口截面平面之間的角度。

圖1 螺旋槽道及其主要幾何參數(shù)的示意圖

       1.2 數(shù)值方法

       采用基于流體體積法的三維CFD模型求解控制方程組，采用幾何重建方法對(duì)氣液界面進(jìn)行重建。該方案使用分段線性方法表示流體之間的界面，使得兩相流的計(jì)算在氣體和液體之間具有清晰的界面^[12]。采用有限體積法求解控制守恒和邊界條件，進(jìn)行了數(shù)值模擬。前期我們已經(jīng)使用了不同的湍流模型來(lái)預(yù)測(cè)螺旋管內(nèi)流動(dòng)的流體力學(xué)和熱特性，為了更精確地預(yù)測(cè)螺旋管內(nèi)兩相流的流體力學(xué)特性[14]和提高了螺旋管內(nèi)流動(dòng)的精度，這里選擇了重整化群k-ε湍流模型快速應(yīng)變流^[17]。

       連續(xù)性方程:

           (1)

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

           (2)

       式中，P、g、F、μ分別為流場(chǎng)中的壓力、重力加速度、作用于系統(tǒng)的物體力和流動(dòng)流體的粘度。

       上面給出的方程依賴于所有相的體積分?jǐn)?shù)，由性質(zhì)ρ和μ隱式表示。在兩相系統(tǒng)中，每個(gè)單元的密度和粘度可以估計(jì)為：

           (3)
           (4)

       通過(guò)求解體積分?jǐn)?shù)的標(biāo)量輸運(yùn)方程，可以追蹤兩相之間的界面：

           (5)

       對(duì)于液膜流動(dòng)，表面張力起著重要的作用。FLUENT中的表面張力模型是Brackbill等人提出的連續(xù)表面力（CSF）模型 ^[18]。在該模型中，表面張力對(duì)VOF計(jì)算的貢獻(xiàn)是動(dòng)量方程的源項(xiàng)，由表面張力產(chǎn)生的體積力可計(jì)算為^[19]：

           (6)

       式中，σ是表面張力系數(shù)，k是根據(jù)單位法向散度定義的自由表面曲率：

           (7)

       1.3 數(shù)值計(jì)算

       本文采用的求解程序是基于商用CFD軟件FLUENT，網(wǎng)格是使用商業(yè)軟件GAMBIT創(chuàng)建的。圖2描繪了用于三維計(jì)算域的網(wǎng)格（六面體網(wǎng)格和楔形網(wǎng)格）。在FLUENT求解過(guò)程中使用非穩(wěn)態(tài)求解器，整個(gè)過(guò)程物理性質(zhì)保持一致的，操作環(huán)境考慮重力作用。邊界條件是入口分別為氣相和液相速度入口，出口采用壓力出口邊界，出口方向變量的擴(kuò)散通量設(shè)為零。在槽道壁面設(shè)置了一個(gè)穩(wěn)定、無(wú)滑移、無(wú)滲透的邊界，并考慮了強(qiáng)化壁面作用處理。壓力插值采用壓力交錯(cuò)選項(xiàng)（PRESTO）格式，壓力-速度耦合采用壓力隱式算子分裂（PISO）格式，動(dòng)量方程采用二階迎風(fēng)差分格式，體積分?jǐn)?shù)方程采用Geo重構(gòu)格式，動(dòng)量、湍流動(dòng)能和湍流耗散率采用二階迎風(fēng)格式。基于兩相流的動(dòng)態(tài)特性，進(jìn)行了時(shí)間步長(zhǎng)為0.0001s的非穩(wěn)態(tài)數(shù)值模擬，連續(xù)性、速度、k和epsilon的收斂準(zhǔn)則為1.0e-5。

       為了更好的獲得螺旋槽道內(nèi)的兩相流特征，首先對(duì)各網(wǎng)點(diǎn)精度進(jìn)了獨(dú)立性驗(yàn)證。數(shù)值模擬的網(wǎng)格測(cè)試考慮了75104、152478、237162和291403個(gè)網(wǎng)格的計(jì)算節(jié)點(diǎn)，網(wǎng)格間距分別為0.0015、0.0012、0.001和0.0095。使用不同網(wǎng)格的預(yù)測(cè)值如表1所示，結(jié)果表明，0.001的區(qū)間節(jié)點(diǎn)比例可以為研究矩形螺旋槽內(nèi)的液相占有率和液膜厚度提供可接受的統(tǒng)計(jì)精度。

（a）    （b）
圖 2 螺旋槽道的網(wǎng)格特征（a）計(jì)算域（b）入口截面
表1 不同網(wǎng)格數(shù)下的預(yù)測(cè)結(jié)果比較

       1.4 數(shù)據(jù)處理

       圖3a示出了從商業(yè)軟件FLUENT 6.3的輸出獲得的橫截面含液率分布的原始圖像，為了更好地測(cè)量液膜厚度，采用數(shù)字圖像處理算法將液膜輪廓與其他流動(dòng)特征分離開來(lái)。通過(guò)設(shè)置含液率閾值，得到了兩相流結(jié)構(gòu)中液相的等表面圖，在本文研究中液相的閾值固定為0.9^[20]。圖3b示出了測(cè)量螺旋槽道外側(cè)的液膜厚度的示意圖，為了更好地從圖片信息中獲取液膜厚度，本文采用了50倍的放大因子進(jìn)行測(cè)量。

（a）                                           （b）
圖 3（a） 截面局部含氣率分布；（b）從氣相分離液膜輪廓后液膜厚度的測(cè)量

2. 結(jié)果與討論

       2.1 模型驗(yàn)證

       VOF方法已廣泛應(yīng)用于界面位置的跟蹤和界面現(xiàn)象的求解，并進(jìn)行了大量的工作^[8,11,17]。對(duì)于螺旋矩形通道內(nèi)的兩相流，關(guān)于液膜厚度的研究成果并不多見(jiàn)。數(shù)值模型的驗(yàn)證是根據(jù)我們以前的研究^[12]中環(huán)形流的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行的。由于離心力的影響，在螺旋矩形通道的外側(cè)突出了環(huán)形兩相流的液膜厚度。為了驗(yàn)證數(shù)值解的準(zhǔn)確性，根據(jù)實(shí)驗(yàn)研究，進(jìn)行了該數(shù)值模擬的結(jié)構(gòu)參數(shù)和兩相質(zhì)量流量的計(jì)算。圖4給出了實(shí)驗(yàn)和數(shù)值模擬所得液膜厚度的比較，圖中還顯示了實(shí)驗(yàn)數(shù)據(jù)的不確定度。實(shí)驗(yàn)數(shù)據(jù)與模擬數(shù)據(jù)的偏差在3.79%以內(nèi)，因此，本研究建議計(jì)算模型有足夠的精確度。

圖 4 實(shí)驗(yàn)和數(shù)值模擬的平均液膜厚度對(duì)比

       2.2 界面特征的演變規(guī)律

       單螺桿膨脹機(jī)槽道外側(cè)的液相分布對(duì)防止氣體泄漏具有重要作用，圖5給出了螺旋槽道不同橫截面處的含液率分布及演變規(guī)律。值得注意的是，當(dāng)η=0時(shí)，由于螺旋槽道的節(jié)距產(chǎn)的扭曲和流體流動(dòng)慣性的影響，流動(dòng)會(huì)發(fā)生扭曲，在y/a=0位置時(shí)可以看到較大的持液率。隨著流動(dòng)的進(jìn)一步發(fā)展，液體被離心力作用移向外側(cè)。同時(shí)，在初始發(fā)展階段，液體被節(jié)距產(chǎn)生的扭轉(zhuǎn)力作用開始向y/a=1方向移動(dòng)。在θ=60o角處，液體分布在水平中心線上接近對(duì)稱，然后，由于二次流運(yùn)動(dòng)的慣性作用，液相的繼續(xù)向y/a=1方向移動(dòng)使得液相分布不對(duì)稱^[13]，這種現(xiàn)象主要是由于螺旋通道入口的扭轉(zhuǎn)效應(yīng)引起的。隨著流動(dòng)的發(fā)展，扭轉(zhuǎn)效應(yīng)弱于曲率效應(yīng)，對(duì)液體脈動(dòng)分布的影響減弱。當(dāng)η=0.9β時(shí)，根據(jù)數(shù)值模擬結(jié)果可以發(fā)現(xiàn)，扭轉(zhuǎn)效應(yīng)對(duì)螺旋槽液相分布的影響基本可以消除，主要受到離心力的作用，螺旋槽的外側(cè)保持較高的持液率。

圖5 螺旋槽道不同橫截面處的含液率分布及演變規(guī)律（a）η=0；（b）η=0.9β

       2.3 液相入射角度對(duì)液膜厚度分布的影響

       從不同η下的等高線變化可以看出η對(duì)液相分布具有明顯的影響，為了進(jìn)一步揭示η對(duì)液相分布的影響，圖6給出了不同η下的液膜厚度變化規(guī)律。從圖6a可以看出，當(dāng)η=0.80β時(shí)，液膜厚度在y/a=0位置明顯較高，因?yàn)?，這是由于η較小，流動(dòng)仍然受到扭轉(zhuǎn)力的明顯影響。隨著流動(dòng)的不斷發(fā)展，在y/a=0位置的液膜厚度先增大后減小，這是螺旋槽道節(jié)距產(chǎn)的扭曲和流體流動(dòng)中的慣性，使得液體更容易在這一區(qū)域積聚。隨著扭轉(zhuǎn)力作用將流體從y/a=0的位置移到y(tǒng)/a=1的位置，使得y/a=1位置的液膜厚度逐漸增大。當(dāng)η=0.85β時(shí)，由圖6b可以看出，液膜厚度在y方向的不對(duì)稱性顯著減小。當(dāng)η=0.90β時(shí)，圖6c表示液膜厚度分布基本與均勻一致。當(dāng)η繼續(xù)增加到0.95β時(shí)，如圖6d所示，在y/a=1區(qū)域出現(xiàn)較厚的膜厚分布。盡管η的值小于β，但由于螺旋槽道產(chǎn)生的較高的二次流導(dǎo)致液體向y/a=1移動(dòng)^[21]。當(dāng)η=β時(shí)，如圖6e表示，y/a=1位置的液膜厚度明顯高于y/a=0位置。對(duì)比不同θ下測(cè)量數(shù)據(jù)之間的薄膜厚度值可以發(fā)現(xiàn)，由于離心力始終在薄膜表面的法向上，隨著流動(dòng)的不斷發(fā)展，外側(cè)壁面薄膜平均厚度逐漸增加^[22]。

(a) η=0.80β                                            (b) η=0.85β

(c) η=0.90β

 (d) η=0.95β                                                      (e) η=β
圖 6 不同進(jìn)口角下液體厚度分布的變化

3. 結(jié)論

       本文借助VOF多相流模型，數(shù)值分析了兩相流體的螺旋槽道內(nèi)的流動(dòng)特性，揭示了氣液兩相界面特征演變規(guī)律及液膜厚度分布規(guī)律，理論分析了液相入射角度對(duì)液膜厚度的影響規(guī)律，主要結(jié)論如下：

       1. VOF多相流模型能夠很好地預(yù)測(cè)螺旋槽道內(nèi)液相分布及液膜厚度變化規(guī)律，數(shù)值模擬結(jié)果和實(shí)驗(yàn)測(cè)量結(jié)果具有很好的一致性。

       2. 當(dāng)液相入射角度為η=0時(shí)，螺旋槽道內(nèi)液相分布呈現(xiàn)明顯的波動(dòng)性，當(dāng)η=0.90β時(shí)，扭轉(zhuǎn)效應(yīng)對(duì)螺旋槽液相分布的影響基本可以忽略不計(jì)。 

       3. 不同入射角度下的液膜厚度對(duì)比可以發(fā)現(xiàn)，當(dāng)η=0.9β時(shí)，螺旋槽道外側(cè)壁面可以獲得比較均勻的液膜厚度分布，當(dāng)η=0.9β時(shí)，y/a=0位置區(qū)域的液膜厚度較高，反之y/a=1位置區(qū)域的液膜厚度較高。

參考文獻(xiàn)

       [1] Wang Wei, Wu Yuting, Ma Chongfang et al, Preliminary experimental study of single screw expander prototype, Appl Therm Eng 2011;31:3684–3688. 
       [2] Ziviani D, Van M Den Broek, Paepe M De, Geometry-Based Modeling of Single Screw Expander for Organic Rankine Cycle Systems in Low-Grade Heat Recovery, Energy Procedia 2014;61:100–103.
       [3] Li Guoqiang, Lei Biao, Wu Yuting et al, Influence of inlet pressure and rotational speed on the performance of high pressure single screw expander prototype, Energy 2018;147:279–285 
       [4] Shen Lili, Wang Wei, Wu Yuting et al, Theoretical and experimental analyses of the internal leakage in single-screw expanders, Int. J. Refrigeration. 86 (2018) 273–281. 
       [5] Murai Y, Yoshikawa S, Toda S, Structure of air-water two-phase flow in helically coiled tubes, Nucl Eng Des 2006;236:94–106. 
       [6] Zhu Hongye, Li Zhaoxu, Yang Xingtuan et al, Flow regime identification for upward two-phase flow in helically coiled tubes, Chem Eng J 2017;308:606–618.
       [7] Zhang Jun, Guo Jun, Gong Daotong et al, An investigation on oil-water separation mechanism inside helical pipes, J Hydrodyn Ser B 2006;18:343– 347.
       [8] Inzoli F, M Colombo, Cammi A et al, CFD study of an air–water flow inside helically coiled pipes, Prog Nucl Energy 2015;85:462–472. 
       [9] da Mota FRM, Pagano DJ, Simulation and experimental study of phase segregation in helical pipes: A new method for flow conditioning, Flow Meas Instrum 2014;35:99–108. 
       [10] Jayakumar JS, Mahajani SM, Mandal JC et al, Thermal hydraulic characteristics of air-water two-phase flows in helical pipes, Chem Eng Res Des 2009;88:501–512. 
       [11] Vashisth S, Nigam KDP, Prediction of flow profiles and interfacial phenomena for two-phase flow in coiled tubes, Chem Eng Process Process Intensif 2009;48:452–463. 
       [12] Liu Xianfei, Xia Guodong, Yang Guang, Experimental study on the characteristics of air–water two-phase flow in vertical helical rectangular channel, Int J Multiph Flow 2015;73:227–237.
       [13] Xia Guodong, Liu Xianfei, Zhai Yuling et al, Single-phase and two-phase flows through helical rectangular channels in single screw expander prototype , J Hydrodyn Ser B 2014;26:114–121. 
       [14] Liu Xianfei, Xia Guodong, Zhai Yuling, Numerical analysis of the two-phase pressure drop and liquid distribution in single-screw expander prototype, Sci. bull 2014;59:4388–4396.
       [15] Xia Guodong, Liu Xianfei, An investigation of two-phase flow pressure drop in helical rectangular channel, Int Commun Heat Mass Transf 2014;54:33–41. 
       [16] Lu Xing, Du Xuepeng, Zeng Min et al Shell-side thermal-hydraulic performances of multilayer spiral-wound heat exchangers under different wall thermal boundary conditions, Appl Therm Eng 2014; 70: 1216–1227. 
       [17] Liu Xianfei, Wang Fang, Li Zhiqiang et al, Parametric investigation of thermal-hydrodynamic performance in the innovative helically coiled heat exchangers in the heat pump system, Energ Buildings, 2020; 216:109961
       [18] Shine SR, Kumar SS, Suresh BN, Numerical study of wave disturbance in liquid cooling film, Propuls Power Res 2013;2:107–118.
       [19] Gu Fang, Liu Chunjiang, Yuan Xiguang et al, CFD simulation of liquid film flow on inclined plates, Chem Eng Technol 2004; 27:1099–1104. 
       [20] Yu Jiawen, Jiang Yiqiang, Cai Weihua et al. Forced convective condensation flow and heat transfer characteristics of hydrocarbon mixtures refrigerant in helically coiled tubes, Int J Heat Mass Transf 2018;124:646–654.
       [21] Zhao Houjian, Li Xiaowei, Wu Xinin, Numerical investigation of supercritical water turbulent flow and heat transfer characteristics in vertical helical tubes, J Supercrit Fluids 2017;127:48–61.
       [22] Kang Zhongtao, Li Xiandong, Mao Xiaobing, Acta Astronautica Experimental investigation on the surface wave characteristics of conical liquid film, Acta Astronaut 2018;149:15––24.

       備注：本文收錄于《建筑環(huán)境與能源》2020年10月刊總第37期（第22屆全國(guó)暖通空調(diào)制冷學(xué)術(shù)年會(huì)文集）。版權(quán)歸論文作者所有，任何形式轉(zhuǎn)載請(qǐng)聯(lián)系作者。