李佳川， 魏英杰， 王聰， 夏維學

(哈爾濱工業(yè)大學 航天學院， 黑龍江 哈爾濱 150001)

0 引言

減小液體中運動物體的阻力是流體動力學中的一個重要問題。較低的阻力意味著更快的速度和更少的能源消耗。減小阻力的方法有很多，例如應用聚合物添加劑[1]、超空泡[2-3]、微氣泡注射[4-5]、疏水涂層[6-8]等。這些方法主要是通過改變物體的表面特性或通過改變物體附近流體的密度來實現(xiàn)減阻功能。

當物體的表面溫度遠高于液體沸點時，周圍液體汽化后的蒸汽不斷匯聚，形成連續(xù)的蒸汽膜，稱為膜態(tài)沸騰過程。由于蒸汽膜的絕熱性會使熱量交換變得很差，物體壁面會達到很高的溫度，甚至有燒毀的危險，在熱工程領(lǐng)域中，一般要避免膜態(tài)沸騰的發(fā)生。最近的研究表明，利用膜態(tài)沸騰形成的蒸汽膜能夠?qū)崿F(xiàn)減阻功能[9-13]，因此膜態(tài)沸騰減阻成為一種較新的減阻方法。

膜態(tài)沸騰產(chǎn)生的蒸汽膜具有滑移效應，類似于嵌入超疏水性壁面的空氣墊。然而超疏水性壁面的減阻效果有限，且在水壓及水流沖刷作用下空氣層很快會消失。持續(xù)加熱物體可以在壁面處形成連續(xù)穩(wěn)定的蒸汽膜，減小摩擦阻力。Saranadhi等[9]研究表明，當雷諾數(shù)在26 100～52 000范圍內(nèi)時，摩擦阻力可以減小80%～90%，減阻效果非常明顯。

膜態(tài)沸騰生成的蒸汽膜，不僅可以減小摩擦阻力，還可以減小壓差阻力。Vakarelski等[10]采用一種沸點和汽化潛熱均低于水的氟化液體，對膜態(tài)沸騰階段球體在該液體中的自由下落過程進行了研究，結(jié)果表明蒸汽膜的存在可以降低球體受到的壓差阻力，使總阻力減小85%. 隨后，Vakarelski等[11-12]采用不同液體對膜態(tài)沸騰球體繞流過程進行了詳細的實驗研究，分析了液體溫度、球體溫度和雷諾數(shù)等對減阻特性的影響。

目前，對膜態(tài)沸騰球體繞流特性的研究還處于初步階段，相關(guān)減阻機理還需要進一步研究。由于膜態(tài)沸騰球體繞流過程涉及相變、湍動、流動分離和多相流動等復雜物理問題，至今還沒有學者對該過程進行數(shù)值仿真研究。Gruncell等[14]對超疏水性球體的繞流過程進行了數(shù)值仿真研究，通過固定一層厚度相等且不變形的空氣層等效疏水特性。對于膜態(tài)沸騰球體繞流過程，不同表面位置蒸汽膜的厚度不同，且存在蒸汽堆積現(xiàn)象，因此用Gruncell等[14]的方法進行等效模擬的準確性較低。

本文采用Mixture多相流模型、剪切壓力傳輸(SST)k-ω湍流模型，耦合蒸發(fā)- 冷凝相變模型，對膜態(tài)沸騰球體繞流過程進行數(shù)值仿真研究。研究膜態(tài)沸騰球體繞流特性，分析雷諾數(shù)對阻力特性的影響，并對減阻機理進行研究。

1 數(shù)值計算方法

1.1 控制方程

為了識別球體表面蒸汽膜，球體壁面附近需要較密的網(wǎng)格，使得整體網(wǎng)格量較大，因此本文將球體繞流問題簡化為二維軸對稱流動問題，總體網(wǎng)格量約為44萬。

采用Mixture多相流模型描述各相流體體積分數(shù)。水和水蒸汽的混合密度為

ρm=ρl(1-αv)+ρvαv，

(1)

式中：ρl、ρv分別為水的密度和水蒸汽的密度；αv為水蒸汽的體積分數(shù)。

質(zhì)量守恒方程為

(2)

式中：t為時間；ui(i= 1, 2)為平面笛卡爾坐標系下的速度分量。

動量守恒方程為

(3)

式中：uj(j=1, 2)為平面笛卡爾坐標系下的速度分量；p為壓力；g為重力加速度；F為表面張力；τij為黏性剪切應力，

(4)

δij為Kronecker符號，μm為混合黏度；uk(k=1，2)為平面笛卡爾坐標系下的速度分量。

能量守恒方程為

(5)

式中：cp為定壓比熱；T為溫度；k為導熱系數(shù)；Φ為相變潛熱產(chǎn)生的能量源項，

(6)

Lv為汽化潛熱，fv=ρvαv/ρm為汽相質(zhì)量分數(shù)。

水汽化成水蒸汽的相間傳遞過程采用蒸發(fā)- 冷凝模型進行模擬，其輸運方程表達式為

(7)

(8)

(9)

cflv與cfvl分別為蒸發(fā)和冷凝系數(shù)，Ts為水的飽和溫度。

球體繞流問題，需要較精確地預測流動分離位置，因此本文湍流模型選用可以較好地模擬分離問題的SSTk-ω模型。

1.2 計算域及網(wǎng)格劃分

計算域、邊界條件和局部網(wǎng)格劃分如圖1所示。本文計算域采用二維軸對稱模型，球體直徑為D，計算域徑向直徑為10D，高度為15D. 球體壁面采用無滑移壁面邊界條件，側(cè)邊壁采用滑移壁面邊界條件，下邊界為速度入口，上邊界為壓力出口。球體位置固定不變，重力方向向下，來流方向向上，入口速度u為定值，可以等效模擬球體在水中勻速自由下落的過程。本文網(wǎng)格采用C-Block劃分策略，為了識別蒸汽膜，球體表面附近網(wǎng)格進行了加密，球體壁面附近第1層網(wǎng)格高度設置為1 μm，滿足SSTk-ω湍流模型所要求的y+<1，壁面附近網(wǎng)格高度增長因子設置為1.05，C-Block外層網(wǎng)格高度增長因子設置為1.1.

圖1 計算域及局部網(wǎng)格劃分Fig.1 Computational domain and local mesh generation

1.3 數(shù)值方法驗證

目前僅有少數(shù)學者針對膜態(tài)沸騰球體水下運動過程進行了實驗研究，本文將數(shù)值仿真得到的膜態(tài)沸騰球體繞流運動阻力系數(shù)與Vakarelski等[11]的膜態(tài)沸騰球體水下自由下落實驗得到的阻力系數(shù)進行對比，以驗證本文數(shù)值方法的準確性。球體水中自由下落過程中，受到水的阻力與運動速度的平方呈正比，球體下落的初始階段速度較小，球體做向下加速運動，當球體速度達到一定值時，水的阻力與球體重力相等，球體將做勻速運動。

球體水中運動的阻力系數(shù)為

Cd=2Fd/(πR2ρlu2)，

(10)

式中：Fd為球體受到的阻力；R為球體半徑；u為球體水中運動速度。

當球體勻速運動時阻力與重力相等時，阻力系數(shù)可以表達為

Cd=[8g(ρ-ρl)R]/(3ρlu2)，

(11)

式中：ρ為球體密度。

根據(jù)Vakarelski等[11]得到的球體勻速運動速度，利用(11)式計算出球體的阻力系數(shù)。

如圖1所示，本文數(shù)值仿真球體位置固定，水相對于球體向上運動，水流速度為u，重力方向向下，可以等效模擬Vakarelski等[11]的球體在水中勻速下落過程。表1為本文數(shù)值仿真與Vakarelski等[11]實驗得到的膜態(tài)沸騰球體水中勻速運動阻力系數(shù)對比，表中雷諾數(shù)Re=ρluD/μl，μl為水的黏度。

表1 數(shù)值與實驗阻力系數(shù)對比

從表1中可以看到，本文數(shù)值仿真選用的球體溫度與Vakarelski等[11]實驗選用的球體溫度不同。Vakarelski等[11]采用超疏水性球體，相對于親水性球體，超疏水性球體可以在較低溫度下形成膜態(tài)沸騰[11,15]。由于同時考慮流動、傳熱、相變和表面潤濕性比較復雜，不易于模擬，本文數(shù)值仿真沒有考慮球體表面潤濕性，為了形成膜態(tài)沸騰，球體溫度設置為450 ℃的較高溫度。Vakarelski等[10-12]的研究表明，球體發(fā)生膜態(tài)沸騰后會繼續(xù)升高球體溫度，球體達到勻速時的速度u變化較小，可以忽略，由(11)式可知，膜態(tài)沸騰球體的阻力系數(shù)隨球體溫度變化較小，可以忽略。因此，可以采用表1所示的不同球體溫度與Vakarelski等[11]的實驗結(jié)果進行對比。在Vakarelski等[11]的研究中，球體水中下落時間小于1.5 s，且由于蒸汽具有較低的導熱系數(shù)，球體溫度降低較小。由于沸騰球體繞流物理過程復雜，本文數(shù)值仿真進行了簡化，僅對定常狀態(tài)下的膜態(tài)沸騰球體繞流進行研究，球體溫度恒為450 ℃.

從表1中可以看到，3種不同Re下的數(shù)值仿真Cd值與Vakarelski等[11]實驗的Cd值最大相差-12.5%. 由于本文數(shù)值計算簡化為二維軸對稱問題，忽略了尾流不對稱性對阻力的影響，會與實際流動過程略有偏差。但阻力系數(shù)范圍與實驗范圍一致，并隨Re的增加而大幅度減小。因此可以采用本文數(shù)值仿真方法對膜態(tài)沸騰球體在水中的運動過程進行研究，并分析相應的減阻機理。

2 計算結(jié)果及分析

2.1 膜態(tài)沸騰球體繞流特性

本文對比分析Re=1.28×105條件下膜態(tài)沸騰球體與普通球體的繞流特性。其中水溫均為95 ℃，膜態(tài)沸騰球體球溫為450 ℃，普通球體球溫為95 ℃. 兩種球體的球體直徑D為20 mm，水流速度u=2 m/s. 計算所得的普通球體阻力系數(shù)為0.48，計算所得的膜態(tài)沸騰球體阻力系數(shù)為0.22(見表1)。

當普通球體流動處于亞臨界Re范圍內(nèi)，即5×103≤Re≤3×105時，阻力系數(shù)約為常數(shù)，Cd≈0.44. 本文計算所得的普通球體阻力系數(shù)在0.44附近，相差約為9.1%. 膜態(tài)沸騰球體的阻力系數(shù)約為普通球體的50%，由此可知膜態(tài)沸騰可以大幅度減小球體受到的阻力。

為了便于分析，對球體表面位置進行規(guī)定，如圖2所示，水流方向向上，球面底端位置P對應角度θ=0°，球面中端位置M對應θ=90°，球面頂端位置N對應θ=180°.

圖2 球體表面位置示意Fig.2 Schematic diagram of surface position of sphere

圖3(a)所示為膜態(tài)沸騰球體繞流的水蒸汽云圖，圖3(b)所示為球體周圍溫度云圖。當球體表面溫度遠高于水的沸點時，汽化成的水蒸汽匯聚到一起，使球體被一層水蒸汽膜包裹，如圖3(a)中M點處的水蒸汽云圖放大部分所示。球體尾部由于旋渦與回流的影響，沿流動方向生成汽泡云。由于球體周圍溫度大部分較低，且溫度范圍為95～450 ℃，跨度較大，圖3(b)僅顯示95～130 ℃部分的溫度云圖，其中溫度大于130 ℃部分顯示的顏色與130 ℃溫度一致。由于水蒸汽的導熱系數(shù)遠小于水，蒸汽膜的存在使高溫球體與周圍水隔開，降低了球體的傳熱能力，溫度場中大于130 ℃分布在球體表面附近的蒸汽膜內(nèi)，球體尾部的汽泡云溫度在130 ℃以下。

圖3 膜態(tài)沸騰球體繞流水蒸汽與溫度云圖Fig.3 Water vapor phase and temperature distribution

圖4所示為普通球體與膜態(tài)沸騰球體的繞流速度v矢量圖。由圖4可見：球體側(cè)面位置附近流速較高，尾部旋渦區(qū)速度較低；由于蒸汽膜的潤滑作用，膜態(tài)沸騰球體繞流速度最大值要大于普通球體；膜態(tài)沸騰球體尾部低速區(qū)比普通球體窄。

圖4 普通球體與膜態(tài)沸騰球體繞流速度矢量圖Fig.4 Velocity fields of ordinary sphere and film boiling sphere

圖5所示為普通球體與膜態(tài)沸騰球體的繞流流線圖。由圖5可見，球體在水中受到的阻力分為摩擦阻力與壓差阻力，摩擦阻力由水的粘滯性引起，壓差阻力由球體上表面和下表面的壓力差產(chǎn)生。對于普通球體，當5×103≤Re≤3×105時，摩擦阻力占總阻力的比例較小，球體主要受壓差阻力的作用。針對球體繞流問題，壓差阻力主要由球體表面發(fā)生流動分離的位置決定，流動分離發(fā)生的越遲緩、越靠近球體尾部，球體繞流受到的壓差阻力就越小。對于普通球體，球體表面流動分離發(fā)生在M點稍前的位置，即分離角度略小于90°[16]，與圖5(a)所示的結(jié)果相吻合。

圖5 普通球體與膜態(tài)沸騰球體繞流流線圖Fig.5 Streamline distributions of flows around ordinary sphere and film boiling sphere

Vakarelski等[10]對膜態(tài)沸騰球體在氟化液體中的自由下落過程進行了研究，發(fā)現(xiàn)膜態(tài)沸騰球體達到勻速時的速度要大于普通球體。結(jié)果表明，膜態(tài)沸騰階段產(chǎn)生的蒸汽膜使流動分離位置向尾部移動，減小了球體受到的壓差阻力，從而減小了球體受到的總阻力。與Vakarelski等[10]的研究相吻合，圖5(b)所示膜態(tài)沸騰球體繞流，分離位置位于球體表面的尾部區(qū)域。相比于普通球體，膜態(tài)沸騰球體尾部旋渦區(qū)較窄，尾部形狀趨于流線型。流線形狀與圖3所示水蒸汽相與溫度相云圖尾部形狀一致。本文的數(shù)值方法可以較好地模擬膜態(tài)沸騰球體繞流問題。

2.2 不同Re膜態(tài)沸騰球體繞流特性

對于普通球體，阻力系數(shù)大小主要與Re有關(guān)。當Re很小(Re<1)時，黏性力起主要作用，阻力系數(shù)可以表達為Cd=24/Re. 隨著Re增大，Cd不存在解析解，Cd隨著Re增大而減小。當5×103≤Re≤3×105時，流動分離點幾乎不發(fā)生改變，Cd幾乎與Re無關(guān)。當Re達到臨界值，即Re≈3.5×105時，邊界層轉(zhuǎn)化為湍流邊界層，使流動分離位置向后移動，阻力系數(shù)減小為Cd≈0.1. 不同于普通球體，Vakarelski等[10-12]的實驗表明，在亞臨界流動范圍內(nèi)，即5×103≤Re≤3×105時，膜態(tài)沸騰球體阻力系數(shù)隨Re增大而減小。

本文采用數(shù)值仿真方法研究亞臨界流動范圍內(nèi)Re對膜態(tài)沸騰球體阻力系數(shù)Cd的影響。Re分別取值5.06×104、1.28×105和2.24×105. 如表1所示，得到的阻力系數(shù)分別為0.38、0.22和0.14，與Vakarelski等[11]的實驗結(jié)果相吻合。這一結(jié)果表明，隨著Re增大，膜態(tài)沸騰球體繞流運動的阻力系數(shù)逐漸減小。在2.1節(jié)已經(jīng)對比分析了Re為1.28×105時普通球體與膜態(tài)沸騰球體的繞流特性，該節(jié)主要分析Re分別取值5.06×104和2.24×105時的膜態(tài)沸騰球體繞流特性，其中球體直徑均為20 mm，水流速度分別為0.79 m/s和3.50 m/s.

圖6～圖8所示分別為兩種Re條件下膜態(tài)沸騰球體繞流的水蒸汽相圖、速度矢量圖和流線圖。由圖6～圖8可以看出，Re越大，球體尾部水蒸汽云越趨于流線型，尾部低速旋渦區(qū)越窄，流動分離位置越靠近球體尾部。

圖6 不同Re膜態(tài)沸騰球體繞流水蒸汽相圖Fig.6 Water vapor phase distributions of flows around film boiling spheres with different Re

圖7 不同Re膜態(tài)沸騰球體繞流速度矢量圖Fig.7 Velocity fields of flows around film boiling spheres with different Re

圖8 不同Re膜態(tài)沸騰球體繞流流線圖Fig.8 Streamline distributions of flows around film boiling spheres with different Re

為了清晰對比不同雷諾數(shù)膜態(tài)沸騰球體繞流的流動分離特性，圖9給出了不同雷諾數(shù)膜態(tài)沸騰球體及常溫球體繞流近壁面流線對比圖。由于在亞臨界流動范圍內(nèi)，常溫球體繞流的流動分離位置幾乎保持一致，圖9僅呈現(xiàn)Re=1.28×105的普通球體近壁面流線。分析圖9可知，普通球體繞流的流動分離位置為θ≈85°，略大于Achenbach[16]實驗所得的θ≈82.5°.Re取值分別為5.06×104、1.28×105和2.24×105時膜態(tài)沸騰球體繞流的流動分離位置分別為θ≈99°、θ≈108°和θ≈120°.膜態(tài)沸騰球體繞流與普通球體相比，流動分離位置偏向尾部，且隨著Re增大，分離角度逐漸增大。在亞臨界流動范圍內(nèi)，阻力大小主要由壓差阻力決定，即由流動分離位置決定。因此，膜態(tài)沸騰繞流受到的阻力要小于普通球體，且隨著Re增大，減阻效果越明顯。

圖9 近壁面流線對比Fig.9 Comparison of streamlines near the sphere surface

2.3 膜態(tài)沸騰球體繞流運動減阻機理

流動分離的條件為逆壓梯度，流動分離發(fā)生的根本原因為流體黏性的存在。在分離區(qū)附近，流體的動能不僅逐漸轉(zhuǎn)變?yōu)閴簭娔?，同時邊界層內(nèi)的粘性流體(特別在臨近壁面處)受到壁面的滯止作用(表現(xiàn)為切應力)，不斷消耗動能。當靠近壁面處流體的動能全部耗盡，流體微團將在分離點處停滯不前，并發(fā)生堆積。分離點后在反向壓差作用下，壁面附近產(chǎn)生回流，排擠上游流體使流動在分離點處發(fā)生分離。Vakarelski等[10-11]的實驗研究指出，膜態(tài)沸騰球體被一層蒸汽膜包裹，減小了切應力對流體的滯止作用，使流動分離點向后移動。

圖10所示為Re=1.28×105時普通球體與膜態(tài)沸騰球體繞流在圖2所示M點附近的速度矢量對比圖。由圖10可見，普通球體壁面附近速度梯度較小，距離壁面較遠位置速度才達到外部主流速度。由于水的黏度遠大于水蒸汽，膜態(tài)沸騰球體在蒸汽膜內(nèi)，速度梯度較大，在距離壁面較近的蒸汽膜外部具有較大的流動速度。且膜態(tài)沸騰球體外層速度與來流速度平行，普通球體速度稍偏離來流方向偏向遠離球面的方向。圖11給出了M點處沿y軸方向流體的x軸方向無量綱速度v(v=vx/u，vx為x軸方向的流動速度)對比，其中來流方向為x軸負方向。球體直徑為20 mm，M點處y軸坐標值為10 mm(位置與方向規(guī)定見圖2)。由圖11可知：普通球體M點位于流動分離點之后近壁面位置發(fā)生回流，因此普通球體近壁面處x軸方向速度為正值(與來流方向相反)；蒸汽膜的存在，使膜態(tài)沸騰球體沿y軸方向流體速度增長較快，且達到的速度最大值略大于普通球體。

圖10 普通球體與膜態(tài)沸騰球體在M點附近速度矢量Fig.10 Velocity vectors of ordinary sphere and film boiling sphere near point M

圖11 M點處沿y軸方向流體的x軸方向無量綱速度對比Fig.11 Comparison of fluid velocities in the direction of x at point M

將蒸汽膜與球體看作整體，膜態(tài)沸騰球體的繞流運動可以等效為流體繞蒸汽膜的運動，則壁面處的無滑移邊界條件可以轉(zhuǎn)換為蒸汽膜處的滑移邊界條件。蒸汽膜邊界處黏滯力的減小使膜態(tài)沸騰球體繞流流動分離點向球體尾部移動，減小了球體受到的阻力。

由2.2節(jié)可知，Re對膜態(tài)沸騰球體繞流的阻力特性有較大影響。Gruncell等[14]通過在球體表面固定一層厚度相等且不變形的氣體層，對超疏水球體的繞流過程進行了數(shù)值仿真研究，結(jié)果表明氣體層的滑移效應使流動分離點向球體尾部移動，極大地減小了球體繞流受到的阻力。對于膜態(tài)沸騰球體的繞流過程，球體表面的蒸汽膜厚度并不處處相等。圖12所示為不同Re膜態(tài)沸騰球體繞流蒸汽相體積分數(shù)為0.5的等值線圖。由圖12可知：蒸汽會在逆壓梯度作用下發(fā)生堆積，使蒸汽膜外部流動發(fā)生分離；類似于普通球體，隨著Re增大，堆積位置會向球體尾部移動；蒸汽膜的減阻作用，會由于蒸汽的堆積而受到抑制。因此膜態(tài)沸騰球體繞流過程中球體的阻力系數(shù)隨著Re增大而減小。

圖12 不同Re水蒸汽相體積分數(shù)為0.5的等值線圖Fig.12 Contour plot of 0.5 vapor fraction for different Re

3 結(jié)論

本文采用數(shù)值仿真方法對亞臨界流動范圍內(nèi)膜態(tài)沸騰球體的繞流特性進行了研究，分析了Re對阻力特性的影響，并對膜態(tài)沸騰減阻機理進行了研究。主要得到以下結(jié)論：

1)本文數(shù)值仿真方法可以較好地模擬膜態(tài)沸騰球體繞流特性，得到的阻力系數(shù)與已有文獻[11]的一致性較好。

2)由于水蒸汽的黏度遠小于水，壁面附近蒸汽膜內(nèi)部的速度梯度較大，相比于普通球體，膜態(tài)沸騰球體繞流在距離壁面較近的水蒸汽膜外部具有較大的流動速度，水蒸汽膜的滑移效應減小了壁面對流體的粘滯力作用，使流動分離點向尾部移動，使球體受到的阻力減小。

3)由于逆壓梯度的存在，膜態(tài)沸騰球體生成的水蒸汽會在球體壁面尾部區(qū)域堆積，隨著Re增大，堆積位置向后移動，延緩流動分離，阻力系數(shù)逐漸減小。