Passive system failures can be categorized into two types: equipment failure and physical process failure (Mezio, 2010). Given that passive systems rely entirely on inherent system attributes, the probability of equipment failure is relatively low, and the primary cause of failure is usually a physical process failure. Due to the lack of operational data for passive systems, uncertainty in data is one of the primary challenges in passive system reliability analysis (Nayak et al., 2008).

The physical processes of passive systems are typically described by mathematical equations like Equations 1, 2. The calculation of the failure probability can be expressed as: Y = y X (1) X = X 1 , X 2 , ⋯ , X n (2)

For a specific task, the system performance function is given by Equation 3: g X = A − y X (3) where A is the failure criterion of the passive safety system; X is the input parameter of the best estimate procedure; y X is the key operating parameter calculated by the program; and g X is the performance function representing the operating state of the system.

For notational convenience, let Y > A denote the failure condition. The performance function can be defined as follows: g X = 0 represents the limit state or failure surface, g X < 0 represents the space of failure events, and g X > 0 represents the space of normal operation. Therefore, the failure probability can be evaluated through the following integral: P f = ∫ g X < 0 f X d X (4) where f X represents the joint probability density of the input parameter X . If all input parameters are independent, Equation 5 will be obtained: f X = f 1 X 1 f 2 X 2 ⋅ ⋅ ⋅ f n X n (5)

Traditionally, the Monte Carlo simulation method is used to estimate the failure probability defined in Equation 4. By establishing a model of the system under study, sampling the basic variables, and conducting a large number of experiments on the model, the value of y X is obtained. Then, according to the selected failure criterion, the g X < 0 cases are calculated and counted, and the ratio of this count to the total number of simulations N is calculated, which is the estimated value of the failure probability P f . The larger the value of N, the closer P f is to the true failure probability. The accuracy of this estimation method can be evaluated by its variance.

For a small number of simulation cycles, the variance may be relatively large. To obtain an acceptable estimate, the sample size should be at least 1 / P f (Xie et al., 2007). However, most passive systems have a very small failure probability, which means that a large number of computational simulations are required (Liu, 2015). Since the system analysis program takes a long time to calculate the results of a single sample, the traditional method is time-consuming. Additionally, for this problem, the system performance function cannot be expressed by a specific mathematical formula, making it impossible to obtain the failure probability simply by solving Equation 4.

To address the aforementioned issue, a rapidly computable surrogate model can be established to simulate the response of the thermal-hydraulic program. This type of model is generally referred to as a Response Surface (RS).

The fundamental idea of the response surface method is as follows: First, sampling is conducted based on the selected input parameters and their distributions, and the corresponding output results are calculated using the best-estimate code. Then, through regression analysis, a “response surface” is constructed to represent the relationship between the output results and the various input parameters.

The response surface method can replace complex thermal-hydraulic programs for calculations, significantly saving computational resources and time, and thus improving the efficiency of reliability analysis. In the application of the response surface method, the variation of the dependent variable is generally determined by multiple independent variables. The essential of this method lies in constructing a mathematical model between the independent and dependent variables, and further investigating the key factors affecting the response surface and their mechanisms, thereby achieving the optimization of the response surface. In this paper, the common polynomial regression and artificial neural network methods are used to construct the response surface.

Figure 1 shows the reliability analysis process using the response surface method.

First, the passive system is identified, and its functions, failure criteria, and performance indicators are clarified, laying the foundation for subsequent simulations. Second, a best-estimate code is developed that complies with internationally recognized guidelines and is consistent with prior information. Third, the relevant parameters affecting system uncertainty are defined, and uncertainty is screened, labeled, and quantified based on the available information. The randomness is described using probability density functions, and the sampling method is determined.

Fourth, parameter sampling is performed, and the sample size is determined based on reliability and confidence requirements. The sampling results are grouped to form an input matrix. Fifth, the input matrix is integrated into the developed model and the best-estimate code is executed. Sixth, the performance indicators of each model are calculated, and an output vector is generated. Seventh, the reliability of the passive safety system is analyzed. Eighth, a new input vector is obtained using the sampling method, and its corresponding performance indicators and output vector are obtained through the response surface (step nine). Subsequent work is to determine whether the calculation converges and to improve the calculation accuracy.

This paper uses VITARS as the best-estimate code to model and calculate. Referring to the overall VITARS coolant system model developed by Lu Guoqing (Liu, 2015), the reactor core, pressurizer, steam generator, and secondary-side passive residual heat removal system relevant to this problem are remodeled based on the requirements of this study.

The passive residual heat removal system studied in this paper is an important measure for preventing severe accidents in the advanced pressurized water reactor. In the event of a loss of feedwater to the steam generator or the station black-out accident combined with the failure of auxiliary feedwater steam-driven pumps, the reactor cooling system is no longer capable of continuously removing heat, resulting in a sharp rise in the reactor temperature and pressure to rise continuously, and thus preventing the normal operation of the residual heat removal system (Wang et al., 2012). In this emergency situation, the PRS system is activated. The height difference between the PRS heat exchanger and the steam generator forms a natural circulation, passively removing the core decay heat and the stored heat of various equipment in the reactor coolant system, further reducing the primary circuit temperature and pressure, and thus avoiding a loss-of-coolant accident that may be caused by primary circuit overpressure. In addition, a cooling water tank is set outside the containment to provide sufficient cooling water for the heat exchanger, ensuring that the PRS system can continue to operate self-sustainably for 72 h after an accident, providing valuable buffer time for accident management and emergency response. The secondary-side passive residual heat removal system is installed on the secondary side of the steam generator in all three coolant loops. The VITARS modeling diagram is shown in Figure 2.

Based on the aforementioned characteristic parameters, the formed debris bed in its natural state exhibits a significant extent of spreading. With an increase in the mass of molten material, the proportion of large-sized fragments in the formed debris bed also increases. Under the constant pressure conditions, as the subcooling of water decreases, the proportion of large-sized fragments in the debris bed increases. Since the molten material is sufficiently cooled in water, the high-temperature impact on the bottom plate of the water tank is relatively small.

In this simulation, to better understand the safety performance of nuclear power plants under extreme conditions and to provide a scientific basis for preventing and responding to similar accidents, we assume that both the auxiliary feedwater system and the emergency core cooling system fail after the accident occurs. Therefore, the residual heat of the reactor will be mainly removed by the PRS system.

The accident sequence is shown in Table 1.

Figures 3, 4 show the water level of the cooling water tank and the peak temperature of the coolant outlet over the entire 72 h after the accident, respectively. Figure 3 shows that the water in the cooling water tank is nearly depleted at around 70 h, while Figure 4 shows that the coolant outlet temperature has not shown an upward trend, indicating that the water volume setting of the PRS cooling water tank is sufficient. In addition, the reactor outlet temperature reaches its peak within 1 h, which means that the analysis of the results within 1 h has already covered the most dangerous situation of the reactor. In subsequent batch calculations, in order to save time, only the reactor parameter response in the first hour is calculated.

During this process, the temperature of the outermost layer of the fuel cladding never exceeds 672K, and the coolant in the core region is always in a liquid state, and the core is not exposed, indicating that the core is always in a safe state. Within 72 h after the accident, the outlet temperature of the reactor coolant can be reduced to nearly 400K, and natural circulation can be maintained in both the primary and secondary sides of the fluid, so it is considered that the reactor coolant system has entered a safe state.

Critical parameters are those that have the greatest impact on system performance. They can generally be divided into operating parameters and structural parameters. For the secondary-side passive residual heat removal system, based on literature research and expert judgment, the input parameters that need to be considered are shown in Table 2.

Since the distribution ranges of general system critical parameters (such as heat transfer area) are not completely symmetrical, their actual parameter ranges may deviate from the standard normal distribution. In this case, a truncated normal distribution is more reasonable, and its distribution is shown in Figure 5.

Therefore, the normal distributions used in this paper are all truncated normal distributions, and the probability density function can be expressed as ψ(μ,σ,a,b;x), and the specific expression is as shown in the Equation 6: ψ μ ¯ , σ ¯ , a , b ; x = 0 x ≤ a ; ϕ μ ¯ , σ ¯ 2 x Φ μ ¯ , σ ¯ 2 ; b − Φ μ ¯ , σ ¯ 2 ; a a < x < b ; 0 b ≤ x . (6) where μ ¯ and σ ¯ are the mean and standard deviation, respectively, of the standard normal distribution; a and b are the lower and upper bounds, respectively, of the truncated normal distribution parameter values.

In the reliability analysis of passive systems, the most important step is to determine the failure criterion of the system. For the PRS system studied in this paper, its key task is to remove the decay heat generated after the reactor shutdown through natural circulation under the accident conditions, and to reduce the coolant temperature and the pressure of the primary cooling system timely and effectively, preventing the occurrence of local boiling in the coolant inside the core, and thus ensuring the integrity of the reactor core cladding (Wang, 2022). Therefore, this study takes the peak temperature of the coolant at the core outlet as the failure criterion of the PRS system. The failure region can be determined as shown in the Equation 7: F = Ω : T p e a k > 618 K (7)

After determining the uncertainty of the critical parameters, Latin Hypercube Sampling (LHS) is employed to sample the critical parameters listed in Table 2. According to the research object, the key parameters mainly come from PRS heat exchangers (PRHX).

LHS draws samples from the entire distribution of random variables and produces more stable and reliable results. In other words, LHS can obtain good sampling results with fewer samples (Huang and Kuang, 2012).

After sampling, a self-written script is used to batch modify the relevant parts of the key parameters in the input cards, and 100 groups of input cards are obtained. Then, the VITARS program is used to calculate these 100 groups of input cards, and the operating data of the reactor is obtained. The output parameters of interest are extracted, and finally, the input parameters and corresponding output parameters of each calculation are grouped as a set for subsequent analysis.

Preliminary batch calculations were conducted, and the results are shown in Figure 6. It was found that within the range of uncertain parameters, the coolant outlet temperature generally reached its peak at around 500 s, with the peak temperature fluctuating slightly around 592 K, which is far below the failure temperature. Therefore, an attempt was made to find the worst-case condition within the range of uncertain parameters to determine whether the PRS system had a probability of failure within the range of uncertainty. After a series of sampling and calculations, the worst-case condition within the range of uncertainty was found, and the coolant temperature response without the PRS system was calculated. A comparison of the two is shown in Figure 7. The results come from small batches simulations using Thermal-Hydraulic programs.

Without PRS system, the coolant outlet temperature began to rise at around 1800 s. This was due to the loss of all feedwater, causing the water in the steam generator to continuously boil while absorbing the heat generated by the core. Finally, the water evaporated at around 1,800 s, and the coolant lost its ultimate heat sink, leading to a temperature increase. The reactor coolant would then begin to boil. However, under the worst-case condition with the PRS system, where the peak temperature of the coolant outlet is highest within the uncertainty parameter range, the friction factor and inlet/outlet resistance coefficient in the heat exchanger tubes were both at their maximum values within the range of uncertainty. Due to tube blockage, the effective flow area was at its minimum value. Due to scaling, the tube thickness increased, and the initial temperature of the heat exchanger water tank was 40°C. Even under these conditions, the PRS system could still effectively remove the heat generated by the core. According to the calculations, the coolant outlet temperature would not rise again within 72 h, and the water in the tank would not be depleted. Therefore, it can be concluded that the PRS is not prone to failure within the range of uncertain parameters and has superior safety.

Therefore, some extreme conditions are superimposed to further test the reliability of PRS system. The coolant temperature decreases due to the opening of the steam generator safety valve and the pressure drop. Therefore, it is feasible to delay the opening of the steam generator safety valve. Considering the maximum allowable pressure of the steam generator, the delay opening time is limited to within 1 min.

Before establishing the response surface, sensitivity analysis is needed to compare the importance of each input parameter to the output parameters. Input parameters with relatively low importance can be directly set to fixed values after sensitivity analysis, thereby reducing the model complexity and minimizing the interference of unnecessary input parameters on the model. In this paper, Spearman’s rank correlation coefficient is selected as the basis for sensitivity analysis.

Spearman’s rank correlation coefficient is used to estimate the strength of the correlation between two random variables and is generally represented by ρ. Its value ranges from [-1, 1]. The stronger the correlation between the two variables, the closer the absolute value of the rank correlation coefficient is to 1. A positive value indicates a positive correlation, a negative value indicates a negative correlation, and a value close to 0 indicates no correlation (Chen et al., 2013).

Sensitivity analysis was performed on the above input parameters and corresponding output parameters, and the results are shown in Figure 8. The output parameter is the failure criterion, which is the peak temperature of the coolant outlet.

Parameters numbered 1, 4, 6, and 9 are the reactor power, PRHX tube diameter, initial water temperature of the cooling water tank, and safety valve delay opening time, respectively. These four parameters have a strong correlation with the output parameters and are selected as the input parameters for response surface training. Given the significant differences in magnitudes among the various feature variables, normalization is necessary to mitigate their impact when used as training data for a response surface model.

The order of the highest-order term of the polynomial response surface can be selected according to the relationship between the response value and the input parameters. In general, the operation of a passive system is more complex, and a first-order model cannot accurately predict its operating results. Therefore, a second-order model is selected to establish the response surface model here (Kirchsteiger and Lavín, 2004).

Using the above 100 groups of data for training and testing, a self-written script is written for regression analysis, and the final result of the quadratic polynomial response surface is obtained as shown in Equation 8: Y X = 477.74 − 0.04 X 0 + 2509.91 X 1 + 0.54 X 2 − 0.01 X 3 − 0.83 X 0 X 1 − 10.61 X 1 X 2 − 0.72 X 1 X 3 + 16434.8 X 1 2 (8) where X 0 is the reactor power; X 1 is the PRHX tube diameter; X 2 is the initial water temperature of the cooling water tank; X 3 is the safety valve delay opening time.

When constructing a multi-layer neural network, a issue that needs attention is the number of hidden layers and the number of nodes in the hidden layer. The hidden layer plays an abstracting role and is responsible for extracting key features from the input parameters, and the number of nodes in the hidden layer directly determines the size and capacity of the network.

To enhance the ability of the neural network to process information, adding hidden layers is an effective way. However, this will also lead to an increase in the complexity of the training process, thus requiring more training data and a longer training period. For the response surface model of this problem, in order to replace the complex thermo-hydraulic program, efforts should be made to reduce the scale of the network structure in order to reduce the length of the training time. Generally, a BP neural network with a three-layer structure is already able to meet the needs of most response surface models. Based on these factors, this study chose to construct a neural network model with only a single hidden layer (Wang et al., 2025).

For the regression problem of this topic, a three-layer BP neural network is selected. The number of neurons in the input layer is 4, the number of neurons in the hidden layer is 9, and the number of neurons in the output layer is 1, which can be expressed as N (4, 9, 1), and the structure is shown in Figure 9. The optimal number of hidden layer neurons is obtained using Bayesian optimization. The principle of Bayesian optimization is based on Bayes’ theorem. First, a probability model of the objective function is constructed, and then this probability model is continuously updated with new sampling points. Based on the model, the next most promising sampling point is calculated, so as to efficiently find the optimal solution of the objective function (Pedroni and Zio, 2017; Lu, 2022).

After automatic parameter tuning using the Bayesian optimization method, the final configuration determined is that the hidden layer uses PReLU, the output layer uses Linear, and the Adam optimizer is used with an initial learning rate of lr = 0.01. For this neural network response surface, the mean squared error (MSE) is 0.008003.

The Adam optimizer used in this article not only achieves the function of adaptive learning rate, but also has simple implementation, high computational efficiency, and low memory requirements, making it an optimizer with excellent performance (Wang et al., 2018; Helton and Davis, 2003; Ye, 2020; Dong et al., 2024).

PReLU, also known as Parametric Rectified Linear Unit, is an improvement on ReLU. As shown in Figure 10, it introduces a learnable slope parameter on the negative half-axis, where represents the output of the previous layer. This design enables neurons to remain in an activated state even when the input is negative, effectively avoiding the “neuron death” problem that may be caused by excessive gradients during the backpropagation process of ReLU. In addition, PReLU also demonstrates many advantages such as high computational efficiency and stable gradient propagation, making it suitable for most deep learning tasks.

After the establishment and training of the response surface model, it is necessary to compare the prediction capabilities of the two response surfaces. By sampling 5,000 times for the input parameters and distributions determined in Section 3.5, new input samples are obtained and substituted into the established quadratic polynomial response surface and neural network response surface for calculation. Using the 100 groups of data generated by VITARS and the 5,000 groups of data output by each response surface, probability distribution diagrams as shown in Figure 11 and cumulative probability distribution diagrams as shown in Figure 12 are plotted.

From Figure 11, it can be seen that similar to the VITARS calculation data, the output results of the two response surfaces approximately follow a normal distribution, and the probability peak appears between 595K and 600K. However, compared with the quadratic polynomial response surface, the neural network can better cover the parameter range of the training set, reflecting the true distribution of the data. For the cumulative probability distribution in Figure 12, due to the scarcity of VITARS data samples, the cumulative distribution has abrupt changes and is not smooth, but the output of the neural network response surface is still very close to the VITARS data, and is sufficiently approximate the true distribution curve. In contrast, the cumulative distribution of the quadratic polynomial response surface lacks a temperature range of about 10K, which is significantly different from the actual distribution. Although the limited sample of the VITARS data results in abrupt changes rather than smooth transitions in its cumulative distribution, the output of the neural network response surface is still more similar to the VITARS data.

To further compare the differences in the fitting capabilities of the two response surface models, the input parameters of the original 100 groups of data are used as the input of the response surface, and the predicted values output by the two response surfaces are compared with the real output of VITARS. In Figure 13, a Parity plot is shown, which is a scatter plot that compares the differences between the model results and the baseline data. A baseline is added to the figure to represent the perfect prediction value of the model. When the predicted point of the model falls on the baseline, it means that the prediction error of the model at that point is 0, and the farther away from this line, the greater the error. In addition, an error interval is also drawn. The predicted values of the neural network response surface are all within the 1% error interval, while only 90% of the predicted results of the quadratic polynomial are within the error interval, which indicates that the prediction capability of the neural network response is also better than that of the quadratic polynomial.

In general, the prediction results of the neural network response surface can well reflect the probability distribution characteristics of the parameters, and compared with the baseline values, the error can be controlled within 1%, which can meet the large-scale and high-precision calculation requirements of the passive system reliability analysis.

Based on the above comparison, the neural network response surface with better prediction performance is adopted as the surrogate model of VITARS. After extracting 5,000, 10,000, 50,000, 100,000, and 200,000 groups of sampling data using Latin hypercube sampling, the trained neural network response surface is used for prediction. To reduce the impact of different results for each sampling, each group of calculations is performed 5 times, and the average value of the results is taken. The failure probability is summarized in Table 3. It can be seen that the failure probability converges at 100,000 calculations, and the failure probability is 0.035%, which is acceptable. The prediction results of 10,000 and 100,000 times are shown in Figure 14.