论文部分内容阅读
DOI：10.16660/j.cnki.1674-098x.2104-5640-5601
摘 要：本文基于格林函数节块法程序TNGFM与子通道程序COBRA-EN开发了三维瞬态核热耦合计算程序，利用压水堆弹棒算例对核热耦合数值稳定性进行了分析，并提出了改进Picard迭代方法。数值模拟结果表明，减小计算时间步长可提高耦合计算数值稳定性，但计算效率较低，在物理计算外迭代中引入近似反馈可有效提高耦合收敛性能并大大减少耦合计算迭代次数。
关键词：瞬态 核热耦合 数值稳定性 计算效率
中图分类号：TL351.1 文献标识码：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 Green’s 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
隨着高性能计算机技术的迅速发展，物理-热工耦合计算趋向于高保真技术发展，采用三维物理-热工耦合分析能更精确模拟瞬态或事故过程，提高瞬态工况下物理及热工参数分布计算的准确性。近年来，国内外针对物理-热工耦合问题开展了大量研究，并开发了多套耦合程序，这些程序多采用算符分裂或Picard迭代等方法求解，容易出现数值不稳定和收敛速度较慢等问题。本文基于格林函数节块法程序TNGFM与子通道程序COBRA-EN开发了瞬态核热耦合程序，对瞬态核热耦合数值稳定性进行了分析，并提出了改进迭代算法改善瞬态核热耦合收敛性能。
1 耦合方法
堆芯瞬态核热耦合需解决耦合实现方式、耦合时间步算法、耦合收敛算法等3个关键问题。
1.1 耦合实现方式
耦合实现方式从数值求解角度可分为强耦合和弱耦合。强耦合是指将物理-热工方程联立求解，如近几年发展的Jacobian-free Newton-Krylov（JFNK）[1]方法即可实现耦合变量的同步更新，但实现较为复杂。目前大多数核热耦合程序如COBAYA4/COBRA-TF[2]、RMC/CTF[3]仍以弱耦合为主，程序独立求解物理-热工方程，通过交换边界数据实现耦合，能充分利用现有程序，实现相对容易。
1.2 耦合时间步算法
耦合时间步算法主要包括物理与热工程序各自计算时间步长、耦合参数传递时间点的选取及耦合程序时间步推进格式。物理与热工程序时间步长的选取通常取决于程序自身的计算精度，较为简便的方法为各程序选取相同的时间步长，并在每个时间步内进行耦合参数的传递。耦合程序时间步推进格式通常分为显式与隐式耦合[4]，显式耦合是指在单个时间步内，物理与热工程序采用一次通过方式计算，不进行迭代，计算量较小但精度较低。隐式耦合则需要物理与热工程序在单个时间步内进行反复迭代直到满足相应收敛准则为止，计算量较大但能保证耦合边界条件一致且稳定性更好。 1.3 耦合收敛算法
传统的核热耦合收敛算法包括OS方法[5]、固定点迭代法（Picard方法）[6]，它们都是先由物理计算得到堆芯内功率分布，将数据传递给热工程序计算得到堆芯内温度、密度等热工参数，再由反馈计算模块更新物理计算所需截面后进行物理计算，如此反复迭代直至计算收敛。该算法由于流程简单较易实现，但也容易出现计算震荡乃至发散的现象，目前主要通过引入松弛因子以保证耦合计算的收敛[7]。
在实际核热耦合计算中，由于初始物理及热工参数均为猜测值，没有必要在每次物理计算后都进行详细的热工分析计算再进行反馈计算，尤其是瞬态计算时此种方式会大大降低计算效率，并可能导致计算不收敛。由此本文提出改进Picard方法，即先近似建立热工参数与截面的函数关系，在物理源迭代的外迭代中引入该近似反馈，在物理计算收敛后进行详细的热工分析计算得到温度、密度等分布，再更新截面参数进行核热耦合计算，如此反复迭代直至计算收敛。
2 计算模型
2.1 耦合方案
本文采用隐式耦合方式开发了三维瞬态核热耦合程序TNGFM-C，在每个时间步内，物理程序求解堆芯功率分布传递给热工分析程序计算得到燃料温度、水温度、水密度等热工参数，然后反馈计算模块根据最新热工参数更新截面传递给物理程序，计算流程如图1所示。
2.2 數值模型
本文选取某PWR寿期初三维弹棒算例进行瞬态核热耦合计算分析，图2为该堆芯控制棒组件布置示意图，计算初始条件为：热态满功率，控制棒组A、B、C、S提出225步，D棒组提出156步，考虑最不利的工况，即最大价值的控制棒组D中的L07弹出堆芯，在瞬态程序中设置该束控制棒弹出速度为1800cm/s，不考虑触发停堆，耦合计算时间步长为0.01s。由于TNGFM为粗网节块法程序，COBRA-EN采用全堆芯分析模型，因此在计算时径向划分均以组件为单位，轴向划分20层，考虑径向及轴向反射层，共划分为4248个物理节块，121个子通道。
3 数值分析
利用TNGFM-C耦合程序对瞬态弹棒工况进行了计算，分析了采用传统Picard方法时耦合计算的数值稳定性，对耦合计算时间步长敏感性进行了分析，并比较了改进Picard方法对耦合收敛性能及计算效率的影响。
3.1 Picard方法
图3为采用传统Picard方法及引入松弛因子时堆芯相对功率随耦合迭代次数的变化。从图中可看出，在单个时间步耦合迭代初期，由于堆芯功率急剧上升，在热工反馈作用下耦合计算存在明显震荡现象，在引入松弛因子后有所改善，但每个时间步内耦合计算仍需多次迭代才能收敛。同时随着时间步的推进，堆芯功率变化程度减缓，在单个时间步内耦合收敛所需迭代次数逐渐减少，但引入松弛因子并未明显减少耦合计算总迭代次数。
3.2 时间步长敏感性分析
在弹棒事故中，在单个时间步内由于引入较大正反应性导致堆芯功率剧烈变化，使得耦合计算出现明显震荡，因此本文分析比较了不同时间步长对核热耦合计算数值稳定性的影响。图4分别为时间步长为0.01s、0.005s、0.002s及0.001s时堆芯相对功率随耦合迭代次数的变化，从图中可看出，时间步长越小，耦合迭代计算过程中功率震荡幅度越小，这是由于时间步长越小则单个时间步内引入的反应性越小，耦合的反馈效应相对越弱，使得单个时间步内耦合收敛所需的迭代次数越少，但时间步长越小整体耦合计算所需时间越长，因此在实际计算中需综合考虑计算稳定性及效率，选择合适的时间步长。
3.3 改进Picard方法
由前文分析可知，采用传统Picard迭代法进行瞬态耦合计算时存在明显震荡现象，引入松弛因子后有所改善，但计算效率未有明显提升，采用较小时间步长虽能提升计算数值稳定性，但计算效率明显降低，因此本文提出改进Picard方法，即在物理计算源迭代的外迭代过程中引入近似反馈计算，物理计算收敛后进行热工及反馈计算。图5为采用改进Picard方法时堆芯相对功率随耦合迭代次数的变化，从图中可看出，由于在物理源迭代中完成了近似耦合计算，物理热工程序间耦合计算迭代次数大大减少，同时计算数值稳定性大大提高。
表1为分别采用传统与改进Picard迭代方法时迭代次数及计算时间对比，从表中可看出，采用传统Picard迭代方法时，物理计算迭代次数较少，但耦合迭代次数较多，采用改进Picard迭代方法时由于在物理计算外迭代中引入近似反馈，物理计算迭代次数较多，但耦合迭代次数大大减少，综合耦合计算时物理热工程序各自计算时间及数据传递时间，采用改进Picard迭代方法耦合计算效率相对较高。
4 结语
本文以格林函数节块法程序TNGFM与子通道程序COBRA-EN为基础开发了三维物理-热工瞬态耦合程序，利用压水堆寿期初弹棒算例对耦合计算数值稳定性进行了分析，结果表明，在瞬态核热耦合计算时，采用传统Picard迭代方法存在明显数值震荡现象，引入松弛因子后有所改善但计算效率无明显提升，同时减小时间步长可提高数值稳定性但计算效率较低，采用改进Picard迭代法能大幅减少耦合迭代次数，提高计算效率，且能改善耦合迭代过程中的数值震荡现象。
参考文献
[1] 周夏峰，李富.基于节块展开法的Jacobian-Free Newton Krylov联立求解物理-热工耦合问题[J].物理学报，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] 廖承奎.三维节块中子动力学方程组的数值解法及物理与热工-水力耦合瞬态过程的数值计算的研究[D].西安：西安交通大学，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.
摘 要：本文基于格林函数节块法程序TNGFM与子通道程序COBRA-EN开发了三维瞬态核热耦合计算程序，利用压水堆弹棒算例对核热耦合数值稳定性进行了分析，并提出了改进Picard迭代方法。数值模拟结果表明，减小计算时间步长可提高耦合计算数值稳定性，但计算效率较低，在物理计算外迭代中引入近似反馈可有效提高耦合收敛性能并大大减少耦合计算迭代次数。
关键词：瞬态 核热耦合 数值稳定性 计算效率
中图分类号：TL351.1 文献标识码：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 Green’s 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
隨着高性能计算机技术的迅速发展，物理-热工耦合计算趋向于高保真技术发展，采用三维物理-热工耦合分析能更精确模拟瞬态或事故过程，提高瞬态工况下物理及热工参数分布计算的准确性。近年来，国内外针对物理-热工耦合问题开展了大量研究，并开发了多套耦合程序，这些程序多采用算符分裂或Picard迭代等方法求解，容易出现数值不稳定和收敛速度较慢等问题。本文基于格林函数节块法程序TNGFM与子通道程序COBRA-EN开发了瞬态核热耦合程序，对瞬态核热耦合数值稳定性进行了分析，并提出了改进迭代算法改善瞬态核热耦合收敛性能。
1 耦合方法
堆芯瞬态核热耦合需解决耦合实现方式、耦合时间步算法、耦合收敛算法等3个关键问题。
1.1 耦合实现方式
耦合实现方式从数值求解角度可分为强耦合和弱耦合。强耦合是指将物理-热工方程联立求解，如近几年发展的Jacobian-free Newton-Krylov（JFNK）[1]方法即可实现耦合变量的同步更新，但实现较为复杂。目前大多数核热耦合程序如COBAYA4/COBRA-TF[2]、RMC/CTF[3]仍以弱耦合为主，程序独立求解物理-热工方程，通过交换边界数据实现耦合，能充分利用现有程序，实现相对容易。
1.2 耦合时间步算法
耦合时间步算法主要包括物理与热工程序各自计算时间步长、耦合参数传递时间点的选取及耦合程序时间步推进格式。物理与热工程序时间步长的选取通常取决于程序自身的计算精度，较为简便的方法为各程序选取相同的时间步长，并在每个时间步内进行耦合参数的传递。耦合程序时间步推进格式通常分为显式与隐式耦合[4]，显式耦合是指在单个时间步内，物理与热工程序采用一次通过方式计算，不进行迭代，计算量较小但精度较低。隐式耦合则需要物理与热工程序在单个时间步内进行反复迭代直到满足相应收敛准则为止，计算量较大但能保证耦合边界条件一致且稳定性更好。 1.3 耦合收敛算法
传统的核热耦合收敛算法包括OS方法[5]、固定点迭代法（Picard方法）[6]，它们都是先由物理计算得到堆芯内功率分布，将数据传递给热工程序计算得到堆芯内温度、密度等热工参数，再由反馈计算模块更新物理计算所需截面后进行物理计算，如此反复迭代直至计算收敛。该算法由于流程简单较易实现，但也容易出现计算震荡乃至发散的现象，目前主要通过引入松弛因子以保证耦合计算的收敛[7]。
在实际核热耦合计算中，由于初始物理及热工参数均为猜测值，没有必要在每次物理计算后都进行详细的热工分析计算再进行反馈计算，尤其是瞬态计算时此种方式会大大降低计算效率，并可能导致计算不收敛。由此本文提出改进Picard方法，即先近似建立热工参数与截面的函数关系，在物理源迭代的外迭代中引入该近似反馈，在物理计算收敛后进行详细的热工分析计算得到温度、密度等分布，再更新截面参数进行核热耦合计算，如此反复迭代直至计算收敛。
2 计算模型
2.1 耦合方案
本文采用隐式耦合方式开发了三维瞬态核热耦合程序TNGFM-C，在每个时间步内，物理程序求解堆芯功率分布传递给热工分析程序计算得到燃料温度、水温度、水密度等热工参数，然后反馈计算模块根据最新热工参数更新截面传递给物理程序，计算流程如图1所示。
2.2 數值模型
本文选取某PWR寿期初三维弹棒算例进行瞬态核热耦合计算分析，图2为该堆芯控制棒组件布置示意图，计算初始条件为：热态满功率，控制棒组A、B、C、S提出225步，D棒组提出156步，考虑最不利的工况，即最大价值的控制棒组D中的L07弹出堆芯，在瞬态程序中设置该束控制棒弹出速度为1800cm/s，不考虑触发停堆，耦合计算时间步长为0.01s。由于TNGFM为粗网节块法程序，COBRA-EN采用全堆芯分析模型，因此在计算时径向划分均以组件为单位，轴向划分20层，考虑径向及轴向反射层，共划分为4248个物理节块，121个子通道。
3 数值分析
利用TNGFM-C耦合程序对瞬态弹棒工况进行了计算，分析了采用传统Picard方法时耦合计算的数值稳定性，对耦合计算时间步长敏感性进行了分析，并比较了改进Picard方法对耦合收敛性能及计算效率的影响。
3.1 Picard方法
图3为采用传统Picard方法及引入松弛因子时堆芯相对功率随耦合迭代次数的变化。从图中可看出，在单个时间步耦合迭代初期，由于堆芯功率急剧上升，在热工反馈作用下耦合计算存在明显震荡现象，在引入松弛因子后有所改善，但每个时间步内耦合计算仍需多次迭代才能收敛。同时随着时间步的推进，堆芯功率变化程度减缓，在单个时间步内耦合收敛所需迭代次数逐渐减少，但引入松弛因子并未明显减少耦合计算总迭代次数。
3.2 时间步长敏感性分析
在弹棒事故中，在单个时间步内由于引入较大正反应性导致堆芯功率剧烈变化，使得耦合计算出现明显震荡，因此本文分析比较了不同时间步长对核热耦合计算数值稳定性的影响。图4分别为时间步长为0.01s、0.005s、0.002s及0.001s时堆芯相对功率随耦合迭代次数的变化，从图中可看出，时间步长越小，耦合迭代计算过程中功率震荡幅度越小，这是由于时间步长越小则单个时间步内引入的反应性越小，耦合的反馈效应相对越弱，使得单个时间步内耦合收敛所需的迭代次数越少，但时间步长越小整体耦合计算所需时间越长，因此在实际计算中需综合考虑计算稳定性及效率，选择合适的时间步长。
3.3 改进Picard方法
由前文分析可知，采用传统Picard迭代法进行瞬态耦合计算时存在明显震荡现象，引入松弛因子后有所改善，但计算效率未有明显提升，采用较小时间步长虽能提升计算数值稳定性，但计算效率明显降低，因此本文提出改进Picard方法，即在物理计算源迭代的外迭代过程中引入近似反馈计算，物理计算收敛后进行热工及反馈计算。图5为采用改进Picard方法时堆芯相对功率随耦合迭代次数的变化，从图中可看出，由于在物理源迭代中完成了近似耦合计算，物理热工程序间耦合计算迭代次数大大减少，同时计算数值稳定性大大提高。
表1为分别采用传统与改进Picard迭代方法时迭代次数及计算时间对比，从表中可看出，采用传统Picard迭代方法时，物理计算迭代次数较少，但耦合迭代次数较多，采用改进Picard迭代方法时由于在物理计算外迭代中引入近似反馈，物理计算迭代次数较多，但耦合迭代次数大大减少，综合耦合计算时物理热工程序各自计算时间及数据传递时间，采用改进Picard迭代方法耦合计算效率相对较高。
4 结语
本文以格林函数节块法程序TNGFM与子通道程序COBRA-EN为基础开发了三维物理-热工瞬态耦合程序，利用压水堆寿期初弹棒算例对耦合计算数值稳定性进行了分析，结果表明，在瞬态核热耦合计算时，采用传统Picard迭代方法存在明显数值震荡现象，引入松弛因子后有所改善但计算效率无明显提升，同时减小时间步长可提高数值稳定性但计算效率较低，采用改进Picard迭代法能大幅减少耦合迭代次数，提高计算效率，且能改善耦合迭代过程中的数值震荡现象。
参考文献
[1] 周夏峰，李富.基于节块展开法的Jacobian-Free Newton Krylov联立求解物理-热工耦合问题[J].物理学报，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] 廖承奎.三维节块中子动力学方程组的数值解法及物理与热工-水力耦合瞬态过程的数值计算的研究[D].西安：西安交通大学，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.