鈉冷快堆非能動停堆機構動導管共軛換熱數值分析

任 逸，喻 宏

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

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

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

動導管位于反應堆堆芯液態鈉出口部位，其作用是為非能動棒驅動機構抓手提供運動通道和導向，液態鈉自非能動組件出口部位沿動導管向上流動，最終流入冷鈉池。

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

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

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

本文采用Trelis前處理軟件對鈉冷快堆非能動停堆機構動導管和其內外兩側的液態鈉建模并生成全域高質量六面體網格，基于開源CFD軟件OpenFOAM中的chtMultiRegionSimpleFoam求解器，解決動導管和其內外兩側的液態鈉構成的共軛換熱問題，得到對動導管數值分析的溫度分布結果，為力學設計提供重要的輸入條件。

1 定義計算域與網格劃分

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

1.1 定義計算域

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

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

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

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

1.2 網格劃分

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

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

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

Inner-Middle-Outer生成六面體網格的數量和質量參數列于表1，其中網格已滿足OpenFOAM內置的checkMesh網格質量檢查工具的要求。由表1可看出，由于Outer區域幾何形狀復雜，切割難度大，使得其網格質量稍低于Inner和Middle。

2 計算模型

2.1 物理模型

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

表1 六面體網格數量和質量Table 1 Quantity and quality of hexahedral mesh

流體域滿足穩態不可壓縮牛頓流體的控制方程。

連續性方程：

(1)

動量守恒方程：

(2)

能量守恒方程：

(3)

固體域滿足傳熱方程：

(4)

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

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

TW|I=TW|II

(5)

qW|I=qW|II

(6)

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

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

2.2 湍流模型

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

(7)

(8)

(9)

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

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

(10)

I=0.16Re-1/8

(11)

(12)

l=0.07L

(13)

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

2.3 物性參數

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

液態鈉密度：

ρ=[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)

液態鈉比定壓熱容：

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

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

(15)

液態鈉動力黏度：

(16)

液態鈉熱導率：

λ=99.5-0.039 1T

(17)

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

316不銹鋼密度：

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

(18)

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

316不銹鋼比定壓熱容：

cp=462+0.134T

(19)

316不銹鋼熱導率：

λ=9.248+0.015 71T

(20)

2.4 邊界條件

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

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

3 數值求解設置

3.1 離散格式設置

3.2 求解和算法設置

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

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

4 計算結果和分析

4.1 網格無關性分析

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

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

表2 用于網格無關性計算的網格設置Table 2 Grid setting for grid independence calculation

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

1) 流體域出口流量

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

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

2) 固體域溫度場

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

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

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

圖8 不同網格下動導管溫度場比較Fig.8 Comparison of temperature distribution of guide tube under different grids

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

4.2 內部流體域軸向溫度分布

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

4.3 動導管軸向溫度分布

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

5 結語

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

圖9 動導管內部流體軸向溫度分布Fig.9 Axial temperature distribution of fluid inside guide tube

圖10 動導管軸向溫度分布Fig.10 Axial temperature distribution of guide tube