基于確定性抽樣的過冷沸騰邊界條件不確定性分析

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

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

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

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

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

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

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

1.1 壁面沸騰模型

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

qw=qC+qQ+qE

(1)

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

(2)

(3)

qE=VdNwρghfgf

(4)

式中：hC為單相對流換熱系數(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 相間模型

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

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

(5)

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

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

(6)

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

湍流耗散力用來描述液相湍流流動(dòng)對汽泡的耗散作用。在湍流耗散力作用下，汽泡將從近壁面區(qū)域被帶走，對徑向空泡份額分布起到展平的作用。湍流耗散力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]對豎直圓管內(nèi)R-12工質(zhì)自下而上的過冷沸騰流動(dòng)進(jìn)行了研究，簡稱DEBORA實(shí)驗(yàn)。實(shí)驗(yàn)中加熱管長3.5 m、圓管內(nèi)徑19.2 mm，實(shí)驗(yàn)測得了空泡份額、液體溫度、汽泡直徑等參數(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)無量綱數(shù)，可反映高壓水的流動(dòng)沸騰換熱特性[3]。模擬中邊界條件輸入包括質(zhì)量流速G、入口溫度T、熱流密度q和壓力p，這些變量在測量過程中由于測量誤差會(huì)引入不確定性，對模擬結(jié)果產(chǎn)生影響。此外，湍流強(qiáng)度I作為速度入口邊界中的邊界值，在充分發(fā)展段的數(shù)值計(jì)算中，若根據(jù)經(jīng)驗(yàn)假定的充分發(fā)展段管長不足以使流體完全充分發(fā)展，可能會(huì)對出口結(jié)果產(chǎn)生影響，進(jìn)而對過冷沸騰段的計(jì)算產(chǎn)生影響，因此湍流強(qiáng)度也作為不確定性源考慮。根據(jù)Fluent手冊[31]推薦和有關(guān)數(shù)值模擬邊界條件不確定性研究中湍流強(qiáng)度的選取[17,34]，模擬中假設(shè)其值為0.05。綜上，本文主要分析這5個(gè)邊界條件作為不確定性源對過冷沸騰徑向參數(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)圓管的對稱性，為節(jié)省計(jì)算資源，截取圓管的二維軸對稱平面作為計(jì)算域。采用壓力-速度耦合算法Coupled進(jìn)行穩(wěn)態(tài)計(jì)算，使用偽瞬態(tài)方法幫助加速收斂，對動(dòng)量、空泡份額、湍動(dòng)能、湍流耗散率、能量采用一階迎風(fēng)的空間離散格式，對梯度采用基于單元體的最小二乘法插值離散格式。模擬實(shí)驗(yàn)中的迭代誤差可通過計(jì)算獲得充分收斂的解來忽略。此外，為降低網(wǎng)格離散對計(jì)算結(jié)果帶來的影響，基于Case4工況針對不同網(wǎng)格(表3)進(jìn)行了網(wǎng)格敏感性分析(圖1，其中r表示徑向位置，R表示圓管半徑)?？煽吹?，mesh3和mesh4、mesh6的結(jié)果基本一致，可認(rèn)為達(dá)到網(wǎng)格無關(guān)解。綜合考慮模擬結(jié)果的準(zhǔn)確性和計(jì)算成本，選取mesh3對DEBORA實(shí)驗(yàn)進(jìn)行數(shù)值模擬，并分析其邊界條件不確定性對結(jié)果的影響。

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

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

3 不確定性分析

3.1 確定性抽樣方法

與蒙特卡羅隨機(jī)抽樣和拉丁超立方抽樣通過抽取大量樣本描述變量概率密度函數(shù)的方式不同，確定性抽樣通過傳遞變量統(tǒng)計(jì)矩的方式抽取樣本，使得樣本數(shù)大幅減少，降低計(jì)算成本[35]。根據(jù)滿足的統(tǒng)計(jì)矩階數(shù)不同，確定性抽樣具有不同抽樣形式。按照順序，前4階統(tǒng)計(jì)矩分別為均值、方差、偏度和平直度。本文采用滿足前4階矩的抽樣形式，并假設(shè)所研究的變量相互獨(dú)立且其概率密度函數(shù)滿足正態(tài)分布，可求出對于變量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ù)目。

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

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

3.2 邊界條件不確定性對過冷沸騰模擬的影響

(15)

(16)

UB=3δφ

(17)

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

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

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

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

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

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

a——空泡份額；b——液體溫度圖4 不確定性源對徑向參數(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)差和不確定度帶寬。

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

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

4 結(jié)論

本文采用確定性抽樣方法，分析了邊界條件不確定性對圓管過冷沸騰徑向參數(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)無量綱數(shù)，可反映高壓水的流動(dòng)沸騰換熱特性。此外，計(jì)算得到了每個(gè)邊界條件不確定性源的影響權(quán)重。得到如下結(jié)論：

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

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

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