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

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

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

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

中圖分類號(hào)：TL351.1? ? ? ? ? ? ? ? ? ? ? ?文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)：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

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

1? 耦合方法

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

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

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

1.2 耦合時(shí)間步算法

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

1.3 耦合收斂算法

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

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

2? 計(jì)算模型

2.1 耦合方案

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

2.2 數(shù)值模型

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

3? 數(shù)值分析

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

3.1 Picard方法

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

3.2 時(shí)間步長敏感性分析

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

3.3 改進(jìn)Picard方法

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

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

4? 結(jié)語

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

參考文獻(xiàn)

[1] 周夏峰，李富.基于節(jié)塊展開法的Jacobian-Free Newton Krylov聯(lián)立求解物理-熱工耦合問題[J].物理學(xué)報(bào)，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é)塊中子動(dòng)力學(xué)方程組的數(shù)值解法及物理與熱工-水力耦合瞬態(tài)過程的數(shù)值計(jì)算的研究[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.