Numerical calculation method for solving neutron transport equation

Abstract

The invention relates to the field of reactor data processing, and particularly discloses a numerical calculation method for solving a neutron transport equation. The method comprises the following steps: 1, describing the geometry of a reactor core of a nuclear reactor by utilizing a construction entity geometric modeling method, and subdividing a curved surface grid along a curved surface structure; 2, directly expanding a neutron flux distribution function in each curved surface grid by using a multivariate polynomial; 3, deriving a neutron transport equation and an integral residual equation of a boundary condition under each curved surface grid; 4, solving geometric integral parameters of each curved surface grid; step 5, establishing an adaptive linear system; and step 6, iteratively calculating and solving a maximum characteristic value and an unknown expansion coefficient.

Description

Technical Field

[0001] This invention relates to the field of reactor data processing, and in particular to a numerical calculation method for solving neutron transport equations. Background Technology

[0002] In real reactors, numerous curved surfaces exist, such as rod-shaped fuel elements, thin-layered flammable poison (IFBA) on cylindrical fuel pellets, and cylindrical containers. How to realistically handle these curved surfaces, especially accurately performing surface integral calculations, is a major challenge for deterministic methods and a primary factor limiting their computational accuracy. Currently, the main deterministic solution techniques that can directly adapt to curved meshes include the Method of Characteristics (MOC) and some Discrete Ordinate (DO) methods based on cylindrical or spherical coordinates. Additionally, there are solution techniques that approximate curved structures using irregular meshes, such as the Variational Nodal Method (VNM), the Discrete Ordinates Transport Method (DNTM), and the Finite Element Method (FEM).

[0003] The method of characteristics (MOC) is an approximate numerical method based on characteristic theory for solving hyperbolic partial differential equations. It's an integral method for solving transport equations. Its core idea is to simultaneously discretize space and angles, integrating the differential form of the neutron transport equation along the trajectory of neutron transport (i.e., characteristic lines), transforming it into a problem of solving a large number of one-dimensional parallel characteristic line equations. Theoretically, as long as the characteristic lines are arranged densely enough to cover the entire computational domain, MOC can be applied to solve neutron transport equations for any geometric problem.

[0004] Discrete ordinate methods based on cylindrical or spherical coordinates can naturally and accurately describe concentric cylindrical or spherical surfaces by dividing the grid in polar coordinates and establishing the difference relationship between the derivative terms of the grids.

[0005] The solution technique for approximating curved surface structures using irregular meshes typically discretizes the space into many triangular or polygonal meshes, and then uses variational principles or weighted residual methods to construct the residual constraints of the governing equations and boundary condition equations under these mesh integrals. The unknown field function is approximated as a polynomial expansion, and the basis functions of the expansion are usually selected from interpolation basis functions such as Lagrange, so that the unknown quantity corresponds to the neutron flux value at each interpolation point (the endpoint of the irregular mesh).

[0006] While the method of characteristics is widely used in high-fidelity physics calculations, it requires repeated online scanning calculations, resulting in high computational costs and multiple iterations. Although acceleration methods have made significant progress and effectively reduced the number of MOC iterations, acceleration bottlenecks still exist (Zhu, K., Kong, B., Zhang, H., Guo, J., & Li, F. (2022). High-Fidelity Neutron Transport Solution of High Temperature Gas-Cooled Reactor by Three-Dimensional Linear Source Method of Characteristics. Nuclear Science and Engineering, 197(6), 1174–1196. https: / / doi.org / 10.1080 / 00295639.2022.2143706). Furthermore, when dealing with regions with fine geometric distributions and drastic changes in cross-section, such as thin arc-shaped fuels and control drum absorbers, the characteristic line spacing needs to be locally reduced to maintain computational accuracy, leading to a surge in computation time by tens of times or more. Moreover, it is difficult to implement and lacks versatility.

[0007] The characteristic line method essentially discretizes space into a cluster of parallel characteristic lines with a certain spacing in different spatial angles. The intersection points of these characteristic lines with the boundaries of any shaped space can be easily calculated, thus determining a characteristic line segment. The rectangle formed by multiplying the length of each characteristic line segment by its spacing can fill any shaped space, and the sum of the areas of these rectangles can approximate the area of ​​any shaped space. However, if the spacing between the characteristic line clusters is too large, errors in area calculation will occur. Furthermore, if a spatial region is particularly small, excessively large spacing between the characteristic line clusters, resulting in no intersection points, will also cause significant calculation errors. Approximation methods using irregular meshes suffer from limitations because the boundaries of irregular meshes are planar, and the Lagrange interpolation basis functions can only describe the function distribution between interpolation point connections, not the curvature of the surface. Therefore, approximation can only be achieved by increasing the number of interpolation points on the surface or by boundary correction.

[0008] Discrete ordinate methods based on cylindrical or spherical coordinates can only describe concentric cylindrical or spherical surfaces due to the limitations of polar coordinates. This limits their application to highly homogenized simplified models and makes them unsuitable for high-fidelity computation.

[0009] Solving techniques that approximate curved surfaces using irregular meshes require sufficiently fine mesh divisions to ensure computational accuracy, thus severely limiting computational efficiency. Furthermore, some discrete ordinate block methods under irregular meshes decompose the three-dimensional neutron transport equations into three one-dimensional equations, which are coupled together using transverse integration leakage techniques. The transverse integration typically employs a "flat leakage" approximation method, which also introduces errors.

[0010] In general, all existing solution techniques on the market have approximate errors in the description of surface geometry. Reducing these errors usually only requires refining the mesh, which drastically increases the computational cost. Summary of the Invention

[0011] As mentioned in the background section, existing methods, when faced with the problem of approximation errors in the geometric description of curved surfaces, typically employ simple and direct methods such as refining the mesh or proposing techniques for approximate mesh boundary corrections until an acceptable level is achieved. This is because deterministic methods mostly involve the integral calculation of neutron flux and neutron flow over mesh elements and boundary surfaces. However, the integration domain of curved surface mesh elements is irregular, and the normal directions on non-curved surfaces are not fixed, making it difficult to perform integral calculations using conventional methods. This limits the direct application of currently available methods or theories to curved surface meshes unless a new solution technique specifically designed for curved surface meshes is proposed, which usually requires a deep understanding and experience with the neutron transport equations and their numerical solution techniques.

[0012] The technical solution provided by this invention is as follows:

[0013] This invention provides a numerical calculation method for solving the neutron transport equation, which includes the following steps:

[0014] Step 1: Describe the nuclear reactor core geometry using the Construct Solid Geometry (CSG) method and mesh the surface along the curved structure;

[0015] Step 2: Expand the neutron flux distribution function directly using a multivariate polynomial in each surface mesh;

[0016] Step 3: Derive the integral residual equations of neutron transport equations and boundary conditions for each surface mesh;

[0017] Step 4: Solve for the geometric integral parameters of each surface mesh;

[0018] Step 5: Establish the solution system, that is, establish a well-posed linear system;

[0019] Step 6: Iteratively calculate and solve for the largest eigenvalue and the unknown expansion coefficients.

[0020] Step 2 specifically involves expanding the neutron angular flux density and source term directly in space using a two-dimensional polynomial, without the need for transverse integral coupling. The angular flux is expanded using a 5th-order polynomial, and the source term is expanded using a 3rd-order polynomial, as shown in equation (1).

[0021]

[0022] In the formula: P k (x,y) are basis functions for the spatial expansion, P0=1, It is the material center of each segment in a coordinate system with the material center of the cell as the origin. When dealing with problems involving uniform material segments, it is also their geometric center; ψ gik S gik These are the expansion coefficients of the basis functions, and are parameters independent of geometric position.

[0023] In step 3, the weighted integral of the residuals of the neutron transport equation and the boundary condition equation within each surface grid is set to 0, and the resulting solution is considered to be an approximate solution to the equation.

[0024] Step 3 specifically involves:

[0025] (1) For the two-dimensional steady-state problem, considering isotropic scattering, the neutron transport equation after energy and angle discretization is obtained, as shown in equation (2):

[0026]

[0027] In the formula: i represents the discrete direction; g represents the energy group; Ω xyi =Ω xi +Ω yi ψ is the spatial angle in a two-dimensional Cartesian coordinate system. gi S represents the angular flux density of neutrons after directional discretization. gi The source term after directional discretization includes fission source term and scattering source term, without considering external sources;

[0028] (2) After discretizing the neutron transport equation in terms of angle and space, perform a weighted integral in each discrete block, and take the form of the weighting function as the same as the expanded basis function (Gallenkin weighted residual method). Then the k′th order weighted residual equation is equation (3):

[0029]

[0030] In particular, the 0th order weighted residual equation is the block integral neutron equilibrium equation, i.e., equation (4):

[0031]

[0032] (3) By defining geometric parameters to simplify the expression of the integral, the k′-th order weighted residual equation is simplified to equation (5):

[0033]

[0034] The geometric parameters are shown in equation (6):

[0035]

[0036] The zeroth-order weighted residual equation can be simplified to equation (7):

[0037]

[0038] In the formula: This is the integral leakage term of neutrons on the m-plane, denoted as L. gim ; That is, the moment of angular flux;

[0039] The integral is defined as in equation (8):

[0040]

[0041] (4) Similarly, the boundary condition is that the neutron flux distribution at the incident surface satisfies the boundary condition, which is Equation (9):

[0042] ψ gi (x,y)-ψ BC (x,y)=0 (9)

[0043] Weighted integration within each discrete block yields the k′-th order weighted residual equation, which is equation (10):

[0044]

[0045] Step 4 is as follows:

[0046] For relatively regular integration regions, the integral prototype is obtained and then the upper and lower limits are substituted for calculation. For integration regions that are more difficult to calculate, numerical calculation methods are used to solve the problem, such as interpolation quadrature formulas or Monte Carlo methods for approximate solutions.

[0047] Step 5 specifically involves:

[0048] For cases where the model boundary conditions or continuity conditions with neighboring nodes need to be satisfied on the incident boundary surface, the weighted residual equation needs to be selected.

[0049] More specifically, the solution system is established as follows:

[0050] For the weighted residual equations of the neutron transport equations within each block, take the three weighted residual equations corresponding to k′=0,1,2 in equation (5); for the weighted residual equations of the boundary conditions on the boundary of each block, consider the case that each block has 4 boundary surfaces, and on average each block has 2 incident surfaces and 2 exit surfaces. Therefore, take the one weighted residual equation corresponding to k′=0 in equation (10) on each incident surface (i.e., the continuity condition of the integral leakage term). In this way, the five expansion coefficients of each energy group of each block correspond to five linearly independent weighted residual equations, and determine the unique solution.

[0051] By combining the five constraint equations for all angles, all energy groups, and all nodes, a large linear system (MS)ψ = 1 / kFψ is formed, where M is the output operator, S is the scattering source operator, F is the fission source operator, and k is the effective multiplication factor.

[0052] Step 6 specifically involves:

[0053] For solving eigenvalue problems, the CSN method adopts the source iteration method, in which the scattering source term is placed on the left side of the equation, and neutrons are artificially divided into different generations. The neutrons of each generation are generated from fission as the starting point, and each source iteration is the process of forming a new generation of neutrons from the previous generation of neutrons after a series of reactions.

[0054] The ratio of the number of neutrons produced by fission in two adjacent iterations is the effective multiplication factor, i.e., equation (11):

[0055]

[0056] in c is the node number, and N is the total number of discrete nodes;

[0057] The source iteration ends when the effective multiplication factor and flux are iterated to meet the convergence criterion, and the eigenvalue problem is solved.

[0058] This invention provides a solution model for the neutron transport equation, which includes a data input module, a data processing module, and a result output module, wherein the data processing module is used to execute the calculation method described above.

[0059] The present invention also provides a non-transitory computer-readable storage medium storing at least one instruction or at least one program segment, said at least one instruction or said at least one program segment being loaded and executed by a processor to implement the method or the model.

[0060] The present invention also provides an electronic device including a processor and the non-transitory computer-readable storage medium.

[0061] The beneficial effects of this invention are as follows: The core of this invention is based on direct multivariate polynomial expansion under a real curved mesh, a concept that differs from all other methods. Existing theories lack transverse integration approximation operations, and this invention demonstrates that it guarantees the neutron conservation condition within the mesh, facilitating the implementation of low-order acceleration methods.

[0062] Specifically, this invention proposes a novel numerical calculation method for neutron transport based on curved surface meshes (i.e., boundaries can be planar or arc surfaces), namely the Curved Surface Nodal method (CSN). This method discretizes neutrons angularly using a discrete ordinate method and spatially divides them directly according to the material composition using a curved surface mesh. The angular flux in each direction is directly expanded in multiple dimensions within the curved surface mesh, and an explicit relationship is established between the interface integral average and the spatial distribution of the internal neutron flux. Then, using weighted residual equations and continuity conditions, a global simultaneous equation is constructed to solve for the expansion coefficients, directly yielding the spatial distribution of angular and ordinate fluxes. Compared to the discrete ordinate method, the spatial solution method of this invention is completely different and directly supports surface geometry. Compared to the MOC method, the spatial solution method of this invention strictly supports surface geometry and is not constrained by the spacing of feature line clusters. Compared to the finite element method, the spatial element partitioning of this method directly supports surface geometry, and the continuity condition between elements is changed from point continuity to surface-average continuity, which can better satisfy the neutron conservation condition within the element. Compared to the nodal method, it is not limited by regular geometric partitioning and is not limited by the additional approximations caused by transverse integration, and is expected to achieve higher accuracy. The technical solution of this invention has practical significance, as its calculation results involve a large number of surface geometries in real reactors, such as rod-shaped fuel elements, coolant channels, irradiation measurement channels, and cylindrical containers. At the same time, in recent years, advanced reactors, especially fast reactors, small modular reactors, and high-flux experimental reactors, have become the focus of global nuclear energy development, exhibiting a design trend of geometric singularities and strong non-uniformity. Core analysis methods based on regular grid partitioning ideas such as rectangles, hexagons, and modularity, which originated from traditional pressurized water reactors, are no longer applicable. There is an urgent need for a solution that can adapt to the requirements of these complex geometries. Attached Figure Description

[0063] Figure 1 The technical flowchart of this invention.

[0064] Figure 2 Geometric model of curved mesh problem.

[0065] Figure 3 Standard flux distribution for surface mesh problems. Among them, (a) standard flux distribution for problem A, (b) standard flux distribution for problem B, and (c) standard flux distribution for problem C. Detailed Implementation

[0066] Example 1: Method Steps of the Invention

[0067] The technical solution of the present invention is as follows Figure 1 As shown.

[0068] Step 1: Describe the reactor core geometry using Constructive Solid Geometry (CSG) and mesh the surface along the curved structure. CSG is the preferred method for many advanced modeling software programs (such as CAD) due to its structured modeling characteristics, which reduce the storage of geometric information and offer greater flexibility. The most important and commonly used technique in CSG is Boolean operation. Simply put, Boolean operation is a three-dimensional set concept. It includes union, complement, and intersection. A union (Boolean union) is the combination of two objects. An intersection (Boolean intersection) is the common part of two objects. Complement (Boolean intersection) is a relative concept, subtracting object B from object A. By specifying different A and B and applying different Boolean operations, arbitrarily complex geometries can be easily obtained.

[0069] Step 2: Expand the neutron flux distribution function directly using a multivariate polynomial in each surface mesh.

[0070] For spatial variables, a constructive polynomial fitting method is adopted. Unlike the traditional transport nodal method, the CSN method directly expands the neutron angular flux density and source term in space using a two-dimensional polynomial without the need for lateral integral term coupling. The angular flux adopts a 5th-order expansion and the source term adopts a 3rd-order expansion, as shown in Equation (1):

[0071]

[0072] In the formula: P k (x,y) are basis functions for the spatial expansion, P0=1, It is the material center of each segment in a coordinate system with the material center of the cell as the origin. When dealing with problems involving uniform material segments, it is also their geometric center; ψ gik S gik These are the expansion coefficients of the basis functions, and are parameters independent of geometric position.

[0073] Step 3: Derive the integral residual equations of neutron transport equations and boundary conditions for each surface mesh.

[0074] Since the approximate polynomial fitting distribution cannot accurately satisfy the differential equations and boundary conditions, i.e., there are residuals, according to the idea of ​​the method of weighted residuals (MWR), the weighted integral of the residuals of the neutron transport equation and boundary condition equation in each surface grid is set to 0, and the resulting solution is considered to be an approximate solution of the equation.

[0075] (1) For the two-dimensional steady-state problem, considering isotropic scattering, the neutron transport equation after energy and angle discretization can be obtained, as shown in equation (2):

[0076]

[0077] In the formula: i represents the discrete direction; g represents the energy group; Ω xyi =Ω xi +Ω yi ψ is the spatial angle in a two-dimensional Cartesian coordinate system. gi S represents the angular flux density of neutrons after directional discretization. gi The source term is the directionally discretized term, which includes fission source term and scattering source term, and does not consider external sources.

[0078] (2) After discretizing the neutron transport equation in terms of angle and space, perform a weighted integral in each discrete block, and take the form of the weighting function as the same as the expanded basis function (Gallenkin weighted residual method). Then the k′th order weighted residual equation is equation (3):

[0079]

[0080] In particular, the 0th order weighted residual equation is the block integral neutron equilibrium equation, i.e., equation (4):

[0081]

[0082] (3) Since the integral in the solution domain (each surface grid block) space is only related to the coordinates, and the part of the flux that is related to the coordinates has been represented by a known finite-order polynomial expansion, the remaining flux expansion coefficients are independent of the integral. Therefore, the integral can be solved relatively easily and can be simplified by defining geometric parameters. The k′-th order weighted residual equation can be simplified to equation (5):

[0083]

[0084] The geometric parameters are shown in equation (6):

[0085]

[0086] The zeroth-order weighted residual equation can be simplified to equation (7):

[0087]

[0088] In the formula: This is the integral leakage term of neutrons on the m-plane, denoted as L. gim ; That is, the moment of angular flux.

[0089] The integral is defined as in equation (8):

[0090]

[0091] (4) Similarly, the boundary condition is that the neutron flux distribution at the incident surface satisfies the boundary condition, which is Equation (9):

[0092] ψ gi (x,y)-ψ BC (x,y)=0 (9)

[0093] Weighted integration within each discrete block yields the k′-th order weighted residual equation, which is equation (10):

[0094]

[0095] Step 4: Solve for the geometric integral parameters of each surface mesh.

[0096] As can be seen from equations (6) and (8), the geometric parameters b, c, and d are all single or double integrals of simple expansion basis functions, and are independent of physical parameters involved in iterative updates, such as neutron flux, source terms, and neutron flux. They can be directly calculated after the model is determined, and can be directly called during the iteration process without repeated calculations. For relatively regular integration regions, the integral prototype can be easily obtained and then calculated by substituting the upper and lower limits. For integration regions that are more difficult to calculate, numerical calculation methods, such as interpolation quadrature formulas or Monte Carlo methods, can be used for approximate solutions. The specific implementation steps will be introduced below using the integral prototype method as an example.

[0097] (1) The integral calculation on the boundary surface can be expressed as the first kind of line integral of the curve multiplied by the height. The formula for the first kind of line integral can be expressed as equation (11):

[0098]

[0099] in y = ψ(t), α ≤ t ≤ β.

[0100] (2) For the calculation of the integral within the node, it can be expressed as a double integral multiplied by the height. The calculation methods of double integrals are mainly classified into two categories: one is to calculate the double integral using rectangular coordinates, which is suitable for integration regions that conform to the X-shape or Y-shape; the other is to calculate the double integral using polar coordinates, which is suitable for integration regions that can be represented by polar coordinates, such as sectors or circles. In specific implementation, the corresponding method can be flexibly selected according to the shape characteristics of the node. Among them, the formula for the X-shape double integral can be expressed as equation (12):

[0101]

[0102] The formula for a Y-type double integral can be expressed as equation (13):

[0103]

[0104] The formula for the double integral in polar coordinates can be expressed as equation (14):

[0105]

[0106] Step 5: Establish the solution system.

[0107] Since the angular flux is expanded using 5th-order basis functions, there are 5 expansion coefficients to be determined for each angle, each node, and each energy group, requiring the identification of 5 constraint equations. As mentioned earlier, a 5th-order weighted residual equation can theoretically be obtained, but the weighted residual method can only obtain an approximate solution distribution within the integration region and cannot guarantee the boundary conditions. Therefore, it is necessary to supplement the boundary condition constraints, i.e., the model boundary conditions or continuity conditions with neighboring nodes must be satisfied on the incident boundary surface. In this case, the weighted residual equation needs to be selected. Usually, the first consideration is to remove higher-order terms. This patent application proposes a solution system establishment scheme:

[0108] For the weighted residual equations of the neutron transport equations within each block, take the three weighted residual equations corresponding to k′=0,1,2 in equation (5); for the weighted residual equations of the boundary conditions on the boundary of each block, consider the case that each block has 4 boundary surfaces, and on average each block has 2 incident surfaces and 2 exit surfaces. Therefore, take the one weighted residual equation corresponding to k′=0 in equation (10) on each incident surface (i.e., the continuity condition of the integral leakage term). In this way, the five expansion coefficients of each energy group of each block correspond to five linearly independent weighted residual equations, and a unique solution can be determined.

[0109] By simultaneously solving the five constraint equations for all angles, all energy groups, and all nodes, a large linear system (MS)ψ = 1 / kFψ can be formed. Here, M is the output operator, S is the scattering source operator, F is the fission source operator, and k is the effective multiplication factor.

[0110] Step 6: Iteratively calculate and solve for the largest eigenvalue and the unknown expansion coefficients.

[0111] For solving the eigenvalue problem, the CSN method adopts the source iteration method, in which the scattering source term is placed on the left side of the equation. The basic principle is to artificially divide neutrons into different generations, with each generation of neutrons starting from the time of fission. Each source iteration is the process of forming a new generation of neutrons from the previous generation of neutrons through a series of reactions. The ratio of the number of neutrons produced by fission in two adjacent iterations is the effective multiplication factor, i.e., equation (15):

[0112]

[0113] in c is the node number, and N is the total number of discrete nodes.

[0114] The source iteration ends when the effective multiplication factor and flux are iterated to meet the convergence criterion, and the eigenvalue problem is solved.

[0115] Verification Example

[0116] The effectiveness of this invention can be verified through the following examples following the steps described in the embodiments:

[0117] Problems with one-dimensional boundaries as curved surfaces (Problem A), problems with one-dimensional interfaces as curved surfaces (Problem B), and cylindrical problems in two-dimensional Cartesian coordinates (Problem C) are presented with geometric models as follows: Figure 2 As shown. All cross-sections are Σ. t =1.0cm -1 , Σ a =0.6cm -1 ,νΣ f =0.2cm -1 , Σ s =0.4cm -1 The boundary condition on the plane is reflection, and the boundary condition on the curved surface is vacuum.

[0118] In this invention, to address the three problems mentioned above, three geometric models are first defined using the solid geometry construction method in step 1, including the origin of the coordinate system, surface equations, and surface half-spaces. Problem B involves two mesh regions; therefore, in addition to describing the two regions separately, it is also necessary to define the neighbor relationships between the meshes. The final geometric information is shown in Table 1.

[0119] Table 1 Geometric Information for the Surface Mesh Problem

[0120]

[0121]

[0122] Then, the polynomial expansion from step 2 is used to approximate the spatial distribution shape of the neutron angular flux in each region.

[0123] Then, according to steps 3-5, a well-posed linear system is established by introducing reasonable constraint equations.

[0124] Finally, the unknown expansion coefficients of the internal corner flux of each segment are solved using the source iteration method described in step 6.

[0125] To ensure computational accuracy, a discrete quadrature group decoupled from the azimuth and polar angles was adopted. The azimuth and polar angle discretization and weighting were set according to the OpenMOC program, and the Tabuchi-Yamamoto 3-polar-angle quadrature group was used for the polar angles. Calculations were performed for azimuth angles of 16, 24, 40, and 48 degrees, with a convergence criterion of 10 for both eigenvalues ​​and flux. -10 The calculation results are all derived from the OpenMC program. The calculation results are shown in Table 2. Figure 3 Standard flux distribution plots were drawn for the three problems.

[0126] Table 2 Calculation results for the surface mesh problem

[0127]

[0128] As shown in Table 2, the CSN method can calculate and solve problems with one-dimensional boundaries, one-dimensional interfaces, and cylinder problems in two-dimensional Cartesian coordinates according to the actual surface structure. It achieves high computational accuracy with a small number of surface meshes, and the computational accuracy increases with the number of discrete angles, which demonstrates the applicability of the CSN method proposed in this invention to surface mesh problems. Figure 3 This demonstrates that CSN can directly obtain the scalar flux distribution within the computational domain, even if the computational domain is an irregular grid. However, it should be noted that the surface in Problem B should satisfy the continuity condition, but the calculated scalar flux is not continuous everywhere. This is because the current boundary condition constraints are only limited to the integral average of the leakage term, which cannot guarantee that the boundary conditions and continuity conditions are satisfied everywhere, which may have a certain impact on the computational accuracy and stability.

Claims

1. A numerical calculation method for solving the neutron transport equation, characterized in that, Includes the following steps: Step 1: Describe the nuclear reactor core geometry using the Construct Solid Geometry (CSG) method and mesh the surface along the curved structure; Step 2: Expand the neutron flux distribution function directly using a multivariate polynomial in each surface mesh; Step 3: Derive the integral residual equations of neutron transport equations and boundary conditions for each surface mesh; Step 4: Solve for the geometric integral parameters of each surface mesh; Step 5: Establish the solution system, that is, establish a well-posed linear system; Step 6: Iteratively calculate and solve for the largest eigenvalue and the unknown expansion coefficients.

2. The numerical calculation method as described in claim 1, characterized in that, Step 2 specifically involves expanding the neutron angular flux density and source term directly in space using a two-dimensional polynomial, without the need for transverse integral coupling. The angular flux is expanded using a 5th-order polynomial, and the source term is expanded using a 3rd-order polynomial, as shown in equation (1). In the formula: P k (x,y) are basis functions for the spatial expansion, P0=1, It is the material center of each segment in a coordinate system with the material center of the cell as the origin. When dealing with problems involving uniform material segments, it is also their geometric center; ψ gik S gik These are the expansion coefficients of the basis functions, and are parameters independent of geometric position.

3. The numerical calculation method as described in claim 1, characterized in that, In step 3, the weighted integral of the residuals of the neutron transport equation and the boundary condition equation within each surface grid is set to 0, and the resulting solution is considered to be an approximate solution to the equation.

4. The numerical calculation method as described in claim 3, characterized in that, Step 3 specifically involves: (1) For the two-dimensional steady-state problem, considering isotropic scattering, the neutron transport equation after energy and angle discretization is obtained, as shown in equation (2): In the formula: i represents the discrete direction; g represents the energy group; Ω xyi =Ω xi +Ω yi ψ is the spatial angle in a two-dimensional Cartesian coordinate system. gi S represents the angular flux density of neutrons after directional discretization. gi The source term after directional discretization includes fission source term and scattering source term, without considering external sources; (2) After discretizing the neutron transport equation in terms of angle and space, perform a weighted integral in each discrete block, and take the form of the weighting function as the same as the expanded basis function (Gallenkin weighted residual method). Then the k′th order weighted residual equation is equation (3): In particular, the 0th order weighted residual equation is the block integral neutron equilibrium equation, i.e., equation (4): (3) By defining geometric parameters to simplify the expression of the integral, the k′-th order weighted residual equation is simplified to equation (5): The geometric parameters are shown in equation (6): The zeroth-order weighted residual equation can be simplified to equation (7): In the formula: This is the integral leakage term of neutrons on the m-plane, denoted as L. gim ; That is, the moment of angular flux; The integral is defined as in equation (8): (4) Similarly, the boundary condition is that the neutron flux distribution at the incident surface satisfies the boundary condition, which is Equation (9): ψ gi (x,y)-ψ BC (x,y)=0 (9) Weighted integration within each discrete block yields the k′-th order weighted residual equation, which is equation (10):

5. The numerical calculation method as described in claim 4, characterized in that, Step 4 specifically involves: For relatively regular integration regions, the integral prototype is obtained and then the upper and lower limits are substituted for calculation. For integration regions that are more difficult to calculate, numerical calculation methods are used to solve them, such as interpolation-type quadrature formulas or Monte Carlo methods for approximate solutions.

6. The numerical calculation method as described in claim 1, characterized in that, Step 5 specifically involves: For cases where the model boundary conditions or continuity conditions with neighboring nodes need to be satisfied on the incident boundary surface, the weighted residual equation needs to be selected. More specifically, the solution system is established as follows: For the weighted residual equations of the neutron transport equations within each block, take the three weighted residual equations corresponding to k′=0,1,2 in equation (5); for the weighted residual equations of the boundary conditions on the boundary of each block, consider the case that each block has 4 boundary surfaces, and on average each block has 2 incident surfaces and 2 exit surfaces. Therefore, take the one weighted residual equation corresponding to k′=0 in equation (10) on each incident surface (i.e., the continuity condition of the integral leakage term). In this way, the five expansion coefficients of each energy group of each block correspond to five linearly independent weighted residual equations, and determine the unique solution. By combining the five constraint equations for all angles, all energy groups, and all nodes, a large linear system (MS)ψ = 1 / kFψ is formed, where M is the output operator, S is the scattering source operator, F is the fission source operator, and k is the effective multiplication factor.

7. The numerical calculation method as described in claim 1, characterized in that, Step 6 specifically involves: For solving eigenvalue problems, the CSN method adopts the source iteration method, in which the scattering source term is placed on the left side of the equation, and neutrons are artificially divided into different generations. The neutrons of each generation are generated from fission as the starting point, and each source iteration is the process of forming a new generation of neutrons from the previous generation of neutrons after a series of reactions. The ratio of the number of neutrons produced by fission in two adjacent iterations is the effective multiplication factor, i.e., equation (11): in c is the node number, and N is the total number of discrete nodes; The source iteration ends when the effective multiplication factor and flux are iterated to meet the convergence criterion, and the eigenvalue problem is solved.

8. A solution model for the neutron transport equation, characterized in that, It includes a data input module, a data processing module, and a result output module, wherein the data processing module is used to perform the calculation method as described in any one of claims 1 to 7.

9. A non-transitory computer-readable storage medium, characterized in that, The storage medium stores at least one instruction or at least one program segment, which is loaded and executed by a processor to implement the method as described in any one of claims 1 to 7, or the model as described in claim 8.

10. An electronic device, characterized in that, Includes a processor and the non-transitory computer-readable storage medium of claim 9.

Images