基于多重網格的多物理耦合程序開發與驗證

李 壯，孫國民，楊子輝，傅 娟，郁 杰

（1.中國科學院合肥物質科學研究院核能安全技術研究所，合肥 230031；2.中國科學技術大學研究生院科學島分院，合肥 230026）

多物理耦合需要建立合適的網格映射，以實現數據傳輸。陳軍等［1］根據反饋效應的強弱，分別在燃料和慢化劑區域使用一對一映射，在包殼區域采用體積權重的網格映射方式，在Linux系統中完成MCNP5（Monte Carlo N-Particle Transport Code System）和STAR-CCM+（STARComputational Continuum Mechanics+）耦合程序的開發工作。Zhang 等［2］基于OpenMC 和FLUENT 探索自適應平衡算法增強并行性能，對蒙特卡羅模型和CFD（Computational Fluid Dynamics）模型使用相同的單元劃分形式，這種一對一映射處理可以簡化數據映射，但對于處理規模較大的模型，建模和網格劃分需要大量的前處理工作。Xu 等［3］基于數值反應堆物理計算程序NECP-X 和CTF（Coolant-Boiling in Rod Arrays-Two Fluids，COBRA-TF）所開發的耦合程序，在軸向上采用相同劃分方式，CTF網格處的值（冷卻劑溫度和密度）直接被傳遞到相應的NECP-X 網格上，在徑向上傳遞溫度和密度信息時，將CTF 模型每一層的4 個網格的信息平均后傳遞給NECP-X 模型相同分層。Weng 等［4］在堆用蒙特卡羅分析程序（Reactor Monte Carlo code，RMC）的結構化網格和有限元分析求解軟件COMSOL 的非結構化網格之間采取結構化網格實現兩個程序之間的通信。Weng 在蒙特卡羅粒子輸運計算中采用結構化網格計數，燃料邊界處的網格會出現計數偏小的問題。本研究解決了傳統數據映射方法前處理煩瑣，結構化計數網格存在計數偏小的問題，提出了基于多重網格的數據映射方法，并基于蒙特卡羅粒子輸運程序OpenMC［5］和CFD 程序OpenFOAM［6］，采用C++編程語言開發外耦合程序MCOF，通過MCOF 程序在耦合過程自動完成數據交互，實現了靈活高效的數據映射。

1 耦合方法

1.1 耦合工具介紹

OpenMC［5］是一款基于蒙特卡羅方法的粒子輸運計算程序。該程序由美國麻省理工學院研發，2012 年末首次對外公布，它支持反應堆及其系統的高保真建模和中子光子耦合模擬計算，經歷多個版本迭代，目前其計算精度已被廣泛認可。

OpenFOAM［6］是一個完全由C++編程的開源的CFD 求解類庫，其面向對象的程序設計支持用戶根據實際問題對源碼進行修改、擴充。研究中對chtMultiRegionSimpleFoam 進行定制化開發，在能量方程中，新增體積功率源項，使該求解器能夠加載OpenMC 計算得到的體積功率。

1.2 OpenMC 和OpenFOAM 的耦合流程

OpenMC 和OpenFOAM 耦合流程如圖1所示，數據交互由耦合程序MCOF 完成。耦合程序MCOF 的工作流程主要有如下步驟：

（1）初始化中子物理模型的溫度分布和密度分布，調用OpenMC 程序進行粒子輸運計算。

（2）從OpenMC 的輸出文件提取計數結果，有fission、nu-fission 和kappa-fission，并運用式（1）將計數結果歸一化為體積功率。

其中，Pcelli是網格i的功率密度，單位是W·m-3；ncelli是在網格i中每次裂變產生的中子數，為nu-fission 和fission 的比值，單位是neutrons/fission；P是反應堆或組件的熱功率，單位是W；kappacelli是網格i的kappa-fission 計數，單位是eV/source；Qcelli是網格i吸收或者釋放的能量，為kappa-fission 和fission 的比值，單位是eV/fission；Vcelli是網格i的體積，單位是m3；Keff是有效增殖因子［7］；0.974 表示在壓水動力堆的燃料部分沉積97.4%的能量［8］。

（3）將歸一化的非均勻分布功率信息加載到OpenFOAM 求解器的能量方程中進行源項更新，并進行熱工水力計算，獲取溫度分布和密度分布。

（4）根據材料卡ID 更新中子物理模型，重復迭代計算，直至溫度分布和功率分布收斂或達到最大迭代次數。

1.3 基于多重網格的數據映射方法

為了解決傳統數據映射方法前處理煩瑣，結構化計數網格存在計數偏小的問題，本研究提出基于多重網格的數據映射方法進行功率信息和密度、溫度信息映射。其中，功率信息映射通過構建中間獨立的均勻結構化網格實現；溫度和密度信息映射的網格劃分與蒙特卡羅模型的單元劃分一致，并通過此網格完成溫度和密度信息映射。

1.3.1 功率信息映射

功率信息映射如圖2 所示，首先，將計數網格中的計數歸一化為體積功率，歸一化公式如式（1）所示；然后構建中間均勻結構化網格（以下簡稱中間網格）映射體積功率；最后，通過mapFields 程序［6］將中間網格上的功率信息加載到熱工水力模型指定區域，如圖2 中熱工水力模型的燃料區域。

圖2 中間獨立結構化網格功率信息映射圖示Fig.2 Independent structured grid for mapping power information

結構化計數網格在燃料邊界處包含較少的可裂變材料，會出現計數偏小的問題，導致歸一化的功率出現巨大偏差。如圖3 所示，以20×20 的網格未修正結果所示，（a）為OpenMC 計算的裂變率結果，可見在燃料的邊界處的計數偏??；（b）為中間網格的體積功率映射結果，可見邊界處的歸一化功率值偏低；導致（c）和（d）中燃料部分邊界處的功率值偏低。計數網格細化可以最大限度逼近燃料邊界，提升計數精度，但在蒙特卡羅粒子輸運計算的粒子總數一定的情況下，更加精細的網格會導致更高的計數不確定性，因為每個計數網格中的結果由少量源粒子確定，使得計數不確定性增加；而要保證精細的網格具有相同的計數精度，就要增加蒙特卡羅計算的粒子數和計算時間成本［9］。圖4 是100×100 的網格裂變率的計數結果，計算中設置200 個批次（batch），每批次108個粒子，100 核并行計算耗時129 分鐘。由圖4 可見，在邊界附近還是會出現鋸齒狀，數據映射出現燃料邊界局部功率偏低的現象。

圖3 20×20 網格未修正結果Fig.3 20×20 Grid uncorrected results

圖4 100×100 網格裂變率的計數結果（a）和映射結果（b）Fig.4 100×100 grid fission rate（a）and mapping result（b）

文章基于改進薩瑟蘭-霍奇曼多邊形裁剪算法（Sutherland-Hodgman algorithm）［10，11］，在燃料的曲線邊界處通過分割處理、逐邊裁剪，自動辨識裁剪窗口，裁剪出目標區域并計算出可裂變材料的體積修正因子，再進行功率歸一化，改善數據映射的結果。與此同時，本文還展開對計數網格無關性驗證的研究，見表1。

表1 計數網格無關性驗證Table 1 Material thermal property parameters

由表1 的結果可見，經過修正，在燃料邊界處不再出現功率值偏低的地方；經過三種網格的計數統計，10×10 網格、20×20 網格和50×50 網格的歸一化功率最高值分別為3.7×107W·m-3、3.8×107W·m-3、3.8×107W·m-3。由于10×10 網格的計數精度相對偏低，50×50網格的蒙特卡羅粒子輸運計算耗時較長，本文采用20×20 網格驗證進行映射方法和耦合程序驗證。

1.3.2 溫度和密度信息映射

在溫度和密度信息映射方面，首先在中子物理模型中顯式地劃分出若干單元，然后熱工水力模型依據中子物理模型的單元劃分，在熱工水力模型中隱式地劃分出網格單元（以下簡稱偽單元），最后將偽單元的溫度和密度均值傳遞到中子物理模型對應單元，更新對應材料卡的溫度和密度。值得注意的是，該信息映射分為徑向映射和軸向映射，在徑向映射時，根據哈爾濱工程大學的Zhang 等［12］的研究，燃料沿徑向的溫度分布對功率分布的影響極小。因此，本文在處理中子物理模型時，燃料徑向上沒有劃分單元，同時考慮減小溫度映射的誤差，冷卻劑在徑向上進行單元劃分，如圖5 所示。軸向映射時，通過在軸向上劃分出多層來進行映射。層數越多，劃分越精細，蒙特卡羅粒子輸運計算消耗計算資源就越大［13］。

圖5 蒙特卡羅模型徑向單元劃分Fig.5 MC model radial cell division

1.3.3 耦合程序收斂判斷

對于耦合程序計算收斂的判斷問題，伊利諾伊大學的研究者［14］和加州大學的研究者［9］認為可以通過判斷有效增殖因子和溫度是否收斂終止程序。有效增殖因子判斷收斂只能判斷出蒙特卡羅粒子輸運計算達到收斂，使得有效增殖因子和溫度判斷收斂不能反映出中子物理和熱工水力的反饋效應，而溫度和功率判斷收斂能夠反映出中子物理模型的功率空間分布對熱工水力模型的溫度空間分布的影響。因此本文選擇溫度和功率作為耦合程序收斂的判據，判斷收斂的計算式如式（2）和式（3）所示。

其中，N為網格總數；n-1 表示前一次耦合計算；n表示當前耦合計算；celli是網格索引；RMS為每個網格前后兩次耦合計算的溫度或功率差的平方與網格總數的比值；ΔT和ΔP分別是溫度收斂因子和功率收斂因子。耦合迭代10次后，當收斂因子ΔT和ΔP都小于1%時，可視為耦合計算達到收斂。

2 程序開發與驗證

2.1 程序開發

MCOF 耦合程序主體采用C++編程語言完成開發。MCOF 耦合程序主要有三個功能模塊，分別為功率映射模塊、材料更新模塊和耦合模塊。功率映射模塊控制蒙特卡羅程序與CFD程序的功率傳遞，材料更新模塊控制CFD 程序和蒙特卡羅程序溫度密度傳遞，耦合模塊主要包含耦合求解中程序的調用策略和收斂方案。

2.2 程序驗證結果與分析

本研究基于多重網格映射方法完成耦合程序開發，參考單棒模型［15］進行驗證。單棒如圖6 所示，綠色部分是包殼，藍色部分是冷卻劑，紅色部分是燃料，模型的幾何和物性參數均采用參考文獻中的數值［15］。

圖6 IPR-R1 TRIGA 單棒模型Fig.6 IPR-R1 TRIGA single rod model

耦合迭代收斂后的功率和溫度沿軸向中心線分布，如圖7 所示，與參考文獻［15］結果的總體相對偏差在3%以下。圖8 所示模型高度的二分之一處沿徑向的功率分布和溫度分布，與文獻總體偏差小于1%。其中功率分布通過多重網格映射方法和OpenFOAM 的mapFields 程序插值后，與參考值符合良好。耦合計算結果充分反映了中子物理和熱工水力的相互作用，體現了反應堆堆芯的反饋機制。耦合迭代計算過程的功率收斂因子和溫度收斂因子如圖9 所示，最終收斂因子分別為0.45%和0.28%。

圖7 沿軸向功率分布和溫度分布Fig.7 Power distribution and temperature distribution in axial direction

圖9 功率收斂因子和溫度收斂因子隨耦合迭代次數的變化Fig.9 Variation of power convergence factor and temperature convergence factor with the number of coupling iterations

3 結語

針對傳統數據映射方法前處理煩瑣，結構化計數網格存在計數偏小等問題，采用多重網格數據映射方法開展耦合計算研究。通過單棒模型對數據映射方法和耦合程序進行耦合計算驗證，與文獻結果進行了對比，功率分布曲線和溫度分布曲線符合良好，驗證了所提出數據映射方法的正確性。