水平管道過冷沸騰換熱的非穩(wěn)態(tài)數(shù)值計算

劉 征,余志毅,2,孫偉華,張 珂

(1.北京理工大學(xué) 機械與車輛學(xué)院,北京 100081;2.水沙科學(xué)與水利水電工程國家重點實驗室(清華大學(xué)),北京 100084)

過冷沸騰是指壁面溫度高于飽和溫度,液體主體溫度低于飽和溫度下所發(fā)生的局部沸騰現(xiàn)象,作為一種強化換熱方式應(yīng)用于內(nèi)燃機、電子元器件、充電電纜、車用電池以及大功率激光武器等重要設(shè)備的散熱[1-7]。有學(xué)者針對過冷沸騰流動換熱進行了大量數(shù)值計算方面的研究,可以分為兩方面:一類是穩(wěn)態(tài)計算,這類計算的計算量較小,一般針對相對復(fù)雜的結(jié)構(gòu),只能獲取各參數(shù)的結(jié)果分布[5-10];另一類是同時反映參數(shù)分布規(guī)律和瞬時特性非穩(wěn)態(tài)計算,側(cè)重于對流動過程參數(shù)的獲取,計算量較大,更能接近實際沸騰過程。魏敬華等[11]對豎直窄縫通道內(nèi)的過冷流動沸騰進行了數(shù)值計算研究,分析了不同壓力和熱流密度下的汽泡動力學(xué)特性;潘良明等[12]建立了過冷流動沸騰凝結(jié)模型,以分析汽泡凝結(jié)過程的流場特性;葛蘇槿[13]從含氣率、物理場、汽泡行為等方面對一傾斜管道內(nèi)過冷沸騰換熱特性進行了數(shù)值研究;Lee等[14-17]針對垂直管道內(nèi)FC-72的過冷沸騰進行了大量研究,首先研究了上升流中的高過冷度核態(tài)沸騰,總結(jié)了換熱系數(shù)沿軸向的變化規(guī)律,提出了一種適用于VOF模型的剪切力計算方法,并研究了高過冷度沸騰在CHF點前后的流動傳熱特性以及微重力條件下的低過冷度沸騰中的汽泡行為和流型演化;Yeo等[18]在不涉及蒸干的情況下成功模擬了微通道中的過冷沸騰,通過與實驗以及經(jīng)驗關(guān)聯(lián)式驗證了數(shù)值計算方法的有效性;Chen等[19]對矩形微通道的過冷沸騰進行了數(shù)值模擬,分析了含氣率和熱流密度對換熱的影響;曹濤濤[20]基于水平通道中的過冷沸騰,提出了確定Lee模型中傳質(zhì)系數(shù)的經(jīng)驗關(guān)聯(lián)式。

現(xiàn)有研究中針對過冷沸騰換熱過程換熱參數(shù)波動和分布規(guī)律方面的研究較少,實際上換熱參數(shù)的波動和分布是氣液相變換熱設(shè)備設(shè)計的重要依據(jù),參數(shù)波動過大產(chǎn)生的熱應(yīng)力會降低換熱設(shè)備的可靠性,而換熱參數(shù)的分布規(guī)律是確定換熱設(shè)備幾何尺寸的重要參考。基于此,針對一宏觀尺度中小直徑水平管道內(nèi)的過冷沸騰換熱過程進行非穩(wěn)態(tài)數(shù)值計算,分析沸騰過程中流動和換熱參數(shù)的分布及波動特性,并對過冷沸騰換熱的規(guī)律進行總結(jié),為新型換熱設(shè)備的相關(guān)設(shè)計提供借鑒經(jīng)驗。

1 物理模型及網(wǎng)格劃分

參考目前應(yīng)用廣泛的換熱設(shè)備(如閉式冷卻系統(tǒng)中的冷卻水套以及板翅式換熱器)的通道高度,利用ICEM_CFD 2020 R2建立了一個高度為10 mm的水平管道二維模型,如圖1所示。管道總長為800 mm,其中入口段200 mm,加熱段100 mm,出口段500 mm。采用下壁面加熱方式,底部加熱面與流固耦合面之間的固體域厚度為0.5 mm。

圖1 物理模型及計算域網(wǎng)格分布

通過ICEM_CFD軟件對計算域劃分4套網(wǎng)格,網(wǎng)格總數(shù)分別為65 950、107 250、225 000、274 000。選擇基準工況,以加熱壁面固體域的平均溫度作為檢驗標準。計算發(fā)現(xiàn),當(dāng)網(wǎng)格總數(shù)在225 000以上時,加熱壁面固體域的平均溫度變化幅度已經(jīng)很小。綜合考慮計算量與計算精度,最終選擇第3套網(wǎng)格,即總數(shù)為225 000的網(wǎng)格,其中,流體域縱向網(wǎng)格數(shù)為65,橫向網(wǎng)格數(shù)為3 600,對加熱段、出口延伸段前部以及壁面附近的網(wǎng)格進行加密,第1層網(wǎng)格高度為0.05 mm,YPlus值為3.78。

2 計算模型及設(shè)置

2.1 計算模型的選取

VOF多相流模型是一種針對兩種或兩種以上流體的界面跟蹤技術(shù),其優(yōu)勢是能夠清晰地捕捉到兩相之間的交界面。為了更準確地捕捉汽泡的生成發(fā)展過程,分析汽泡的運動規(guī)律,選擇VOF模型對過冷沸騰換熱過程進行計算。

湍流模型選擇SSTk-ω模型,湍流阻力因子設(shè)置為10,以對壁面處的汽泡生成發(fā)展的非穩(wěn)態(tài)過程進行更好的描述[21]。

此外,相變模型是過冷沸騰數(shù)值計算所涉及的重要模型。由Lee在1980年提出的半隱式傳熱傳質(zhì)模型——Lee模型[22],是目前應(yīng)用最為廣泛的相變模型。該模型將相變過程簡化為飽和溫度與實際溫度的溫差驅(qū)動的傳質(zhì)過程,其形式簡單且能夠模擬沸騰和凝結(jié)過程的全階段,Lee模型的基本方程如下[21-22]:

連續(xù)方程為

(1)

Lee模型相變過程質(zhì)量傳遞源項求解方式為

(2)

能量傳遞源項為質(zhì)量源項與汽化潛熱的乘積,表示為

Q=mhfg

(3)

式中:α為相體積分數(shù);ρ為密度,kg/m3;v為速度矢量;Q為能量傳遞源項,J/(s·m3);m為質(zhì)量傳遞源項,kg/(s·m3);hfg為汽化潛熱,J/kg;T為流體溫度,K;Tsat為飽和溫度,K;r為相變傳質(zhì)系數(shù),s-1,其取值需要根據(jù)具體工況進行調(diào)整,本研究所采取的方式為結(jié)合計算關(guān)聯(lián)式的結(jié)果對傳質(zhì)系數(shù)的取值進行若干次試算,最終在3種不同熱流密度(150、200、250 kW/m2)下,汽化傳質(zhì)系數(shù)分別確定為0.1、0.5、2.0。鑒于氣相密度較低,一般認為凝結(jié)傳質(zhì)系數(shù)應(yīng)當(dāng)取較大值,Kim等[23]認為凝結(jié)傳質(zhì)系數(shù)應(yīng)當(dāng)在汽化傳質(zhì)系數(shù)的基礎(chǔ)上乘以液體密度與氣體密度的比值,以此為基礎(chǔ),近似將汽化傳質(zhì)系數(shù)與凝結(jié)傳質(zhì)系數(shù)的比確定為1/1 600,故凝結(jié)傳質(zhì)系數(shù)的取值分別為160、800、3 200;下標v和l分別代表氣相和液相,下標e和c分別代表汽化和凝結(jié)過程。

此外,基于核沸騰起始需要一定的過熱條件,通過UDF功能將Bergles提出的沸騰起始點計算關(guān)聯(lián)式[24]嵌入由Lee模型構(gòu)建的傳熱傳質(zhì)源項表達式中,計算沸騰起始階段所需過熱度,作為計算沸騰核化前單相液體換熱的補充,計算關(guān)聯(lián)式為

(4)

式中:p為絕對壓力,Pa;qONB為沸騰起始點熱流密度,W/m2;TONB為該熱流密度下所對應(yīng)的沸騰起始點溫度,K;TONB-Tsat即所在位置達到沸騰起始點所需過熱度,沸騰起始點溫度通過所在位置相鄰接觸面的局部熱流密度以及所在位置的局部壓力求得。

基于此,將成核位置限制在靠近加熱面的近壁區(qū),當(dāng)流體存在一定過熱,但局部區(qū)域不存在相界面時,認為該區(qū)域處于沸騰核化前的單相階段,由于此時將要發(fā)生的汽化過程位于液相內(nèi)部,近似以沸騰起始點溫度作為此時的汽化溫度,將式(2)的質(zhì)量傳遞源項修正為

(5)

當(dāng)沸騰過程存在相界面或發(fā)生氣相的凝結(jié)過程時,質(zhì)量源項求解仍使用式(2)。最后,通過Fluent中的UDF功能,將當(dāng)前質(zhì)量源項的求解結(jié)果引入計算。綜合考慮式(4)和(5),質(zhì)量源項可以表達為

(6)

2.2 邊界條件和其他設(shè)置

邊界條件方面,入口邊界流速為0.4 m/s,過冷度為5 K;出口壓力設(shè)置為常壓101 325 Pa;加熱壁面設(shè)置3種熱流密度,分別為150、200、250 kW/m2;表面張力的求解采用Brackbill等提出的CSF表面張力模型[25]。

壓力速度的耦合采用SIMPLEC算法。對于標量方程各項的空間離散,壓力采用PRESTO!格式,體積分數(shù)采用幾何重構(gòu)格式,能量和動量采用二階迎風(fēng)格式,湍動能和湍流耗散率采用一階迎風(fēng)格式;為保證計算的收斂性,采用可變時間步長,范圍為5×10-6～5×10-4s,限制全局庫朗數(shù)為1以下,每個時間步的最大迭代次數(shù)設(shè)置為50,計算總時長設(shè)置為1.2 s。

3 計算結(jié)果分析

3.1 結(jié)果合理性驗證

數(shù)值計算結(jié)果的可靠性通過過冷流動沸騰換熱經(jīng)驗關(guān)聯(lián)式驗證。常用的沸騰換熱關(guān)聯(lián)式中適用于過冷沸騰的關(guān)聯(lián)式包括Chen關(guān)聯(lián)式、Moles-Shaw關(guān)聯(lián)式、Shah關(guān)聯(lián)式的修正形式、Gungor-Winterton關(guān)聯(lián)式以及Liu-Winterton關(guān)聯(lián)式等。Chen[26]關(guān)聯(lián)式作為最早的疊加關(guān)聯(lián)式,將流動沸騰換熱看作是液相強制對流換熱和核態(tài)沸騰換熱兩種機制共同作用,通過兩種機制的線性疊加求解流動沸騰過程的對流換熱系數(shù)。此后在Chen關(guān)聯(lián)式的基礎(chǔ)上衍生出若干更為精確的沸騰傳熱關(guān)聯(lián)式,如Gungor-Winterton關(guān)聯(lián)式、Liu-Winterton關(guān)聯(lián)式等。Gungor-Winterton[27]關(guān)聯(lián)式在Chen關(guān)聯(lián)式基礎(chǔ)上針對核態(tài)沸騰換熱過程的對流換熱系數(shù)求解采用了更精確的Cooper[28]關(guān)聯(lián)式,針對過冷沸騰考慮了核態(tài)沸騰換熱和強制對流換熱驅(qū)動溫差不同做出改進,提升了求解精度;Liu-Winterton[29]關(guān)聯(lián)式與Gungor-Winterton關(guān)聯(lián)式類似,區(qū)別為Liu-Winterton關(guān)聯(lián)式采用了非線性疊加方式。適用范圍方面,Shah[30]過冷沸騰關(guān)聯(lián)式的驗證數(shù)據(jù)多數(shù)集中在普朗特數(shù)1以下,雷諾數(shù)10 000以上; 而Gungor-Winterton關(guān)聯(lián)式與Liu-Winterton關(guān)聯(lián)式采用了相同的數(shù)據(jù)庫,適用于雷諾數(shù)為568.9～875 000、普朗特數(shù)0.83～9.1的工況范圍[29],可覆蓋本文數(shù)值計算的工況范圍。基于以上,選擇了Gungor-Winterton關(guān)聯(lián)式對數(shù)值計算進行驗證,其計算過冷沸騰對流換熱系數(shù)的表達式如下:

q=hl(Tw-Tb)+Shpool(Tw-Tsat)

(7)

(8)

(9)

(10)

F=1+24 000Bo1.16

(11)

(12)

式中:q為熱流密度,W/m2;hl為單相對流換熱系數(shù),W/(m2·K);hpool為沸騰換熱系數(shù),W/(m2·K);F為對流換熱強化因子,過冷沸騰不考慮該系數(shù),取值為1;S為沸騰抑制因子;kl為液體熱導(dǎo)率,W/(m·K);d為水力直徑,m;Rel為液相雷諾數(shù);Prl為液相普朗特數(shù);Bo為沸騰數(shù);G為質(zhì)量流速,kg/(m2·s);pr為對比壓力;M為相對分子質(zhì)量。

3種熱流密度條件下,數(shù)值計算與實驗關(guān)聯(lián)式計算結(jié)果對比如表1所示。可以看出,盡管存在一定偏差,但數(shù)值計算與關(guān)聯(lián)式計算結(jié)果具有相同的趨勢,認為當(dāng)前所計算結(jié)果是合理的。

表1 數(shù)值計算與Gungor-Winterton關(guān)聯(lián)式計算結(jié)果對比

3.2 過冷沸騰流場的演變與換熱特性波動

以熱流密度q為250 kW/m2的工況為例,流場在計算時間t為0.5 s時刻后到達相對穩(wěn)定狀態(tài)。圖2為該工況下0.7～1.2 s過冷沸騰流場中氣相分布的變化及局部區(qū)域矢量圖。可以看出,由于下壁面的加熱,汽泡不斷形成并隨主流向下游運動。總體上看,加熱段的前半部分汽泡較細、較少,處于核化、脫離和初步生長階段,而后半部分汽泡明顯增大,氣泡的合并和消亡現(xiàn)象更明顯。伴隨著汽泡的核化、生長、脫離、消亡等過程,壁面與流體間的換熱受到干擾。如圖2所示,由0.9和1.2 s時刻壁面附近的局部流場圖可見,汽泡的生長和運動等行為對壁面附近流動造成強烈的擾動,并形成明顯的漩渦。這些擾動使得管道換熱面的換熱性能得到強化的同時,將造成換熱特性的波動。

圖2 熱流密度250 kW/m2工況加熱段氣相分布隨時間t的變化及部分局部放大區(qū)矢量圖

圖3為與圖2相對應(yīng)的溫度分布。可以看出,從加熱段進口到加熱段約中間位置,熱邊界層沿流動方向逐漸增厚,此后該厚度無明顯增加,在一均值上下波動。基于此特點,可將管內(nèi)流動沸騰分為熱邊界層發(fā)展段與核沸騰主導(dǎo)段,250 kW/m2工況下兩者的分界線約在橫坐標X等于250 mm處。

在熱邊界層發(fā)展段,加熱表面放出的熱量部分用于加熱壁面附近流體,使得該區(qū)域流體的溫度逐漸升高,另一部分則通過液體的汽化過程被吸收;當(dāng)熱邊界層溫度達到飽和溫度,邊界層內(nèi)熱傳遞逐漸趨向穩(wěn)定,進入核沸騰主導(dǎo)段,此時,主流區(qū)的溫升以及汽化過程使得壁面附近流體溫度隨位置變化相對較小,壁面所放出的熱量主要通過汽化過程被吸收。同時,對比圖2和3可知,核沸騰主導(dǎo)段的流場隨時間的波動更為劇烈,表明沸騰過程增加了流場中的不穩(wěn)定性。

圖3 熱流密度250 kW/m2工況加熱段溫度分布隨時間t的變化

圖4分別為加熱區(qū)域各部分對流換熱系數(shù)隨時間的波動,每個樣本點為0.05 s時間段內(nèi)該區(qū)域的平均值。對比發(fā)現(xiàn),各段對流換熱系數(shù)隨時間均存在波動性,熱邊界層發(fā)展段前部的換熱系數(shù)波動幅值較小,而其下游區(qū)域換熱系數(shù)波動幅值沿流動方向整體呈現(xiàn)增加趨勢。造成前者波動趨勢的原因是熱邊界層發(fā)展段前部的壁面溫度與近壁面液相溫度較低,顯熱帶走了較多的熱量,受到汽泡核化以及生長脫離過程影響較小,換熱系數(shù)波動不明顯,另一方面該段壁溫相對較低,汽泡脫離直徑較小,使得汽泡生長產(chǎn)生熱阻所造成的影響更小,該段各部分的對流換熱系數(shù)平均波動幅值小于核沸騰主導(dǎo)段;而下游區(qū)域近壁區(qū)溫度升高造成壁面換熱特性受汽泡生長脫離周期的影響明顯,對流換熱系數(shù)產(chǎn)生明顯波動,受到汽泡脫離直徑增加的影響,波動幅值整體較大。如200～220 mm段在0.5 s記錄時間段內(nèi)的波動幅值為414.16 W/(m2·K),而>260～280 mm段以及>280～300 mm段的波動幅值分別為1 196.33和815.93 W/(m2·K),>220～300 mm整體平均波動幅值為794.14 W/(m2·K),約為熱邊界層發(fā)展段前部波動幅值的2倍,且波動規(guī)律相似,表明汽泡的生長運動加劇了對流換熱系數(shù)的波動,使得入口區(qū)域以外的波動規(guī)律具有相似特性。

3.3 熱流密度對過冷沸騰流場的影響

圖5～7分別為t為1.0 s時刻3個工況下的溫度分布、氣相分布以及速度分布,三者之間相互影響。如圖5所示,加熱段熱流密度的變化直接影響加熱段近壁區(qū)溫度分布,不同熱流密度下,壁面附近液體溫度升高速度不同,沸騰發(fā)生的位置和強度相應(yīng)不同,對管內(nèi)速度場產(chǎn)生不同影響。如圖6所示,隨著熱流密度增加,壁面附近出現(xiàn)明顯汽泡的位置提前。q為150 kW/m2工況下壁面附近存在的汽泡數(shù)量很少,汽泡的生長、合并與脫離現(xiàn)象很不明顯;與之相對,q為250 kW/m2的工況下汽泡活動范圍得以擴大,其在加熱段前部即產(chǎn)生細泡狀流型,而后由于壁面溫度的升高以及汽泡的生長合并使得其直徑不斷增大。同時,各工況下產(chǎn)生的汽泡對近壁區(qū)域流場產(chǎn)生不同程度的擾動,汽化過程造成的膨脹使得下游流速不同程度的升高,如圖7所示,熱流密度越大,產(chǎn)生的汽泡越多,近壁區(qū)域流場速度分布越不規(guī)則,下游速度越大,其中,流體速度用v表示。

圖5 各工況加熱段溫度分布

圖6 各工況加熱段氣相分布

圖7 各工況加熱段速度分布

3.4 過冷沸騰對流換熱系數(shù)隨位置的變化

由于當(dāng)前工況水平管道的過冷沸騰是從單相逐漸轉(zhuǎn)變?yōu)閮上?加熱面的對流換熱系數(shù)隨流動位置變化,不同工況下對流換熱系數(shù)隨位置的變化存在一定差異。

圖8為3個工況下,局部平均對流換熱系數(shù)隨位置的變化。與換熱的波動特性類似,造成各局部平均換熱系數(shù)變化的原因同樣包含熱邊界層沿流動方向的發(fā)展以及沸騰對傳熱的強化兩種機制。

圖8 0.7～1.2 s對流換熱系數(shù)h沿水平坐標X的變化

對比圖5可以看出,熱邊界層發(fā)展段對流換熱系數(shù)沿流動方向不斷下降,熱流密度越大,該區(qū)域?qū)α鲹Q熱系數(shù)沿流動方向變化越快,發(fā)展所需通道的長度越短,如圖8所示,150 kW/m2工況下該段長度為80 mm,200 kW/m2工況下為60 mm,250 kW/m2工況下為50 mm。這是由于熱邊界層處于發(fā)展階段時,發(fā)展初期的熱邊界層很薄,近壁區(qū)液體溫度較低且與主流區(qū)溫度相近,由于溫度場的連續(xù)性,壁面溫度較低且與主流區(qū)的溫差較小,單相導(dǎo)熱及對流換熱強度較高,隨著壁面對附近流體的加熱,近壁區(qū)域熱邊界層增厚,熱阻相應(yīng)變大,換熱系數(shù)下降。熱流密度越大,近壁區(qū)域的局部過熱層產(chǎn)生越快,因此,熱邊界層發(fā)展所需的通道長度越短。

隨著熱邊界層發(fā)展到一定程度,近壁區(qū)域的過熱層內(nèi)部發(fā)生沸騰,流場中開始出現(xiàn)明顯的汽泡,沸騰對傳熱的強化作用不斷增強,對流換熱系數(shù)由逐漸下降轉(zhuǎn)變?yōu)橼呌诜€(wěn)定,進入核沸騰主導(dǎo)段。由于汽泡在各段的生長階段不同,進入核沸騰主導(dǎo)段后對流換熱系數(shù)隨位置的變化存在一定小范圍浮動,表明該階段熱邊界層發(fā)展對傳熱的抑制減緩,核沸騰對換熱的強化作用增強。

4 結(jié) 論

1)短通道內(nèi)過冷沸騰流場的演變及換熱特性波動受熱邊界層發(fā)展和沸騰兩種機制的影響,據(jù)此可將短通道加熱區(qū)域分為熱邊界層發(fā)展段和核沸騰主導(dǎo)段。熱邊界層發(fā)展段的前部換熱系數(shù)波動幅度較小,原因是入口附近的近壁區(qū)溫度較低,受沸騰過程影響較小;熱邊界層發(fā)展段后部及核沸騰主導(dǎo)段在壁面附近形成過熱液體層,受沸騰影響較大,換熱系數(shù)的波動幅度增大。

2)熱流密度可促進水平管內(nèi)流動參數(shù)和換熱參數(shù)沿流動方向的變化。一方面,熱流密度增加,近壁區(qū)溫度、平均含氣率及流動速度沿流動方向的提升越快,沸騰越劇烈;另一方面,熱流密度的增加使得對流換熱系數(shù)隨熱邊界層發(fā)展達到相對穩(wěn)定所需的距離越短。

綜上,針對高熱流密度工況下的過冷流動沸騰換熱,可以采取適當(dāng)減小通道長度的方式,以獲得相對更大的對流換熱系數(shù),同時減小沸騰引起的流場參數(shù)波動和不穩(wěn)定性。