CN116579205A - A Calculation Method of PWR Nuclear-Heat Coupling - Google Patents

Abstract

The invention discloses a calculation method for nuclear thermal coupling of a pressurized water reactor. The scheme may include: establishing a neutron diffusion equation group and a thermal hydraulic equation group of the pressurized water reactor core to be analyzed, and combining the neutron diffusion equation group and The thermal hydraulic equations are coupled to obtain the nuclear thermal coupling equations; the format of the preprocessing matrix and the residual equations is dynamically adjusted; the preprocessing matrix and the residual equations are dynamically updated; the JFNK combined with the preprocessing matrix is used to solve the nuclear thermal coupling Residual equations; judging the convergence of JFNK solution and adjusting the simulation step size; judging whether the simulation process is over.

Description

A Calculation Method of PWR Nuclear-Heat Coupling

technical field

The invention relates to the technical field of nuclear reactor cores, in particular to a calculation method for nuclear-thermal coupling of a pressurized water reactor.

Background technique

Reactor-nuclear-thermal coupling technology can combine the complex neutron diffusion equations and thermal-hydraulic equations in the reactor to complete the simulation of core neutrons, thermal-hydraulic and other physical fields. Currently, the traditional coupling technology is mainly through the Operator Splitting (OS) method or the Picard iteration method. In the OS method, the partial solutions of the neutron diffusion equation module and the thermal-hydraulic equation module are used as new estimates of the boundary conditions of the other side, and are alternately calculated in each time step until the end of the simulation. Due to the lag in updating the physical field parameters in the transient calculation, this method requires a small coupling time step, otherwise it will cause a large error. The Picard iterative method requires neutron diffusion equations and thermal-hydraulic equations to iteratively converge in each time step before advancing to the next time step. The Picard iterative method is a method that supports full implicit coupling, which can ensure that the convergence is completed within each time step, and the time step size can be increased appropriately, but there are problems of slow calculation speed and weak calculation stability.

The Jacobian-Free Newton–Krylov (JFNK) method is a strong coupling algorithm, which puts the diffusion equations and thermal hydraulic equations to be coupled in a large matrix for unified solution, so as to Guarantee the unified update of parameters during the iterative solution process, thereby improving the convergence speed and reducing the iterative error. Generally speaking, compared with the traditional operator decomposition method or Picard iteration method, the JFNK method has faster convergence speed, higher convergence accuracy, and better calculation stability. Therefore, the use of JFNK method is beneficial to improve the calculation efficiency of nuclear thermal coupling simulation of PWR. However, in the transient simulation, when the PWR coolant phase changes, if the two-phase fluid equation is used to be compatible with the single-phase fluid equation, it may appear ill-conditioned, or even degenerate into an indeterminate equation system, resulting in distortion. However, if the phase transition occurs, the corresponding equations are listed according to the specific situation, and then the simultaneous solution can solve the problems of ill-conditioned equations and indeterminate equations. However, when the JFNK method itself solves large-scale matrices simultaneously, its calculation process is relatively complicated, and its expansion ability is not strong. When the number and format of the physical field equations change, the entire matrix to be solved will change, which is a challenge for the JFNK algorithm. Therefore, how to use the JFNK method to solve the two-phase nuclear thermal coupling of the PWR stably and efficiently is a technical problem to be solved.

Contents of the invention

The invention provides a calculation method for nuclear heat coupling of a pressurized water reactor, which is used to overcome at least one technical problem existing in the prior art.

An embodiment of the present invention provides a calculation method for nuclear thermal coupling of a pressurized water reactor, including:

Step S1, establishing the neutron diffusion equation group and the thermal hydraulic equation group of the pressurized water reactor core to be analyzed, coupling the neutron diffusion equation group and the thermal hydraulic equation group to obtain the nuclear thermal coupling equation group ; Using the extended node method to discretize the neutron diffusion equations to obtain the first grid set sum, and using the sub-channel method to discretize the thermal and hydraulic equations to obtain the second grid set and; establish the mapping relationship between the first grid set and the second grid set; the first grid set and the second grid set are based on the PWR The core is obtained by using the same grid division method;

According to the boundary conditions and initial conditions of the physical field of the pressurized water reactor core, initialize the parameters of the physical field of the core, and utilize the few-group cross-section library initialization to include transport cross-section, scattering cross-section, absorption cross-section, fission cross-section, middle The physical quantity including the reaction cross section such as sub-generation cross section;

The initial time step is calculated according to the following formula, which is calculated as follows:

Among them, the symbol Δt represents the maximum allowable time step; the symbol _Δxi represents the unit grid size, the unit is m; the symbol U _i represents the maximum velocity in the unit grid, the unit is m/s; the symbol η represents the empirical coefficient greater than 0 ;

Step S2, dynamically adjust the format of the preprocessing matrix and residual equations, specifically including:

Step S21, read the physical parameters at the previous moment relative to the moment to be analyzed, calculate the process parameters, and store the process parameters in the computer memory;

Step S22, judging the phase state of each grid in the first grid set and the second grid set, for the first arbitrary grid, if the first arbitrary grid is liquid phase, then the two-phase residual equations corresponding to the grid degenerate into liquid phase residual equations, and the vapor phase physical quantity in the corresponding vector to be solved is deleted; if the first arbitrary grid is a vapor phase, then the grid The two-phase residual equations corresponding to the grid degenerate into vapor-phase residual equations, and the liquid-phase physical quantities in the corresponding unsolved vectors are deleted; if the first arbitrary grid is two-phase, the two-phase residual equations and The vector to be solved remains unchanged; the phase state parameters of all grids in the first grid set and the second grid set are traversed, the number of equations to be solved and physical quantities to be solved is counted, and the corresponding preprocessing matrix is determined and the scale and generation form of residual equations;

Step S3, dynamically updating the preprocessing matrix and residual equations, specifically including:

Step S31, reading the physical parameters and process parameters at the previous moment from the memory;

Step S32, judging the phase state of each grid in the first grid set and the second grid set, for the second arbitrary grid, if the second arbitrary grid is two phase, then fill in the grid’s residual equations and vectors to be solved based on the two-phase physical parameter database; if the second arbitrary grid is liquid phase, fill in the grid’s residual equation based on the liquid phase physical parameter Group and vectors to be solved, if the second arbitrary grid is a vapor phase, fill in the residual equations and the vectors to be solved of the grid based on the vapor phase physical parameter database, until the residual equations and the vectors to be solved are completed;

Step S4, combining the JFNK of the preprocessing matrix to solve the nuclear thermal coupling residual equation group;

Step S5, judging the convergence of the JFNK solution and adjusting the simulation step size;

Step S6, judging whether the simulation process ends.

Preferably, if the second arbitrary grid is two-phase, fill in the grid's residual equations and vectors to be solved based on the two-phase physical parameter database, specifically including:

Based on physical parameters including neutron flux density and two-phase pressure, temperature, flow velocity, and cross-sectional gas fraction read out from the memory at the previous moment, calculate neutron flux density , macroscopic reaction cross-section, total heat power density, and process parameters including two-phase heat transfer coefficient, wall heat power, and wall resistance coefficient.

Preferably, if the second arbitrary grid is a liquid phase, fill in the grid's residual equations and vectors to be solved based on the liquid phase physical parameter database, specifically including:

Based on the physical parameters including neutron flux density, liquid phase pressure, temperature, and flow velocity read from the memory at the last moment, calculate neutron flux density, macroscopic reaction cross section, Process parameters including the total heat power density, heat transfer coefficient of the liquid phase, wall heat power, and wall resistance coefficient.

Preferably, if the second arbitrary grid is a vapor phase, fill in the grid's residual equations and vectors to be solved based on the vapor phase physical parameter database, specifically including:

Based on the physical parameters including neutron flux density, vapor phase pressure, temperature, and flow rate read from the memory about the previous moment, calculate neutron flux density, macroscopic reaction cross section, Process parameters including the total heat power density, heat transfer coefficient of the vapor phase, wall heat power, and wall resistance coefficient.

Preferably, when the phase state is two phases, the nuclear thermal coupling equations are:

x＝[Φ,C,Q,P,α _g ,u _g ,u _l ,w _g ,w _l ,h _g ,h _l ] ^T

Among them, f _Φ represents the neutron flux residual equation group; f _C represents the delayed neutron concentration residual equation group; f _Q represents the fission energy residual equation group; f _P represents the vapor phase mass residual equation group; Represents the liquid phase mass residual equation group; /> Represents the vapor phase axial momentum residual equations; /> Represents the liquid phase axial momentum residual equations; /> Represents the residual equations of the vapor phase transverse momentum; Represents the liquid phase lateral momentum residual equations; /> Represents the vapor phase energy residual equations; /> Indicates the liquid phase energy residual equations; the superscript n+1 is the variable value of the current time step, and the superscript n is the variable value of the previous time step; Δt is the time step, and the unit is s; Φ is the neutron Flux density, the unit is cm ^-2 s ^-1 ; v is the neutron velocity, the unit is cm/s; C is the density of the delayed neutron precursor nucleus, the unit is cm ^-3 ; J is the neutron flux density, The unit is cm ^-2 s ^-1 ; a ₁ , a ₂ , a ₃ , and a ₄ are the correlation coefficients of the extended nodal method; λ _d is the decay constant of the delayed neutron precursor nucleus, and the unit is t ^-1 ; Σ _f is The fission cross section is in cm ^-1 ; β and k _eff are the proportion of delayed neutrons and the effective multiplication factor respectively; E _f is the energy produced by a single fission in J; Q is the energy produced by fission in W; V is The volume of the grid, the unit is m ³ ; A _u and A _w are the flow areas perpendicular to the axial and transverse velocities, the unit is m ² ; ρ _g , ρ _l are the densities of the vapor phase and the liquid phase, respectively , the unit is kg/m ³ ; α _g , α _l are the vapor content and liquid content of the section respectively; h _g , h _l are the specific enthalpy of the vapor phase and the liquid phase respectively, and the unit is kJ/kg; u _g , u _l is the axial flow velocity of vapor phase and liquid phase respectively, the unit is m/s; w _g , w _l are the transverse flow velocity of vapor phase and liquid phase respectively, m/s; Γ is the phase of vaporization from liquid phase Variable mass flow rate, the unit is kg/s; P _J , P _J+1 are the axial upstream and downstream grid pressures, the unit is Pa; P _ii , P _jj are the horizontal upstream and downstream grid pressures, the unit is is Pa; τ _wg,u , τ _wl,u are the wall resistances of the axial vapor phase and liquid phase, respectively, in kg m/s ² ; τ _wg,w , τ _wl,w are the transverse vapor phase and liquid The wall resistance of phase, the unit is kg _· m _/ s ² ; τ _gl is the resistance between vapor and liquid, the unit is kg·m/s ² ; Total heat, in W; n _j represents the total number of adjacent grids in the axial direction; n _gap represents the total number of adjacent grids in the lateral direction.

One embodiment of this specification can at least achieve the following beneficial effects:

(1) The PWR nuclear-thermal coupling calculation method of the present invention can efficiently and stably calculate the two-phase nuclear-thermal coupling problem of the PWR, and avoid problems such as ill-conditioned and uncertain solutions.

(2) The present invention utilizes the JFNK method to carry out two-phase nuclear thermal coupling simultaneous solution to the PWR core, and each physical field variable is updated simultaneously in the JFNK iteration, compared with traditional OS methods, methods such as Picard, has higher Convergence accuracy, faster convergence speed and better stability.

(3) The dynamic Jacobi preprocessing method and the dynamic residual equation group method of the present invention dynamically adjust the scale and format of the preprocessing matrix and the residual equation group according to the phase state in each discrete grid, avoid infinitesimal quantities, and reduce pathological conditions , to eliminate indefinite equations and avoid indefinite solutions.

(4) The present invention proposes an intelligent calculation time step scheme, fully utilizes the fully implicit iterative feature supported by the JFNK algorithm, increases the simulation time step, reduces the number of simulations, and greatly improves calculation efficiency.

(5) The present invention does not need to limit the discrete form of the neutron diffusion equation and thermal-hydraulic equations, but only needs to use discrete equations to describe the physical variable relationship between cells to construct the corresponding Jacobi preconditioning matrix and difference equations.

Description of drawings

In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings that need to be used in the description of the embodiments or the prior art. Obviously, the accompanying drawings in the following description are only These are some embodiments of the present invention. Those skilled in the art can also obtain other drawings based on these drawings without creative work.

Fig. 1 represents the method flow diagram of the computing method of a kind of pressurized water reactor nuclear thermal coupling provided by the present invention;

Fig. 2 represents the flow chart of dynamically adjusting the preprocessing matrix and the format of equations to be solved in a kind of calculation method of PWR nuclear thermal coupling provided by the present invention;

Fig. 3 shows a flow chart of dynamically updating the preprocessing matrix and the equations to be solved in a calculation method for nuclear thermal coupling of a PWR provided by the present invention.

Detailed ways

In order to make the purpose, technical solutions and advantages of one or more embodiments of this specification more clear, the following will clearly and completely describe the technical solutions of one or more embodiments of this specification in conjunction with specific embodiments of this specification and corresponding drawings . Apparently, the described embodiments are only some of the embodiments in this specification, not all of them. Based on the embodiments in this specification, all other embodiments obtained by persons of ordinary skill in the art without making creative efforts belong to the protection scope of one or more embodiments in this specification.

It should be understood that although the terms first, second, third, etc. may be used in this application document to describe various information, such information should not be limited to these terms. These terms are only used to distinguish information of the same type from one another.

As stated above, current system simulations are difficult to extend to high-dimensional, multi-physics objects due to the traditional modeling methods of piping nodes and thermal components. Therefore, how to efficiently model and couple each object in an appropriate way to form a multi-object, multi-scale, multi-physics system is a current research hotspot. The classic coupling methods include operator decomposition method (OS) and Picard iteration method. Using the traditional OS coupling method may lead to non-convergence of the physical field parameters within the time step, making the fine model calculate inaccurate or even wrong results. However, the Picard method can make the physical field parameters converge in each step, but the calculation amount is large, the convergence stability is poor, and the calculation efficiency is low. Therefore, a system simulation method with high convergence speed and high stability compatible with multiple objects, multiple scales, and multiple physical fields is worth studying. The invention provides a calculation method for nuclear heat coupling of a pressurized water reactor, which is used to solve the technical problems in the prior art.

As shown in Figure 1, the flow of the simulation method may include the following steps:

Step S1: Initialize equations, time step and physical field parameters:

Step S11: Initialize the equation group. The invention can combine neutron diffusion equations and thermal hydraulic equations into a large nuclear-heat coupled equations. Taking the typical core modeling of PWR as an example, the neutron diffusion equations can be discretized using the extended block method to obtain the first grid set sum, and the thermal and hydraulic equations can be processed using the sub-channel method The discretization process obtains the sum of the second grid set, and then establishes a mapping relationship between the first grid set and the second grid set. When performing specific grid division, the first grid set sum and the second grid set sum can be obtained based on the same grid division method for the pressurized water reactor core, so as to facilitate the establishment of the discrete network of the extended nodal method The mapping relationship between the grid and the discrete grid of sub-channels. Wherein, when the phase state is two phases, the nuclear thermal coupling equation is as follows:

x＝[Φ,C,Q,P,α _g ,u _g ,u _l ,w _g ,w _l ,h _g ,h _l ] ^T

Among them, f _Φ represents the neutron flux residual equation group; f _C represents the delayed neutron concentration residual equation group; f _Q represents the fission energy residual equation group; f _P represents the vapor phase mass residual equation group; Represents the liquid phase mass residual equation group; /> Represents the vapor phase axial momentum residual equations; /> Represents the liquid phase axial momentum residual equations; /> Represents the residual equations of the vapor phase transverse momentum; Represents the liquid phase lateral momentum residual equations; /> Represents the vapor phase energy residual equations; /> Indicates the liquid phase energy residual equations; the superscript n+1 is the variable value of the current time step, and the superscript n is the variable value of the previous time step; Δt is the time step, and the unit is s; Φ is the neutron Flux density, the unit is cm ^-2 s ^-1 ; v is the neutron velocity, the unit is cm/s; C is the density of the delayed neutron precursor nucleus, the unit is cm ^-3 ; J is the neutron flux density, The unit is cm ^-2 s ^-1 ; a ₁ , a ₂ , a ₃ , and a ₄ are the correlation coefficients of the extended nodal method; λ _d is the decay constant of the delayed neutron precursor nucleus, and the unit is t ^-1 ; Σ _f is The fission cross section is in cm ^-1 ; β and k _eff are the proportion of delayed neutrons and the effective multiplication factor respectively; E _f is the energy produced by a single fission in J; Q is the energy produced by fission in W; V is The volume of the grid, the unit is m ³ ; A _u and A _w are the flow areas perpendicular to the axial and transverse velocities, the unit is m ² ; ρ _g , ρ _l are the densities of the vapor phase and the liquid phase, respectively , the unit is kg/m ³ ; α _g , α _l are the vapor content and liquid content of the section respectively; h _g , h _l are the specific enthalpy of the vapor phase and the liquid phase respectively, and the unit is kJ/kg; u _g , u _l is the axial flow velocity of vapor phase and liquid phase respectively, the unit is m/s; w _g , w _l are the transverse flow velocity of vapor phase and liquid phase respectively, m/s; Γ is the phase of vaporization from liquid phase Variable mass flow rate, the unit is kg/s; P _J , P _J+1 are the axial upstream and downstream grid pressures, the unit is Pa; P _ii , P _jj are the horizontal upstream and downstream grid pressures, the unit is is Pa; τ _wg,u , τ _wl,u are the wall resistances of the axial vapor phase and liquid phase, respectively, in kg m/s ² ; τ _wg,w , τ _wl,w represent the lateral vapor phase and liquid phase resistance The wall resistance of phase, the unit is kg _· m _/ s ² ; τ _gl is the resistance between vapor and liquid, the unit is kg·m/s ² ; Total heat, in W; n _j represents the total number of adjacent grids in the axial direction; n _gap represents the total number of adjacent grids in the lateral direction.

Step S12: Initializing physical field parameters. According to the boundary conditions and initial conditions of the core physical field, the parameters of the core physical field are initialized, and then physical quantities such as the reaction cross section are initialized by using the few-group macroscopic cross-section library. The boundary conditions and initial conditions of the physical field of the core are generally physical quantities such as flow velocity, pressure, temperature, and specific enthalpy. In this example, the flow velocity and specific enthalpy at the core inlet and the pressure at the core outlet are selected. The small group macroscopic reaction cross-section generally includes transport cross-section, scattering cross-section, absorption cross-section, fission cross-section, neutron generation cross-section and so on. The small-group macroscopic cross-section library can query the corresponding cross-section data according to the information of fluid density, temperature, boron concentration and fuel temperature.

Step S13: Initialize the time step. Calculate the initial time step according to the Courant number, and store it in memory for subsequent calculation calls. The formula for calculating the time step is as follows:

Among them, Δt is the maximum allowed time step; Δx _i is the unit grid size, m; U _i is the maximum velocity in the unit grid, m/s; η is the empirical coefficient, greater than 0, when the equations are fully implicit format, it can take a value of 3 or even higher.

Step S2: Dynamically adjust the format of the preprocessing matrix and residual equations, as shown in Figure 2, which may specifically include:

Step S21: Read the solution vector parameters at the previous moment, and calculate the process parameters, and store them in memory for other steps to read. Because nuclear thermal coupling is a multi-physics simulation, considering the number of grids and the number of equations to be solved, if all the physical quantities are used as the number of solutions to be solved and written as large vectors to be solved, the residual equation and preprocessing matrix will be very large. Huge and difficult to calculate. However, in the present invention, the heat conduction equations of solid structural parts such as fuel rods, control rod conduits, and positioning grids are used as process parameter equations, which is beneficial to ensure calculation accuracy while reducing calculation amount and improving calculation efficiency. Because the process parameter equations will be updated and solved after each inexact Newton iteration of JFNK nuclear thermal coupling, and participate in the next inexact Newton iteration until convergence, the vector to be solved and the process parameters have been updated, satisfying the full implicit coupling requirements, the calculation accuracy is small. Solving the process parameter equations independently will help reduce the scale of the nuclear-thermal coupling equations, greatly reduce the amount of JFNK nuclear-thermal coupling calculations, and improve computational efficiency.

Step S22: Determine the composition format of the preprocessing matrix and residual equations. Due to the change of the phase state in each cell in the transient calculation, the number and format of the residual equations will change, and the corresponding vectors to be solved will change. As a result, the preprocessing matrix also changes. Therefore, it is necessary to traverse the phase state parameters of all grids, count the number of equations and physical quantities to be solved, and determine the size and generation form of the corresponding preprocessing matrix and residual equations.

In this step, the phase state of each grid in the first grid set and the second grid set can be judged. For the first arbitrary grid, if the first arbitrary grid is a liquid phase, then The two-phase residual equations corresponding to this grid degenerate into liquid-phase residual equations, and the vapor phase physical quantity in the corresponding vector to be solved is deleted; The phase residual equations degenerate into vapor phase residual equations, and the liquid phase physical quantities in the corresponding vectors to be solved are deleted; if the first arbitrary grid is two-phase, the two-phase residual equations and the vectors to be solved remain unchanged ; Traverse the first grid set and the phase state parameters of all grids in the second grid set and count the number of equations to be solved and physical quantities to be solved, and determine the corresponding preprocessing matrix and residual equations size and form.

More specifically, the present invention can judge which phase state the grid is in according to the numerical results of the cross-sectional vapor content of each grid in the previous step, and dynamically adjust the corresponding equations. If it is two-phase, the two-phase residual equations and the vector to be solved remain unchanged. If it is a liquid phase, the two-phase residual equations corresponding to the grid degenerate into a liquid phase residual equations, corresponding to the deletion of the vapor phase physical quantities in the vector to be solved, as shown in the following formula:

Similarly, if it is a vapor phase, the two-phase residual equations corresponding to the grid degenerate into vapor phase residual equations, corresponding to the deletion of liquid phase physical quantities in the vector to be solved, as shown in the following formula:

The present invention adopts the Jacobi matrix as the JFNK preprocessing matrix, and the preprocessing matrix is generated by partial derivatives of the discrete equations of all grids and the variables to be solved, and its block matrix form is as follows:

When a grid phase state is single-phase, its residual equation is only a single-phase equation. When forming the Jacobi matrix, use the single-phase residual equation to obtain partial derivatives of its independent variables to obtain partial derivatives, and place the partial derivatives in corresponding positions of the Jacobi matrix.

It can be seen that when there is a single-phase grid in the present invention, since the two-phase equations degenerate into single-phase equations, the overall residual equations and preprocessing matrix scale will be reduced, thus reducing the amount of calculation. can increase the calculation speed.

Step S3: Dynamically update the preprocessing matrix and residual equations, as shown in Figure 3, which may specifically include:

Step S31: Read the solution vector parameters, process parameters and time step at the last moment from the memory.

Step S32: Update the parameters of the residual equations according to the format of the residual equations generated in the second step. Using the solution vector parameters, process parameters and time steps read from the memory at the previous moment, the residual equations are updated.

Specifically, as shown in Figure 3, if any second grid in the sum of the first grid set and the sum of the second grid is two-phase, fill in the residual of the grid based on the two-phase physical parameter database Equations and vectors to be solved, which can specifically include:

Based on physical parameters including neutron flux density and two-phase pressure, temperature, flow velocity, and cross-sectional gas fraction read out from the memory at the previous moment, calculate neutron flux density , macroscopic reaction cross-section, total heat power density, and process parameters including two-phase heat transfer coefficient, wall heat power, and wall resistance coefficient.

If the second arbitrary grid is a liquid phase, fill in the grid's residual equations and vectors to be solved based on the liquid phase physical parameter database, which may specifically include:

Based on the physical parameters including neutron flux density, liquid phase pressure, temperature, and flow velocity read from the memory at the last moment, calculate neutron flux density, macroscopic reaction cross section, Process parameters including the total heat power density, heat transfer coefficient of the liquid phase, wall heat power, and wall resistance coefficient.

If the second arbitrary grid is a vapor phase, fill in the grid's residual equations and vectors to be solved based on the vapor phase physical parameter database, which may specifically include:

Based on the physical parameters including neutron flux density, vapor phase pressure, temperature, and flow rate read from the memory about the previous moment, calculate neutron flux density, macroscopic reaction cross section, Process parameters including the total heat power density, heat transfer coefficient of the vapor phase, wall heat power, and wall resistance coefficient.

Step S33: Update the preprocessing matrix according to the format of the preprocessing matrix generated in the second step. The preprocessing matrix is updated by using the solution vector parameters, process parameters and time steps read from the memory at the previous moment. This example uses the Jacobi matrix as the preprocessing matrix. The Jacobi matrix is a large matrix formed by taking partial derivatives of all the elements in the solution vector for each residual equation system. Its format is as follows:

The Jacobi matrix is used in the preprocessing calculation of the Krylov subspace in JFNK to improve the calculation efficiency of the Krylov subspace method. If the left preprocessing is used, the Jacobi matrix is incompletely decomposed by LU and then inverted, and then left multiplied on the matrix to be solved.

Step S4: Combining with the preprocessing matrix, using JFNK to solve the residual equations of nuclear thermal coupling includes:

Step S41: Non-exact Newton iteration. The non-exact Newton iteration method is the outer loop of the JFNK iterative solution. The non-exact Newton iteration includes: constructing the Jacobi matrix-vector product linear equation system, updating the iterative solution, and judging the convergence of the iterative process. Among them, the update of iterative decoupling includes the update of the Krylov subspace solution and the update of the process parameters. Iterative convergence judgment refers to judging whether the residual of the equation meets the convergence requirements in the outer loop, and exits the loop if the requirements are met, and continues the outer loop if the requirements are not met. Constructing the Jacobi matrix-vector product linear equation system refers to transforming the nonlinear equation system in the form of an inexact Newton equation, and taking the difference of the Jacobi matrix-vector product to form a linear equation system for iterative solution of the Krylov subspace. The inexact Newton equation looks like this:

Among them, F(x _k ) is the function value of the kth iteration, F′(x _k ) is the Jacobi matrix value of the kth iteration, is the iteration step size to be solved, and η _k is a mandatory item. However, since it is difficult to solve the Jacobi matrix, and the Jacobi matrix in the inexact Newton equation appears in the form of matrix-vector product with the iterative step size, the JFNK method adopts the differential approximation as follows:

Among them, ε is a small constant. In this way, the complex inexact Newton equations are transformed into linear equations, which is convenient for the next step to use the Krylov subspace method to efficiently solve the equations.

Step S42: Krylov subspace iteratively solves the system of linear equations. In the process of solving JFNK, the Krylov subspace iteration method is an inner loop algorithm of JFNK iterative solution, which is used to solve large linear equations. The present invention uses the preprocessing matrix in the third step and combines the Krylov subspace iteration method to solve the iteration step size stably and efficiently. Record the left preprocessing matrix as M, and the non-exact Newton equation after left preprocessing is:

The above formula is in an iteration step, except for the iteration step size is unknown, all other numbers are known. Using the Krylov subspace iteration method can complete the calculation of the iteration step with a high convergence speed.

Step S5: Judging the convergence of the JFNK solution and adjusting the simulation step size includes:

Judgment of JFNK solution convergence. During the JFNK solving process, if the preset number of iterations (for example, in a technical solution of a possible embodiment, the number of iterations can be set to 1,000 times) does not converge, it is determined that the solution of the equation does not converge. When the equation calculation does not converge, divide the current simulation time step by 2, update the time step value in the memory, and re-execute step 3 (dynamically update the preprocessing matrix and residual equations) to update the preprocessing matrix and residual equations. If the calculation of the equation converges, it is considered that the calculation of the physical field parameters at this simulation moment is completed, and the next step is performed.

Step S6: Judging the end time of the simulation, and making corresponding adjustments to the simulation process, which may specifically include:

Judgment of the simulation end time: compare the current simulation time with the required simulation time, if the required simulation time is not reached, calculate the simulation step size according to the Courant number, and execute the second step (determine the formation). The program terminates when the required simulation time is reached.

In the technical solution of this application, when using the JFNK method to solve the aforementioned large-scale nuclear-thermal coupled equations process, physical quantities such as neutron field, temperature field, pressure field, and flow velocity field will be updated simultaneously in each non-exact Newton iteration step, compared to The OS method and Picard method, which solve the physical fields separately, have higher calculation accuracy of nuclear thermal coupling. Moreover, compared with the first-order convergence speed of the OS method and the Picard method, the JFNK method adopted in the present invention has a second-order convergence speed and faster calculation speed. In addition, the present invention dynamically adjusts the corresponding equations according to the phase state of the discrete grid. Compared with using two-phase equations to model the phase states of each grid without distinction, the present invention has fewer equations and a smaller amount of calculation, and It can avoid ill-conditioned problems when calculating single-phase physical fields with two-phase equations.

Those skilled in the art can understand that the accompanying drawing is only a schematic diagram of an embodiment, and the modules or processes in the accompanying drawing are not necessarily necessary for implementing the present invention.

Those skilled in the art can understand that: the modules in the device in the embodiment can be distributed in the device in the embodiment according to the description in the embodiment, and can also be changed and located in one or more devices different from the embodiment. The modules in the above embodiments can be combined into one module, and can also be further split into multiple sub-modules.

Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that: it can still be Modifications are made to the technical solutions described in the foregoing embodiments, or equivalent replacements are made to some of the technical features; these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims (5)

1. The calculation method of the pressurized water reactor nuclear thermal coupling is characterized by comprising the following steps:

step S1, establishing a neutron diffusion equation set and a thermal hydraulic equation set of a pressurized water reactor core to be analyzed, and coupling the neutron diffusion equation set and the thermal hydraulic equation set to obtain a nuclear thermal coupling equation set; discretizing the neutron diffusion equation set by using an expansion joint block method to obtain a first grid set sum, and discretizing the thermal hydraulic equation set by using a sub-channel method to obtain a second grid set sum; establishing a mapping relation between the first grid set and the second grid set; the first grid set and the second grid set are obtained based on the same grid division mode of the pressurized water reactor core;

initializing parameters of a physical field of the reactor core according to boundary conditions and initial conditions of the physical field of the pressurized water reactor core, and initializing physical quantities including a transport section, a scattering section, an absorption section, a fission section, a neutron production section and other reaction sections by using a few-group section library;

the initial time step is calculated according to the following formula:

wherein the symbol Δt represents the maximum time step allowed; symbol Deltax _i Representing the cell grid size in m; symbol U _i Representing the maximum speed in the cell grid in m/s; the symbol η represents an empirical coefficient greater than 0;

step S2, dynamically adjusting the format of the preprocessing matrix and the residual equation set, and specifically comprising the following steps:

s21, reading each physical parameter at the previous moment relative to the moment to be analyzed, calculating a process parameter, and storing the process parameter into a computer memory;

step S22, judging the phase states of each grid in the first grid set and the second grid set, and for any first grid in the first grid set, if any first grid is a liquid phase, degrading a two-phase residual equation set corresponding to the grid into a liquid-phase residual equation set, and deleting vapor phase physical quantities in corresponding vectors to be solved; if any one of the first grids is a vapor phase, degrading a two-phase residual equation set corresponding to the grid into a vapor phase residual equation set, and deleting liquid phase physical quantity in a corresponding vector to be solved; if the first arbitrary grid is two phases, the two-phase residual equation set and the vector to be solved are unchanged; traversing phase parameters of all grids in the first grid set and the second grid set, counting the number of equations to be solved and physical quantities to be solved, and determining the scale and the generation form of a corresponding preprocessing matrix and residual equation set;

step S3, dynamically updating a preprocessing matrix and a residual equation set, wherein the method specifically comprises the following steps:

step S31, reading all physical parameters and process parameters of the previous moment from a memory;

step S32, judging the phase states of each grid in the first grid set and the second grid set, for a second arbitrary grid in the first grid set and the second grid set, if the second arbitrary grid is two-phase, filling a residual equation set and a vector to be solved of the grid based on a two-phase physical parameter database, if the second arbitrary grid is liquid phase, filling the residual equation set and the vector to be solved of the grid based on the liquid phase physical parameter database, if the second arbitrary grid is vapor phase, filling the residual equation set and the vector to be solved of the grid based on the vapor phase physical parameter database until the residual equation set and the vector to be solved are completed;

s4, solving a nuclear thermal coupling residual equation set by combining JFNK of a pretreatment matrix;

s5, judging the convergence of JFNK solution and adjusting the simulation step length;

and S6, judging whether the simulation process is finished.

2. The method for calculating nuclear thermal coupling of pressurized water reactor according to claim 1, wherein if the second arbitrary grid is two-phase, filling out a residual equation set and a vector to be solved of the grid based on a two-phase physical parameter database, specifically comprising:

based on the physical parameters including neutron flux density, pressure, temperature, flow rate and section gas content of the two phases at the previous moment read from the memory, process parameters including neutron flux density, macroscopic reaction section, total thermal power density, heat exchange coefficient of the two phases, wall thermal power and wall resistance coefficient are calculated.

3. The method for calculating nuclear thermal coupling of pressurized water reactor according to claim 1, wherein if the second arbitrary grid is a liquid phase, filling out a residual equation set and a vector to be solved of the grid based on a liquid phase physical parameter database, specifically comprising:

based on the physical parameters including neutron flux density, liquid phase pressure, temperature and flow velocity read from the memory at the previous moment, process parameters including neutron flux density, macroscopic reaction section, total heat power density, liquid phase heat exchange coefficient, wall heat power and wall resistance coefficient are calculated.

4. The method for calculating nuclear thermal coupling of pressurized water reactor according to claim 1, wherein if said second arbitrary grid is a vapor phase, filling out a residual equation set and a vector to be solved of the grid based on a vapor phase physical parameter database, specifically comprising:

based on the physical parameters including neutron flux density, vapor phase pressure, temperature and flow rate read from the memory at the previous time, process parameters including neutron flux density, macroscopic reaction section, total thermal power density, vapor phase heat exchange coefficient, wall thermal power and wall resistance coefficient are calculated.

5. The method of claim 1, wherein when the phase is two phases, the set of nuclear thermal coupling equations is:

x＝[Φ,C,Q,P,α _g ,u _g ,u _l ,w _g ,w _l ,h _g ,h _l ] ^T

wherein f _Φ Representing a neutron flux residual equation set; f (f) _C Representing a delayed neutron concentration residual equation set; f (f) _Q Representing a system of residual equations of the fracture energy; f (f) _P Representing a set of vapor phase mass residual equations;representing a liquid phase mass residual equation set; />Representing a set of vapor phase axial momentum residual equations; />Representing a liquid phase axial momentum residual equation set; />Representing a set of vapor phase lateral momentum residual equations; />Representing a liquid phase lateral momentum residual equation set; />Representing a set of vapor phase energy residual equations; />Representing a liquid phase energy residual equation set; the superscript n+1 is the variable value of the current time step, and the superscript n is the variable value of the previous time step; Δt is the time step in s; phi is neutron flux density in cm ^-2 s ^-1 The method comprises the steps of carrying out a first treatment on the surface of the v is neutron velocity in cm/s; c is the density of the delayed neutron precursor nuclei in cm ^-3 The method comprises the steps of carrying out a first treatment on the surface of the J is neutron flux density in cm ^-2 s ^-1 ；a ₁ 、a ₂ 、a ₃ 、a ₄ Is the correlation coefficient of the expansion block method; lambda (lambda) _d Is the decay constant of the delayed neutron precursor nucleus, and has the unit of t ^-1 ；Σ _f Is a fissile section in cm ^-1 ；β、k _eff The neutron fraction and the effective multiplication factor are delayed respectively; e (E) _f Is energy generated by single fission, and the unit is J; q is fission generating energy, and the unit is W; v is the volume of the mesh in m ³ ；A _u 、A _w Flow areas orthogonal to the axial and transverse velocities, respectively, in m ² ；ρ _g 、ρ _l The density of the vapor phase and the liquid phase, respectively, in kg/m ³ ；α _g 、α _l The cross section steam content and the cross section liquid content are respectively; h is a _g 、h _l Specific enthalpy of vapor phase and liquid phase, respectively, in kJ/kg; u (u) _g 、u _l The axial flow rates of the vapor phase and the liquid phase are respectively expressed in m/s; w (w) _g 、w _l The lateral flow rates, m/s, of the vapor and liquid phases, respectively; Γ is the phase change mass flow in kg/s of the liquid phase vaporized into the vapor phase; p (P) _J 、P _J+1 The pressures of the axially upstream and downstream grids, respectively, in Pa; p (P) _ii 、P _jj The pressures of the transverse upstream and downstream grids are Pa; τ _wg,u 、τ _wl,u Wall resistance in kg.m/s for the axial vapor and liquid phases, respectively ² ；τ _wg,w 、τ _wl,w Wall drag in kg.m/s for the transverse vapor and liquid phases, respectively ² ；τ _gl Is the resistance of the vapor phase and the liquid phase, and the unit is kg.m/s ² ；q _g 、q _l The total heat transferred to the vapor phase and the liquid phase in the grid is W; n is n _j Representing the total number of adjacent grids in the axial direction; n is n _gap Representing the total number of adjacent grids in the transverse direction.