張 翔，彭敏俊，叢騰龍，李孝佳，陳衣然

(哈爾濱工程大學(xué) 核安全與仿真技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001)

過(guò)冷沸騰具有換熱系數(shù)高的特點(diǎn)，但隨著熱流密度的增加，可能達(dá)到沸騰臨界，導(dǎo)致壁溫飛升損壞設(shè)備[1]。近年來(lái)隨著計(jì)算機(jī)技術(shù)的發(fā)展和對(duì)兩相流機(jī)理的認(rèn)識(shí)，計(jì)算流體力學(xué)(CFD)方法越來(lái)越多地應(yīng)用于反應(yīng)堆熱工水力和流體流動(dòng)傳熱的研究，并已得到認(rèn)可[2]，因此，基于CFD方法的過(guò)冷沸騰流動(dòng)研究引起廣泛關(guān)注，包括基準(zhǔn)實(shí)驗(yàn)?zāi)M[3-5]、汽泡動(dòng)力學(xué)研究[6-8]、模型對(duì)比驗(yàn)證[9-11]、棒束通道流動(dòng)傳熱[12-13]等。

在上述的CFD模擬中，邊界條件、模型經(jīng)驗(yàn)系數(shù)等均作為確定值輸入，計(jì)算結(jié)果無(wú)法反映在實(shí)際邊界條件測(cè)量和模型系數(shù)簡(jiǎn)化過(guò)程中引入的不確定性帶來(lái)的影響。然而Bestion等[14]指出，在沸騰流動(dòng)的數(shù)值模擬中，邊界條件不確定性和模型不確定性的影響不可忽略。因此，需在過(guò)冷沸騰的CFD模擬中引入不確定性分析，提高結(jié)果的可信度。Badillo等[15]和Pro?ek等[16]分別利用多項(xiàng)式混沌展開(kāi)和隨機(jī)抽樣對(duì)單相GEMIX攪混實(shí)驗(yàn)的入口邊界條件進(jìn)行了不確定性分析。Dunn等[17]和Hedberg等[18]分別利用拉丁超立方抽樣和確定性抽樣對(duì)湍流模型在單相突擴(kuò)流中的不確定性進(jìn)行了分析。Zhang等[19]和Cong等[20]針對(duì)Bartolomei過(guò)冷沸騰實(shí)驗(yàn)進(jìn)行了數(shù)值模擬，并分別對(duì)入口邊界條件和沸騰模型進(jìn)行了不確定性分析，但研究中采用的拉丁超立方抽樣法(分別需要740、1 040個(gè)樣本數(shù)據(jù))產(chǎn)生的樣本數(shù)較多，計(jì)算成本較高，不適合應(yīng)用于復(fù)雜問(wèn)題的計(jì)算，且二者實(shí)驗(yàn)數(shù)據(jù)均為截面平均值在軸向的分布，均未研究不確定性傳遞對(duì)沸騰流動(dòng)參數(shù)徑向分布的影響。

本文利用Fluent軟件構(gòu)建數(shù)學(xué)物理模型，對(duì)DEBORA實(shí)驗(yàn)[21]進(jìn)行數(shù)值模擬，采用確定性抽樣方法對(duì)邊界條件輸入?yún)?shù)進(jìn)行抽樣計(jì)算，獲得主要徑向參數(shù)分布的期望值和置信區(qū)間，并與實(shí)驗(yàn)數(shù)據(jù)對(duì)比驗(yàn)證，研究不確定性對(duì)過(guò)冷沸騰流動(dòng)徑向參數(shù)的影響。之后進(jìn)一步分析得到每個(gè)不確定性源的影響權(quán)重，為模型改進(jìn)方向提供參考。

1 數(shù)學(xué)物理模型

基于歐拉兩流體模型，對(duì)過(guò)冷沸騰中的兩相分別建立質(zhì)量、動(dòng)量和能量守恒方程，結(jié)合標(biāo)準(zhǔn)k-ε湍流模型、RPI壁面沸騰模型、相間作用力以及一些輔助模型對(duì)豎直圓管內(nèi)過(guò)冷沸騰流動(dòng)建立數(shù)值模型。

1.1 壁面沸騰模型

RPI壁面沸騰模型是由Kurul等[22]提出的，用于計(jì)算壁面向流體傳遞的熱流。該模型將總熱流密度qw分為單相對(duì)流熱流密度qC、驟冷換熱熱流密度qQ和蒸發(fā)熱流密度qE3個(gè)部分，即：

qw=qC+qQ+qE

(1)

qC=hC(Tw-Tl)(1-Ab)

(2)

(3)

qE=VdNwρghfgf

(4)

式中：hC為單相對(duì)流換熱系數(shù)；Tw為壁面溫度；Tl、kl、cp,l、ρl為液體溫度、導(dǎo)熱率、比熱容、密度；ρg為氣體密度；hfg為蒸發(fā)潛熱；Vd為汽泡體積，可根據(jù)汽泡脫離直徑Dw計(jì)算；Ab、f和Nw分別為汽泡影響面積、汽泡脫離頻率和汽化核心密度。

壁面沸騰模型中選用的具體子模型如下：汽泡脫離直徑Dw，Tolubinski-Kostanchuk模型[23]；汽泡影響面積Ab，Delvalle-Kenning模型[24]；汽泡脫離頻率f，Cole模型[25]；汽化核心密度Nw，Lemmert-Chawla模型[26]。

1.2 相間模型

兩相間由于速度差等因素存在相互作用力，主要包括曳力、升力和湍流耗散力，這些相間作用力通過(guò)在動(dòng)量方程中引入源項(xiàng)來(lái)表示。

作用于汽泡上的曳力FD由兩相相對(duì)運(yùn)動(dòng)引入，主要由局部流動(dòng)參數(shù)(包括空泡份額αg、汽泡直徑dg、汽泡速度vg和液體速度vl)計(jì)算得到：

(5)

式中，CD為曳力系數(shù)，可由Schiller-Naumann模型[27]計(jì)算得到。

在豎直圓管內(nèi)過(guò)冷沸騰流動(dòng)中，由于液相在徑向的速度梯度不均勻，會(huì)對(duì)其中的汽泡產(chǎn)生升力FL作用，作用方向垂直于流動(dòng)方向，表示為：

(6)

式中，CL為升力系數(shù)，可由Moraga模型[28]計(jì)算得到。

湍流耗散力用來(lái)描述液相湍流流動(dòng)對(duì)汽泡的耗散作用。在湍流耗散力作用下，汽泡將從近壁面區(qū)域被帶走，對(duì)徑向空泡份額分布起到展平的作用。湍流耗散力FTD由Burns模型[29]計(jì)算：

(7)

式中：Kgl為相間動(dòng)量交換系數(shù)；Dl為液體耗散標(biāo)量；CTD為湍流耗散力系數(shù)；σgl為耗散普朗特?cái)?shù)。

近壁面區(qū)域的相間質(zhì)量交換率mE和主流區(qū)的相間質(zhì)量交換率m分別由式(8)、(9)計(jì)算。

(8)

(9)

式中，ql和qg分別為兩相界面到液體和汽泡的能量傳遞，可由式(10)、(11)計(jì)算。

ql=hgl(Tg-Tl)

(10)

qg=αgρgcp,g(Tsat-Tg)/δt

(11)

式中：hgl為體積傳熱系數(shù)，由Ranz-Marshall模型[30]計(jì)算；Tg和cp,g分別為氣體溫度和比熱容；δt為時(shí)間尺度，取值為0.05[31]。

2 DEBORA實(shí)驗(yàn)工況及數(shù)值建模

Garnier等[21]對(duì)豎直圓管內(nèi)R-12工質(zhì)自下而上的過(guò)冷沸騰流動(dòng)進(jìn)行了研究，簡(jiǎn)稱DEBORA實(shí)驗(yàn)。實(shí)驗(yàn)中加熱管長(zhǎng)3.5 m、圓管內(nèi)徑19.2 mm，實(shí)驗(yàn)測(cè)得了空泡份額、液體溫度、汽泡直徑等參數(shù)的徑向分布，可用于數(shù)據(jù)驗(yàn)證。為排除單個(gè)工況的計(jì)算結(jié)果巧合，本文選取了4個(gè)實(shí)驗(yàn)工況[3,32-33](表1)。實(shí)驗(yàn)工況下，2.62 MPa的R-12具有與15.7 MPa高壓水相近的沸騰相關(guān)無(wú)量綱數(shù)，可反映高壓水的流動(dòng)沸騰換熱特性[3]。模擬中邊界條件輸入包括質(zhì)量流速G、入口溫度T、熱流密度q和壓力p，這些變量在測(cè)量過(guò)程中由于測(cè)量誤差會(huì)引入不確定性，對(duì)模擬結(jié)果產(chǎn)生影響。此外，湍流強(qiáng)度I作為速度入口邊界中的邊界值，在充分發(fā)展段的數(shù)值計(jì)算中，若根據(jù)經(jīng)驗(yàn)假定的充分發(fā)展段管長(zhǎng)不足以使流體完全充分發(fā)展，可能會(huì)對(duì)出口結(jié)果產(chǎn)生影響，進(jìn)而對(duì)過(guò)冷沸騰段的計(jì)算產(chǎn)生影響，因此湍流強(qiáng)度也作為不確定性源考慮。根據(jù)Fluent手冊(cè)[31]推薦和有關(guān)數(shù)值模擬邊界條件不確定性研究中湍流強(qiáng)度的選取[17,34]，模擬中假設(shè)其值為0.05。綜上，本文主要分析這5個(gè)邊界條件作為不確定性源對(duì)過(guò)冷沸騰徑向參數(shù)計(jì)算結(jié)果的影響。假設(shè)5個(gè)邊界條件相互獨(dú)立且其概率密度函數(shù)服從正態(tài)分布，假設(shè)的輸入變量在3σ區(qū)間內(nèi)的誤差列于表2，即5個(gè)輸入變量低于此誤差的概率為99.74%。

表2 輸入變量在3σ區(qū)間內(nèi)的誤差Table 2 Error of input variable in 3σ region

考慮到實(shí)驗(yàn)圓管的對(duì)稱性，為節(jié)省計(jì)算資源，截取圓管的二維軸對(duì)稱平面作為計(jì)算域。采用壓力-速度耦合算法Coupled進(jìn)行穩(wěn)態(tài)計(jì)算，使用偽瞬態(tài)方法幫助加速收斂，對(duì)動(dòng)量、空泡份額、湍動(dòng)能、湍流耗散率、能量采用一階迎風(fēng)的空間離散格式，對(duì)梯度采用基于單元體的最小二乘法插值離散格式。模擬實(shí)驗(yàn)中的迭代誤差可通過(guò)計(jì)算獲得充分收斂的解來(lái)忽略。此外，為降低網(wǎng)格離散對(duì)計(jì)算結(jié)果帶來(lái)的影響，基于Case4工況針對(duì)不同網(wǎng)格(表3)進(jìn)行了網(wǎng)格敏感性分析(圖1，其中r表示徑向位置，R表示圓管半徑)?？煽吹剑琺esh3和mesh4、mesh6的結(jié)果基本一致，可認(rèn)為達(dá)到網(wǎng)格無(wú)關(guān)解。綜合考慮模擬結(jié)果的準(zhǔn)確性和計(jì)算成本，選取mesh3對(duì)DEBORA實(shí)驗(yàn)進(jìn)行數(shù)值模擬，并分析其邊界條件不確定性對(duì)結(jié)果的影響。

表3 網(wǎng)格信息Table 3 Mesh information

圖1 網(wǎng)格敏感性分析Fig.1 Grid sensitivity analysis

3 不確定性分析

3.1 確定性抽樣方法

與蒙特卡羅隨機(jī)抽樣和拉丁超立方抽樣通過(guò)抽取大量樣本描述變量概率密度函數(shù)的方式不同，確定性抽樣通過(guò)傳遞變量統(tǒng)計(jì)矩的方式抽取樣本，使得樣本數(shù)大幅減少，降低計(jì)算成本[35]。根據(jù)滿足的統(tǒng)計(jì)矩階數(shù)不同，確定性抽樣具有不同抽樣形式。按照順序，前4階統(tǒng)計(jì)矩分別為均值、方差、偏度和平直度。本文采用滿足前4階矩的抽樣形式，并假設(shè)所研究的變量相互獨(dú)立且其概率密度函數(shù)滿足正態(tài)分布，可求出對(duì)于變量q，滿足前4階矩的抽樣點(diǎn)分別為2個(gè)周邊點(diǎn)qi1、qi3以及1個(gè)中心點(diǎn)qi2，并賦予相應(yīng)的權(quán)重因子w[18,36]。

(12)

(13)

式中：μ和σ分別為平均值和標(biāo)準(zhǔn)差；下標(biāo)i表示第i個(gè)不確定性源。

為合并所有中心點(diǎn)，中心點(diǎn)的權(quán)重因子需重新計(jì)算，計(jì)算公式如下：

(14)

式中，N為不確定性源的數(shù)目。

通過(guò)上述方法，針對(duì)表2中列舉的帶有不確定性的5個(gè)邊界條件抽取樣本，得到的樣本及權(quán)重因子列于表4。

表4 抽樣樣本及權(quán)重Table 4 Sample and weight

3.2 邊界條件不確定性對(duì)過(guò)冷沸騰模擬的影響

(15)

(16)

UB=3δφ

(17)

圖2 邊界條件不確定性對(duì)徑向參數(shù)的影響Fig.2 Effect of boundary condition uncertainty on radial parameter

圖2中空泡份額沿壁面向管道中心方向，呈逐漸下降趨勢(shì)，這是由于流體在壁面受熱蒸發(fā)產(chǎn)生汽泡，部分汽泡在相間力的作用下向管道中心運(yùn)動(dòng)，而管道主流區(qū)域處于過(guò)冷狀態(tài)，汽泡運(yùn)動(dòng)過(guò)程中會(huì)受到冷凝作用。流道中心的空泡份額置信區(qū)間相比于壁面附近較小，說(shuō)明和壁面附近相比，管道中心空泡份額的計(jì)算受邊界條件不確定性的影響較低。相反，根據(jù)徑向3σ區(qū)間的帶寬分布來(lái)看，壁面附近氣體速度和汽泡直徑的不確定度帶寬達(dá)到最小值，說(shuō)明壁面附近氣體速度和汽泡直徑的計(jì)算對(duì)邊界條件不確定性的敏感度較低。此外，邊界條件不確定性對(duì)液體溫度和液體速度的影響在徑向上分布較為均勻，置信區(qū)間的帶寬變化不大。

壁面軸向位置的最大空泡份額是判斷臨界熱流現(xiàn)象發(fā)生的重要參數(shù)[39-40]，該參數(shù)的準(zhǔn)確計(jì)算對(duì)預(yù)測(cè)臨界熱流密度具有重要意義。然而隨著邊界條件不確定性的引入，基于確定單一輸入值的CFD方法得到壁面最大空泡份額，并將其用于預(yù)測(cè)臨界熱流密度的計(jì)算方式，不足以得到足夠可靠的預(yù)測(cè)值。在考慮邊界條件不確定性基礎(chǔ)上，計(jì)算得到了壁面最大空泡份額期望值和3σ置信區(qū)間，如圖3所示。Case1、Case2、Case3和Case4等4個(gè)工況的置信區(qū)間帶寬分別為0.087、0.089、0.090和0.089，與各自期望值相比，所占的百分比在19%～24%之間。因此，從模擬結(jié)果來(lái)看，邊界條件不確定性在壁面最大空泡份額計(jì)算中的影響不可忽略，也說(shuō)明在今后臨界熱流密度模擬預(yù)測(cè)中引入不確定性分析具有一定必要性。

圖3 邊界條件不確定性對(duì)壁面最大空泡份額的影響Fig.3 Effect of boundary condition uncertainty on maximum void fraction at wall

3.3 邊界條件不確定性源的影響權(quán)重

基于不確定性源相互獨(dú)立假設(shè)，每個(gè)不確定性源對(duì)過(guò)冷沸騰徑向參數(shù)的影響權(quán)重Ci可通過(guò)式(18)計(jì)算，計(jì)算結(jié)果如圖4所示，圖中不同徑向位置用不同符號(hào)表示，為避免結(jié)果的偶然性，針對(duì)表1中的4個(gè)工況分別進(jìn)行不確定性源的影響權(quán)重計(jì)算。

a——空泡份額；b——液體溫度圖4 不確定性源對(duì)徑向參數(shù)的影響權(quán)重Fig.4 Uncertainty source contribution on radial parameter

(18)

UBi=3δφi

(19)

(20)

式中，δφi和UBi分別為第i個(gè)不確定性源的標(biāo)準(zhǔn)差和不確定度帶寬。

對(duì)于徑向空泡份額(圖4a)，4個(gè)工況中不確定性源的影響權(quán)重計(jì)算結(jié)果基本相同，液體入口溫度的不確定性起主要影響，在3個(gè)測(cè)量位置其影響權(quán)重分別約為30%、37%和42%。另外，湍流強(qiáng)度邊界不確定性的影響權(quán)重幾乎可忽略不計(jì)，這是因?yàn)楸疚目紤]了流體的充分發(fā)展，將充分發(fā)展段出口參數(shù)作為加熱沸騰段的入口邊界，湍流強(qiáng)度邊界影響足夠小說(shuō)明加熱沸騰段的入口流體已充分發(fā)展。隨著由圓管中心到管道壁面的徑向位置的改變，可看到液體入口溫度和運(yùn)行壓力不確定性的影響權(quán)重逐漸減小，而質(zhì)量流速和壁面熱流密度不確定性的影響逐漸升高。值得注意的是，對(duì)于壁面空泡份額(r/R=1)，液體入口溫度和壁面熱流密度不確定性的影響權(quán)重均達(dá)到約30%，因此，在依據(jù)壁面最大空泡份額預(yù)測(cè)臨界熱流密度的模擬實(shí)驗(yàn)中，液體入口溫度和壁面熱流密度須謹(jǐn)慎處理。

與空泡份額不同，在不確定性源對(duì)液體溫度影響權(quán)重的計(jì)算中，不同工況的計(jì)算結(jié)果略有差異，但影響權(quán)重隨徑向位置變化的趨勢(shì)相同(圖4b)。這是由于4個(gè)工況實(shí)驗(yàn)條件不同，Case1、Case2、Case4、Case3的出口平衡含氣率依次上升(表1)。從圖4b可看出，隨著由流道中心到四周的徑向位置的改變，液體入口溫度、壁面熱流密度和質(zhì)量流速不確定性的影響權(quán)重均逐漸減小，而運(yùn)行壓力不確定性的影響權(quán)重逐漸增加。這說(shuō)明在圓管中心，液體溫度受入口流體溫度不確定性的影響較大，而靠近壁面區(qū)域的液體溫度主要受運(yùn)行壓力的影響。此外，隨著不同工況出口平衡含氣率的上升，運(yùn)行壓力不確定性對(duì)壁面附近液體溫度的影響也逐漸增大，對(duì)于4個(gè)工況中出口平衡含氣率最高的Case3，運(yùn)行壓力的影響權(quán)重甚至達(dá)到約70%。

4 結(jié)論

本文采用確定性抽樣方法，分析了邊界條件不確定性對(duì)圓管過(guò)冷沸騰徑向參數(shù)分布的影響。研究選用流動(dòng)工質(zhì)為R-12的DEBORA實(shí)驗(yàn)的4個(gè)工況作為模擬的基準(zhǔn)實(shí)驗(yàn)。實(shí)驗(yàn)工況下，2.62 MPa的R-12具有與高壓水相近的沸騰相關(guān)無(wú)量綱數(shù)，可反映高壓水的流動(dòng)沸騰換熱特性。此外，計(jì)算得到了每個(gè)邊界條件不確定性源的影響權(quán)重。得到如下結(jié)論：

1) 通過(guò)建立過(guò)冷沸騰模擬的數(shù)值模型，得到了邊界條件不確定性對(duì)徑向參數(shù)分布的影響趨勢(shì)。與實(shí)驗(yàn)值相比，本文建立的模型可較好地模擬DEBORA實(shí)驗(yàn)中的空泡份額和流體溫度徑向分布，但對(duì)于氣相速度和汽泡直徑的計(jì)算存在一定偏差。

2) 通過(guò)計(jì)算邊界條件不確定性影響下壁面最大空泡份額的期望值和置信區(qū)間，發(fā)現(xiàn)邊界條件不確定性在壁面最大空泡份額計(jì)算中的影響不可忽略。

3) 通過(guò)對(duì)比分析邊界條件不確定性源的影響權(quán)重，發(fā)現(xiàn)邊界條件中流體入口溫度的不確定性對(duì)徑向空泡份額的計(jì)算影響較大，另外越靠近壁面位置，壁面熱流密度的不確定性對(duì)空泡份額計(jì)算的影響越大。而運(yùn)行壓力和流體入口溫度的不確定性是影響徑向液體溫度計(jì)算的主要因素。