基于OpenFOAM的液態(tài)金屬鉛鉍三維流動換熱特性數(shù)值模擬研究

何少鵬，王明軍,*，章 靜，田文喜，蘇光輝，秋穗正

(1.西安交通大學(xué) 動力工程多相流國家重點實驗室，陜西 西安 710049；2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

液態(tài)金屬鉛鉍快堆作為第4代核能系統(tǒng)的主要堆型之一，有自然循環(huán)能力強、中子經(jīng)濟(jì)性好等諸多優(yōu)點。但液態(tài)金屬鉛鉍由于特殊的熱物性質(zhì)(如低普朗特數(shù)、高熱導(dǎo)率等)，其流動邊界層與熱邊界層已不在一個量級，流動換熱特性已與水、空氣等常見工質(zhì)有很大不同[1]，這給液態(tài)金屬鉛鉍反應(yīng)堆的研究設(shè)計帶來了許多困難。故亟需針對液態(tài)金屬鉛鉍開發(fā)更為精確的湍流換熱模型以進(jìn)行高精度數(shù)值模擬。

在主流商用CFD軟件中，通常采取定值湍流普朗特數(shù)的SED模型對液態(tài)金屬流動換熱問題進(jìn)行模擬求解。該模型基于雷諾比擬原理，即假設(shè)熱邊界層與流動邊界層具有相似性，其成立的條件為Pr接近于1。對于空氣、水等常見工質(zhì)，使用湍流普朗特數(shù)為0.85～1的SED模型往往能得到工程上可接受的結(jié)果。但針對液態(tài)金屬等低普朗特數(shù)流體開展的一些CFD模擬研究表明，定值湍流普朗特數(shù)的SED模型已不能達(dá)到很好的預(yù)測精度。Kawamura等[2]利用直接數(shù)值模擬研究了雷諾數(shù)與普朗特數(shù)對液態(tài)金屬流動換熱的影響，結(jié)果表明采用定值湍流普朗特數(shù)的SED模型已不能很好地與實驗關(guān)聯(lián)式吻合。Cheng等[3]通過CFD模擬研究了單通道內(nèi)液態(tài)金屬的熱工水力現(xiàn)象，發(fā)現(xiàn)當(dāng)SED模型的湍流普朗特數(shù)由0.9增加到1.5時，預(yù)測努塞爾數(shù)將會下降約30%。Chen等[4]使用一種新網(wǎng)格標(biāo)記分離方法研究了繞絲棒束組件三維流動換熱特性，結(jié)果表明繞絲會引起橫流以及軸向速度周期振蕩，影響熱工水力特性。Manservisi等[5]針對液態(tài)金屬開發(fā)了一種更為精確的k-ε-kθ-εθ四因子模型。該模型是一種顯式的代數(shù)熱通量模型，充分考慮了熱邊界層與流動邊界層的差異性，引入脈動溫度方均值及其耗散率兩個新參數(shù)，建立了湍流熱擴散率與脈動溫度方均值間的數(shù)值關(guān)系，能實現(xiàn)對液態(tài)金屬湍流換熱的精確模擬。Manservisi等[6-7]使用四因子模型在三角形柵元通道和四邊形柵元通道內(nèi)進(jìn)行了數(shù)值模擬與驗證，均取得了較好的預(yù)測結(jié)果。

考慮到主流商業(yè)軟件代碼閉源，封裝嚴(yán)格，不利于用戶對模型進(jìn)行修改和開發(fā)，而開源計算流體力學(xué)軟件OpenFOAM采用C++語言編程，具有編程環(huán)境自由、模型開發(fā)簡便等優(yōu)點，廣泛應(yīng)用于能源化工、海洋船舶等領(lǐng)域?；陂_源平臺進(jìn)行模型添加和求解器開發(fā)對實現(xiàn)更高精度的液態(tài)金屬流動換熱CFD模擬有重要意義。故本文基于k-ε-kθ-εθ四因子模型在開源程序OpenFOAM中開發(fā)針對液態(tài)金屬的自定義求解器，闡述數(shù)值模擬中所使用的數(shù)學(xué)物理模型，包括守恒方程與四因子模型方程，并對帶繞絲棒束通道中液態(tài)金屬鉛鉍的流動換熱問題進(jìn)行模擬研究。

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

1.1 守恒方程

液態(tài)金屬反應(yīng)堆中液態(tài)金屬冷卻劑在正常工況下的流動可視為單相、不可壓縮的湍流流動。因此，經(jīng)過不可壓縮流簡化及雷諾時均處理后，動量、質(zhì)量、能量守恒方程如式(1)所示。

(1)

式中：U為速度平均值；T為溫度平均值；U′為湍流脈動速度；T′為湍流脈動溫度；cp為比定壓熱容；λ為熱導(dǎo)率。

(2)

Prt=νt/αt

(3)

湍流普朗特數(shù)給定后(本文分別取針對液態(tài)金屬推薦的湍流普朗特數(shù)1.5和2.0進(jìn)行模擬計算)，借由湍流模型求得的局部湍流黏度νt，便可進(jìn)一步求解出湍流熱擴散系數(shù)αt。至此，方程組(1)得以封閉，方程組理論可解。

1.2 k-ε-kθ-εθ四因子模型

1) 低雷諾數(shù)修正的k-ε模型

標(biāo)準(zhǔn)k-ε模型是一種較為常見的二方程湍流模型，適合模擬充分發(fā)展的湍流流動，對低雷諾數(shù)過渡流和近壁區(qū)域計算結(jié)果不夠理想。在四因子模型中，應(yīng)用了針對低雷諾數(shù)修正的k-ε方程。與標(biāo)準(zhǔn)k-ε模型相似，引入了Boussinesq假設(shè)，認(rèn)為湍流雷諾應(yīng)力與牛頓剪切引力相似，均與應(yīng)變呈正比，雷諾應(yīng)力張量可寫為湍流黏度與速度梯度的乘積，即：

(4)

式中，τt為湍流雷諾應(yīng)力。

與標(biāo)準(zhǔn)k-ε方程相似，低雷諾數(shù)修正的k-ε方程也定義了2個新的變量，即湍動能k及其耗散率ε，如式(5)所示。

(5)

與標(biāo)準(zhǔn)k-ε模型不同的是，本模型定義的湍流黏度為：

νt=Cμkτlu

(6)

式中，Cμ為模型常數(shù)。為更加準(zhǔn)確地模擬流體的近壁面行為，該模型定義了局部動力學(xué)特征時間τlu[9-10]，由4個參數(shù)計算得到，即：

τlu=f1μA1μ+f2μA2μ

(7)

式(7)中4個參數(shù)均模化為歸一化壁面距離Rδ與湍流雷諾數(shù)Ret的函數(shù)，如式(8)所示。

(8)

歸一化壁面距離Rδ和湍流雷諾數(shù)Ret定義如下：

(9)

式中：δ為各網(wǎng)格單元到壁面的距離；ν為層流運動黏性系數(shù)。

為計算局部湍流黏度νt，需得到式(7)中定義的局部動力學(xué)特征時間τlu，即需要解出歸一化壁面距離和湍流雷諾數(shù)，進(jìn)而需解出湍動能k及其耗散率ε。修正后的k-ε方程如式(10)所示。

(10)

式中：

(11)

(12)

關(guān)于低雷諾數(shù)修正的k-ε模型的更多闡釋可參考文獻(xiàn)[11]。

2)kθ-εθ模型

與k-ε模型類似，kθ-εθ模型同樣定義了2個新變量：脈動溫度方均值kθ及其耗散率εθ。

(13)

與SED模型不同的是，在四因子模型中將湍流熱擴散率αt定義為式(14)所示的3個參數(shù)的乘積。

αt=Cθkτlθ

(14)

式中：Cθ為模型常數(shù)；τlθ為局部熱力學(xué)特征時間，用來更加精確地描述湍流流動中的換熱特性[5]。

τlθ=f1θB1θ+f2θB2θ

(15)

式中，f1θ、f2θ、B1θ與B2θ4個參數(shù)被定義為歸一化壁面距離Rδ、湍流雷諾數(shù)Ret等變量的函數(shù)。

(16)

(17)

式中，R為熱力學(xué)特征時間與動力學(xué)特征時間的比值。

(18)

綜上所述，為計算湍流熱擴散系數(shù)αt，需得到局部熱力學(xué)特征時間τlθ，進(jìn)一步需解出脈動溫度方均值kθ及其耗散率εθ，總結(jié)kθ-εθ方程組如下：

(19)

式中：Cd1=1、CP1=0.925、CP2=0.9均為模型常數(shù)；Pk定義見式(11)；Pθ與Cd2定義見式(20)。

(20)

2 帶繞絲棒束通道CFD模擬

2.1 幾何模型

本文選取德國卡爾斯魯厄理工學(xué)院(KIT)盒間流實驗[13]中的帶繞絲棒束組件作為CFD模擬和模型驗證對象。該棒束組件為六邊形，包含有7根帶繞絲燃料棒，每根棒表面施以均勻熱流密度以模擬燃料棒對鉛鉍合金流體加熱。從組件進(jìn)口向出口看，7根繞絲順時針旋轉(zhuǎn)并纏繞在7根棒上。該帶繞絲棒束組件具體幾何模型與截面視圖如圖1所示，主要幾何尺寸如下：P=20.5 mm、D=16.0 mm、d=4.4 mm、W=20.8 mm、FF=61.05 mm、L=262 mm、Dh=10.34 mm。

圖1 帶繞絲棒束幾何模型及其截面視圖Fig.1 Geometric model and sectional view of wire-wrapped rod bundle

2.2 計算方法

本文選取繞絲螺距長度為262 mm的1個通道作為物理模型。進(jìn)口溫度均為573 K，施加速度入口邊界條件。出口為固定壓力條件，組件壁和繞絲表面設(shè)為絕熱壁面，棒束表面施加均勻的熱流密度。模擬中考慮了液態(tài)鉛鉍物性隨溫度的變化，梯度項離散采用Gauss linear格式，擴散項采用Bounded Gauss linear格式，其余使用一階迎風(fēng)差分格式，分別使用SED模型與四因子模型進(jìn)行計算，具體模擬工況列于表1。

表1 模擬工況Table 1 Simulated condition

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

如圖1所示，本文選取組件出口截面上1條垂直線段A，以沿線段A的冷卻劑溫度變化特性及速度分布為對象獲得網(wǎng)格獨立性解，如圖2所示。從計算結(jié)果可看出，當(dāng)網(wǎng)格數(shù)為166萬、196萬和234萬時，沿線段A的冷卻劑速度、溫度變化相差不大，故最終選取166萬網(wǎng)格數(shù)的網(wǎng)格為網(wǎng)格獨立解開展后續(xù)數(shù)值模擬研究。

圖2 網(wǎng)格無關(guān)性分析Fig.2 Grid independence analysis

2.4 計算結(jié)果分析

帶繞絲棒束通道內(nèi)各截面處的速度、溫度云圖如圖3所示。結(jié)果表明，由于繞絲對流體的攪混作用，各截面處速度分布呈現(xiàn)非對稱、不均勻特征，且沿軸向發(fā)生劇烈變化。由溫度云圖還可看出，通道內(nèi)不對稱速度分布引起了不對稱溫度分布。通道中心截面橫流速度場矢量圖如圖4所示，由于繞絲的攪混作用，橫流速度場呈現(xiàn)出旋轉(zhuǎn)特點，且橫流旋轉(zhuǎn)方向與繞絲旋轉(zhuǎn)方向一致。

圖3 各截面速度與溫度云圖Fig.3 Velocity and temperature contour of several sections

圖4 中間截面(z=131 mm)橫流速度矢量圖Fig.4 Cross flow velocity vector graph of middle section (z=131 mm)

本文共選取3個典型局部通道(圖5)，即中心通道、邊通道、角通道進(jìn)行模型驗證。選取截面位置為z=262 mm，即出口截面。針對液態(tài)金屬建立的Mikityuk換熱實驗關(guān)聯(lián)式[14]和Ushakov換熱實驗關(guān)聯(lián)式[15]如式(21)、(22)所示，式中X為柵距棒徑比，X=P/D。

圖5 局部通道示意圖Fig.5 Local channel schematic diagram

1) Mikityuk實驗關(guān)聯(lián)式：

Nu=0.047(1-e-3.8(X-1))(Pe0.77+250)

30≤Pe≤5 000;1.1≤X≤1.95

(21)

2) Ushakov實驗關(guān)聯(lián)式：

1≤Pe≤4 000;1.3≤X≤2.0

(22)

關(guān)聯(lián)式(21)、(22)的實驗條件與本文模擬條件之間的對比列于表2。由表2可知，本文模擬條件與關(guān)聯(lián)式在棒束柵格幾何形狀、柵距棒徑比X、Pe等方面具有較好的一致性。盡管本文模擬幾何的X未在Ushakov關(guān)聯(lián)式推薦范圍內(nèi)，但文獻(xiàn)[14]將X為1.15～1.95的實驗數(shù)據(jù)與Ushakov關(guān)聯(lián)式進(jìn)行了對比，結(jié)果表明，對于X<1.3的實驗數(shù)據(jù)，Ushakov關(guān)系式也能符合很好。因此在本文模擬條件下，Mikityuk關(guān)聯(lián)式與Ushakov關(guān)聯(lián)式具有一定的可靠度。

表2 關(guān)聯(lián)式實驗條件與模擬條件的對比Table 2 Comparison of correlations’ experimental condition and simulation condition

使用湍流普朗特數(shù)取1.5和2.0的SED模型以及k-ε-kθ-εθ四因子模型的計算結(jié)果示于圖6，圖中M關(guān)聯(lián)式為Mikityuk關(guān)聯(lián)式預(yù)測值，U關(guān)聯(lián)式為Ushakov關(guān)聯(lián)式預(yù)測值。對比結(jié)果表明，四因子模型與實驗關(guān)聯(lián)式符合良好，在本文所模擬條件下能更為精確地預(yù)測繞絲棒束通道中的Nu。

圖6 中心通道、邊通道、角通道內(nèi)努塞爾數(shù)Fig.6 Nusselt numbers of center, edge and corner channels

圖7為使用SED模型(Prt=1.5)和四因子模型得到的沿圖1中直線A的湍流黏度νt與湍流熱擴散率αt的分布曲線。

a——SED模型(Prt=1.5)；b——四因子模型圖7 沿線段A的湍流黏度和湍流熱擴散率Fig.7 Turbulent viscosity and turbulent thermal diffusivity along line A

從圖7a可看出，使用SED模型時，νt與αt的變化曲線具有高度一致性，兩者呈現(xiàn)一定的等比例關(guān)系。這是由于SED模型基于雷諾比擬原理，直接使用式(3)進(jìn)行湍流熱擴散率αt的計算，且由式(3)可知，νt與αt間的比例系數(shù)即為定值的湍流普朗特數(shù)1.5。然而對于低普朗特數(shù)的液態(tài)金屬鉛鉍，其動量擴散與能量擴散具有不同的行為特征[16]，此時使用假設(shè)熱邊界層與流動邊界層相似的雷諾比擬原理往往會造成較大的誤差。

與SED模型不同，四因子模型考慮了熱邊界層與流動邊界層的差異，引入了脈動溫度方均值及其耗散率2個參數(shù)來描述熱邊界層的能量耗散行為，構(gòu)造控制方程(式(19))求解湍流熱擴散率αt，從而避免了使用雷諾比擬原理。從圖7b可看出，νt和αt間并不呈完全的等比例關(guān)系，尤其在棒束壁面附近，表征動量擴散的νt與表征能量擴散的αt呈不同的變化趨勢，說明四因子模型能較好地模擬動量擴散與能量擴散的不一致性，更加準(zhǔn)確地描述液態(tài)金屬鉛鉍流場中的熱擴散行為。

上述結(jié)果表明：四因子模型在預(yù)測繞絲棒束通道內(nèi)的Nu、描述液態(tài)鉛鉍流場動量和能量擴散行為方面相比傳統(tǒng)SED模型表現(xiàn)出了一定優(yōu)勢，本文基于四因子模型開發(fā)的適用于低普朗特數(shù)液態(tài)金屬流動傳熱模擬的OpenFOAM求解器可較好地預(yù)測該模擬工況下繞絲棒束通道中液態(tài)金屬鉛鉍流動換熱現(xiàn)象。

3 結(jié)論

本文基于開源CFD平臺OpenFOAM開發(fā)了適用于液態(tài)金屬模擬的k-ε-kθ-εθ四因子模型自定義求解器，對帶繞絲棒束通道中液態(tài)鉛鉍流動換熱現(xiàn)象進(jìn)行了模擬研究。計算中充分考慮了熱邊界層與流動邊界層的差異性，采用針對低雷諾數(shù)修正的k-ε模型求解湍流黏度，使用1組kθ-εθ偏微分方程求解湍流熱擴散率，實現(xiàn)了對液態(tài)金屬流動換熱問題較為精確的數(shù)值模擬。計算得到了帶繞絲棒束通道中液態(tài)金屬鉛鉍速度、溫度等重要熱工水力參數(shù)三維分布，深入分析了繞絲對棒束通道內(nèi)液態(tài)金屬流動換熱過程的影響規(guī)律，對比分析了四因子模型和傳統(tǒng)SED模型的模擬結(jié)果，證明了四因子模型在預(yù)測液態(tài)金屬鉛鉍流動傳熱過程中具有一定優(yōu)勢，上述模型和求解器能較好地模擬液態(tài)金屬鉛鉍在繞絲棒束中的流動換熱問題。