Nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraint

Abstract

The invention discloses a nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraint, and relates to the field of nuclear reactor xenon oscillation prediction. The problems that in the prior art, xenon oscillation prediction is difficult, a physical model method can only be suitable for a specific reactor type, and a data driving method needs high calculation cost are solved. The invention establishes a xenon oscillation agent model with physical constraint based on data assimilation. The agent model is divided into two parts, namely a parameter search part and a xenon oscillation process curve prediction part. The parameter searching part refers to that the model searches global optimal parameters of multiple groups of working conditions through a data assimilation method under the limitation of physical constraints. According to the xenon oscillation process curve prediction part, the xenon oscillation process under the unknown working condition is predicted through the fitting relation between the parameters. The method is also suitable for the technical field of agent models for predicting the xenon oscillation curve.

Description

Technical Field

[0001] The present invention relates to the technical field of nuclear reactor xenon oscillation prediction, and in particular to a nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints. Background Art

[0002] Xenon oscillations typically occur in large thermal neutron reactors with high neutron flux rates, such as large pressurized water reactors (PWRs). Radial xenon oscillations converge, and their stability increases with core burnup. However, the stability of axial xenon oscillations decreases with increasing core burnup or core power. Therefore, axial xenon oscillations are a worthy research topic in large thermal neutron reactors. In actual PWRs, xenon oscillations can cause periodic oscillations in the axial power deviation ΔI. When the power deviation ΔI exceeds ±3% FP (full power), a master control alarm is triggered. When the power deviation exceeds ±5% FP, the cumulative timer must not exceed 1 hour within 12 consecutive hours. Furthermore, xenon oscillations can shift the position of reactor heat pipes and alter the power density crest factor, leading to localized increases in power levels and temperatures. If uncontrolled, these oscillations can even cause fuel element meltdowns. Xenon oscillations can also cause alternating temperature fluctuations in the core, exacerbating thermal stresses in core materials and causing premature material failure. Therefore, in order to ensure the safe and stable operation of the unit, xenon oscillation must be discovered early and correct intervention measures must be taken.

[0003] Xenon oscillations involve complex nonlinear processes, including the interaction between neutron flux, iodine-135, and xenon-135, which are difficult to describe using simple linear models. During reactor operation, the concentration of xenon-135 cannot be directly measured, and the xenon distribution can only be inferred through indirect methods (such as neutron flux measurements), which increases the uncertainty of predictions. Therefore, effective prediction methods are needed. By accurately predicting xenon oscillations, measures can be taken in advance to suppress the oscillations and avoid potential safety risks. Accurately predicting and effectively controlling the xenon oscillation process is of great significance for ensuring reactor safety and optimizing operation.

[0004] Among them, the physical model method proposed in the prior art is based on two groups of one-dimensional diffusion equations and the concentration variation equations of iodine and xenon. It uses Fourier sine series to expand the spatial distribution of xenon, iodine, and flux. The axial difference parameter is introduced to transform the equation into a non-homogeneous equation. This method proposes an analytical solution for the time-varying variation of the reactor axial power shape index (ASI) during xenon oscillation. Once the cross-sectional data under equilibrium xenon conditions and the corresponding physical constants are known under the corresponding operating conditions, the reactor axial power shape index (ASI) at the initial time and its corresponding derivative are input to predict the xenon oscillation process.

[0005] However, the physical model approach has the following drawbacks: Initializing the equilibrium xenon requires a sufficient amount of pre-measured data to fit the cosine function, and core conditions cannot be changed to validate this data. This approach has limited applicability to operating conditions and is generally applicable only to specific reactor types.

[0006] A hybrid dynamic prediction scheme combining a data-driven approach proposed in the prior art with a data assimilation method and dynamic mode decomposition (DMD) is used to predict the full-core power distribution caused by xenon oscillations in the HPR1000 reactor. This coupled scheme is divided into two stages: (1) Reconstruction stage: In the initial stage, only the IDB-DA method is used for physical field reconstruction. (2) Prediction and reconstruction stage: In the subsequent stage, the detector node power values ​​predicted by DMD and the IDB-DA method are combined to predict the full-core power distribution.

[0007] However, the disadvantages of the data-driven method are: (1) The implementation of the data assimilation method is relatively complex, requiring the selection of appropriate background physical fields and hyperparameters such as regularization factors, and requiring a large amount of experimental data and adjustments. When the number of detectors configured in the core is limited, insufficient data will lead to increased prediction errors. This method is a data-driven method and is not subject to physical constraints, so the credibility of the predicted results is relatively poor. (2) Although the DMD method has high computational efficiency, the IDB-DA method may require high computational costs when processing large-scale physical fields. (3) The performance of data assimilation is highly dependent on the quality of the observation data. If there is large noise or error in the observation data, it may affect the accuracy of the prediction results. Summary of the Invention

[0008] The present invention addresses the problems in the prior art such as the difficulty in predicting xenon oscillations, the problem that the physical model method can only be applied to specific reactor types, and the problem that the data-driven method requires high computational costs.

[0009] To solve the above technical problems, the present invention is achieved through the following technical solutions:

[0010] The present invention proposes a nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints, the method comprising the following steps:

[0011] S1. Select five xenon transient operating conditions with the same reactor power level and different degrees of sudden drop. Set the duration of each xenon transient step to 1 hour. Calculate the change of the reactor axial power shape index (ASI) and the corresponding physical parameters over time for each operating condition under 38 xenon transient steps, and output the corresponding results.

[0012] S2. Setting up a proxy model, inputting the change of the reactor axial power shape index ASI over time and the corresponding physical parameters for each group of operating conditions in S1 into the proxy model;

[0013] The physical parameters include the reactor core height H, the decay constants of iodine and xenon nuclides, the diffusion coefficients D1 and D2 of the full reactor fast group and thermal group, and the absorption cross sections Σ a1 ,Σ a2 , fission cross section of full reactor fast group and hot group Σ f1 ,Σ f2 , the removal cross section of the whole stack Σ R , the number of neutrons produced per fission in the reactor ν and the reaction cross section σ of the nuclide xenon in the thermal group Xe ; and calculate the average value of the above parameters after averaging over the entire stack;

[0014] S3. Define the search parameters for the yield of iodine and xenon γ I , γ Xe , the average nuclear concentrations of iodine and xenon Average neutron flux of the whole reactor fast group and thermal group and the power coefficient α in the reactor p The value range of

[0015] S4. Executing the parameter search part of the proxy model, the proxy model performs a global search and dimensionality reduction of parameters for the five input working conditions, and performs a local correction of the search parameters for each input working condition;

[0016] S5, judging whether the relative standard deviation of the final output search parameters of the five groups of working conditions is within 5%, if so, proceeding to S6; if not, repeating the process of S4;

[0017] S6. The proxy model draws the ASI-t curve corresponding to each set of working conditions in combination with the search parameter value finally output by S4, and performs time-space translation and scaling transformation based on the original ASI-t curve under the working condition;

[0018] S7. Determine whether the improvement rate η of the MSE for each group of working conditions before and after the temporal and spatial translation correction and the scaling correction is greater than 10%;

[0019] S8, execute the xenon oscillation process curve prediction part of the agent model, combine the input working conditions and the output parameters to establish the power transient value ΔP and E Pd (0), α p and And the linear interpolation expressions of ΔP and the time translation value Δt, the spatial translation value Δy and the scaling value α between any two working conditions;

[0020] S9, input the ΔP value of the working condition to be predicted, and obtain the E of the working condition to be predicted through the fitting relationship obtained in S8. Pd (0), α pThe predicted values ​​of the parameters and the time translation value Δt, space translation value Δy and scaling value α are input into the ASI calculation model to obtain the ASI-t curve under the predicted working conditions.

[0021] Furthermore, a preferred embodiment is provided, in which the method for calculating the change of the reactor axial power shape index ASI over time and the corresponding physical parameters of each group of operating conditions in S1 is as follows:

[0022]

[0023] Where, P B With P T They represent the power values ​​of the lower and upper parts of the reactor axis respectively.

[0024] Furthermore, a preferred embodiment is provided, wherein S2 further includes a step of calculating the ASI value of the xenon transient operating condition determination curve at time t=0, and the ASI value of the xenon transient operating condition determination curve at time t=0 is named E P (0) and the derivative E Pd (0), the derivative value is the ASI change during the first xenon transient step, and its calculation expression is:

[0025]

[0026] Furthermore, a preferred embodiment is provided, in which the search parameters iodine and xenon yield γ are defined in S3. I , γ Xe , the average nuclear concentrations of iodine and xenon Average neutron flux of the whole reactor fast group and thermal group and the power coefficient α in the reactor p The value range does not exceed 0.01 times and 100 times the actual measurement value inside the reactor.

[0027] Furthermore, a preferred embodiment is provided, wherein S4 also includes the agent model for the average neutron flux of the full reactor fast group and the thermal group that varies with time and space. with the average nuclear concentrations of iodine and xenon The steps of introducing the axial Fourier expansion and axial difference parameters,

[0028] The expressions are:

[0029]

[0030] E i =b i,2 (t)-b i,2 (∞)

[0031] Where i=1, 2, 3, 4 represent Four variables, n represents the order of the Fourier expansion coefficient, b i,2 (∞) represents the value of the second-order Fourier expansion coefficient of the variable under equilibrium conditions; E i Axial difference function representing the coefficients of the second-order Fourier expansion of the variable.

[0032] Furthermore, a preferred embodiment is provided, in which the method for performing local correction of search parameters for each input working condition in S4 is: performing global search and local correction of the parameters of the working condition based on a differential evolution-SLSQP hybrid algorithm.

[0033] Furthermore, a preferred embodiment is provided, in which the method for performing the time-space translation transformation and scaling transformation based on the original ASI-t curve under the working condition in S6 is:

[0034] The expression that defines the transformation function is:

[0035]

[0036] Among them, Δt, Δy, and α represent the time transformation value of the curve respectively. The default positive value is the left translation and the spatial transformation value. The default positive value is the left translation and scaling transformation value. With 1 as the base value, each value in the curve is multiplied by the scaling value to obtain the scaled new curve. is the ASI value after transformation;

[0037] Δt∈[-2,2],|Δy|≤0.1max(|y1|),α∈[0.5,1.5]

[0038] Furthermore, a preferred embodiment is provided, in which the method for determining in S7 whether the improvement rate η of the MSE of each group of working conditions before and after the temporal and spatial translation correction and the scaling correction is greater than 10% is:

[0039]

[0040] If so, the corrected time shift value Δt, spatial shift value Δy and scaling value α are output; if not, the output results of the time shift value Δt and spatial shift value Δy are 0, and the scaling value α is 1.

[0041] Solution 3: A computer device includes a memory and a processor, wherein the memory stores a computer program. When the processor runs the computer program stored in the memory, the processor executes any one of the methods described in Solution 1.

[0042] Solution 4: A computer-readable storage medium storing a computer program, wherein the computer program, when executed by a processor, implements the steps of the method described in any one of Solution 1.

[0043] The present invention is beneficial in that:

[0044] The nuclear reactor xenon oscillation prediction method described in this paper, based on data assimilation and physical constraints, establishes a proxy model suitable for predicting xenon oscillation curves. Combining the advantages of physical modeling and data-driven approaches, it is highly adaptable to reactor types and can produce prediction results without extensive computation. It also predicts the evolution of xenon oscillations within a certain range. This method combines the advantages of physical modeling and data-driven approaches, and the results obtained are of certain reference value.

[0045] The proxy model of the present invention automatically reduces the dimensionality of the parameters required for input, thereby reducing the complexity of subsequent variable fitting and reducing a certain amount of calculation.

[0046] The proxy model of the present invention has low computational cost and high computational efficiency. It takes an average of 2 to 3 minutes to obtain search parameter results corresponding to five groups of working conditions, and less than 10 seconds to obtain the results of the subsequent parameter fitting and prediction parts.

[0047] The proxy model described in the present invention can predict the evolution trend of xenon oscillation at any point within the operating condition range through several operating conditions, which can greatly simplify the calculation amount of core power capability analysis and provide a certain reference for the power change during the xenon oscillation process in the core.

[0048] The present invention is also applicable to the technical field of proxy models for predicting xenon oscillation curves. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 This is a flow chart of the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in embodiment 1.

[0050] Figure 2 Schematic diagram of the original ASI-t curves of the five calculation conditions at the initial power level of 100% FP described in the eleventh embodiment.

[0051] Figure 3 This is a comparison diagram of the ASI-t curves of the predicted working condition of the initial power level 100% FP described in the eleventh embodiment.

[0052] Among them, (a) is the comparison diagram of the predicted curve and the actual curve when the initial power level is 100% FP and suddenly drops to 65% FP, (b) is the comparison diagram of the predicted curve and the actual curve when the initial power level is 100% FP and suddenly drops to 55% FP, (c) is the comparison diagram of the predicted curve and the actual curve when the initial power level is 100% FP and suddenly drops to 45% FP, (d) is the comparison diagram of the predicted curve and the actual curve when the initial power level is 100% FP and suddenly drops to 35% FP, (e) is the comparison diagram of the predicted curve and the actual curve when the initial power level is 100% FP and suddenly drops to 25% FP, and (f) is the comparison diagram of the predicted curve and the actual curve when the initial power level is 100% FP and suddenly drops to 15% FP.

[0053] Figure 4 Schematic diagram of the original ASI-t curve of the five calculation conditions at the initial power level 90% FP described in the eleventh embodiment.

[0054] Figure 5 This is a comparison diagram of the ASI-t curves of the predicted working condition of the initial power level 90% FP described in the eleventh embodiment.

[0055] Among them, (a) is the comparison diagram of the predicted curve and the actual curve when the initial power level of 90% FP suddenly drops to 65% FP, (b) is the comparison diagram of the predicted curve and the actual curve when the initial power level of 90% FP suddenly drops to 55% FP, (c) is the comparison diagram of the predicted curve and the actual curve when the initial power level of 90% FP suddenly drops to 45% FP, (d) is the comparison diagram of the predicted curve and the actual curve when the initial power level of 90% FP suddenly drops to 35% FP, (e) is the comparison diagram of the predicted curve and the actual curve when the initial power level of 90% FP suddenly drops to 25% FP, and (f) is the comparison diagram of the predicted curve and the actual curve when the initial power level of 90% FP suddenly drops to 15% FP. DETAILED DESCRIPTION

[0056] In order to make the purpose, technical solutions and advantages of the implementation methods of this application clearer, the technical solutions in the implementation methods of this application will be clearly and completely described below in combination with the drawings in the implementation methods of this application. Obviously, the described implementation methods are only part of the implementation methods of this application, not all of the implementation methods.

[0057] Embodiment 1: This embodiment provides a method for predicting nuclear reactor xenon oscillations based on data assimilation and physical constraints, the method comprising the following steps:

[0058] S1. Select five xenon transient operating conditions with the same reactor power level and different degrees of sudden drop. Set the duration of each xenon transient step to 1 hour. Calculate the change of the reactor axial power shape index (ASI) and the corresponding physical parameters over time for each operating condition under 38 xenon transient steps, and output the corresponding results.

[0059] S2. Setting up a proxy model, inputting the change of the reactor axial power shape index ASI over time and the corresponding physical parameters for each group of operating conditions in S1 into the proxy model;

[0060] The physical parameters include the reactor core height H, the decay constants of iodine and xenon nuclides, the diffusion coefficients D1 and D2 of the full reactor fast group and thermal group, and the absorption cross sections Σ a1 ,Σ a2 , fission cross section of full reactor fast group and hot group Σ f1 ,Σ f2 , the removal cross section of the whole stack Σ R , the number of neutrons produced per fission in the reactor ν and the reaction cross section σ of the nuclide xenon in the thermal group Xe ; and calculate the average value of the above parameters after averaging over the entire stack;

[0061] S3. Define the search parameters for the yields of iodine and xenon γ I , γ Xe , the average nuclear concentrations of iodine and xenon Average neutron flux of the whole reactor fast group and thermal group and the power coefficient α in the reactor p The value range of

[0062] S4. Execute the parameter search part of the proxy model. The proxy model performs a global search and dimensionality reduction of parameters for the five input working conditions, and performs a local correction of the search parameters for each input working condition.

[0063] S5, judging whether the relative standard deviation of the final output search parameters of the five groups of working conditions is within 5%, if so, proceeding to S6; if not, repeating the process of S4;

[0064] S6. The proxy model draws the ASI-t curve corresponding to each set of working conditions in combination with the search parameter value finally output by S4, and performs time-space translation and scaling transformation based on the original ASI-t curve under the working condition;

[0065] S7. Determine whether the improvement rate η of the MSE for each group of working conditions before and after the temporal and spatial translation correction and the scaling correction is greater than 10%;

[0066] S8, execute the xenon oscillation process curve prediction part of the agent model, combine the input working conditions and the output parameters to establish the power transient value ΔP and E Pd (0), α p and And the linear interpolation expressions of ΔP and the time translation value Δt, the spatial translation value Δy and the scaling value α between any two working conditions;

[0067] S9, input the ΔP value of the working condition to be predicted, and obtain the E of the working condition to be predicted through the fitting relationship obtained in S8. Pd (0), α p The predicted values ​​of the parameters and the time translation value Δt, space translation value Δy and scaling value α are input into the ASI calculation model to obtain the ASI-t curve under the predicted working conditions.

[0068] Embodiment 2: This embodiment further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Embodiment 1. The method for calculating the change of the reactor axial power shape index ASI over time and the corresponding physical parameters for each group of operating conditions in S1 is as follows:

[0069]

[0070] Where, P B With P T They represent the power values ​​of the lower and upper parts of the reactor axis respectively.

[0071] Implementation method 3: This implementation method further limits the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in implementation method 2. S2 also includes the step of calculating the ASI value of the xenon transient operating condition determination curve at time t = 0, and the ASI value of the xenon transient operating condition determination curve at time t = 0 is named E P (0) and the derivative E Pd (0), the derivative value is the ASI change during the first xenon transient step, and its calculation expression is:

[0072]

[0073] Implementation 4: This implementation further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Implementation 1. In S3, the search parameters iodine and xenon yield γ are defined. I , γ Xe , the average nuclear concentrations of iodine and xenon Average neutron flux of the whole reactor fast group and thermal group and the power coefficient α in the reactor p The value range does not exceed 0.01 times and 100 times the actual measurement value inside the reactor.

[0074] Implementation 5. This implementation further limits the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Implementation 1. S4 also includes the agent model for the average neutron flux of the full reactor fast group and thermal group that varies with time and space. with the average nuclear concentrations of iodine and xenon The steps of introducing the axial Fourier expansion and axial difference parameters,

[0075] The expressions are:

[0076]

[0077] E i =b i,2 (t)-b i,2 (∞)

[0078] Where i=1, 2, 3, 4 represent Four variables, n represents the order of the Fourier expansion coefficient, b i,2 (∞) represents the value of the second-order Fourier expansion coefficient of the variable under equilibrium conditions; E i Axial difference function representing the coefficients of the second-order Fourier expansion of the variable.

[0079] Implementation method six. This implementation method further limits the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in implementation method five. The method for local correction of search parameters for each input working condition in S4 is: global search and local correction of the parameters of the working condition based on the differential evolution-SLSQP hybrid algorithm.

[0080] Embodiment 7. This embodiment further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Embodiment 2. The method for performing spatiotemporal translation and scaling transformation based on the original ASI-t curve under this operating condition in S6 is as follows:

[0081] The expression that defines the transformation function is:

[0082]

[0083] Among them, Δt, Δy, and α represent the time transformation value of the curve respectively. The default positive value is the left translation and the spatial transformation value. The default positive value is the left translation and scaling transformation value. With 1 as the base value, each value in the curve is multiplied by the scaling value to obtain the scaled new curve. is the ASI value after transformation;

[0084] Δt∈[-2,2], |Δy|≤0.1max(|y1|), α∈[0.5,1.5].

[0085] Embodiment 8. This embodiment further limits the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Embodiment 2. In S7, the method for determining whether the improvement rate η of the MSE for each group of operating conditions before and after the spatiotemporal translation correction and scaling correction is greater than 10% is as follows:

[0086]

[0087] If so, the corrected time shift value Δt, spatial shift value Δy and scaling value α are output; if not, the output results of the time shift value Δt and spatial shift value Δy are 0, and the scaling value α is 1.

[0088] Embodiment 9: A computer device includes a memory and a processor, wherein the memory stores a computer program. When the processor runs the computer program stored in the memory, the processor executes the method described in any one of embodiments 1 to 8.

[0089] Embodiment 10: A computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps of the method described in any one of embodiments 1 to 8 are implemented.

[0090] Implementation 11: The examples provided in this implementation are used to explain the above implementations 1 to 10, and specifically include the following contents:

[0091] The purpose of this implementation is to establish a xenon oscillation proxy model with physical constraints based on data assimilation. The proxy model consists of two parts: a parameter search component and a xenon oscillation process curve prediction component. The parameter search component uses data assimilation to search for globally optimal parameters for multiple operating conditions within the constraints of physical constraints. The xenon oscillation process curve prediction component uses the fitted relationship between parameters to predict the xenon oscillation process under unknown operating conditions.

[0092] This embodiment specifically includes the following steps:

[0093] S1: Select five sets of xenon transient operating conditions with different degrees of sudden drop in reactor power level. Set the duration of each xenon transient step to 1 hour. Based on the traditional nuclear design program, calculate the change of the reactor axial power shape index (ASI) over time and the corresponding physical parameters for each set of operating conditions under 38 xenon transient steps, and output the corresponding results. The calculation expression of ASI is:

[0094]

[0095] Where, P B With P T They represent the power values ​​of the lower and upper parts of the reactor axis respectively.

[0096] S2: Input the change of reactor axial power shape index ASI over time and the corresponding physical parameters for each set of operating conditions into the proxy model.

[0097] The required physical parameters include the reactor core height H, the decay constants of iodine and xenon nuclides, the diffusion coefficients D1 and D2 of the full reactor fast group and thermal group, and the absorption cross sections Σ of the full reactor fast group and thermal group. a1 ,Σ a2 , fission cross section of full reactor fast group and hot group Σ f1 ,Σ f2 , the removal cross section of the whole stack Σ R , the number of neutrons produced per fission in the reactor ν and the reaction cross section σ of the nuclide xenon in the thermal group Xe .

[0098] These parameters are the average values ​​obtained after the whole stack is averaged. At the same time, it is also necessary to determine the ASI value of the curve at time t=0 based on the calculated xenon transient conditions (named E P (0)) and the derivative E Pd (0), the derivative value is the ASI change during the first xenon transient step, and its calculation expression is:

[0099]

[0100] S3: Define the dimension and value range of the search parameters. The number of search parameters is the dimension of the search parameters, which is generally taken as the yield of iodine and xenon γ I , γ Xe , the average nuclear concentrations of iodine and xenon Average neutron flux of the whole reactor fast group and thermal group and the power coefficient α in the reactor p The upper and lower limits of the search parameter range are set at most 0.01 times and 100 times the actual measured value in the reactor, respectively. The specific value range of the parameter needs to be adjusted based on the results of the model algorithm.

[0101] S4: Execute the proxy model's parameter search. After running the proxy model, it will perform a global parameter search based on the input ASI-t curve and the value range of the search parameter. The search results obtained by the input can be used to draw a curve similar to the original ASI-t curve of the corresponding working condition based on the formula in the model, achieving the desired effect.

[0102] The ASI-t calculation function of the proxy model is derived based on the reactor physics two-group one-dimensional diffusion model and the physical model of iodine and xenon concentration changes, as shown in Expression 3-6:

[0103]

[0104]

[0105] The proxy model is used for the average neutron flux of the full reactor fast group and thermal group with time and space variations. with the average nuclear concentrations of iodine and xenon The axial Fourier expansion and axial difference parameter are introduced, and their expressions are:

[0106]

[0107] E i =b i,2 (t)-b i,2 (∞)

[0108] Where i=1, 2, 3, 4 represent Four variables, n represents the order of the Fourier expansion coefficient, b i,2 (∞) represents the value of the second-order Fourier expansion coefficient of the variable under equilibrium conditions. E i Axial difference function representing the coefficients of the second-order Fourier expansion of the variable.

[0109] This leads to the derivation of the homogeneous equations for the four axial difference parameters:

[0110] 0=C 11 E1(t)+C 12 E2(t)

[0111] 0=C 21 E1(t)+C 22 E2(t)+C 24

[0112]

[0113] Where,

[0114]

[0115]

[0116] Finally, the relationship between the axial power shape index ASI and time is obtained through Laplace transform:

[0117] E P (t) = exp(Bt)[E P (0)cos(ωt)+qsin(ωt)]

[0118] Where, E P That is, the reactor axial power shape index ASI, and there is

[0119]

[0120]

[0121] Therefore, the ASI-t curve of the proxy model will be drawn based on formula (13).

[0122] Let x vector represent the combination of seven search parameter values, y j,1 is the input original ASI-t curve value, y j,2 (x) is the curve value obtained after inputting the seven search parameter values. MSE represents the mean square error between the input original ASI-t curve value and the calculated new curve value. Its expression is as follows:

[0123]

[0124] Where A represents the number of samples, which in the model represents the number of xenon transient steps.

[0125] The differential evolution-SLSQP hybrid algorithm of the surrogate model searches for the x vector to find the optimal x vector value that minimizes the MSE. The principle of the hybrid algorithm is as follows:

[0126] 1. Population initialization

[0127]

[0128] Among them, NP is the population number, which is 50 in this model and can be adjusted according to the model results. n with u n They represent the minimum and maximum values ​​in the range of the nth parameter in the x vector value.

[0129] 2. Mutation Strategy

[0130] After the population is initialized, the mutation strategy is used to mutate the population. The best / 1 mutation strategy is selected, and its expression is as follows:

[0131]

[0132] F(G)=0.5+1.0·e -0.02G

[0133] Among them, G is the maximum evolutionary generation, which is 150 in this model. is the value of the variant individual, is the optimal individual value, and F(G) is the scaling factor, which changes with the maximum evolutionary generations.

[0134] 3. Adaptive crossover operation

[0135] After mutation, the overall crossover between the mutation vector and the target vector is performed, and the expression is as follows:

[0136]

[0137] Among them, CR represents the crossover operator, which is 0.7 in this model. Represents the mth population, the nth search parameter, and the individual value of the G+1th population.

[0138] 4. Dynamic range shrinkage

[0139] During the population evolution process, the parameter range is constantly shrinking, and its expression is as follows:

[0140]

[0141] in, and Represents the average value of all populations under the nth search parameter of the Gth generation, l' n with u' n The upper and lower limits of the new search parameter value range respectively.

[0142] 5. Individual local correction

[0143] The optimal variable vector x obtained through the above process best Neighborhood constraints are imposed to obtain better results for each working condition. The constraints are as follows:

[0144]

[0145] st||xx best || ∞ ≤0.1(u' n -l′ n )

[0146] The constraint equation is solved using the BFGS quasi-Newton method.

[0147] After inputting the physical parameters and search parameter ranges under the corresponding working conditions, the model will perform global search and local correction of the parameters of the five working conditions based on the differential evolution-SLSQP hybrid algorithm. Under the action of the algorithm, the parameters that are insensitive to the model (usually γ I , γ Xe 、 ) will get the same result, achieving the effect of dimensionality reduction. In other words, these parameters corresponding to the five groups of working conditions only need the same value to reproduce the original curve. The hybrid algorithm performs sensitive parameters (i.e. α p ) is corrected locally to reduce the mean square error (MSE) and obtain search results that are more consistent with the original working condition data.

[0148] S5: Determine whether the relative standard deviation of the final output search parameters of the five groups of working conditions is within 5%. If so, proceed to the next step; if not, repeat the S4 process.

[0149] S6: The proxy model draws the ASI-t curve under the corresponding working condition based on the search parameter value finally output by S4, and performs spatiotemporal translation and scaling transformation based on the original ASI-t curve under the working condition.

[0150] The expression that defines the transformation function is as follows:

[0151]

[0152] Among them, Δt, Δy, and α represent the time transformation value of the curve (the default positive value is a left shift), the spatial transformation value (the default positive value is a left shift), and the scaling transformation value (with 1 as the base value, each value in the curve is multiplied by the scaling value to obtain a scaled new curve). is the ASI value after transformation.

[0153] The range of the three change values ​​is:

[0154] Δt∈[-2,2],|Δy|≤0.1max(|y1|),α∈[0.5,1.5]

[0155] The mean square error expression between the transformed curve value and the original curve value is:

[0156]

[0157] The initial value of the new curve obtained by the above transformation at t=0 is E P (0) and the effect of further reducing the MSE.

[0158] S7: Determine whether the improvement rate η of the MSE for each set of working conditions before and after the temporal and spatial translation correction and scaling correction is greater than 10%. The calculation formula for η is as follows:

[0159]

[0160] If yes, the time shift value Δt, the space shift value Δy and the scaling value α are output. If no, the time shift value Δt, the space shift value Δy are output as 0, and the scaling value α is output as 1.

[0161] S8: Execute the xenon oscillation process curve prediction part. In order to predict the relationship between the ASI and time under the xenon transient operating conditions with the same reactor power level suddenly dropping to other degrees, it is necessary to introduce a new variable power transient value ΔP to represent the magnitude of the reactor power level suddenly dropping under this condition (for example, if the reactor power level suddenly drops from 100% FP to 40% FP, the value of ΔP is 60). Then establish the power transient value ΔP and E Pd (0), α pThe relationship between ΔP and E Pd (0), α p and ΔP and The relationships between linear and exponential changes are approximated, so corresponding fitting equations can be established to calculate all relevant parameters. For the time shift Δt, spatial shift Δy, and scaling value α output by S7, linear fitting between two working conditions is used to make predictions.

[0162] S9: Input the ΔP value of the working condition to be predicted, and obtain the E of the working condition to be predicted through the fitting relationship obtained in S8 Pd (0), α p The predicted values ​​of the parameters and time shift Δt, space shift Δy and scaling α are input into the ASI calculation model to obtain the ASI-t curve under the predicted working conditions.

[0163] Example 1: Select the initial power level of 100% FP at the initial stage of HPR1000 reactor loading, and five groups of xenon transient operating conditions with power transients of ΔP of 30, 40, 60, 80, and 90. Calculate the ASI-t curves of the five groups of operating conditions and the corresponding initial physical parameters.

[0164] Table 1 Initial physical parameter values

[0165]

[0166]

[0167] Using the proxy model to perform parameter search, we get the following results:

[0168] Table 2 Final values ​​of search parameters

[0169]

[0170] Table 3 Transformation parameter values

[0171]

[0172] The fitting relationship obtained from Table 2 is as follows:

[0173] E Pd (0)=5.69526E-7-4.02001E-7*ΔP

[0174]

[0175] Δt, α, and Δy are obtained using linear interpolation based on the known power transient ΔP.

[0176] By inputting ΔP=35, 45, 55, 65, 75, 85 and calculating the corresponding parameters, the prediction curve of ASI-t is drawn as shown in Figure 3.

[0177] Example 2: Select the initial power level of 90% FP at the initial stage of HPR1000 reactor loading, and the power transient is ΔP of 30, 40, 60, 70, and 80. Use the traditional nuclear design program to calculate and obtain the ASI-t curves of the five conditions and the corresponding initial physical parameters, as shown in the following example: Figure 4 and as shown in Table 4.

[0178] Table 4 Initial physical parameter values

[0179] Initial physical parameters variable value unit H 410.8505 cm <![CDATA[λ I ]]> 2.92E-05 1 / s <![CDATA[λ Xe ]]> 2.10E-05 1 / s <![CDATA[σ Xe ]]> 8.28784E-19 <![CDATA[cm 2 ]]> <![CDATA[D1]]> 1.15244 <![CDATA[cm -1 ]]> <![CDATA[D2]]> 0.31457 <![CDATA[cm -1 ]]> <![CDATA[Σ a1 ]]> 0.00607 <![CDATA[cm -1 ]]> <![CDATA[Σ a2 ]]> 0.06231 <![CDATA[cm -1 ]]> <![CDATA[Σ R ]]> 0.06396 <![CDATA[cm -1 ]]> <![CDATA[Σ f1 ]]> 0.00139 <![CDATA[cm -1 ]]> <![CDATA[Σ f2 ]]> 0.02801 <![CDATA[cm -1 ]]> ν 1.51965 <![CDATA[E P (0)]]> 0.08961 <![CDATA[E Pd (0)(90%FP-10%FP)]]> -3.16905E-5 <![CDATA[E Pd (0)(90%FP-20%FP)]]> -2.74449E-5 <![CDATA[E Pd (0)(90%FP-30%FP)]]> -2.33321E-5 <![CDATA[E Pd (0)(90%FP-50%FP)]]> -1.54077E-5 <![CDATA[E Pd (0)(90%FP-60%FP)]]> -1.15401E-5

[0180] Using the proxy model to perform parameter search, we get the following results:

[0181] Table 5 Final values ​​of search parameters

[0182]

[0183]

[0184] Table 6 Transformation parameter values

[0185] Transformation results 90-10 (% FP) 90-20 (% FP) 90-30 (% FP) 90-50 (% FP) 90-60 (% FP) Δt 0.00E+00 0.00E+00 -2.47E-02 5.53E-01 1.13E+00 α 1.00E+00 1.00E+00 9.68E-01 8.84E-01 8.89E-01 Δy 0.00E+00 0.00E+00 7.79E-04 3.87E-02 5.51E-02

[0186] The fitting relationship obtained from Table 5 is as follows:

[0187] E Pd (0)=6.79992E-7-4.02487E-7*ΔP

[0188]

[0189] Δt, α, and Δy are obtained using linear interpolation based on the known power transient ΔP.

[0190] By inputting ΔP=25,35,45,55,65,75 and calculating the corresponding parameters, the prediction curve of ASI-t is drawn as shown in Figure 5 shown.

[0191] Those skilled in the art will understand that the above description is only a preferred embodiment of the present invention, and the features described in the various embodiments and / or claims of the present disclosure may be combined or coupled in various ways, even if such a combination or coupling is not explicitly described in the present disclosure. It is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art may still modify the technical solutions described in the aforementioned embodiments or make equivalent substitutions for some of the technical features therein. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention shall be included in the scope of protection of the present invention.

[0192] Although preferred embodiments of the present invention have been described, those skilled in the art may make additional changes and modifications to these embodiments once they are aware of the basic inventive concepts. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the present invention. Obviously, those skilled in the art may make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, the present invention is intended to include such changes and modifications as fall within the scope of the claims and their equivalents.

Claims

1. A nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints, characterized by: The method comprises the following steps: S1. Select five xenon transient operating conditions with the same reactor power level and different degrees of sudden drop. Set the duration of each xenon transient step to 1 hour. Calculate the change of the reactor axial power shape index (ASI) and the corresponding physical parameters over time for each operating condition under 38 xenon transient steps, and output the corresponding results. S2. Setting up a proxy model, inputting the change of the reactor axial power shape index ASI over time and the corresponding physical parameters for each group of operating conditions in S1 into the proxy model; The physical parameters include the reactor core height H, the decay constants of iodine and xenon nuclides, the diffusion coefficients D1 and D2 of the full reactor fast group and thermal group, and the absorption cross sections Σ a1 ,Σ a2 , fission cross section of full reactor fast group and hot group Σ f1 ,Σ f2 , the removal cross section of the whole stack Σ R , the number of neutrons produced per fission in the reactor ν and the reaction cross section σ of the nuclide xenon in the thermal group Xe ; and calculate the average value of the above parameters after averaging over the entire stack; S3. Define the search parameters for the yields of iodine and xenon γ I , γ Xe , the average nuclear concentrations of iodine and xenon Average neutron flux of the whole reactor fast group and thermal group and the power coefficient α in the reactor p The value range of S4. Execute the parameter search part of the proxy model. The proxy model performs a global search and dimensionality reduction of parameters for the five input working conditions, and performs a local correction of the search parameters for each input working condition. S5, judging whether the relative standard deviation of the final output search parameters of the five groups of working conditions is within 5%, if so, proceeding to S6; if not, repeating the process of S4; S6. The proxy model draws the ASI-t curve corresponding to each set of working conditions in combination with the search parameter value finally output by S4, and performs time-space translation and scaling transformation based on the original ASI-t curve under the working condition; S7. Determine whether the improvement rate η of the MSE for each group of working conditions before and after the temporal and spatial translation correction and the scaling correction is greater than 10%; S8, execute the xenon oscillation process curve prediction part of the agent model, combine the input working conditions and the output parameters to establish the power transient value ΔP and E Pd (0), α p and And the linear interpolation expressions of ΔP and the time translation value Δt, the spatial translation value Δy and the scaling value α between any two working conditions; S9, input the ΔP value of the working condition to be predicted, and obtain the E of the working condition to be predicted through the fitting relationship obtained in S8. Pd (0), α p The predicted values ​​of the parameters and the time translation value Δt, space translation value Δy and scaling value α are input into the ASI calculation model to obtain the ASI-t curve under the predicted working conditions.

2. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints according to claim 1, characterized in that: The change of reactor axial power shape index ASI with time and the corresponding physical parameters for each group of operating conditions in S1 are as follows: Where, P B With P T They represent the power values ​​of the lower and upper parts of the reactor axis respectively.

3. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints according to claim 1, characterized in that: S2 also includes the step of calculating the ASI value of the xenon transient working condition determination curve at time t=0, and the ASI value of the xenon transient working condition determination curve at time t=0 is named E P (0) and the derivative E Pd (0), the derivative value is the ASI change during the first xenon transient step, and its calculation expression is:

4. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints according to claim 1, characterized in that: The search parameters for iodine and xenon yields γ are defined in S3. I , γ Xe , the average nuclear concentrations of iodine and xenon Average neutron flux of the whole reactor fast group and thermal group and the power coefficient α in the reactor p The value range does not exceed 0.01 times and 100 times the actual measurement value inside the reactor.

5. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints according to claim 1, characterized in that: S4 also includes the proxy model for the average neutron flux of the full reactor fast group and thermal group for time- and space-varying variables with the average nuclear concentrations of iodine and xenon The steps of introducing the axial Fourier expansion and axial difference parameters, The expressions are: E i =b i,2 (t)-b i,2 (∞) Where i=1, 2, 3, 4 represent Four variables, n represents the order of the Fourier expansion coefficient, b i,2 (∞) represents the value of the second-order Fourier expansion coefficient of the variable under equilibrium conditions; E i Axial difference function representing the coefficients of the second-order Fourier expansion of the variable.

6. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints according to claim 1, characterized in that: The method for performing local correction of search parameters for each input working condition in S4 is: performing global search and local correction of the parameters of the working condition based on the differential evolution-SLSQP hybrid algorithm.

7. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints according to claim 1, characterized in that: The method for performing time-space translation and scaling transformation based on the original ASI-t curve under this working condition in S6 is: The expression that defines the transformation function is: Among them, Δt, Δy, and α represent the time transformation value of the curve respectively. The default positive value is the left translation and the spatial transformation value. The default positive value is the left translation and scaling transformation value. With 1 as the base value, each value in the curve is multiplied by the scaling value to obtain the scaled new curve. is the ASI value after transformation; Δt∈[-2,2], |Δy|≤0.1max(|y1|), α∈[0.5,1.5].

8. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints according to claim 1, characterized in that: In S7, the method for judging whether the improvement rate η of the MSE of each group of working conditions before and after the temporal and spatial translation transformation correction and the scaling correction is above 10% is as follows: If so, the corrected time shift value Δt, spatial shift value Δy and scaling value α are output; if not, the output results of the time shift value Δt and spatial shift value Δy are 0, and the scaling value α is 1.

9. A computer device comprising a memory and a processor, characterized in that A computer program is stored in the memory. When the processor runs the computer program stored in the memory, the processor executes the method according to any one of claims 1 to 8.

10. A computer-readable storage medium, characterized in that The computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 8 are implemented.

Images