A Reduced Order Model Based on ANN-POD Algorithm for Steady-State Neutronics and Thermal-Hydraulics Coupling Problem

. Te neutronics and thermal-hydraulics (N/TH) coupling behavior analysis is a key issue for nuclear power plant design and safety analysis. Due to the high-dimensional partial diferential equations (PDEs) derived from the N/TH system, it is usually time consuming to solve such a large-scale nonlinear equation by the traditional numerical solution method of PDEs. To solve this problem, this work develops a reduced order model based on the proper orthogonal decomposition (POD) and artifcial neural networks (ANNs) to simulate the N/TH coupling system. In detail, the POD method is used to extract the POD modes and corresponding coefcients from a set of full-order model results under diferent boundary conditions. Ten, the backpropagation neural network (BPNN) is utilized to map the relationship between the boundary conditions and POD coefcients. Terefore, the physical felds under the new boundary conditions could be calculated by the predicated POD coefcients from ANN and POD modes from snapshot. In order to assess the performance of an ANN-POD-based reduced order method, a simplifed pressurized water reactor model under diferent inlet coolant temperatures and inlet coolant velocities is utilized. Te results show that the new reduced order model can accurately predict the distribution of the physical felds, as well as the efective multiplication factor in the N/TH coupling nuclear system, whose relative errors are within 1%.


Introduction
Nuclear reactor is a complex multiphysics system, resulting from the coupling between neutronics, fuid dynamics, heat transfer, and mechanics. Te neutronics and thermalhydraulics coupling system has attracted signifcant research attention due to its great infuence on the reactor power distribution and nuclear safety. Several coupling algorithms have been developed to solve this large-scale neutronics/thermal-hydraulics nonlinear system. Picard iteration is the most common method since it could fully reuse the existing codes for individual physics [1][2][3]. Specifcally, in Picard iteration, the coupling boundary information is exchanged between diferent physical felds to consider the coupling efect, while each physical feld is still solved by the original existing code. However, the outer iterations are required among all the physical felds whose convergence rate is only linear. Recently, several attempts have been made to accelerate the Picard iteration, such as the residual balance method [4] and Anderson acceleration [5]. Another widely used multiphysics solver is Jacobian-free Newton-Krylov method (JFNK), where all the physical felds are solved together to achieve the 2nd order convergence rate [6][7][8][9]. However, both the Picard iteration and JFNK method have to solve a large-scale nonlinear algebraic equation system. When simulating the same physical model under diferent boundary conditions, this large-scale system should be solved repeatedly, usually leading to expensively computational cost.
In order to solve this problem, the method of reduced order models (ROMs) based on the proper orthogonal decomposition (POD) has emerged as a promising strategy to efciently simulate the N/TH coupling system, as well as other issues in nuclear engineering community, such as the eigenvalue problem [10,11], the neutron transport problem [12], and the multigroup neutron difusion problem [13]. Vergari and German use the POD-based ROM method to solve the N/TH coupling problem, where the full-order model (FOM) is projected onto the POD basis to produce the low-dimensional system [14,15]. Here, by the projection operator, the partial diferential equations (PDEs) about physical felds have been transformed to ordinary diferential equations (ODEs) about the corresponding coefcient vectors which still have to be solved by numerical methods, called as the POD-Galerkinbased ROM method. Several attempts have been made to avoid solving the equations by using the artifcial neural network (ANN) and have been applied in nuclear reactor problems including the abnormal event diagnosis [16,17] and system state prediction [18,19]. Although both the POD-Galerkin method and ANN method have been successfully applied in nuclear reactor felds, there are relatively few reports of the POD-ANN method [20,21] which is fully data driven and fully equation solver free. Furthermore, amongst the studies reviewed, research about the application of the POD-ANN method in the N/TH coupling problem has not been made yet.
Terefore, in this paper, a nonintrusive POD-ANN method for the N/TH coupling problem has been developed. Te POD method is employed to capture the characteristic spatial structure and modes of physical felds and eigenvalues under diferent boundary conditions by the snapshot method. Ten, the backpropagation neutral network is selected to map the relationship between the boundary conditions and the POD coefcients. Finally, the trained ANN model is utilized to predict the POD coefcients with new boundary conditions, and then the distributions of coupled multiphysics felds could be reconstructed.
Te reminder of the paper is organized as follows: Section 2 provides a brief introduction of the coupling model and the POD-ANN-based ROM method used in this study. Section 3 presents the application of the proposed method to a simplifed N/TH coupling model. Finally, the conclusion of this work is summarized in Section 4.

Conservation Equations.
A simplifed pressurized water reactor model in a simplfed cylinder r-z two-dimensional water reactor model has been considered in this work. Te detailed information about this model could be found in [22]. Since the objective of this work is to develop a method for the fast-running N/TH coupling problem in the nuclear reactor, a simplifed neutron and thermal-hydraulics coupling model is selected here for the purpose of easy to study. Please note that the coupling relationship and values of the parameters are considered as practical as possible. Te results of this work could have the potential to indicate the performance of the POD-ANN method for the complicated and practical neutronics and thermal-hydraulics issue. From the mathematical perspective, the neutronics and thermal-hydraulics system can be described by fve nonlinear partial diferential equations.
Neutron difusion equation with constraint equation: where P is the fxed total reactor power and Ω is the computational domain. Solid porous medium energy conservation equation: Coolant continuity equation: Coolant momentum conservation equation: Coolant energy conservation equation: Te pressure Poisson equation used in this work is derived from the mass conservation equation (4) and the momentum equation (5): Te neutron cross sections in equation (1) are the functions depending on the fuel temperature and the coolant temperature. Te heat conduction of the solid fuel in equation (2) depends on the fuel temperature itself. Te meaning of each term is listed in Table 1. Te expression of the parameters can be seen in [22].

Numerical Discretization Method.
Te fnite volume method is utilized in this work to discretize the nonlinear PDEs. Particularly, the convection term in equation (6) is discretized with the upwind diference scheme, and the central diference scheme is utilized for difusion terms in these equations. Besides, the staggered grid mesh is employed to solve the coolant fow feld distribution and coolant temperature feld distribution. After discretization, the matrix form of the nonlinear equations is expressed by the following equation [22].
In equation (8), the term Y is the unknowns, A is the coefcient matrix, and the submatrix blocks of which are dependent on the unknowns Y. Te information about the boundary conditions and the source term is included in Q, which is not a constant vector obviously. Te Picard iteration is utilized to solve N/TH equation system (8), where the coupling information, such as the fssion power, fuel temperature, and coolant temperature, is exchanged among diferent physical felds. Specifcally, the efective multiplication factor k eff in the neutron difusion equation is solved numerically with the power iteration method.

POD-ANN-Based ROM Method
Te ROM method aims to replace the original full-order problem with an alternation of the lower order but accurate model. Te POD method is a widely used ROM method due to its broad applicability to linear and nonlinear systems. In this paper, a full-order N/TH coupling model described in Section 2.1.2 is considered. Te coupling equation system can be expressed by a function f(x, δ): R d ⟶ R N , where x is the domain of the space and N is the size of x, which depends on the number of the discretization points of the partial diferential equations. In this work, the magnitude for N is 10 3 . δ is the parameters which equal to boundary conditions in this work. Te basic mechanic of POD is to fnd a set of deterministic function u i (x) called POD basis which can best approximate the concerning function f(x, δ) in the least square sense. It is the optimization problem: where k i (δ) is the POD coefcient and ||•|| is the standard Euclidean.
Te Snapshot method is one way to achieve this purpose. Let S be a matrix consisting of the snapshots of the solution. In this work, the snapshots y i � f(x, δ i ) ∈ R N×L , i � 1 · · · L are collected from the numerical simulation of the full-order model constructed in Section 2.1.1 with a series of diferent boundary conditions. L is the number of the snapshots, whose value is usually much smaller than that of N.
We apply the singular value decomposition to S, and the basic vectors are the left singular vectors of S and the frst r eigenvector U � u 1 (x), u 2 (x) · · · u r (x) are chosen to be the POD basis vectors. Te number of the POD basis r is determined by the predefned energy ration c: With the POD basis, the original function can be reconstructed by

ANN Backpropagation Neutral Network.
Te backpropagation neutral network (BPNN) is a powerful artifcial neural network, which has been successfully applied in engineering communities. A typical BPNN model usually consists of input layers, hidden layers, and output layers. Here are the two main processes in the BPNN: forward propagation of information and error backpropagation. Except for the input layers, for each neuron in the output layers and hidden layers, the activation function transformed the outputs and the weight sum of the previous layer with a bias added. Te activation functions of the hidden layer and the output layer used in this work are "tanh" and "identity," respectively.
Te accuracy of the BPNN highly depends on the quantity of hidden layers and quantity of neurons at each hidden layer. Teoretically, a BPNN with one hidden layer can approach to any continuous function in the closed interval. Although increasing the number of hidden layers has the potential to improve accuracy, too many hidden layers complicate the network, increasing the training time and may lead to the overftting issue. At present, the trialand-error method is the most common and practical way to determine the number of hidden layers and number of nodes in the hidden layers. Te basic structure of one hidden layer BPNN used in this work is elucidated in Figure 1. Te four training processes of the BPNN are included as follows: frstly, the connection weights are initialized randomly within the given range; then, the predicted value is generated after a series of operations on the input data; after that, the connection weights are adjusted according to the error between the predicted value and the expected value; the Levenberg-Marquardt method is utilized in this work as the training function, which is a combination of the Grand method and Gauss-Newton method; such process continues until the error in norm 2 reaches the convergence criterion or reaches the maximum number of iterations.

POD-ANN Algorithm.
Based on the combination of POD with the ANN-BPNN, a computational model to solve the N/TH coupling problem is developed in this paper. Te implementation of this new approach is shown in Figure 2. A full-order N/TH coupling model, as described in Section 2.1, is calculated under diferent perturbed boundary conditions to establish a library of the required results. Ten, the POD method is utilized to extract the features of the physical feld based on the simulation results, where the POD basis and corresponding coefcients are obtained. After that, the BPNN is used to map the relationship between the boundary conditions and POD coefcients. Finally, the trained BPNN is utilized to predict the POD coefcients as outputs, with new boundary conditions as inputs. Te predicted POD modes and POD basis are then used to reconstruct the physical felds, and the deviation between the POD-ANN predicted values and numerical solution of N/TH coupling equations values is employed to evaluate the performance of the POD-ANN method.

Problem Description and Sampling Design.
A simplifed cylinder r-z two-dimensional pressurized water reactor model with a radius of 1.6 meters and a height of 4.4 meters is used to evaluate the performance of the proposed POD-ANN method. Te neutronics and thermalhydraulics coupling behavior is described by equations (1)- (7). Te vacuum boundary condition is defned at the outface boundaries for the neutron difusion equation, and the adiabatic boundary condition is defned at the boundary for solid energy conservation. Te coolant pressure at outlet boundary is set to be 15.5 MPa, and the inlet boundary is set to be fxed velocity and fxed temperature boundary condition. Te boundary conditions are listed in Table 2. Te reactor power is set to be 200 MW in equation (2). Figure 3 shows the numerical calculation result of physics distribution under the boundary condition of a coolant inlet velocity of 4 m/s and inlet temperature of 270℃.
A series of diferent inlet temperatures and speeds of the coolant is selected to generate the snapshots. Te inlet coolant temperature ranges from 280°C to 300°C, and the inlet coolant velocity is within the range of 4 m/s to 6 m/s. Te intervals of the coolant temperature and coolant velocity are 2℃ and 0.2 m/s, respectively. Terefore, there are 121 cases which are determined by the computational model described in Section 2.1.2, constructing a library of snapshots as training test. In order to verify the prediction accuracy of the model for new inputs, for the testing cases, the inlet coolant temperature ranges from 281℃ to 299℃ and the inlet coolant velocity is within the range of 4.1 m/s to 5.9 m/s.Te intervals of the coolant temperature and coolant velocity are 2℃ and 0.2 m/s for testing cases, whose total number is 100. Te distributions of the training and testing sample points are presented in Figure 4. Te locations of the testing sample are totally diferent from those of the training sample.

POD Mode Features.
Te POD method is performed on the snapshots generated by the 121 training sets to obtain the POD modes which include the important features of diferent physical felds in the N/TH coupling system. Te decay in the eigenvalue of the correlation matrices for each physical quantity is presented in Figure 5. It can be seen that  Science and Technology of Nuclear Installations the eigenvalues decay rapidly, which means only a few POD modes are sufcient to represent the characteristics of the physics felds. Te number of the POD basis is determined by the energy ration c in equation (11). In this work, c � 1 × 10 − 7 , indicating that the error between the numerical results and the reconstructed results is less than 0.1% in norm 2. For the neutron fux feld, only the frst three POD modes are required to satisfy the criterion of energy ration c. Te dimensions of the reduced spaces for the approximation of the diferent physical quantity are summarized in Table 3, and the distributions of the POD bases of each physical quantity are shown in Figures 6(a)-6(e). By comparing with the physics distribution in Figure 3, it can be seen that the frst or frst two POD basis could capture the key features of the snapshots. We take the frst three POD basis of the neutron fux as an example, Figure 6(a) shows that the frst order POD basis, providing the basic spatial distribution of the neutron fux. Te frst order POD basis is almost axial symmetry along the z axis, and the maximum is achieved at the center of the z axis and circle center in the r axis, satisfying the theoretical prediction of the neutron fux distribution. Te second order POD basis makes some adjustments among the z axis, which is approximately positive at the low part and negative at the upper part due to the temperature feedback efect. Compared with the frst two POD basis with one and two extreme points, respectively, the third POD basis has three extreme points which could further consider the local detailed distribution of the neutron fux feld. Te POD bases of other physical felds have the similar phenomenon. Please note that the efective multiplication factor, which is a key parameter for steadystate N/TH calculation, is also considered here.   Science and Technology of Nuclear Installations

Artifcial Neural Network Structure.
Te POD coefcients of high order change drastically with inlet velocity u 0 and inlet coolant temperature T 0 . As an example, Figure 7 shows the 3rd POD coefcient of fuel temperature under diferent inlet velocities. In this work, the BPNN is utilized to calculate the POD coefcients under diferent coolant inlet parameters. Please note that since the complicated relationship between POD coefcient and inlet conditions, the simple interpolation methods, such as the bilinear interpolation, are incompetent in such cases. Six BPNNs are constructed for the six physical quantities, respectively, where the inlet coolant velocity and temperature are the inputs and the sampling dataset of the POD coefcient are the outputs. Table 4 shows the settings of the training parameters of each BPNN. For all of the BPNNs, one hidden layer is used, and the optimized number of neurons at the hidden layer is determined for each BPNN after several attempts. Te Levenberg-Marquardt method is used as the training function due to its high efciency for small and medium-size problems. Te performance of the network is evaluated by the mean square error (MSE) between the target value and the predicted value of the outputs. Te term n is the number of the outputs, that is, the POD coefcients. Te information about residual convergence of each network is presented in Figure 8. For neutron fux, there are only a few tens of epochs to achieve the predefned stop criterion, while it requires about several hundreds of epochs for fuel solid temperature, coolant temperature, coolant velocity, coolant pressure, and efective multiplication factor, perhaps due to the relatively complicated  relationship between the POD coefcients and inputs (coolant velocity and temperature).
Te approximate felds are reconstructed based on the predicted POD coefcients from the ANN and POD bases for the training samples to evaluate the performance of the POD-ANN model. Ten, the POD-ANN model is further applied to the testing samples to assess the generalization ability of the POD-ANN model. Here, the deviation between reconstructed felds and the numerical results of PDEs in Section 2.1.2 is determined by two conventional evaluation metrics, the mean relative error (MRE), and the relative L 2 error (l 2 ), which is defned as where y pre i,j is the predicted results obtained from the POD-ANN model and y rel i,j is the numerical solution of PDEs. N is the total number of the discrete points, and n is the sample size which is 121 in the training cases and 100 in the testing cases.
Te maximum MRE and the mean MRE of each physical feld are listed in Table 5, and those of maximum l 2 error and the mean l 2 error are presented in Table 6. Te MRE and l 2 of each feld are under 1%, and the mean MRE of the efective multiplication factor reached 10 − 5 magnitude, which indicates the proposed model accurately approximate the steady-state N/TH coupling problem. Among the six physical quantities considered in this work, the MRE and l 2 error of coolant velocity in the testing cases are the largest, almost one order of magnitude larger than other physical quantities. Tis might be contributed by the relatively small change of coolant velocity between inlet and outlet. Moreover, the computational cost is also analyzed, as shown in Table 7, and the computational efciency of the POD-ANN method is nearly 80 times higher than that of the traditional Picard iteration method. Furthermore, Figure 9 shows the spatial distribution of the MRE of the training and testing cases of each physical feld except for the efective multiplication factor. As shown in Figures 9(a) and 9(b), the maximum MRE of the predicted neutron fux and fuel temperature is achieved near the coolant outlet, which might be contributed by the smaller values and relatively larger gradients of these two physical quantities. Please note that, due to the low level of the neutron fux in these regions, the relatively larger predicted errors have minor impact on the power level. Te MREs for the temperature and velocity and pressure of the coolant shown in Figures 9(c)-9(e) are also small, whose value are less than 0.1%.
In order to explore the prediction accuracy varied with input boundary conditions, the l 2 error distribution of each case, including both the training and testing samples, is shown in Figure 10. Te l 2 error of the training and testing   Science and Technology of Nuclear Installations        cases has similar distribution. Te l 2 error is relatively uniform and distributed within a small range, along with the input variables. Moreover, the performance of the POD-ANN model is also evaluated for the situations whose inlet coolant temperature and velocity are linear distribution. In detail, the average value of the inlet parameter, such as inlet coolant temperature and velocity, is kept as a constant, but the slope of the linear distribution is changed. Te slopes of u 0 and T 0 change from −1 to 0, and the sampling interval is 0.1. Terefore, the 121 samples are generated and randomly divided into training (70) and testing (51) sets. Te POD bases are then extracted, and the corresponding coefcients are calculated. Te distributions of the frst POD basis are shown in Figure 11. Te MRE and l 2 error of each physical feld in the testing cases under linear distribution inlet parameters are shown in Tables 8 and 9, respectively. Te max MRE and l 2 error of all the physical felds are within 1%, which indicates that the POD-ANN method also has the capability of accurate prediction for the complicated boundary conditions.

Conclusion
An efcient reduced order model based on the POD and ANN for the neutronics and thermal-hydraulics coupling problem is developed in this work. In order to evaluate its performance, a simplifed cylinder r-z two-dimensional pressurized water reactor model is considered with diferent inlet temperature and speed of coolant. In detail, the POD approach is employed to extract the characteristic of the physical felds by snapshots. Te spatial distribution of the POD bases shows that the frst or frst two POD bases capture the main features of the snapshots. Te higher order POD bases could further consider the local detailed distribution. Te BPNN is constructed with the boundary conditions as input and the POD coefcients as output. Te trained BPNNs are then utilized to predict the POD coefcients with new boundary conditions. To evaluate the accuracy of the POD-ANN-based reduced order model, the MRE and l 2 error between the predicted results from the reduced order model and the numerical solution of PDEs are analyzed. Results show the MRE and l 2 error of each feld are under 1% which indicates that the proposed approach can accurately predict the distribution of the physical felds, as well as the efective multiplication factor k eff , for this steady-state N/TH coupling issue. Te computational time of the POD-ANN method is nearly 80 times less than that of the traditional Picard iteration method. Furthermore, the relatively large MRE appears at the outlet, perhaps due to a small value or large gradient in these regions.
In this work, the BPNN is used to reduce the training time, but the main drawback of the BPNN is the overftting problem. In order to solve this issue, the deep learning neural network's plan to be used and applied to the practical multiphysical model in the future [23].

Data Availability
Te data used to support the fndings of this study are available from the corresponding author on request.

Conflicts of Interest
Te authors declare that they have no conficts of interest.