Abstract

The 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 differential 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 artificial neural networks (ANNs) to simulate the N/TH coupling system. In detail, the POD method is used to extract the POD modes and corresponding coefficients from a set of full-order model results under different boundary conditions. Then, the backpropagation neural network (BPNN) is utilized to map the relationship between the boundary conditions and POD coefficients. Therefore, the physical fields under the new boundary conditions could be calculated by the predicated POD coefficients from ANN and POD modes from snapshot. In order to assess the performance of an ANN-POD-based reduced order method, a simplified pressurized water reactor model under different inlet coolant temperatures and inlet coolant velocities is utilized. The results show that the new reduced order model can accurately predict the distribution of the physical fields, as well as the effective multiplication factor in the N/TH coupling nuclear system, whose relative errors are within 1%.

1. Introduction

Nuclear reactor is a complex multiphysics system, resulting from the coupling between neutronics, fluid dynamics, heat transfer, and mechanics. The neutronics and thermal-hydraulics coupling system has attracted significant research attention due to its great influence 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 [13]. Specifically, in Picard iteration, the coupling boundary information is exchanged between different physical fields to consider the coupling effect, while each physical field is still solved by the original existing code. However, the outer iterations are required among all the physical fields 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 fields are solved together to achieve the 2nd order convergence rate [69]. 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 different 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 efficiently 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 diffusion 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 differential equations (PDEs) about physical fields have been transformed to ordinary differential equations (ODEs) about the corresponding coefficient vectors which still have to be solved by numerical methods, called as the POD-Galerkin-based ROM method. Several attempts have been made to avoid solving the equations by using the artificial 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 fields, 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.

Therefore, in this paper, a nonintrusive POD-ANN method for the N/TH coupling problem has been developed. The POD method is employed to capture the characteristic spatial structure and modes of physical fields and eigenvalues under different boundary conditions by the snapshot method. Then, the backpropagation neutral network is selected to map the relationship between the boundary conditions and the POD coefficients. Finally, the trained ANN model is utilized to predict the POD coefficients with new boundary conditions, and then the distributions of coupled multiphysics fields could be reconstructed.

The 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 simplified N/TH coupling model. Finally, the conclusion of this work is summarized in Section 4.

2. Numerical Methods

2.1. Neutronics and Thermal-Hydraulics Coupling Model
2.1.1. Conservation Equations

A simplified pressurized water reactor model in a simplfied cylinder r-z two-dimensional water reactor model has been considered in this work. The 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 simplified 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. The 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 five nonlinear partial differential equations.

Neutron diffusion equation with constraint equation:where P is the fixed 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:

The pressure Poisson equation used in this work is derived from the mass conservation equation (4) and the momentum equation (5):

The neutron cross sections in equation (1) are the functions depending on the fuel temperature and the coolant temperature. The heat conduction of the solid fuel in equation (2) depends on the fuel temperature itself. The meaning of each term is listed in Table 1. The expression of the parameters can be seen in [22].

2.1.2. Numerical Discretization Method

The finite volume method is utilized in this work to discretize the nonlinear PDEs. Particularly, the convection term in equation (6) is discretized with the upwind difference scheme, and the central difference scheme is utilized for diffusion terms in these equations. Besides, the staggered grid mesh is employed to solve the coolant flow field distribution and coolant temperature field distribution. After discretization, the matrix form of the nonlinear equations is expressed by the following equation [22].

In equation (8), the term is the unknowns, is the coefficient matrix, and the submatrix blocks of which are dependent on the unknowns . The information about the boundary conditions and the source term is included in , which is not a constant vector obviously. The Picard iteration is utilized to solve N/TH equation system (8), where the coupling information, such as the fission power, fuel temperature, and coolant temperature, is exchanged among different physical fields. Specifically, the effective multiplication factor in the neutron diffusion equation is solved numerically with the power iteration method.

2.2. POD-ANN-Based ROM Method
2.2.1. POD Method

The ROM method aims to replace the original full-order problem with an alternation of the lower order but accurate model. The 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. The coupling equation system can be expressed by a function , where is the domain of the space and N is the size of , which depends on the number of the discretization points of the partial differential equations. In this work, the magnitude for N is . is the parameters which equal to boundary conditions in this work. The basic mechanic of POD is to find a set of deterministic function called POD basis which can best approximate the concerning function in the least square sense. It is the optimization problem:where is the POD coefficient and is the standard Euclidean.

The Snapshot method is one way to achieve this purpose. Let be a matrix consisting of the snapshots of the solution. In this work, the snapshots are collected from the numerical simulation of the full-order model constructed in Section 2.1.1 with a series of different boundary conditions. is the number of the snapshots, whose value is usually much smaller than that of .

We apply the singular value decomposition to S, and the basic vectors are the left singular vectors of and the first eigenvector are chosen to be the POD basis vectors. The number of the POD basis is determined by the predefined energy ration :

With the POD basis, the original function can be reconstructed by

2.2.2. ANN Backpropagation Neutral Network

The backpropagation neutral network (BPNN) is a powerful artificial 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. The activation functions of the hidden layer and the output layer used in this work are “tanh” and “identity,” respectively.

The accuracy of the BPNN highly depends on the quantity of hidden layers and quantity of neurons at each hidden layer. Theoretically, 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 overfitting issue. At present, the trial-and-error method is the most common and practical way to determine the number of hidden layers and number of nodes in the hidden layers. The basic structure of one hidden layer BPNN used in this work is elucidated in Figure 1. The four training processes of the BPNN are included as follows: firstly, 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.

2.2.3. 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. The 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 different perturbed boundary conditions to establish a library of the required results. Then, the POD method is utilized to extract the features of the physical field based on the simulation results, where the POD basis and corresponding coefficients are obtained. After that, the BPNN is used to map the relationship between the boundary conditions and POD coefficients. Finally, the trained BPNN is utilized to predict the POD coefficients as outputs, with new boundary conditions as inputs. The predicted POD modes and POD basis are then used to reconstruct the physical fields, 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.

3. Numerical Results

3.1. Problem Description and Sampling Design

A simplified 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. The neutronics and thermal-hydraulics coupling behavior is described by equations (1)–(7). The vacuum boundary condition is defined at the outface boundaries for the neutron diffusion equation, and the adiabatic boundary condition is defined at the boundary for solid energy conservation. The coolant pressure at outlet boundary is set to be 15.5 MPa, and the inlet boundary is set to be fixed velocity and fixed temperature boundary condition. The boundary conditions are listed in Table 2. The 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 different inlet temperatures and speeds of the coolant is selected to generate the snapshots. The inlet coolant temperature ranges from 280°C to 300°C, and the inlet coolant velocity is within the range of to . The intervals of the coolant temperature and coolant velocity are and , respectively. Therefore, 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 to and the inlet coolant velocity is within the range of to .The intervals of the coolant temperature and coolant velocity are and for testing cases, whose total number is 100. The distributions of the training and testing sample points are presented in Figure 4. The locations of the testing sample are totally different from those of the training sample.

3.2. POD Mode Features

The POD method is performed on the snapshots generated by the 121 training sets to obtain the POD modes which include the important features of different physical fields in the N/TH coupling system. The decay in the eigenvalue of the correlation matrices for each physical quantity is presented in Figure 5. It can be seen that the eigenvalues decay rapidly, which means only a few POD modes are sufficient to represent the characteristics of the physics fields. The number of the POD basis is determined by the energy ration in equation (11). In this work, , indicating that the error between the numerical results and the reconstructed results is less than 0.1% in norm 2. For the neutron flux field, only the first three POD modes are required to satisfy the criterion of energy ration . The dimensions of the reduced spaces for the approximation of the different 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 first or first two POD basis could capture the key features of the snapshots. We take the first three POD basis of the neutron flux as an example, Figure 6(a) shows that the first order POD basis, providing the basic spatial distribution of the neutron flux. The first 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 flux distribution. The 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 effect. Compared with the first 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 flux field. The POD bases of other physical fields have the similar phenomenon. Please note that the effective multiplication factor, which is a key parameter for steady-state N/TH calculation, is also considered here.

3.3. Artificial Neural Network Structure

The POD coefficients of high order change drastically with inlet velocity u0 and inlet coolant temperature T0. As an example, Figure 7 shows the 3rd POD coefficient of fuel temperature under different inlet velocities. In this work, the BPNN is utilized to calculate the POD coefficients under different coolant inlet parameters. Please note that since the complicated relationship between POD coefficient 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 coefficient 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. The Levenberg–Marquardt method is used as the training function due to its high efficiency for small and medium-size problems. The performance of the network is evaluated by the mean square error (MSE) between the target value and the predicted value of the outputs. The term n is the number of the outputs, that is, the POD coefficients. The information about residual convergence of each network is presented in Figure 8. For neutron flux, there are only a few tens of epochs to achieve the predefined stop criterion, while it requires about several hundreds of epochs for fuel solid temperature, coolant temperature, coolant velocity, coolant pressure, and effective multiplication factor, perhaps due to the relatively complicated relationship between the POD coefficients and inputs (coolant velocity and temperature).

The approximate fields are reconstructed based on the predicted POD coefficients from the ANN and POD bases for the training samples to evaluate the performance of the POD-ANN model. Then, 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 fields 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 error (), which is defined aswhere is the predicted results obtained from the POD-ANN model and 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.

The maximum MRE and the mean MRE of each physical field are listed in Table 5, and those of maximum error and the mean error are presented in Table 6. The MRE and of each field are under 1%, and the mean MRE of the effective multiplication factor reached 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 error of coolant velocity in the testing cases are the largest, almost one order of magnitude larger than other physical quantities. This 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 efficiency 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 field except for the effective multiplication factor. As shown in Figures 9(a) and 9(b), the maximum MRE of the predicted neutron flux 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 flux in these regions, the relatively larger predicted errors have minor impact on the power level. The 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 error distribution of each case, including both the training and testing samples, is shown in Figure 10. The error of the training and testing cases has similar distribution. The 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. The slopes of u0 and T0 change from −1 to 0, and the sampling interval is 0.1. Therefore, the 121 samples are generated and randomly divided into training (70) and testing (51) sets. The POD bases are then extracted, and the corresponding coefficients are calculated. The distributions of the first POD basis are shown in Figure 11. The MRE and error of each physical field in the testing cases under linear distribution inlet parameters are shown in Tables 8 and 9, respectively. The max MRE and error of all the physical fields are within 1%, which indicates that the POD-ANN method also has the capability of accurate prediction for the complicated boundary conditions.

4. Conclusion

An efficient 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 simplified cylinder r-z two-dimensional pressurized water reactor model is considered with different inlet temperature and speed of coolant. In detail, the POD approach is employed to extract the characteristic of the physical fields by snapshots. The spatial distribution of the POD bases shows that the first or first two POD bases capture the main features of the snapshots. The higher order POD bases could further consider the local detailed distribution. The BPNN is constructed with the boundary conditions as input and the POD coefficients as output. The trained BPNNs are then utilized to predict the POD coefficients with new boundary conditions. To evaluate the accuracy of the POD-ANN-based reduced order model, the MRE and error between the predicted results from the reduced order model and the numerical solution of PDEs are analyzed. Results show the MRE and error of each field are under 1% which indicates that the proposed approach can accurately predict the distribution of the physical fields, as well as the effective multiplication factor , for this steady-state N/TH coupling issue. The 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 overfitting 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

The data used to support the findings of this study are available from the corresponding author on request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by The National Key R&D Program of China no. 2022YFB1903000, National Natural Science Foundation of China no. 12275150, Beijing Natural Science Foundation no. 1212012, and Research Project of China National Nuclear Corporation.