基于CFD的安全殼內(nèi)浮力驅(qū)動自然對流特性分析

鄧豪放王安慶呂續(xù)艦

1（南京理工大學(xué)能源與動力工程學(xué)院 南京 210094）

2（中國電力科學(xué)研究院有限公司 北京 100192）

研究封閉空腔內(nèi)自然對流換熱現(xiàn)象，對于核反應(yīng)堆安全殼安全十分重要［1］。相關(guān)研究結(jié)果對安全殼內(nèi)氣體流動換熱計算模型的開發(fā)和推廣具有重要意義。封閉空腔內(nèi)的自然對流及換熱現(xiàn)象是傳熱學(xué)、流體力學(xué)等學(xué)科研究的重要課題之一，在核科學(xué)等工程領(lǐng)域有著廣泛的應(yīng)用，因此國內(nèi)外諸多學(xué)者都對此開展了實(shí)驗和數(shù)值模擬研究。Barakos等［2］使用控制體積法研究了差異加熱空腔內(nèi)基準(zhǔn)問題；Ampofo等［3］測量了0.75 m×0.75 m×1.5 m空腔內(nèi)自然對流實(shí)驗的基準(zhǔn)數(shù)據(jù)；堯俊等［4］利用矩形封閉空腔開展了底部弧形面自然對流換熱實(shí)驗；王煒波等［5］進(jìn)行了單側(cè)加熱矩形通道內(nèi)混合對流換熱的實(shí)驗研究。除實(shí)驗之外，研究人員還使用計算流體力學(xué)（Computational Fluid Dynamics，CFD）軟件對此進(jìn)行了數(shù)值模擬研究［6?9］。

在封閉空腔對流換熱基準(zhǔn)問題研究基礎(chǔ)上，各研究機(jī)構(gòu)通過建立安全殼模型來研究安全殼內(nèi)的傳熱現(xiàn)象及氣體分布等問題。其中德國Becker科技公司 建 立 的THAI（Thermal-hydraulics，Hydrogen，Aerosols and Iodine）系列實(shí)驗裝置有較好的代表性。TH21實(shí)驗是此系列實(shí)驗中的一個基礎(chǔ)性實(shí)驗，該實(shí)驗研究了對容器壁進(jìn)行不同溫度加熱和冷卻作用而產(chǎn)生的自然對流現(xiàn)象，為后續(xù)實(shí)驗研究奠定了基礎(chǔ)［10?11］。由于THAI實(shí)驗為相關(guān)計算程序的開發(fā)和驗證提供了大量實(shí)驗數(shù)據(jù)，國內(nèi)外諸多學(xué)者都對此進(jìn)行了大量研究。胡效明等［12］使用CFX程序?qū)HAI HM-2實(shí)驗進(jìn)行了驗證，分析了安全殼壁面的換熱和蒸汽冷凝情況下氦氣（氫氣代替工質(zhì)）的流動擴(kuò)散分布；Kelm等［13］使用OpenFOAM定制求解器“containmentFOAM”對安全殼中H2與空氣混合分布過程進(jìn)行研究，得到了比通用CFD軟件更優(yōu)的模擬結(jié)果；Zhang等［14］利用GASFLOW-MPI中三種湍流模型對安全殼內(nèi)氣體分層展開模擬研究，并與TH20實(shí)驗進(jìn)行對比，得出大渦模擬能準(zhǔn)確反映實(shí)驗現(xiàn)象；Hu等［15］利用GASFLOW-MPI程序?qū)H22實(shí)驗進(jìn)行模擬和分析，清楚地捕獲了安全殼內(nèi)氦氣的濃度和溫度分布；Dehbi等［16］利用非定常RANS模型對TH24自然循環(huán)實(shí)驗進(jìn)行計算，證明了氣體輻射對自然對流及氣體混合過程的影響作用。

由于核反應(yīng)堆事故時安全殼內(nèi)氣體行為具有較高的復(fù)雜性，氣體組分眾多，并伴隨著輻射和換熱現(xiàn)象，需要對其機(jī)理進(jìn)一步研究。為此，本文針對核反應(yīng)堆安全殼內(nèi)浮力驅(qū)動的自然對流換熱現(xiàn)象，采用CFD方法進(jìn)行數(shù)值研究。首先確認(rèn)了CFD方法對模擬封閉空腔內(nèi)自然對流和換熱現(xiàn)象的適用性和準(zhǔn)確性，隨后針對THAI-TH21模型開展計算分析，并與實(shí)驗數(shù)據(jù)進(jìn)行了比較分析，研究了安全殼內(nèi)輻射換熱及浮力驅(qū)動的自然對流現(xiàn)象。研究結(jié)果作為安全殼內(nèi)氣體復(fù)雜流動換熱的基礎(chǔ)，可進(jìn)一步分析多組分氣體的實(shí)際流動換熱現(xiàn)象。

1 計算模型

1.1 計算方法

本文基于有限體積方法，將流體域劃分成為若干個控制體單元，在每個控制體單元內(nèi)使用離散化方法求解瞬時Navier-Stokes方程（N-S方程），即質(zhì)量、能量和動量守恒方程，求解得到各點(diǎn)參數(shù)在空間內(nèi)的分布。三個守恒方程分別如下所示［17?18］。

質(zhì)量守恒方程：

動量守恒方程：

能量守恒方程：

式中：ρ為流體密度，kg?m?3；t為時間，s；u為流速，m?s-1；τˉ為粘性應(yīng)力張量，N?m?2；p為壓力，Pa；g為重力加速度，m?s-2；h為流體的比焓，J?kg-1；K為流體動能；α為導(dǎo)熱系數(shù)，W?（m?K）-1；Srad為輻射源項，W?m-3。

焓與溫度的關(guān)系為：

式中：cp為比定壓熱容，J?（kg·K）-1；T為流體溫度，K。

在本研究中，差異加熱空腔和THAI-TH21模型均采用剪切應(yīng)力SSTk-ω湍流模型，并激活Full Buoyancy Effects（浮力效應(yīng)）。壓力插值項為Body Force Weighted，干空氣介質(zhì)密度模型為Boussinesq，均采用Couple求解器開展計算。

1.2 輻射模型驗證

首先對Discrete Ordinates（DO）輻射模型計算封閉空腔內(nèi)的輻射換熱現(xiàn)象準(zhǔn)確性進(jìn)行評估。所用模型為簡單立方體模型（圖1），其中6個面的表面發(fā)射率和溫度值如表1所示。將立方體的面凈輻射熱流量（Φr）作為驗證參數(shù)，如果壁面是輻射的凈接收者，則Φr為正值；如果壁面是輻射的凈發(fā)射者，則Φr為負(fù)值。通過將多表面系統(tǒng)輻射傳熱的角系數(shù)法［17］計算得到解析解和用DO輻射模型得到的數(shù)值解進(jìn)行比較，驗證輻射模型的適用性。兩種方法得到的結(jié)果如表1所示，從中可以得出，用DO輻射模型計算得到的Φr值與解析解相比，相對誤差最大為3.79%。因此，DO輻射模型可以較好地模擬封閉空腔內(nèi)的輻射換熱現(xiàn)象。

圖1 輻射熱流計算的簡單立方體模型Fig.1 Simple cube model for calculating radiant heat flow

表1 輻射模型驗證Table 1 Radiation model validation

為進(jìn)一步探究封閉空間內(nèi)輻射效應(yīng)對換熱的影響，選擇差異加熱空腔內(nèi)自然對流的實(shí)驗基準(zhǔn)數(shù)據(jù)［3］作為模擬結(jié)果的驗證。差異加熱空腔（圖2）的熱壁面和冷壁面分別保持在323.15 K和283.15 K，頂部和底部的壁面與周圍環(huán)境絕熱。此模型用于分析二維情況下空腔內(nèi)的自然流動及輻射換熱現(xiàn)象，研究結(jié)果將應(yīng)用于后續(xù)THAI-TH21模型的模擬中。用于差異加熱空腔和TH21模型的干空氣熱物理性質(zhì)如表2所示，鋼結(jié)構(gòu)壁面和絕熱壁面的材料物性如表3所示。

表3 鋼和絕熱壁面材料特性Table 3 Material properties of steel and insulation wall

圖2 差異加熱空腔示意圖Fig.2 Schematic diagram of differentially heated cavity

表2 干空氣熱物理性質(zhì)Table 2 Material properties of dry air

1.2.1 近壁面分析處理

對于自然對流，通常有以下兩種方法來模擬近壁面區(qū)域：一是不使用壁面函數(shù)，對壁面附近網(wǎng)格進(jìn)行加密處理，完全解析近壁面邊界層；二是將第一個網(wǎng)格節(jié)點(diǎn)設(shè)置在粘性底層之外（y+>6），使用壁面函數(shù)對近壁面進(jìn)行求解。相比而言，第一種方法更為精確，但會大大增加網(wǎng)格數(shù)量，使計算成本增加；第二種方法降低了對近壁面處網(wǎng)格的要求，提高了計算效率，因此應(yīng)用較為廣泛。

本文分別采用這兩種不同方法，對y+足夠小的網(wǎng)格（y+<2），完全解析邊界層；對第一個網(wǎng)格節(jié)點(diǎn)位于粘性底層之外的網(wǎng)格（y+>2），采用標(biāo)準(zhǔn)壁面函數(shù)對近壁面進(jìn)行求解。

定義局部努塞爾數(shù)［3］如下：

式中：L為特征長度；T為溫度，K；ΔT為冷熱壁面溫差。

根據(jù)不同模擬方法得出的局部努塞爾數(shù)和實(shí)驗數(shù)據(jù)對比如圖3所示。從圖3中可知，采用標(biāo)準(zhǔn)壁面函數(shù)進(jìn)行計算時，局部努塞爾數(shù)誤差較大。因此，對于自然對流問題，若要精準(zhǔn)捕捉近壁面流動及溫度分布情況，必須保證y+<2來對近壁面邊界層進(jìn)行求解。

圖3 不同y+下局部努塞爾數(shù)對比圖Fig.3 Comparison of local Nusselt number under different y+

1.2.2 計算模型驗證

為驗證計算模型的可靠性，分別計算了差異加熱空腔一半高度處（y=0.375 m）的速度和溫度分布，并與實(shí)驗數(shù)據(jù)進(jìn)行比較，如圖4所示?？梢园l(fā)現(xiàn)，速度分布和溫度分布相對于空腔中心近似是中心對稱的，在水平方向上，兩個等溫壁面附近的溫度分布和速度分布存在著一個急劇變化。從x=0.075~0.675 m范圍內(nèi)，空腔內(nèi)在遠(yuǎn)離壁面區(qū)域流體幾乎處于靜止?fàn)顟B(tài)。在遠(yuǎn)離壁面的差異加熱空腔中部，計算結(jié)果與實(shí)驗數(shù)據(jù)吻合良好，在近壁面附近，實(shí)驗數(shù)據(jù)與數(shù)值模擬結(jié)果變化趨勢一致，數(shù)值上有一定的誤差，但誤差均在10%以下?？傮w而言，此計算方法對空腔內(nèi)的溫度及速度分布可以作出較為準(zhǔn)確的模擬，即可將此模型及計算方法用于后續(xù)安全殼的計算研究。

圖4 差異加熱空腔中部速度(a)和溫度(b)分布Fig.4 Velocity(a)and temperature(b)distribution in the middle of differentially heated cavity

1.2.3輻射作用影響分析

由于差異加熱空腔模型為絕熱封閉空腔，熱壁面上的壁面平均熱流密度值均有冷壁面上與其數(shù)值相等，符號相反的熱流密度值相對應(yīng)（圖5），此外腔體表面輻射換熱效應(yīng)會抑制自然對流，對流效應(yīng)的減弱主要由輻射效應(yīng)進(jìn)行補(bǔ)償，這與Sharma等［19］的研究結(jié)果是一致的。

圖5 差異加熱空腔冷熱壁面平均熱流密度Fig.5 The average wall heat flux of differential heating cavity

為進(jìn)一步了解輻射作用的影響，將有無輻射作用下空腔中部溫度值與實(shí)驗值進(jìn)行對比，對比圖如圖6所示，可以得出，模擬計算得到的空腔中部溫度分布結(jié)果在遠(yuǎn)離壁面處與實(shí)驗數(shù)據(jù)吻合良好，在靠近壁面處存在一定的誤差。經(jīng)分析，引起誤差的原因為邊界條件設(shè)置問題，此算例在模擬中，冷熱壁面的發(fā)射率被設(shè)置為恒定值（ε=0.3），但實(shí)驗中并未對壁面發(fā)射率作具體說明，除此之外，實(shí)驗中頂部和底部也不是模擬計算時設(shè)置的那樣完全絕熱，因此近壁面處會存在較大誤差。從有無輻射情況下封閉空腔內(nèi)的溫度分布云圖（圖7）可以看出，在這兩種情形下，封閉空腔內(nèi)部的溫度都呈層狀分布，與無輻射的情況相比，在考慮輻射作用的情況下，腔體在垂直方向上的溫差更小，空腔內(nèi)的溫度分布更加均勻。

圖6 差異加熱空腔中部溫度分布Fig.6 The temperature distribution in the middle of differentially heated cavity

圖7 差異加熱空腔溫度云圖Fig.7 The temperature contour of differentially heated cavity

2 計算結(jié)果分析

2.1 核反應(yīng)堆安全殼模型

本文采用的核反應(yīng)堆安全殼模型是THAI實(shí)驗?zāi)Ｐ停?1］，該模型具有工業(yè)規(guī)模的多隔間，常用于研究在嚴(yán)重事故期間核反應(yīng)堆安全殼中的氣體流動、溫度分布及傳熱現(xiàn)象。THAI實(shí)驗?zāi)Ｐ偷闹黧w部分是一個60 m3的不銹鋼容器，高9.2 m，直徑3.2 m，如圖8所示。裝置外部殼體裝有恒溫夾套，用于保持壁面的恒溫條件，裝置的內(nèi)部部件可以更換，可用于研究不同條件下的實(shí)驗現(xiàn)象。

圖8 THAI實(shí)驗裝置幾何參數(shù)和CFD模型Fig.8 Geometric parameters and CFD model of THAI facility

對于TH21模型，為簡化計算，忽略了內(nèi)筒的支撐結(jié)構(gòu)。在此模型中，設(shè)定加熱和冷卻壁面為等溫邊界條件，其余外壁面為混合邊界條件。THAI網(wǎng)格由六面體結(jié)構(gòu)網(wǎng)格組成，在靠近容器壁面、內(nèi)筒壁面和冷熱壁面的交界處均做加密處理（圖9），網(wǎng)格加密方法與差異加熱空腔模型壁面附近網(wǎng)格加密方法一致，以便更好捕捉壁面處的流動和溫度分布。

圖9 TH21核反應(yīng)堆模型網(wǎng)格Fig.9 Diagram of TH21 reactor model grid

在進(jìn)行TH21模型網(wǎng)格無關(guān)性研究時，如與實(shí)際模型設(shè)置保持完全一致，計算成本和時間將大幅提高。因此對計算模型進(jìn)行合理簡化，簡化后的計算模型僅設(shè)置三個不同的壁面條件：壁面1（安全殼內(nèi)表面）、壁面2（內(nèi)筒內(nèi)表面）和壁面3（內(nèi)筒外表面）。其中壁面1和壁面2溫度T=393 K，發(fā)射率ε=0.5，壁面3溫度T=367 K，發(fā)射率ε=0.1。以壁面3壁面輻射熱流量Φ為參數(shù)，來研究TH21模型的網(wǎng)格無關(guān)性。分別計算了148W、212W和302W三種網(wǎng)格數(shù)量下的壁面3平均熱流密度值Φ，結(jié)果如圖10所示。其中212W網(wǎng)格相對于148W網(wǎng)格的誤差為1.259%，302W網(wǎng)格相對于212W網(wǎng)格的誤差為0.101%，可認(rèn)為212W網(wǎng)格數(shù)量達(dá)到收斂狀態(tài)。在本研究中，TH21模型均采用212W網(wǎng)格數(shù)量進(jìn)行后續(xù)計算。

圖10 TH21模型網(wǎng)格無關(guān)性驗證Fig.10 Verification of grid-independence of TH21 model

2.2 安全殼內(nèi)溫度和壓力分布特性

在TH21實(shí)驗中［20］，安全殼內(nèi)空氣在40 500 s內(nèi)從初始大氣壓力（101 325 Pa）升高至124 000 Pa。整個實(shí)驗過程分為兩個階段：壓力和溫度上升的升溫升壓階段和壓力達(dá)到穩(wěn)定的準(zhǔn)靜態(tài)階段。在升溫升壓階段初期，數(shù)值模擬得到的壓力值與實(shí)驗值幾乎相同；在升溫升壓階段后期，數(shù)值模擬得到的壓力值逐漸高于實(shí)驗值，到達(dá)準(zhǔn)靜態(tài)階段后，計算壓力值與實(shí)驗值差值趨于穩(wěn)定。具體原因為，在升溫升壓階段初期，氣體溫度較低，壁面散熱產(chǎn)生的影響較小，隨著氣體溫度升高，實(shí)際實(shí)驗中散熱量增加，使得計算得到的壓力值逐漸高于實(shí)驗值。數(shù)值模擬得到的準(zhǔn)靜態(tài)階段壓力值為126 000 Pa，與實(shí)驗值相比，最大誤差為1.61%。整個過程壓力變化曲線與實(shí)驗值對比如圖11所示。由于實(shí)驗時實(shí)驗裝置壁面存在熱量損失，熱量損失具體數(shù)值未知，因此無法在計算時對散熱量作出精確設(shè)置，所以在模擬時未設(shè)置熱量損失，這也是模擬得到的穩(wěn)態(tài)壓力值高于實(shí)驗值的原因。從準(zhǔn)靜態(tài)階段壓力云圖（圖11）可以看出，整個容器內(nèi)壓力呈層狀分布，上高下低，氣體對安全殼上部壁面產(chǎn)生的壓力載荷較大。

圖11 壓力隨時間變化曲線及準(zhǔn)靜態(tài)階段壓力云圖Fig.11 The pressure-time curves and the pressure contour of quasi-static stage

在TH21實(shí)驗中，裝置內(nèi)氣體的流動是由浮力驅(qū)動的，y方向上的氣體密度不同是產(chǎn)生浮力的原因，密度不同是由壁面溫度不同引起的，因此應(yīng)關(guān)注TH21裝置內(nèi)不同半徑處溫度值隨高度的變化。根據(jù)實(shí)驗時裝置內(nèi)熱電偶的分布選取了兩條距中心線不同距離的豎直線來監(jiān)測溫度值分布，分別是內(nèi)筒外側(cè)和外壁面內(nèi)側(cè)的中間位置（r=1.15）和接近外壁面內(nèi)側(cè)的近壁面處（r=1.50），計算得到的這兩條直線上溫度值隨高度變化與實(shí)驗值對比曲線如圖12所示。從圖12可以得出，在靠近冷壁面的位置（y>6.6 m），CFD模擬計算結(jié)果與實(shí)驗數(shù)據(jù)吻合較好。在冷熱壁面交界高度處，由于冷熱壁面溫差引起工質(zhì)的流動摻混，此高度處溫度值產(chǎn)生“突變”。在靠近熱壁面的位置（y=1.8~6.6 m），計算結(jié)果與實(shí)驗結(jié)果溫度值分布趨勢類似，但在溫度數(shù)值上存在誤差，相對誤差均在2%以下。誤差的原因為實(shí)驗時熱電偶未對輻射效應(yīng)進(jìn)行標(biāo)定，使得熱電偶在測量有輻射作用時的溫度值出現(xiàn)了誤差。輻射效應(yīng)對熱電偶的影響主要表現(xiàn)在熱壁面附近，所以在熱壁面附近模擬得到的溫度值與實(shí)驗值存在誤差。值得注意的是，在熱壁面附近的環(huán)形空間內(nèi)，不同半徑處的溫度值出現(xiàn)了“交叉”，結(jié)合圖13速度流線圖可知原因為環(huán)形空間內(nèi)存在較弱的對流回路，使得回路附近不同半徑處的溫度值產(chǎn)生“交叉”。在模擬達(dá)到準(zhǔn)靜態(tài)階段后，處理得到z平面溫度云圖如圖12所示。從準(zhǔn)靜態(tài)階段溫度云圖中可以得出，準(zhǔn)靜態(tài)階段時TH21安全殼內(nèi)溫度接近于均勻狀態(tài)，高溫區(qū)分布在熱壁面和內(nèi)筒外壁面附近，低溫區(qū)分布在冷壁面和安全殼最底部。

圖12 不同半徑處溫度隨高度變化曲線及準(zhǔn)靜態(tài)階段溫度云圖Fig.12 Curve of temperature variation with height at different radius and temperature contour of quasi-static stage

圖13為不同高度處速度隨半徑變化曲線以及z平面速度流線圖。從速度隨變化曲線可知，在不同高度處，熱壁面附近流體流速均是同一高度處最高的。在h=3.5 m和4.9 m時，由于內(nèi)筒壁面的輻射作用，內(nèi)筒外表面的溫度要高于流體溫度，因此內(nèi)筒外壁面附近和熱壁面附近流體流速均高于環(huán)形空間內(nèi)流體流速。隨著h的增加，熱壁面對流體的加熱作用逐漸累積，空氣溫度也隨之增高，流速也逐漸升高，因此4.9 m高度處流速高于3.5 m處流速。h=6.4 m時，由于冷熱流體的交匯在交界處產(chǎn)生回路，因此在半徑0.6~1.3 m之間，y方向流速無明顯變化。從流線圖可以看出，由于壁面的溫差作用，在熱壁面附近產(chǎn)生上升流動，在冷壁面附近產(chǎn)生下降流動。熱壁面附近的上升流動和冷壁面的附近的下降流動在交界處形成回路。流場存在兩個主要的對流回路，由上升流動和下降流動交匯產(chǎn)生。在熱壁面和內(nèi)筒之間的環(huán)形空間內(nèi)，存在一些較弱的對流回路。在TH21模型的頂部和底部空間中，只存在較弱的流動現(xiàn)象，氣體流動可忽略不計。

圖13 不同高度處速度隨半徑變化曲線及準(zhǔn)靜態(tài)階段流線圖Fig.13 Curve of velocity variation with radius at different heights and stream pattern of quasi-static stage

3 結(jié)語

本文基于CFD方法對封閉空間內(nèi)的自然對流換熱現(xiàn)象以及TH21實(shí)驗進(jìn)行了計算分析，并與實(shí)驗結(jié)果進(jìn)行了對比，得到了以下主要結(jié)論：

1）對DO輻射模型的驗證和對差異加熱空腔模型分析證明了DO輻射模型對計算封閉空間內(nèi)自然對流及換熱現(xiàn)象的適用性。對差異加熱空腔的分析中得出在模擬自然對流現(xiàn)象時要保證近壁面y+<2，完全解析近壁面邊界層，并且在計算研究時應(yīng)開啟輻射模型，考慮輻射效應(yīng)對流動換熱的影響。

2）在TH21實(shí)驗中，由于壁面熱量損失數(shù)值未知而無法精確模擬，導(dǎo)致得到的穩(wěn)態(tài)壓力值高于實(shí)驗值約為1.61%。對準(zhǔn)靜態(tài)階段TH21裝置內(nèi)溫度分析得出，在冷壁面附近區(qū)域，工質(zhì)氣體（干空氣）的溫度分布與實(shí)驗數(shù)據(jù)吻合良好，在熱壁面附近計算結(jié)果與實(shí)驗數(shù)據(jù)存在誤差，最大誤差小于2%。

3）輻射效應(yīng)會使內(nèi)筒外壁面溫度升高，使得內(nèi)筒附近氣體流速增加。在熱壁面和內(nèi)筒之間的環(huán)形空間內(nèi)，存在一些較弱的對流回路，環(huán)形空間內(nèi)對流回路為同一高度不同半徑處溫度值產(chǎn)生“交叉”的原因。

作者貢獻(xiàn)聲明鄧豪放：仿真計算，初稿撰寫；王安慶：仿真模型調(diào)試；呂續(xù)艦：課題指導(dǎo)，稿件修改。