蘇興康，顧 龍,3,*

(1.中國科學(xué)院 近代物理研究所，甘肅 蘭州 730000；2.中國科學(xué)院大學(xué) 核科學(xué)與技術(shù)學(xué)院，北京 100049；3.蘭州大學(xué) 核科學(xué)與技術(shù)學(xué)院，甘肅 蘭州 730000)

由于鉛及鉛鉍共晶合金(LBE)熱工水力性能良好、低壓下沸點(diǎn)高、化學(xué)惰性強(qiáng)，鉛冷快堆(LFR)被選為第4代核能系統(tǒng)的主力堆型之一[1]。中國加速器驅(qū)動嬗變研究裝置(CiADS)將LBE冷卻快堆作為次臨界堆堆型[2-3]。三角形排列棒束是LFR燃料組件常見的形式，研究LBE冷卻棒束的熱工水力特性對LFR安全性和經(jīng)濟(jì)性具有重要意義。在不考慮燃料棒相互作用的情況下，三角形棒束可簡化為用六邊形傳熱單元等效的環(huán)管[4-6]。因此，環(huán)管內(nèi)的湍流換熱研究有助于分析三角形棒束熱工水力。LBE的不透明、腐蝕性特點(diǎn)增加了LBE環(huán)境中傳熱實(shí)驗(yàn)的難度[7-8]，需借助計(jì)算流體力學(xué)(CFD)數(shù)值模擬LBE在環(huán)形管道內(nèi)的湍流換熱傳熱過程。然而，LBE的分子普朗特?cái)?shù)(Pr≈0.01～0.03)較低，導(dǎo)致雷諾時(shí)均方法(RANS)中的雷諾比擬假設(shè)會影響LBE的傳熱計(jì)算[9]。即用于模擬傳統(tǒng)流體(Pr≈1)的恒定湍流普朗特?cái)?shù)(Prt)，可能會影響LBE傳熱的評估[10]。在過去的40年里，為提高LBE的計(jì)算精度，國內(nèi)外研究了用于封閉能量方程的兩方程kθ-εθ模型。用于封閉動量方程的兩方程k-ε湍流模型和用于封閉能量方程的兩方程kθ-εθ傳熱模型統(tǒng)稱為四方程k-ε-kθ-εθ模型。

文獻(xiàn)[11-21]對四方程模型的系數(shù)函數(shù)和近壁熱湍流行為進(jìn)行了大量的研究，發(fā)展了多種適用于不同Pr流體的四方程模型。近年來，Manservisi等[22-24]在前人研究的基礎(chǔ)上，發(fā)展了新的四方程LBE湍流換熱模型，并利用直接數(shù)值模擬(DNS)數(shù)據(jù)和LBE傳熱實(shí)驗(yàn)數(shù)據(jù)對模型進(jìn)行了驗(yàn)證，在平板、管道、三角形和四邊形棒束通道內(nèi)的計(jì)算結(jié)果較好。文獻(xiàn)[25-30]在Manservisi四方程模型的基礎(chǔ)上，提出了比耗散率k-ω-kθ-ωθ和對數(shù)比耗散率k-Ω-kθ-Ωθ形式的四方程模型，并進(jìn)一步計(jì)算和驗(yàn)證了LBE在平面、管道、后臺階和六邊形棒束中的湍流換熱，預(yù)測精度良好。

盡管人們對低Pr流體的可靠計(jì)算工具愈加感興趣，但商業(yè)代碼仍較缺乏，基于四方程模型的LBE環(huán)管內(nèi)的熱工水力研究較少。因此，本文基于開源CFD軟件OpenFOAM，開發(fā)四方程求解器4eqnFoam?；?eqnFoam，計(jì)算LBE在平板和圓管內(nèi)的湍流換熱，并與DNS數(shù)據(jù)和實(shí)驗(yàn)關(guān)系式進(jìn)行對比驗(yàn)證。最后，基于該求解器，計(jì)算LBE在環(huán)管內(nèi)的流動傳熱，獲得不同加熱條件下LBE的關(guān)鍵參量，為計(jì)算鉛鉍流體湍流換熱提供參考。

1 湍流換熱模型

環(huán)管幾何模型示于圖1。內(nèi)、外壁面分別受恒熱流qi、qo加熱。邊界條件劃分為入口Γi、出口Γo、外壁面Γw,o和內(nèi)壁面Γw,i。Ri和Ro分別為內(nèi)外半徑。

圖1 同心環(huán)管示意圖Fig.1 Sketch of concentric annulus

對于不可壓縮、常物性、無重力的充分發(fā)展雷諾時(shí)均方程組，有：

(1)

(2)

(3)

(4)

其中，ux為沿流動方向的速度分量。

(5)

其中，k和νt分別為湍動能和湍流運(yùn)動黏度。采用引入了Kolmogorov速度尺度uε=(vε)1/4的Abe兩方程k-ε湍流模型計(jì)算νt[17]。ε為湍流耗散率。

(6)

其中，αt為湍流熱擴(kuò)散系數(shù)。選用Manservisi兩方程kθ-εθ傳熱模型計(jì)算αt[22-24]。該模型引入了湍流溫度動能kθ和湍流溫度耗散率εθ，能較好地計(jì)算LBE湍流傳熱。

2 數(shù)值求解器

OpenFOAM是基于C++編寫的開源CFD類庫[31-32]?；谏鲜鐾牧鱾鳠崮Ｐ秃蚈penFOAM，開發(fā)了4eqnFoam求解器。如圖2所示，4eqnFoam主要由動量方程UEqn.H、能量方程TEqn.H、SIMPLE壓力泊松方程pEqn.H、k-ε湍流模型VEqn.H、kθ-εθ換熱模型HEqn.H等子程序組成。前處理時(shí)，采用GAMBIT軟件劃分網(wǎng)格得到mesh文件，通過fluentMeshToFoam工具導(dǎo)入OpenFOAM中得到網(wǎng)格文件polyMesh。在前處理時(shí)，初始及邊界條件數(shù)據(jù)、網(wǎng)格及流體物性數(shù)據(jù)、離散及求解格式數(shù)據(jù)分別存放于0、constant、system文件夾中。求解器執(zhí)行計(jì)算時(shí)，前處理數(shù)據(jù)被求解器相關(guān)子程序調(diào)用。基于有限體積法離散并求解穩(wěn)態(tài)控制方程。對流項(xiàng)采用高斯迎風(fēng)格式、擴(kuò)散項(xiàng)采用高斯線性格式。矩陣選用高斯賽德爾進(jìn)行迭代求解，殘差收斂條件為：Max|Qn+1/Qn-1|<10-10，其中Q代表ui、T、p、k、ε、kθ、εθ，n為迭代步數(shù)。

圖2 4eqnFoam求解器在OpenFOAM中的框架Fig.2 Framework of 4eqnFoam solver in OpenFOAM

3 計(jì)算結(jié)果

3.1 數(shù)值驗(yàn)證(三維平板)

圖3 平板無量綱速度、無量綱雷諾應(yīng)力、無量綱溫度和無量綱雷諾熱流Fig.3 Non-dimensional velocity,non-dimensional Reynolds stress,non-dimensional temperature and non-dimensional Reynolds heat flux of plane

3.2 數(shù)值驗(yàn)證(圓管，Ri=0)

無量綱速度u+和無量綱雷諾應(yīng)力τ12沿圓管壁面無量綱法向距離y+的分布如圖4所示。由圖4a可見，工況B、D與管道壁面速度分布規(guī)律吻合良好，線性區(qū)滿足u+=y+，且對數(shù)區(qū)滿足u+=2.5lny++5，相對誤差均在4%以內(nèi)。由圖4b可見，工況A無量綱雷諾應(yīng)力分布與工況DNS-A吻合良好，相對誤差在2%以內(nèi)。

圖4 圓管無量綱速度和無量綱雷諾應(yīng)力Fig.4 Non-dimensional velocity and non-dimensional Reynolds stress of pipe

圖5 圓管無量綱溫度和無量綱溫度脈動Fig.5 Non-dimensional temperature and non-dimensional temperature fluctuation of pipe

圓管無量綱雷諾熱流〈uj′θ′〉+的分布如圖6a所示，工況A的雷諾熱流分布預(yù)測結(jié)果與DNS-A分布趨勢相近，整體預(yù)測結(jié)果高于DNS結(jié)果，最大相對誤差約為20%。努塞爾數(shù)Nu隨派克利數(shù)Pe的分布示于圖6b，將努塞爾數(shù)計(jì)算結(jié)果與Kirillov推薦的圓管傳熱關(guān)系式(Nu=4.5+0.018Pe0.8，25

3.3 環(huán)管內(nèi)流動傳熱

1)環(huán)管Ⅰ：外壁絕熱內(nèi)壁恒熱流加熱(qi=qw，qo/qi=0)

表1 環(huán)管Ⅰ計(jì)算工況Table 1 Calculation condition of annulus Ⅰ

圖7為環(huán)管Ⅰ無量綱速度u/ub和無量綱雷諾應(yīng)力τ12分別沿外壁離開到內(nèi)壁的無量綱距離y/d的分布。工況A的速度場、雷諾應(yīng)力場均與工況DNS-A吻合良好。最大速度出現(xiàn)在y/d=1.10處，零雷諾應(yīng)力出現(xiàn)在y/d=1.14附近。

圖7 環(huán)管Ⅰ無量綱速度和無量綱雷諾應(yīng)力Fig.7 Non-dimensional velocity and non-dimensional Reynolds stress in annulus Ⅰ

圖8 不同雷諾數(shù)下的環(huán)管Ⅰ無量綱溫度和無量綱溫度脈動Fig.8 Non-dimensional temperature and non-dimensional temperature fluctuation for different Re in annulus Ⅰ

圖9a、b分別為不同雷諾數(shù)下環(huán)管Ⅰ的無量綱雷諾熱流〈uj′θ′〉+和無量綱熱擴(kuò)散系數(shù)α+=αt/α(α+用于表征湍流熱流量與分子熱流量的比值)。由圖9a可見，雷諾數(shù)越大，無量綱雷諾熱流峰值越大。由圖9b可見，在近壁區(qū)，α+?1，分子導(dǎo)熱作用占主導(dǎo)地位。當(dāng)LBE的流動雷諾數(shù)較小時(shí)，全流場α+均小于1。只有當(dāng)雷諾數(shù)增加到一定程度時(shí)，如工況F，開始出現(xiàn)轉(zhuǎn)折點(diǎn)α+≥1，湍流熱擴(kuò)散作用在主流區(qū)開始強(qiáng)于LBE的分子導(dǎo)熱作用。

環(huán)管Ⅰ中的湍流普朗特?cái)?shù)Prt=vt/αt分布如圖9c所示。主流區(qū)的Prt分布相對近壁區(qū)平滑。隨著雷諾數(shù)的增加，平均湍流普朗特?cái)?shù)〈Prt〉逐漸減小，當(dāng)雷諾數(shù)增加到一定程度時(shí)，如工況D、E和F，〈Prt〉變化不明顯。此外，選取了多個外壁絕熱內(nèi)壁加熱的環(huán)形管道內(nèi)液態(tài)金屬傳熱關(guān)系式對比Nu計(jì)算結(jié)果，如文獻(xiàn)[38-42]的傳熱關(guān)系式。在低Pe下4eqnFoam預(yù)測的Nu結(jié)果與仇子鋮[39]和Jaeger[40]關(guān)系式符合良好，中等Pe下與朱鋒杰[38]和石雙凱[41]關(guān)系式符合良好，而高Pe下與張貴勤[42]傳熱關(guān)系式符合良好(圖9d)。

圖9 不同Re下環(huán)管Ⅰ的無量綱雷諾熱流、無量綱熱擴(kuò)散系數(shù)、湍流普朗特?cái)?shù)和努塞爾數(shù)Fig.9 Non-dimensional Reynolds heat flux,non-dimensional thermal diffusivity,turbulent Prandtl number and Nusselt number for different Re in annulus Ⅰ

2)環(huán)管Ⅱ：內(nèi)外壁面熱流密度相等(qi=qo=qw)

令qi=qo=qw，環(huán)管Ⅱ其他計(jì)算工況與表1一致。內(nèi)外壁面熱流密度相等時(shí)，θ+的最大值和零梯度值出現(xiàn)在最大速度值對應(yīng)位置(y/d=1.10)附近(圖10a)。相比于外壁絕熱時(shí)，無量綱溫度脈動在外壁附近多出1個峰值，整體波動變化情況更加劇烈，且兩峰之間的谷值位置均偏向內(nèi)壁(圖10b)。無量綱雷諾熱流的分布曲線和雷諾應(yīng)力分布曲線均不對稱(圖10c)。零雷諾熱流的位置出現(xiàn)在y/d=1.07附近，且不隨雷諾數(shù)變化。隨著雷諾數(shù)的增加，平均湍流普朗特?cái)?shù)〈Prt〉逐漸減小，當(dāng)雷諾數(shù)增加到一定程度時(shí)，如工況D、E和F，〈Prt〉變化不明顯(圖10d)。

圖10 不同Re下環(huán)管Ⅱ的無量綱溫度、無量綱溫度脈動、無量綱雷諾熱流和湍流普朗特?cái)?shù)Fig.10 Non-dimensional temperature,non-dimensional temperature fluctuation,non-dimensional Reynolds heat flux and turbulent Prandtl number for different Re in annulus Ⅱ

3.4 四方程模型求解器性能分析及適應(yīng)性評價(jià)

當(dāng)前，4eqnFoam求解器是基于OpenFOAM V6版本特性編寫的，編譯集成在Ubuntu 18.04系統(tǒng)環(huán)境中。串行運(yùn)行下，求解不同幾何和網(wǎng)格算例時(shí)，求解器所消耗的計(jì)算時(shí)間和內(nèi)存列于表2。測試處理器為Intel(R)Core(TM)i7-9700 CPU @ 3.00 GHz，系統(tǒng)運(yùn)行內(nèi)存為16 GB @ 2667 MHz，測試計(jì)算平臺是8核心8線程的HP ProDesk 680 G4 MT。

表2 4eqnFoam求解器性能分析Table 2 Performance analysis of 4eqnFoam solver

從平板、圓管和環(huán)管的計(jì)算驗(yàn)證結(jié)果來看，本研究開發(fā)的4eqnFoam求解器具備了初步計(jì)算不可壓縮、常物性和無重力條件下LBE湍流換熱過程的能力，并能獲得雷諾熱流、溫度脈動和湍流普朗特?cái)?shù)等結(jié)果。由于驗(yàn)證過程中只考慮了Pr=0.025常物性下LBE在簡單幾何中的對流傳熱，對于更廣范圍內(nèi)的低Pr流體、復(fù)雜幾何結(jié)構(gòu)(后臺階、棒束通道等)以及變物性條件下4eqnFoam的適應(yīng)性和計(jì)算精度，還需進(jìn)一步研究。

4 結(jié)論

為獲得LBE在環(huán)管內(nèi)不同加熱條件下的傳熱現(xiàn)象，基于開源CFD程序OpenFOAM，開發(fā)了四方程模型求解器4eqnFoam。基于4eqnFoam求解器，將LBE在平板和圓管內(nèi)的傳熱計(jì)算結(jié)果與DNS和實(shí)驗(yàn)關(guān)系式進(jìn)行了對比，驗(yàn)證了4eqnFoam的有效性。將模擬計(jì)算擴(kuò)展到LBE在內(nèi)壁加熱外壁絕熱、內(nèi)外壁面等熱流加熱的環(huán)管內(nèi)，獲得了不同雷諾數(shù)下LBE充分發(fā)展的無量綱溫度、溫度脈動、雷諾熱流、熱擴(kuò)散系數(shù)和湍流普朗特?cái)?shù)等分布。結(jié)果表明：隨著雷諾數(shù)的增加，兩種加熱條件下LBE在環(huán)管內(nèi)的無量綱溫度脈動、雷諾熱通量、湍流熱擴(kuò)散作用變大，而平均湍流普朗特?cái)?shù)減?。豢拷诿娴耐牧髌绽侍?cái)?shù)較湍流核心處的變化劇烈；當(dāng)雷諾數(shù)增加到一定程度時(shí)，平均湍流普朗特?cái)?shù)變化不明顯。基于四方程模型求解器，可為計(jì)算LBE熱工水力現(xiàn)象提供更多參考。