任 逸，喻 宏

(中國原子能科學(xué)研究院，北京 102413)

鈉冷快堆作為第4代核能系統(tǒng)中最成熟的快堆技術(shù)，在全世界已有超過420堆·年的經(jīng)驗[1]，我國已完成中國實驗快堆(CEFR)的建造，目前正在進行中國示范快堆(CFR600)的設(shè)計。無保護失流事故(ULOF)是一種發(fā)生概率極小的嚴重事故工況，它是指失流事故時反應(yīng)堆保護系統(tǒng)因機械故障或操作失誤而不能正常工作，僅能依靠反應(yīng)堆固有安全機制被動響應(yīng)的嚴重事故工況[2]。為提高反應(yīng)堆固有安全性，實現(xiàn)ULOF中控制反應(yīng)性的基本安全功能，示范快堆設(shè)計了采用非能動機制的停堆系統(tǒng)，即液體懸浮式非能動停堆機構(gòu)。

液體懸浮式非能動停堆機構(gòu)由兩部分構(gòu)成，分別是非能動組件和非能動棒驅(qū)動機構(gòu)。非能動組件由固定的組件套筒和可移動的非能動棒構(gòu)成，在反應(yīng)堆處于正常運行工況時，非能動棒懸浮于上工位；當發(fā)生失流事故時，非能動棒自動下落，插入堆芯實現(xiàn)停堆[3]。

動導(dǎo)管位于反應(yīng)堆堆芯液態(tài)鈉出口部位，其作用是為非能動棒驅(qū)動機構(gòu)抓手提供運動通道和導(dǎo)向，液態(tài)鈉自非能動組件出口部位沿動導(dǎo)管向上流動，最終流入冷鈉池。

動導(dǎo)管內(nèi)外兩側(cè)的環(huán)境差異較大，最大溫差約為160 ℃。這樣的溫差所引起的熱應(yīng)力對動導(dǎo)管材料本身是一種考驗，在長期運行時可能會使動導(dǎo)管壽命降低、斷裂韌性下降。而動導(dǎo)管一旦斷裂，將會出現(xiàn)嚴重的事故工況。這種在局部區(qū)域因較大溫差導(dǎo)致的熱疲勞問題已成為核電廠重點關(guān)注的安全問題。例如2008年10月，瑞典奧斯卡港核電廠(Oskarshamn Nuclear Power Plant)3號機組(沸水堆)停堆換料期間在控制棒焊接接頭處發(fā)現(xiàn)了裂縫，其業(yè)主公司Tinoco等[4]進行廣泛研究后，發(fā)現(xiàn)裂縫是由在控制棒和上部管道之間的環(huán)形間隙的冷層流(60 ℃)和旁路流動的熱湍流(276 ℃)熱混合導(dǎo)致的熱疲勞產(chǎn)生的。

由于動導(dǎo)管兩側(cè)的熱邊界條件無法預(yù)先給定，而是受到液態(tài)鈉和壁面之間相互作用的制約，因此動導(dǎo)管和內(nèi)外兩側(cè)的液態(tài)鈉構(gòu)成一個典型的共軛換熱問題，也稱為耦合傳熱問題[5]。本工作采用開源CFD軟件OpenFOAM中的chtMultiRegionSimpleFoam求解器求解該共軛換熱問題。

chtMultiRegionSimpleFoam是OpenFOAM中較成熟且經(jīng)驗證的穩(wěn)態(tài)共軛換熱求解器，已成功應(yīng)用于核能安全分析領(lǐng)域。2013年瑞典皇家理工學(xué)院(KTH)Ignacio[6]基于該求解器，對導(dǎo)致瑞典奧斯卡港核電廠3號機組控制棒破裂的低頻溫度波動現(xiàn)象進行了數(shù)值模擬研究；2017年，比利時馮卡門流體動力學(xué)研究所(VKI)Koloszar等與比利時核能研究中心(SCK-CEN)Keijers等[7]基于該求解器共同開發(fā)了用于MYRRHA熱工水力行為數(shù)值分析平臺——MyrrhaFoam；2016年，意大利米蘭理工大學(xué)核工程系Pini等[8]在該求解器的基礎(chǔ)上添加了內(nèi)熱源項，開發(fā)了新的求解器chtSourceMultiRegionFoam，研究有分布式內(nèi)熱源存在的熔鹽回路的自然循環(huán)問題，且建造了實驗臺架對計算結(jié)果進行了對比驗證。

本文采用Trelis前處理軟件對鈉冷快堆非能動停堆機構(gòu)動導(dǎo)管和其內(nèi)外兩側(cè)的液態(tài)鈉建模并生成全域高質(zhì)量六面體網(wǎng)格，基于開源CFD軟件OpenFOAM中的chtMultiRegionSimpleFoam求解器，解決動導(dǎo)管和其內(nèi)外兩側(cè)的液態(tài)鈉構(gòu)成的共軛換熱問題，得到對動導(dǎo)管數(shù)值分析的溫度分布結(jié)果，為力學(xué)設(shè)計提供重要的輸入條件。

1 定義計算域與網(wǎng)格劃分

定義計算域與網(wǎng)格劃分是CFD計算的第1步，特別是對多個液體、固體域組成的復(fù)雜共軛換熱問題，識別、選取、定義合理的計算域并對應(yīng)劃分高質(zhì)量的網(wǎng)格尤為重要。本工作采用Csimoft公司開發(fā)的先進前處理軟件Trelis進行幾何建模與網(wǎng)格劃分。

1.1 定義計算域

圖1為通過Trelis建立的幾何模型，除簡化少量用于加工目的的倒角以外，基本保持原設(shè)計尺寸和形狀不變。在圖1的基礎(chǔ)上，需進一步確定并生成相應(yīng)的流體和固體計算域。動導(dǎo)管作為中間固體域(定義為Middle)，其內(nèi)側(cè)流體域(定義為Inner)是從非能動組件操作頭進入的液態(tài)鈉，外側(cè)流體域(定義為Outer)是其周圍1圈6個燃料組件操作頭所形成的流體域。外側(cè)流體域的選擇可滿足動導(dǎo)管共軛換熱問題的研究和求解要求，同時可減少網(wǎng)格數(shù)、降低網(wǎng)格復(fù)雜度，提高求解效率。最終形成Inner-Middle-Outer 3部分構(gòu)成的計算域。

圖1 正常運行工況下動導(dǎo)管 及其周圍組件操作頭Fig.1 Guide tube and its surrounding assemblies’ operating heads under normal operation condition

Inner-Middle-Outer所構(gòu)成的整體計算域的規(guī)模為：直徑181.5 mm，高度580 mm。其中中間固體計算域Middle與圖1所示動導(dǎo)管完全相同，對于流體計算域Inner和Outer，通過在對應(yīng)圓柱區(qū)域中去除圖1所示的固體部分得到，具體可通過Trelis中Boolean操作的Substract方法實現(xiàn)，生成的Inner和Outer計算域如圖2所示。

圖2 Inner(a)和Outer(b)流體域的生成Fig.2 Generation of fluid domains of Inner (a) and Outer (b)

1.2 網(wǎng)格劃分

高質(zhì)量網(wǎng)格的生成是CFD計算的前提條件，也是CFD工作中人工工作量最大的部分[9]，在本工作中約占總工作量的70%。本工作采用Trelis對整個計算域生成六面體網(wǎng)格。雖然六面體網(wǎng)格的生成耗時多、難度大，但其相比四面體網(wǎng)格有更好的精度、更高的效率。在六面體網(wǎng)格生成過程中，幾何切割方案的制定是高質(zhì)量網(wǎng)格生成的基礎(chǔ)和難點，特別是OpenFOAM對網(wǎng)格質(zhì)量的要求較高，如對網(wǎng)格的偏斜率、非正交度參數(shù)等指標有要求。

而要針對復(fù)雜幾何沿流動方向進行“掃掠”生成高質(zhì)量的全域結(jié)構(gòu)化/非結(jié)構(gòu)化六面體網(wǎng)格，必須對幾何結(jié)構(gòu)整體進行系統(tǒng)分析，從CFD計算的全局全過程考慮制定幾何切割方案。例如圖3所示Outer的切割方案中切割數(shù)約為300塊。在此過程中Trelis和OpenFOAM需反復(fù)迭代，直到網(wǎng)格質(zhì)量滿足OpenFOAM的網(wǎng)格檢查要求。

圖3 Inner(a)、Middle(b)和Outer(c)切割方案Fig.3 Decomposition schemes of Inner (a), Middle (b) and Outer (c)

Inner-Middle-Outer生成六面體網(wǎng)格的數(shù)量和質(zhì)量參數(shù)列于表1，其中網(wǎng)格已滿足OpenFOAM內(nèi)置的checkMesh網(wǎng)格質(zhì)量檢查工具的要求。由表1可看出，由于Outer區(qū)域幾何形狀復(fù)雜，切割難度大，使得其網(wǎng)格質(zhì)量稍低于Inner和Middle。

2 計算模型

2.1 物理模型

一般的換熱問題均會給定邊界的溫度分布、熱流密度或這兩者的線性組合，而共軛換熱邊界條件不能事先給定，如何構(gòu)成定解問題并且合理地映射邊界值是共軛換熱問題的難點。chtMultiRegionSimpleFoam采用分區(qū)求解法求解共軛換熱問題，即逐個求解各計算域的控制方程，再對相鄰的邊界的溫度(T)和熱流密度(q)耦合，重復(fù)此過程直到所有變量收斂[5]。

表1 六面體網(wǎng)格數(shù)量和質(zhì)量Table 1 Quantity and quality of hexahedral mesh

流體域滿足穩(wěn)態(tài)不可壓縮牛頓流體的控制方程。

連續(xù)性方程：

(1)

動量守恒方程：

(2)

能量守恒方程：

(3)

固體域滿足傳熱方程：

(4)

式中：xi為方向坐標；ui為速度；ρ為流體密度；p為流體壓強；μ為流體動力黏度；h為流體比焓；λl為流體導(dǎo)熱系數(shù)；Tl為流體溫度；Sh為流體內(nèi)熱源；λs為固體導(dǎo)熱系數(shù)；Ts為固體溫度。

在相鄰計算域的交界面，對TW和qW進行耦合，例如對于Ⅰ區(qū)和Ⅱ區(qū)的邊界需滿足[5]：

TW|I=TW|II

(5)

qW|I=qW|II

(6)

此外，對于本工作中耦合邊界面的處理有其特殊性，雖然Middle規(guī)則的幾何形狀可生成全域的結(jié)構(gòu)化六面體網(wǎng)格，但由于Outer和Inner的不規(guī)則結(jié)構(gòu)，如Inner在貫穿在斜面上的3個圓柱形流域，使得六面體網(wǎng)格生成時的掃掠方向受限，導(dǎo)致邊界出現(xiàn)了不能保角映射的非結(jié)構(gòu)化網(wǎng)格，如圖4所示。針對此問題本文采用了OpenFOAM中的任意網(wǎng)格界面(AMI)技術(shù)。AMI技術(shù)可使每個面接受鄰近的面和其交叉部分的權(quán)重的貢獻，可確保在兩個相鄰邊界面間插值且不會出現(xiàn)不連續(xù)問題。

圖4 非保角映射邊界網(wǎng)格示意圖 (Inner-Middle-Outer)Fig.4 Boundary mesh diagram of nonconformal mapping (Inner-Middle-Outer)

2.2 湍流模型

雷諾時均(RANS)模擬方法中的k-ε兩方程模型是工程中最常用的湍流模型之一，在液態(tài)金屬管道流計算中已應(yīng)用且經(jīng)過實驗驗證。中國原子能科學(xué)研究院李淞等[10]應(yīng)用標準k-ε模型對CEFR燃料組件管腳流量分配進行了數(shù)值模擬，且模擬結(jié)果得到實驗驗證。本文采用標準k-ε兩方程模型[11]，湍流脈動動能k和湍流脈動動能的耗散率ε的穩(wěn)態(tài)控制方程如下：

(7)

(8)

(9)

其中，u′i、u′j為速度脈動量。方程組包含5個可調(diào)整的常量，標準k-ε模型采用綜合性數(shù)據(jù)得到的常數(shù)值可廣泛用于各種湍流，Cμ=0.09，σk=1.00，σε=1.30，C1ε=1.44，C2ε=1.92。

此外，關(guān)于Inner和Outer的湍流邊界條件，即湍流入口的k和ε，采用如下常用經(jīng)驗公式[11]計算：

(10)

I=0.16Re-1/8

(11)

(12)

l=0.07L

(13)

其中：Re為主流區(qū)的雷諾數(shù)；U為入口流速；I為湍流強度；L為當量直徑。

2.3 物性參數(shù)

Inner和Outer計算域均為液態(tài)鈉，其物性參數(shù)計算公式[12-13]參考國際原子能機構(gòu)(IAEA)和比利時國家核能研究中心(SCK·CEN)發(fā)布的官方數(shù)據(jù)，具體如下所示。

液態(tài)鈉密度：

ρ=[0.896 606 79+0.516 134 3×

(T×10-3)-1.829 721 8×(T×10-3)2+

2.201 624 7×(T×10-3)3-1.397 563 4×

(T×10-3)4+0.448 668 94×(T×10-3)5-

0.057 963 628×(T×10-3)6]×103

(14)

液態(tài)鈉比定壓熱容：

cp=(38.12-0.69T-2-1.949 3×

10-2T+1.024×10-5T2)/22.99×1 000

(15)

液態(tài)鈉動力黏度：

(16)

液態(tài)鈉熱導(dǎo)率：

λ=99.5-0.039 1T

(17)

Middle采用316不銹鋼，其物性參數(shù)計算公式[12]如下。

316不銹鋼密度：

ρ=8 084-0.420 9T-3.894×10-5T2

(18)

計算時認為動導(dǎo)管密度不變，取動導(dǎo)管平均溫度477 ℃(750 K)時的值7 746 kg/m3作為動導(dǎo)管的密度。

316不銹鋼比定壓熱容：

cp=462+0.134T

(19)

316不銹鋼熱導(dǎo)率：

λ=9.248+0.015 71T

(20)

2.4 邊界條件

本工作計算的工況是正常運行工況，即非能動棒懸浮于上工位。處于正常運行工況的動導(dǎo)管內(nèi)外組件操作頭的結(jié)構(gòu)和入口冷卻劑參數(shù)分別如圖1和圖5所示。冷卻劑入口邊界值包括給定的溫度值和流量值，設(shè)置為質(zhì)量流量邊界條件，該輸入值取自設(shè)計值[14]，屬于獨立的邊界條件；冷卻劑出口設(shè)置為自由流出邊界條件；固體部分均設(shè)置為固壁邊界條件。

圖5 正常運行工況下動導(dǎo)管 內(nèi)外組件操作頭入口冷卻劑參數(shù)Fig.5 Coolant parameters for inlet of inner and outer assemblies’ operating heads of guide tube under normal operation condition

3 數(shù)值求解設(shè)置

3.1 離散格式設(shè)置

3.2 求解和算法設(shè)置

chtMultiRegionSimpleFoam求解器采用Simple算法求解流體域。因壓力場在求解中收斂速度較慢[6]，故設(shè)置多重網(wǎng)格求解器求解壓力場以加快收斂速度。為提高求解的穩(wěn)定性，采用了亞松弛技術(shù)，松弛因子設(shè)置為0.3。此外，各場量內(nèi)迭代求解的殘差設(shè)置為1×10-7。

此外，本工作采用計算域分解法進行并行計算，OpenFOAM中的scotch分解方法是一種自動分解方法，其劃分原則是使每塊分割后的網(wǎng)格邊界面最小化。本工作分別將各流體、固體域劃分為6個計算，即6核并行計算。

4 計算結(jié)果和分析

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

網(wǎng)格無關(guān)性分析是CFD驗證工作中最重要的活動之一[15]。為進行網(wǎng)格無關(guān)性分析，通過保持相同的幾何切割方案和網(wǎng)格質(zhì)量、僅改變網(wǎng)格劃分間隔的方法形成由疏到密的3套網(wǎng)格，具體網(wǎng)格設(shè)置和最細網(wǎng)格尺寸列于表2。

Inner-Middle-Outer在3套網(wǎng)格下的比焓殘差(相對殘差)如圖6所示。圖中的橫線為1×10-4，對于比焓殘差，認為低于1×10-4即收斂?？煽闯觯S著網(wǎng)格的加密，3套網(wǎng)格計算至5 000外迭代步均達到較好收斂。選取3套網(wǎng)格的5 000迭代步時Inner和Outer的出口流量及Middle的溫度分布分別對流體域和固體域進行網(wǎng)格無關(guān)性分析。

表2 用于網(wǎng)格無關(guān)性計算的網(wǎng)格設(shè)置Table 2 Grid setting for grid independence calculation

圖6 不同疏密網(wǎng)格下的Inner(a)、Middle(b)和Outer(c)比焓殘差Fig.6 Specific enthalphy residuals of Inner (a), Middle (b) and Outer (c) under different grids

1) 流體域出口流量

不同疏密網(wǎng)格下Inner和Outer的出口流量列于表3。由表3可看出，隨著網(wǎng)格的加密，Inner和Outer出口流量變化極小，可認為mesh3的流體域網(wǎng)格已滿足網(wǎng)格無關(guān)性要求。

表3 不同網(wǎng)格下的出口流量比較Table 3 Comparison of outlet flows under different grids

2) 固體域溫度場

Inner、Middle和Outer整體溫度分布結(jié)果如圖7所示。

由圖7可看出，動導(dǎo)管外壁面的溫差小于內(nèi)壁面的溫差，因此選取動導(dǎo)管溫度場梯度大的內(nèi)壁面，即Middle和Inner的交界面，高度選取4.38 cm，對不同疏密網(wǎng)格下圍繞該高度一周的溫度場計算結(jié)果進行比較，具體位置如圖8a所示。由于Middle外部受到相對高溫的Outer流體的加熱、內(nèi)部受到相對低溫的Inner流體的冷卻，特別是Inner中非能動組件操作頭的3束分流對Middle的冷卻，形成如圖8b所示的波動的溫度分布?？煽闯觯瑥氖杈W(wǎng)格mesh1到密網(wǎng)格mesh3，其溫度場分布差異趨于收斂，mesh3的固體域網(wǎng)格亦滿足網(wǎng)格無關(guān)性要求。由圖8可得到動導(dǎo)管在周向0.07 m范圍內(nèi)溫差為97 ℃。

圖7 Inner(a)、Middle(b)和Outer(c)溫度分布Fig.7 Temperature distributions of Inner (a), Middle (b) and Outer (c)

圖8 不同網(wǎng)格下動導(dǎo)管溫度場比較Fig.8 Comparison of temperature distribution of guide tube under different grids

由以上對于不同疏密網(wǎng)格的流體域出口流量和固體域溫度場的計算結(jié)果可得出，網(wǎng)格的劃分方案滿足網(wǎng)格無關(guān)性的要求，以下選擇mesh3網(wǎng)格的計算結(jié)果進行分析。

4.2 內(nèi)部流體域軸向溫度分布

為驗證選取的計算域高度合理，有必要對動導(dǎo)管內(nèi)外兩側(cè)的流體軸向溫度進行分析。本文選取動導(dǎo)管內(nèi)部流體(Inner)軸向溫度進行研究，結(jié)果如圖9所示。由靠近出口處的h=420 mm和h=520 mm兩圖可看出，流體已經(jīng)過較充分的混合，在出口的徑向和軸向的溫度變化均很小(溫差約為20 ℃)，可判斷選取的計算域高度合理。

4.3 動導(dǎo)管軸向溫度分布

選取動導(dǎo)管3個軸向方向——Line A、Line B和Line C的溫度分布如圖10所示。由圖10可看出，動導(dǎo)管在軸向0.17 m高度范圍內(nèi)的溫差為80 ℃，且動導(dǎo)管軸向溫度自下而上先急劇降低后再略有升高。前者是因動導(dǎo)管內(nèi)非能動組件操縱頭出口冷鈉對動導(dǎo)管的沖刷冷卻占主導(dǎo)作用，使動導(dǎo)管溫度急劇降低；后者是由于動導(dǎo)管外部的燃料組件出口的熱鈉對動導(dǎo)管的傳熱占主導(dǎo)作用。

5 結(jié)語

本工作通過Trelis前處理軟件生成全域高質(zhì)量六面體網(wǎng)格，基于開源CFD軟件OpenFOAM中的chtMultiRegionSimpleFoam求解器的應(yīng)用，計算鈉冷快堆非能動停堆棒驅(qū)動機構(gòu)動導(dǎo)管和其內(nèi)外兩側(cè)的液態(tài)鈉構(gòu)成的共軛換熱問題，對網(wǎng)格無關(guān)性和計算域高度選取進行了驗證和深入分析，得到了較為可靠的動導(dǎo)管數(shù)值分析的溫度分布結(jié)果，得出動導(dǎo)管內(nèi)壁面周向0.07 m范圍內(nèi)溫差為97 ℃，動導(dǎo)管軸向0.17 m高度范圍內(nèi)的溫差為80 ℃。在這樣窄的范圍內(nèi)該溫差將會產(chǎn)生較大的熱應(yīng)力，數(shù)值分析結(jié)果將為動導(dǎo)管力學(xué)設(shè)計提供重要的輸入條件。

圖9 動導(dǎo)管內(nèi)部流體軸向溫度分布Fig.9 Axial temperature distribution of fluid inside guide tube

圖10 動導(dǎo)管軸向溫度分布Fig.10 Axial temperature distribution of guide tube