李宗坤,范廣銘,*,曾曉波,嚴一鳴,馬富兵

(1.哈爾濱工程大學 核科學與技術學院,黑龍江 哈爾濱 150001;2.哈爾濱工程大學 黑龍江省核動力裝置性能與設備重點實驗室,黑龍江 哈爾濱 150001)

由于相變的存在,過冷沸騰作為一種高效的傳熱方式被廣泛應用,因此許多學者對過冷沸騰的傳熱計算開展了相關研究。隨著傳熱學的發(fā)展,過冷沸騰傳熱計算已逐漸從純經(jīng)驗模型過渡到半經(jīng)驗模型,再過渡到如今的機理模型,其中汽泡動力學和壁面熱流分配模型在過冷沸騰熱流密度計算中的應用受到廣泛關注與研究。

在采取壁面熱流分配機理模型計算壁面熱流密度時,汽泡脫離直徑、浮升直徑等參數(shù)對壁面熱流密度計算具有很大的影響。Klausner等[1]和Zeng等[2]的汽泡力平衡模型是目前應用較多的求解汽泡相關參數(shù)的模型。在最近的研究中,Akand等[3-4]分別對豎直加熱壁面和以IVR(implementation of in-vessel retentio)為背景的傾斜加熱壁面的汽泡受力分析開展了相關研究,并對汽泡浮升直徑和臨界熱流密度(CHF)進行了計算,計算得到的汽泡浮升直徑的平均誤差約為9.89%,CHF的平均誤差約為7.2%?；谄輨恿W研究,王暢等[5]建立了窄矩形通道內(nèi)的過冷流動沸騰傳熱系數(shù)計算關系式,計算誤差可控制在±30%以內(nèi)。Ren等[6-8]分別對汽化核心密度、窄矩形通道汽泡滑移、結(jié)合汽泡動力學的汽泡脫離直徑和汽泡浮升直徑等開展了研究。

對于壁面熱流分配機理模型,經(jīng)典RPI模型是廣為接受的機理模型,該模型最初應用于池式沸騰的熱流密度計算,后來進一步發(fā)展到過冷流動沸騰中。在RPI模型的基礎上,Sateesh等[9]考慮了豎直加熱壁面的汽泡滑移效應,Yeoh等[10]基于垂直環(huán)形通道中向上強迫對流過冷沸騰開發(fā)了壁面熱流分配機理模型,閆美月[11]針對矩形通道內(nèi)強迫對流過冷沸騰建立了壁面熱流分配機理模型,以達到計算矩形通道內(nèi)壁面熱流密度的目的。

目前,壁面熱流分配模型多數(shù)以汽泡力平衡模型為基礎,并且在目前所公開的文獻中,多數(shù)研究集中在垂直或水平平板加熱工況。然而,實際情況中存在一些特殊的棒束通道,如傾斜安裝的換熱器、海洋核動力浮動平臺或海洋條件下船舶核動力裝置中的棒束通道,棒束長時間處于傾斜狀態(tài)[12-13]。在傾斜條件下,由于汽泡受力方向的改變,傾斜棒束壁面汽泡可能會直接脫離加熱壁面而不存在滑移現(xiàn)象,也可能出現(xiàn)汽泡繞加熱棒曲面滑移的現(xiàn)象。然而,壁面熱流分配機理模型很大程度上依賴于汽泡相關參數(shù),在計算壁面熱流密度的過程中,傾斜條件對汽泡行為以及汽泡參數(shù)會產(chǎn)生一定影響,進而對壁面熱流分配模型的計算產(chǎn)生很大影響。因此,針對傾斜條件下棒束壁面的傳熱計算,有必要結(jié)合汽泡動力學和壁面熱流分配模型開展壁面熱流密度計算的進一步研究。本文主要進行傾斜條件下汽泡力平衡模型的建立與計算,在力平衡模型基礎上建立壁面熱流分配模型,并與實驗數(shù)據(jù)進行對比驗證。

1 汽泡力平衡模型

根據(jù)Klausner經(jīng)典汽泡受力理論[1],對傾斜條件下棒束的上加熱壁面進行汽泡受力分析,如圖1所示,其中汽泡在x方向和y方向的受力方程可表示為:

圖1 傾斜加熱棒表面汽泡受力分析

∑Fx=Fs,x+Fsl+Fh+

Fcp+Fbsinθinc+Fdu,x

∑Fy=Fs,y+Fqs+Fbcosθinc+Fdu,y

(1)

式中:Fs為表面張力,N;Fsl為剪切升力,N;Fh為水動力壓力,N;Fcp為接觸壓力,N;Fb為浮力,N;θinc為傾斜角度,(°);Fdu為非對稱生長力,N;Fqs為流動曳力(準穩(wěn)態(tài)曳力),N;下標x代表垂直于流向的方向,y代表平行于流動方向的方向。

金頔等[14]和肖仁杰等[15]也曾把傾斜角度的影響考慮在汽泡受力分析上,但對于汽泡脫離或浮升的條件,均是x方向或y方向的力打破平衡,即當∑Fx>0或∑Fy>0時,汽泡開始脫離或浮升。但傾斜條件的存在會導致近壁面汽泡行為發(fā)生改變。一方面,實驗現(xiàn)象拍攝結(jié)果表明,在傾斜上加熱壁面的汽泡,滑移距離非常短,并且多數(shù)汽泡會沿浮力方向(豎直向上)直接脫離加熱壁面,如圖2所示;另一方面,相比于汽泡尺寸,相鄰加熱棒壁面對近壁面汽泡行為影響較小,即在浮力引起汽泡脫離的方向上無限制汽泡脫離的其余壁面,且根據(jù)Zeng等[2]的研究,浮力在汽泡脫離過程中占據(jù)重要地位。因此,對于傾斜上加熱壁面汽泡脫離判定條件定為浮力方向(豎直向上)受力總和大于0,即∑Fg>0時,汽泡脫離加熱壁面,如式(2)所示:

圖2 汽泡脫離示意圖

∑Fg=Fb+Fqscosθinc+(Fsl+Fh+Fcp)·

Fs,xsinθinc+Fs,ycosθinc>0

(2)

式中,θi為汽泡傾斜角。

根據(jù)Zeng等[2]的研究,相比于其余汽泡受力的量級,水動力壓力和接觸壓力的量級遠小于其余汽泡受力量級。因此,在處理傾斜加熱壁面汽泡受力時,水動力壓力和接觸壓力忽略不計,即Fh、Fcp→0。簡化后的方程如式(3)所示:

∑Fg=Fb+Fqscosθinc+

Fs,xsinθinc+Fs,ycosθinc>0

(3)

1) 浮力

浮力是汽泡在重力場中由于汽液密度差引起的力,計算公式如式(4)所示:

(4)

式中:Vb為汽泡體積,m3;rb為汽泡生長半徑,m;ρl為液相密度,kg/m3;ρv為汽相密度,kg/m3;g為重力加速度,m/s2。

2) 非對稱生長力

非對稱生長力是汽泡由于非對稱生長而導致的力,本文采用Klausner方法計算,如式(5)所示:

(5)

式中,r′b和r″b分別為汽泡生長速度和汽泡生長加速度。沿汽泡傾斜角(θi)方向分解為x、y方向的對稱生長力,其中汽泡傾斜角根據(jù)Ren等[8]的研究取15°。

汽泡生長速度和汽泡生長加速度的計算需要汽泡尺寸和時間的相對關系,即汽泡生長控制方程。目前,Zuber[16]汽泡生長模型應用較廣泛,因此,本文采用Zuber模型,如式(6)所示:

(6)

式中:b為考慮汽泡非球狀的經(jīng)驗修正系數(shù),Zeng的研究表明b的范圍為1～1.73,而且b取1時,其實驗汽泡浮升直徑和脫離直徑符合最好;Ja為雅可比數(shù);α為熱擴散率;t為汽泡生長時間,s;kl為液相熱導率,W/(m·℃);cpl為比定壓熱容,kJ/(kg·℃);Tw為加熱壁面溫度,℃;Tsat為飽和水溫度,℃;γ為汽化潛熱,kJ/kg。根據(jù)現(xiàn)有研究以及最終實驗結(jié)果,本文在計算傾斜工況時b取1。對R(t)分別求一階導數(shù)和二階導數(shù)即可獲得汽泡生長速度和汽泡生長加速度。

3) 剪切升力

由于加熱壁面法向速度梯度的存在而導致汽泡在加熱壁面法向方向上受到剪切升力。本文采用Mei等[17]推導的剪切升力公式,如式(7)所示:

(7)

式中:Gs為剪切速率;Reb為汽泡雷諾數(shù)。

(8)

Reb=2ρ1urrb/μl

(9)

式中:ur為汽泡質(zhì)心處液相速度與汽泡速度之差,在汽泡脫離前,認為汽泡速度為0,ur即為汽泡質(zhì)心處液相速度;μl為動力黏度,Pa·s。本文采用Reichardt[18]提出的公式計算近壁面液相速度,如式(10)所示:

(10)

式中:u(x)為壁面法向距離x處的液相流速,m/s;x+為無量綱壁面距離;u*為摩擦速度;δ、c和χ為經(jīng)驗常數(shù),分別取值0.4、7.4和11。計算時,取法向距離為汽泡半徑時的近壁面液相速度。

x+和u*計算公式如式(11)所示:

(11)

式中:vl為運動黏度,m2/s;τw為壁面剪應力,kg/(m·s2)。

(12)

式中:u為主流速度,m/s;λ為摩擦系數(shù),由Blasius公式計算。

λ=0.316 4Re-0.25

(13)

4) 流動曳力

由于流體與汽泡存在速度差而引起的阻礙汽泡在流體中發(fā)生相對運動所受到的力,稱為流動曳力。采用Mei推導的公式進行計算,如式(14)所示:

Fqs=6πvlρlurrb·

(14)

5) 表面張力

根據(jù)Klausner的研究,表面張力計算公式如式(15)所示:

(15)

式中:dw為汽泡與壁面的接觸直徑,m;σ為表面張力系數(shù),N/m;α和β分別為汽泡的前、后接觸角,rad。根據(jù)Yun等[19]和Jia等[20]的研究,dw=0.134rb,α和β分別取45°和36°。

以時間為變量,利用Zuber汽泡生長模型控制汽泡尺寸,計算得到汽泡受力。根據(jù)傾斜工況下汽泡受力的不同與相關計算,依據(jù)判別標準判斷汽泡是否脫離,進而得到汽泡脫離直徑(dd)和汽泡生長時間(tg)等汽泡相關參數(shù)。

2 壁面熱流分配機理模型

在過冷沸騰傳熱特性研究過程中,已開發(fā)出許多過冷沸騰傳熱計算模型,主要包含3種壁面熱流密度的計算形式:純經(jīng)驗關系式、壁面熱流密度劃分半經(jīng)驗關系式、壁面熱流密度分配機理模型[21]。本文主要針對第3種機理模型給出傾斜加熱棒束壁面熱流分配模型。

對于傾斜條件下棒束的上加熱壁面,只存有汽泡脫離現(xiàn)象,因此總壁面熱流可分為蒸發(fā)熱流密度qme(kW/m2)、淬冷熱流密度qtc(kW/m2)和單相對流熱流密度qc(kW/m2)。

1) 蒸發(fā)熱流密度

蒸發(fā)熱流是由汽泡底部微液層和過熱液體層蒸發(fā)而產(chǎn)生的,一方面由汽泡以汽化潛熱的形式帶離加熱壁面,另一方面用于抵抗汽泡頂部冷凝效應以維持汽泡形態(tài)與生長。根據(jù)徐平昊[22]的相關研究,汽泡頂部冷凝效應相對較小并對其進行了忽略處理,Sateesh等[9]對汽泡頂部冷凝效應也進行了忽略處理,因此本文忽略汽泡頂部冷凝效應,壁面蒸發(fā)熱流密度計算如式(16)所示:

(16)

式中:dd為汽泡受力分析計算中得到的汽泡脫離直徑,m;nA為汽化核心密度,m-2;f為汽泡脫離頻率,s-1。

(17)

式中:tg為生長時間,即汽泡從產(chǎn)生到脫離的時間,s;tw為等待時間,即同一成核位點汽泡脫離后至下一汽泡出現(xiàn)的時間,s。在上述汽泡受力分析計算時,可得到汽泡生長時間,根據(jù)van Stralen等[23]的表示,即tw=3tg,Sateesh等[9]在進行模型建立時,也采用了上述相關性。

2) 淬冷熱流密度

由于汽泡脫離后,汽泡周圍冷流體填補到汽泡脫離位置,直接與加熱壁面接觸而發(fā)生瞬態(tài)導熱產(chǎn)生傳熱。根據(jù)Sateesh等[9]的研究,其計算公式如式(18)所示:

(18)

式中:ΔTb=Tw-Tl,Tl為流體溫度,℃;Ab為脫離汽泡影響面積,m2;K為經(jīng)驗常數(shù),用來對汽泡投影面積進行修正得到汽泡影響面積,在Yeoh等[24]的研究成果中K的取值在1.8～2之間,本文中K取2.0。

3) 單相對流熱流密度

在不受汽泡影響的區(qū)域,會發(fā)生單相強制對流,計算公式如式(19)所示:

qc=hc(1-Ab)ΔTb

(19)

式中,hc為單相對流傳熱系數(shù),根據(jù)Sateesh和Chuang等[25]的研究,本文采取Churchill-Chu[26]公式計算:

(20)

式中:Gr為格拉曉夫數(shù);Pr為普朗特數(shù);D為當量直徑。

因此,最終總熱流密度qs為:

qs=qme+qtc+qc

(21)

輔助汽化核心密度計算選取Chuang等[25]考慮帶有傾斜角度影響的加熱平板汽化核心密度計算公式,即:

nA=-3 400+5 470ΔTsat+3.64θinc

(22)

式中,ΔTsat=Tw-Tsat,℃。

3 計算流程

計算流程如圖3所示,給定初始時刻t0,依據(jù)Zuber汽泡生長模型,計算相關汽泡尺寸。

圖3 計算流程圖

針對傾斜工況,計算汽泡受力,依據(jù)汽泡脫離判別標準,判斷汽泡受力是否達到脫離條件,若達到,受力計算結(jié)束,輸出dd和tg,否則增加一步時間步長Δt,利用Zuber模型控制汽泡生長,繼續(xù)計算汽泡受力,直至達到汽泡脫離判別標準,輸出dd和tg。

利用汽泡受力計算結(jié)果——dd、tg和tw等參數(shù)作為壁面熱流分配模型輸入項,根據(jù)壁面熱流分配機理模型,計算各熱流密度分量,并得到傾斜工況下最終壁面熱流密度計算值。

4 結(jié)果與分析

本文利用已有實驗裝置[27],開展傾斜加熱棒束過冷流動沸騰實驗。驗證實驗數(shù)據(jù)范圍如下:傾斜角度為10°～20°;熱流密度為50.63～86.08 kW/m2;流體過冷度為7.34～18.54 ℃;壁面過熱度為7.53～11.80 ℃;實驗段連通大氣,因此實驗壓力為當?shù)卮髿鈮?質(zhì)量流速約為37.95～38.86 kg/(m2·s)。

定義評價指標,即平均相對誤差(MRE)如下:

(23)

驗證結(jié)果如圖4所示,根據(jù)以前的研究結(jié)果[28],熱流密度最大相對不確定度僅為0.31%,誤差較小?？梢钥闯?熱流密度實驗值變化趨勢與熱流密度計算值變化趨勢基本一致,93.75%數(shù)據(jù)計算誤差在±30%以內(nèi),81.25%數(shù)據(jù)計算誤差在±20%以內(nèi),平均相對誤差為13.07%,該模型對傾斜條件下棒束通道過冷流動沸騰熱流密度具備一定的計算能力。

圖4 模型計算熱流密度與實驗熱流密度比較

但是計算模型依舊存在些許偏差,可以看出,當流體過冷度較高時,計算偏差相對較大,而當流體過冷度較低時,計算偏差相對較小?？梢灶A見的是,當流體過冷度更大時,模型計算值將會更低,偏差將會更大。相同的結(jié)果也出現(xiàn)在Chuang[25]和肖波齊[29]等的實驗與模型中,這說明模型本身存在計算能力上的不足與欠缺。

為了探究模型計算值和實驗值偏差的原因,選取模型計算值良好的工況——傾斜10°、壁面過熱度10.69 ℃、熱流密度75.95 kW/m2工況作為主要參考工況,對主要參數(shù)壁面溫度進行進一步分析。計算過程中,將所有工況壁溫設置為參考工況壁溫,保持其余熱工參數(shù)不變,計算結(jié)果如圖5所示?？梢?在參考工況熱流密度之前,由于參考工況壁溫高于對應實驗工況壁溫,因此導致模型計算的熱流密度高于原熱流密度計算值;在參考工況熱流密度之后,由于參考工況壁溫低于對應實驗工況壁溫,因此導致模型計算的熱流密度低于原熱流密度計算值?？梢?壁溫變化在模型計算中起著重要作用,在實驗工況范圍內(nèi),遵從高壁溫對應高熱流密度這一物理規(guī)律。

圖5 壁溫影響分析

為進一步研究過冷沸騰區(qū)間壁溫變化趨勢對模型計算的影響,繪制了如圖6所示的典型壁溫變化趨勢圖,當熱流密度處于35 kW/m2時,壁溫變化趨勢發(fā)生轉(zhuǎn)折,同時結(jié)合實驗現(xiàn)象觀察,加熱棒壁面也開始出現(xiàn)汽泡,即認為此處為過冷沸騰起始點(ONB點)?？梢园l(fā)現(xiàn),在ONB點后,雖然壁溫變化趨勢開始由單相時的直線發(fā)生變化,但壁溫并不是立即變平,而是緩慢上升一段區(qū)間,然后再逐漸變得十分平緩。而根據(jù)圖5可知,在模型計算過程中,同一工況壁溫的升高會導致模型計算熱流密度的升高?？梢?在高過冷度過冷沸騰區(qū),壁溫依舊在緩慢上升,相對于低過冷度過冷沸騰區(qū)變化十分平緩的壁溫,整體水平要低,因此導致模型計算值較實驗值低。而在低過冷度過冷沸騰區(qū),壁溫變化十分平緩,模型計算值與實驗值符合良好。模型對于低過冷度過冷沸騰具備良好的計算能力,對于高過冷度過冷沸騰的計算仍值得進一步探究。

圖6 典型沸騰曲線

由諸多汽化核心密度公式可知,汽化核心密度與壁面過熱度有很大相關性。如圖6所示,相比于低過冷度過冷沸騰,高過冷度過冷沸騰中,在變化同等熱流密度的情況下,壁面過熱度變化程度更大。根據(jù)式(22)可知,更大程度的壁面過熱度變化導致汽化核心密度變化程度更大,而汽化核心密度對于模型計算熱流密度影響很大。選取傾斜10°、壁面過熱度10.69 ℃、熱流密度75.95 kW/m2工況作為主要參考工況對汽化核心密度影響進行進一步分析。計算過程中,保持參考工況各參數(shù)不變,在實驗工況范圍內(nèi),改變汽化核心密度進行相關計算,結(jié)果如圖7所示。可見,當汽化核心密度低于參考工況汽化核心密度時,計算熱流密度偏低。汽化核心密度低于參考工況汽化核心密度越多,計算熱流密度低于參考工況熱流密度就越多。汽化核心密度高于參考工況汽化核心密度時同理。

圖7 汽化核心密度影響分析

相比于低過冷度過冷沸騰區(qū),高過冷度過冷沸騰區(qū)壁溫變化更為明顯,且壁溫整體偏低,進而導致汽化核心密度變化更為明顯,且汽化核心密度整體偏低。而較低的汽化核心密度會導致較低的模型熱流密度計算值,造成高過冷度過冷沸騰區(qū)熱流密度計算值存在一定程度的偏低。因此,在高過冷度過冷沸騰區(qū),模型計算熱流密度偏差相對較大,而在低過冷度過冷沸騰中,模型計算熱流密度與實驗值符合較好。汽化核心密度對于壁面過熱度的依賴性對模型計算很重要,汽化核心密度與壁面過熱度之間的相關變化趨勢以及計算關系對模型計算影響很大。

5 結(jié)論

本文結(jié)合汽泡動力學與壁面熱流分配機理模型,對傾斜加熱棒束過冷沸騰壁面熱流密度進行了研究,得到如下結(jié)論。

1) 在傾斜條件下棒束上加熱壁面的汽泡受力分析過程中,根據(jù)實驗現(xiàn)象,制定“浮力方向受力打破平衡,汽泡脫離”的判定準則,從而計算得到汽泡脫離直徑dd、汽泡生長時間tg和汽泡等待時間tw等在壁面熱流分配機理模型中所需要的相關汽泡參數(shù)。

2) 對相同工況同時建立了汽泡力平衡模型和壁面熱流分配模型,并與實驗值進行了對比驗證。計算結(jié)果表明,以力平衡模型為基準,壁面熱流分配模型計算結(jié)果在一定范圍內(nèi)具有良好的計算能力,平均相對誤差為13.07%,誤差線±30%以內(nèi)涵蓋93.75%的數(shù)據(jù),±20%以內(nèi)涵蓋81.25%的數(shù)據(jù)。

3) 耦合汽泡力平衡模型的壁面熱流分配模型在低過冷度過冷沸騰區(qū)具備十分良好的計算能力;在高過冷度過冷沸騰區(qū)域,計算偏差雖然在可接受范圍之內(nèi),但相比于低過冷度過冷沸騰區(qū)偏差更大且計算值相比于實驗值偏低,這是由于過冷沸騰不同區(qū)域本身特性所導致。相比于低過冷度過冷沸騰,高過冷度過冷沸騰整體偏低且較為明顯的壁溫變化趨勢導致的汽化核心密度的改變是模型計算值偏低的主要原因。