A reactor transport criticality simulation method and system based on Krylov-Schur algorithm

Abstract

A reactor transport criticality simulation method and system based on the Krylov-Schur algorithm, relating to the field of nuclear reactor core design and safety, and specifically to a reactor transport criticality simulation method and system based on the Krylov-Schur algorithm. The purpose of the present invention is to realize the transport criticality calculation and solve the maximum eigenvalue using the Krylov-Schur method, and to improve the convergence accuracy of the maximum eigenvalue and the eigenvector. Technical points: The steady-state neutron transport equation is discretized and deformed in the form of an operator, the energy variable E in the steady-state neutron transport equation is approximated by multiple groups, the angle variable is discretized by vertical scale to obtain a coupled equation expressed in the form of an operator, a standard eigenvalue problem is given, and the eigenvalues ​​and eigenvectors of the matrix P in the standard eigenvalue problem are solved. The method of the present invention is used for reactor transport criticality simulation.

Description

Technical Field

[0001] The present invention relates to the fields of reactor transport criticality simulation, nuclear reactor core design and safety, and in particular to a reactor transport criticality simulation method and system based on a Krylov-Schur algorithm. Background Art

[0002] The core problem of reactor criticality calculation is to solve the maximum eigenvalue and corresponding eigenvector of the neutron transport equation. In addition, high-order eigenvalues ​​and eigenvectors play an important role in the analysis of reactor stability, neutron flux oscillation and other issues and in reactor monitoring. Traditionally, the power iteration (PI) (CN110717275A) method is used in the field of nuclear reactor physics to solve reactor critical parameters. However, when the dominance ratio is close to 1, the efficiency of the power iteration method decreases, and it can only solve the maximum eigenvalue, and cannot obtain high-order eigenvalues ​​and eigenvectors.

[0003] In recent years, researchers have begun to explore the use of more advanced numerical iterative methods to solve reactor criticality problems: Verdu et al. [1] The implicit restart Arnoldi method (IRAM) was used to solve multiple eigenvalues ​​of the neutron diffusion equation. The Three Mile Island nuclear power plant reactor was used for testing, and the performance of IRAM was tested against that of subspace iteration. The test found that compared with the subspace iteration method, the IRAM method showed better computational performance and was more reliable and stable, and was not affected by eigenvalue degradation. Warsa et al. [2] In the linear discontinuous finite element program, the IRAM eigenvalue calculation was implemented based on ARPACK, and IRAM was initialized using PI. The performance of the algorithm was verified in detail and compared with the performance of power iteration. A sensitivity analysis was performed on the parameters that may affect IRAM, such as the number of PIs used for initialization and the spatial dimension of the inner iteration Krylov iteration. They found that: compared with PI, IRAM has better performance and is more stable; for problems with a dominance ratio close to 1, IRAM can increase the calculation speed by 10-100 times; the larger the dominance ratio, the greater the advantage of IRAM over PI; the existence of upper scattering will reduce the performance advantage of IRAM over PI to a certain extent. However, all the tests by Warsa et al. were completed in a serial environment, and no detailed study was conducted on the performance of IRAM in a parallel environment. Nicolo Abrate et al. [3]The computational stability and efficiency of Filtered Power Method (FPM), IRAM, Sub-space Iteration (SSI) and other methods in reactor criticality calculations based on the diffusion equation were analyzed and studied. Various factors affecting the efficiency of iterative methods, such as the selection of initial iteration vectors, subspace dimensions, and convergence criteria, were analyzed in detail. Finally, they pointed out that considering various factors such as computational stability, robustness, ease of use, and computational efficiency, IRAM may be the best choice. However, their research was only conducted on the diffusion equation. In addition, they inferred that IRAM may not be suitable for parallel computing because Gram-Schmidt orthogonalization is a serial process, but no parallel performance tests were conducted on IRAM. Munoz-Cobo et al. [4] The NODAL-LAMBDA program was developed based on two groups of improved coarse-grid finite difference methods and IRAM to estimate the high-order lambda eigenvalues ​​and high-order harmonics of the reactor under different operating conditions. Numerical experiments show that the NODAL-LAMBDA program has achieved a good compromise between computational efficiency and computational accuracy, and can efficiently calculate the high-order eigenvalues ​​of the reactor under unstable events with high accuracy. Morató et al. [5] The 1-D and 2-D neutron transport equations are discretized using the discrete ordinate (SN) method and the finite difference method, and the SLEPc mathematical library is used to solve the neutron transport equations. [6] The Krylov-Schur method in the paper solves multiple eigenvalues ​​and eigenfunctions of the discretized linear equation system and develops the FORTRAN program n-DOTEC. The correctness and computational accuracy of n-DOTEC are tested using multiple benchmark problems. However, the program can only calculate 1-D and 2-D problems. Bernal et al. [7] The TRIVAC Program [8] The Krylov-Schur method in SLEPc is also used [9] Solve high-order eigenvalues ​​and eigenvectors. TRIVAC uses the Raviart-Thomas method and the Raviart-Thomas-Schneider method to discretize the neutron diffusion equation under Cartesian and Hexagonal geometry respectively. Bernal et al. compared the Krylov-Schur method with the Hotelling deflation technique and found that the Krylov-Schur method is more reliable and stable, but how to use the Krylov-Schur algorithm to effectively simulate reactor transport criticality has not given a feasible solution.

[0004] The prior art with document number CN110717275A discloses a three-dimensional neutron flux numerical simulation method for a pressurized water reactor core. The method first divides the three-dimensional pressurized water reactor core to be simulated into several layers along the axial direction, and establishes a two-dimensional neutron transport model for each layer based on the characteristic line method; then divides the three-dimensional pressurized water reactor core to be simulated into several strips based on the radial grid element geometry, and establishes a one-dimensional neutron transport model for each strip based on the discrete ordinate method; and iterates the residual model composed of the two-dimensional neutron transport model and the one-dimensional neutron transport model by the JFNK method until convergence, and obtains the neutron flux distribution of the three-dimensional pressurized water reactor core. Compared with the prior art, the invention converts the two-dimensional neutron transport model and the one-dimensional neutron transport model into residual models and solves them simultaneously, has a second-order convergence speed, and has good iterative solution stability. It can be used for transport module calculation of numerical reactors, improves the efficiency of numerical reactor transport calculation, increases calculation stability, and saves the nuclear time cost generated by numerical calculation. Summary of the invention

[0005] The technical problems to be solved by the present invention are:

[0006] The purpose of the present invention is to provide a reactor transport criticality simulation method and system based on the Krylov-Schur algorithm, so as to realize the transport criticality calculation and solve the maximum eigenvalue by using the Krylov-Schur method, and improve the convergence accuracy of the maximum eigenvalue and the eigenvector.

[0007] The technical solution adopted by the present invention to solve the above technical problems is:

[0008] A reactor transport criticality simulation method based on Krylov-Schur algorithm, the implementation process of the method is:

[0009] Step 1: Discretize the steady-state neutron transport equation and transform it into an operator form

[0010] The steady-state neutron transport equation is expressed as:

[0011] [Ω·▽+Σ t (r)]ψ(r,E,Ω)=Q(r,E,Ω) (1)

[0012] where r is a spatial variable, Ω is an angular variable, ψ(r,E,Ω) is the angular neutron flux, Σ t (r) is the total cross section, Q(r,E,Ω) is the total neutron source intensity; ▽ is the gradient operator;

[0013] By performing multi-group approximation on the energy variable E in equation (1) and discretizing the angle variable Ω with discrete ordinates, we can obtain the following coupled equations:

[0014] [Ω m ▽+Σtg (r)]ψ g,m (r) = Q g,m (r),g=1,…,G,m=1,…,M (2)

[0015] Where: m is the angle direction index, g is the energy group index, G is the maximum number of energy groups, ψ g,m (r) means: the neutron angular flux density of the mth angle of the gth group at r, Σ t,g (r) means the neutron scattering cross section (neutron total reaction cross section) of the g group at r;

[0016] Total neutron source Q g,m (r) is calculated by the following formula:

[0017]

[0018] Among them: λ is the effective multiplication coefficient, g' is the energy group index of the non-g group, Σ s,g'→g (r) is the neutron scattering cross section from the g'th to the g group at r, ν is the number of neutrons produced in each fission, Σ f,g' (r) is the neutron fission cross section of the g'th group at r, χ g is the fission spectrum of the g-th group, φg'(r) is the scalar flux of the g'-th group at r and is calculated as follows:

[0019]

[0020] where w m is the mth angle Ω m The weight of

[0021] The equation (2) can be expressed in operator form as:

[0022]

[0023] Where: L = (Ω·▽+Σ) is the input operator, D is the discrete angle to moment operator, M is the moment to discrete angle operator, S is the scattering operator, F is the fission operator; Σ is the total reaction cross section matrix; ψ is the neutron angular flux density matrix;

[0024] Step 2: Construct standard eigenvalue formula

[0025] Equation (5) can be rearranged into the following generalized eigenvalue problem:

[0026]

[0027] Where: I is the unit matrix, and the matrix A consists of G×G blocks, each of which represents an operator of energy group g:

[0028]

[0029] The symbols with a bar represent operators within the group. is the intra-group moment to discrete angle operator, is the inverse of the input operator L in the group, S gg' is the scattering operator from the g'th group to the gth group, is the discrete angle-to-moment operator within the group;

[0030] Equation (6) is further transformed into a standard eigenvalue problem:

[0031] λψ=Pψ . (8)

[0032] in

[0033] Equation (8) will serve as the starting point for the subsequent discussion of the solution to the transport criticality problem;

[0034] Step 3: Calculate the eigenvalues ​​and eigenvectors of P based on the Krylov-Schur algorithm

[0035] The Arnoldi process generates an orthogonal basis for the Krylov subspace and generates an upper Hessenberg matrix H, which represents the effect of P on this subspace; the Arnoldi process establishes the following relationship:

[0036]

[0037] Where: V k is the basis vector of the Krylov subspace, H k is the upper Hessenberg matrix, is the kth unit vector, h k+1,k Yes H k The last entry of the kth column vector of the matrix, v k+1 is the next Arnoldi vector orthogonal to the current Krylov subspace, and the eigenvalues ​​and eigenvectors of P are calculated based on the aforementioned parameter relationship;

[0038] For H m Perform Schur decomposition to obtain an orthogonal matrix (if H m is a real number) or a unitary matrix (if H m is a complex number)Q and an upper triangular matrix S;

[0039] The Ritz values ​​of the eigenvalues ​​of approximating P can be found on the diagonal of S; if these Ritz values ​​have converged, the Krylov-Schur iteration terminates; otherwise, the Ritz pairs (eigenvalues ​​and corresponding eigenvectors) are sorted and the desired eigenvalues ​​and corresponding vectors are selected; V m and Hm It is updated accordingly and an additional Arnoldi step is performed to expand the Krylov subspace;

[0040] Schur decomposition is used to extract Ritz pairs from the Hessenberg matrix and is also used to truncate the subspace, retaining important components (i.e., eigenvalues ​​and their corresponding vectors), and then restart the iterative process to ensure efficient convergence without significantly losing information from previous iterations.

[0041] Furthermore, the effective incremental coefficient λ is the first element on the diagonal of the matrix S.

[0042] Furthermore, the discrete angle to moment operator D, the moment to discrete angle operator M, the scattering operator S, and the fission operator F are determined as follows:

[0043] The calculation method of the discrete angle to moment operator D is:

[0044]

[0045] Where, Ω0,…,Ω M are discrete angles; w0,…,w M ; and is the mth angle, (l,n)th even-order and odd-order spherical harmonics;

[0046] The moment to discrete angle operator M is calculated as:

[0047]

[0048] In the formula, the symbols have the same meaning as above;

[0049] The scattering operator S is calculated as:

[0050]

[0051] where an element in S [S] gg' represents the scattering operator from the g'th group to the gth group and is calculated as follows:

[0052]

[0053] Among them, Σ s,g'→g (r) is the neutron scattering cross section from the g'th group to the gth group at r;

[0054] The fission operator F is:

[0055]

[0056] Among them, Σ f,g'(r) is the neutron fission cross section of the g'th group at r, χ g is the fission spectrum of the gth group, and ν is the number of neutrons in each fission.

[0057] A reactor transport criticality simulation system based on Krylov-Schur algorithm, the system has a program module corresponding to the steps of the above technical solution, and executes the steps of the reactor transport criticality simulation method based on Krylov-Schur algorithm when running.

[0058] A computer-readable storage medium stores a computer program, wherein the computer program is configured to implement the steps of the reactor transport criticality simulation method based on the Krylov-Schur algorithm when called by a processor.

[0059] The present invention has the following beneficial technical effects:

[0060] In reactor transport criticality calculations, the power iteration (PI) method is usually used to solve the maximum eigenvalue and the corresponding eigenvector, but the efficiency of power iteration decreases when the dominance ratio is close to 1. This paper explores the use of the Krylov-Schur method for transport criticality calculations to solve the maximum eigenvalue, and implements it in the large-scale parallel discrete ordinate transport program Marvin, and verifies the correctness and efficiency using the Takeda benchmark model 1. Preliminary numerical results show that when calculating the maximum eigenvalue and eigenvector, the Krylov-Schur method and the PI method have comparable convergence accuracy, but the calculation time increases by up to 30%. In the future, the convergence performance of the Krylov-Schur method will be further tested in depth, and research on the calculation of high-order eigenvalues ​​and eigenvectors using this method will be carried out. The keywords of the present invention are: criticality calculation; power iteration method; IRAM method; Krylov-Schur method.

[0061] The present invention adopts the Krylov-Schur method instead of the power iteration method to calculate the maximum eigenvalue and eigenvector of the transport equation, and implements it in the large-scale SN transport program Marvin independently developed by this research group. The convergence performance of the Krylov-Schur method is preliminarily tested using the Takeda benchmark model 1, and the expected technical effect is achieved. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 This is a screenshot of the pseudocode program corresponding to the Krylov-Schur method flowchart;

[0063] Figure 2 This is a schematic diagram of the principle of Takeda benchmark model 1;

[0064] Figure 3This is a flow chart of the Krylov-Schur method for neutron transport criticality calculations. DETAILED DESCRIPTION

[0065] Combined with Figure 1-3 The implementation of the reactor transport criticality simulation method based on the Krylov-Schur algorithm of the present invention is described as follows:

[0066] The steady-state neutron transport equation can be written as:

[0067] [Ω·▽+Σ t (r)]ψ(r,E,Ω)=Q(r,E,Ω) (1)

[0068] where r is a spatial variable, Ω is an angular variable, ψ(r,Ω) is the angular neutron flux, Σ t (r) is the total cross section, Q(r,E,Ω) is the total neutron source intensity; ▽ is the gradient operator;

[0069] By performing multi-group approximation on the energy variable E in equation (1) and discretizing the angle variable Ω with discrete ordinates, we can obtain the following coupled equations:

[0070] [Ω m ▽+Σ tg (r)]ψ g,m (r) = Q g,m (r), g = 1, ..., G, m = 1, ..., M (2) where m is the angle direction index, g is the energy group index, G is the maximum number of energy groups, and the total neutron source Q g,m (r) is calculated by the following formula:

[0071]

[0072] Where λ is the effective multiplication coefficient, g' is the energy group index of the non-g group, Σ s,g'→g (r) is the neutron scattering cross section from the g'th to the g group at r, ν is the number of neutrons produced in each fission, Σ f,g' (r) is the neutron fission cross section of the g'th group at r, χ g is the fission spectrum of the g-th group, φg'(r) is the scalar flux of the g'-th group at r and is calculated as follows:

[0073]

[0074] where w m is the mth angle Ω m The weight of .

[0075] Equation (2) can be expressed in operator form as:

[0076]

[0077] Where L = (Ω·▽+Σ) is the input operator, D is the discrete angle to moment operator, M is the moment to discrete angle operator, S is the scattering operator, and F is the fission operator;

[0078] Equation (5) can be rearranged into the following generalized eigenvalue problem:

[0079]

[0080] Note that I is the identity matrix and the matrix A consists of G×G blocks, each of which represents an operator of energy group g:

[0081]

[0082] The symbols with a bar represent operators within the group. is the intra-group moment to discrete angle operator, is the inverse of the input operator L in the group, S gg' is the scattering operator from the g'th group to the gth group, is the discrete angle-to-moment operator within the group;

[0083] Equation (6) is further transformed into a standard eigenvalue problem:

[0084] λψ=Pψ . (8)

[0085] in

[0086] Equation (8) will serve as the starting point for our subsequent discussion of solutions to transport criticality problems.

[0087] The iterative process of the Krylov-Schur method is as follows Figure 3 As shown, its core is the Arnoldi process, which generates an orthogonal basis for the Krylov subspace and generates an upper Hessenberg matrix H, which represents the effect of P on this subspace. The Arnoldi process establishes the following relationship:

[0088]

[0089] Where V m is the basis vector of the Krylov subspace, H m is the upper Hessenberg matrix, is the mth unit vector, h m+1,m Yes H m The last entry of the mth column vector of the matrix, v m+1 is the next Arnoldi vector that is orthogonal to the current Krylov subspace. This relationship allows the Krylov-Schur method to efficiently compute the eigenvalues ​​and eigenvectors of P in a small and manageable subspace.

[0090] Next, for H m Perform Schur decomposition to obtain an orthogonal matrix (if H m is a real number) or a unitary matrix (if H m is a complex number) Q and an upper triangular matrix S. The Ritz values ​​(the eigenvalues ​​of approximating P) can be found on the diagonal of S. If these Ritz values ​​have converged, the Krylov-Schur iteration terminates. Otherwise, the Ritz pairs (eigenvalues ​​and corresponding eigenvectors) are sorted and the desired eigenvalues ​​and corresponding vectors are selected. V m and H m is updated accordingly, and additional Arnoldi steps are performed to expand the Krylov subspace.

[0091] Schur decomposition is a key part of the Krylov-Schur method. It not only helps to extract Ritz pairs from the Hessenberg matrix, but also provides a systematic way to truncate the subspace, retain the important components (i.e. eigenvalues ​​and their corresponding vectors), and then restart the iterative process. This ensures efficient convergence without significantly losing information from previous iterations.

[0092] from Figure 3As can be seen in the figure, the difference between the Krylov-Schur method and the traditional Arnoldi method is the presence of a restart process; the difference from the explicit restart Arnoldi is that the restart process is implicit; the difference from the traditional IRAM is that the implicit shifted QR decomposition process is not used. In the traditional implicit restart Arnoldi method, after the Hessenberg matrix is ​​generated, the IRAM algorithm selects one or several specific shift values ​​(shifts) based on the Ritz value of the previous iteration. These shift values ​​are used in the subsequent QR step to select the required eigenvalues; then, for each shift, the algorithm simulates the shift and QR decomposition by performing a series of QR operations on the Hessenberg matrix instead of directly calculating A-σI; after all the selected shifts are completed, the new Hessenberg matrix will be used to establish a reduced-size Krylov subspace that retains all information related to the required "Ritz pairs"; finally, based on this reduced Krylov subspace, the Arnoldi iteration is continued to complete the calculation. The Krylov-Schur method first uses a series of orthogonal transformations to rearrange the Schur matrix S and put the required eigenvalues ​​at the head of S; then, after selecting a set of required "Ritz" pairs, it truncates the Krylov subspace and only retains the parts corresponding to these "Ritz pairs"; finally, it continues to perform the Arnoldi process on the basis of the reduced Krylov subspace to expand the subspace.

[0093] Numerical verification of the method of the present invention is carried out:

[0094] The present invention uses Takeda benchmark model 1 for preliminary verification (such as Figure 2 Model 1 is a simplified PWR core with two groups, which is based on the Kyoto University critical device KUCA. The model includes two operating conditions: 1) control rod insertion; 2) control rod withdrawal. The model size is 50cm×50cm×50cm, and includes three areas: fuel, reflector, control rod or cavity area.

[0095] The calculation of the working condition 1 (i.e., the control rod is proposed) was carried out in this paper, and the calculation grid size was 1cm×1cm×1. The comparison methods are the power iteration method (PI) and the Krylov-Schur method. In the PI method, the eigenvalue convergence limit is 1.0E-5; in the Krylov-Schur method, the Arnoldi iteration residual convergence limit is 1.0E-3. The eigenvalue calculation results of working condition 1 are shown in Table 1, and the flux calculation results of the Krylov-Schur method are shown in Table 2. From the results in the table, it can be seen that the function of the Krylov-Schur critical calculation solver developed in this paper is correct. The eigenvalue calculation results differ from the results of the PI solver by 10pcm and the difference from the reference solution is 50pcm; the relative error of the average flux in each zone and the MC reference solution is within 1%. In terms of calculation time, the calculation time of Krylov-Schur is similar to that of PI, but it must be pointed out that PI can only give the maximum eigenvalue and the maximum eigenvector, but the Krylov-Schur solver can give multiple eigenvalues ​​at the same time.

[0096] Table 1 Takeda benchmark model 1 PI and Krylov-Schur method calculation results of the lifting condition

[0097]

[0098] Table 2 Takeda benchmark problem model 1 Krylov-Schur method flux calculation results and relative errors with MC benchmark solutions

[0099]

[0100]

[0101] The article also calculated the working condition 2 (i.e. control rod insertion). At this time, the convergence limit of the Arnoldi iteration residual was 1.0E-4, and the other parameters were consistent with working condition 1. The eigenvalue calculation results of working condition 2 are shown in Table 3, and the flux calculation results of the Krylov-Schur method are shown in Table 4. It can also be seen from the results in the table that the Krylov-Schur solver developed in this paper is correct, the eigenvalue results are consistent with the PI results, and the relative error between the average flux in each zone and the reference solution is less than 5%. In terms of calculation time, Krylov-Schur increases 30% compared with PI, from 33.6 seconds to 45 seconds. However, it must also be pointed out that PI can only calculate the maximum eigenvalue, but Krylov-Schur can give multiple eigenvalue results at the same time.

[0102] Table 3 Calculation results of PI and Krylov-Schur method for Takeda benchmark model 1 under the condition of rod insertion

[0103]

[0104] Table 4. Krylov-Schur method flux calculation results for Takeda benchmark model 1 and relative errors with MC benchmark solutions

[0105]

[0106]

[0107] This study further investigates the difference in parallel computing performance between the discrete ordinate transport critical calculation based on the Krylov-Schur method and the traditional PI method. The model benchmark uses the Takeda benchmark model 1. At this point, in order to fully study the parallel performance of the two methods, this study increases the number of discrete grids to 200×200×200, that is, 8 million grids. Different numbers of processors (CPUs) were used for testing, namely 4, 16, 32, 64, 128, 256, 512 and 1024. The evaluation standard used is the internationally accepted standard "Parallel Computing Efficiency (PCE)"∈, which is defined as:

[0108]

[0109] Among them, S is the speedup ratio, P is the number of CPUs used in the current test, and P0 is the number of CPUs used in the initial test.

[0110] The speedup ratio S is defined as:

[0111]

[0112] Where T is the calculation time of the current test, and T0 is the calculation time of the initial test.

[0113] The test results are as follows:

[0114] Table 5 Comparison of PCE results of Krylov-Schur and PI discrete ordinate transport critical calculations

[0115]

[0116]

[0117] As can be seen from the table, compared with the PI discrete ordinate transport critical calculation, the Krylov-Schur discrete ordinate transport critical calculation method shows better parallel computing performance. From 16 cores to 1024 cores, the PCE of the Krylov discrete ordinate transport critical calculation method is always higher than that of the PI discrete ordinate transport critical calculation method, proving that this method has very high parallel scalability and has better application potential in future high-performance large-scale parallel computing.

[0118] Summary: This paper studies the application of the Krylov-Schur method in discrete ordinate transport critical calculations, and uses Takeda benchmark model 1 to conduct a preliminary study on the performance of this method. Numerical results show that for the maximum eigenvalue and eigenvector, the convergence accuracy of the Krylov-Schur method is comparable to that of the PI method, and the calculation time is the same as or 30% higher than that of the PI method; but the Krylov-Schur method can obtain multiple eigenvalues ​​and eigenvectors at the same time. In addition, this study also used Takeda benchmark model 1 to test the parallel performance of the method, and found that the method has better parallel performance than the PI method and has great potential for large-scale parallel application. In future research work, we will further study the efficiency of the Krylov-Schur method, and explore the use of the Krylov-Schur method to obtain high-order eigenvalues ​​and eigenvectors for SN transport critical calculations and their application in nuclear reactor physics calculations.

[0119] In summary, the method proposed in the present invention solves the practical technical problem of the present invention. The method of the present invention has been verified through simulation experiments and practical applications, and the technical effects and practicality claimed by the present invention have been verified.

[0120] The algorithm (method) proposed in the present invention is the underlying technical core of the present invention, and various products can be derived based on the algorithm.

[0121] Based on the algorithm (method) proposed in the present invention, a reactor transport criticality simulation system based on the Krylov-Schur algorithm is developed using a programming language. The system has program modules corresponding to the steps of the above-mentioned technical solution, and executes the steps in the reactor transport criticality simulation method based on the Krylov-Schur algorithm during operation.

[0122] The computer program of the developed system (software) is stored on a computer-readable storage medium, and the computer program is configured to implement the steps of the reactor transport criticality simulation method based on the Krylov-Schur algorithm when called by a processor, that is, the present invention is materialized on a carrier to become a computer program product.

[0123] Various implementations of the systems and techniques described herein can be realized in digital electronic circuit systems, integrated circuit systems, dedicated ASICs (application specific integrated circuits), computer hardware, firmware, software, and / or combinations thereof. These various implementations can include: being implemented in one or more computer programs that can be executed and / or interpreted on a programmable system including at least one programmable processor, which can be a special purpose or general purpose programmable processor that can receive data and instructions from a storage system, at least one input device, and at least one output device, and transmit data and instructions to the storage system, the at least one input device, and the at least one output device.

[0124] The computer programs (also referred to as programs, software, software applications, or codes) of the present invention include machine instructions for programmable processors, and these computer programs can be implemented using high-level procedural and / or object-oriented programming languages, and / or assembly / machine languages. As used herein, the terms "machine-readable medium" and "computer-readable medium" refer to any computer program product, device, and / or device (e.g., disk, optical disk, memory, programmable logic device PLD) for providing machine instructions and / or data to a programmable processor, including a machine-readable medium that receives machine instructions as machine-readable signals. The term "machine-readable signal" refers to any signal for providing machine instructions and / or data to a programmable processor.

[0125] It should be understood that the various forms of processes shown above can be used to reorder, add or delete steps. For example, the steps recorded in this application can be executed in parallel, sequentially or in different orders, as long as the expected results of the technical solution disclosed in this application can be achieved, they are all within the scope of protection of the present invention.

[0126] The references cited in the present invention are listed as follows:

[0127] [1].G.Verdu,R.Miro,D.Ginestar,V.Vidal.The implicit restarted Arnoldi method,an efficient alternative to solve the neutron diffusion equation.Ann.Nucl.Energy,1999,26:579-593.

[0128] [2].J.S.Warsa,T.A.Wareing,J.E.Morel J.M.McGhee,R.B.Lehoucq.Krylov subspaceiterations for deterministic k-Eigenvalue calculations.Nucl.Sc.Eng.2004,147:26-42.

[0129] [3].Nicolo Abrate,Giovanni Bruna,Sandra Dulla,Piero Ravetto.Assessmentof numerical methods for the evaluation of higher-order harmonics indiffusion theory.Annals of Nuclear Energy,2019,128:455-470.

[0130] [4].J.L.Munoz-Cobo,R.Miró,Aaron Wysocki,A.Soler.3D calculation of thelambda eigenvalues and eigenmodes of the two-group neutron diffusion equationby coarse-mesh nodal methods.Progress in Nuclear Energy.2019,110:393-409.

[0131] [5].S.Morato,A.Bernal,R.Miro,Jose E.Roman,G.Verdu.Calculation ofmodes of the multi-group neutron transport equation using the discreteordinates and Finite Difference Method.Annals of Nuclear Energy,2020,137:107077.

[0132] [6].Hernandez V,Roman JE,Vidal V.SLEPc:a scalable and flexibletoolkit for the solution of eigenvalue problems.ACM Trans Math Software.2005;31:351–362.

[0133] [7].A.Bernal,A.Hebert,J.E.Roman,R.Miro,G.Verdu.A Krylov-Schursolution of the eigenvalue problem for the neutron diffusion equationdiscretized with the Raviart-Thomas method.Journal of Nuclear Science andTechnology,2017,54(10):1085-1094.

[0134] [8].Hébert A.Auser guide for TRIVAC version 4.Montréal(Canada):Institut de génie nucleaire, Polytechnique de Montréal;2014.(Report no.;IGE-2293).

[0135] [9].Stewart GW.AKrylov-Schur algorithm for large eigenproblems.SIAM JMatrix Anal Appl.2002;23:601–614.

[0136]

[10] .G.Zhang,Z.Li,Marvin:Aparallel three-dimensional transport codebased on the discrete ordinates method for reactor shielding calculations.Progress inNuclear Energy,2021,137:103786.

Claims

1. A reactor transport criticality simulation method based on Krylov-Schur algorithm, characterized in that: The implementation process of the method is: Step 1: Discretize the steady-state neutron transport equation and transform it into an operator form The steady-state neutron transport equation is expressed as: where r is a spatial variable, Ω is an angular variable, ψ(r,E,Ω) is the angular neutron flux, Σ t (r) is the total cross section, Q(r,E,Ω) is the total neutron source intensity; is the gradient operator; By performing multi-group approximation on the energy variable E in equation (1) and discretizing the angle variable Ω with discrete ordinates, we can obtain the following coupled equations: Where: m is the angle direction index, g is the energy group index, G is the maximum number of energy groups, ψ g,m (r) means: the neutron angular flux density of the mth angle of the gth group at r, Σ t,g (r) means the neutron scattering cross section (neutron total reaction cross section) of the g group at r; Total neutron source Q g,m (r) is calculated by the following formula: Among them: λ is the effective multiplication coefficient, g' is the energy group index of the non-g group, Σ s,g'→g (r) is the neutron scattering cross section from the g'th to the g group at r, ν is the number of neutrons produced in each fission, Σ f,g' (r) is the neutron fission cross section of the g'th group at r, χ g is the fission spectrum of the g-th group, φg'(r) is the scalar flux of the g'-th group at r and is calculated as follows: where w m is the mth angle Ω m The weight of Equation (2) can be expressed in operator form as: in: is the input operator, D is the discrete angle to moment operator, M is the moment to discrete angle operator, S is the scattering operator, and F is the fission operator; Σ is the total reaction cross section matrix; ψ is the neutron angular flux density matrix; Step 2: Construct standard eigenvalue formula Equation (5) can be rearranged into the following generalized eigenvalue problem: Where: I is the unit matrix, and the matrix A consists of G×G blocks, each of which represents an operator of energy group g: The symbols with a bar represent operators within the group. is the intra-group moment to discrete angle operator, is the inverse of the input operator L in the group, S gg' is the scattering operator from the g'th group to the gth group, is the discrete angle-to-moment operator within the group; Equation (6) is further transformed into a standard eigenvalue problem: λψ=Pψ. (8) in Equation (8) will serve as the starting point for the subsequent discussion of the solution to the transport criticality problem; Step 3: Calculate the eigenvalues ​​and eigenvectors of P based on the Krylov-Schur algorithm The Arnoldi process generates an orthogonal basis for the Krylov subspace and generates an upper Hessenberg matrix H, which represents the effect of P on this subspace; the Arnoldi process establishes the following relationship: Where: V k is the basis vector of the Krylov subspace, H k is the upper Hessenberg matrix, is the kth unit vector, h k+1,k Yes H k The last entry of the kth column vector of the matrix, v k+1 is the next Arnoldi vector orthogonal to the current Krylov subspace, and the eigenvalues ​​and eigenvectors of P are calculated based on the aforementioned parameter relationship; For H m Perform Schur decomposition to obtain an orthogonal matrix (if H m is a real number) or a unitary matrix (if H m is a complex number)Q and an upper triangular matrix S; The Ritz values ​​of the eigenvalues ​​of approximate P can be found on the diagonal of S; if these Ritz values ​​have converged, the Krylov-Schur iteration will terminate; otherwise, the Ritz eigenvalues ​​and corresponding eigenvectors are sorted and the required eigenvalues ​​and corresponding vectors are selected; V m and H m It is updated accordingly and an additional Arnoldi step is performed to expand the Krylov subspace; Schur decomposition is used to extract Ritz pairs from the Hessenberg matrix and also to truncate the subspace, retaining the important components, i.e., eigenvalues ​​and their corresponding vectors, and then restart the iterative process to ensure efficient convergence without significantly losing information from previous iterations.

2. The reactor transport criticality simulation method based on the Krylov-Schur algorithm according to claim 1, characterized in that: The effective multiplication coefficient λ is the first element on the diagonal of the matrix S.

3. A reactor transport criticality simulation method based on Krylov-Schur algorithm according to claim 1 or 2, characterized in that: The discrete angle to moment operator D, the moment to discrete angle operator M, the scattering operator S, and the fission operator F are determined as follows: The calculation method of the discrete angle to moment operator D is: Where, Ω0,…,Ω M are discrete angles; w0,…,w M ; and is the mth angle, (l,n)th even-order and odd-order spherical harmonics; The moment to discrete angle operator M is calculated as: The scattering operator S is calculated as: where an element in S [S] gg' represents the scattering operator from the g'th group to the gth group and is calculated as follows: Among them, Σ s,g'→g (r) is the neutron scattering cross section from the g'th group to the gth group at r; The fission operator F is: Among them, Σ f,g' (r) is the neutron fission cross section of the g'th group at r, χ g is the fission spectrum of the gth group, and ν is the number of neutrons in each fission.

4. A reactor transport criticality simulation system based on Krylov-Schur algorithm, characterized in that: The system has a program module corresponding to the steps of any one of claims 1 to 3 above, and executes the steps of the reactor transport criticality simulation method based on the Krylov-Schur algorithm when running.

5. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and the computer program is configured to implement the steps of the reactor transport criticality simulation method based on the Krylov-Schur algorithm according to any one of claims 1 to 3 when called by a processor.

Images