何 雯，趙陳儒，薄涵亮

(清華大學(xué) 核能與新能源技術(shù)研究院，北京 100084)

在壓水堆和一體化小型水堆的熱工設(shè)計(jì)中，不但允許堆芯冷卻劑發(fā)生過冷沸騰，還允許在少量冷卻劑通道中發(fā)生飽和沸騰，其目的是在一定的系統(tǒng)壓力下，提高壓水堆堆芯出口處的冷卻劑溫度，以及一體化小型水堆的自然循環(huán)能力，從而改善整個(gè)核電站的熱效率。因此，研究過冷流動(dòng)沸騰中的流動(dòng)和傳熱特性對(duì)反應(yīng)堆的安全運(yùn)行和經(jīng)濟(jì)性都具有重要意義[1]。然而，流動(dòng)沸騰過程非常復(fù)雜，涉及到質(zhì)量、動(dòng)量和能量的交換，且伴隨著大量氣泡的產(chǎn)生、生長(zhǎng)、脫離、運(yùn)動(dòng)和冷凝，以往關(guān)于流動(dòng)沸騰過程的研究主要集中于宏觀實(shí)驗(yàn)、理論研究和數(shù)值模擬。隨著氣泡動(dòng)力學(xué)研究的發(fā)展，揭示流動(dòng)沸騰微觀機(jī)理成為可能，而氣泡動(dòng)力學(xué)研究成果也為過冷流動(dòng)沸騰的機(jī)理揭示和理論模型發(fā)展奠定了寶貴的基礎(chǔ)。

工程上通常采用兩流體模型來描述反應(yīng)堆中的沸騰兩相流動(dòng)，能較直觀地獲得兩相的宏觀特性。李小暢[2]基于兩流體模型對(duì)壓水堆棒束子通道的兩相流動(dòng)進(jìn)行了模擬。張小英等[3]采用兩流體模型針對(duì)核反應(yīng)堆蒸汽發(fā)生器二次側(cè)的單元通道過冷段進(jìn)行了數(shù)值模擬，獲得了通道內(nèi)空泡份額、溫度、速度等參數(shù)。此外由于兩流體模型中封閉子模型較多，如壁面熱流密度分割模型、湍流模型、氣液之間界面力模型等，很多研究者對(duì)其敏感性進(jìn)行了分析。Cheung等[4]通過改變壁面熱流密度模型中的氣泡脫離直徑、脫離頻率和活化核心密度模型，耦合兩流體模型對(duì)過冷流動(dòng)沸騰進(jìn)行了數(shù)值模擬。Yeoh等[5]采用群體平衡方法(PBM)根據(jù)氣泡在主流流場(chǎng)的運(yùn)動(dòng)行為建立了氣泡密度平衡方程，耦合兩流體模型，對(duì)實(shí)際主流中存在多粒徑的氣液兩相流進(jìn)行了數(shù)值模擬。Mali等[6]對(duì)不同相界面模型和湍流模型的準(zhǔn)確度進(jìn)行了評(píng)估，然后基于兩流體模型對(duì)過冷流動(dòng)沸騰進(jìn)行了模擬。然而，盡管數(shù)值模擬應(yīng)用廣泛，但為了描述兩相界面處的質(zhì)量、能量和動(dòng)量交換，數(shù)值模擬通常需要較多的相界面方程。當(dāng)應(yīng)用到大規(guī)模工程計(jì)算時(shí)，這些方程將增加模擬所需的計(jì)算量，在計(jì)算過程中也容易導(dǎo)致發(fā)散。

基于實(shí)驗(yàn)得到的經(jīng)驗(yàn)關(guān)系式是預(yù)測(cè)過冷流動(dòng)沸騰的另一種方法。如，Cioncolini等[7]基于2 673個(gè)實(shí)驗(yàn)數(shù)據(jù)提出一種空泡份額的預(yù)測(cè)方法；Cai等[8]則基于1 118個(gè)垂直圓管的實(shí)驗(yàn)數(shù)據(jù)對(duì)現(xiàn)有關(guān)系式進(jìn)行評(píng)價(jià)，同樣提出一種預(yù)測(cè)空泡份額的經(jīng)驗(yàn)關(guān)系式。然而，這些關(guān)系式受限于所采用的實(shí)驗(yàn)數(shù)據(jù)，適用范圍有限。因此，基于理論分析獲得的機(jī)理模型是預(yù)測(cè)沸騰過程的一種更重要的方法。Levy[9]基于質(zhì)量含氣率和熱平衡值的關(guān)系，提出一種預(yù)測(cè)空泡份額的方法；Zeitoun[10]認(rèn)為當(dāng)壁面氣泡受到的蒸發(fā)作用和冷凝作用被打破時(shí)，氣泡脫離壁面進(jìn)入主流，提出了一種預(yù)測(cè)凈蒸氣產(chǎn)生點(diǎn)(OSV點(diǎn))的理論模型；Okawa[11]同樣提出一種預(yù)測(cè)OSV點(diǎn)的方法，認(rèn)為該點(diǎn)的出現(xiàn)與泡狀流向彈狀流的轉(zhuǎn)變有關(guān)。然而，盡管理論模型較經(jīng)驗(yàn)關(guān)系式具有更廣的適用范圍，但與經(jīng)驗(yàn)關(guān)系式相同，對(duì)過冷沸騰的預(yù)測(cè)都是一維流動(dòng)，即僅考慮流場(chǎng)在軸向上的質(zhì)量、能量和動(dòng)量交換，而忽略了徑向上的氣泡行為和傳熱傳質(zhì)，進(jìn)而降低了這些關(guān)系式或模型的精準(zhǔn)性。因此，有必要提出一種預(yù)測(cè)過冷流動(dòng)沸騰的新方法，該方法與一維模型相比能獲得更多的兩相參數(shù)，又能較數(shù)值模型更穩(wěn)定、精確、計(jì)算量更小。

沸騰過程中的傳熱傳質(zhì)與氣泡行為息息相關(guān)，對(duì)沸騰過程的預(yù)測(cè)離不開對(duì)氣泡行為的準(zhǔn)確描述。杜靜宇[12]基于氣泡邊界層的概念，將流場(chǎng)劃分為主流和邊界層兩個(gè)區(qū)域，基于氣泡動(dòng)力學(xué)子模型對(duì)氣泡邊界層進(jìn)行建模，重新建立了沸騰傳熱機(jī)理模型，用于沸騰曲線和臨界熱流密度的預(yù)測(cè)，并給出了氣泡從邊界層進(jìn)入主流時(shí)的速度和尺寸分布，通過數(shù)值模擬將該邊界條件與主流流體模型進(jìn)行了初步單向耦合。周毓佳[13]則主要對(duì)主流區(qū)域內(nèi)的氣泡動(dòng)力學(xué)特性進(jìn)行建模和分析，建立了泡狀流下單氣泡在流場(chǎng)中的受力和運(yùn)動(dòng)模型，并基于歐拉-拉格朗日框架，將氣泡邊界層模型單向耦合到主流流場(chǎng)的數(shù)值模擬中，對(duì)過冷流動(dòng)沸騰下主流區(qū)域的多氣泡運(yùn)動(dòng)特性進(jìn)行了分析。何雯等[14-16]基于現(xiàn)有實(shí)驗(yàn)數(shù)據(jù)，對(duì)邊界層區(qū)域內(nèi)氣泡動(dòng)力學(xué)參數(shù)的影響因素和概率分布進(jìn)行了分析，包括脫離直徑、活化核心密度和脫離頻率，對(duì)這些參數(shù)的現(xiàn)有經(jīng)驗(yàn)關(guān)系式進(jìn)行評(píng)價(jià)，并針對(duì)脫離直徑和脫離頻率提出了新的預(yù)測(cè)模型。關(guān)舜然[17]基于非均勻多尺寸模型，在歐拉-歐拉框架下對(duì)絕熱兩相流動(dòng)工況中的氣泡聚并和破碎進(jìn)行了研究，在歐拉-拉格朗日框架下，基于氣泡動(dòng)力學(xué)微觀描述和氣泡邊界層及氣泡運(yùn)動(dòng)模型，對(duì)過冷沸騰中主流氣泡的聚并、破碎對(duì)氣泡運(yùn)動(dòng)和尺寸分布的影響進(jìn)行分析，建立了綜合考慮氣泡聚并、破碎、冷凝的離散氣泡運(yùn)動(dòng)模型，并運(yùn)用于數(shù)值模擬中?？梢?，目前對(duì)于氣泡在過冷流動(dòng)沸騰中的行為已有較成熟的認(rèn)識(shí)，而將流場(chǎng)劃分為主流和邊界層兩個(gè)區(qū)域的概念不僅能簡(jiǎn)化流場(chǎng)的計(jì)算，還能獲得更多徑向上的流場(chǎng)信息。然而，兩個(gè)流場(chǎng)間的傳熱傳質(zhì)過程目前僅能通過數(shù)值模擬進(jìn)行單向耦合，兩個(gè)區(qū)域內(nèi)流場(chǎng)的傳熱傳質(zhì)和氣泡行為還有待進(jìn)一步研究。

因此，本文基于氣泡動(dòng)力學(xué)和氣泡邊界層模型，提出一種預(yù)測(cè)過冷流動(dòng)沸騰的新的理論模型，在徑向上將流場(chǎng)劃分為多個(gè)區(qū)域，采用氣液兩相流動(dòng)的分相模型描述不同區(qū)域內(nèi)的氣泡行為和流動(dòng)傳熱傳質(zhì)，利用一套準(zhǔn)二維質(zhì)量、能量和動(dòng)量方程對(duì)流場(chǎng)進(jìn)行數(shù)學(xué)建模，進(jìn)而獲得不同區(qū)域內(nèi)的流場(chǎng)信息；基于所獲得的邊界層流場(chǎng)信息，可確定ONB點(diǎn)和OSV點(diǎn)。

1 理論分析

1.1 物理模型

過冷流動(dòng)沸騰起始于單相過冷液體，沿軸向方向當(dāng)壁溫超過液體的飽和溫度時(shí)，第1個(gè)氣泡開始在壁面核化，此時(shí)，流場(chǎng)到達(dá)沸騰起始點(diǎn)(ONB點(diǎn))，該點(diǎn)為流場(chǎng)從單相轉(zhuǎn)變?yōu)閮上嗟霓D(zhuǎn)折點(diǎn)。ONB點(diǎn)之后，流場(chǎng)進(jìn)入高過冷沸騰段，此時(shí)不斷有氣泡在壁面核化并逐漸長(zhǎng)大，但由于過冷度仍較高，這些氣泡通常僅附著在壁面，既不滑移也不浮升，空泡份額也較低[8]。隨著流場(chǎng)的逐漸升溫，當(dāng)氣泡長(zhǎng)到足夠大時(shí)，受曳力的影響將離開核化點(diǎn)，部分沿著壁面滑移，部分直接離開壁面進(jìn)入主流[18]。通常將氣泡離開核化點(diǎn)的位置定義為OSV點(diǎn)，該點(diǎn)之后流場(chǎng)進(jìn)入低過冷沸騰段，此時(shí)主流內(nèi)開始出現(xiàn)氣泡，空泡份額會(huì)出現(xiàn)明顯上升，直至流場(chǎng)溫度達(dá)到飽和狀態(tài)，即Tsat點(diǎn)。由此可見，過冷流動(dòng)沸騰可劃分為高過冷沸騰段和低過冷沸騰段，不同階段流場(chǎng)情況和氣泡行為存在較大的區(qū)別，需要進(jìn)行單獨(dú)分析。

基于氣泡邊界層的思想，本文仍將流場(chǎng)劃分為主流和邊界層兩個(gè)區(qū)域，采用氣泡脫離直徑作為邊界層的厚度[12]。氣泡達(dá)到該尺寸后開始離開核化點(diǎn)，部分進(jìn)入主流，而沿著壁面滑移的氣泡的尺寸和形狀基本不變[18]，因此會(huì)形成一個(gè)附著在壁面的氣泡層。氣泡邊界層與常采用的溫度邊界層和速度邊界不同，前者定義為流體與壁面的溫度差達(dá)到主流流體與壁面的溫度差的99%處到壁面的距離，主要反映邊界層內(nèi)的對(duì)流換熱效果，后者定義為沿壁面切向的流動(dòng)速度達(dá)到自由來流速度的99%處到壁面的距離，主要反映邊界層內(nèi)動(dòng)量和質(zhì)量的損失情況。三者的相對(duì)大小受工況的影響較大，需要后續(xù)進(jìn)行深入研究。本文以恒熱流下垂直圓管內(nèi)過冷流動(dòng)沸騰為例進(jìn)行分析，各階段兩個(gè)區(qū)域內(nèi)的氣泡行為如圖1所示。

圖1 恒熱流下垂直圓管內(nèi)過冷流動(dòng)沸騰示意圖Fig.1 Schematic diagram of subcooled flow boiling in vertical heated tube under constant heat flux

準(zhǔn)確預(yù)測(cè)ONB點(diǎn)對(duì)反應(yīng)堆運(yùn)行安全至關(guān)重要，因?yàn)榭张莘蓊~的變化會(huì)對(duì)反應(yīng)性造成較大的波動(dòng)[8]。通常情況下，當(dāng)壁面附近的流場(chǎng)能滿足氣泡生長(zhǎng)所需要的能量時(shí)，ONB點(diǎn)出現(xiàn)，因此，確定ONB點(diǎn)相當(dāng)于是確定滿足氣泡生長(zhǎng)的流場(chǎng)環(huán)境。基于能量平衡，氣泡的生長(zhǎng)通常受3部分能量的影響，即過熱液層的蒸發(fā)、微液膜的蒸發(fā)和過冷流體的冷凝。Raj等[19]發(fā)現(xiàn)，當(dāng)過冷沸騰氣泡在壁面穩(wěn)定存在時(shí)，過熱液層通常小于氣泡直徑，即氣泡底端為過熱流體，而頂端仍為過冷流體，此時(shí)氣泡受到的蒸發(fā)作用和冷凝作用平衡?？紤]到流動(dòng)沸騰流場(chǎng)非常復(fù)雜，尤其是近壁面，且溫度分布很難準(zhǔn)確確定，因此，本模型假定在ONB點(diǎn)和OSV點(diǎn)處，邊界層內(nèi)溫度線性分布，當(dāng)氣泡穩(wěn)定附著在壁面時(shí)，其中心位置溫度等于飽和溫度。對(duì)于ONB點(diǎn)氣泡對(duì)應(yīng)的尺寸，受裝置測(cè)量精度的限制，通常能測(cè)到的最小氣泡尺寸為0.01 mm，因此假定ONB點(diǎn)處出現(xiàn)的第1個(gè)氣泡的尺寸為0.01 mm，中心位置溫度為TONB，而OSV點(diǎn)對(duì)應(yīng)的氣泡尺寸即為脫離直徑(Dd)。若氣泡邊界層的平均溫度為Tb，則可通過ONB點(diǎn)氣泡中心溫度TONB和壁面溫度Tw及氣泡邊界層的平均溫度Tb來確定ONB點(diǎn)，具體公式如下：

(1)

這樣，對(duì)于高過冷沸騰段，邊界層平均溫度低于飽和溫度，處于過冷狀態(tài)；對(duì)于OSV點(diǎn)，邊界層平均溫度等于飽和溫度，處于飽和狀態(tài)；對(duì)于低過冷沸騰段，邊界層平均溫度大于飽和溫度，處于過熱狀態(tài)；而在Tsat點(diǎn)，過熱度達(dá)到最大值。顯然，隨著氣泡邊界層平均溫度的逐漸上升，邊界層內(nèi)空泡份額逐漸上升。

1.2 基本假設(shè)

為簡(jiǎn)化計(jì)算，假定流場(chǎng)為穩(wěn)態(tài)流場(chǎng)且忽略截面徑向壓力波動(dòng)和軸向壓力變化對(duì)物性的影響。此外，根據(jù)前文的分析，對(duì)ONB點(diǎn)和OSV點(diǎn)重新定義，認(rèn)為兩個(gè)點(diǎn)的位置需要同時(shí)滿足氣泡尺寸分別達(dá)到0.01 mm和脫離直徑，且氣泡周圍流場(chǎng)的平均溫度等于飽和溫度兩個(gè)條件。若認(rèn)為過冷流動(dòng)沸騰段邊界層內(nèi)平均溫度和空泡份額均呈線性增加，則確定ONB點(diǎn)和OSV點(diǎn)的空泡份額以及飽和點(diǎn)對(duì)應(yīng)的平均溫度和空泡份額就非常重要。對(duì)于ONB點(diǎn)，由于氣泡數(shù)量少且尺寸小，可認(rèn)為空泡份額為0。而對(duì)于OSV點(diǎn)，本文假定邊界層空泡份額等于0.3，原因在于Okawa[11]認(rèn)為OSV點(diǎn)后空泡份額的快速上升與氣泡的聚并有關(guān)，基于現(xiàn)有的流型轉(zhuǎn)變機(jī)理，認(rèn)為此時(shí)邊界層內(nèi)氣泡呈四面體晶格排列，對(duì)應(yīng)的空泡份額為0.3。而對(duì)于飽和點(diǎn)，目前尚未見有研究對(duì)該點(diǎn)氣泡在壁面空間排列的分析，在課題組目前的研究中，通過分析不同邊界層空泡份額對(duì)模型準(zhǔn)確性的影響發(fā)現(xiàn)，邊界層空泡份額取0.5時(shí)準(zhǔn)確度最高，此時(shí)氣泡在壁面近似呈均勻的緊密排列，因此本文取該值作為飽和點(diǎn)對(duì)應(yīng)的值。此外，由于此時(shí)邊界層內(nèi)氣泡沿壁面滑動(dòng)，對(duì)流場(chǎng)進(jìn)行擾動(dòng)，此時(shí)邊界層內(nèi)溫度分布不能再假定為線性分布，而應(yīng)取對(duì)數(shù)分布，經(jīng)過敏感性分析，此時(shí)取邊界層過熱度等于壁面過熱度的1/3，準(zhǔn)確度更高，具體研究成果將在后續(xù)文章中詳述。

綜上，該模型所有的基本假設(shè)如下：1) 流場(chǎng)為穩(wěn)態(tài)流場(chǎng)；2) 截面徑向壓力波動(dòng)和軸向壓力變化對(duì)物性的影響忽略不計(jì)；3) 氣泡邊界層的厚度等于氣泡脫離直徑，氣泡僅在邊界層內(nèi)產(chǎn)生和生長(zhǎng)；4) ONB點(diǎn)的氣泡尺寸等于0.01 mm，且周圍流場(chǎng)的平均溫度等于飽和溫度；5) OSV點(diǎn)的氣泡尺寸等于脫離直徑，且周圍流場(chǎng)的平均溫度等于飽和溫度；6) 過冷沸騰邊界層內(nèi)的空泡份額呈線性上升，ONB點(diǎn)空泡份額為0，OSV點(diǎn)空泡份額為0.3，飽和點(diǎn)空泡份額為0.5；7) 過冷流動(dòng)沸騰邊界層內(nèi)的平均溫度呈線性上升，ONB點(diǎn)溫度過冷，OSV點(diǎn)溫度等于飽和溫度，飽和點(diǎn)溫度等于壁面過熱度的1/3。

2 數(shù)學(xué)描述

2.1 基本方程

由于低過冷沸騰段的氣泡行為更復(fù)雜，因此以該段為例展示數(shù)學(xué)模型的建立過程，單相段和高過冷沸騰段的數(shù)學(xué)模型可由該模型簡(jiǎn)化得到。與主流相比，邊界層內(nèi)氣液兩相的流速較小，兩者之間的速度差可忽略，再加上邊界層內(nèi)兩相的溫度一致，因此該區(qū)域可基于均相流模型描述。相反，由于主流內(nèi)氣液兩相同時(shí)存在溫度差和速度差，因此該區(qū)域需要基于分相流模型描述，被二次劃分為主流液體區(qū)域和主流氣體區(qū)域，具體如圖2所示。各區(qū)域之間由于氣泡或液體的運(yùn)動(dòng)會(huì)存在質(zhì)量、能量和動(dòng)量的交換。

圖2 低過冷流動(dòng)沸騰段數(shù)學(xué)模型Fig.2 Mathematical model for slightly subcooled flow boiling

對(duì)于質(zhì)量交換，由于氣泡在壁面不斷產(chǎn)生并攜帶部分液體進(jìn)入主流，因此從邊界層進(jìn)入主流的為氣液兩相，相反，由于氣泡進(jìn)入主流后隨流體以更高的流速運(yùn)動(dòng)，因此僅考慮液相從主流返回邊界層進(jìn)行下一輪蒸發(fā)?；诖耍瑢?duì)于邊界層區(qū)域、主流液體區(qū)域和主流氣體區(qū)域，質(zhì)量守恒方程分別如下：

M′>b-Mb=Mcb-Mbc

(2)

M′>cl-Mcl=Mbc(1-xb)-Mcb+Mcon

(3)

M′>cg-Mcg=Mbcxb-Mcon

(4)

其中：Mb、Mcl和Mcg分別為邊界層、主流液體區(qū)域和主流氣體區(qū)域的質(zhì)量流量；Mcb和Mbc分別為主流進(jìn)入邊界層的質(zhì)量流量和邊界層進(jìn)入主流的質(zhì)量流量；Mcon為冷凝質(zhì)量流量。

對(duì)于能量交換，壁面的熱量首先傳遞給邊界層，用于氣泡的蒸發(fā)，而后再通過與主流的質(zhì)量交換將能量傳遞到主流，3個(gè)區(qū)域的能量方程如下：

qwξwn+McbHcl-MbcHb-Qc=

M′>bHb-MbHb

(5)

Mbc(1-xb)Hbl-McbHcl+HcgMcon+Qc=

M′>clH′>cl-MclHcl

(6)

MbcxbHbg-HcgMcon=M′>cgHcg-McgHcg

(7)

其中：qw為熱流密度；ξw為管道周長(zhǎng)；xb為邊界層區(qū)域質(zhì)量含氣率；Hcl和Hbl分別為主流和邊界層液體的焓；Hcg和Hbg分別為主流和邊界層氣體的焓(Hcg=Hbg)；Hb為邊界層混合物的焓；Qc為邊界層和主流液體之間的導(dǎo)熱量。由于兩個(gè)區(qū)域存在溫差，以圓管為例，有：

(8)

-Ab(p′-p)+τiζin-ρbgAbn-τw,vζwn+

McbUcl-MbcUb=M′>bU′>b-MbUb

(9)

-Acl(p′-p)-τiζin+F-ρclgAcln+

Mbc(1-xb)Ub+MconUcg-McbUcl=

M′>clU′>cl-MclUcl

(10)

-Acg(p′-p)-F-ρcggAcgn+MbcxbUb-

MconUcg=M′>cgU′>cg-McgUcg

(11)

其中：p為系統(tǒng)壓力；ζi為主流與邊界層交界面的周長(zhǎng)；τi、τw,v分別為主流液體與邊界層之間、邊界層與壁面之間的應(yīng)力；F為主流氣液兩相之間的曳力。

2.2 方程求解

對(duì)于低過冷沸騰段，未知數(shù)有9個(gè)，分別為系統(tǒng)壓力(p)、主流空泡份額(αc)、3個(gè)區(qū)域的流速(Ub、Ucl、Ucg)、主流液體溫度(Tcl)、冷凝量(Mcon)以及主流和邊界層的徑向交換量(Mbc，Mcb)，其中冷凝量通過冷凝方程求得。而對(duì)于控制方程，由于Hcg=Hbg，因此方程(4)和(7)一致，即控制方程為8個(gè)，因此，在獲得3個(gè)力的條件下利用Matlab可對(duì)應(yīng)求解。然而，力的求解需要獲得交界面處的速度梯度，τi和τw,v可通過壁面速度分布公式獲得，而由于主流氣體在不斷增多且形狀變化較大，F(xiàn)的確定非常困難。因此，求解時(shí)將主流氣液兩區(qū)域的動(dòng)量方程合并，消去F，補(bǔ)充滑速比公式用于反映兩區(qū)域的流速關(guān)系，具體求解步驟可參考He等[20]的研究。而對(duì)于單相段和高過冷沸騰段，由于主流不存在氣體，該區(qū)域仍基于均相流，因此控制方程只有6個(gè)，未知數(shù)減少為7個(gè)，分別為主流和邊界層區(qū)域的流速和溫度、壓力、徑向交換量。因此，通過補(bǔ)充壁面速度分布公式可求解。

2.3 封閉方程

1) 氣泡脫離直徑

氣泡脫離直徑Dd直接影響氣泡邊界層的厚度，目前，利用無量綱數(shù)將脫離直徑各影響因素綜合考慮在內(nèi)的經(jīng)驗(yàn)關(guān)系式法被廣泛采用。杜靜宇[12]采用大量實(shí)驗(yàn)數(shù)據(jù)，對(duì)現(xiàn)有經(jīng)驗(yàn)關(guān)系式進(jìn)行評(píng)價(jià)，然后重新提出一套將力平衡結(jié)合在內(nèi)的經(jīng)驗(yàn)關(guān)系式。但何雯等[14]發(fā)現(xiàn)該關(guān)系式的準(zhǔn)確性隨著普朗克數(shù)的增大而減小，因此對(duì)該關(guān)系式進(jìn)行了修正，得到的新關(guān)系式準(zhǔn)確度更高、適用范圍也更廣。因此，本模型采用該關(guān)系式計(jì)算脫離直徑，其數(shù)學(xué)形式如式(12)所示。該公式的適用范圍為：水力直徑1～42.4 mm，壁面過熱度0.5～50.1 K，液體過冷度2～50.1 K，質(zhì)量流量67～1 927 kg/(m2·s)，工質(zhì)包括水、HFE-301和FC-87。

(12)

其中，σ、μl、cp,l、λl和ilg分別為表面張力、液體動(dòng)力黏度、液體比定壓熱容、液體導(dǎo)熱系數(shù)和汽化潛熱。

2) 滑速比公式

滑速比(S)定義為氣液兩相流速之比?；俦鹊臏y(cè)定非常困難，原因在于流動(dòng)沸騰中相變不停發(fā)生，且氣體形狀不斷發(fā)生變化，因此目前關(guān)于滑速比的研究及相應(yīng)預(yù)測(cè)模型不多[21]。由于氣液兩相的流速差通常隨著氣體含量的增多而變大，因此認(rèn)為滑速比受質(zhì)量含氣率x的影響較大，其表示單位時(shí)間內(nèi)流過通道某一截面的兩相流體總質(zhì)量中氣相所占的比例?；诖?，本文提出新的預(yù)測(cè)模型(式(13))，具體過程將在后續(xù)文章中詳述。

(13)

得到滑速比后，空泡份額和質(zhì)量含氣率可通過式(14)進(jìn)行轉(zhuǎn)換，對(duì)于邊界層，S=1。

(14)

3) 主流和邊界層平均流速

與管道尺寸相比，氣泡邊界層通常很薄，對(duì)于垂直管道內(nèi)的流動(dòng)沸騰，氣泡脫離直徑的最可幾值為0.162 mm[12]。因此，可假定這個(gè)區(qū)域內(nèi)的流體速度呈線性分布，此時(shí)邊界層區(qū)域的平均流速(Ub)為兩區(qū)域交界面處流速(Ul)的1/2[20]，即：

Ub=Ul/2

(15)

而交界面處的流速可采用適用于單相湍流的對(duì)數(shù)分布公式計(jì)算：

(16)

(17)

其中：Ac和A分別為主流區(qū)域和管道截面的面積；G為質(zhì)量流速；f為摩擦系數(shù)，是雷諾數(shù)(Re)的函數(shù)。

(18)

一旦Ub確定，則主流液體區(qū)域的平均流速可由下式計(jì)算：

Ucl=(GA-UbρbAb)/(ρclAcl+SρcgAcg)

(19)

而主流氣體區(qū)域的流速則通過滑速比公式得到。對(duì)于單相段和高過冷段，主流沒有氣體，則主流區(qū)域的平均流速為：

Uc=(GA-UbρbAb)/ρcAc

(20)

4) 冷凝公式

對(duì)于低過冷沸騰段，氣泡進(jìn)入主流后由于溫差會(huì)被冷凝，冷凝換熱量計(jì)算公式如下：

dQcon(z)=kcon(Tcg-Tcl)dz

(21)

其中：Tcg和Tcl分別為主流氣液兩相的溫度；kcon為冷凝系數(shù)，受多個(gè)參數(shù)影響，如流速、水力直徑、熱流密度等，采用Rouhani關(guān)系式[22]計(jì)算，適用范圍為1.9～23.8 MPa下的水，熱流密度為18～1 200 kW/m2。

(22)

其中：Cs為加熱周長(zhǎng)；α(z)為不同軸向位置處的截面空泡份額。一旦冷凝換熱量確定，冷凝的質(zhì)量流量Mcon可通過下式確定：

(23)

5) 黏性力

不同區(qū)域間由于速度差會(huì)產(chǎn)生黏性力，該力與流體黏性和界面處速度梯度有關(guān)，其中邊界層的速度梯度可通過壁面速度分布公式獲得：

(24)

3 模型驗(yàn)證和應(yīng)用

3.1 模型驗(yàn)證

空泡份額是流動(dòng)沸騰中的一個(gè)重要兩相參數(shù)，對(duì)換熱量、壓降和流動(dòng)不穩(wěn)定性等均有影響。因此，本文采用垂直圓管內(nèi)截面平均空泡份額的實(shí)驗(yàn)數(shù)據(jù)對(duì)模型進(jìn)行驗(yàn)證，工質(zhì)為水，由于低過冷沸騰段空泡份額更明顯且呈快速增長(zhǎng)，因此主要對(duì)這一階段的空泡份額進(jìn)行對(duì)比。Cai等[8]曾對(duì)預(yù)測(cè)過冷流動(dòng)沸騰的現(xiàn)有經(jīng)驗(yàn)關(guān)系式進(jìn)行了對(duì)比和評(píng)價(jià)，包括OSV點(diǎn)的確定、質(zhì)量含氣率的計(jì)算以及質(zhì)量含氣率與空泡份額的轉(zhuǎn)換關(guān)系等，基于大量實(shí)驗(yàn)數(shù)據(jù)得到一套用于預(yù)測(cè)垂直圓管內(nèi)過冷流動(dòng)沸騰的經(jīng)驗(yàn)關(guān)系式的組合，因此，該模型也作為本模型的對(duì)比之一。為定量反映模型的準(zhǔn)確性，定義絕對(duì)誤差(MAE)如下：

(25)

圖3為模型與Ferrell[23]和Rouhani[24]實(shí)驗(yàn)數(shù)據(jù)的對(duì)比，其中橫坐標(biāo)為熱平衡干度(xth=(h-hl,sat)/ilg)，也稱為熱平衡含氣率。可發(fā)現(xiàn)，本模型能對(duì)大多數(shù)工況下的空泡份額進(jìn)行準(zhǔn)確預(yù)測(cè)，MAE分別為29.1%和17.7%。與Cai模型[8]相比，本模型與實(shí)驗(yàn)數(shù)據(jù)的增長(zhǎng)趨勢(shì)更相近，且準(zhǔn)確度更高。原因在于，Cai模型假定OSV點(diǎn)處空泡份額為0，即忽略此時(shí)邊界層包含的氣泡數(shù)量，因此會(huì)造成模型在起始段與實(shí)驗(yàn)數(shù)據(jù)存在差異，此外，Cai模型采用經(jīng)驗(yàn)關(guān)系式確定OSV點(diǎn)的具體位置，該關(guān)系式的誤差同樣會(huì)對(duì)模型的預(yù)測(cè)結(jié)果造成較大的影響。圖4為模型預(yù)測(cè)與Mali等[6]實(shí)驗(yàn)中空泡份額和流體溫度的準(zhǔn)確性對(duì)比。圖4顯示，模型對(duì)流體溫度預(yù)測(cè)的準(zhǔn)確度高，MAE為0.3%，對(duì)空泡份額預(yù)測(cè)的準(zhǔn)確性雖不如Cai模型，但仍與實(shí)驗(yàn)數(shù)據(jù)相差不大，MAE為15.8%。

圖3 低過冷流動(dòng)沸騰段空泡份額預(yù)測(cè)值與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比Fig.3 Comparisons of predicted void fraction with experimental data in slightly subcooled flow boiling

圖4 低過冷流動(dòng)沸騰段空泡份額和液體溫度預(yù)測(cè)值與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比Fig.4 Comparisons of predicted void fraction and temperature with experimental data in slightly subcooled flow boiling

綜上，通過與垂直圓管內(nèi)空泡份額和流體溫度實(shí)驗(yàn)數(shù)據(jù)的對(duì)比，可證明本模型在預(yù)測(cè)過冷流動(dòng)沸騰中具有較高的準(zhǔn)確性，空泡份額整體誤差為22.3%。由于實(shí)驗(yàn)數(shù)據(jù)有限，本文驗(yàn)證的工況如下：壓力為0.827～4.5 MPa的水、質(zhì)量流速為520～1 440 kg/(m2·s)、熱流密度為243～888 kW/m2、水力直徑為6.1～15.4 mm。盡管驗(yàn)證范圍有限，但本模型主要依靠機(jī)理分析，普適性很廣，本模型的適用性主要受封閉方程中經(jīng)驗(yàn)關(guān)系式的適用范圍影響，但影響有限，因此可認(rèn)為本模型有應(yīng)用到更大工況范圍的潛力。后續(xù)基于更多的實(shí)驗(yàn)數(shù)據(jù)可對(duì)這些經(jīng)驗(yàn)關(guān)系式以及本模型的適用性進(jìn)行進(jìn)一步確定和評(píng)價(jià)。

3.2 模型應(yīng)用

為展現(xiàn)模型在工程應(yīng)用上的意義，本文以某一體化自然循環(huán)小型水堆燃料元件加熱通道工質(zhì)為例，將其應(yīng)用到反應(yīng)堆堆芯出口溫度恰好等于飽和溫度的一種微沸騰工況，預(yù)測(cè)其過冷流動(dòng)沸騰的狀態(tài)，其棒束結(jié)構(gòu)如圖5所示，水自下向上流動(dòng)，流速為0.5 m/s，系統(tǒng)壓力為6.5 MPa，對(duì)應(yīng)流體的飽和溫度為554 K，流體入口過冷度為35 K，熱流密度為189 kW/m2，管道長(zhǎng)度為2.1 m，由于流場(chǎng)對(duì)稱，因此選擇其中1/4柵格進(jìn)行計(jì)算。圖6為兩相參數(shù)沿通道方向的變化，需要說明的是，對(duì)于該工況，在低過冷沸騰開始階段，由冷凝公式計(jì)算得到的冷凝量大于邊界層進(jìn)入主流的氣體含量，造成起始階段主流空泡份額為負(fù)，而出現(xiàn)該現(xiàn)象的原因可能來自冷凝公式的誤差。為避免該現(xiàn)象的出現(xiàn)，計(jì)算該工況時(shí)對(duì)冷凝公式乘以0.3進(jìn)行修正，而在不同工況下的修正情況將在后續(xù)文章中進(jìn)一步討論。結(jié)果顯示，整個(gè)過冷沸騰段的長(zhǎng)度占總長(zhǎng)度的52.6%，ONB點(diǎn)和OSV點(diǎn)的位置分別距離通道出口1 764 mm和1 300 mm，這兩個(gè)點(diǎn)的確定對(duì)燃料元件通道內(nèi)流場(chǎng)的分析，如不同運(yùn)行工況的影響、反應(yīng)性的影響、自然循環(huán)能力的增強(qiáng)等都非常重要，對(duì)反應(yīng)堆微沸騰工況安全運(yùn)行和經(jīng)濟(jì)性也具有重要意義。此外，3個(gè)區(qū)域的空泡份額、流速以及系統(tǒng)壓力的變化也分別示于圖6b～d。值得說明的是，邊界層流速在OSV點(diǎn)前后出現(xiàn)明顯變化，原因在于OSV點(diǎn)前模型采用壁面速度分布公式獲得邊界層流速，并假定此時(shí)邊界層內(nèi)速度呈線性分布，而OSV點(diǎn)后由于控制方程的增多，模型不再采用該公式和假設(shè)。

圖5 燃料元件通道內(nèi)的管束結(jié)構(gòu)及流道Fig.5 Tube bundle structure and flow passage in fuel element channel

圖6 燃料元件通道內(nèi)過冷流動(dòng)沸騰兩相參數(shù)的變化Fig.6 Variations of two-phase flow parameter in subcooled flow boiling in fuel element channel

4 結(jié)論

本文基于氣泡動(dòng)力學(xué)理論提出了一套用于預(yù)測(cè)過冷流動(dòng)沸騰的邊界層模型，該模型在徑向上將流場(chǎng)劃分為多個(gè)區(qū)域，通過一組準(zhǔn)二維質(zhì)量、能量和動(dòng)量方程，將不同區(qū)域內(nèi)的氣泡行為和區(qū)域間的傳熱傳質(zhì)考慮在內(nèi)，其中主流區(qū)域主要基于分相流模型，邊界層基于均相流模型。利用獲得的邊界層流場(chǎng)信息，本文提出了一種確定ONB點(diǎn)和OSV點(diǎn)的新方法，與垂直管道內(nèi)的空泡份額和流體溫度的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了驗(yàn)證對(duì)比，絕對(duì)誤差分別為22.3%和0.3%，并成功應(yīng)用于某一體化自然循環(huán)小型水堆燃料元件通道流體的計(jì)算。

隨著研究的進(jìn)一步深入，本模型將考慮更多的氣泡動(dòng)力學(xué)過程，并將推廣至飽和沸騰，獲得更大范圍流動(dòng)沸騰的計(jì)算分析，可為壓水堆和一體化小型水堆的熱工設(shè)計(jì)和安全分析提供有效的計(jì)算分析工具。