A fully coupled system simulation method for reactors

Abstract

The present invention discloses a reactor fully coupled system simulation method, including: analyzing an object model, establishing a black box equation group; analyzing the black box equation group, generating a block matrix equation group; initializing various physical field parameters; preprocessing the block matrix equation group to accelerate JFNK solution; JFNK solution of the block matrix equation group; judging convergence, and adjusting the corresponding simulation step size; judging the simulation end time, and adjusting the corresponding simulation process. The present invention supports system simulation of various types of reactors, and the reactor fully coupled system simulation method supports multi-object, multi-scale, and multi-physical field coupling solution, which can solve the defects of the previous simulation software modeling model single, scale single, and less physical fields, supports black box solution, and does not need to understand the specific calculation process inside the black box model during use. The black box model only needs to provide the required input and output data, has strong expansion compatibility, and supports other software to participate in simulation as a black box model.

Description

Technical Field

[0001] The invention relates to the technical field of reactor simulation, and in particular to a system simulation method for a fully coupled reactor. Background Art

[0002] The reactor system is a complex multi-physics, multi-scale, and multi-object coupled system. The current mainstream system simulation programs at home and abroad, such as RELAP5, TRACE, etc., use rough pipeline nodes and thermal components to complete one-dimensional system simulation. This type of system simulation program has the defects of a single modeling model, a single scale, and fewer physical fields. With the development of computer technology and simulation technology, the simulation of various parts of the reactor system has become more refined and accurate. For example, CFD technology can be used to construct three-dimensional cores, steam generators, turbines and other structures, and accurately analyze the flow field, stress and other information of the equipment. The use of refined nuclear thermal coupling technology can simultaneously consider the interactive effects of the neutron field and the flow field, and calculate the power and flow field distribution with high fidelity. However, these refined and high-fidelity simulations are limited to a certain device or a small system, and do not support the simulation of the entire reactor system.

[0003] Due to the traditional modeling methods of pipeline nodes and thermal components, current system simulation is difficult to expand to high-dimensional and multi-physics field objects. Therefore, how to efficiently model and couple each object in its own appropriate way to form a multi-object, multi-scale, and multi-physics field system is a current research hotspot. Classic coupling methods include operator decomposition (OS) and Picard iteration. The use of traditional OS coupling methods may cause the physical field parameters within the time step to not converge, resulting in inaccurate or even wrong results calculated by the detailed model. The Picard method can make the physical field parameters converge within 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 that is compatible with multiple objects, multiple scales, and multiple physical fields is worth studying. Summary of the invention

[0004] The present invention provides a reactor fully coupled system simulation method to overcome at least one technical problem existing in the prior art.

[0005] An embodiment of the present invention provides a reactor fully coupled system simulation method, comprising:

[0006] Step S1, determining an object to be simulated in a reactor, and establishing several black box models of the object to be simulated based on the simulation purpose of the object to be simulated; establishing a function mapping relationship between input parameters and output parameters of each black box model in the several black box models; listing the obtained mapping relationship of each black box model in a matrix form to obtain a black box equation group corresponding to the several black box models;

[0007] Step S2, analyzing the black box equation group to generate an associated equation group corresponding to the black box equation group, wherein the associated equation group is used to describe the associated relationship between the input and output variables of each black box model, including boundary conditions, space constraints, and physical field equation constraints, so as to associate each black box model;

[0008] Combining the black box equation group and the association equation group to obtain a block matrix equation group, wherein the block matrix equation group is used to describe the association relationship between the black box models having the association relationship;

[0009] Step S3, assigning initial values ​​of various physical field parameters in the object to be simulated under steady-state conditions of the reactor system according to the physical field boundary conditions and initial conditions of each model in the object to be simulated;

[0010] Step S4, using the JFNK method to solve the block matrix equations, specifically includes:

[0011] Step S41, an outer loop calculation process based on the iterative solution of the JFNK method, including: constructing a Jacobi matrix-vector product linear equation group, updating the iterative solution and judging the iterative convergence; wherein, the construction of the Jacobi matrix-vector product linear equation group refers to deforming the nonlinear equation group in the form of an inexact Newton equation, and differentiating the Jacobi matrix-vector product to form a linear equation group for iterative solution in the Krylov subspace; the iterative solution update includes updating the inexact Newton iterative solution and the process parameters; the iterative convergence judgment refers to judging whether the equation residual meets the convergence requirement in the loop, if it meets the requirement, exiting the loop, if it does not meet the requirement, continuing the loop;

[0012] Step S42, using the Krylov subspace iteration method as an inner loop to solve the linear equation system;

[0013] Step S5, judging the convergence;

[0014] Step S6: Determine whether the simulation process is finished.

[0015] Preferably, if the object to be simulated is a core in the reactor, the simulation purpose of the object to be simulated is to establish several black box models of the object to be simulated, specifically including:

[0016] Taking the refinement requirements, computing resources, and error requirements when simulating the core as constraints, a three-dimensional neutron diffusion equation group and a three-dimensional two-phase thermal hydraulic equation group are used for discretization and modeling to calculate the fuel assembly power and coolant physical field respectively;

[0017] The three-dimensional neutron diffusion equations are as follows:

[0018]

[0019] The superscript n+1 represents the value of the variable at the next moment; the superscript n represents the value of the variable at this moment; the symbol Δt represents the time step, unit: s; the symbol Φ represents the neutron flux density, unit: cm -2 s -1 ; Symbol J represents neutron flux density, unit: cm -2 s -1 ; The symbol v represents the neutron movement speed, unit: cm / s; The symbol C represents the density of the delayed neutron precursor nucleus, unit: cm -3 ; Symbols a1, a2, a3, a4 are the coefficients of the correlation coefficient of the extended nodal method and the physical quantities including unit cell scale, core reaction cross section, core temperature, coolant temperature, and coolant density; symbol λ d It represents the decay constant of the delayed neutron precursor nucleus, unit: t -1 ; Symbol Σ f Indicates the fission cross section, unit: cm -1 ; The symbol β represents the delayed neutron fraction; the symbol k eff Indicates effective proliferation factor; symbol E f Indicates the energy produced by each fission, unit: J; symbol Q indicates the heat produced by unit fission, unit: W;

[0020] The three-dimensional two-phase thermal hydraulic equations are shown as follows:

[0021]

[0022] Among them, the symbol ρ g Indicates the density of the vapor phase, unit: kg / m 3 ; Symbol ρ l Indicates the density of the liquid phase, unit: kg / m 3 ; Symbol α g Indicates the cross-sectional vapor content, symbol α l Indicates the liquid content of the cross section; symbol h g Indicates the specific enthalpy of the vapor phase, unit: kJ / kg; symbol h l Indicates the specific enthalpy of the liquid phase, unit: kJ / kg; symbol u g Indicates the axial velocity of the vapor phase, unit: m / s; symbol u l Indicates the axial velocity of the liquid phase, unit: m / s; symbol w g Indicates the transverse velocity of the vapor phase, unit: m / s; symbol w lIt represents the lateral flow velocity of the liquid phase, unit: m / s; the symbol k represents the proportion of thermal power transferred from the fuel wall to the liquid phase; the symbol Γ represents the phase change mass flow rate of liquid phase vaporization to vapor phase, unit: kg / s; the symbol P represents the pressure, unit: Pa.

[0023] Preferably, the step S6 specifically includes: comparing the current simulation time with the required simulation time, if the required simulation time is not reached, calculating the simulation step length according to the Courant number, and continuing the calculation; if the required simulation time is reached, the program is terminated; wherein the Courant number calculation formula is as follows:

[0024]

[0025] where Δt is the maximum allowed time step; Δx i is the unit grid size, in cm; symbol U i is the maximum velocity in the unit grid, m / s; the symbol η represents an empirical coefficient, and in black box coupling, the empirical coefficient is greater than 0 and less than 1.

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

[0027] (1) The method of the present invention supports simulation of various types of reactor systems. The reactor full-coupled system simulation method of the present invention supports multi-object, multi-scale, and multi-physical field coupling solutions, which can solve the defects of previous simulation software such as single modeling model, single scale, and fewer physical fields.

[0028] (2) The present invention utilizes the JFNK strong coupling solver to ensure high convergence speed and high stability of coupling. Compared with traditional weak coupling methods such as the OS method or the Picard method, the JFNK method has higher convergence accuracy, faster convergence speed and better stability.

[0029] (3) The present invention supports black box solving. During use, it is not necessary to understand the specific calculation process inside the black box model. The black box model only needs to provide the required input and output data.

[0030] (4) The present invention has strong expansion compatibility and supports any software that can provide input and output interfaces to be embedded in the algorithm as a black box, which is conducive to users using mature models or software to participate in reactor system simulation calculations. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0032] Figure 1 A schematic diagram showing a flow chart of a reactor fully coupled system simulation method provided by the present invention;

[0033] Figure 2 The figure shows a JFNK solution flow chart in a reactor fully coupled system simulation method provided by the present invention. DETAILED DESCRIPTION

[0034] In order to make the purpose, technical solutions and advantages of one or more embodiments of this specification clearer, the technical solutions of one or more embodiments of this specification will be clearly and completely described below in combination with the specific embodiments of this specification and the corresponding drawings. Obviously, the described embodiments are only part of the embodiments of this specification, not all of the embodiments. Based on the embodiments in this specification, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of one or more embodiments of this specification.

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

[0036] As stated above, the current system simulation is difficult to expand to high-dimensional, multi-physics field objects due to the traditional pipeline node and thermal component modeling methods. Therefore, how to efficiently model and couple each object in its own appropriate way to form a multi-object, multi-scale, multi-physics field system is a current research hotspot. Classic coupling methods include operator decomposition method (OS) and Picard iteration method. The use of traditional OS coupling method may cause the physical field parameters within the time step to not converge, so that the fine model calculates inaccurate or even wrong results. The use of Picard method can make the physical field parameters converge within each step, but the amount of calculation 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 that is compatible with multiple objects, multiple scales, and multiple physical fields is worth studying. The present invention provides a system simulation method for fully coupled reactors to solve the technical problems in the prior art.

[0037] like Figure 1 As shown, the simulation method process may include the following steps:

[0038] Step S1: Analyze the object model and establish a black box equation group, which may specifically include:

[0039] Step S11: Analyze the object and establish a black box model of the simulation object. Taking a pressurized water reactor as an example, the simulation object can be equipment or systems such as the core, pressurizer, steam generator, main pump, turbine, condenser, deaerator, regenerator, valve, pipeline, etc., or it can be a component unit of the equipment or system. For different equipment, various adaptive methods are used to establish models in combination with the system simulation requirements. Taking the core as an example, after comprehensively considering the conditions such as refinement requirements, computing resources, and error requirements, the three-dimensional neutron diffusion equations and the three-dimensional two-phase thermal-hydraulic equations are used for discretization and modeling, and the fuel assembly power and coolant physical field are calculated respectively. Among them, the three-dimensional neutron diffusion equations are shown as follows:

[0040]

[0041] The superscript n+1 represents the value of the variable at the next moment; the superscript n represents the value of the variable at this moment; the symbol Δt represents the time step, unit: s; the symbol Φ represents the neutron flux density, unit: cm -2 s -1 ; Symbol J represents neutron flux density, unit: cm -2 s -1 ; The symbol v represents the neutron movement speed, unit: cm / s; The symbol C represents the density of the delayed neutron precursor nucleus, unit: cm -3 ; Symbols a1, a2, a3, a4 are the coefficients of the correlation coefficient of the extended nodal method and the physical quantities including unit cell scale, core reaction cross section, core temperature, coolant temperature, and coolant density; symbol λ d It represents the decay constant of the delayed neutron precursor nucleus, unit: t -1 ; Symbol Σ f Indicates the fission cross section, unit: cm -1 ; The symbol β represents the delayed neutron fraction; the symbol k eff Indicates effective proliferation factor; symbol E f It represents the energy produced by each fission, unit: J; the symbol Q represents the heat produced by unit fission, unit: W.

[0042] The three-dimensional two-phase thermal hydraulic equations are shown as follows:

[0043]

[0044] Among them, the symbol ρ g Indicates the density of the vapor phase, unit: kg / m 3 ; Symbol ρ l Indicates the density of the liquid phase, unit: kg / m 3 ; Symbol αg Indicates the cross-sectional vapor content, symbol α l Indicates the liquid content of the cross section; symbol h g Indicates the specific enthalpy of the vapor phase, unit: kJ / kg; symbol h l Indicates the specific enthalpy of the liquid phase, unit: kJ / kg; symbol u g Indicates the axial velocity of the vapor phase, unit: m / s; symbol u l Indicates the axial velocity of the liquid phase, unit: m / s; symbol w g Indicates the transverse velocity of the vapor phase, unit: m / s; symbol w l It represents the lateral flow velocity of the liquid phase, unit: m / s; the symbol k represents the proportion of thermal power transferred from the fuel wall to the liquid phase; the symbol Γ represents the phase change mass flow rate of liquid phase vaporization to vapor phase, unit: kg / s; the symbol P represents the pressure, unit: Pa.

[0045] Step S12: Describe the mapping relationship of the black box model.

[0046] As a system simulation method supporting black box models, the present invention does not need to know the specific modeling method or internal solution process of each object. However, the present invention requires each black box object to indicate that after inputting what kind of undetermined parameters or boundary conditions, it can return (output) the boundary conditions or physical quantities required by other black box models, as shown in the following formula:

[0047]

[0048] Among them, the symbol y i , (i=1,n) indicates the boundary conditions or physical quantities that need to be returned, the symbol n indicates the total number of returned variables; the symbol x i , (i=1,m) represents the unknown parameters or boundary conditions that need to be input, and the symbol m represents the total number of input variables.

[0049] Step S13: Generate a black box equation group.

[0050] The mapping relationship of each black box model is listed in matrix form to generate a black box equation group, as shown in the following formula:

[0051]

[0052]

[0053] ...

[0055] Among them, the subscripts A, B, and C refer to the names of the black box objects, and the black box equation group has no limit on the number of black box equations. Taking the core as an example, the black box equation group of the three-dimensional neutron diffusion equation group is in the following form:

[0056]

[0057] Among them, the input block vector x Nu and returns the block vector y Nu They are:

[0058] x Nu =[α g ,ρ l ,T F ,T l ] T

[0059] y Nu =[Q] T

[0060] That is, after inputting physical quantities such as the cross-sectional vapor content, coolant density, internal fuel temperature, and coolant temperature of each cell into the black box model of the three-dimensional neutron diffusion equation group, the black box model independently solves the three-dimensional neutron diffusion equation group and finally returns the heat generated by the fission of each cell.

[0061] Assuming that the boundary conditions and initial conditions of the flow field except the thermal power are fixed, the black box equations of the three-dimensional two-phase thermal hydraulic equations are as follows:

[0062]

[0063] Among them, the input block vector x Sub and returns the block vector y Sub They are:

[0064] x Sub =[Q] T

[0065] y Sub =[α g ,ρ l ,T F ,T l ] T

[0066] That is, after inputting the heat generated by fuel fission into the black box model of the three-dimensional two-phase thermal-hydraulic equations, the black box model independently solves the three-dimensional two-phase thermal-hydraulic equations, and finally returns the physical quantities such as the cross-sectional vapor content, coolant density, fuel internal temperature, and coolant temperature of each unit cell.

[0067] Step S2: Analyze the black box equations to generate a block matrix equation system, which may specifically include:

[0068] Step S21: Analyze the black box equations and generate the association equations. The association equations describe the relationship between the input and output variables of each black box object, such as boundary conditions, spatial constraints, physical field equation constraints, etc., and can associate each black box object. According to the input and output variables provided by each black box object, the association equations can be expressed as block matrices and block vectors as follows:

[0069] G AB (x A ,y A ,x B ,y B )=0

[0070]

[0071]

[0072]

[0073]

[0074] Among them, x A ,x B are the input block vectors of objects A and B respectively; A ,y B are the output block vectors of objects A and B respectively; the associated equation group G AB A block matrix consisting of the boundary conditions, spatial constraints, physical field equation constraints, etc. between the black box models. If there are physical quantities related between multiple black box objects, it is only necessary to expand the correlation function to the corresponding dimension. For example, the correlation equations of three black box objects are shown as follows:

[0075] G ABC (x A ,y A ,x B ,y B ,x C ,y C )=0

[0076]

[0077]

[0078]

[0079]

[0080]

[0081]

[0082] Taking the core as an example, the three-dimensional neutron diffusion equations are used to calculate the heat generation of each unit fuel, and the three-dimensional two-phase thermal hydraulic equations are used to calculate the coolant physical field. There is a correlation between the three-dimensional neutron diffusion equations and the three-dimensional two-phase thermal hydraulic equations, and the corresponding correlation equations are as follows:

[0083] x Nu -y Sub =0

[0084] x Sub -y Nu =0

[0085] Among them, the symbol x Nu 、Symbol x Sub are the input block vectors of the aforementioned three-dimensional neutron diffusion equations and three-dimensional two-phase thermal hydraulic equations; symbol y Nu 、Symbol y Sub They are the return block vectors of the aforementioned three-dimensional neutron diffusion equations and three-dimensional two-phase thermal-hydraulic equations, respectively.

[0086] Step S22: Generate a block matrix equation group.

[0087] The block matrix equation group is composed of a black box equation group and an associated equation group, and is used to describe the black box model and the relationship between the black box models. The present invention combines the black box equation group generated in step S1 and the associated equation group generated in step S21 to obtain a block matrix equation group, which is as follows:

[0088]

[0089]

[0090]

[0091] G AB (x A ,y A ,x B ,y B )=0

[0092] G BC (x B ,y B ,x C ,y C )=0

[0093] G AC (x A ,y A ,x C ,y C )=0

[0094] This set of equations should flexibly adjust the associated equations according to the specific coupling relationship between the black boxes and ensure that the entire block matrix equation has a specific solution.

[0095] Taking the core as an example, the black box equations of the three-dimensional neutron diffusion equations and the three-dimensional two-phase thermal hydraulic equations and the input and output variables x between the black box equations Nu 、x Sub ,y Nu ,y Sub The resulting set of associated equations forms a block matrix equation system as follows:

[0096]

[0097] Step S3: Initializing various physical field parameters includes:

[0098] According to the physical field boundary conditions and initial conditions of each model in the system, the initial values ​​of each physical field parameter are assigned when the reactor system is in steady state. Among them, the boundary conditions and initial conditions are generally physical quantities such as mass flow, pressure, temperature, specific enthalpy, neutron flux distribution, power, etc. In the present invention, because the system simulates each physical field to form a closed loop and satisfies the conservation equation, the initial values ​​of each physical field do not have to be strictly consistent with the physical field parameters of the system in steady state, and the assigned initial value can be roughly equal to the actual steady-state parameter.

[0099] Optionally, step S4 may also be included, namely: preprocessing of the block matrix equation group to accelerate the JFNK solution, including:

[0100] Step S41: differential approximation of partial derivatives of black box model equations. Since the system simulation coupling object is a black box model, the partial derivative expression of the black box equation cannot be accurately known. Therefore, the present invention uses a differential method to approximate the partial derivatives of the black box equation, and the first-order difference is shown as follows:

[0101]

[0102] Where F is a black box function, x i is an independent variable of the black box function.

[0103] Step S42: Block matrix preprocessing, the differential approximation of each black box equation and the partial derivatives of the associated equation are used to form a Jacobi matrix for preprocessing. For example, the Jacobi matrix of the three black box equations is as follows:

[0104]

[0105] The Jacobi matrix is ​​used for the preprocessing calculation of the Krylov subspace in JFNK. When using it, the Jacobi matrix is ​​inverted after incomplete LU decomposition and then applied to the equation matrix to complete the preprocessing process. This can reduce the number of iterations and memory overhead of the Krylov subspace to speed up the JFNK solution.

[0106] Step S5: JFNK solves the block matrix equations including:

[0107] Step S51: Inexact Newton iteration. The inexact Newton iteration method is the outer loop of the JFNK iterative solution. The inexact Newton iteration includes: constructing the Jacobi matrix-vector product linear equations, updating the iterative solution and judging the iterative convergence. Among them, the iterative solution update includes the updating of the inexact Newton iterative solution and the process parameters. Iterative convergence judgment refers to judging whether the equation residual meets the convergence requirements in the loop. If it meets the requirements, the loop is exited, and if it does not meet the requirements, the loop continues. Constructing the Jacobi matrix-vector product linear equations refers to deforming the nonlinear equations in the form of inexact Newton equations, and differentiating the Jacobi matrix-vector products to form a linear equations for Krylov subspace iterative solution. The inexact Newton equations are shown below:

[0108]

[0109] Among them, F(x k ) is the kth iteration function value, F′(x k ) is the Jacobi matrix value of the kth iteration, is the iteration step length to be solved, η k is a mandatory term. However, since the Jacobi matrix is ​​difficult to solve and the Jacobi matrix in the inexact Newton equation appears in the form of matrix-vector product with the iteration step, the JFNK method uses differential approximation as follows:

[0110]

[0111] Among them, ε is a small constant, which converts the complex inexact Newton equations into a system of linear equations, making it easier to use the Krylov subspace method to efficiently solve the system of equations in the next step.

[0112] Step S52: Krylov subspace iteration solves the linear equations. In the JFNK solution process, the Krylov subspace iteration method is the inner iteration loop algorithm of the JFNK iterative solution, which is used to solve large linear equations. If there is no preprocessing matrix, the Krylov subspace iteration can be used directly to solve the equations. If there is a preprocessing matrix, such as the left preprocessing matrix is ​​denoted as M, the inexact Newton equation after left preprocessing is:

[0113]

[0114] In the above formula, in one iteration step, except for the iteration step length Except for the unknown number, all other numbers are known. The Krylov subspace iteration method can be used to complete the calculation of the iterative step with a higher convergence speed.

[0115] Step S6: Determine the convergence and make corresponding adjustments, including:

[0116] Judgment of the convergence of the JFNK solution: During the JFNK solution process, if the preset number of iterations (for example, in a possible embodiment, the number of iterations can be set to 10,000 times) does not converge, it is considered that the equation solution does not converge. When the equation calculation does not converge, the current simulation time step can be divided by 2, and step S4 is re-executed to generate a block preprocessing matrix. If the equation calculation converges, it is considered that the physical field parameter calculation at this simulation moment has been completed, and the next step is performed.

[0117] Step S7: Determine the end of the simulation and make corresponding adjustments, which may specifically include:

[0118] Determination of 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 length according to the Courant number and execute step S4 again. If the required simulation time is reached, the program terminates. The Courant number calculation formula is as follows:

[0119]

[0120] The symbol Δt is the maximum allowed time step; the symbol Δx i is the unit grid size, unit: cm; symbol U i It represents the maximum velocity in the unit grid, unit: m / s; the symbol η represents the empirical coefficient. In black box coupling, the value of this empirical coefficient is generally greater than 0 and less than 1.

[0121] Those skilled in the art can understand that the accompanying drawings are only schematic diagrams of an embodiment, and the modules or processes in the accompanying drawings are not necessarily required to implement the present invention.

[0122] 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 of the embodiment, or can be changed accordingly and located in one or more devices different from the embodiment. The modules in the above embodiment can be combined into one module, or can be further divided into multiple sub-modules.

[0123] 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 it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A reactor fully coupled system simulation method, characterized in that: The simulation method comprises: Step S1, determining an object to be simulated in a reactor, and establishing several black box models of the object to be simulated based on the simulation purpose of the object to be simulated; establishing a function mapping relationship between input parameters and output parameters of each black box model in the several black box models; listing the obtained mapping relationship of each black box model in a matrix form to obtain a black box equation group corresponding to the several black box models; Step S2, analyzing the black box equation group to generate an associated equation group corresponding to the black box equation group, wherein the associated equation group is used to describe the associated relationship between the input and output variables of each black box model, including boundary conditions, space constraints, and physical field equation constraints, so as to associate each black box model; Combining the black box equation group and the association equation group to obtain a block matrix equation group, wherein the block matrix equation group is used to describe the association relationship between the black box models having the association relationship; Step S3, assigning initial values ​​of various physical field parameters in the object to be simulated under steady-state conditions of the reactor system according to the physical field boundary conditions and initial conditions of each model in the object to be simulated; Step S4, using the JFNK method to solve the block matrix equations, specifically includes: Step S41, an outer loop calculation process based on the iterative solution of the JFNK method, including: constructing a Jacobi matrix-vector product linear equation group, updating the iterative solution and judging the iterative convergence; wherein, the construction of the Jacobi matrix-vector product linear equation group refers to deforming the nonlinear equation group in the form of an inexact Newton equation, and differentiating the Jacobi matrix-vector product to form a linear equation group for iterative solution in the Krylov subspace; the iterative solution update includes updating the inexact Newton iterative solution and the process parameters; the iterative convergence judgment refers to judging whether the equation residual meets the convergence requirement in the loop, if it meets the requirement, exiting the loop, if it does not meet the requirement, continuing the loop; Step S42, using the Krylov subspace iteration method as an inner loop to solve the linear equation system; Step S5, judging the convergence; Step S6: Determine whether the simulation process is finished.

2. The reactor fully coupled system simulation method according to claim 1, characterized in that: If the object to be simulated is the core of the reactor, the simulation purpose of the object to be simulated is to establish several black box models of the object to be simulated, specifically including: Taking the refinement requirements, computing resources, and error requirements when simulating the core as constraints, a three-dimensional neutron diffusion equation group and a three-dimensional two-phase thermal hydraulic equation group are used for discretization and modeling to calculate the fuel assembly power and coolant physical field respectively; The three-dimensional neutron diffusion equations are as follows: The superscript n+1 represents the value of the variable at the next moment; the superscript n represents the value of the variable at this moment; the symbol Δt represents the time step, unit: s; the symbol Φ represents the neutron flux density, unit: cm -2 s -1 ; Symbol J represents neutron flux density, unit: cm -2 s -1 ; The symbol v represents the neutron movement speed, unit: cm / s; The symbol C represents the density of the delayed neutron precursor nucleus, unit: cm -3 ; Symbols a1, a2, a3, a4 are the coefficients of the correlation coefficient of the extended nodal method and the physical quantities including unit cell scale, core reaction cross section, core temperature, coolant temperature, and coolant density; symbol λ d It represents the decay constant of the delayed neutron precursor nucleus, unit: t -1 ; Symbol Σ f Indicates the fission cross section, unit: cm -1 ; The symbol β represents the delayed neutron fraction; the symbol k eff Indicates effective proliferation factor; symbol E f Indicates the energy produced by each fission, unit: J; symbol Q indicates the heat produced by unit fission, unit: W; The three-dimensional two-phase thermal hydraulic equations are shown as follows: Among them, the symbol ρ g Indicates the density of the vapor phase, unit: kg / m 3 ; Symbol ρ l Indicates the density of the liquid phase, unit: kg / m 3 ; Symbol α g Indicates the cross-sectional vapor content, symbol α l Indicates the liquid content of the cross section; symbol h g Indicates the specific enthalpy of the vapor phase, unit: kJ / kg; symbol h l Indicates the specific enthalpy of the liquid phase, unit: kJ / kg; symbol u g Indicates the axial velocity of the vapor phase, unit: m / s; symbol u l Indicates the axial velocity of the liquid phase, unit: m / s; symbol w g Indicates the transverse velocity of the vapor phase, unit: m / s; symbol w l It represents the lateral flow velocity of the liquid phase, unit: m / s; the symbol k represents the proportion of thermal power transferred from the fuel wall to the liquid phase; the symbol Γ represents the phase change mass flow rate of liquid phase vaporization to vapor phase, unit: kg / s; the symbol P represents the pressure, unit: Pa.

3. The reactor fully coupled system simulation method according to claim 1, characterized in that: The step S6 specifically includes: comparing the current simulation time with the required simulation time, if the required simulation time is not reached, calculating the simulation step length according to the Courant number, and continuing the calculation; if the required simulation time is reached, the program is terminated; wherein the Courant number calculation formula is as follows: where Δt is the maximum allowed time step; Δx i is the unit grid size, in cm; symbol U i is the maximum velocity in the unit grid, m / s; the symbol η represents an empirical coefficient, and in black box coupling, the empirical coefficient is greater than 0 and less than 1.

Images