湍流模型對瞬態(tài)流固耦合換熱數(shù)值模擬影響分析

丁文杰，郭海兵，黃洪文，郭斯茂

(中國工程物理研究院 核物理與化學(xué)研究所，四川 綿陽 621900)

聚變-裂變混合能源堆(Z-FFR)利用電磁內(nèi)爆加熱并壓縮氘氚靶丸發(fā)生聚變，產(chǎn)生聚變中子驅(qū)動次臨界包層實現(xiàn)能量放大。由于聚變中子源的脈沖特性，98.7%的能量在0.01 s 內(nèi)沉積于次臨界包層的燃料中[1]，燃料瞬時溫度階躍最大可達(dá)300 K[2]。Z-FFR的瞬態(tài)釋能特性與傳統(tǒng)穩(wěn)態(tài)釋能系統(tǒng)存在較大差異，給次臨界包層的熱工水力設(shè)計帶來未知和不確定因素。已有部分學(xué)者對功率瞬態(tài)變化條件下的流動傳熱問題展開研究。1960年，Sparrow等[3]通過對湍流模型求解近似解獲得了圓管湍流流動的瞬態(tài)換熱響應(yīng)結(jié)果。1967年，Soliman等[4]通過考慮湍流邊界層的效應(yīng)獲得了平板溫度變化的解析解，但該結(jié)果的換熱系數(shù)高出實驗值的50%。1983年，Steward等[5]開展了液氦在受到壁面階躍熱流量時瞬態(tài)強制對流換熱響應(yīng)的實驗研究。2005年，Bhowmik等[6]對電加熱片的瞬態(tài)對流換熱進(jìn)行了實驗研究，并根據(jù)實驗數(shù)據(jù)擬合了計算瞬態(tài)換熱系數(shù)的經(jīng)驗關(guān)系式。Liu等[7-9]對氦氣及氮氣外掠加熱平板、水平圓柱時功率指數(shù)上升的瞬態(tài)強制對流換熱進(jìn)行了理論與實驗研究，并根據(jù)理論與實驗研究的結(jié)果擬合了計算瞬態(tài)換熱系數(shù)的經(jīng)驗關(guān)系式。

以上研究通過理論或?qū)嶒灥姆椒ǐ@得了瞬態(tài)對流傳熱過程中的響應(yīng)，但這些文獻(xiàn)均未給出不同湍流模型對數(shù)值模擬結(jié)果的影響。本文基于ANSYS Fluent程序和Liu等[7-9]的實驗數(shù)據(jù)，對瞬態(tài)對流換熱過程中的湍流模型進(jìn)行敏感性分析。

1 數(shù)值模型

1.1 幾何建模

熱負(fù)載功率指數(shù)提升實驗平臺由日本神戶大學(xué)Liu等搭建，用于分析高溫氣冷堆功率瞬變過程中的對流換熱影響。該平臺實驗段為一水平放置、內(nèi)徑20 mm的圓形流道，流道內(nèi)的流體可采用氦氣或氮氣，在流道中心位置水平放置加熱元件，該加熱元件為薄平板狀或細(xì)圓柱狀。以薄平板狀加熱元件為例(圖1)，加熱元件為一厚0.1 mm、寬4.0 mm、長50 mm的金屬鉑片，為防止實驗過程中加熱元件發(fā)生移動，金屬鉑片兩端均采用長為45 mm的銅片固定。

圖1 實驗段示意圖Fig.1 Schematic of test section

根據(jù)上述描述，采用GAMBIT軟件對該實驗段進(jìn)行建模。當(dāng)加熱元件為細(xì)圓柱狀時，可直接根據(jù)軸對稱原則建立二維模型。當(dāng)加熱元件為薄平板狀時，固體過薄導(dǎo)致網(wǎng)格劃分量大，不利于瞬態(tài)計算，為此，分別建立三維模型與簡化二維模型計算相同的穩(wěn)態(tài)工況，兩種模型計算的壁面平均換熱系數(shù)相對偏差為2.1%，即加熱元件為薄平板狀時同樣可采用二維模型。薄平板狀加熱元件的幾何模型如圖2所示。

圖2 薄平板狀加熱元件的幾何模型Fig.2 Geometric model of thin plate heating element

1.2 物性參數(shù)選取

1) 流體熱物性

實驗段出口壓力為0.5 MPa，且實驗過程中保持不變，故流體熱物性參數(shù)可看作是僅隨溫度變化的函數(shù)。通過文獻(xiàn)[10-14]可獲得氮氣與氦氣在0.5 MPa、300～1 000 K下的熱物性參數(shù)，進(jìn)而采用多項式擬合的方法可獲得兩種氣體熱物性參數(shù)的擬合關(guān)系式，供ANSYS Fluent程序分析使用。

氮氣熱物性擬合關(guān)系式如下：

ρN2=-1.479 8e-8T3+3.766 56e-5T2-

0.034 01T+12.801 67

(1)

cp,N2=-3.535 35e-10T3+ 8.679 65e-7T2-

4.642 14e-4T+ 1.120 24

(2)

λN2=-7.857 14e-9T2+

6.724e-5T+6.577 38e-3

(3)

ηN2=-1.333 33e-5T2+

0.050 76T+4.063 1

(4)

式中：ρ為密度，kg/m3；cp為比定壓熱容，kJ·kg-1·K-1；λ為熱導(dǎo)率，W·m-1·K-1；η為動力黏性系數(shù)，μPa·s；T為溫度，K。

氦氣的比定壓熱容取常數(shù)5.19 kJ·kg-1·K-1，其余熱物性參數(shù)擬合關(guān)系式如下：

ρHe=-2.095 96e-9T3+5.343 07e-6T2-

0.004 83T+1.820 48

(5)

λHe=-7.202 38e-8T2+

3.842 3e-4T+0.048 09

(6)

ηHe=-8.333 33e-6T2+

0.048 19T+6.284 52

(7)

2) 固體熱物性

在計算過程中，假定固體材料的熱物性不隨溫度發(fā)生變化。根據(jù)文獻(xiàn)[15]可獲得金屬鉑與銅的熱物性參數(shù)，如表1所列。

表1 固體材料物性參數(shù)Table 1 Property parameter of solid material

1.3 湍流模型與網(wǎng)格無關(guān)性分析

ANSYS Fluent程序提供了大量的湍流模型，如k-ε模型、k-ω模型等，其中每一大類的湍流模型又包含若干個子類模型。考慮到本文所研究的對象為簡單的壁面流動，并不考慮強旋流動和彎曲壁面的影響，故選取4種典型的湍流模型進(jìn)行分析，包括標(biāo)準(zhǔn)k-ε模型、標(biāo)準(zhǔn)k-ω模型、過渡SST模型以及LRRIP雷諾應(yīng)力模型(RSM模型)。

壁面函數(shù)法是基于所計算問題的壁面附近黏性支層以外區(qū)域的無量綱速度與溫度分別服從對數(shù)律分布的假設(shè)而獲得的[16]。在功率指數(shù)瞬變過程中，靠近加熱壁面附近的流體密度將隨溫度瞬變而改變，引起近壁處的流場發(fā)生變化，破壞了無量綱速度與溫度的對數(shù)律分布假定。為此，通過控制邊界層網(wǎng)格y+的大小，確保所有湍流模型的壁面模型均采用強化壁面處理而非壁面函數(shù)法。建立5套網(wǎng)格劃分方案，采用不同湍流模型進(jìn)行穩(wěn)態(tài)模擬，結(jié)果示于圖3，可見，當(dāng)y+在1左右時，各湍流模型模擬的穩(wěn)態(tài)換熱系數(shù)幾乎保持不變，即具有網(wǎng)格無關(guān)解。

圖3 不同網(wǎng)格劃分方案下的穩(wěn)態(tài)換熱系數(shù)Fig.3 Steady-state heat transfer coefficient in different grids

2 計算結(jié)果與分析

2.1 瞬態(tài)過程中的流固耦合換熱特性

在功率指數(shù)增長實驗過程中，加熱元件被直流電源加熱，采用控制系統(tǒng)對熱功率進(jìn)行控制并測量其兩端的電壓與電流，進(jìn)而根據(jù)計算的電阻值反推加熱元件的體平均溫度。加熱元件的體功率瞬態(tài)變化如下：

Q(t)=Q0exp(t/τ)

(8)

式中：t為瞬態(tài)時間；τ為指數(shù)增長周期；Q(t)為t時刻加熱元件的體功率密度；Q0為初始時刻的體功率密度。

加熱元件的等效表面熱流密度采用下式計算：

(9)

式中：q(t)為t時刻加熱元件的等效表面熱流密度；δ為平板加熱元件的厚度(當(dāng)加熱元件為圓柱狀時為該圓柱的半徑)；ρs、cs、Ts分別為加熱元件的密度、比熱容和體平均溫度。

通過式(9)獲得加熱元件的等效表面熱流密度，該等效熱流密度的最大不確定度為2.4%[10]。

采用k-ε模型模擬氦氣外掠平板流動的功率、溫升瞬態(tài)變化，結(jié)果示于圖4。其模擬工況為：流體速度10 m/s、氣體入口溫度313 K、功率指數(shù)增長周期分別為92 ms和735 ms。加熱元件的溫升與主流流體的溫升分別代表加熱元件溫度、主流流體溫度與氣體入口溫度的差值。從圖4可看出，在加熱元件的體功率指數(shù)增長過程中，其等效表面熱流密度、體平均溫度以及流體溫度也隨之呈指數(shù)形式增長。根據(jù)式(9)可計算出：當(dāng)指數(shù)增長周期為92 ms時，加載到金屬鉑上的熱量約有20%將會被冷卻流體以對流的形式帶走，而剩下的近80%的熱量仍留在固體內(nèi)，使固體域的溫度升高；當(dāng)指數(shù)增長周期為735 ms時，約有60%的熱量被流體帶走，而僅有約40%的熱量留在固體內(nèi)。不難發(fā)現(xiàn)，瞬態(tài)加熱過程與穩(wěn)態(tài)時熱量全部傳遞到流體的過程存在明顯的不同：在加熱功率指數(shù)增長的過程中，加熱功率一部分被流體以對流形式帶走，而另一部分繼續(xù)留在固體內(nèi)，使固體升溫。流體對流換熱帶走的熱量占加熱功率的比值與加熱功率的提升速率有關(guān)，即提升速率越快，對流換熱帶走的熱量占比越低。

當(dāng)模擬工況中采用氦氣工質(zhì)，且流體速度為10 m/s、氣體入口溫度為313 K、功率指數(shù)增長周期為92 ms時，不同瞬態(tài)時刻沿金屬鉑中心位置法線方向的流體速度分布與溫度分布示于圖5(y為網(wǎng)格中心點距金屬鉑中線位置的距離)。從圖5a可看出，流體速度沿法線方向的分布隨時間而變化，這與采用定熱物性假設(shè)而獲得的速度分布不變有明顯的差異。此外，隨著時間的增加，靠近加熱壁面流體的流速不斷減小，而遠(yuǎn)離加熱壁面流體的流速不斷增加，這是由于壁面溫度不斷升高進(jìn)而使靠近壁面的流體熱膨脹，且流體黏性增大所導(dǎo)致的。從圖5b可看出，在靠近加熱壁面1 mm處的流體溫度沿法向方向劇烈變化，而遠(yuǎn)離加熱壁面的流體溫度仍保持入口溫度。

圖4 不同指數(shù)增長周期下的傳熱特性Fig.4 Heat transfer characteristic in different periods

圖5 不同時刻的流體速度與溫度分布Fig.5 Velocity and temperature distribution at different time

2.2 湍流模型對瞬態(tài)流固耦合換熱的影響

1) 氮氣外掠平板流動實驗

使用不同湍流模型對氮氣外掠平板流動實驗(實驗工況為：流體速度4 m/s、氣體入口溫度313 K、功率指數(shù)增長周期174 ms)進(jìn)行模擬，將模擬的固體域溫度以及根據(jù)式(9)計算的等效表面熱流密度與實驗結(jié)果進(jìn)行對比，以分析不同湍流模型對流固耦合換熱的影響，結(jié)果示于圖6。從圖6a可看出，采用SST模型模擬的加熱元件體溫度略高于其他3種模型的模擬結(jié)果，且4種模型的模擬結(jié)果整體上均與實驗數(shù)據(jù)相吻合。從圖6b可看出，不同湍流模型模擬的等效表面熱流密度變化趨勢一致，但4種模型的模擬結(jié)果均低于實驗值，其中k-ε模型與實驗數(shù)據(jù)最接近，而SST模型與實驗數(shù)據(jù)相差最大。

圖6 氮氣外掠平板流動的模擬結(jié)果與實驗數(shù)據(jù)的對比Fig.6 Comparison between numerical results and experimental data for nitrogen gas flowing over plate

為整體考慮湍流模型的性能，需進(jìn)一步比較各模型的均方根偏差E。定義E為：

(10)

式中：qCAL為根據(jù)數(shù)值模擬結(jié)果獲得的等效表面熱流密度；qEXP為根據(jù)實驗數(shù)據(jù)獲得的等效表面熱流密度；n為實驗數(shù)據(jù)點個數(shù)。

根據(jù)式(10)計算各模型的均方根偏差，結(jié)果示于圖7。從圖7可看出，SST模型的均方根偏差最大，接近28.3%。本次模擬結(jié)果中，均方根偏差最小的為k-ε模型，RSM模型次之，SST模型最大。

2) 氦氣外掠平板實驗

使用4種湍流模型對氦氣外掠平板流動實驗(實驗工況為：流體速度10 m/s、氣體入口溫度313 K、功率指數(shù)增長周期92 ms)進(jìn)行模擬，并與實驗結(jié)果進(jìn)行對比，結(jié)果示于圖8。從圖8可看出，4種湍流模型計算的等效表面熱流密度均低于實驗值，其中k-ε模型與RSM模型與實驗值較為接近，而k-ω模型及SST模型與實驗值偏差較大。

圖7 各模型模擬氮氣外掠平板流動的均方根偏差Fig.7 Root mean square deviation of different turbulence models for nitrogen gas flowing over plate

各模型的均方根偏差示于圖9。從圖9可看出，采用4種湍流模型獲得的等效表面熱流密度均與實驗值有較大偏差，均方根偏差均超過20%。但相比于其他3種湍流模型，k-ε模型仍是均方根偏差最小的模型。

圖8 氦氣外掠平板流動的模擬結(jié)果與實驗值的對比Fig.8 Comparison between numerical results and experimental data for helium gas flowing over plate

圖9 各模型模擬氦氣外掠平板流動的均方根偏差Fig.9 Root mean square deviation of different turbulence models for helium gas flowing over plate

3) 氦氣外掠圓柱實驗

使用4種湍流模型對氦氣外掠圓柱流動實驗(實驗工況為：流體速度10 m/s、氣體入口溫度313 K、功率指數(shù)增長周期1.5 s)進(jìn)行模擬，并與實驗結(jié)果進(jìn)行對比，結(jié)果示于圖10。從圖10可看出，4種湍流模型計算的等效表面熱流密度相互間差異較小，且均與實驗值相吻合。

各模型的均方根偏差示于圖11。從圖11可看出，采用4種湍流模型獲得的等效表面熱流密度與實驗值的均方根偏差均低于16%。其中，k-ε模型與RSM模型的均方根偏差較小，SST模型次之，k-ω模型最大。

4) 數(shù)值模擬與實驗值的對比分析

本研究選取的4種湍流模型，根據(jù)求解變量的差異又可分為ε類湍流模型與ω類湍流模型。其中，標(biāo)準(zhǔn)k-ε模型與RSM模型屬于ε類湍流模型，而標(biāo)準(zhǔn)k-ω模型與過渡SST模型則屬于ω類湍流模型。

圖10 氦氣外掠圓柱流動的模擬結(jié)果與實驗數(shù)據(jù)對比Fig.10 Comparison between numerical results and experimental data for helium gas flowing over cylinder

圖11 各模型模擬氦氣外掠圓柱流動的均方根偏差Fig.11 Root mean square deviation of different turbulence models for helium gas flowing over cylinder

在ANSYS Fluent程序求解湍流流動的能量守恒方程時，一般認(rèn)為速度場和溫度場可比擬，并通過引入湍流普朗特數(shù)Prt計算湍流導(dǎo)熱系數(shù)，即：

(11)

式中：μt為湍流黏性系數(shù)；λt為湍流導(dǎo)熱系數(shù)。

研究過程中Prt取ANSYS Fluent程序的默認(rèn)值0.85。根據(jù)式(11)可知，湍流模型計算的湍流黏性系數(shù)越高，湍流導(dǎo)熱系數(shù)越大，則湍流換熱越強，等效表面熱流密度越大。故而不同湍流模型計算的等效表面熱流密度差異直接取決于湍流模型計算的湍流黏性系數(shù)。

3組實驗數(shù)據(jù)與數(shù)值模擬結(jié)果的對比結(jié)果(圖6b、8、10)表明：采用不同湍流模型計算的等效表面熱流密度均低于實驗值，這可能是由于瞬態(tài)加熱過程中壁面附近流體發(fā)生擾動強化了湍流換熱，而該強化作用在數(shù)值模擬中并未考慮。均方根偏差計算表明，ε類湍流模型與實驗值吻合度更高，明顯優(yōu)于ω類湍流模型，這是由于兩類湍流模型求解的湍流黏性系數(shù)不同所致。

3 夾持銅片對實驗結(jié)果的影響

由于實驗過程中金屬鉑的壁面熱流密度無法直接測量，Liu等[7-9]通過加熱功率與金屬鉑的溫度進(jìn)行計算，間接獲得其值。然而實驗過程中金屬鉑不僅對流體進(jìn)行加熱，同時與兩側(cè)的夾持銅片發(fā)生熱傳導(dǎo)，若直接采用數(shù)值模擬所獲得的表面熱流密度與實驗值進(jìn)行對比，將會因為忽略與銅片的導(dǎo)熱而產(chǎn)生較大誤差。為此，在上述數(shù)值模擬與實驗結(jié)果的對比過程中，表面熱流密度均采用等效表面熱流密度，即通過式(9)計算獲得。根據(jù)氦氣外掠平板流動實驗(實驗工況為：流體速度10 m/s，氣體入口溫度313 K，功率指數(shù)增長周期92 ms)的數(shù)值模擬結(jié)果，金屬鉑真實表面熱流密度與等效熱流密度的對比示于圖12。可看出，兩者存在約20%的相對偏差。即當(dāng)前的熱負(fù)載功率指數(shù)提升實驗平臺可通過間接比較等效熱流密度的方法對瞬態(tài)條件下對流傳熱特性進(jìn)行一定程度的驗證，但等效熱流密度計算方法無法剔除夾持銅片導(dǎo)熱引入的誤差。為更好地驗證瞬態(tài)條件下的對流傳熱特性，有必要對該實驗平臺進(jìn)行一定的改進(jìn)。

圖12 金屬鉑的真實表面熱流密度與等效表面熱流密度對比Fig.12 Comparison between actual surface heat flux and equivalent surface heat flux of platinum

4 結(jié)論

通過ANSYS Fluent程序，針對功率指數(shù)提升實驗平臺建立分析模型，選取標(biāo)準(zhǔn)k-ε模型、標(biāo)準(zhǔn)k-ω模型、過渡SST模型及RSM模型4種湍流模型，通過控制y+<1確保所有模型均采用強化壁面處理，對瞬態(tài)過程中的流固耦合換熱進(jìn)行了分析，得到以下結(jié)論。

1) 熱負(fù)載指數(shù)提升過程中，熱功率一部分用于對流換熱而另一部分仍留在固體內(nèi)，且對流換熱占熱功率的比值與熱功率的提升速率有關(guān)，熱功率提升速率越高，對流換熱的占比越低。

2) 瞬態(tài)加熱過程中，加熱壁面附近流體由于溫度的瞬變進(jìn)而引起密度與黏性的變化，其流體速度與溫度分布將不斷發(fā)生變化，即定熱物性假設(shè)與壁面函數(shù)法不適用于流體熱物性隨溫度顯著變化的瞬態(tài)對流傳熱特性的分析。

3) 采用4種湍流模型模擬的等效表面熱流密度均低于實驗值。均方根偏差計算表明，ε類湍流模型的模擬結(jié)果與實驗值吻合度較高，且明顯優(yōu)于ω類湍流模型。

4) 當(dāng)前的熱負(fù)載功率指數(shù)提升實驗平臺可通過間接比較等效熱流密度的方法對瞬態(tài)條件下的對流傳熱特性進(jìn)行一定程度的驗證，但等效熱流密度計算方法無法剔除夾持銅片導(dǎo)熱引入的誤差。為更好地驗證瞬態(tài)條件下的對流傳熱特性，有必要對該實驗平臺進(jìn)行一定的改進(jìn)。