模块式高温气冷堆(high-temperature gas-cooled reactor,HTGR)是具有固有安全特性的第四代先进球床式核动力反应堆[1-2]。堆芯在严重事故工况下仅需依靠自然的热传递就能将堆芯余热有效载出。近年来,高温球床传热成为国内外学者的研究焦点[3]。
德国、南非和中国先后开展了大型球床高温传热实验[4-6],测量了不同温度下的球床平均有效导热系数,结果指出高温(超过700 ℃)下球床堆辐射换热发挥主导作用。早期颗粒辐射研究主要集中在辐射传输方程上,将研究对象简化为均匀的吸收、发射和散射介质,主要应用在大气、粉尘和气溶胶等稀相颗粒系统。针对稠密的高温球床,Gusarov[7]指出辐射传输方程预测结果显著低于球床辐射等效导热系数的实验值。
在微观层次研究中,离散单元法广泛应用在球床球流运动、堆积和搅混中[2, 8]。由于球床堆堆积结构复杂,颗粒间的辐射换热受到周围球的阻挡,直接计算辐射换热角系数极为烦琐。在计算流体力学-离散单元法(computational fluid dynamics-discrete element method, CFD-DEM) 仿真中常采用近距离截断模型[9-10],但难以描述远程颗粒对之间的辐射特征。Johnson等[11]提出的径向距离近似方案仍无法反映随机堆积下的球床局部热辐射特征。
机器学习的发展为解决高温球床辐射传热提供了新的研究视角,以梯度提升决策树(gradient boosting decision tree, GBDT)[12]算法为代表的人工智能方法能够高效捕获复杂传热问题的非线性特征[13]。基于以上背景,本研究利用机器学习中的GBDT算法,结合高温球床堆中颗粒间辐射规律,训练开发颗粒间辐射角系数智能回归预测模型。
1 机器学习模型针对典型的10 MW高温气冷实验堆(HTR-10)堆芯[14],利用离散单元法获取球床随机堆积结构,如图 1所示。在高温堆球床的颗粒尺度辐射模型中,任意2个燃料球间的净辐射换热量Qr,ij可表示为
$ Q_{\mathrm{r}, i j}=\varepsilon \sigma_{\mathrm{b}} A_i X_{i j}\left(T_j^4-T_i^4\right). $ | (1) |
![]() |
图 1 球床随机堆积结构示意图 |
其中:ε、σb和Ai分别表示颗粒表面发射率、黑体辐射常数和颗粒i的表面积,高温堆中颗粒发射率取0.8,球形燃料元件设计直径为60 mm;Ti和Tj分别表示颗粒i和j表面的绝对温度;Xij表示受到其他周围颗粒阻挡下的颗粒间角系数。
此外,基于连续介质动力学假设,大型高温球床宏观辐射有效导热系数为
$ k_{\mathrm{r}}=\varepsilon \sigma T^3 \cdot \frac{2}{3} \rho \bar{A} \frac{1}{M} \sum\limits_i \sum\limits_j X_{i j} r_{i j}^2 . $ | (2) |
其中:ρ、rij和M分别表示球床中燃料球平均数密度、颗粒间球心距和颗粒总数量;
可以看出,在球床辐射换热模型中无论是微观颗粒尺度计算还是宏观分析,仅Xij需要根据相对堆积结构计算,存在辐射换热关系的颗粒不仅包括周围相互接触的颗粒,还包括距离较远的颗粒。通过计算发现,需要考虑周围3~4倍直径范围内约110个颗粒的辐射作用。颗粒间角系数计算极为复杂,利用传统数值方法在中央处理器(central processing unit,CPU)上计算1×104个颗粒间角系数矩阵约需1.6×107 s,难以用于分析超过42×104个燃料球的大型高温堆球床。此时燃料球间热辐射角系数的机器学习计算模型可表示为
$ X_n=B_{0, n}+\sum\limits_{k=1}^N f_{k, n} . $ | (3) |
其中:n表示单次计算中颗粒球的数量,即n=2表示需要计算无其他颗粒阻挡的2个颗粒间的角系数,n=5表示需要考虑周围3个颗粒对2个颗粒间角系数的影响;B0,n表示角系数半解析基础计算模型;
基础计算模型主要基于角系数的解析理论。利用角系数的定义,目前仅针对无阻挡的2个等径颗粒球情况存在显式解析解,具体可表示为
$ \begin{gathered} X_0\left(r_{12}\right)=\frac{1}{2}- \\ \frac{1}{6 \pi}\left[\left(4-r_{12}^2\right) K\left(\frac{2}{r_{12}}\right)+\left(r_{12}^2+4\right) E\left(\frac{2}{r_{12}}\right)\right] . \end{gathered} $ | (4) |
其中:X0表示颗粒间的角系数;
$ \begin{gathered} X_0\left(r_{12}\right)=\frac{1}{2}- \\ \frac{1}{3 \pi r_{12}}\left[\left(r_{12}^2-4\right) K\left(\frac{r_{12}}{2}\right)+\left(r_{12}^2+4\right) E\left(\frac{r_{12}}{2}\right)\right] . \end{gathered} $ | (5) |
虽然解析模型存在非初等函数,但椭圆积分函数是许多编程语言的内置函数。此角系数的解析结果可用于验证数值模型的精度。计算2个颗粒球间的角系数只需考虑基础计算模型,无需额外的机器学习模型预测,即X2=B0,2=X0(r12)。
针对3个颗粒球情况(n=3),角系数的主要影响因素包括球心距r12和外部颗粒球的阻挡,如图 2所示。角系数X12的半解析基础计算模型为
$ B_{0, 3}=X_0\left(r_{12}\right) \cdot g\left(s_3\right) \text {. } $ | (6) |
![]() |
图 2 球床中多个颗粒相对位置阻挡示意图 |
其中:g(s)为颗粒阻挡函数;s3表示外部第3个颗粒中心到颗粒1和2连线的无量纲距离,即s3= x/d,d表示颗粒直径,x表示距离矢量。同样,s4和s5表示外部第4和第5个颗粒到连线的无量纲距离。当s3≥1时,外部颗粒距离较远,不参与颗粒1和2之间的辐射换热,即g=1;当s3=0时,外部颗粒刚好在颗粒1和2中间,完全阻挡了两者之间的热辐射换热,即g=0。结合光线追踪数值计算结果(见图 3),阻挡函数可近似拟合为
$ g(s)=5.1358 s^5-12.023 s^4+6.6396 s^3+1.248 s^2 . $ | (7) |
![]() |
图 3 角系数基础计算模型阻挡函数 |
该函数描述了外部颗粒对辐射换热角系数的阻挡作用,在计算中无需考虑较远的周围颗粒的影响,具体评估详见第2节。若周围存在多个颗粒,阻挡作用更为复杂,难以用简单表达式准确计算。这里近似取乘积作用,对于n≥4,多球阻挡的情况(见图 2),角系数X12的基础计算模型调整为
$ B_{0, n}=X_0\left(r_{12}\right) \cdot\left[\prod\limits_{m=3}^n g\left(s_m\right)\right]^{1.4} . $ | (8) |
其中:sm表示外部每一个颗粒(0≤sm≤1)到颗粒1和2连线的距离,m为外部颗粒编号,经排序后有sm < sm+1。
高温球床辐射换热的基础计算模型可以捕捉角系数的主要成分,包括核心的沿着球心距离衰减的规律和周围颗粒的平均阻挡作用。由于只采用了解析表达式,计算速度快,在精度需求不高的情况下可直接采用基础计算模型获取球床颗粒间角系数的近似值。GBDT模型用于修正角系数基础模型并进一步提高计算精度,其模型结构为
$ F_N=\sum\limits_{k=1}^N f_{k, n}=f_{N, n}+\sum\limits_{k=1}^{N-1} f_{k, n}=f_{N, n}+F_{N-1} . $ | (9) |
计算中逐一添加回归树,而第N棵树用于预测前N-1棵树预测值与真实结果的残差,保证精度逐渐提升,避免模型退化问题。模型中的输入参数是各颗粒的几何位置坐标,输出量是角系数真实值与基础计算模型预测结果的差值,即
$ C_m=X_n\left(\boldsymbol{x}_m\right)-B_{0, n}\left(\boldsymbol{x}_m\right) . $ | (10) |
为便于统一分析和对比,对Cm进行标准化预处理,即
$Y_m=\frac{C_m-\bar{C}_m}{\sigma}. $ | (11) |
其中,
$\mathrm{MSE}=\frac{1}{M} \sum\limits_{m=1}^M\left(Y_m-\hat{Y}_m\right)^2 $ | (12) |
此时GBDT模型的回归系数Rf2可进一步简化为
$ R_f^2=1-\frac{\sum\left(Y_m-\hat{Y}_m\right)^2}{\sum\left(Y_m-\bar{Y}\right)^2}=1-\mathrm{MSE} . $ | (13) |
本研究利用开源机器学习平台LightGBM[12]开发GBDT模型,而训练和优化中采用自动机器学习平台FLAML[15]自动搜索最佳模型训练参数。
2 结果与讨论 2.1 球床角系数数据集本研究基于模块式高温气冷堆示范项目HTR-PM堆芯[16]球床离散堆积结构(见图 4)获取用于训练机器学习模型的角系数数据集。球床直径3 m,随机填充约42×104个直径为60 mm的燃料球。球床局部燃料球可与大量颗粒球产生辐射换热,任意颗粒对间的角系数均会受到多个周围颗粒阻挡, 故在数据集处理中删除周围不参与颗粒间辐射换热、sm>1的颗粒球,进而建立含有不同颗粒数量的子数据集。
![]() |
图 4 HTR-PM堆芯球床 |
对于机器学习模型的输出值,利用光线追踪算法[17],在英伟达GeForce RTX 3090显卡平台上并行计算每种工况下的热辐射角系数。计算中光线数量均固定为1.0×108,在无阻挡的2个颗粒间角系数数值解与解析解的最大绝对偏差为7.8×10-6,如图 5所示。
![]() |
图 5 HTR-PM球床中2个颗粒间角系数 |
由于高温球床中燃料颗粒密集堆积,周围颗粒对角系数阻挡显著。对比数值结果发现,距离大于8倍颗粒半径(r12≥8R)的颗粒对间辐射换热极为微弱(ΔX < 10-4),角系数近似为0。在颗粒数量超过12时,最外部颗粒对角系数的影响极小(ΔX < 10-4),故实际计算中可近似忽略n≥12的外部颗粒影响。综合以上规律,结合HTR-PM球床离散堆积特性,分别建立了n=3,4,…,11的多球工况,且r12 < 8R的9个角系数子数据集。每个子数据集的输入维度固定,输出维度均为1,HTR-PM球床共1.66×107条角系数工况,具体信息如表 1所示。
n | 输入维度 | 数据量 |
3 | 9 | 14 156 |
4 | 12 | 99 349 |
5 | 15 | 499 567 |
6 | 18 | 1 388 908 |
7 | 21 | 1 691 701 |
8 | 24 | 1 777 568 |
9 | 27 | 2 555 826 |
10 | 30 | 3 755 957 |
11 | 33 | 4 864 188 |
2.2 模型复杂度与回归系数
大数据集中的80%数据用于训练模型,剩余20%用于验证。为避免梯度提升决策树模型过拟合,需要合理限制模型的复杂度。梯度提升决策树模型复杂度可定义为
$ C=N \cdot L. $ | (14) |
其中: L表示每棵树叶子数量,C表示模型复杂度。在n=3,4,…,11时分别讨论FLAML计算精度,选取复杂度低且精度高的训练模型。FLAML获取的模型复杂度如图 6所示,在n=3时不同参数组合共554个模型,对比发现在回归树数量为22棵,每棵树24个叶子节点时,MSE为0.033 5。此时模型复杂度为528,是兼顾模型精度和复杂度的最优模型,继续增加复杂度时模型精度提升不明显。随着颗粒数量增加,在n=7时,最优模型中回归树数量增加到1 929棵,每棵树10个叶子节点,MSE降低到0.007 3。在最优的模型中需要更多棵回归树以达到更高的精度,在n=11时回归树数量达到8 549,叶子节点数仅为7,MSE为0.023 8。不同情况下的模型参数与回归系数如表 2所示,其中RB2表示基础计算模型的回归系数。
![]() |
图 6 梯度提升决策树模型的复杂度 |
n | RB2 | N | L | Rf2 |
3 | 0.999 | 22 | 24 | 0.971 |
4 | 0.998 | 103 | 18 | 0.990 |
5 | 0.998 | 628 | 15 | 0.992 |
6 | 0.997 | 1 197 | 13 | 0.992 |
7 | 0.997 | 1 929 | 10 | 0.992 |
8 | 0.994 | 4 165 | 8 | 0.994 |
9 | 0.981 | 2 425 | 15 | 0.991 |
10 | 0.954 | 5 578 | 10 | 0.987 |
11 | 0.914 | 8 549 | 7 | 0.978 |
利用基础计算模型直接计算辐射角系数,角系数大数据集上的回归系数在颗粒数目较少(n < 9)时精度较高,回归系数大于0.994,而在n≥9时精度有所下降。最优的梯度提升决策树模型用于预测角系数真实值与基础模型的差值,其预测回归系数Rf2均在0.978以上。综合基础预测模型和机器学习模型,角系数预测回归系数在n=3,4,…,11时均大于0.999,同等硬件(英特尔酷睿i9-10900K)上实现了平均13.5×104倍加速。
3 结论本文针对高温气冷堆球床辐射换热过程开展机器学习模型研究,基础计算模型和GBDT模型用于高效预测考虑周围阻挡条件下的颗粒间角系数。
基础计算模型基于角系数的显式解析理论,合理描述了角系数随球心距的衰减和周围颗粒的平均阻挡作用,可快速计算球床堆中辐射角系数的核心主导部分。
利用HTR-PM球床随机堆积结构,结合光线追踪技术建立了高温堆球床高精度角系数大数据集,共包含1.66×107条角系数工况,最大绝对误差约7.8×10-6。角系数大数据集很好地覆盖了各种常规局部结构。
GBDT模型用于优化并高效提升角系数预测精度,训练后的模型能够反映颗粒间辐射角系数计算规则。通过限制模型的复杂度避免了机器学习模型过拟合现象。综合机器学习模型角系数回归系数超过0.999,实现了13.5×104倍加速效果,可用于在颗粒尺度上分析大规模高温球床辐射换热规律。
[1] |
张作义, 吴宗鑫, 王大中, 等. 我国高温气冷堆发展战略研究[J]. 中国工程科学, 2019, 21(1): 12-19. ANG Z Y, WU Z X, WANG D Z, et al. Development strategy of high temperature gas cooled reactor in China[J]. Strategic Study of CAE, 2019, 21(1): 12-19. (in Chinese) |
[2] |
姜胜耀, 桂南, 杨星团, 等. 球床式高温气冷堆球流及球床辐射的理论与实验研究进展[J]. 原子能科学技术, 2019, 53(10): 1918-1929. JIANG S Y, GUI N, YANG X T, et al. Theoretical and experimental research progress on pebble flow and effective thermal conductivity in pebble-type HTGR[J]. Atomic Energy Science and Technology, 2019, 53(10): 1918-1929. (in Chinese) |
[3] |
PENG Z B, DOROODCHI E, MOGHTADERI B. Heat transfer modelling in discrete element method (DEM)-based simulations of thermal processes: Theory and model development[J]. Progress in Energy and Combustion Science, 2020, 79: 100847. DOI:10.1016/j.pecs.2020.100847 |
[4] |
ROUSSEAU P G, DU TOIT C G, VAN ANTWERPEN W, et al. Separate effects tests to determine the effective thermal conductivity in the PBMR HTTU test facility[J]. Nuclear Engineering and Design, 2014, 271: 444-458. DOI:10.1016/j.nucengdes.2013.12.015 |
[5] |
WU Y Y, REN C, YANG X T, et al. Repeatable experimental measurements of effective thermal diffusivity and conductivity of pebble bed under vacuum and helium conditions[J]. International Journal of Heat and Mass Transfer, 2019, 141: 204-216. DOI:10.1016/j.ijheatmasstransfer.2019.06.071 |
[6] |
任成, 杨星团, 李聪新, 等. 高温堆球床等效导热系数测量实验加热系统设计及可行性验证[J]. 清华大学学报(自然科学版), 2014, 54(8): 1068-1072. REN C, YANG X T, LI C X, et al. Heating system design and validation for the pebble bed effective thermal conductivity experiment in a high temperature gas-cooled reactor[J]. Journal of Tsinghua University (Science and Technology), 2014, 54(8): 1068-1072. (in Chinese) |
[7] |
GUSAROV A V. Statistical approach to radiative transfer in the heterogeneous media of thin-wall morphology-Ⅱ: Applications[J]. Journal of Heat Transfer, 2019, 141(1): 012701. DOI:10.1115/1.4040958 |
[8] |
GUI N, HUANG X L, YANG X T, et al. HTR-PM-based 3D pebble flow simulation on the effects of base angle, recirculation mode and coefficient of friction[J]. Annals of Nuclear Energy, 2020, 143: 107442. DOI:10.1016/j.anucene.2020.107442 |
[9] |
MEHRABIAN R, SHIEHNEJADHESAR A, SCHARLER R, et al. Multi-physics modelling of packed bed biomass combustion[J]. Fuel, 2014, 122: 164-178. DOI:10.1016/j.fuel.2014.01.027 |
[10] |
ZHOU Z Y, YU A B, ZULLI P. Particle scale study of heat transfer in packed and bubbling fluidized beds[J]. AIChE Journal, 2009, 55(4): 868-884. DOI:10.1002/aic.11823 |
[11] |
JOHNSON E F, TARI I, BAKER D. Radiative heat transfer in the discrete element method using distance based approximations[J]. Powder Technology, 2021, 380: 164-182. |
[12] |
KE G L, MENG Q, FINLEY T, et al. LightGBM: A highly efficient gradient boosting decision tree[C]// Proceedings of the 31st International Conference on Neural Information Processing Systems. Long Beach, USA: ACM, 2017: 3149-3157.
|
[13] |
HUGHES M T, KINI G, GARIMELLA S. Status, challenges, and potential for machine learning in understanding and applying heat transfer phenomena[J]. Journal of Heat Transfer, 2021, 143(12): 120802. |
[14] |
GAO Z Y, SHI L. Thermal hydraulic calculation of the HTR-10 for the initial and equilibrium core[J]. Nuclear Engineering and Design, 2002, 218(1-3): 51-64. |
[15] |
WANG C, WU Q, WEIMER M, et al. FLAML: A fast and lightweight AutoML library[C]// Proceedings of 4th Machine Learning and Systems Conference. San Jose, USA, 2021: 434-447.
|
[16] |
ZHANG Z Y, DONG Y J, LI F, et al. The Shandong Shidao Bay 200 MWe high-temperature gas-cooled reactor pebble-bed module (HTR-PM) demonstration power plant: An engineering and technological innovation[J]. Engineering, 2016, 2(1): 112-118. |
[17] |
HOWELL J R, DAUN K J. The past and future of the Monte Carlo method in thermal radiation transfer[J]. Journal of Heat Transfer, 2021, 143(10): 100801. |