數值堆熱工流體程序CVR-PACA驗證及典型應用

王明軍，鞠浩然，趙民富，李偉卿，劉天才，胡長軍，楊 文，田文喜，秋穗正，蘇光輝

(1.西安交通大學 核科學與技術學院，陜西 西安 710049；2.中國原子能科學研究院 反應堆工程技術研究所，北京 102413；3.北京科技大學 計算機與通信工程學院，北京 100083)

數值堆技術是當前核能領域的新興技術，是各核能大國相互爭奪的核科技戰略高地之一。數值堆技術基于高精度數值仿真技術、大數據應用技術和高速數據傳輸技術，應用先進的反應堆物理、熱工水力、安全分析、燃料性能、材料性能、結構力學、耦合接口等計算軟件，面向某一反應堆實體系統進行全方位、全周期的數字化模擬，并采用可視化及人機交互技術直觀展示模擬預測結果，完成反應堆內各物理場現象的研究，是未來核能系統的先進研究平臺。當前時期，美國、歐盟均提出了各自的數值堆技術路線，并大力發展相應的關鍵技術：歐洲主要完成了以SALOME開源軟件平臺為基礎的反應堆多物理場耦合模擬技術，拓展了平臺集成、模型開發、耦合技術、不確定性分析及驗證等方法，對輕水堆安全分析、運行和工程設計給出綜合性解決方案；美國能源部則通過啟動NEAMS計劃，用于部署和發展先進核能技術，計劃支持開發的多種高精度計算分析程序發展迅速，并已在各相關科研領域產生了一定的影響。此外，日本、韓國等核能強國也積極參與到數值堆相關技術的研究工作中。

為實現由核能大國向核能強國的跨越，我國緊跟國際核能科技發展前沿，投入了大量資源進行數值堆相關的專業軟件研發和應用工作。其中，依托國家重點研發計劃項目“數值反應堆原型系統開發及示范應用”的資助，開發了中國數值反應堆原型系統CVR1.0[1]，系統主要包含熱工、物理、材料、燃料、結構力學5種專業工具及所依托的耦合平臺。其中，熱工水力分析工具包括高精細CFD計算工具CVR-PACA及堆芯子通道分析工具CVR-PASA。本文詳細介紹CVR-PACA程序在反應堆堆芯熱工分析方面的典型應用情況——各類棒束通道內湍流流場的分析建模方法、預測結果分析及模型驗證等應用，驗證依托高精度譜元方法的大渦模擬(LES)模型及非穩態雷諾時均(URANS)模型的適用性，多角度展示CVR-PACA在數值堆精細化熱工分析方面的多場景應用能力和計算精度，定量證明該程序具有出任CVR1.0系統熱工計算模塊重要組成部分的能力。

1 數值方法及湍流模型簡介

1.1 譜元方法

譜元方法是一種高精度偏微分方程離散數值方法，是譜方法和有限元方法的融合。譜元法采用高斯多項式函數插值的方式獲得譜元內部的場參數的精細分布，能夠極大降低截斷誤差，提高計算精度。CVR-PACA采用勒讓德正交多項式作為插值多項式基底：

(1)

(D(u(x,y,z)|Ωk),ν)=(f,ν)

(2)

式中：D(u)=f為所求解的偏微分方程；(·,·)為內積算符；ν為試函數空間。當采用伽遼金方法時，由于試函數空間ν及基函數空間為同一正交函數空間，因而要求計算網格必須為全六面體網格，以保障函數空間的正交性。這一要求為空間離散網格劃分過程提出了較大挑戰。

1.2 時間分裂法求解N-S方程

不可壓縮流動N-S方程的一般形式為：

(3)

為敘述方便，將動量方程簡化為以下形式：

(4)

非穩態項應用k階向后差分格式離散：

L(un+1)+fn+1

(5)

式中，bj為k階向后差分格式各離散項前系數。

非線性項N(u)采用顯式外推格式：

(6)

式中，aj為顯式外推格式各離散項前系數。

線性項L(u)進行以下變形：

(7)

式中：μ為流體動力黏度；I為單位矩陣。

外力源項f采用顯式計算方法，式(4)可整理為：

(8)

式中：F(un)為本時刻已知項；除外力源項f外，上標n+1為本時刻待求項；上標小于等于n為已知項。式(8)可采用時間分裂法分3步完成求解：

(9)

(10)

(11)

式(9)為非線性步，用于求解時間相關項。式(10)為壓力步，通過變形獲得以下壓力泊松方程，利用譜元方法數值求解壓力對流場的作用：

(12)

(13)

采用譜元方法結合速度邊界條件即可求得當前步速度分布，進而完成時間分裂法求解瞬態N-S方程的步進計算過程。

1.3 湍流模型

目前，CVR-PACA能夠支持高通濾波大渦模擬(HPF-LES)模型及URANS SSTk-ω模型。其中HPF-LES模型在文獻[2]中有較為詳細的介紹。濾波后的無量綱化N-S方程形式為:

(14)

(15)

URANS模型方面選擇SSTk-ω模型用于計算驗證，該模型已廣泛應用于攪混翼棒束通道流場預測[6-9]，并獲得了較高的認可度。文獻[10]中詳細說明了相關模型的實現過程。

2 典型應用算例

2.1 單棒光棒通道

為驗證HPF-LES模型對棒間湍流攪混現象的預測能力，本文對韓國SOFEL實驗段的單根光棒稠密柵元進行建模[11]，預測流道中的湍流流動攪混特性。為使計算建模區域長度包含4個完整的子通道間湍流脈動擬序結構，保障計算結果的正確性[12]，研究將建模區域長度設為40Dh，Dh為流道的水力直徑。軸向設置周期性邊界使計算區域封閉。計算預測的速度場分布如圖1所示。可以看到當流場達到充分發展流態時，相鄰子通道間的窄縫處存在較為明顯的橫向脈動現象。

a——猝發過程；b——充分發展流態圖1 猝發及充分發展時的湍流無量綱速度場Fig.1 Turbulent dimensionless velocity field during burst and fully developed condition

算例對節徑比P/D=1.03、1.06，雷諾數Re=20 500條件下的湍流脈動現象進行了定量分析，并與實驗結果進行對比，結果如圖2所示。可以直觀看到，HPF-LES模型能準確預測湍流脈動流場細節，獲取各脈動峰，脈動頻率計算方面也符合較好。為實現對湍流脈動現象的進一步分析，采用湍流脈動系數β定量評估子通道間的湍流攪混能力：

a——節徑比1.03，Re=20 500；b——節徑比1.06，Re=20 500圖2 子通道間窄縫中心處的橫向湍流脈動速度瞬時值Fig.2 Transverse transient velocity at center of narrow gap between adjacent sub-channels

(16)

研究整理了P/D=1.06條件下β隨Re的變化，并與理論推導值和實驗測量值進行對比，結果如圖3所示。可以看到，該結果定量驗證了HPF-LES模型在預測單棒通道內湍流流場方面的精確程度，仿真計算結果與實驗值及理論推導規律符合良好，采用相同的模型參數可外推至其他相似幾何條件用于準確預測相關湍流脈動現象。

圖3 湍流脈動系數與Re的關系Fig.3 Relationship between turbulent mixing coefficient and Reynolds number

2.2 3×3稠密柵元光棒棒束通道

為進一步分析稠密柵元光棒棒束通道內更細致的湍流攪混現象，研究利用3×3棒束通道對不同類型子通道間的攪混現象進行定量分析。計算所采用的網格生成策略與2.1節相似：近壁面第1層網格高度嚴格滿足y+<1，軸向長度為40Dh，流動出入口設為周期性邊界。P/D=1.06、Re=15 000時計算獲得的湍流流場瞬時值如圖4所示。

1——中心通道間脈動速度提取點；2——邊通道間脈動速度提取點；3——邊通道與中心通道間脈動速度提取點；4——角通道與邊通道間脈動速度提取點圖4 3×3光棒棒束子通道內的無量綱速度場分布Fig.4 Dimensionless velocity field distribution in sub-channel of 3×3 bare rod bundle

由圖4中縱截面可看出，不同子通道間的湍流脈動速度大小差異較為明顯。為定性研究該現象，選取圖4中標出的4處窄縫中心位置提取橫向瞬時速度，用于定量對比。圖5示出P/D=1.06時不同Re條件下不同位置湍流脈動速度瞬時值的對比。對比同一子圖各曲線可知，湍流脈動速度的幅值和頻率受Re影響顯著，可采用監測脈動速度的方式初步判斷流體流動的狀態(層流、轉捩、充分發展湍流流態)；對比不同子圖中的相同Re條件下的曲線可知，湍流脈動速度幅值和頻率受監測點兩側子通道幾何形狀影響顯著。在中心通道間的湍流脈動現象最為顯著，臨近角通道相連窄縫處的脈動攪混能力最弱。

a——點1處；b——點2處；c——點3處；d——點4處圖5 各監測點的脈動速度瞬時值Fig.5 Transient turbulent mixing velocity at different monitoring points

研究采用快速傅里葉變換(FFT)方法對上述脈動速度進行了頻域分析，獲得了相應的脈動峰值特征頻率，P/D=1.06、Re=20 500時的結果如圖6所示。由圖6可見：脈動峰值受幾何位置影響顯著，監測點距離角通道越近，脈動峰值越低。原因在于，相比于中心通道內的流體渦旋，邊角通道內的渦旋受幾何限制更為顯著，通道較為局促，流動摩擦阻力較大，渦旋含能較低，渦旋旋轉頻率減弱。但由于幾何形狀難以采用準確的數學參數予以反映，因而目前可直接采用修正系數的方法反映幾何形狀對脈動頻率的影響。

a——點1處；b——點2處；c——點3處；d——點4處圖6 各監測點脈動特征頻率Fig.6 Characteristic frequency at different monitoring points

對于湍流脈動引起的子通道間攪混效應，同樣采用無量綱量湍流脈動系數β進行定量分析，P/D=1.06時的仿真預測結果列于表1。由表1可見，幾何條件對湍流攪混效應同樣具有較為顯著的影響，角通道間隙處的攪混系數較中心通道處降低25%左右。因此，在采用計算獲得的湍流攪混系數指導子通道程序計算子通道間橫向攪混計算時，與脈動頻率類似，可適當引入修正系數，提高其建模精度；或可直接采用角通道臨近位置所獲得的脈動系數結果進行計算，以保證計算過程的保守性。

表1 湍流脈動系數隨不同監測點和Re的分布Table 1 Turbulent mixing coefficient distribution case at different monitoring points and Re

2.3 帶攪混翼棒束通道

為驗證CVR-PACA在預測攪混翼下游湍流流場的精確程度，本文應用該程序對攪混翼棒束通道基準題Matis-H進行了建模計算[13]。幾何建模方面，計算域選取了中心3×3棒束所圍成的2×2子通道作為建模對象；周向采用4對周期性邊界用于計算域的封閉。速度入口邊界設定方面，本次計算測試了簡易速度邊界設置方案對計算結果的影響。簡易速度邊界的具體設置方案如下：1) 定位格架上游光棒長度不進行加長處理，入口段光棒長度僅設置為1Dh；2) 速度入口面上的速度分布采用層流流動的速度分布，具體為：

(17)

網格劃分方面，研究采用網格分區劃分結合網格分裂的方法建立復雜結構的全六面體網格。分區方案如圖7a所示。為滿足大渦模擬湍流模型計算要求，在生成邊界層網格時，需嚴格滿足近壁面第1層網格高度y+<1的條件。最終建立的六面體網格的部分結果如圖7b、c所示。

a——網格劃分策略；b——定位格架條帶部分網格；c——定位格架下游光棒部分網格圖7 3×3棒束六面體計算網格Fig.7 Hex mesh of 3×3 rod bundle

不同模型攪混翼下游不同位置處時均速度場及脈動速度場計算結果與實驗值的對比如圖8、9所示。由圖8，圖9a、b可看到，在時均場預測方面，HPF-LES模型及URANS SSTk-ω模型均與實驗值較為接近，均能較好地預測撕裂式攪混翼下游的時均速度分布；然而，脈動場方面，由圖9c、d看到，URANS模型的預測能力明顯不足，HPF-LES模型準確預測了脈動速度分布的峰值位置，在脈動場預測方面具備一定能力。

為驗證簡易速度邊界設置方法配合HPF-LES模型的聯合計算方案對攪混翼下游流場的預測能力，將該方案計算結果與基準題sp05算例[13]的計算結果進行了定量對比。sp05算例與本研究采用相同的幾何建模方案，不同點在于sp05算例采用RSM-SSG RANS模型完成計算，入口段采用延長光棒長度的方式制造充分發展湍流流態支撐計算。圖9示出本文研究與sp05算例的流場對比，可以看處，HPF-LES模型在時均速度場預測方面更為準確，脈動速度場預測方面，HPF-LES模型同樣能夠更好地捕捉相關細節，預測相關峰值，因而HPF-LES模型配合簡易速度邊界設置的聯合計算方案具有一定可行性。

a——攪混翼下游0.5Dh處y方向速度分量；b——攪混翼下游1.0Dh處x方向速度分量；c——攪混翼下游1.0Dh處z方向速度分量；d——攪混翼下游4.0Dh處z方向速度分量圖8 時均速度場綜合對比Fig.8 Comprehensive comparison of time-averaged velocity

a——攪混翼下游1.0Dh處時均速度z分量；b——攪混翼下游4.0Dh處時均速度z分量；c——攪混翼下游1.0Dh處脈動速度z分量；d——攪混翼下游4.0Dh處脈動速度z分量圖9 不同入口處理方案下攪混翼下游流場的對比Fig.9 Comparison of downstream flow field of mixing vane under different inlet treatment schemes

為進一步探究HPF-LES模型與簡易入口設置方案匹配的具體原因，截取了定位格架處的速度剖面進行研究，如圖10所示。由圖10可知：采用HPF-LES模型計算時，定位格架的下沿作為產生的前臺階幾何形狀能夠猝發湍流流態的形成，從而快速促進湍流流態的發展；同時，湍流在流經定位格架條帶結構時，湍流流態迅速發展，流經攪混翼附近時已轉化為充分發展狀態；對比而言，由于URANS模型很大程度上均勻化了流場細節，難以利用條帶結構促進湍流渦的形成，因而不能很好地猝發湍流流型，獲得的脈動細節也不夠充分。

a——HPF-LES模型；b——URANS SST k-ω 模型圖10 HPF-LES模型及URANS模型的瞬時無量綱速度分布Fig.10 Transient dimensionless velocity distribution calculated by HPF-LES model and URANS model

圖11示出采用HPF-LES模型及URANS模型計算得到的攪混翼下游湍動能云圖及橫向速度矢量分布的預測結果，可以看到，HPF-LES模型預測結果顯示出了湍動能的分布規律：在與軸向垂直的橫截面上，湍動能一般集中分布在存在大型渦旋的子通道中央，同時可觀察到湍動能隨小渦旋在子通道間遷移的特性；軸向上，隨著橫截面到攪混翼的距離增加，流動的混亂程度減弱，湍動能也逐漸減弱。URANS模型不能很好預測這一分布規律，需采用合理的入口條件支撐計算。

a——HPF-LES模型，0.5Dh；b——URANS模型，0.5Dh；c——HPF-LES模型，4Dh；d——URANS模型，4Dh圖11 攪混翼下游不同截面處湍動能及橫向時均速度矢量分布Fig.11 Turbulent kinetic energy and time-averaged transverse velocity profiles at different cross sections of mixing vane downstream

2.4 帶繞絲棒束通道

研究KALLA-KIT 19棒束組件實驗，對帶繞絲棒束通道內的熱工現象進行了計算分析，以驗證建模方案的可行性和精確性。實驗采用液態鉛鉍合金作為工質。繞絲樣式采用簡化的垂直連接方式，其詳細幾何參數參見文獻[14]。計算采用網格分裂方案實現六面體網格生成，最終的網格劃分結果如圖12所示。

a——四面體網格總覽；b——截面網格細節圖12 19棒繞絲組件網格示意圖Fig.12 Mesh scheme of 19 wire-wrapped rod bundles

計算采用無量綱化N-S方程進行。為還原KALLA-KIT實驗采用的定熱流密度的加熱方案，本文采用固定入口流體溫度和壁面熱流密度的方式設定邊界條件，用于支持能量方程求解。壁面熱流密度的無量綱化方面采用下式實現：

(18)

式中：qw為壁面絕對熱流密度；ρ流體密度；T0為特征溫度；v0為特征速度，本算例設為入口平均速度；q′w為無量綱化壁面熱流密度；cp為流體比定壓熱容。為簡化假設，設定熱流密度均勻布置于燃料棒與繞絲的表面。

圖13、14分布示出19棒繞絲算例的速度與壓力瞬時值分布及壓力分布的細節。壓力細節方面，由于繞絲沿各燃料棒順時針繞動，相應的導流特性使各截面上的壓力峰值區域也沿軸向按順時針變化；另外由于采用URANS模型計算，一定程度上反映了瞬時流場的流動細節，因而壓力場不盡光滑平順。

圖13 19棒繞絲組件的無量綱速度(a)及無量綱壓力(b)分布Fig.13 Dimensionless velocity (a) and dimensionless pressure (b) profiles in 19 wire-wrapped rod bundles

壓力分布：a——z=40Dh；b——z=80Dh；c——z=120Dh圖14 19棒繞絲組件的無量綱壓力分布細節Fig.14 Dimensionless pressure distribution detail in 19 wire-wrapped rod bundles

圖15示出z=20Dh處的瞬時溫度分布與橫向速度矢量分布。由圖15可見，截面內溫度分布受流場影響顯著。由于子通道內部存在渦旋，致使通道內部溫度分布較為均勻；同時，繞絲的位置變化通過影響速度壓力分布，最終對溫度分布及峰值點存在影響：由于溫度梯度減小，液態金屬導熱輸熱能力減弱，致使繞絲與燃料棒焊接位置處的傳熱性能更為惡化，在繞絲的背離面，由于子通道內壓力較低，使得局部流動受阻，背離面焊接處換熱惡化，因而出現溫度峰值。迎風面由于壓力升高，促進局部流體流動，因而導致子通道內溫度梯度增大，改善換熱，局部溫度未達到峰值水平。

圖15 19棒繞絲組件的無量綱溫度與橫向速度矢量分布Fig.15 Dimensionless temperature and transverse velocity vector profiles in 19 wire-wrapped rod bundles

上述流場流動特性也可由圖16直觀觀察：相鄰子通道1和子通道2的流線起始端由紅色圈及橘黃色圈標注。其中子通道1為背離面子通道，子通道2為迎風面子通道。可以看到迎風面子通道流體在黃色方框處受繞絲擠壓有所偏移，在下游處有部分流體壓入子通道2；同時由于子通道2中存在渦旋，對流動有一定的約束作用，使得紅色圈內的流線初始呈直線分布(深藍色流線簇)；隨著流體流動，由于紅色箭頭所指繞絲的背離，形成局部真空，子通道2中的部分流體被吸入該繞絲的背離子通道中，同時由于子通道2中的繞絲旋入子通道1中，對流體產生了擠壓，因而在后段部分更多的流體被壓入上層子通道中，使得流線分布范圍更加寬泛，流體在子通道間的攪混更為復雜充分。

圖16 19棒繞絲組件的局部流場流線示意圖Fig.16 Scheme of local streamline in 19 wire-wrapped rod bundles

綜合分析可知，采用瞬態URANS模型計算能獲得繞絲棒束內更多的流動和換熱細節，為驗證設計的熱工特性提供更加豐富的數據，更為精準地指導相關設計的優化改進。

3 結論

本文采用CVR-PACA對單棒光棒通道湍流流場、3×3光棒棒束湍流流場、Matis-H壓水堆棒束通道基準題、19棒帶繞絲組件通道湍流流場進行了仿真計算，驗證了基于高精度譜元數值方法的HPF-LES模型及URANS SSTk-ω模型的適用性和計算精度。模擬過程采用網格分裂技術，結合混合網格生成技術實現了帶攪混翼棒束通道、帶繞絲棒束通道等復雜幾何結構的全六面體網格劃分，實現了利用高精度數值方法對各類堆芯組件內流動換熱現象預測的典型應用。通過計算分析，可得到以下結論。

1) 對照實驗結果，HPF-LES模型能精確計算光棒稠密柵元棒束通道內的湍流流場，精細預測子通道間的湍流攪混現象。該現象受Re及幾何條件影響顯著，計算獲得的相關結論可用于子通道程序的精細化模型開發及修正。

2) 配合高精細劃分的全六面體網格，HPF-LES模型及URANS模型均可較好地預測帶有攪混翼的定位格架棒束通道下游的時均流場。對于湍流脈動流場，HPF-LES模型受入口段影響較小，可以準確預測脈動速度峰值，但URANS模型受入口段流態影響顯著，因此采用具有充分發展湍流流態的入口對于URANS模型計算極為重要。

3) 研究實現了繞絲棒束的全六面體網格的生成策略，采用URANS模型預測了繞絲幾何條件下的流動換熱現象，基于高精度CFD方案的URANS模型能獲得子通道內的精細化壓力分布和子通道內的旋渦分布情況，獲得的流場結果與溫度場分布規律契合較好。通過觀察流線能夠更為直觀地了解繞絲棒束通道內的攪混特性。

通過上述分析方法能夠獲得更精確的流動和換熱細節，為驗證設計的熱工特性提供更加豐富的數據支撐，為方案優化改進提供定量反饋。