摘要：為解決傳統壓水堆輻照監督理論計算中多次物理建模、粒子輸運計算及效率低的問題，基于蒙特卡羅(MC)方法與離散縱標(SN)方法耦合的正向權重一致共軛驅動重要性抽樣(FW-CADIS)方法，建立壓力容器輻照監督管中子注量計算的加速方法。在某CPR1000機組中開展計算精度和速度的影響因素研究，驗證該方法在不同堆芯參數下的適用性并給出SN輸運模擬參數建議值;通過該機組驗證與確認，結果表明：與直接MC方法相比，中子注量率計算的品質因子(FOM)提高約95~181倍;中子注量計算結果與實測值的相對偏差不超過8%。該方法可顯著提高計算效率，同時滿足工程應用的精度要求。

　　關鍵詞：中子注量;輻照監督;壓水堆;正向權重一致共軛驅動重要性抽樣(FW-CADIS);蒙特卡羅(MC)

　　論文《壓水堆輻照監督理論計算的加速方法及驗證》發表在《核動力工程》，版權歸《核動力工程》所有。本文來自網絡平臺，僅供參考。

　　1 引言

　　壓水堆運行過程中，反應堆壓力容器(RPV)長期受中子輻照會產生輻照損傷，影響結構完整性和使用壽命[1-2]。輻照監督作為評估壓力容器輻照損傷的關鍵手段，通過在壓力容器內壁布置輻照監督管，監測中子注量并推斷壓力容器的實際輻照水平，為設備壽命評估提供依據[3-4]。

　　傳統輻照監督理論計算主要依賴蒙特卡羅(MC)粒子輸運方法，需對每根監督管在每個燃料循環的輻照情況單獨建模計算，存在多次物理建模、粒子輸運計算量大、計算效率低等問題[5-6]。離散縱標(SN)方法計算速度快、易收斂，但精度相對較低[7-8]。正向權重一致共軛驅動重要性抽樣(FW-CADIS)方法通過耦合MC與SN方法，利用SN方法快速獲取減方差參數，可顯著降低MC計算的統計方差，提高收斂速度[9-10]。

　　本文基于FW-CADIS方法，結合堆芯核設計程序PCM、SN程序JSNT和MC程序JMCT，建立壓水堆輻照監督理論計算的加速方法，通過CPR1000機組驗證其準確性和高效性，為工程應用提供技術支撐。

　　2 理論基礎與實現過程

　　2.1 FW-CADIS方法原理

　　FW-CADIS方法核心是通過SN方法計算正向和共軛中子注量率分布，生成MC計算所需的減方差參數(源偏倚因子、權窗因子)，從而加速MC計算收斂，具體步驟如下：

　　1. SN正向輸運計算：建立SN物理模型，開展正向中子輸運計算，獲取粗略的正向中子注量率分布(phi(r, E, Omega))，并根據式(4)確定共軛源分布(q^+(r, E, Omega))：

　　[q^{+}(r, E, Omega)=frac{sigma(r, E, Omega)}{iint phi(r, E, Omega) sigma(r, E, Omega) d E d Omega} ag{4}]

　　式中，(sigma(r, E, Omega))為截面分布。

　　2. SN共軛輸運計算：基于共軛源分布(q^+(r, E, Omega))開展SN共軛輸運計算，獲得共軛中子注量率分布(phi^+(r, E, Omega))，并通過式(5)、(6)計算源偏倚因子(hat{q}(r, E, Omega))和權窗因子(w(r, E, Omega))：

　　[hat{q}(r, E, Omega)=frac{phi^{+}(r, E, Omega) q(r, E, Omega)}{iiint phi^{+}(r, E, Omega) q(r, E, Omega) d r d E d Omega} ag{5}]

　　[w(r, E, Omega)=frac{iiint phi^{+}(r, E, Omega) q(r, E, Omega) d r d E d Omega}{phi^{+}(r, E, Omega)} ag{6}]

　　式中，(q(r, E, Omega))為原始源分布。

　　3. MC正向輸運計算：將減方差參數代入MC物理模型，開展正向中子輸運計算，實現結果快速收斂。

　　2.2 輻照監督計算理論基礎

　　壓力容器輻照監督計算屬于純外源問題，通過堆芯功率分布構造源分布，采用MC程序固定源模式計算輻照監督管和壓力容器的中子注量率分布，結合運行歷史得到累積中子注量，具體計算過程如下：

　　1. 累積中子注量計算：

　　[Psi_{calc}=int phi_{calc}left(t_{i} ight) d t_{i} ag{7}]

　　式中，(Psi_{calc})為累積中子注量理論計算值，(phi_{calc}(t_i))為第(i)個燃料循環的中子注量率計算值，(t_i)為第(i)個燃料循環的輻照時間。

　　2. 超前因子計算：

　　[L=frac{Psi_{ISC, calc}}{Psi_{RPV, calc}} ag{8}]

　　式中，(L)為超前因子，(Psi_{ISC, calc})、(Psi_{RPV, calc})分別為輻照監督管和壓力容器的累積中子注量理論計算值。

　　3. 壓力容器實際中子注量推斷：

　　[Psi_{RPV, mes}=frac{Psi_{ISC, mes}}{L} ag{9}]

　　式中，(Psi_{ISC, mes})為輻照監督管累積中子注量實測值，(Psi_{RPV, mes})為壓力容器累積中子注量實際值。

　　工程中需計算能量大于0.1 MeV和1 MeV的中子注量，用于輻照損傷評估。

　　2.3 加速計算方法實現過程

　　基于FW-CADIS的輻照監督加速計算流程分為4個步驟，如圖1所示：

　　[圖1 加速計算方法實現流程]

　　1. 堆芯中子學參數計算及物理建模：基于堆芯參數和實測通量圖數據，使用PCM程序進行燃耗-輸運耦合計算，得到各循環不同燃耗下的原子密度、水密度、硼濃度、功率分布及裂變中子參數，用于構建粒子輸運模擬的物理模型和源分布(空間分布、能譜)。

　　2. SN正向計算：建立通用SN幾何模型(無需精細建模，采用典型幾何參數)，堆芯物理參數采用循環燃耗加權平均值替代，具體加權公式如下：

　　[M=frac{sum_{i}left(t_{i} × frac{sum_{j} B_{i j} m_{i, j}}{sum_{j} B_{i j}} ight)}{sum_{i} t_{i}} ag{10}]

　　式中，(M)為加權平均參數，(t_i)為第(i)個循環的輻照時間，(B_{ij})為第(i)個循環第(j)個燃耗步的燃耗深度，(m_{i,j})為第(i)個循環第(j)個燃耗步的參數值。

　　3. SN共軛計算：基于SN正向計算結果，構造共軛源分布，開展SN共軛輸運計算，生成MC計算所需的減方差參數。

　　4. MC正向計算：將減方差參數代入MC模型，開展正向中子輸運計算，得到輻照監督管中子注量率分布。

　　3 輻照監督加速計算與驗證

　　3.1 反應堆模型

　　基于某CPR1000機組的真實堆芯物理參數、幾何參數和材料參數，使用JSNT/JMCT程序建立1/4堆芯粒子輸運物理模型(選取含2根輻照監督管的象限)，反應堆徑向部件從內到外依次為燃料組件、圍板、吊籃、熱屏、輻照監督管、壓力容器、保溫層。堆芯燃料組件為AFA3G型(共157個)，采用均勻化模型;輻照監督管采用精細化建模，包含包容器、殼體、底架、力學性能試樣及上中下三個中子探測器。

　　[圖2 CPR1000的計算模型剖面圖]

　　(a)X-Y平面剖面圖;(b)X-Z平面剖面圖

　　3.2 影響因素分析

　　3.2.1 堆芯物理參數影響分析

　　選取水密度、燃耗分布、功率分布、能譜、硼濃度等典型堆芯物理參數，開展敏感性分析，比較加速方法與直接MC方法的中子注量率計算結果，如表1所示。結果表明：在典型堆芯參數變化范圍內，兩種方法的計算偏差均在3倍統計誤差(3σ)以內，證明加速方法在不同堆芯狀態下具有準確性和無偏性。

　　表1 堆芯參數敏感性分析

　　| 參數類型 | 變化類型 | 統計誤差/% | | 計算偏差/% |

　　| | | 加速方法 | 直接MC方法 | |

　　| 水密度變化 | -0.2 g/cm³ | 0.52 | 1.04 | 0.53 |

　　| | -0.1 g/cm³ | 0.57 | 1.16 | -0.56 |

　　| | 0 g/cm³ | 0.62 | 1.30 | -0.87 |

　　| | 0.1 g/cm³ | 0.69 | 1.43 | 0.01 |

　　| | 0.2 g/cm³ | 0.74 | 1.56 | 1.05 |

　　| 燃耗分布 | 第1循環 | 0.63 | 1.35 | 1.18 |

　　| | 第2循環 | 0.63 | 1.36 | 1.90 |

　　| | 第3循環 | 0.62 | 1.35 | 2.33 |

　　| | 第4循環 | 0.62 | 1.35 | 2.54 |

　　| 功率分布 | 第1循環 | 0.62 | 1.30 | -0.87 |

　　| | 第2循環 | 0.64 | 1.33 | -1.95 |

　　| | 第3循環 | 0.63 | 1.29 | 0.32 |

　　| | 第4循環 | 0.63 | 1.30 | -0.48 |

　　| 能譜 | 第1循環 | 0.62 | 1.30 | -0.87 |

　　| | 第2循環 | 0.62 | 1.31 | 0.19 |

　　| | 第3循環 | 0.62 | 1.29 | -1.23 |

　　| | 第4循環 | 0.62 | 1.28 | -1.46 |

　　| | 瓦特譜 | 0.63 | 1.31 | -0.60 |

　　| 硼濃度變化 | -200 ppm | 0.62 | 1.30 | 1.80 |

　　| | -100 ppm | 0.62 | 1.30 | 0.10 |

　　| | 0 ppm | 0.62 | 1.30 | -0.72 |

　　| | 100 ppm | 0.62 | 1.29 | 0.21 |

　　| | 200 ppm | 0.62 | 1.30 | 0.83 |

　　(注：1 ppm=10??)

　　3.2.2 SN模擬參數影響分析

　　SN模擬的離散角度和空間離散尺寸直接影響減方差參數質量，進而影響MC計算效果：

　　1. 離散角度影響：選取S2、S4、S8、S16四種離散角度，空間離散尺寸均為2 cm，分析對MC收斂速度和計算結果的影響(表2)。結果表明：不同離散角度下，中子注量率偏差均在3σ以內;離散角度對SN計算時間影響顯著，S2離散角度的SN計算時間僅148 s，建議優先選擇S2~S4離散角度以節約計算時間。

　　表2 不同角度離散數對計算結果的影響

　　| 角度離散數 | SN計算時間/s | 模擬粒子數為10? | | 模擬粒子數為10? | |

　　| | | 統計誤差/% | 中子注量率偏差/% | 統計誤差/% | 中子注量率偏差/% |

　　| S2 | 148 | 4.76 | 0.08 | 0.68 | 0.49 |

　　| S4 | 279 | 4.49 | 0.51 | 0.64 | -0.11 |

　　| S8 | 718 | 4.14 | -5.56 | 0.62 | 0.06 |

　　| S16 | 2378 | 4.23 | -3.89 | 0.63 | -0.84 |

　　2. 空間離散尺寸影響：選取2、3、5、10、15、20 cm六種空間離散尺寸，角度離散數為S16，分析對MC收斂速度和計算結果的影響(表3)。結果表明：空間離散尺寸越小，MC收斂速度越快，但SN計算時間越長;空間離散尺寸過大(如20 cm)會導致減方差效果下降;建議選擇5~10 cm的空間離散尺寸，兼顧SN計算效率和MC收斂速度。

　　表3 不同空間離散尺寸對計算結果的影響

　　| 空間離散尺寸/cm | 總網格數 | SN計算時間/s | 統計誤差/% | 中子注量率偏差/% |

　　| 2 | 4.0×10? | 19789 | 0.60 | 0.38 |

　　| 3 | 1.2×10? | 7894 | 0.62 | -0.34 |

　　| 5 | 2.2×10? | 2378 | 0.63 | -0.84 |

　　| 10 | 2.8×10? | 2407 | 0.81 | 0.04 |

　　| 15 | 1.2×10? | 2388 | 1.17 | -2.21 |

　　| 20 | 3.5×10³ | 2841 | 3.87 | -3.09 |

　　3.3 驗證與確認

　　在CPR1000機組中開展加速方法的驗證與確認，對比加速方法與直接MC方法的計算效率和精度：

　　1. 計算效率對比：直接MC方法模擬粒子數為10?，CPU核數80個;加速方法中MC模擬粒子數為10?，SN計算角度離散數S4、空間離散尺寸5 cm。結果表明(表4)：加速方法的FOM提高95~181倍，統計誤差均低于1%，滿足工程要求(統計誤差<5%);直接MC方法需約1 h收斂至統計誤差5%，而加速方法僅需30 s以內(圖5)。

　　表4 兩種方法快中子注量理論計算值的對比

　　| 燃料循環 | 模擬方法 | 統計誤差/% | 計算時間/h | FOM | 計算偏差/% | FOM比值 |

　　| 第1循環 | 直接MC方法 | 0.47 | 94.16 | 6.06 | 0.63 | 150 |

　　| | 加速方法 | 0.60 | 0.38 | 907.68 | - | - |

　　| 第2循環 | 直接MC方法 | 0.34 | 133.86 | 8.27 | 1.28 | 113 |

　　| | 加速方法 | 0.60 | 0.37 | 935.59 | - | - |

　　| 第3循環 | 直接MC方法 | 0.33 | 229.45 | 5.16 | 0.30 | 181 |

　　| | 加速方法 | 0.60 | 0.37 | 932.42 | - | - |

　　| 第4循環 | 直接MC方法 | 0.46 | 67.20 | 8.72 | 0.90 | 95 |

　　| | 加速方法 | 0.60 | 0.43 | 826.38 | - | - |

　　[圖5 直接MC方法與加速方法的收斂速度對比]

　　2. 計算精度對比：將累積中子注量理論計算值與實測值對比(表5)，加速方法的計算偏差不超過8%，符合工程級精度要求。

　　表5 累積中子注量理論計算值與實測值的對比

　　| 探測器位置 | 能量范圍 | 計算值與實測值偏差/% | |

　　| | | 直接MC方法 | 加速方法 |

　　| 上 | E>1.0 MeV | 6 | 5 |

　　| | E>0.1 MeV | 2 | 1 |

　　| 中 | E>1.0 MeV | 6 | 8 |

　　| | E>0.1 MeV | 3 | 4 |

　　| 下 | E>1.0 MeV | 8 | 6 |

　　| | E>0.1 MeV | 7 | 7 |

　　4 結論

　　本文基于FW-CADIS方法，耦合MC與SN方法，建立了壓水堆輻照監督理論計算的加速方法，通過CPR1000機組驗證得到以下結論：

　　1. 加速方法的減方差參數在不同堆芯物理參數(水密度、燃耗分布、功率分布等)下具有廣泛適用性，無需重復SN建模和減方差參數計算，可直接用于不同堆芯狀態的MC加速計算。

　　2. 給出SN模擬參數建議值：角度離散數選取S2~S4，空間離散尺寸選取5~10 cm，兼顧計算效率和精度。

　　3. 與直接MC方法相比，加速方法的FOM提高95~181倍，計算時間顯著縮短;中子注量計算結果與實測值偏差不超過8%，滿足工程應用精度要求。

　　該方法可有效解決傳統輻照監督計算效率低的問題，尤其適用于多循環、多迭代設計場景，為壓水堆輻照監督理論計算提供高效可行的技術方案。

　　參考文獻

　　[1] 能源行業核電標準化技術委員會. 輕水冷卻反應堆壓力容器輻照監督：NB/T 20220-2013[S]. 北京：核工業標準化研究所，2013：1-25.

　　[2] 張斌，鄭君蕭，李曉靜，等. 反應堆壓力容器快中子注量計算不確定性分析[J]. 核技術，2018，41(2)：020605.

　　[3] 許懷金. RPV快中子注量率計算評估方法研究[D]. 北京：華北電力大學(北京)，2022.

　　[4] 鄭君蕭. 基于SN方法的反應堆壓力容器快中子注量率計算方法研究[D]. 北京：華北電力大學(北京)，2017.

　　[5] 夏春梅，梅其良，丁謙學，等. DORT程序進行RPV中子注量率計算的可靠性驗證[J]. 核科學與工程，2016，36(3)：329-334.

　　[6] VASILIEV A, FERROUKHI H, ZIMMERMANN M A, et al. Development of a CASMO-4/SIMULATE3/MCNPX calculation scheme for PWR fast neutron fluence analysis and validation against RPV scraping test data[J]. Annals of Nuclear Energy, 2007, 34(8)：615-627.

　　[7] 靳忠敏，陳義學，石生春，等. 基于MCNP的壓力容器快中子注量率計算參數敏感性分析[J]. 原子能科學技術，2011，45(2)：195-199.

　　[8] 石秀安，蘇耿華，包鵬飛. 壓水堆輻照監督管中子注量計算方法改進研究[J]. 核科學與工程，2018，38(4)：585-589.

　　[9] 鄧理鄰，呂煥文，譚怡，等. 輻照監督管中子注量率精細化模型計算方法研究[J]. 核動力工程，2013，34(S1)：84-86.

　　[10] 劉巧鳳，韓靜茹. CAP1400反應堆壓力容器快中子注量獨立審核計算[J]. 中國核電，2018，11(4)：462-465.

　　[11] TIEP N H, HOANG S M T, HARTANTO D, et al. Investigation of the VVER-1000 reactor pressure vessel neutron fluence and displacement per atom using MCNP6[J]. Radiation Physics and Chemistry, 2020, 177：109141.

　　[12] NGUYEN H T, NHU V H P, NGUYEN M T. Calculation of neutron and gamma fluences on VVER reactor pressure vessel[J]. Nuclear Science and Technology, 2016, 6(4)：18-25.

　　[13] ARDEKANI S F G, HADAD K. Monte Carlo evaluation of neutron irradiation damage to the VVER1000 RPV[J]. Nuclear Energy and Technology, 2017, 3(2)：73-80.

　　[14] 李剛，張寶印，鄧力，等. 蒙特卡羅粒子輸運程序JMCT研制[J]. 強激光與粒子束，2013，25(1)：158-162.

　　[15] 鄧力，李剛，張寶印，等. JMCT蒙特卡羅中子-光子輸運程序全堆芯pin-by-pin模型的模擬[J]. 原子能科學技術，2014，48(6)：1061-1066.

　　[16] 鄧力，李瑞，丁謙學，等. 基于JMCT秦山核電廠一期反應堆屏蔽計算與分析[J]. 核動力工程，2021，42(2)：173-179.

　　[17] 李剛，鄧力，張寶印，等. BEAVRS基準模型熱零功率狀態的JMCT分析[J]. 物理學報，2016，65(5)：052801.

　　[18] 劉雄國，鄧力，胡澤華，等. 基于VENUS-Ⅲ國際基準模型的JMCT程序驗證[J]. 計算物理，2016，33(5)：570-580.

　　[19] 楊超，程湯培，鄧力，等. 三維并行程序JSNT對HBR-2裝置的屏蔽計算與分析[J]. 原子能科學技術，2019，53(2)：250-255.

　　[20] 張廣春，程湯培，鄧力，等. 基于JSNT程序的壓水堆屏蔽計算[J]. 強激光與粒子束，2017，29(3)：036020.

　　[21] 史濤. 蒙特卡羅粒子輸運問題中的減方差方法研究[D]. 綿陽：中國工程物理研究院，2018.

　　[22] 王超，鄧力. 基于MCNP多群伴隨計算與正算分析的比較[J]. 原子核物理評論，2012，29(4)：388-394.

　　[23] WAGNER C J, PEPLOW E D, MOSHER W S. FW-CADIS method for global and regional variance reduction of Monte Carlo radiation transport calculations[J]. Nuclear Science and Engineering, 2014, 176(1)：37-57.