Calculation method for neutron transport equation of nuclear reactor

Abstract

The invention discloses a calculation method for a nuclear reactor neutron transport equation, and belongs to the field of nuclear reactor core design and safety. According to the method, the problem that the convergence speed of an inner iteration process of a transport equation and a low-order equation in a traditional PI-CMFD method is low is solved. The method specifically comprises the following steps: step 1, reading a coarse mesh model, a fine mesh model, material information and boundary conditions of a nuclear reactor; step 2, establishing a numerical discrete model of a steady-state multi-group transport equation according to the fine grid model read in the step 1; 3, after boundary constraint and material constraint are applied to the established numerical discrete model of the steady-state multi-group transport equation, a discrete equation set is constructed by transforming a test function in the model; and 4, solving the constructed discrete equation set based on the grid mapping relation between the coarse grid model and the fine grid model by adopting a PI method and a CMFD method to obtain a neutron angular flux value and a characteristic value of the neutron transport equation. The method can be applied to calculation of the neutron transport equation of the nuclear reactor.

Description

Technical Field

[0001] The invention belongs to the field of nuclear reactor core design and safety, and in particular relates to a calculation method for a neutron transport equation of a nuclear reactor. Background Art

[0002] The neutron transport equation describes the specific behavior of neutrons in a nuclear reactor system and is the basis for the design and research of a nuclear reactor system. There are currently two mainstream methods for solving the neutron transport equation: the power iteration method (PI method) and the Newton-Kronov method (JFNK method) that does not require explicit construction of the Jacobian matrix. Although the JFNK method has a faster convergence speed in most cases, it requires a better initial value, otherwise it may converge to an erroneous result. In general, an iterative acceleration method is introduced to improve the convergence speed of the PI method. The coarse mesh finite difference acceleration method (CMFD) is an efficient acceleration technique. Some studies have shown that a reasonably optimized PI method accelerated by CMFD (PI-CMFD) even achieves a faster convergence speed than the JFNK method.

[0003] The PI method is divided into external iteration and internal iteration during the calculation process. The external iteration is mainly used to update the eigenvalues ​​and fission sources, while the internal iteration is mainly responsible for updating the scattering source size and flux size of neutrons in different energy groups. CMFD will also construct a low-order neutron conservation eigenvalue equation. The iterative solution process of this equation is similar to that of the transport equation, and the PI method can also be used. The internal iteration process of the PI method generally uses the source iteration method, but this process converges slowly in strong scattering problems, and the CMFD method mainly accelerates the external iteration process.

[0004] In summary, in order to further improve the convergence speed of transport calculations, it is necessary to accelerate the internal iteration process of the transport equations and low-order equations in the PI-CMFD method. Summary of the invention

[0005] The purpose of the present invention is to solve the problem of slow convergence speed of the inner iteration process of the transport equation and the low-order equation in the traditional PI-CMFD method, and to propose a calculation method for the neutron transport equation of the nuclear reactor.

[0006] The technical solution adopted by the present invention to solve the above technical problems is: a calculation method for the neutron transport equation of a nuclear reactor, the method specifically comprising the following steps:

[0007] Step 1: read the coarse grid model, fine grid model, material information and boundary conditions of the nuclear reactor, and construct a grid mapping relationship between the coarse grid model and the fine grid model;

[0008] Step 2: Establish a numerical discrete model of the steady-state multi-group transport equation based on the fine grid model read in step 1;

[0009] Step 3, after applying boundary constraints and material constraints to the established numerical discrete model of steady-state multi-group transport equations according to material information and boundary conditions, a discrete equation group is constructed by transforming the test function in the numerical discrete model;

[0010] Step 4: Use the PI method and the CMFD method, and solve the constructed discrete equations based on the grid mapping relationship between the coarse grid model and the fine grid model to obtain the neutron angular flux value and eigenvalue of the neutron transport equation.

[0011] Furthermore, the numerical discrete model of the steady-state multi-group transport equation is established based on the discrete ordinate discontinuous finite element method.

[0012] Furthermore, the numerical discrete model of the steady-state multi-group transport equation is specifically:

[0013]

[0014] in, represents the migration behavior of neutrons, represents the neutron consumption behavior, represents the scattering behavior of neutrons, Represents the fission behavior of neutrons.

[0015] Furthermore, the migration behavior of the neutrons Neutron consumption behavior Neutron scattering behavior and the fission behavior of neutrons Specifically:

[0016]

[0017] in, represents the migration of neutrons on the boundary surface, represents the migration of neutrons between adjacent units, represents the migration of neutrons within the cell;

[0018]

[0019]

[0020]

[0021] in, is the spatial position coordinate; is the jth test function; ▽ is the gradient operator; D is the set of all units of the fine grid model, is the collection of fine mesh boundary cell faces of the nuclear reactor; in the collection yes The outward normal vector of the fine mesh cell surface at ; on the set Γ of the fine mesh internal cell surfaces of the nuclear reactor, yes Satisfaction The element face normal vector of ; and w m They are the discrete quadrature groups in the discrete ordinate method. The discrete angles and the angle integral weights corresponding to the discrete angles, M is the number of discrete angles, Indicates the energy group g and the discrete direction The magnitude of the neutron angular flux under ; Represented on the internal element faces of the fine mesh The interruption, is the neutron angular flux of the upstream unit in the direction of neutron motion,

[0022]

[0023]

[0024]

[0025]

[0026]

[0027]

[0028] in, Is in the energy group g, position The total cross-section of the material at; L is the order of Legendre expansion, l = 0, 1, ..., L; Is in the energy group g, position The neutron generation cross section of the material at is in position The probability that the energy group of the new neutron produced when the neutron undergoes fission reaction with the material is exactly g; It's location The l-th-order Legendre moment of the scattering cross section from energy group g' to g at , g' = 1, 2, ..., G, G is the total number of energy groups; is the discrete angle The l-order p-th spherical harmonic function, p = -l, -l+1, ..., l; is in the energy group g′ and the discrete direction The magnitude of the neutron angular flux underm is the discrete angle The corresponding angle integral weight; It's location The magnitude of the neutron flux at energy group g'; is the i′th basis function in element e in the discontinuous finite element method, ψ g,m,i′ yes is the expansion coefficient of , and N is the number of basis functions in unit e.

[0029] Furthermore, the specific process of applying boundary constraints to the numerical discrete model of the steady-state multi-group transport equation is:

[0030] make is the reflection boundary, is a vacuum boundary, then

[0031]

[0032]

[0033]

[0034] For vacuum boundaries The incident neutron flux is 0, that is, when hour Then we have:

[0035]

[0036] For reflective boundaries When two discrete directions and satisfy hour, and exist There is a reflection relationship about the reflection boundary;

[0037] like but yes The reflection angle is like but yes The reflection angle is Then we have:

[0038]

[0039] Furthermore, the specific process of applying material constraints to the numerical discrete model of the steady-state multi-group transport equation is:

[0040] Let D = D1∪D2, and after applying material constraints on material D1, become become become After applying material constraints on material D2, become become become but

[0041]

[0042]

[0043]

[0044]

[0045]

[0046]

[0047] Let D = D aniso ∪D iso , D aniso represents the material region where anisotropic scattering is considered and the Legendre expansion order is L, D iso represents the material area where anisotropic scattering is not considered, then:

[0048]

[0049]

[0050]

[0051] In the formula, is the case when only isotropic scattering is considered. Abbreviation of .

[0052] Furthermore, the matrix operator form of the discrete equation system is:

[0053]

[0054] Where L includes the behaviors represented by equations (3) to (7); F represents the fission behavior of neutrons; ψ is ψ in equation (11) g,m,i′ The vector composed.

[0055] Furthermore, the specific process of step 4 is as follows:

[0056] Step 401, initialize the outer iteration number n=0, and initialize the eigenvalue and solution vector:

[0057] k0 =k0 (27)

[0058]

[0059] Among them, k 0 is the initialized eigenvalue, ψ 0 is the initialized solution vector;

[0060] Step 402: Update the fission source Q of the neutron transport equation using the eigenvalue and solution vector n :

[0061]

[0062] Among them, |·| represents the modulus of the vector;

[0063] Step 403: According to the fission source Q n The inner iteration equation is established, and the Krylov subspace method is used to solve equation (30):

[0064]

[0065] Among them, ψ n+1 is the updated solution vector, and then according to ψ n+1 Update feature values:

[0066] k n+1 =|Fψ n+1 | (31)

[0067] Among them, k n+1 is the updated feature value;

[0068] Let the number of outer iterations of the PI method for solving the neutron transport equation be n = n + 1;

[0069] Step 404: According to ψ n Determine whether the convergence condition is met;

[0070] If the convergence condition is met, the calculation process ends and outputs ψ n and k n ;

[0071] Otherwise, continue to execute step 405;

[0072] Step 405: Using ψ n and k n To update the parameters of the coarse grid model, that is, to calculate equations (32) to (39), and then to construct a low-order diffusion approximation equation based on the updated coarse grid model parameters;

[0073] The specific process of step 405 is as follows:

[0074]

[0075]

[0076]

[0077]

[0078]

[0079]

[0080]

[0081]

[0082] Among them, V c (I) is the volume of coarse mesh unit I; is the average neutron flux size of fine mesh unit i in energy group g, i represents the fine mesh unit covered by coarse mesh unit I; yes Initial value; is the neutron flux of coarse grid unit I in energy group g; V i is the volume of fine mesh unit i; is the total cross-section of the material at unit I in the coarse mesh under energy group g; Δx(I) is the width of coarse mesh unit I in the direction of the line connecting the centers of the two; Δx(J) is the width of coarse mesh unit J in the direction of the line connecting the centers of the two; yes Initial value; is the neutron flux size of the coarse mesh unit J under energy group g; is the neutron flux diffusion coefficient obtained based on finite differences; is the neutron flux correction factor obtained using the transport neutron flux, It's location The 0th-order Legendre moment of the scattering cross section from energy group g' to g at ;

[0083] The intermediate variable Defined as:

[0084]

[0085]

[0086]

[0087] in, is the outward normal vector of the coarse mesh unit I on the common surface, is the external normal vector of the coarse mesh unit J on the common surface; Sc (I, J) is the area of ​​the common surface of two adjacent coarse mesh units I and J;

[0088] The low-order diffusion approximation equation is constructed based on the updated coarse-grid model parameters and the finite difference method:

[0089]

[0090]

[0091]

[0092] in, is the net neutron flow from coarse mesh unit I to coarse mesh unit J; is the scattering cross section of the coarse-grid element I from energy group g' to g; is the neutron flux of coarse grid unit I in energy group g'; k c is the eigenvalue of the low-order diffusion approximation equation of the coarse grid model; is the neutron generation cross section of the material at coarse mesh unit I; is the probability that the energy group of the new neutron produced when the neutron undergoes fission reaction with the material at the coarse mesh unit I is exactly g;

[0093] Substituting equations (44) and (45) into equation (43), we obtain the finite difference discrete model of the coarse grid model:

[0094]

[0095]

[0096]

[0097]

[0098] in, is the vector consisting of the neutron fluxes of the coarse-grid units of energy group g, The first element in is the neutron flux size of the coarse mesh unit I, is the vector consisting of the sub-scalar fluxes of the coarse-grid units of energy group g′, yes The element in row I and column J in yes The element in row I and column J in yes The element in the I-th row and J-th column, J adjacent to I means that the coarse mesh unit I and the coarse mesh unit J are adjacent;

[0099] The finite difference discrete model is written in the form of equation (50):

[0100]

[0101]

[0102]

[0103]

[0104]

[0105] Step 406, initialize the inner layer iteration number l′=0, and initialize the neutron flux vector and eigenvalue of the low-order diffusion approximation equation:

[0106] k c,0 =k n (55)

[0107]

[0108]

[0109] Step 407: Use the coarse-grid neutron flux vector Φ c,l' And the fission source term on the right side of the eigenvalue update equation (50):

[0110]

[0111] Step 408: c and k c To update:

[0112] L c Φ c,l′+1 =Q c,l′ (59)

[0113] k c,l′+1 =|F c Φ c,l′+1 | (60)

[0114] Among them, Φ c,l′+1 is the updated coarse-grid neutron flux vector, k c,l′+1 is the updated feature value;

[0115] Step 409, set the number of PI method outer iterations of the coarse grid model to l′=l′+1;

[0116] Determine whether the convergence condition is met, if the convergence condition is met, execute step 410, if not, return to execute step 407;

[0117] Step 410: Calculate ψ n+1 / 2 and kn+1 / 2 :

[0118]

[0119] k n+1 / 2 =k c,final (62)

[0120] in, and k c,final They are the neutron flux size and eigenvalue size after the PI method of the coarse grid model converges;

[0121]

[0122] Step 411: Using ψ n+1 / 2 and k n+1 / 2 Return to step 402.

[0123] Furthermore, the convergence condition in step 404 is:

[0124]

[0125]

[0126] Among them, ε ψ is the convergence limit of the solution vector, ε k is the convergence limit of the eigenvalue of the transport equation.

[0127] The beneficial effects of the present invention are:

[0128] The method of the present invention uses the PI method accelerated by the CMFD method to solve the transport equation and has the following advantages:

[0129] The discontinuous finite element method (DFEM) is used to discretize the spatial variables in the transport equation, which has good compatibility with unstructured grids.

[0130] The transport equation and the low-order equations of the CMFD acceleration method are solved using the PI method, and the inner iteration equations of both are solved using the Krylov subspace method. Compared with the traditional source iteration method, the Krylov method can further improve the convergence speed and solution accuracy of the inner iteration.

[0131] When solving linear problems of the form Ax=b, the Krylov solver is fully utilized. It only requires the product of matrix A and vector v (v is a vector of the same length as x), without caring about the characteristics of the calculation process. The complex matrix-vector product in the transport calculation is replaced by an equivalent integral calculation.

[0132] Considering that the main composition of the matrix A of different transport problems in actual simulation is basically the same, its composition is divided into items according to their respective integral calculation formulas and implemented through modular programming. Then, in actual calculations, it is called by region and boundary as needed. This not only improves the functional extensibility of the program, but also increases the code reuse rate. BRIEF DESCRIPTION OF THE DRAWINGS

[0133] Figure 1 It is a flow chart of a calculation method for a nuclear reactor neutron transport equation of the present invention;

[0134] Figure 2 It is a schematic diagram of a fine mesh;

[0135] Figure 3 It is a rough network diagram. DETAILED DESCRIPTION

[0136] Specific implementation method 1: Combination Figure 1 The present embodiment is described as follows. The present embodiment describes a method for calculating the neutron transport equation for a nuclear reactor, and the method specifically comprises the following steps:

[0137] Step 1: Figure 2 and Figure 3 As shown, the coarse grid model, fine grid model, material information and boundary conditions of the nuclear reactor are read, and a grid mapping relationship between the coarse grid model and the fine grid model is established;

[0138] Step 2: Establish a numerical discrete model of the steady-state multi-group transport equation based on the fine grid model read in step 1;

[0139] Step 3, after applying boundary constraints and material constraints to the established numerical discrete model of steady-state multi-group transport equations according to material information and boundary conditions, a discrete equation group is constructed by transforming the test function in the numerical discrete model;

[0140] Step 4: Use the PI method and the CMFD method, and solve the constructed discrete equations based on the grid mapping relationship between the coarse grid model and the fine grid model (the grid mapping relationship means which fine grids are covered by the coarse grid for a coarse grid) to obtain the neutron angular flux value and eigenvalue of the neutron transport equation.

[0141] Specific implementation method 2: This implementation method is different from the specific implementation method 1 in that the numerical discrete model of the steady-state multi-group transport equation is established based on the discrete ordinate discontinuous finite element method (SN-DFEM method).

[0142] The other steps and parameters are the same as those in the first embodiment.

[0143] The steady-state SN multigroup transport equation is:

[0144]

[0145] Multiply both ends of the above equation by the test function Then, by performing volume integration on the entire fine mesh, we can obtain the following formula:

[0146]

[0147] in,

[0148]

[0149]

[0150]

[0151]

[0152] because so

[0153]

[0154] And using Gaussian divergence theorem:

[0155]

[0156] Where, e is the fine mesh unit, is the element face set of the fine mesh element e, is the outward normal vector of the fine mesh element e on the element surface. So,

[0157]

[0158] Because the internal cell surface of each fine mesh is shared by two adjacent fine mesh cells, the two adjacent fine mesh cells each have The direction of the neutron movement is exactly the opposite. The adjacent fine mesh unit of this fine mesh unit located downstream of the neutron movement direction So we have:

[0159]

[0160]

[0161]

[0162]

[0163] in, is the collection of fine mesh boundary cell faces of the nuclear reactor; in the collection superior, yes The outward normal vector of the fine mesh cell surface at ; on the set Γ of the fine mesh internal cell surfaces of the nuclear reactor, yes Satisfaction The element face normal vector.

[0164] Apply the upwind format (i.e., use replace then:

[0165]

[0166] so:

[0167]

[0168] For the discontinuous finite element method, the neutron angular flux Can be used in any unit e To approximate, ψ g,m,i yes The expansion coefficient, N is the number of basis functions in unit e;

[0169]

[0170] Based on the above theory, a numerical discrete model of the steady-state multi-group transport equation of the specific implementation method three is established.

[0171] Specific implementation method three: This implementation method is different from specific implementation method one or two in that the numerical discrete model of the steady-state multi-group transport equation is specifically:

[0172]

[0173] in, represents the migration behavior of neutrons, represents the neutron consumption behavior, represents the scattering behavior of neutrons, Represents the fission behavior of neutrons.

[0174] The other steps and parameters are the same as those in the first or second embodiment.

[0175] Specific implementation method 4: This implementation method is different from the specific implementation methods 1 to 3 in that the migration behavior of the neutrons Neutron consumption behavior Neutron scattering behavior and the fission behavior of neutrons Specifically:

[0176]

[0177] in, represents the migration of neutrons on the boundary surface, represents the migration of neutrons between adjacent units, represents the migration of neutrons within the cell;

[0178]

[0179]

[0180]

[0181] in, is the spatial position coordinate; is the jth test function; ▽ is the gradient operator; D is the set of all units of the fine grid model, is the collection of fine mesh boundary cell faces of the nuclear reactor; in the collection superior, yes The outward normal vector of the fine mesh cell surface at ; on the set Γ of the fine mesh internal cell surfaces of the nuclear reactor, yes Satisfaction The element face normal vector of ; and w m They are the discrete quadrature groups in the discrete ordinate method (SN method) The discrete angles and the angle integral weights corresponding to the discrete angles, M is the number of discrete angles, Indicates the energy group g and the discrete direction The magnitude of the neutron angular flux under ; Represented on the internal element faces of the fine mesh The interruption, is the neutron angular flux of the upstream unit in the direction of neutron motion,

[0182]

[0183]

[0184]

[0185]

[0186]

[0187]

[0188] in, Is in the energy group g, position The total cross-section of the material at; L is the order of Legendre expansion, l = 0, 1, ..., L; Is in the energy group g, position The neutron generation cross section of the material at is in position The probability that the energy group of the new neutron produced when the neutron undergoes fission reaction with the material is exactly g; It's location The l-th-order Legendre moment of the scattering cross section from energy group g' to g at , g' = 1, 2, ..., G, G is the total number of energy groups; is the discrete angle The l-order p-th spherical harmonic function, p = -l, -l+1, ..., l; is in the energy group g′ and the discrete direction The magnitude of the neutron angular flux under m is the discrete angle The corresponding angle integral weight; It's location The magnitude of the neutron flux at energy group g'; is the i′th basis function in element e in the discontinuous finite element method, the neutron angular flux Can be used in any unit e To approximate, ψ g,m,i′ yes is the expansion coefficient of , and N is the number of basis functions in unit e.

[0189] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0190] Specific implementation method 5: This implementation method is different from any one of the specific implementation methods 1 to 4 in that the specific process of applying boundary constraints to the numerical discrete model of the steady-state multi-group transport equation is as follows:

[0191] make is the reflection boundary, is a vacuum boundary, then

[0192]

[0193]

[0194]

[0195] For vacuum boundaries The incident neutron flux is 0, that is, when hour Then we have:

[0196]

[0197] For reflective boundaries When two discrete directions and satisfy hour, and exist There is a reflection relationship about the reflection boundary;

[0198] like but yes The reflection angle is like but yes The reflection angle is like Then no one can find a solution that satisfies the reflection relation. At this time, the boundary surface integral must also be 0, so:

[0199]

[0200] The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

[0201] Different boundary conditions are expressed mathematically as follows: The constraints imposed are different. These boundary conditions will be imposed on the discrete model by substituting these constraints into equation (3). and vacuum boundary Take the nuclear reactor as an example.

[0202] In actual calculations, there may be other types of boundary conditions or more than two boundaries at the same time. Different boundary conditions will only cause the program to impose different boundary constraints on different boundary surfaces, which will lead to different substitutions. After the constraint is set, the integral calculation formula corresponding to formula (3) changes, such as formula (15) and formula (16) in this implementation mode; and the more the number of boundary conditions, The more subsets it is divided into (similar to and ), each subset represents a boundary condition, It can also be expanded into more integral terms over different subsets of boundaries (similar to and ), but this does not affect the calculation process of the algorithm and will not be listed in detail here.

[0203] Specific implementation method 6: This implementation method is different from the specific implementation methods 1 to 5 in that the specific process of applying material constraints to the numerical discrete model of the steady-state multi-group transport equation is as follows:

[0204] In different material regions, the influence of material properties on the problem will be reflected in equations (6) to (8) and constraints.

[0205] For equation (6) and equation (8), and The constraint is relatively simple. Let D = D1 ∪ D2, and after applying the material constraint on material D1, become become become After applying material constraints on material D2, become become become but

[0206]

[0207]

[0208]

[0209]

[0210]

[0211]

[0212] For equation (7), since the Legendre expansion order of the scattering source considered on different materials may be different, its constraint application is more complicated. Here, let D = D aniso ∪D iso , D aniso represents the material region where anisotropic scattering is considered and the Legendre expansion order is L, D iso represents the material region where anisotropic scattering is not considered (i.e. the Legendre expansion order is 0), then:

[0213]

[0214]

[0215]

[0216] In the formula, is the case when only isotropic scattering is considered. Abbreviation of .

[0217] The other steps and parameters are the same as those in Specific Implementations 1 to 5.

[0218] Similar to the boundary conditions, if there are more materials in D, D can be further subdivided into different subsets according to the different material constraints, and then similar to equations (17), (20), and (23) in this embodiment, and Expanded into the integral sum over multiple subsets. If there are multiple materials with and without anisotropic scattering, D aniso and D iso Then expand it like D, and the integral processing method after expansion is similar to equation (17) and equation (20) in this embodiment.

[0219] Specific implementation method 7: This implementation method is different from any one of specific implementation methods 1 to 6 in that the matrix operator form of the discrete equation group is:

[0220]

[0221] Wherein, L includes the behaviors represented by equations (3) to (7), indicating the reaction processes such as neutron migration, consumption, and scattering; F represents the fission behavior of neutrons (i.e., equation (8)); ψ is ψ in equation (11) g,m,i′ The vector composed.

[0222] The other steps and parameters are the same as those in Specific Embodiments 1 to 6.

[0223] For the nonlinear problem of formula (26), the PI method can be used to solve it. The PI method can be divided into external iteration and internal iteration in the calculation process. In the internal iteration, the present invention uses the Krylov subspace method to replace the traditional source iteration method to calculate the internal iteration equation. In the external iteration, if the transport calculation result of this external iteration is determined to not meet the convergence condition, the CMFD acceleration method is executed once.

[0224] Specific implementation eight: This implementation is different from specific implementations one to seven in that the specific process of step 4 is as follows:

[0225] Step 401, initialize the number of outer iterations n=0, and initialize the eigenvalues ​​and solution vectors (using the solution vectors, the neutron angular flux results of the transport equation can be obtained):

[0226] k 0 =k0 (27)

[0227]

[0228] Among them, k 0 is the initialized eigenvalue, ψ 0 is the initialized solution vector, k0 and ψ0 are the initial values ​​and initial vectors given based on experience;

[0229] Step 402: Update the fission source Q of the neutron transport equation using the eigenvalue and solution vector n :

[0230]

[0231] Among them, |·| represents the modulus of the vector;

[0232] Step 403: According to the fission source Q n The inner iteration equation (Equation (30)) is established, and the Krylov subspace method is used to solve Equation (30):

[0233] Lψ n+1 =Q n (30)

[0234] Among them, ψ n+1 is the updated solution vector, and then according to ψ n+1 Update feature values:

[0235] k n+1 =|Fψ n+1 | (31)

[0236] Among them, k n+1 is the updated feature value;

[0237] Let the number of outer iterations of the PI method for solving the neutron transport equation be n = n + 1;

[0238] Step 404: According to ψ n Determine whether the convergence condition is met;

[0239] If the convergence condition is met, the calculation process ends and outputs ψ n and k n ;

[0240] Otherwise, continue to execute step 405;

[0241] Step 405: Using ψ n and k n To update the parameters of the coarse grid model, that is, to calculate equations (32) to (39), and then to construct a low-order diffusion approximation equation based on the updated coarse grid model parameters;

[0242] The specific process of step 405 is as follows:

[0243]

[0244]

[0245]

[0246]

[0247]

[0248]

[0249]

[0250]

[0251] Among them, V c (I) is the volume of coarse mesh unit I; First, use the latest ψ n Calculated by formula (11) and formula (9) Then Substituting into equation (33) the average neutron flux of fine mesh unit i in energy group g is calculated, where i represents the fine mesh unit covered by coarse mesh unit I (there is a strict nesting relationship between the coarse mesh and the fine mesh, that is, each coarse mesh unit just perfectly covers one or more fine mesh units, such as Figure 1 and Figure 2 shown); yes Initial value; is the neutron flux of coarse grid unit I in energy group g; V i is the volume of fine mesh unit i; is the total cross-section of the material at unit I in the coarse mesh under energy group g; Δx(I) is the width of coarse mesh unit I in the direction of the line connecting the centers of the two; Δx(J) is the width of coarse mesh unit J in the direction of the line connecting the centers of the two; yes Initial value; is the neutron flux size of the coarse mesh unit J under energy group g; is the neutron flux diffusion coefficient obtained based on finite differences; is the neutron flux correction factor obtained using the transport neutron flux, It's location The 0th-order Legendre moment of the scattering cross section from energy group g' to g at ;

[0252] The other symbols with “—” top line and “c” superscript are the physical parameters of the coarse mesh model. On the boundary surface, since the coarse mesh unit I has no adjacent units, and The calculation formula of will be different from equation (38) and equation (39). The specific form of the calculation formula is related to the boundary conditions, but this will not affect other calculation processes, so it will not be elaborated here.

[0253] In order to maintain conservation with the neutron transport equation, the intermediate variable Defined as:

[0254]

[0255]

[0256]

[0257] in, is the outward normal vector of the coarse mesh unit I on the common surface, is the external normal vector of the coarse mesh unit J on the common surface; S c (I, J) is the area of ​​the common surface of two adjacent coarse mesh units I and J;

[0258] The low-order diffusion approximation equation is constructed based on the updated coarse-grid model parameters and the finite difference method:

[0259]

[0260]

[0261]

[0262] in, is the net neutron flow from coarse mesh unit I to coarse mesh unit J; is the scattering cross section of the coarse-grid element I from energy group g' to g; is the neutron flux of coarse grid unit I in energy group g'; k c is the eigenvalue of the low-order diffusion approximation equation of the coarse grid model; is the neutron generation cross section of the material at coarse mesh unit I; is the probability that the energy group of the new neutron produced when the neutron undergoes fission reaction with the material at the coarse mesh unit I is exactly g;

[0263] Substituting equations (44) and (45) into equation (43), we obtain the finite difference discrete model of the coarse grid model:

[0264]

[0265]

[0266]

[0267]

[0268] in, is the vector consisting of the neutron fluxes of the coarse-grid units of energy group g, The first element in is the neutron flux size of the coarse mesh unit I, is the vector consisting of the sub-scalar fluxes of the coarse-grid units of energy group g′, yes The element in row I and column J in yes The element in row I and column J in yes The element in the I-th row and J-th column, J adjacent to I means that the coarse mesh unit I and the coarse mesh unit J are adjacent;

[0269] The finite difference discrete model is written in the form of equation (50) (similar to equation (26)):

[0270]

[0271]

[0272]

[0273]

[0274]

[0275] Step 406: Since the equation form of equation (50) is consistent with that of equation (26), the PI iteration process similar to that of equation (26) can be used to solve equation (50), and the Krylov method is also used in the inner iteration of PI:

[0276] Initialize the inner iteration number l′=0, and initialize the neutron flux vector and eigenvalue of the low-order diffusion approximation equation:

[0277] k c,0 =k n (55)

[0278]

[0279]

[0280] Step 407: Use the coarse-grid neutron flux vector Φ c,l' And the fission source term on the right side of the eigenvalue update equation (50):

[0281]

[0282] Step 408: c and k c To update:

[0283] L c Φ c,l′+1 =Q c,l′ (59)

[0284] k c,l′+1 =|F c Φ c,l′+1 | (60)

[0285] Among them, Φ c,l′+1 is the updated coarse-grid neutron flux vector, k c,l′+1 is the updated feature value;

[0286] Equation (59) is similar to the linear equation system of equation (30), and can also be solved using the Krylov method. However, because L c The form is relatively simple, and L c The dimension size is much smaller than L, so the matrix L can be directly calculated c and vector v c The product of gets the residual vector required by the Krylov solver.

[0287] Step 409, set the number of PI method outer iterations of the coarse grid model to l′=l′+1;

[0288] Determine whether the convergence condition is met, if the convergence condition is met, execute step 410, if not, return to execute step 407;

[0289] The convergence condition in this step is similar to the convergence condition in step 404, namely:

[0290]

[0291] Step 410: Calculate ψ n+1 / 2 and k n+1 / 2 :

[0292]

[0293] k n+1 / 2 =k c,final (62)

[0294] in, and k c,final are the neutron flux and eigenvalue of the coarse grid model after the PI method converges. is the position in the area covered by coarse network unit I The neutron angular flux of the nth outer iteration of the transport calculation at ; is the neutron angular flux after acceleration using the convergence results of the coarse grid model PI method; k n+1 / 2 It is the eigenvalue after the convergence result of the PI method using the coarse grid model is accelerated;

[0295] Depend on have

[0296]

[0297] because is the vector ψ n+1 / 2 The elements in the above formula can be used to obtain the accelerated ψ n+1 / 2 ;

[0298] Step 411: Using ψ n+1 / 2 and k n+1 / 2 Return to step 402.

[0299] Using ψ n+1 / 2 and k n+1 / 2 Updating the fission source term in the transport calculation in step 402 can have the effect of accelerating the iterative process of the neutron transport calculation.

[0300] The other steps and parameters are the same as those in Specific Embodiments 1 to 7.

[0301] Formula (30) is a system of linear equations, which can be written as the following more common linear equation system form:

[0302] Ax=b (64)

[0303] Among them, for the transport equation:

[0304] A=L (65)

[0305] x=ψ n+1 (66)

[0306] b=Q n (67)

[0307] In order to solve equation (64) using the Krylov subspace method, it is necessary to provide the calculated result of Av to the Krylov solver, where v is an intermediate vector with the same length as x generated by the Krylov method during the solution process. If T1, T2, T3, R and S are used to represent the behaviors represented by equations (3) to (7), then for the transport equation:

[0308] L=T1+T2+T3+R+S (68)

[0309] Av=Lv=T1v+T2v+T3v+Rv+Sv (69)

[0310] The following describes how to implement T1, T2, T3, R, S and F. Taking operator R as an example, substituting equation (11) into equation (6) yields:

[0311]

[0312] Since the matrix system (Eq. (26)) is obtained by continuously transforming the test function in Eq. (1) So any row R in the matrix R y Operations on vector v (R y v) is equivalent to: When the test function is When , the element v in the vector v is g,m,i As ψ g,m,i Substituting into equation (70) for calculation, as shown in equations (71) and (72), the subscripts "g, m, i" correspond to the position of a certain element of vectors v and ψ, but the specific element position depends on the operation when assembling the matrix system:

[0313]

[0314]

[0315] The other operators are similar, except that the calculation formula substituted into them is changed from formula (6) to one of formulas (3) to (5) and (7) to (8) (or the derivatives of these formulas after substituting boundary constraints and material constraints, such as formula (15) and formula (16). Therefore, the product operations of matrix operators such as T1, T2, T3, R, S and F with vector v can be replaced by integral calculations such as formula (72), thereby avoiding the storage of these complex matrices.

[0316] In practical applications, due to different physical conditions, the various items in equation (69) may be slightly different in different problems during actual execution. Among them, the most typical is that due to different boundary conditions, the calculation behavior corresponding to T1v is composed of equation (15), equation (16) or other unlisted boundary condition calculation formulas independently (when there is only one boundary condition) or on different boundary surfaces (when there are multiple boundary conditions). Because the Krylov solver generally only requires the calculation result of Av, and does not care about the specific calculation process, a modular programming method can be used, that is, various integral calculation behaviors that may be included are modularly programmed, such as implementing equations (15) and (16), equations (24) and (25) in different calculation modules. In actual calculations, the program will call different calculation modules for instantiation in different material regions and boundary surfaces according to the requirements of different physical conditions, so as to implement the functions of T1, T2, T3, R, S and F in different regions and boundary surfaces, and then accumulate the calculation results of these calculation modules to obtain the result vector of Av. Moreover, when there are new functional requirements (corresponding to a new integral calculation formula), you only need to write a new integral module based on the new integral calculation formula, and then introduce the module into the calculation to add new functions. This not only improves the functional extensibility of the program, but also increases the code reuse rate. A large number of integral modules can be reused in different problems.

[0317] Specific implementation method 9: This implementation method is different from any one of specific implementation methods 1 to 8 in that the convergence condition in step 404 is:

[0318]

[0319]

[0320] Among them, ε ψ is the convergence limit of the solution vector, ε k is the convergence limit of the eigenvalue of the transport equation.

[0321] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0322] The above calculation examples of the present invention are only used to explain the calculation model and calculation process of the present invention in detail, and are not intended to limit the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still within the scope of protection of the present invention.

Claims

1. A calculation method for the neutron transport equation of a nuclear reactor, characterized in that: The method specifically comprises the following steps: Step 1: read the coarse grid model, fine grid model, material information and boundary conditions of the nuclear reactor, and construct a grid mapping relationship between the coarse grid model and the fine grid model; Step 2: Establish a numerical discrete model of the steady-state multi-group transport equation based on the fine grid model read in step 1; Step 3, after applying boundary constraints and material constraints to the established numerical discrete model of steady-state multi-group transport equations according to material information and boundary conditions, a discrete equation group is constructed by transforming the test function in the numerical discrete model; Step 4: Use the PI method and the CMFD method, and solve the constructed discrete equations based on the grid mapping relationship between the coarse grid model and the fine grid model to obtain the neutron angular flux value and eigenvalue of the neutron transport equation.

2. A calculation method for the neutron transport equation of a nuclear reactor according to claim 1, characterized in that: The numerical discrete model of the steady-state multi-group transport equation is established based on the discrete ordinate discontinuous finite element method.

3. A calculation method for the neutron transport equation of a nuclear reactor according to claim 2, characterized in that: The numerical discrete model of the steady-state multi-group transport equation is specifically: in, represents the migration behavior of neutrons, represents the neutron consumption behavior, represents the scattering behavior of neutrons, Represents the fission behavior of neutrons.

4. A calculation method for the neutron transport equation of a nuclear reactor according to claim 3, characterized in that: The migration behavior of the neutron Neutron consumption behavior Neutron scattering behavior and the fission behavior of neutrons Specifically: in, represents the migration of neutrons on the boundary surface, represents the migration of neutrons between adjacent units, represents the migration of neutrons within the cell; in, is the spatial position coordinate; is the jth test function; is the gradient operator; D is the set of all units of the fine grid model, is the collection of fine mesh boundary cell faces of the nuclear reactor; in the collection superior, yes The outward normal vector of the fine mesh cell surface at ; on the set Γ of the fine mesh internal cell surfaces of the nuclear reactor, yes Satisfaction The element face normal vector of ; and w m They are the discrete quadrature groups in the discrete ordinate method. The discrete angles and the angle integral weights corresponding to the discrete angles, M is the number of discrete angles, Indicates the energy group g and the discrete direction The magnitude of the neutron angular flux under ; Represented on the internal element faces of the fine mesh The interruption, is the neutron angular flux of the upstream unit in the direction of neutron motion, in, Is in the energy group g, position The total cross-section of the material at; L is the order of Legendre expansion, l = 0, 1, ..., L; Is in the energy group g, position The neutron generation cross section of the material at is in position The probability that the energy group of the new neutron produced when the neutron undergoes fission reaction with the material is exactly g; It's location The l-th-order Legendre moment of the scattering cross section from energy group g' to g at , g' = 1, 2, ..., G, G is the total number of energy groups; is the discrete angle The l-order p-th spherical harmonic function, p = -l, -l+1, ..., l; is in the energy group g′ and the discrete direction The magnitude of the neutron angular flux under m is the discrete angle The corresponding angle integral weight; It's location The magnitude of the neutron flux at energy group g'; is the i′th basis function in element e in the discontinuous finite element method, ψ g,m,i′ yes is the expansion coefficient of , and N is the number of basis functions in unit e.

5. A calculation method for the neutron transport equation of a nuclear reactor according to claim 4, characterized in that: The specific process of applying boundary constraints to the numerical discrete model of the steady-state multi-group transport equation is as follows: make is the reflection boundary, is a vacuum boundary, then For vacuum boundaries The incident neutron flux is 0, that is, when hour Then we have: For reflective boundaries When two discrete directions and satisfy hour, and exist There is a reflection relationship about the reflection boundary; like but yes The reflection angle is like but yes The reflection angle is Then we have:

6. A calculation method for the neutron transport equation of a nuclear reactor according to claim 5, characterized in that: The specific process of applying material constraints to the numerical discrete model of the steady-state multi-group transport equation is as follows: Let D = D1∪D2, and after applying material constraints on material D1, become become become After applying material constraints on material D2, become become become but Let D = D aniso ∪D iso , D aniso represents the material region where anisotropic scattering is considered and the Legendre expansion order is L, D iso represents the material area where anisotropic scattering is not considered, then: In the formula, is the case when only isotropic scattering is considered. Abbreviation of .

7. A method for calculating the neutron transport equation for a nuclear reactor according to claim 6, characterized in that: The matrix operator form of the discrete equation system is: Where L includes the behaviors represented by equations (3) to (7); F represents the fission behavior of neutrons; ψ is ψ in equation (11) g,m,i′ The vector composed.

8. A calculation method for the neutron transport equation of a nuclear reactor according to claim 7, characterized in that: The specific process of step 4 is as follows: Step 401, initialize the number of outer iterations n=0, and initialize the eigenvalues ​​and solution vectors: k 0 =k0 (27) Among them, k 0 is the initialized eigenvalue, ψ 0 is the initialized solution vector; Step 402: Update the fission source Q of the neutron transport equation using the eigenvalue and solution vector n : Among them, · represents the modulus of the vector; Step 403: According to the fission source Q n The inner iteration equation is established, and the Krylov subspace method is used to solve equation (30): Lψ n+1 =Q n (30) Among them, ψ n+1 is the updated solution vector, and then according to ψ n+1 Update feature values: k n+1 =Fψ n+1 (31) Among them, k n+1 is the updated feature value; Let the number of outer iterations of the PI method for solving the neutron transport equation be n = n + 1; Step 404: According to ψ n Determine whether the convergence condition is met; If the convergence condition is met, the calculation process ends and outputs ψ n and k n ; Otherwise, continue to execute step 405; Step 405: Using ψ n and k n To update the parameters of the coarse grid model, that is, to calculate equations (32) to (39), and then to construct a low-order diffusion approximation equation based on the updated coarse grid model parameters; The specific process of step 405 is as follows: Among them, V c (I) is the volume of coarse mesh unit I; is the average neutron flux size of fine mesh unit i in energy group g, i represents the fine mesh unit covered by coarse mesh unit I; yes Initial value; is the neutron flux of coarse grid unit I in energy group g; V i is the volume of fine mesh unit i; is the total cross-section of the material at unit I in the coarse mesh under energy group g; Δx(I) is the width of coarse mesh unit I in the direction of the line connecting the centers of the two; Δx(J) is the width of coarse mesh unit J in the direction of the line connecting the centers of the two; yes Initial value; is the neutron flux size of the coarse mesh unit J under energy group g; is the neutron flux diffusion coefficient obtained based on finite differences; is the neutron flux correction factor obtained using the transport neutron flux, It's location The 0th-order Legendre moment of the scattering cross section from energy group g' to g at ; The intermediate variable Defined as: in, is the outward normal vector of the coarse mesh unit I on the common surface, is the external normal vector of the coarse mesh unit J on the common surface; S c (I, J) is the area of ​​the common surface of two adjacent coarse mesh units I and J; The low-order diffusion approximation equation is constructed based on the updated coarse-grid model parameters and the finite difference method: in, is the net neutron flow from coarse mesh unit I to coarse mesh unit J; is the scattering cross section of the coarse-grid element I from energy group g' to g; is the neutron flux of coarse grid unit I in energy group g'; k c is the eigenvalue of the low-order diffusion approximation equation of the coarse grid model; is the neutron generation cross section of the material at coarse mesh unit I; is the probability that the energy group of the new neutron produced when the neutron undergoes fission reaction with the material at the coarse mesh unit I is exactly g; Substituting equations (44) and (45) into equation (43), we obtain the finite difference discrete model of the coarse grid model: in, is the vector consisting of the neutron fluxes of the coarse-grid units of energy group g, The first element in is the neutron flux size of the coarse mesh unit I, is the vector consisting of the sub-scalar fluxes of the coarse-grid units of energy group g′, yes The element in row I and column J in yes The element in row I and column J in yes The element in the I-th row and J-th column, J adjacent to I means that the coarse mesh unit I and the coarse mesh unit J are adjacent; The finite difference discrete model is written in the form of equation (50): Step 406, initialize the inner layer iteration number l′=0, and initialize the neutron flux vector and eigenvalue of the low-order diffusion approximation equation: k c,0 =k n (55) Step 407: Use the coarse-grid neutron flux vector Φ c,l ' and the fission source term on the right side of the eigenvalue update equation (50): Step 408: c and k c To update: 50 c Φ c,l′+1 =Q c,l′ (59) k c,l′+1 =F c F c,l′+1 (60) Among them, Φ c,l′+1 is the updated coarse-grid neutron flux vector, k c,l′+1 is the updated feature value; Step 409, set the number of PI method outer iterations of the coarse grid model to l′=l′+1; Determine whether the convergence condition is met, if the convergence condition is met, execute step 410, if not, return to execute step 407; Step 410: Calculate ψ n+1 / 2 and k n+12 : k n+12 =k c,final (62) in, and k c,final They are the neutron flux size and eigenvalue size after the PI method of the coarse grid model converges; Step 411: Using ψ n+1 / 2 and k n+12 Return to step 402.

9. A calculation method for the neutron transport equation of a nuclear reactor according to claim 8, characterized in that: The convergence condition in step 404 is: Among them, ε ψ is the convergence limit of the solution vector, ε k is the convergence limit of the eigenvalue of the transport equation.

Images