基于四方程模型的液態(tài)鉛鉍共晶合金后臺(tái)階湍流混合對(duì)流數(shù)值模擬

中圖分類號(hào)：TL334 文獻(xiàn)標(biāo)志碼：A DOI： 10.19907/j.0490-6756.240255

Numerical simulation of turbulent mixed convection over a backward-facing step for liquid lead-bismush eutectic based on four-equation model

DONGWen-Meng，ZHANGKe-Fan，CHENHong-Li，CHENZhao2 （1.SchoolofNuclear Science and Technology，Universityof Science and Technology of China，Hefei 230027， China；2.China Nuclear Power Technology Research Institute Co.，Ltd.，Shenzhen 518028，China）

Abstract： Liquid lead-bismuth eutectic （LBE）has important applications in next-generation reactors due to its excellnt physical properties，but its low Prandtl number invalidates the similarity assumption between velocity and temperature fields in turbulent numerical simulations.Inthis work，afour-equation model wasused to simulate a vertical backward-facing step flow.By applying a uniform heat flux at the wall，the flow ofLBE under forced and mixed convection was studied.In allthecases，the Reynolds number is maintained at 48O5 and the expansion ratio is fixed at1.5.The calculation results were compared with Direct Numerical Simulation（DNS）results， and the wallfriction coeffcient，velocity，turbulent kinetic energy，Reynolds stress，temperature，and Nusselt number （ Nu） were provided for Richardson numbers （ Ri） of O and O.2. The results show that buoyancy exerts asignificant influence on the LBE flow field.For buoyancy-assisted LBE flow，the size of the recirculation region isreduced，and the reattachment length is shortened. In the case of Ri=0.2 ，the corner vortex dominates the main vortex，causing the main vortex to detach from the wall and resulting in a positive surface friction coeficient.The obtained results align well with direct numerical simulation results，validating that the four-equation model can accurately simulate the mixed convection of liquid LBE.

Keywords：Liquid lead-bismuth eutectic；Four-equation model;Mixed convection；Backward-facing step

1引言

以煤炭為主的傳統(tǒng)化石燃料能源結(jié)構(gòu)現(xiàn)已無(wú)法適應(yīng)現(xiàn)行經(jīng)濟(jì)的快速發(fā)展，并且產(chǎn)生了嚴(yán)重的環(huán)境問題.綜合考慮現(xiàn)有的煤炭?jī)?chǔ)量以及尚未被發(fā)現(xiàn)的煤炭資源，總體來(lái)講可供人類使用約250年.同時(shí)，燃燒化石燃料會(huì)向空氣中排放大量的有害氣體以及煙塵，這些有害氣體與煙塵會(huì)對(duì)環(huán)境造成嚴(yán)重的破壞.因此，發(fā)展新型清潔能源是非常有必要的，目前，水電的電力不夠穩(wěn)定，太陽(yáng)能和風(fēng)能受限于地理位置，無(wú)法大規(guī)模推廣應(yīng)用，因此無(wú)法短時(shí)間內(nèi)在總電力裝機(jī)容量中占有較大的份額.在剩下的清潔能源中，核電是目前唯一能大規(guī)模工業(yè)應(yīng)用并替代化石燃料的能源.

截至2018年7月，第4代先進(jìn)核能系統(tǒng)國(guó)際論壇組織（GenerationIVInternationalForum，GIF）確定了6種第4代先進(jìn)核能技術(shù)的候選堆型，其中的2種快中子反應(yīng)堆使用鈉、鉛、鉛鉍共晶合金（Lead-BismuthEutectic，LBE）等液態(tài)金屬作為冷卻劑．這主要是因?yàn)榭熘凶臃磻?yīng)堆的堆芯體積更小，堆芯功率密度更大，也因此對(duì)冷卻劑快速換熱載熱能力的要求更高，使得熱堆中常用的冷卻劑已不能滿足要求，而液態(tài)金屬以其優(yōu)良的傳熱特性成為快堆中首選的冷卻劑.除上述領(lǐng)域外，液態(tài)金屬在半導(dǎo)體生產(chǎn)和太陽(yáng)能熱電等領(lǐng)域均有著非常廣泛的應(yīng)用前景.因此，對(duì)液態(tài)金屬開展全面深入的研究十分必要．

液態(tài)金屬在對(duì)流換熱方面與水、空氣等普通流體區(qū)別很大.對(duì)于普朗特?cái)?shù) （Pr） 約為1的水、空氣等流體，一般可以認(rèn)為這些流體在壁面附近的速度邊界層和溫度邊界層厚度近似相等.由于速度邊界層內(nèi)部可以認(rèn)為是層流狀態(tài)，因此分子熱傳導(dǎo)作用也主要在層流區(qū)域內(nèi)；而在遠(yuǎn)離速度邊界層的湍流區(qū)域內(nèi)，分子熱傳導(dǎo)的作用較弱，以湍流擴(kuò)散作用為主.相比之下，液態(tài)金屬由于 Pr 較低，因此可以得出其溫度邊界層要厚于速度邊界層，使得分子熱傳導(dǎo)作用同時(shí)影響層流區(qū)域與湍流區(qū)域，其湍流換熱表現(xiàn)與傳統(tǒng)流體差異很大.這引發(fā)的問題在于：當(dāng)采用雷諾時(shí)均法（ReynoldsAveragedNavier-Stokes，RANS）對(duì)流體進(jìn)行數(shù)值模擬時(shí)，通常需借助湍流普朗特?cái)?shù) （Pr_t） 來(lái)求解湍流熱流項(xiàng).該參數(shù)用于表征速度邊界層與溫度邊界層的厚度比值．對(duì)于傳統(tǒng)流體，采用 Pr_t= 0.85的假設(shè)可以很好地滿足計(jì)算要求；而對(duì)于液態(tài)金屬，因其速度場(chǎng)與溫度場(chǎng)的不相似性，該假設(shè)不再適用.現(xiàn)有數(shù)值研究表明，液態(tài)金屬的 Pr 并非常數(shù)，其值會(huì)隨幾何構(gòu)型、流動(dòng)工況等因素發(fā)生變化.因此在使用RANS方法計(jì)算液態(tài)金屬的湍流換熱時(shí)，需要使用更復(fù)雜的模型[1-3].基于簡(jiǎn)單梯度擴(kuò)散假設(shè)（Simple GradientDiffusion Hypoth-esis，SGDH），將雷諾熱流 轉(zhuǎn)化為湍流熱擴(kuò)散系數(shù) α_t 的計(jì)算，建立溫度脈動(dòng) k_θ 及其耗散率 ε_θ 兩方程模型來(lái)計(jì)算 α_ι 這種使用 k^-ε^-k_θ^-ε_θ 的形式同時(shí)計(jì)算湍流運(yùn)動(dòng)黏度與湍流熱擴(kuò)散系數(shù)的模型稱為四方程模型．此外，低 Pr 流體更容易受到浮力的影響，使情況變得復(fù)雜[4].

本文使用四方程模型研究垂直后臺(tái)階上液態(tài)LBE的湍流混合對(duì)流行為，一個(gè)典型的后臺(tái)階幾何結(jié)構(gòu)如圖1所示，其包含了復(fù)雜的流動(dòng)分離和再附著等流動(dòng).后臺(tái)階作為一種重要的流動(dòng)力學(xué)模型已被研究了半個(gè)多世紀(jì)[5].針對(duì) Pr?1 的液態(tài)金屬開展的數(shù)值研究并不多，首先針對(duì)液態(tài)金屬在后臺(tái)階中的行為開展研究的是Niemann和Frohlich[6-8].他們選取液態(tài)鈉 （Pr=0.0088） 作為工作流體，在入口流動(dòng)雷諾數(shù) Re_h=4805 的條件下使用直接數(shù)值模擬（DirectNumericalSimulation，DNS）方法模擬了強(qiáng)制對(duì)流和混合對(duì)流，分析了浮力對(duì)于換熱、再附著區(qū)等方面的影響，獲得了重要的數(shù)據(jù).文獻(xiàn)[9，10]分別使用代數(shù)熱通量模型與四方程模型，就液態(tài)鈉在后臺(tái)階上的混合對(duì)流行為與DNS數(shù)據(jù)進(jìn)行了對(duì)比驗(yàn)證，其結(jié)果與DNS數(shù)據(jù)均擬合得較好．文獻(xiàn)［11]使用DNS方法研究了以LBE為代表的液態(tài)金屬在不同幾何和浮力條件下的換熱過程，幾何結(jié)構(gòu)選取圓管和后臺(tái)階，提供的精確數(shù)據(jù)可供后期不同的湍流熱通量模型驗(yàn)證和開發(fā).盡管液態(tài)鈉與液態(tài)LBE同屬于液態(tài)金屬，但是鈉屬于堿金屬，其密度和熔點(diǎn)相對(duì)較低，傳熱能力更優(yōu)秀；而液態(tài)LBE屬于重金屬，相比于堿金屬類別，其密度和沸點(diǎn)相對(duì)更高，傳熱能力弱一些，具有更強(qiáng)的腐蝕性.因此有必要針對(duì)液態(tài)LBE開展相應(yīng)的數(shù)值模擬研究.

本文采用四方程模型，將與湍流熱擴(kuò)散系數(shù)相關(guān)的浮力項(xiàng)引人k-ωSST（ShearStressTransport）模型，計(jì)算了液態(tài)LBE在后臺(tái)階上的換熱數(shù)據(jù)，并將計(jì)算結(jié)果與DNS數(shù)據(jù)進(jìn)行了比較[12].

2 湍流模型

使用RANS方法模擬液態(tài)LBE在垂直后臺(tái)階上的湍流流動(dòng)，基于不可壓縮流假設(shè)和Boussinesq近似來(lái)模擬浮力，控制方程如下.

其中， u_i 是速度分量， g_i 是重力項(xiàng)， u 是運(yùn)動(dòng)黏度， α 是熱擴(kuò)散系數(shù).

根據(jù)渦黏性模型（EddyViscosityModel，EVM）以及簡(jiǎn)單梯度擴(kuò)散假設(shè)，雷諾應(yīng)力項(xiàng)以及雷諾熱流項(xiàng)的表達(dá)式如下所示．

四方程模型主要由計(jì)算雷諾應(yīng)力的湍流模型和計(jì)算雷諾熱流密度的換熱模型組成.湍流模型使用浮力修正的 kω SST模型.kωSST模型由Menter13提出，是一種被廣泛應(yīng)用的模型.標(biāo)準(zhǔn)k^-ω 模型是低雷諾數(shù) Re 模型，適用于低雷諾數(shù)的流動(dòng)，即邊界層較厚并可以解析粘性子層.e模型在自由流動(dòng)區(qū)域內(nèi)表現(xiàn)較好，但是在邊界層附近需要使用經(jīng)驗(yàn)函數(shù).這些函數(shù)本質(zhì)上是針對(duì)平板邊界層流動(dòng)得出的，因此在一些結(jié)構(gòu)中不太準(zhǔn)確．而k^-ω 模型不需要這些函數(shù)，從而可以提供更好的精度.另一方面， k-ω 模型的常見問題是對(duì) k，ω 等湍流參數(shù)的初始值過于敏感，而 k^-ε 模型在這些方面不那么容易受影響， k^-ω SST模型通過混合函數(shù)來(lái)使得模型可以在二者中自由切換，在距壁面較近的位置使用比耗散率模型，而在距壁面較遠(yuǎn)的位置使用耗散率模型.通過這種方式，使得 k-ω SST模型不需要邊界層函數(shù).其公式如下.

其中， k 是湍動(dòng)能， 是湍動(dòng)能產(chǎn)生項(xiàng)， s 是渦量，y 是距壁面的距離， F₁ 是混合函數(shù)， G_b 是由于浮力影響引起的湍動(dòng)能產(chǎn)生項(xiàng)[14]， β_T 是液態(tài)LBE的熱膨脹系數(shù)， α_t 是湍流熱擴(kuò)散系數(shù)，默認(rèn)模型系數(shù)如表1所示.

φ=F₁φ₁+（1-F₁）φ₂

湍流換熱模型如下[15]，湍流熱擴(kuò)散系數(shù) α_t 的模型通常建立在已有的湍流運(yùn)動(dòng)黏度 u_t 模型的基礎(chǔ)上，同時(shí)引入分子普朗特?cái)?shù) Pr 效應(yīng)和混合時(shí)間尺度 R

其中， R 是熱湍流時(shí)間尺度 τ_θ 與動(dòng)力時(shí)間尺度 τ_u 之比，即 R=τ_θ/τ_u ；熱湍流時(shí)間尺度 τ_θ=k_θ/ε_θ;R_t 為湍流雷諾數(shù)： R_δ 為特征雷諾數(shù).

為了獲得溫度脈動(dòng) k_θ 及其耗散率 ε_θ 的分布，還需要對(duì)其建立微分輸運(yùn)方程.方程模型系數(shù)如表2所示．

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

為便于本文結(jié)果與DNS數(shù)據(jù)的比較，本文采取了與DNS研究相同的后臺(tái)階幾何結(jié)構(gòu).該結(jié)構(gòu)參考了文獻(xiàn)[16]的后臺(tái)階實(shí)驗(yàn)結(jié)構(gòu)，如圖2所示．圖2為一個(gè)垂直的后臺(tái)階結(jié)構(gòu)，其中臺(tái)階高度為 h 入口在臺(tái)階上游2h處，入口高度為2h，出口高度為3h，膨脹比為1.5，臺(tái)階后加熱段長(zhǎng)度為 20h，h= 0.04m ，重力沿 x 軸負(fù)方向，本文設(shè)置計(jì)算域?yàn)槎S，因此展向?qū)挾葹镺.參考DNS中的處理方式，在入口處給定充分發(fā)展的湍流入口條件和固定入口溫度，上下壁面使用無(wú)滑移壁面條件，除下壁面加熱段給定恒熱流密度外，其余均被視為絕熱墻.

工作流體為液態(tài)LBE，溫度設(shè)置為546.3K，即入口溫度，在這個(gè)溫度下，LBE的普朗特?cái)?shù)為Pr=0.025 除了普朗特?cái)?shù)外，定義本算例的其他2個(gè)無(wú)量綱數(shù)為雷諾數(shù)和理查森數(shù).雷諾數(shù)基于人□平均速度 U_b 和臺(tái)階高度 h 計(jì)算， Re_b=U_bh/ν= 4805.理查森數(shù)定義為 Ri=Gr/Re_b²=gβΔT/U_b² 特征溫差 ΔT 由加熱段施加的恒熱流 定義， ，其中λ是LBE的導(dǎo)熱系數(shù)．為與DNS數(shù)據(jù)進(jìn)行對(duì)比，本文研究了 Ri=0 的強(qiáng)迫對(duì)流情況與 Ri=0.2 的混合對(duì)流情況，其中在研究混合對(duì)流時(shí)額外在出口后增加 10h 的區(qū)域以減少出口對(duì)流場(chǎng)的影響.

4 結(jié)果與分析

4.1 網(wǎng)格無(wú)關(guān)性分析

首先驗(yàn)證后臺(tái)階的網(wǎng)格無(wú)關(guān)性，如表3所示，本文使用3種不同數(shù)量的網(wǎng)格.強(qiáng)制對(duì)流下，后臺(tái)階的一個(gè)典型特征是發(fā)生在臺(tái)階后的流動(dòng)分離和再附著，使用再附著點(diǎn)的流向位置 評(píng)估數(shù)值結(jié)果對(duì)網(wǎng)格的敏感性，并給出了不同網(wǎng)格結(jié)果相對(duì)于最精細(xì)網(wǎng)格的相對(duì)誤差.可以看出，Mesh2與Mesh3的結(jié)果相差小于 1% ，因此，本工作中進(jìn)行的計(jì)算采用Mesh2.

4.2流場(chǎng)

4.2.1強(qiáng)迫對(duì)流在強(qiáng)迫對(duì)流情況下，計(jì)算得到的平均流向速度云圖如圖3所示，范圍從 x/h= -2.0 至 x/h=20 ：圖4給出了臺(tái)階附近的流線分布圖，可以觀察到RANS模擬較好地再現(xiàn)了后臺(tái)階構(gòu)型流動(dòng)的典型特征，即在臺(tái)階后發(fā)生流動(dòng)分離與流動(dòng)重新附著，并在下游進(jìn)一步發(fā)展成附著示.曲線有2個(gè)明顯的與零線的交點(diǎn)，其中由負(fù)到正的交點(diǎn)稱為再附著點(diǎn)，表示主流重新接觸壁面的位置.計(jì)算得到的再附著點(diǎn)長(zhǎng)度為7.50，DNS得到的長(zhǎng)度為7.01.在回流區(qū)內(nèi)，臺(tái)階和壁面之間存在一個(gè)較小的逆時(shí)針再循環(huán)區(qū)域，稱為角渦.在DNS計(jì)算中，角渦在垂向占據(jù)整個(gè)臺(tái)階高度；而在RANS計(jì)算中，角渦僅占約0.6個(gè)臺(tái)階高度.

圖6和圖7給出了區(qū)域內(nèi)不同位置的流向和垂向平均速度分布圖.圖8和圖9給出了區(qū)域內(nèi)不同位置的雷諾應(yīng)力及湍動(dòng)能分布.這些圖中包含對(duì)應(yīng)位置的DNS參考數(shù)據(jù).速度場(chǎng)的計(jì)算數(shù)據(jù)與參考數(shù)據(jù)吻合較好，湍流場(chǎng)曲線的峰值位置預(yù)測(cè)較準(zhǔn)確，但峰值略低于DNS的值.由圖6可見，由于臺(tái)階處流道面積突然擴(kuò)張，導(dǎo)致主流的流向速度在臺(tái)階處減小至0，并沿下游方向逐漸恢復(fù)，壁面附近出現(xiàn)回流，速度峰值約為主流的 20% ，且壁面附近速度逐漸減小，直到再附著點(diǎn)位置速度由負(fù)轉(zhuǎn)正.由圖9可見， x/h=1 處的湍動(dòng)能峰值與DNS數(shù)據(jù)相差較大，這可能是因RANS低估了角渦導(dǎo)致的.

4.2.2混合對(duì)流混合對(duì)流情況下計(jì)算得到的平均速度云圖如圖10所示.圖11給出了臺(tái)階后的流線分布圖，與強(qiáng)迫對(duì)流相比，這種情況下，流場(chǎng)的特點(diǎn)是由于浮力引起的流體加速在加熱面上形成了射流.此外，受浮力影響，回流區(qū)與主渦的縮小，相應(yīng)的角渦顯著增加，甚至超過主渦尺寸．從圖12可以看出，壁面摩擦系數(shù)均為正值，表明主渦已脫離壁面.

圖13和圖14給出了區(qū)域內(nèi)不同位置的流向和垂向平均速度分布圖.從圖中可以看出，浮力對(duì)于流場(chǎng)產(chǎn)生了明顯的影響.流向速度的峰值在下游部分從流動(dòng)中心位置轉(zhuǎn)向下壁面，同時(shí)為保證流量守恒，上壁面附近的流速逐漸減小，對(duì)于垂向速度分量，由于剪切層變小，導(dǎo)致流體需要更快地再附著到壁面，因此剪切層附近的垂向速度相較于強(qiáng)迫對(duì)流情況明顯變大，在恢復(fù)區(qū)內(nèi)浮力對(duì)于垂向速度分量的影響并不明顯．

受浮力的影響，最明顯的區(qū)別發(fā)生在下壁面．對(duì)于強(qiáng)迫對(duì)流，雷諾應(yīng)力的剪切應(yīng)力分量 ′v^′gt; 在整個(gè)下壁面附近是負(fù)數(shù)；而對(duì)于混合對(duì)流，在x/hgt;11 后發(fā)現(xiàn)剪切應(yīng)力在下壁面附近改變了符號(hào)，導(dǎo)致下壁面附近形成了新邊界層，通過剪切應(yīng)力也可以觀察到壁面射流與速度減小的外流之間的邊界.一般來(lái)說，在強(qiáng)迫對(duì)流下，湍流剪切應(yīng)力隨著距離臺(tái)階距離增加而減小，而在混合對(duì)流中，湍流剪切應(yīng)力則由于浮力的作用而增大.

圖15和圖16給出了區(qū)域內(nèi)不同位置的雷諾應(yīng)力及湍動(dòng)能分布．與平均速度場(chǎng)類似，湍流場(chǎng)也湍動(dòng)能表征湍流的強(qiáng)度，它受浮力影響導(dǎo)致的行為與雷諾應(yīng)力類似，在剪切層中迅速增長(zhǎng)并衰減，峰值向下壁面靠近；在恢復(fù)區(qū)內(nèi)，湍動(dòng)能逐漸減弱，但在浮力的影響下近壁面的湍流逐漸開始恢復(fù).在當(dāng)前研究的浮力條件！ （Ri=0.2） 下，臺(tái)階后的湍流強(qiáng)度整體上弱于強(qiáng)迫對(duì)流，從文獻(xiàn)[7]對(duì)液態(tài)鈉的混合對(duì)流 （Ri=0.338） 研究中可以觀察到，在靠近出口的位置湍流強(qiáng)度明顯大于強(qiáng)迫對(duì)流.

四方程模型預(yù)測(cè)了更低的動(dòng)量交換，因此再循環(huán)長(zhǎng)度略大于DNS結(jié)果，隨著離再循環(huán)區(qū)域距離的增加，結(jié)果與DNS結(jié)果之間的一致性提高.綜上所述，可以認(rèn)為四方程模型能夠較好地預(yù)測(cè)混合對(duì)流情況下的流場(chǎng).

4.3 溫度場(chǎng)

4.3.1強(qiáng)迫對(duì)流平均溫度的分布圖如圖17所示，在強(qiáng)迫對(duì)流下，壁面溫度最高的位置位于再循環(huán)區(qū)內(nèi)，從流道中心到下壁面方向溫度劇烈升高.沿下壁面的溫度分布和 Nu 數(shù)分布如圖18和圖19所示，對(duì)于強(qiáng)迫對(duì)流，由于被加熱流體的反向流動(dòng)減少了換熱，其速度相對(duì)較小，因此溫度在緊鄰臺(tái)階的位置達(dá)到最大值，并持續(xù)降低直到再附著點(diǎn)，在再附著點(diǎn)處溫度達(dá)到最低值.在隨后的恢復(fù)區(qū)由于壁面處恒定的熱流密度，壁面溫度呈線性增加.

四方程模型在 x/h=1 以外的位置都與DNS數(shù)據(jù)擬合得較好，壁面溫度曲線和 Nu 數(shù)曲線在回流區(qū)的拐點(diǎn)均比DNS數(shù)據(jù)更向上游.這可能是因?yàn)槟Ｐ偷凸懒私菧u的大小導(dǎo)致的．4.3.2混合對(duì)流混合對(duì)流的平均溫度分布如圖20所示.圖21和圖22給出了沿下壁面的平均溫度和 Nu 數(shù)分布圖.

由于液態(tài)LBE導(dǎo)熱性能很強(qiáng)，因此浮力引起的壁面速度場(chǎng)變化對(duì)于溫度場(chǎng)的影響相對(duì)較小.強(qiáng)迫對(duì)流和混合對(duì)流在恢復(fù)區(qū)內(nèi)的平均溫度 相差小于0.2.但是在回流區(qū)內(nèi)，由于浮力導(dǎo)致的主渦減小和角渦增大，壁面溫度的峰值向下游移動(dòng)，并且峰值大小顯著降低．觀察 Nu 數(shù)分布圖， Nu 數(shù)在回流區(qū)的最小值位置對(duì)應(yīng)于壁面平均溫度的最大值位置，相應(yīng)的最大值點(diǎn)對(duì)應(yīng)壁面平均溫度最小值的位置.浮力增強(qiáng)了角渦的對(duì)流效果，削弱了主渦向上游傳遞熱量，因此回流區(qū)的換熱整體得到提升.混合對(duì)流下主渦的消失，使得在下壁面的流向速度得到提前恢復(fù)，促進(jìn)了對(duì)流作用，因此恢復(fù)區(qū)的傳熱相較于強(qiáng)迫對(duì)流也得到了進(jìn)一步增強(qiáng)．

相比于強(qiáng)迫對(duì)流，混合對(duì)流情況下通過四方程模型得到的曲線與DNS數(shù)據(jù)擬合得更好.這可能是因?yàn)橹鳒u脫離壁面使得壁面上的渦結(jié)構(gòu)簡(jiǎn)單化，而在恢復(fù)區(qū)內(nèi)，無(wú)論是強(qiáng)迫對(duì)流還是混合對(duì)流都與DNS結(jié)果擬合得較好.因此，四方程模型可以較好地預(yù)測(cè)混合對(duì)流下的溫度場(chǎng).

5總結(jié)

液態(tài)金屬作為快中子反應(yīng)堆中優(yōu)秀的冷卻劑，其普朗特?cái)?shù) Pr 較低，溫度邊界層更厚，與傳統(tǒng)流體的湍流換熱表現(xiàn)相差很大，使用湍流數(shù)值計(jì)算中常用的湍流普朗特?cái)?shù)模型進(jìn)行計(jì)算會(huì)產(chǎn)生較大的誤差．本文針對(duì)垂直后臺(tái)階結(jié)構(gòu)，選取液態(tài)LBE作為工作流體，使用雷諾時(shí)均模擬中的四方程模型，針對(duì)強(qiáng)迫對(duì)流與混合對(duì)流兩種工況，分別進(jìn)行數(shù)值模擬研究并與直接數(shù)值模擬結(jié)果進(jìn)行對(duì)比，以研究四方程模型在復(fù)雜工況下計(jì)算液態(tài)LBE流動(dòng)傳熱的能力：

本文分別研究了 Ri=0 的強(qiáng)迫對(duì)流與 Ri=0.2 的混合對(duì)流情況下的數(shù)值結(jié)果與DNS結(jié)果之間的差異，對(duì)于強(qiáng)迫對(duì)流，四方程模型略微高估了主渦的大小；計(jì)算得到的平均流向速度與垂向速度與DNS結(jié)果擬合得較好；在雷諾剪切應(yīng)力與湍動(dòng)能方面預(yù)測(cè)得到的曲線峰值位置與DNS結(jié)果一致，但在臺(tái)階處的峰值大小略低；計(jì)算得到的壁面溫度峰值與努塞爾數(shù)峰值相較DNS結(jié)果略微向上游偏移.對(duì)于混合對(duì)流，四方程模型準(zhǔn)確地預(yù)測(cè)了主渦隨浮力增強(qiáng)與壁面分離的行為；平均流向速度和垂向速度同強(qiáng)迫對(duì)流情況類似，與DNS數(shù)據(jù)擬合得較好；湍動(dòng)能與雷諾剪切應(yīng)力曲線在恢復(fù)區(qū)內(nèi)與DNS相差很小，準(zhǔn)確預(yù)測(cè)了浮力增加所導(dǎo)致的流體域內(nèi)湍流強(qiáng)度下降以及在出口附近的射流；在溫度場(chǎng)方面也準(zhǔn)確預(yù)測(cè)了平均溫度峰值向下游的移動(dòng)．上述對(duì)比證明四方程模型可以較好地模擬液態(tài)LBE在復(fù)雜結(jié)構(gòu)下的混合對(duì)流傳熱：

本研究所采用的是各向同性四方程模型，所使用的標(biāo)量擴(kuò)散系數(shù)并不適合預(yù)測(cè)2種湍流熱流分量，導(dǎo)致湍流熱通量的流向分量被低估，后續(xù)可以通過改進(jìn)近似或考慮完整的各向異性效應(yīng)，以改進(jìn)現(xiàn)有的四方程模型.

參考文獻(xiàn)：

[1] ChengX，TakNI.Investigation on turbulent heat transferto lead-bismuth eutectic flows in circular tubes for nuclear applications[J].Nucl Eng Des， 2006，236（4）：385.

[2] Cheng X，Tak NI.CFD analysis of thermalhydraulic behavior of heavy liquid metals in subchannels[J].NuclEngDes，2006，236（18）：1874.

[3] Groetzbach G.Challenges in low-Prandtl number heat transfer simulation and modelling[J]. Nucl Eng Des，2013，264：41.

[4] HongB，ArmalyBF，Chen TS.Laminarmixed convection in a duct with a backward-facing step： Theeffects of inclination angle and Prandtl number[J].IntJHeatMassTran，1993，36（12）：3059.

[5] Moore T WF.Some experiments on the reattachmentofa laminar boundary layerseparating from a rearward facingstep ona flatplateaerofoil[J].JR Aeronaut Soc，1960，64（599）：668.

[6] NiemannM，F(xiàn)rohlich J.Directnumerical simulation of turbulent heat transfer behind a backward -facing step at low Prandtl number[J].Proc Appl Math Mech，2014，14（1）：659.

[7] NiemannM，F(xiàn)rohlich J.Buoyancy-affected backwardfacing step flow with heat transfer at low Prandtl number[J].IntJHeatMassTran，2O16，1Ol：1237.

[8] Niemann M，F(xiàn)rohlich J.Buoyancy effects on turbulent heat transfer behind a backward-facing step in liquid metal flow[M]//GrigoriadisD，GeurtsB，KuertenH，etal.Directandlarge-eddy simulation X. ERCOFTAC Series，vol 24.Cham：Springer， 2017：513.

[9] DeSantis A，Shams A.Application ofan algebraic turbulent heat flux model to a backward facing step flow at low Prandtl number [J].Ann Nucl Energy， 2018，117：32.

[10]Da Via R，Manservisi S.Numerical simulation of forced and mixed convection turbulent liquid sodium flow over a vertical backward facing step with a four parameter turbulence model [J]. Int JHeat Mass Tran，2019，135：591.

[11] GeZH.Direct numerical simulation study of turbulentheattransferin liquidmetals[D].Hefei：Univer sity of Science and technology of China，2Ol8.[葛志 浩.液態(tài)金屬湍流換熱的直接數(shù)值模擬研究[D]. 合肥：中國(guó)科學(xué)技術(shù)大學(xué)，2018.]

[12]Zhao P，WangC，Ge Z，et al.DNS of turbulent mixed convection over a vertical backward-facing stepfor lead-bismuth eutectic[J].IntJHeatMass Tran，2018，127：1215.

[13]Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications [J].AIAA J， 1994，32（8）：1598.

[14] Van Maele K，MerciB.Application of two buoyancymodified k-εturbulencemodelsto differenttypesof buoyant plumes[J].Fire SafetyJ，2006，41（2）：122.

[15] SuXK.Numerical study of liquid lead-bismuth turbulent heat transfer based on isotropic four-equation modeling[D].Beijing：University ofChineseAcademyof Sciences（Institute of modern physics，ChineseAcademy of Sciences），2022.[蘇興康.基于各 向同性四方程模型的液態(tài)鉛鉍湍流換熱數(shù)值研 究[D].北京：中國(guó)科學(xué)院大學(xué)（中國(guó)科學(xué)院近代物 理研究所），2022.]

[16]KasagiN，Matsunaga A.Three-dimensional particletrackingvelocimetrymeasurementof turbulence statisticsand energy budget in a backward-facing step flow[J].IntJHeatFluidFl，1995，16（6）：477.

（責(zé)任編輯：于白茹）