0.   引　言

在实际反应堆运行过程中，现有的反应堆监测系统针对堆内瞬态变化所探测到的数据往往是有限的^[1]，若要对反应堆瞬态进行镜像模拟，实现对未被探测器监测范围覆盖的关键数据监测，特别是对于堆内中子注量率的监测，需要基于这类有限的实测数据点和反应堆基本参数完成全堆关键参数的还原。

针对中子注量率的求解，传统反应堆物理数值计算一般采用确定论或蒙特卡罗方法来求解中子扩散方程，在面对较为复杂的几何边界条件和截面参数时，则需要做更多的假设，以实现对解空间做空间离散和时间离散后在离散结点上使用准静态方法等手段进行数值迭代计算。虽然大量实践证明，对于热中子反应堆，应用较为简单的双群理论就可以得到满足工程需要的结果^[2]，但是对于高保真的交互型反应堆模拟技术，例如反应堆数字孪生技术^[3-5]而言，满足其精度需求仍具有一定的挑战。

近年来，一种基于物理信息驱动的神经网络算法（PINN），其能够以少量的训练数据得到满足物理约束的模型^[6]的特性引起了反应堆物理分析领域学者的深度关注。刘东^[7]等使用PINN对多维中子扩散方程进行求解，并结合一系列加速收敛措施以及有效增殖系数（k_eff）并行搜索技术，对方程进行求解，取得了可观的精度结果。Wang^[8]等针对中子扩散方程的不光滑性引入了保守型PINN（cPINN），通过在不同材料性质子域中设置额外的守恒约束，实现了非均匀介质内中子扩散方程的求解。Yang^[9]等提出了一种数据驱动的PINN（DEPINN）用于求解中子扩散特征值问题，利用了少量来自实验的先验数据实现了神经网络训练精度和效率的提高。Prantikos^[10]等以接入反应堆数字孪生技术为目标，利用PINN求解点堆动力学下的中子扩散方程，明确了PINN相对传统数值求解方法的优势在于不需要耗费大量时间进行复杂网格划分，并提供在不规则问题和高维问题上更高效求解的途径。Elhareef^[11]等基于PINN对单能群和多能群扩散方程进行求解，其结果在两能群问题上的误差比在定源问题上大近15倍。而对于PINN，Kharazmi^[12]等基于子域Petrov-Galerkin方法进一步开发了具有更高精度的变分残差物理信息神经网络（hp-VPINN）。

现有研究普遍在以下两个方面存在一定的挑战，一是在精确度上，目前只做到了传统反应堆物理数值计算方法的同等量级甚至是低于传统方法的精度；二是在超参数设置上，普遍采用通用设置或是经验设置进行处理，并仅以敏感性分析作为设置值合理性的解释。针对这两方面挑战，本文基于hp-VPINN，结合k_eff智能搜索与反演等创新型关键技术，构建适用于求解中子扩散方程的hp-VPINN模型，并采用鲸鱼优化算法（WOA）^[13]的超参数自寻优技术实现网络结构优化，对单能单介质、单能多介质、双群多介质的中子扩散方程由简到繁进行求解，以验证各技术的有效性以及对各中子扩散方程的适用性。

1.   中子扩散方程物理模型

核反应堆物理分析中，在中子输运理论上作中子散射各向同性的近似假设消除角度自变量，并对连续的能量进行离散和平均化处理，再根据计算域中子守恒原则，可以得到多群的稳态中子扩散方程。

稳态情况下，更需要关注临界时的状态，而中子反应不一定都能达到临界状态，消除上述方程关于时间的非稳态项后还需引入$ {k}_{\text{eff}} $使反应达到临界状态，此时表示为如下的特征值方程：

式中，${\phi _g}$、$ {D}_{g} $、$ {\chi _g} $分别为第$ g $能群的中子注量率、扩散系数、能谱系数；${{{\varSigma}} _{\text{{r}},g}}$、${{{\varSigma}} _{\text{{a}},g}}$、$ {{{\varSigma}} _{\text{{f}},g}} $、${{{\varSigma}} _{\text{{s}},{g^\prime } \to g}}$分别表示第$ g $群中子的移出截面、吸收截面、裂变截面以及从第$ {g^\prime } $能群散射到第$ g $能群的散射截面；${v_{g^\prime }}$为相应的裂变产生的中子数；${S_g}$为外部第$ g $群中子源项。

当能量为单能且无外源时，瞬态和稳态中子扩散方程分别可写为：

式中$ ， L $为扩散长度；$V$为中子速率；${k_\infty }$为无限增殖系数；t为时间变量；${\boldsymbol{r}}$为空间变量。式（2）~式（3）可统一描述为如下形式的二阶微分方程：

式中，x为方程变量。

2.   求解中子扩散方程的hp-VPINN基本原理

PINN在传统神经网络的基础上，将表征物理信息的控制方程以及定解条件融入到损失函数中进行梯度下降，利用人工神经网络强大的非线性拟合能力，使其输出在求解域内满足控制方程并在边界上满足边界条件的偏微分方程解${{\phi _{{\rm{NN}}}}}\left( {\boldsymbol{x}} \right)$^[6]。hp-VPINN是在PINN的基础上结合子域Petrov-Galerkin方法^[14]发展而来，与一般的Galerkin方法相比，其允许试函数与权函数来自不同的函数空间，本研究采用人工神经网络${{\phi _{{\rm{NN}}}}}$作为试函数、式（5）所示勒让德多项式P(x)相邻两项之差作为基来构建K维权函数空间。权函数$w\left( x \right)$乘以式（4）并在整个求解域上进行积分，得到残差R。

基于子域Petrov-Galerkin方法的hp-VPINN还可对求解域进行分解实现求解域划分子域尺寸h的精细化，并能将变分残差投影到由高阶多项式组成的p维空间，提高求解精度，hp-VPINN还可以使用分部积分公式对二阶偏微分方程弱形式进行推导，降低对神经网络求导的阶数，使得神经网络训练更加灵活。具体步骤将在第4节算例中进行展示。

3.   基于hp-VPINN求解中子扩散方程的关键技术

3.1   多区域并行收敛技术

实际反应堆中往往采用分区布置，一般可分为增殖区与包裹区，不同区域采用不同材料，此时描述中子扩散过程的方程参数亦不相同，因此需要对整个堆芯进行空间分区域求解。通过基于子域Petrov-Galerkin方法的hp-VPINN框架可以很好地对求解域进行分解，并在每个子域上采用不同的损失函数进行训练，达到多区域并行收敛的效果。

3.2   区域划分加密技术

在求解中出现梯度较大的区域，对该区域进行区域划分加密，类似于有限元中的网格加密计算，从而提高神经网络训练结果精度。

3.3   $ {k}_{\text{eff}} $智能搜索以及反演技术

对于式（2）和式（3），只有当$ {k}_{\text{eff}} $等于确定数值时方程才有解，要精确求出$ {k}_{\text{eff}} $以及对应的特征函数$ \phi\left(\boldsymbol{\boldsymbol{r}}\right) $是非常困难的。若在hp-VPINN求解方程过程中采用源迭代方法循环往复地训练神经网络，会造成求解效率低的问题。为此基于开发的$ {k}_{\text{eff}} $智能搜索技术对瞬态中子扩散方程进行计算，若取一段时长$\tau $，在$t > \tau $时，裂变反应达到稳态，则有$ \phi\left(\boldsymbol{r},t\right) $对时间的偏导数等于零，将这一条件加入损失函数中，从而在网络训练过程中根据$ \phi_{\text{NN}}\left(\boldsymbol{r},t\right) $对时间的偏导数值不断对$ {k}_{\text{eff}} $进行调整，直到搜索到正确值。而对于稳态中子扩散方程，开发了$ {k}_{\text{eff}} $反演技术，将少量的从某些位置采样到的已知数据作为神经网络训练集，使hp-VPINN输出的解在这些区域上符合这些数据，从而通过方程反演出未知$ {k}_{\text{eff}} $，并获得方程的解。

3.4   方程组并行收敛技术

在对中子扩散方程进行能量平均化处理成G个能群后，会得到由G个偏微分方程构成的偏微分方程组，从而增大了hp-VPINN的求解难度。将每个偏微分方程单独作为hp-VPINN的损失项，从而实现神经网络求解多个耦合方程的并行收敛。

3.5   神经网络超参数优化技术

深度神经网络通常包含一系列超参数^[15]，例如神经网络的深度、宽度以及学习率等，而hp-VPINN额外增加了定解条件损失项相对于方程损失项的权重等超参数。对于超参数的设定，目前没有成熟可靠的理论可循，但其对神经网络的训练效率及精度有显著的影响^[16]，本文采用鲸鱼优化算法^[13]对超参数进行智能优化，其已被证明在求解精度和收敛速度方面均优于粒子群优化算法、差分进化算法、引力搜索算法和灰狼优化算法。

构建采用鲸鱼优化算法的神经网络超参数优化模型，流程如图1所示。通过使用超参数样本数据对神经网络进行训练，寻找超参数最优值。若训练结果中计算的适应度是最大值，则将其设为本次的最优结果，并与全局最优适应度进行对比，若该值大于全局最优适应度值，则进行替换。

图  1  基于鲸鱼优化算法的神经网络超参数优化流程

Figure  1.  Workflow Diagram of Neural Network Hyperparameter Optimization Based on Whale Optimization Algorithm

下载: 全尺寸图片 幻灯片

4.   方程求解与技术验证

4.1   一维平板方程中子扩散方程

考虑临界状态下一维平板方程：${\nabla ^2}\phi \left( x \right) + B_{\text{g}}^2 \cdot \phi \left( x \right) = 0$，$x \in \left[ { - 1,1} \right]$，$\phi \left( { \pm 1} \right) = 0$，临界状态下几何曲率$B_{\text{g}}^2 = {\pi ^2}$，并给定$\phi \left( 0 \right) = 0.5$，此时的解析解为$\phi \left( x \right) = 0.5{\text{cos}}\left( {\pi x} \right)$，构造损失函数l：

式中，$ {\tau }_\text{b} $、$ {\tau }_\text{i} $分别为边界条件损失项$ {l}_\text{b} $、初始条件损失项$ {l}_\text{i} $相对控制方程损失项$ {l}_\text{f} $的权重。

将原来的强形式设为形式1：$ {l}_{\text{f}}^{\left(0\right)} $，对$ {l}_\text{f} $采用分部积分并结合多项式的紧支撑特性可得到另外两种弱形式（形式2、3）：

3种形式使用的网络结构均为[1,80,80,80,80,1]、学习率为0.001，取100个高斯积分点、25阶勒让德多项式作为试函数，求解结果见图2，其均方误差分别为$8.71 \times {10^{ - 7}}$、$3.42 \times {10^{ - 6}}$、$2.16 \times {10^{ - 6}}$，其中形式1结果最准确，但采用形式2、3降低了对神经网络求导的阶数，提高了神经网络训练速度，经过实验验证，3种形式迭代5000步的时间分别为51.65、32.89、23.32 s，形式2和形式3迭代速度比形式1快两倍左右。从各形式损失函数的下降趋势可以发现，形式1与形式2的损失函数下降程度较大，而形式3在第350步左右即收敛不再下降。接下来的求解均采用形式1进行计算。

图  2  相对误差分布以及损失函数下降趋势

Figure  2.  Relative Error Distribution and Loss Function Decline Trend

下载: 全尺寸图片 幻灯片

采用域分解方法将求解域分为两段区域分别训练，分别在各相应子域采用勒让德多项式（其他区域为0）的分段函数作为权函数，因此构造出3个Case算例。求解结果如图3所示，图中黑色虚线为区域分界线（下同），在区域分解后，所使用的权函数能够在对应子域上使神经网络输出的解很好地符合物理方程解析解。该算例说明hp-VPINN的域分解方法能够很好地应用于中子扩散方程。

图  3  子域中使用勒让德多项式作为权函数时的求解结果

Figure  3.  Solution Results Using Legendre Polynomial as Weight Function in Subdomain

下载: 全尺寸图片 幻灯片

4.2   一维球坐标中子扩散方程

考虑一个半径为1（包括外推距离）的球形裸堆，若令$\phi \left( 0 \right) = \pi $，此时的中子扩散方程解析解为：

分别将求解域划分为[−1,−2/3]、[−2/3,2/3]、[2/3,1]三个区域，以及在梯度较大处进行加密划分：[−1,−2/3]、[−2/3,−1/3]、[−1/3,1/3]、[1/3,2/3]、[2/3,1]五个区域，计算所选参数以及网络结构同4.1节。从图4可以观察到，在梯度较大处进行局部加密后，相对误差出现了明显下降，加密前后的均方误差分别为$1.1 \times {10^{- 7}}$、$5.9 \times {10^{- 8}}$，加密前后分别训练5000步的耗时分别为61.64、103.55 s。该算例说明对区域划分进行加密可以有效提高hp-VPINN求解精度。

图  4  求解域分区加密前后的相对误差分布

Figure  4.  Relative Error Distribution before and after Mesh Refinement in the Solution Domain

下载: 全尺寸图片 幻灯片

4.3   非均匀介质内中子扩散方程

对于存在常源的一维单群双介质问题，一维空间被2种材料划分为3个区域，在$x = 0$处存在恒定点源向两侧发射中子。此时的中子扩散方程可写为：

式中，下标1、2为不同材料的区分；a为不同材料交界面处的位置。在$ 0 < \left| x \right| < \dfrac{{\left| a \right|}}{2} $处为材料1，其余处为材料2。

边界条件为：

（1）对称面中子源J(x)：

（2）当$x \to \infty $时，${\phi _2} = 0$，即：

（3）在介质交界面处中子注量率连续，即：

（4）在介质交界面处中子流密度连续，即：

当$x > 0$时，中子扩散方程的解析解为：

将求解域设置为$x \in \left[ {0,1} \right]$，并将区域划分为$\left[ {0,a/2} \right]$、$\left[ {a/2,1} \right]$。运用hp-VPINN技术求解该问题需要构造合适的损失函数，对于边界条件1，构造：

对于边界条件2，构造：

对于边界条件3，其为使用tanh函数作为激活函数的全连接神经网络，是一个连续函数，故无需设置专门的损失项。

对于边界条件4，构造：

控制方程损失项如4.1小节所述方法进行设置。取a=1 cm、D₁=D₂=1 cm、扩散长度${L_1} = 1\;{\text{cm}}$、${L_2}$分别为0.5、1.0、1.5 cm三个算例对中子扩散方程进行求解，结果及相对误差如图5所示，其均方误差分别为$2.12 \times {10^{ - 6}}$、$2.64 \times {10^{ - 8}}$、$6.76 \times {10^{ - 8}}$，说明所提出的方法对多介质中子扩散方程同样适用。网络结构及计算参数与4.1小节相同，迭代3000步运行时间为22.12 s。该算例说明hp-VPINN对非均匀介质下的中子扩散方程同样适用。

图  5  非均匀介质内中子扩散方程的求解结果与相对误差

Figure  5.  Solution Results and Relative Error of Neutron Diffusion Equation in Inhomogeneous Media

下载: 全尺寸图片 幻灯片

4.4   一维瞬态中子扩散方程

求解一维瞬态中子扩散方程，如式（2）所示，瞬态一维平板中子扩散方程解析解见文献[2]。

本文结合临界时物理特点开发的$ {k}_{\text{eff}} $智能搜索技术将$ {k}_{\text{eff}} $初始化为1，$t > \tau $时，有$\dfrac{{\partial \phi \left( {{\boldsymbol{r}},t} \right)}}{{\partial t}} = 0$，因此构造相应损失项如下：

式中，${n_r}$为$t = {t_j}$时空间采样点数量；n_t为时间采样点数量；$ t_j > \tau,j=1,2,3,\cdots ,n_t $。

设置2个算例对该方法进行验证，方程参数为$V = 2200\;{\text{m}}/{\text{s}}$、D=0.00211 m、L²=2103.7 m²、a=2 m，算例1、2的k_∞分别为1.0041、1.001，初始条件分别为cos(πx/a)−0.4[cos(2πx/a)+1]、cos(2πx/a)+1。求解结果如图6所示，均方误差均在${10^{ - 6}}$量级，通过智能搜索技术搜索到2个算例的$ {k}_{\text{eff}} $相对于参考值0.9995811、1.0035790的绝对误差分别为0.02pcm和0.01pcm（1pcm=10⁻⁵），迭代8000步用时为352.88 s。该算例说明基于hp-VPINN技术开发的k_eff智能搜索技术能够准确求出k_eff。

图  6  一维瞬态中子扩散方程求解结果

Figure  6.  Solution Results of 1D Transient Neutron Diffusion Equation

下载: 全尺寸图片 幻灯片

4.5   多群多区中子扩散方程

求解多群多区中子扩散方程，对反演技术及超参数优化技术（鲸鱼优化算法）进行验证。采用2D-TWIGL基准题^[17]的二维矩形堆芯结构，包含增殖区和包裹区，包括1个稳态计算和2个瞬态计算。堆芯形状为边长160 cm的正方形，包含3个区域。由于对称性，采用四分之一堆芯进行计算，堆芯布置及其参数如图7、表1所示。

图  7  2D-TWIGL基准题二维矩形堆芯结构示意图

Figure  7.  Schematic Diagram of 2D-TWIGL Benchmark Problem　　　　

下载: 全尺寸图片 幻灯片

表  1  分区截面参数

Table  1.  Zonal Cross-section Parameters

区域	g	${D_g}/{\text{cm}}$	$ \mathrm{\mathit{\Sigma}}_{\mathrm{a},g}/\text{c}\text{m}^{-1} $	$ v\mathit{\mathrm{\mathit{\Sigma}}}_{\mathrm{f},g}/\text{c}\text{m}^{-1} $	${{{\varSigma}} _{1 \to 2}}/{\text{c}}{{\text{m}}^{ - 1}}$
①,②	1	1.4	0.01	0.007	0.01
①,②	2	0.4	0.015	0.2
③	1	1.3	0.008	0.003	0.01
③	2	0.5	0.05	0.06

下载: 导出CSV 

| 显示表格

采用hp-VPINN结合参数反演技术对该基准题的稳态问题进行计算。首先采用KOMODO程序^[18]进行求解得到堆芯40×40中子注量率数据，再通过拉丁超立方采样选取全堆1%的数据，结合控制方程组以及边界条件等物理信息进行训练。对于复杂的方程组以及分区布置，为了使神经网络能够精确收敛、达到最好的性能，本文融合鲸鱼优化算法对其超参数进行寻优，设置人工鲸鱼数量为5、最大迭代次数为10次，并设置神经网络最大迭代次数为10000，训练结束条件为方程损失项小于10⁻⁴且$ {k}_{\text{eff}} $绝对误差小于20pcm，鲸鱼优化算法适应度设置为与KOMODO程序计算结果均方误差以及$ {k}_{\text{eff}} $均方误差成反比。图8展示了鲸鱼优化算法迭代过程中，每次迭代的最大归一化适应度变化，最终得到最优参数为学习率$4.9 \times {10^{ - 4}}$、数据损失项和边界条件损失项相对于方程损失项权重分别为8、9，神经网络结构为6层隐藏层[56,40,62,89,97,79]，耗时为1142.23 s。求解结果以及将对角线处归一化的求解结果分别如图9、10所示，与文献[19]参考结果进行对比，结果一致。反演$ {k}_{\text{eff}} $为0.913222，相对于参考值0.913214的绝对误差仅0.8pcm。

图  8  适应度随迭代次数的变化

Figure  8.  Variation of Fitness with Iteration Times

下载: 全尺寸图片 幻灯片

图  9  不同中子群求解结果

Figure  9.  Solution Results of Different Neutron Groups

下载: 全尺寸图片 幻灯片

图  10  归一化求解结果与参考解对角线分布

Figure  10.  Diagonal Distribution of Normalized Solution Results and Reference Solution

下载: 全尺寸图片 幻灯片

以上结果说明，基于hp-VPINN技术的中子扩散方程计算方法在全堆1%中子注量率数据基础上，通过鲸鱼优化算法实现了超参数自动优化，对全堆中子注量率分布实现了高精度求解。在反应堆物理分析上，该方法可为反应堆数字孪生技术提供一种较低实际数据依赖的全堆中子注量率分布还原手段，使依据现有反应堆内中子注量率实测数据反演反应堆实际运行过程中无法直接测量的全堆芯相关数据成为可能。

5.   结　论

本文基于hp-VPINN对中子扩散方程进行求解，并结合鲸鱼优化算法实现了神经网络超参数的自动优化。通过多个算例进行验证，证明了各种技术的有效性以及对各种中子扩散方程的适用性，并得出以下结论：

（1）hp-VPINN模型在正向求解单能单介质、单能多介质下的中子扩散方程时仅使用物理信息驱动即可达到良好的精度。

（2）$ {k}_{\text{eff}} $智能搜索技术可以高效代替源迭代搜索到方程的$ {k}_{\text{eff}} $并求解方程。

（3）反向求解多群多介质的中子扩散方程时，hp-VPINN模型采用鲸鱼优化算法搜索最佳超参数促进了神经网络的收敛，实现了仅采用全堆芯1%的数据即可精准反演${k_{\text{eff}}}$并重建全域中子注量率分布，反演结果绝对误差小于10⁻⁵。

本文提出方法对于先进反应堆模拟技术，例如反应堆数字孪生技术，提供了一条少量训练数据输入且较高精度输出的中子注量率求解可靠途径。