A fast calculation method for nuclear reactor thermal hydraulic loop system based on model reduction

Abstract

The application discloses a kind of based on model reduction's nuclear reactor thermal loop system fast calculation method, using full-order model to solve full-order data, generate full-order matrix, solve POD base, reconstruct temperature field and velocity field using POD base, solve coefficient array, verify reconstruction solution using POD base;The accuracy difference between the reconstruction solution after reduction and the full-order data obtained by directly solving the full-order model can be controlled within 1%, the calculation degree of freedom changes from (variable number × grid number) to the sum of the POD order of variable, the calculation precision is high, and the calculation speed is fast;Under single working condition, the time required for existing commercial software to calculate full-order system loop model is usually about 2 minutes, and the time of reconstruction result by POD base after reduction is about 10 ms, and the acceleration ratio of the entire loop system model can reach 10 4 The order of magnitude;For a single component, the calculation time of existing full-order model is usually in the order of magnitude of 10s, and the calculation time of reconstruction solution by POD is about 2ms, and the acceleration ratio is about 10 3 The order of magnitude.

Description

Technical Field

[0001] This invention relates to the field of modeling and solution methods for nuclear reactor loop systems, and more particularly to a rapid calculation method for nuclear reactor thermal loop systems based on model order reduction. Background Technology

[0002] Currently, the nuclear reactors in operation in China are mainly third-generation reactors, represented by pressurized water reactors; while fourth-generation fast reactors are also developing rapidly in China, such as lead-bismuth fast reactors, sodium-cooled fast reactors, and high-temperature gas-cooled reactors.

[0003] In third-generation and fourth-generation reactors, the primary, secondary, and tertiary loop systems contain most of the reactor's components. The primary loop system mainly includes the reactor core, steam generator, main pump, and various types of pipes, valves, and chambers. The secondary loop system mainly consists of the secondary side of the steam generator, turbine, condenser, deaerator, and other facilities. The tertiary loop system includes cooling towers, cooling water pools, cooling water pumps, and associated coolant and cooling pipes.

[0004] Modeling the nuclear reactor loop system is crucial for simulating and predicting its operational behavior and performance, and is a vital step in reactor design, construction, and operation. The model involves knowledge of heat conduction, fluid mechanics, reaction kinetics, and control theory. By establishing a suitable model, the thermodynamic characteristics, safety, and economics of the reactor under different operating conditions can be evaluated. The modeling process can be achieved through mathematical equations and computational simulations.

[0005] The components included in the nuclear reactor loop system model are interconnected by pipes of different types and sizes, and the flow is controlled by many valves; in these modeling models, the governing equations are generally composed of the continuity equation, the momentum conservation equation, and the energy conservation equation.

[0006] For example, the continuity equation: Where ρ represents fluid density, t represents time, u represents fluid velocity along the x-axis, v represents fluid velocity along the y-axis, and w represents fluid velocity along the z-axis.

[0007] For example, the energy conservation equation: Where U is the velocity vector, η is the dynamic viscosity of the fluid, p is the pressure within the infinitesimal volume, and S u S v S w These are the generalized source terms of the three momentum equations.

[0008] For example, the momentum conservation equation:

[0009]

[0010] Where h is the specific enthalpy of the fluid, λ is the thermal conductivity, T is the temperature of the fluid, and Φ is the portion of mechanical energy converted into heat energy due to viscosity, called the dissipation function; S h It serves as the internal heat source for the fluid.

[0011] In the field of engineering, in order to improve the computational efficiency of the thermal system loop of large reactors, the above governing equations are often simplified into a one-dimensional form:

[0012] 1) Continuity equation:

[0013] 2) Energy conservation equation:

[0014] 3) Momentum conservation equation:

[0015] Full-order modeling is a numerical computation method that models the heat transfer process as mathematical equations and solves these equations using numerical methods. Commonly used numerical computation methods for heat transfer include the finite difference method, the finite element method, the finite volume method, and computational fluid dynamics. Full-order models can solve the aforementioned governing equations using these numerical computation methods.

[0016] Ideally, detectors or temperature measuring devices in the nuclear reactor loop would feed back data such as temperature and pressure to the system model, allowing the reactor loop system model to make accurate predictions about subsequent changes in the reactor based on real-time field data.

[0017] However, current models of nuclear reactor loop systems have large state space dimensions and numerical computational degrees of freedom, resulting in slow computation speed, low computational efficiency, and an inability to respond in real time to sudden changes in reactor operating conditions.

[0018] More importantly, the reactor primary loop system model can assess or simulate different accidents and anomalies, such as coolant leaks or fuel rod failures, using real-time field data to predict their impact on the reactor and take corresponding safety measures. This process can be summarized as: field data feedback - reactor primary loop system model calculation - providing field response. However, the actual situation is far from ideal. Although the calculation speed of the reactor primary loop system model has improved compared to the detailed three-dimensional models of individual components, the degree of freedom in numerical calculation is still very large, and its calculation speed is still relatively slow, making it difficult to achieve a good real-time response, especially for some transient situations. It is difficult to accurately predict the thermal evolution characteristics of the system and conduct accident analysis. If the calculation efficiency is improved by simplifying the model, such as reducing the number of grids, using simpler control equations, or using fitting formulas, then the simulation accuracy cannot be guaranteed.

[0019] In summary, current nuclear reactor system loop models often require several minutes or even tens of minutes to complete the full-order model calculation under transient and steady-state conditions. They cannot calculate the required accurate data instantly, and computational accuracy and efficiency cannot be guaranteed simultaneously. Summary of the Invention

[0020] To address the aforementioned technical problems, this invention provides a rapid calculation method for nuclear reactor thermal loop systems based on model order reduction, which offers high calculation accuracy and fast calculation speed.

[0021] The technical solution of this invention is as follows: A rapid calculation method for a nuclear reactor thermal loop system based on model order reduction, which first uses a full-order model to solve for full-order data, then generates a full-order matrix, then solves for the POD basis, then reconstructs the temperature and velocity fields using the POD basis, then solves for the coefficient array, and finally verifies the reconstructed solution using the POD basis; wherein, when solving for the coefficient array, the expressions for the corresponding POD basis reconstructed solutions for a single set of coefficient arrays are as follows:

[0022] Among them, a p,j,k b p,j,k Represents the variable value of each node reconstructed from the POD basis, j∈[1,NX];

[0023] Construct a residual function E that satisfies the governing equations. The convection term is represented by a first-order upwind scheme, and the diffusion term is represented by a center-split scheme. The expression for the corresponding residual function is as follows:

[0024]

[0025]

[0026] Among them, Ea E b This represents the residual function of variables a and b at the central node;

[0027] At the boundary nodes, the first type of boundary conditions are used for the inlet nodes, and the second type of boundary conditions are used for the outlet nodes.

[0028] The expressions for the residual functions of the momentum conservation equations at the inlet and outlet are:

[0029]

[0030]

[0031] The residual function expressions of the energy conservation equations at the inlet and outlet are:

[0032]

[0033]

[0034] The overall residual function E is expressed as follows:

[0035] E = E a +E b +E a,in +E a,out +E b,in +E b,out ;

[0036] There exists a unique set of C m,a With C m,b Make the residual function E after using the least squares method p The minimum value is the extreme point; solve for E. p The expression for the extreme point is:

[0037]

[0038]

[0039] Construct the equation LC = R, and the expressions for the left-hand matrix L and the right-hand matrix C are as follows:

[0040]

[0041]

[0042] in, l represents an element within L; in the process of solving for the minimum residual, it will be related to C. m Place irrelevant terms on the right-hand side of the equation, which will be related to C. m The relevant terms are placed on the left side of the equation.

[0043] The aforementioned rapid calculation method for nuclear reactor thermal loop systems based on model order reduction includes the following: When using a full-order model to solve for full-order data, the fluid modeled in the full-order model is a compressible fluid, with flow velocity as a parameter, and its governing equations are as follows:

[0044] Continuity equation:

[0045] Energy conservation equation:

[0046] Momentum conservation equation:

[0047] Where T is the fluid temperature, q is the output heat, P is the fluid pressure, and S is the source term; the fluid density ρ, the fluid velocity u along the x-axis, and the fluid's specific heat capacity at constant pressure C are selected. p Using the fluid properties, such as the hydrodynamic viscosity μ, as parameters, parameter variations under different operating conditions are set; the numerical heat transfer discretization scheme is selected; temperature and pressure are selected as variables, and the temperature and pressure fields under different flow velocities are solved using numerical algorithms.

[0048] The aforementioned method for rapid calculation of nuclear reactor thermal loop systems based on model order reduction includes a method where, when generating the full-order matrix, the full-order matrix of the temperature field is A. Where m is the number of grid cells and n is the number of flow velocity values.

[0049] A = [α1, α2, ..., α] k , …, α n ],

[0050] Where, α k Indicates the flow velocity as u k The velocity field distribution at time k∈[1,n]

[0051] α k =[a 1,k ,…,a h,k ,…,a NX,k ] T ,

[0052] Among them, a h,k This represents the value of variable a in the h-th grid when the parameter takes the k-th value, where h∈[1,NX].

[0053] Similarly, the full-order matrix of the pressure field is B. Where m is the number of grid cells and n is the number of flow velocity values.

[0054] B = [β1, β2, ..., β] k , ..., β n ],

[0055] Where, β k Indicates the flow velocity as u k Temperature field distribution at time k∈[1,n]

[0056] β k =[b 1,k ,…,b h,k ,…,b NX,k ] T ,

[0057] Among them, b h,k This represents the value of variable b in the h-th grid when the parameter takes the k-th value, where h∈[1,NX].

[0058] The aforementioned method for rapid calculation of nuclear reactor thermal loop systems based on model order reduction includes the following step: when generating the full-order matrix, the full-order matrix A of the temperature field and the full-order matrix B of the pressure field are merged to generate a new matrix Z as the new full-order matrix.

[0059] The aforementioned fast calculation method for nuclear reactor thermal loop systems based on model order reduction, wherein: when solving for the POD basis, the POD basis is related to the square matrix S a or S b The relationship between eigenvalues ​​and eigenvectors:

[0060] S a =A T A,

[0061] S b =B T B;

[0062] Let {λ1,λ2,…,λ} N} is a square matrix S a or S b The eigenvectors are such that λ1>λ2>…>λ N If the value is greater than 0, then its corresponding eigenvector is {X1,X2,…,X...} N}, corresponding to POD base P m The expression is:

[0063]

[0064]

[0065] The corresponding POD order N p The expression is:

[0066]

[0067]

[0068] Where I(NP) represents the energy percentage of the POD base.

[0069] The aforementioned fast calculation method for nuclear reactor thermal loop systems based on model order reduction, wherein: when solving for the POD basis, the matrix S is used... z The relationship between eigenvalues ​​and eigenvectors is S z =Z T Z, corresponding to POD base P m The expression is:

[0070]

[0071] POD order N p The expression is:

[0072]

[0073] The fast calculation method for nuclear reactor thermal loop system based on model order reduction, wherein the expression for the variable values ​​of each node reconstructed from the POD basis when reconstructing the temperature and velocity fields is as follows:

[0074]

[0075]

[0076] Where C is the coefficient array, N p Here, m represents the coefficient number in the coefficient array, i represents the grid number, and k represents the parameter number.

[0077] The aforementioned rapid calculation method for nuclear reactor thermal loop systems based on model order reduction, wherein: when reconstructing the temperature and velocity fields, POD-based P m Divided into P m,a With P m,b , and [P m,1 ,…P m,NX ,P m,NX+1 ,…,P m,2NX ] T =[P m,a,1 ,…P m,a,NX ,P m,b,1 ,…,P m,b,NX ] T .

[0078] The fast calculation method for nuclear reactor thermal loop system based on model order reduction is described above, wherein: when solving the coefficient array, the diffusion term adopts the Quick format instead of the center split format.

[0079] The aforementioned fast calculation method for nuclear reactor thermal loop systems based on model order reduction, wherein: when verifying the reconstructed solution, any number within the range of u is arbitrarily chosen, u = x, x ∈ [min(u k ), max(u k The full-order model is used to calculate the full-order data; and the reconstructed solution is calculated using the POD basis. Residuals are then calculated for the values ​​of each variable at each node.

[0080]

[0081]

[0082] And take the maximum value among them:

[0083] Compared with existing technologies, the fast calculation method for nuclear reactor thermal loop systems based on model order reduction provided by this invention can control the accuracy difference between the solution reconstructed using POD basis and the full-order data obtained directly from the full-order model to within 1% after order reduction. The computational degrees of freedom change from (number of variables × number of grids) to the sum of the POD orders of the variables. This not only results in high computational accuracy but also high computational speed. Under a single operating condition, existing commercial software typically takes about 2 minutes to calculate the full-order system loop model, while after POD order reduction, the time for POD basis reconstruction is about 10 ms, achieving a speedup of 10 for the entire loop system model. 4 The computation time for a single component is on the order of magnitude of 10 seconds for existing full-order models, while the computation time for reconstructing the solution from the POD is around 2 milliseconds, resulting in a speedup of approximately 10^5 seconds. 3 The order of magnitude. Attached Figure Description

[0084] The accompanying drawings described herein are for illustrative purposes only and are not intended to limit the scope of the invention in any way; the shapes and proportions of the components in the drawings are merely illustrative and are intended to aid in understanding the invention, and are not intended to specifically limit the shapes and proportions of the components of the invention; those skilled in the art, under the guidance of this invention, can select various possible shapes and proportions to implement the invention according to specific circumstances.

[0085] Figure 1 This is a flowchart of the rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to the present invention;

[0086] Figure 2 This is a comparison chart of the full-order model calculation data and the basic reconstruction solution curve of POD after the POD order reduction of a single heat exchange pipeline according to the present invention;

[0087] Figure 3This is an error curve of temperature at different nodes after POD reduction of a single heat exchange pipe according to the present invention;

[0088] Figure 4 This is a comparison chart of the calculation data of the full-order model of variable a after POD order reduction for a single heat transfer component and the POD basis reconstruction solution curve.

[0089] Figure 5 This is a comparison chart of the calculation data of the full-order model of variable b after POD order reduction for a single heat transfer component and the POD basis reconstruction solution curve.

[0090] Figure 6 This is a graph showing the error between the POD reconstruction solution of variable a at each grid and the full-order data calculated by the full-order model after POD order reduction for a single heat transfer component according to the present invention.

[0091] Figure 7 This is a graph showing the error between the POD reconstruction solution and the full-order data calculated by the full-order model at each grid after the POD order reduction of a single heat transfer component according to the present invention. Detailed Implementation

[0092] The specific embodiments and examples of the present invention will be described in detail below with reference to the accompanying drawings. The specific embodiments described are only used to explain the present invention and are not intended to limit the specific embodiments of the present invention.

[0093] Before performing Proper Orthogonal Decomposition (POD) reduction on a nuclear reactor loop system model, it is necessary to isolate each component in the nuclear reactor loop system, know the governing equations of each component, and understand the data transfer between the various components and pipeline models. Then, POD reduction calculations are performed on each component and pipeline separately, ensuring that the transmission interfaces and parameters are consistent with the full-order model. The specific steps and operations for reducing the order of a component in the nuclear reactor loop system model can be performed either by using commercial software to model and numerically solve the full-order model, or by using self-programmed languages ​​such as C++ to numerically solve the full-order model.

[0094] like Figure 1 As shown, the present invention provides a rapid calculation method for a nuclear reactor thermal loop system based on model order reduction, comprising the following steps:

[0095] Step S110: Solve for the full-order data using the full-order model;

[0096] Step S120: Generate a full-order matrix;

[0097] Step S130: Solve for the POD basis;

[0098] Step S140: Reconstruct the temperature and velocity fields using the POD basis;

[0099] Step S150: Solve for the coefficient array;

[0100] Step S160: Verify the reconstructed solution using the POD basis.

[0101] After order reduction, the accuracy difference between the solution reconstructed using POD basis and the full-order data obtained directly from the full-order model can be controlled within 1%. The computational degrees of freedom change from (number of variables × number of grids) to the sum of the POD orders of the variables. At the same time, the time to reconstruct the full-order data under a single working condition will also remain on the order of about 10ms.

[0102] Taking the one-dimensional steady-state equation as an example, in step S110, the fluid in the full-order model is a compressible fluid, and the flow velocity is a parameter. Its governing equations are as follows:

[0103] Continuity equation:

[0104] Energy conservation equation:

[0105] Momentum conservation equation:

[0106] Where T is the fluid temperature, q is the output heat, P is the fluid pressure, and S is the source term, including frictional pressure drop and gravitational pressure drop terms; a variable can be selected as the parameter, such as fluid density ρ, fluid velocity along the x-axis u, and fluid specific heat capacity at constant pressure C. p Fluid properties such as fluid dynamic viscosity μ are used as parameters, and parameter variations are set under different operating conditions; assuming flow velocity is a parameter, the flow velocity takes a total of n values, with the value being u. k k∈(1,…,n), the number of grids is NX; any commonly used numerical heat transfer discretization scheme can be selected; assuming that the convection term is selected as the first-order upwind scheme and the diffusion term is selected as the center split scheme, temperature and pressure are selected as variables, and then the temperature field and pressure field under different flow velocities are solved by numerical algorithm, and these temperature field and pressure field under different flow velocities are the full-order data.

[0107] The so-called first-order upwind scheme refers to a numerical method commonly used to solve partial differential equations. It is based on Taylor series expansion and forward difference to discretize time and space, thereby obtaining a numerical solution. The characteristics of the first-order upwind scheme are simple computation, but low accuracy, and it is prone to numerical dissipation and instability when dealing with problems with strong convection terms. Therefore, in practical applications, the first-order upwind scheme is usually replaced by higher-order numerical methods.

[0108] In step S120, the full-order matrix of the temperature field is A. in, The notation is for a matrix (the same applies below), where m×n represents m rows and n columns (the same applies below), m is the number of grid cells, and n is the number of flow velocity values.

[0109] A = [α1, α2, ..., α] k , …, α n ],

[0110] Where, α k Indicates the flow velocity as u k The velocity field distribution at time k∈[1,n]

[0111] α k =[a 1,k ,…,a h,k ,…,a NX,k ] T ,

[0112] Among them, a h,k This represents the value of variable a in the h-th grid when the parameter takes the k-th value, where h∈[1,NX].

[0113] Similarly, the full-order matrix of the pressure field is B. Where m is the number of grid cells and n is the number of flow velocity values.

[0114] B = [β1, β2, ..., β] k , ..., β n ],

[0115] Where, β k Indicates the flow velocity as u k Temperature field distribution at time k∈[1,n]

[0116] β k =[b 1,k ,…,b h,k ,…,b NX,k ] T ,

[0117] Among them, b h,k This represents the value of variable b in the h-th grid when the parameter takes the k-th value, where h∈[1,NX].

[0118] When solving for the POD basis in step S130,

[0119] According to POD construction theory, the definition of a POD basis is related to the square matrix S. a or S b The relationship between eigenvalues ​​and eigenvectors:

[0120] S a =AT A,

[0121] S b =B T B;

[0122] If {λ1,λ2,…,λ} N} is a square matrix S a or S b The eigenvectors are such that λ1>λ2>…>λ N If the value is greater than 0, then its corresponding eigenvector is {X1,X2,…,X...} N Therefore, the corresponding POD base P m :

[0123]

[0124]

[0125] According to POD theory, the magnitude of an eigenvalue represents the importance of the corresponding eigenvector; in fact, the eigenvectors of first-order eigenvalues ​​have a critical proportion to the distribution of both temperature and pressure fields, thus allowing the determination of the POD order N. p :

[0126]

[0127]

[0128] Among them, I(N) P This is called the energy percentage of the POD base, which can be customized according to different full-order matrices.

[0129] In step S140, when reconstructing the temperature and velocity fields, the expression for the variable values ​​of each node reconstructed from the POD basis is as follows:

[0130]

[0131]

[0132] Where C is the coefficient array, N p The dimension is denoted by m, where m represents the coefficient number in the coefficient array, i represents the grid number, and k represents the parameter number. Different POD bases correspond to different coefficient matrices, and their dimensions are also different.

[0133] In step S150, during the process of solving the coefficient array, for the solution of a single set of coefficient arrays, the difference between the full-order data calculated by the POD basis reconstruction solution and the full-order model can be minimized. This is calculated using a mathematical method for finding the minimum value. Assume the POD basis reconstruction solutions are as follows:

[0134] Among them, a p,j,k b p,j,k Represents the variable value of each node reconstructed from the POD basis, j∈[1,NX];

[0135] Assuming the residual function is E, then E needs to satisfy the governing equations. When constructing the residual function, the convection term is selected using the first-order upwind scheme, and the diffusion term adopts the center-split scheme. The calculation process of the corresponding residual function is as follows:

[0136]

[0137]

[0138] Among them, E a E b The sum of the squared residuals of the full-order data of variables a and b at the central node and the POD basis reconstruction solution is called the residual function of the central node.

[0139] At the boundary nodes, the residuals can be calculated by substituting into the governing equations for different types of boundary conditions. For example, the first type of boundary conditions are used for the inlet nodes, while the second type of boundary conditions are used for the outlet nodes.

[0140] The residual functions of the momentum conservation equations at the inlet and outlet are expressed as follows:

[0141]

[0142]

[0143] The residual function expressions of the energy conservation equations at the inlet and outlet are:

[0144]

[0145]

[0146] The overall residual function E is expressed as follows:

[0147] E = E a +E b +E a,in +E a,out +E b,in +E b,out ;

[0148] To minimize E, simply make E p According to the least squares method, there must exist a unique set of C. m,a With C m,b Make E p Minimum,E p Let be the residual function after applying the least squares method, and let the minimum value be the extreme point;

[0149] Solve for E p The expression for the extreme point is:

[0150]

[0151]

[0152] For E, that is, there exists (N) p,a +N p,b There are ) equations, and C in each equation m Each has a corresponding coefficient; through numerous numerical examples and mathematical derivation, this extreme point is the minimum point, thus the equation LC = R can be constructed, where the expressions for the left matrix L and the right matrix C are respectively:

[0153]

[0154]

[0155] in, l represents an element within L; in the process of solving for the minimum residual, it will be related to C. m Place irrelevant terms on the right-hand side of the equation, which will be related to C. m The relevant terms are placed on the left side of the equation.

[0156] In step S160, when verifying the reconstructed solution, any number within the range of u is chosen, u = x, x ∈ [min(u k ), max(u k The full-order model is used to calculate the full-order data; and the reconstructed solution is calculated using the POD basis. Residuals are then calculated for the values ​​of each variable at each node.

[0157]

[0158]

[0159] And take the maximum value among them:

[0160] Under normal circumstances, E m,a With E m,b The difference should be less than 1%.

[0161] Compared to existing technologies, current commercial software typically takes about 2 minutes to compute a full-order system loop model under a single operating condition. However, after performing POD order reduction, the time to reconstruct the result from the POD basis is about 10 ms, resulting in a speedup of 10 for the entire loop system model. 4The computation time for a single component is on the order of magnitude of 10 seconds for existing full-order models, but the computation time for reconstructing the solution using POD is around 2 milliseconds, resulting in a speedup of approximately 10^5 seconds. 3 The order of magnitude.

[0162] For example, for a single heat exchanger pipe, after POD order reduction, the comparison results between its full-order model calculation data and the POD basic reconstruction solution curve are as follows: Figure 2 As shown, the temperature error results at different nodes are as follows: Figure 3 As shown, after the nuclear reactor loop system model is reduced in order by POD, its state space dimension will generally be reduced to about 20 dimensions. For a single component, the calculation time for reconstructing the solution by POD is about 2ms, and the error can be basically controlled to below 1%.

[0163] The rapid calculation method for nuclear reactor thermal loop systems based on model order reduction of the present invention can also be modified for steps S120, S130 and S140 in the foregoing embodiments as follows:

[0164] In step S120, after obtaining the full-order matrix A of the temperature field and the full-order matrix B of the pressure field, the two are merged to generate a new matrix Z as the new full-order matrix.

[0165] When solving for the POD basis in step S130, the matrix S is used. z The relationship between eigenvalues ​​and eigenvectors then becomes S z =Z T Z, and the corresponding POD base P m The expression then becomes:

[0166]

[0167] Therefore, the order N of POD p The expression can be represented as:

[0168]

[0169] During step S140, when reconstructing the temperature and velocity fields, the POD basis P at this time m Divided into P m,a With P m,b And: [P m,1 ,…P m,NX ,P m,NX+1 ,…,P m,2NX ] T =[P m,a,1 ,…P m,a,NX ,P m,b,1 ,…,P m,b,NX ] T;

[0170] The expression for the variable value of each node reconstructed from the POD basis remains unchanged, but the coefficient matrix corresponding to different POD bases is the same.

[0171] Furthermore, in the rapid calculation method for nuclear reactor thermal loop systems based on model order reduction of this invention, when solving the coefficient array in step S150 of the aforementioned embodiment, the same discretization scheme as the full-order model is used when constructing the residual function E. In most cases, for loop system models built by some commercial software, only the form of the governing equations is known, and the specific discretization scheme of each equation is unknown. However, in POD order reduction, step S150 requires equation discretization when constructing the residual function, which requires the use of a discretization scheme. If a different discretization scheme than the full-order model is used, for example, the first-order upwind scheme is used for the convection term in the full-order model, while the Quick scheme is used for the diffusion term when constructing the residual function, the error of the obtained result will not be too large. For specific examples, please refer to Embodiment 2 below.

[0172] Example 1: Rapid Calculation of POD Order Reduction for Steady-State Model of Primary Side of Spiral Steam Generator

[0173] Taking the helical tube steam generator in the primary loop system of the fourth-generation fast reactor lead-bismuth fast reactor as an example, its fluid is compressible fluid helium gas, and the governing equations of the primary side steady-state model of the helical tube steam generator are the momentum conservation equation and the energy conservation equation.

[0174] Momentum conservation equation:

[0175] In the formula, ρ is the density of helium (kg / m³). 3 ), u is the flow velocity of helium (m / s), P is the pressure of helium (Pa), μ is the dynamic viscosity of helium (Pa·s), and g is the acceleration due to gravity (m / s²). 2 ), where θ is the inclination angle of the helical tube, f is the friction coefficient, and d is the diameter of the helical tube (m).

[0176] Energy conservation equation:

[0177] In the formula, C p The specific heat capacity of helium at constant pressure (J / (kg·m³)) 3 T is the temperature of helium (K), λ is the thermal conductivity, and q is the heat output to the secondary side (W / m). 3 ).

[0178] For helium, the specific physical properties selected are as follows:

[0179]

[0180] C p =5194

[0181] μ = 3.674·10 -7 ·T 0.7

[0182]

[0183] The grid number is set to NX, and the helium flow rate is set as a parameter with 11 possible values ​​as follows:

[0184] As shown in Table 1:

[0185] 1 -2.0 2 -2.6 3 -3.2 4 -3.8 5 -4.4 6 -5.0 7 -5.6 8 -6.2 9 -6.8 10 -7.4 11 -8.0

[0186] Table 1

[0187] The corresponding variables and their corresponding codes are set as shown in Table 2 below:

[0188] ρT a λT b ρ c μ d P e

[0189] Table 2

[0190] Temperature and pressure were selected as variables, and then numerical algorithms were used to solve for the temperature and pressure fields of helium at different flow rates. Full-order data were calculated, and the full-order data were used to generate a full-order matrix M.

[0191] M = [α1, α2, ..., α] k , …, α 11 ];

[0192] Where M∈m×n, m=5*NX, n=11; α k Indicates the flow velocity as u k Velocity field distribution at time,

[0193] α k =[a 1,k ,…,a h,k ,…,a NX,k ,b 1,k ,…,b h,k ,…,b NX,k ,c 1,k ,…,c h,k ,…,c NX,k ,k∈[1,11],d 1,k ,…,d h,k ,…,d NX,k ,e 1,k ,…,e h,k ,…,e NX,k ] T ;

[0194] Among them, a h,k This represents the value of variable 'a' in the h-th grid when the parameter takes the k-th value; b h,k This represents the value of variable b in the h-th grid when the parameter takes the k-th value; c h,k This represents the value of variable c in the h-th grid when the parameter takes the k-th value; d h,k This represents the value of variable d in the h-th grid when the parameter takes the k-th value; e h,k This represents the value of variable e in the h-th grid when the parameter takes the k-th value.

[0195] Next, following the POD basis solution steps in the aforementioned scheme, the POD basis of the full-order matrix M is obtained:

[0196] S m =M T M;

[0197] The POD basis of the full-order matrix M is the same as that of the square matrix S. m The relationship between eigenvalues ​​and eigenvectors; if {λ1,λ2,…,λ... N} is a square matrix S m The eigenvalues ​​satisfy λ1>λ2>…>λ N If the value is greater than 0, its corresponding eigenvector is {X1,X2,…,X}. N}, then the POD basis is defined as:

[0198]

[0199] Set the critical percentage value and solve for the eigenvalues ​​and eigenvectors of the critical percentage to determine the order N of the POD. p :

[0200]

[0201] The above variables are then reconstructed using the POD basis:

[0202]

[0203]

[0204]

[0205]

[0206]

[0207] Next, the coefficient array C is solved to minimize the residual function E. Here, the same coefficient matrix is ​​used before different variables. Since the flow velocity is negative, the inlet node is located at number NX and the outlet node is located at number 1.

[0208] The residual functions of the momentum conservation equation and the residual functions on the boundary are expressed as follows:

[0209]

[0210]

[0211]

[0212] The residual functions of the energy conservation equation and the residual functions on the boundary are expressed as follows:

[0213]

[0214]

[0215]

[0216] The overall residual function E is expressed as follows:

[0217] E = E 1,in +E 1,out +E1+E 2,in +E 2,out +E2.

[0218] The coefficient array C that minimizes the residual function is obtained using the least squares method; then, u = 5.6 is selected, the corresponding coefficient matrix is ​​calculated, and the POD reconstruction solution is obtained. The result is then calculated using a full-order model; finally, the two sets of data are compared. The calculation formula is as follows:

[0219]

[0220]

[0221]

[0222]

[0223]

[0224]

[0225] E i,k,x This represents the error of variable x in the i-th grid when the parameter takes the k-th value; the average error and maximum error of variable x are shown in Table 3 below:

[0226] a <![CDATA[7.9531×10 -6 ]]> <![CDATA[2.2979×10 -5 ]]> b 0.0500 0.1304 c 0.0435 0.1092 d 0.0085 0.0318 e <![CDATA[7.2316×10 -5 ]]> <![CDATA[1.4344×10 -4 ]]>

[0227] Table 3

[0228] The maximum error between the full-order model calculation and the reduced-order model calculation in Table 3 is 0.1304%. The degrees of freedom of the full-order model are 5*NX, and the degrees of freedom of the reduced-order model are N. p The computational efficiency has been significantly improved; it has been verified that the computation time for the full-order model is 14 seconds, and the computation time for the POD basis reconstruction solution is 2 ms, achieving a speedup of 10. 4 Order of magnitude.

[0229] Example 2: Error Comparison of Full-Order Model and POD Reduced-Order Model in a Single Heat Transfer Component Using Different Discretization Schemes

[0230] For a single heat transfer component, the energy conservation equation is:

[0231]

[0232] In the full-order model calculation, the convection term uses the Quick scheme, and the diffusion term uses the central split scheme. After integrating and discretizing the governing equations, we can obtain:

[0233]

[0234] Its physical properties are as follows:

[0235]

[0236]

[0237] The number of grids is set to NX, and the flow rate is set as a parameter. The 11 flow rate values ​​are shown in Table 4 below:

[0238] 1 -2.0 2 -2.6 3 -3.2 4 -3.8 5 -4.4 6 -5.0 7 -5.6 8 -6.2 9 -6.8 10 -7.4 11 -8.0

[0239] Table 4 (same as Table 1)

[0240] And ρT and λT are treated as variables, with the corresponding symbols shown in Table 5 below:

[0241] ρT a λT b

[0242] Table 5

[0243] The inlet boundary is subject to boundary conditions of the first kind, and the outlet boundary is subject to boundary conditions of the second kind. The finite volume method is used to numerically solve the full-order model, calculate the full-order data, and generate a full-order matrix M from the full-order data.

[0244] M = [α1, α2, ... α] k , …, α 11 ];

[0245] Where M∈m×n, m=2*NX, n=11; α k Indicates the flow velocity as u k The velocity field distribution at time k∈[1,11], α k =[a 1,k ,…,a h,k ,…,a NX,k ,b 1,k ,…,b h,k ,…,b NX,k ] T ;

[0246] Among them, a h,k This represents the value of variable 'a' in the h-th grid when the parameter takes the k-th value; b h,k This represents the value of variable b in the h-th grid when the parameter takes the k-th value.

[0247] Next, following the POD basis solution steps in the aforementioned scheme, the POD basis of the full-order matrix M is obtained:

[0248] S m =M T M;

[0249] The POD basis of the full-order matrix M is the same as that of the square matrix S. m The relationship between eigenvalues ​​and eigenvectors; if {λ1,λ2,…,λ... N} is a square matrix S m The eigenvalues ​​satisfy λ1>λ2>…>λ N If the value is greater than 0, its corresponding eigenvector is {X1,X2,…,X}. N}, then the POD basis is defined as:

[0250]

[0251] Set the critical percentage value and solve for the eigenvalues ​​and eigenvectors of the critical percentage to determine the order N of the POD. p :

[0252]

[0253] The above variables are then reconstructed using the POD basis:

[0254]

[0255]

[0256] Next, the coefficient array C is solved to minimize E. In this residual matrix, the discretization scheme of the convection term at the center node adopts the Quick scheme, the convection term at the second inlet node adopts the first-order upwind scheme, the boundary conditions at the inlet node are first-type boundary conditions, and the boundary conditions at the outlet node are second-type boundary conditions. The corresponding residual function expressions are as follows.

[0257] The residual function expression for the first node at the inlet is:

[0258]

[0259] The residual function expression for the second node at the inlet is:

[0260]

[0261] The residual function expression for the central node is:

[0262]

[0263] The residual function expression for the exit node is:

[0264]

[0265] The overall residual function expression is:

[0266] E = E in,1 +E in,2 +E medium +E out .

[0267] The least squares method is used to find the residual array C that minimizes the coefficient array; then, u = 5.6 is selected, the corresponding coefficient matrix is ​​calculated, and the POD reconstruction solution is obtained. The full-order result is then calculated using the full-order model; finally, the two sets of data are compared. The calculation formula is as follows:

[0268]

[0269]

[0270] E i,k,x This represents the error of variable x in the i-th grid when the parameter takes the k-th value; a comparison of the full-order model calculation data of variable a with the POD basis reconstruction solution curve is shown below. Figure 4 As shown, the comparison between the full-order model calculation data of variable b and the POD basis reconstruction solution curve is as follows: Figure 5 As shown; the error curves between the POD reconstruction solution and the full-order data calculated by the full-order model at each grid for variable a are as follows. Figure 6 As shown, the maximum error is 1.6680 × 10⁻⁶. -4%; The error curves between the POD reconstruction solution and the full-order data calculated by the full-order model at each grid for variable b are shown below. Figure 7 As shown, the maximum error is 0.5834%.

[0271] The results above show that in the current example, when the discrete scheme used for the numerical calculation of the full-order model is different from the discrete scheme used in the process of reconstructing the solution using the POD basis, the maximum error is 0.5834%, which is within an acceptable range.

[0272] After the nuclear reactor loop system model is reduced in order by POD, its state space dimension will generally be reduced to about 20 dimensions. For a single component, the calculation time for reconstructing the solution by POD is about 2ms, and the error can be basically controlled to below 1%.

[0273] It should be noted that the rapid calculation method for nuclear reactor thermal loop system based on model reduction of the present invention can simultaneously ensure calculation accuracy and calculation efficiency. It does not implement simple functions with overly complicated steps, nor does it use conventional or simple features for combination or stacking. Therefore, it conforms to the common sense of technical improvement and has practicality.

[0274] Any content not described in detail in this specification is prior art known to those skilled in the art.

[0275] It should be understood that the above description is only a preferred embodiment of the present invention and is not sufficient to limit the technical solution of the present invention. For those skilled in the art, within the spirit and principles of the present invention, additions, subtractions, substitutions, transformations or improvements can be made based on the above description, and all such additions, subtractions, substitutions or improvements should fall within the protection scope of the appended claims of the present invention.

Claims

1. A rapid calculation method for nuclear reactor thermal loop systems based on model order reduction, characterized in that, First, the full-order model is used to solve the full-order data, then the full-order matrix is ​​generated. Next, the POD basis is solved, and then the temperature and velocity fields are reconstructed using the POD basis. Then, the coefficient array is solved, and finally, the reconstructed solution is verified using the POD basis. Specifically, when solving the coefficient array, the expressions for the corresponding POD basis reconstructed solution for a single set of coefficients are as follows: A p =[a p,1,k ,…,a p,NX,k ] T B p =[b p,1,k ,…,b p,NX,k ] T , where a p,j,k b p,j,k Represents the variable value of each node reconstructed from the POD basis, j∈[1,NX]; Construct a residual function E that satisfies the governing equations. The convection term is represented by a first-order upwind scheme, and the diffusion term is represented by a center-split scheme. The expression for the corresponding residual function is as follows: Among them, E a E b This represents the residual function of variables a and b at the central node; At the boundary nodes, the first type of boundary conditions are used for the inlet nodes, and the second type of boundary conditions are used for the outlet nodes. The expressions for the residual functions of the momentum conservation equations at the inlet and outlet are: The residual function expressions of the energy conservation equations at the inlet and outlet are: The overall residual function E is expressed as follows: E=E a +E b +E a,in +E a,out +E b,in +E b,out ; There exists a unique set of C m,a With C m,b Make the residual function E after using the least squares method p The minimum value is the extreme point; solve for E. p The expression for the extreme point is: Construct the equation LC = R, and the expressions for the left-hand matrix L and the right-hand matrix C are as follows: in, l represents an element within L; in the process of solving for the minimum residual, it will be related to C. m Place irrelevant terms on the right-hand side of the equation, which will be related to C. m The relevant terms are placed on the left side of the equation.

2. The rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to claim 1, characterized in that: When using a full-order model to solve full-order data, the fluid modeled is a compressible fluid with velocity as a parameter, and its governing equations are as follows: Continuity equation: Energy conservation equation: Momentum conservation equation: Where T is the fluid temperature, q is the output heat, P is the fluid pressure, and S is the source term; the fluid density ρ, the fluid velocity u along the x-axis, and the fluid's specific heat capacity at constant pressure C are selected. p Using the fluid properties, such as the hydrodynamic viscosity μ, as parameters, parameter variations under different operating conditions are set; the numerical heat transfer discretization scheme is selected; temperature and pressure are selected as variables, and the temperature and pressure fields under different flow velocities are solved using numerical algorithms.

3. The rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to claim 1, characterized in that: When generating the full-order matrix, the full-order matrix of the temperature field is A. Where m is the number of grid cells and n is the number of flow velocity values. A＝[α1，α2，…，α k ，…，α n ]， Where, α k Indicates the flow velocity as u k The velocity field distribution at time k∈[1,n] a k =[a 1,k ,…,a h,k ,…,a NX,k ] T , Among them, a h,k This represents the value of variable a in the h-th grid when the parameter takes the k-th value, where h∈[1,NX]. Similarly, the full-order matrix of the pressure field is B. Where m is the number of grid cells and n is the number of flow velocity values. B＝[β1，β2，…，β k ，…，β n ], Where, β k Indicates the flow velocity as u k Temperature field distribution at time k∈[1,n] b k =[b 1,k ,…,b h,k ,…,b NX,k ] T , Among them, b h,k This represents the value of variable b in the h-th grid when the parameter takes the k-th value, where h∈[1,NX].

4. The rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to claim 1, characterized in that: When generating the full-order matrix, the full-order matrix A of the temperature field and the full-order matrix B of the pressure field are merged to generate a new matrix Z as the new full-order matrix.

5. The rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to claim 1, characterized in that: When solving for the POD basis, the POD basis is related to the square matrix S. a or S b The relationship between eigenvalues ​​and eigenvectors: S a =A T A, S b =B T B; Let {λ1,λ2,…,λ} N } is a square matrix S a or S b eigenvectors and λ1>λ2>…>λ N If the value is greater than 0, then its corresponding eigenvector is {X1,X2,…,X...} N }, corresponding to POD base P m The expression is: The corresponding POD order N p The expression is: Where I(NP) represents the energy percentage of the POD base.

6. The rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to claim 1, characterized in that: When solving for the POD basis, with respect to the square matrix S z The relationship between eigenvalues ​​and eigenvectors is S z =Z T Z, corresponding to POD base P m The expression is: POD order N p The expression is:

7. The rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to claim 1, characterized in that: When reconstructing the temperature and velocity fields, the expression for the variable values ​​of each node reconstructed from the POD basis is as follows: Where C is the coefficient array, N p Here, m represents the coefficient number in the coefficient array, i represents the grid number, and k represents the parameter number.

8. The rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to claim 1, characterized in that: When reconstructing the temperature and velocity fields, POD-based P m Divided into P m,a With P m,b , and [P m,1 ,…P m,NX ,P m,NX+1 ,…,P m,2NX ] T =[P m,a,1 ,…P m,a,NX ,P m,b,1 ,…,P m,b,NX ] T .

9. The rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to claim 1, characterized in that: When solving the coefficient array, the diffusion term uses the Quick format instead of the center split format.

10. The rapid calculation method for nuclear reactor thermal loop system based on model order reduction according to claim 1, characterized in that: When verifying the reconstructed solution, any number within the range of u is chosen, u = x, x ∈ [min(u k ),max(u k The full-order data is calculated using the full-order model; the reconstructed solution is calculated using the POD basis, and the residuals of each variable at each node are calculated. And take the maximum value among them:

Images