趙波 曹欣榮 李權(quán) 齊飛鵬 黃永忠 馬強 陳浩 陳飛飛

DOI：10.16660/j.cnki.1674-098x.2104-5640-5601

摘? 要：本文基于格林函數(shù)節(jié)塊法程序TNGFM與子通道程序COBRA-EN開發(fā)了三維瞬態(tài)核熱耦合計算程序，利用壓水堆彈棒算例對核熱耦合數(shù)值穩(wěn)定性進行了分析，并提出了改進Picard迭代方法。數(shù)值模擬結(jié)果表明，減小計算時間步長可提高耦合計算數(shù)值穩(wěn)定性，但計算效率較低，在物理計算外迭代中引入近似反饋可有效提高耦合收斂性能并大大減少耦合計算迭代次數(shù)。

關(guān)鍵詞：瞬態(tài)? 核熱耦合? 數(shù)值穩(wěn)定性? 計算效率

中圖分類號：TL351.1? ? ? ? ? ? ? ? ? ? ? ?文獻標(biāo)識碼：A文章編號：1674-098X（2021）05（a）-0090-05

Research on Numerical Stability of Transient Neutronics/Thermal-Hydraulic Coupling in Reactor Core

ZHAO Bo1，2? CAO Xinrong2? LI Quan1? QI Feipeng1? HUANG Yongzhong1? MA Qiang1? CHEN Hao1? CHEN Feifei1

（1. Science and Technology on Reactor System Design Technology Laboratory， Nuclear Power Institute of? China， Chengdu， Sichuan Province， 610213? China; 2. Harbin Engineering University， Harbin， Heilongjiang Province， 150001? China）

Abstract： A three-dimensional transient neutronics/thermal-hydraulics coupling calculation program is developed based on Greens function nodal method program TNGFM and subchannel thermal-hydraulic analysis program COBRA-EN. The numerical stability of neutronics/thermal-hydraulics coupling is analyzed with a PWR rod ejection example and an improved Picard method is proposed.The numerical simulation results show that reducing the time step can improve the numerical stability of coupling computation with a low efficiency. Introducing approximate feedback into external iteration of neutronics can effectively improve the coupling convergence performance and greatly reduce the number of iterations of coupling computation.

Key Words： Transient; Neutronics/thermal-hydraulics coupling; Numerical stability; Computational efficiency

隨著高性能計算機技術(shù)的迅速發(fā)展，物理-熱工耦合計算趨向于高保真技術(shù)發(fā)展，采用三維物理-熱工耦合分析能更精確模擬瞬態(tài)或事故過程，提高瞬態(tài)工況下物理及熱工參數(shù)分布計算的準(zhǔn)確性。近年來，國內(nèi)外針對物理-熱工耦合問題開展了大量研究，并開發(fā)了多套耦合程序，這些程序多采用算符分裂或Picard迭代等方法求解，容易出現(xiàn)數(shù)值不穩(wěn)定和收斂速度較慢等問題。本文基于格林函數(shù)節(jié)塊法程序TNGFM與子通道程序COBRA-EN開發(fā)了瞬態(tài)核熱耦合程序，對瞬態(tài)核熱耦合數(shù)值穩(wěn)定性進行了分析，并提出了改進迭代算法改善瞬態(tài)核熱耦合收斂性能。

1? 耦合方法

堆芯瞬態(tài)核熱耦合需解決耦合實現(xiàn)方式、耦合時間步算法、耦合收斂算法等3個關(guān)鍵問題。

1.1 耦合實現(xiàn)方式

耦合實現(xiàn)方式從數(shù)值求解角度可分為強耦合和弱耦合。強耦合是指將物理-熱工方程聯(lián)立求解，如近幾年發(fā)展的Jacobian-free Newton-Krylov（JFNK）[1]方法即可實現(xiàn)耦合變量的同步更新，但實現(xiàn)較為復(fù)雜。目前大多數(shù)核熱耦合程序如COBAYA4/COBRA-TF[2]、RMC/CTF[3]仍以弱耦合為主，程序獨立求解物理-熱工方程，通過交換邊界數(shù)據(jù)實現(xiàn)耦合，能充分利用現(xiàn)有程序，實現(xiàn)相對容易。

1.2 耦合時間步算法

耦合時間步算法主要包括物理與熱工程序各自計算時間步長、耦合參數(shù)傳遞時間點的選取及耦合程序時間步推進格式。物理與熱工程序時間步長的選取通常取決于程序自身的計算精度，較為簡便的方法為各程序選取相同的時間步長，并在每個時間步內(nèi)進行耦合參數(shù)的傳遞。耦合程序時間步推進格式通常分為顯式與隱式耦合[4]，顯式耦合是指在單個時間步內(nèi)，物理與熱工程序采用一次通過方式計算，不進行迭代，計算量較小但精度較低。隱式耦合則需要物理與熱工程序在單個時間步內(nèi)進行反復(fù)迭代直到滿足相應(yīng)收斂準(zhǔn)則為止，計算量較大但能保證耦合邊界條件一致且穩(wěn)定性更好。

1.3 耦合收斂算法

傳統(tǒng)的核熱耦合收斂算法包括OS方法[5]、固定點迭代法（Picard方法）[6]，它們都是先由物理計算得到堆芯內(nèi)功率分布，將數(shù)據(jù)傳遞給熱工程序計算得到堆芯內(nèi)溫度、密度等熱工參數(shù)，再由反饋計算模塊更新物理計算所需截面后進行物理計算，如此反復(fù)迭代直至計算收斂。該算法由于流程簡單較易實現(xiàn)，但也容易出現(xiàn)計算震蕩乃至發(fā)散的現(xiàn)象，目前主要通過引入松弛因子以保證耦合計算的收斂[7]。

在實際核熱耦合計算中，由于初始物理及熱工參數(shù)均為猜測值，沒有必要在每次物理計算后都進行詳細的熱工分析計算再進行反饋計算，尤其是瞬態(tài)計算時此種方式會大大降低計算效率，并可能導(dǎo)致計算不收斂。由此本文提出改進Picard方法，即先近似建立熱工參數(shù)與截面的函數(shù)關(guān)系，在物理源迭代的外迭代中引入該近似反饋，在物理計算收斂后進行詳細的熱工分析計算得到溫度、密度等分布，再更新截面參數(shù)進行核熱耦合計算，如此反復(fù)迭代直至計算收斂。

2? 計算模型

2.1 耦合方案

本文采用隱式耦合方式開發(fā)了三維瞬態(tài)核熱耦合程序TNGFM-C，在每個時間步內(nèi)，物理程序求解堆芯功率分布傳遞給熱工分析程序計算得到燃料溫度、水溫度、水密度等熱工參數(shù)，然后反饋計算模塊根據(jù)最新熱工參數(shù)更新截面?zhèn)鬟f給物理程序，計算流程如圖1所示。

2.2 數(shù)值模型

本文選取某PWR壽期初三維彈棒算例進行瞬態(tài)核熱耦合計算分析，圖2為該堆芯控制棒組件布置示意圖，計算初始條件為：熱態(tài)滿功率，控制棒組A、B、C、S提出225步，D棒組提出156步，考慮最不利的工況，即最大價值的控制棒組D中的L07彈出堆芯，在瞬態(tài)程序中設(shè)置該束控制棒彈出速度為1800cm/s，不考慮觸發(fā)停堆，耦合計算時間步長為0.01s。由于TNGFM為粗網(wǎng)節(jié)塊法程序，COBRA-EN采用全堆芯分析模型，因此在計算時徑向劃分均以組件為單位，軸向劃分20層，考慮徑向及軸向反射層，共劃分為4248個物理節(jié)塊，121個子通道。

3? 數(shù)值分析

利用TNGFM-C耦合程序?qū)λ矐B(tài)彈棒工況進行了計算，分析了采用傳統(tǒng)Picard方法時耦合計算的數(shù)值穩(wěn)定性，對耦合計算時間步長敏感性進行了分析，并比較了改進Picard方法對耦合收斂性能及計算效率的影響。

3.1 Picard方法

圖3為采用傳統(tǒng)Picard方法及引入松弛因子時堆芯相對功率隨耦合迭代次數(shù)的變化。從圖中可看出，在單個時間步耦合迭代初期，由于堆芯功率急劇上升，在熱工反饋作用下耦合計算存在明顯震蕩現(xiàn)象，在引入松弛因子后有所改善，但每個時間步內(nèi)耦合計算仍需多次迭代才能收斂。同時隨著時間步的推進，堆芯功率變化程度減緩，在單個時間步內(nèi)耦合收斂所需迭代次數(shù)逐漸減少，但引入松弛因子并未明顯減少耦合計算總迭代次數(shù)。

3.2 時間步長敏感性分析

在彈棒事故中，在單個時間步內(nèi)由于引入較大正反應(yīng)性導(dǎo)致堆芯功率劇烈變化，使得耦合計算出現(xiàn)明顯震蕩，因此本文分析比較了不同時間步長對核熱耦合計算數(shù)值穩(wěn)定性的影響。圖4分別為時間步長為0.01s、0.005s、0.002s及0.001s時堆芯相對功率隨耦合迭代次數(shù)的變化，從圖中可看出，時間步長越小，耦合迭代計算過程中功率震蕩幅度越小，這是由于時間步長越小則單個時間步內(nèi)引入的反應(yīng)性越小，耦合的反饋效應(yīng)相對越弱，使得單個時間步內(nèi)耦合收斂所需的迭代次數(shù)越少，但時間步長越小整體耦合計算所需時間越長，因此在實際計算中需綜合考慮計算穩(wěn)定性及效率，選擇合適的時間步長。

3.3 改進Picard方法

由前文分析可知，采用傳統(tǒng)Picard迭代法進行瞬態(tài)耦合計算時存在明顯震蕩現(xiàn)象，引入松弛因子后有所改善，但計算效率未有明顯提升，采用較小時間步長雖能提升計算數(shù)值穩(wěn)定性，但計算效率明顯降低，因此本文提出改進Picard方法，即在物理計算源迭代的外迭代過程中引入近似反饋計算，物理計算收斂后進行熱工及反饋計算。圖5為采用改進Picard方法時堆芯相對功率隨耦合迭代次數(shù)的變化，從圖中可看出，由于在物理源迭代中完成了近似耦合計算，物理熱工程序間耦合計算迭代次數(shù)大大減少，同時計算數(shù)值穩(wěn)定性大大提高。

表1為分別采用傳統(tǒng)與改進Picard迭代方法時迭代次數(shù)及計算時間對比，從表中可看出，采用傳統(tǒng)Picard迭代方法時，物理計算迭代次數(shù)較少，但耦合迭代次數(shù)較多，采用改進Picard迭代方法時由于在物理計算外迭代中引入近似反饋，物理計算迭代次數(shù)較多，但耦合迭代次數(shù)大大減少，綜合耦合計算時物理熱工程序各自計算時間及數(shù)據(jù)傳遞時間，采用改進Picard迭代方法耦合計算效率相對較高。

4? 結(jié)語

本文以格林函數(shù)節(jié)塊法程序TNGFM與子通道程序COBRA-EN為基礎(chǔ)開發(fā)了三維物理-熱工瞬態(tài)耦合程序，利用壓水堆壽期初彈棒算例對耦合計算數(shù)值穩(wěn)定性進行了分析，結(jié)果表明，在瞬態(tài)核熱耦合計算時，采用傳統(tǒng)Picard迭代方法存在明顯數(shù)值震蕩現(xiàn)象，引入松弛因子后有所改善但計算效率無明顯提升，同時減小時間步長可提高數(shù)值穩(wěn)定性但計算效率較低，采用改進Picard迭代法能大幅減少耦合迭代次數(shù)，提高計算效率，且能改善耦合迭代過程中的數(shù)值震蕩現(xiàn)象。

參考文獻

[1] 周夏峰，李富.基于節(jié)塊展開法的Jacobian-Free Newton Krylov聯(lián)立求解物理-熱工耦合問題[J].物理學(xué)報，2016（65）：092801-1-092801-8.

[2] N.Garcia-Herranz， D.Cuervo， et al.Multiscale neutronics/thermal-hydraulics coupling with COBAYA4 code for pin-by-pin PWR transient analysis[J].Nuclear Engineering and Design，2017，321：38-47.

[3] Yugao Ma， Shichang Liu，et al. RMC/CTF multiphysics solutions to VERA core physics benchmark problem 9 [J].Annals of Nuclear Energy，2019，133：837-852.

[4] 廖承奎.三維節(jié)塊中子動力學(xué)方程組的數(shù)值解法及物理與熱工-水力耦合瞬態(tài)過程的數(shù)值計算的研究[D].西安：西安交通大學(xué)，2002.

[5] Wang Xianmeng，Wu Mingyu et al.Multi-physics coupling simulation in virtual reactors[J].Simulation：Transactions of the Society for Modeling and Simulation International，2019：1-16.

[6] Jincheng Wang，Qin Wang，et al.Review on neutronic/thermal-hydraulic coupling simulation methods for nuclear reactor analysis[J].Annals of Nuclear Energy，2020，137：11-12.

[7] Juanjuan Guo，Shichang Liu，et al.Versatility and stabilization improvements of full core neutronics/thermal-hydraulics coupling between RMC and CTF[J].Nuclear Engineering and Design，2018，332：88-98.