Pin-by-pin level steady-state physical and thermal coupling analysis method for liquid metal cooled fast reactor

Abstract

The application discloses a liquid metal cooled fast reactor pin-by-pin level steady-state physical-thermal coupling analysis method, which is based on an OpenFOAM platform, three-dimensional steady-state neutron transport equations and three-dimensional CFD control equations are coupled under a sub-channel grid, fast and high-fidelity physical-thermal coupling calculation is realized, and three-dimensional physical and thermal phenomena such as neutron space scattering, cross flow and backflow are accurately analyzed. The steps include the following: firstly, OpenMC fine modeling is adopted to generate a few-group constant library; secondly, geometric modeling and grid division are carried out; then, initial physical fields and boundary conditions are set, and physical-thermal internal iteration solving is respectively carried out, neutron flux and coolant temperature field are gradually converged through physical-thermal external iteration, and finally, steady-state physical-thermal coupling characteristics are obtained. The application provides an efficient and accurate pin-by-pin level coupling analysis method for the physical-thermal characteristics of a liquid metal cooled fast reactor, and provides effective guidance for design optimization and safety analysis.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of high-fidelity advanced nuclear reactor multi-physics coupling simulation and analysis, in particular to a pin-by-pin level steady-state physical thermal coupling analysis method for liquid metal cooled fast reactor. BACKGROUND

[0002] In the operation of liquid metal cooled fast reactor, there is a strong coupling relationship between the core neutron physics and thermal hydraulic, which has a key influence on the nuclear reactor design and safe operation. In the liquid metal cooled fast reactor, there are also different characteristics from the pressurized water reactor. First, the fast neutron energy spectrum has a wide range of variation, and the good thermal performance of liquid metal makes the temperature difference of reactor inlet and outlet larger, resulting in more obvious changes in nuclear reaction cross section. Second, fast neutrons have a long free path, which leads to the possibility that local power changes in the reactor may affect the power distribution of multiple components. Third, the liquid metal cooled fast reactor adopts a special structure of pool type or hexagonal assembly, and there are three-dimensional thermal hydraulic phenomena such as thermal stratification, cross flow, backflow, and reverse flow. Due to the above reasons, it is essential to develop a high-fidelity three-dimensional pin-by-pin level physical thermal coupling analysis program to deeply understand the physical thermal interaction, thereby optimizing the reactor design and improving the safety of nuclear reactors.

[0003] In the existing pin-by-pin level physical thermal coupling analysis program for liquid metal cooled fast reactor, the neutron physics is mostly described by point reactor or diffusion equation, and the thermal hydraulic is mostly described by sub-channel method, which cannot reflect the scattering of neutrons in space, as well as the key three-dimensional phenomena such as coolant cross flow, backflow, and reverse flow. With the development of computers, future research will need to integrate more accurate neutron physics models, more complex thermal hydraulic calculations, and comprehensive three-dimensional CFD simulations to comprehensively improve the design accuracy and operational safety of fast reactors. SUMMARY

[0004] In order to overcome the problems existing in the prior art, the purpose of the present application is to provide a liquid metal cooled fast reactor pin-by-pin level steady-state physical thermal coupling analysis method, which is based on the finite volume framework OpenFOAM of a unified platform, realizes data sharing in memory for neutron physics calculation and thermal hydraulic calculation. In the present application, the three-dimensional neutron transport equation is used to describe the physical phenomenon in the neutron physics part, and the three-dimensional CFD control equation is used to describe the physical phenomenon in the thermal hydraulic part, which can reflect the physical action behavior with high fidelity. In the overall solution of the present application, the OpenMC fine modeling is used to generate the constant library and interpolation calling of the few-group constant; the same grid (i.e. sub-channel grid division) is used to simplify the calculation, and the parameters are transmitted between different modules through physical field mapping and boundary condition mapping; the Picard iteration coupling is used. The present application can provide a fast and high-fidelity pin-by-pin level coupling analysis method for analyzing the physical thermal coupling characteristics of the liquid metal cooled fast reactor, and provide suggestions and guidance for the design and safety characteristic analysis of the liquid metal cooled fast reactor.

[0005] In order to achieve the above purpose, the present application adopts the following technical scheme:

[0006] A liquid metal cooled fast reactor pin-by-pin level steady-state physical thermal coupling analysis method, comprising the following steps:

[0007] Step 1: using OpenMC to perform pin-by-pin fine modeling on the structure and material of the liquid metal cooled fast reactor, and extracting the homogenized group constant extraction region according to the actual material and structure of the reactor core, extracting the few-group constant under different global temperature and coolant density conditions through Monte Carlo calculation, including neutron fraction χ g , energy release fission cross section (κΣ f ) g , total cross section, i.e. transport correction cross section effective fission neutron yield cross section (νΣ f ) g and scattering cross section Σ s,g'→g , and processing them into a few-group constant library readable by OpenFOAM through post-processing;

[0008] Step 2: using a geometry modeling software and a grid division software to perform geometry modeling and grid division on the liquid metal cooled fast reactor, wherein the component region is divided into a sub-channel form of grid, and the structure material region is freely divided into a grid according to the requirement; the homogenized group constant extraction region is named according to the OpenMC calculation;

[0009] Step 3: setting the initial physical field and boundary condition: the initial physical field includes neutron flux field distribution and neutron flux density second moment φ2 0 field distribution, coolant temperature T0 field distribution, coolant density p 0 field distribution and coolant velocity U 0 field distribution; physical calculation outer boundary set vacuum boundary, thermal calculation inlet boundary set coolant mass flow, the rest are according to the default setting;

[0010] Step 4: Enter the physical thermal coupling outer iteration, according to the region of the few-group constant and the initial or last iteration thermal calculation in the grid as the global temperature of the coolant temperature field distribution and the coolant density field distribution, interpolation calculation, update the few-group constant;

[0011] Step 5: Enter the neutron physics part inner iteration, solve the three-dimensional steady-state neutron transport equation by using the simplified spherical harmonic function SP3 method; the specific as follows:

[0012] Step 5-1: Calculate the scattering source term S scatter and fission source term S fission ;

[0013] Step 5-2: Solve the three-dimensional steady-state neutron transport equation after approximation by using the simplified spherical harmonic function SP3 method;

[0014] Step 5-3: Calculate the effective multiplication coefficient k eff ;

[0015] Step 5-4: Determine whether the neutron flux density and the effective multiplication coefficient k eff converge; if it converges, end the neutron physics inner iteration, enter step 6; if it does not converge, return to step 5-1;

[0016] Step 6: According to the neutron flux field distribution obtained by physical calculation and the design reference power of liquid metal cooled fast reactor, the actual power field distribution of liquid metal cooled fast reactor is calculated, which is used as the heat source term of thermal hydraulic calculation for updating;

[0017] Step 7: Update the sub-channel grid information according to the coolant velocity field distribution and the coolant density field distribution obtained by the last calculation;

[0018] Step 7-1: First, distinguish the inner channel through the number of adjacent faces of the microelement, and then distinguish the edge channel and the corner channel through the calculation result value of the sum of the normal vector of all faces of the microelement, avoid the error caused by marking the channel, and improve the calculation efficiency;

[0019] Step 7-2: According to the velocity field obtained by the last thermal hydraulic calculation, update the information including flow area A, wetted perimeter P, hydraulic diameter D h , Reynolds number Re; these information will be calculated through the resistance model, liquid metal heat transfer model, sub-channel shape factor transformation model and turbulent mixing model to obtain the resistance second-order tensor K and turbulent viscosity coefficient μt , turbulent heat conduction coefficient K t , subchannel grid shape factor A

[0020] Step 8: enter the thermal-hydraulic part iteration, using the SIMPLE algorithm to solve the three-dimensional CFD control equation;

[0021] The three-dimensional CFD control equation is as follows:

[0022] Mass conservation equation:

[0023]

[0024] Momentum conservation equation:

[0025]

[0026] Energy conservation equation:

[0027]

[0028] Where: p is the coolant density, U is the coolant velocity, p is the coolant pressure, g is the gravity coefficient, h is the specific enthalpy of the coolant, Q is the heat source term;

[0029] The SIMPLE algorithm decouples the velocity field and the pressure field for iterative solution, after adding the energy conservation equation, the order of the internal iteration equation solution is, energy conservation equation, velocity prediction equation, pressure Poisson equation, velocity correction equation; When the velocity field, pressure field and temperature field converge, the thermal part iteration is ended, and the coolant temperature field distribution, coolant velocity field distribution, coolant density field distribution, coolant pressure field distribution are obtained;

[0030] Step 9: judge whether the physical-thermal coupling outer iteration converges, if the maximum value of the two iterations is less than 0.01K, the effective multiplication coefficient k eff is less than 1e-5, the physical-thermal coupling calculation is ended; If not, return to step 4;

[0031] Step 10: perform physical-thermal coupling characteristic analysis, including analyzing and determining the necessity of physical-thermal coupling research, optimizing the design of liquid metal cooled fast reactor, improving the safety of liquid metal cooled fast reactor; Analyze the fuel Doppler temperature feedback and coolant density feedback phenomena of liquid metal cooled fast reactor; How the thermal parameters affect the cross section of the core and the structural material, thereby affecting the three-dimensional neutron physics effect of the neutron behavior; How neutron physics affects the distribution of thermal parameters, and analyze the three-dimensional thermal-hydraulic phenomena of liquid metal cooled fast reactor under the physical feedback, including cross flow, reverse flow and backflow.

[0032] The step 1 is specifically as follows:

[0033] Step 1-1: Firstly, pin-by-pin fine modeling of liquid metal cooled fast reactor geometry and materials in OpenMC, including fuel, air gap, cladding, coolant, reflector, core structure support material;

[0034] Step 1-2: Divide and homogenize the low enrichment fuel region, high enrichment fuel region, axial reflector, radial reflector, control assembly region, external support structure material region, and coolant downcomer region;

[0035] Step 1-3: Set the energy group structure division and number, set the Monte Carlo calculation particle iteration number, set different global temperatures and coolant densities for Monte Carlo calculation and homogenization of few-group constant extraction;

[0036] Step 1-4: After OpenMC generates the few-group constant library, post-process it into an OpenFOAM readable few-group constant library according to different global temperatures and coolant densities.

[0037] The step 2 is specifically as follows:

[0038] Step 2-1: Use a geometry modeling software to establish a geometric model;

[0039] Step 2-2: The component region adopts a sub-channel prism mesh division method, with the fuel rod center located at the intersection of the prisms; the radial reflector, structural material, and coolant downcomer adopt an arbitrary free mesh division method;

[0040] Step 2-3: Divide the same region for the geometry and name it according to the OpenMC homogenization extracted few-group constant region, and in step 4, the group constant interpolation calculation will be performed according to the group constant region and the global temperature and coolant density in the grid.

[0041] In the step 4,

[0042] The interpolation calling formula is as follows:

[0043]

[0044] Where, the influencing factors include global temperature and coolant density, the reference point is global temperature T0 and coolant density ρ0, the temperature change point is global temperature T1 and coolant density ρ0, and the density change point is global temperature T0 and coolant density ρ1.

[0045] The step 5 is specifically as follows:

[0046] Step 5-1: The calculation formula of fission source term S fission and scattering source term S scatter is as follows:

[0047]

[0048] S scatter =∑ s,g'→g,i φ 0,g',i

[0049] where χ g,i is the fraction of neutrons of energy group g in cell i, (νΣ f ) g',i is the effective fission neutron yield cross section of energy group g' in cell i, φ 0,g',i is the neutron flux density of energy group g' in cell i, ∑ s,g'→g,i is the scattering cross section of energy group g' to energy group g in cell i;

[0050] Step 5-2: After simplification by using the simplified spherical harmonics SP3 method, the steady-state three-dimensional neutron transport equation is written as:

[0051]

[0052] where D g,i =1 / (3∑ t,g,i ), ∑ r,g,i =∑ t,g,i -∑ s,g→g,i

[0053] where D g,i is the diffusion coefficient of energy group g in cell i, ∑ t,g,i is the transport cross section of energy group g in cell i, φ 0,g,i is the neutron flux density of energy group g in cell i, φ 2,g,i is the neutron flux density second moment of energy group g in cell i, ∑ r,g,i is the removal cross section of energy group g in cell i, ∑ s,g→g,i is the scattering cross section of energy group g to energy group g in cell i;

[0054] The boundary condition uses the approximate simplified MASHARK boundary condition:

[0055]

[0056] where φ * i =φ 0,i +2φ 2,i , φ * i is the adjoint neutron flux density in cell i, φ 0,i is the neutron flux density in cell i, φ 2,i is the neutron flux density second moment in cell i, α i is the ratio of incoming neutron flow to outgoing neutron flow in the control body, when α i =0, it is a vacuum boundary condition, when αi = 1, full reflection boundary condition;

[0057] Step 5-3: Calculate the effective multiplication factor k eff , the neutron flux density ratio of the current iteration step and the last iteration step;

[0058] Step 5-4: Determine whether the neutron flux density and the effective multiplication factor k eff converge; if the neutron flux density of the two iterations is less than 1e-4 and the effective multiplication factor k eff of the two iterations is less than 1e-5, then go to Step 6; if not, return to Step 5-1;

[0059] The step 6 is specifically as follows:

[0060] Step 6-1: The actual power distribution of the liquid metal cooled fast reactor is calculated according to the following formula:

[0061]

[0062] In the formula, P i is the actual power distribution of the liquid metal cooled fast reactor, P0 is the design reference power of the liquid metal cooled fast reactor, (vΣ f ) g',i is the effective fission neutron yield cross section of energy group g in grid i, V i is the volume of grid i;

[0063] Step 6-2: The heat source term of the thermal hydraulic calculation is calculated according to the following formula:

[0064]

[0065] In the formula, q V is the heat source term, V i,fluid is the actual coolant fluid volume in grid i after removing the solid fuel.

[0066] The step 7 is specifically as follows:

[0067] Step 7-1: First, distinguish the inner channel by the number of adjacent faces of the microelement body, the number of adjacent faces of the control body of the inner channel is 5, while the number of adjacent faces of the control body of the edge channel and the corner channel is 6; then distinguish the edge channel and the corner channel by the calculation result value of the sum of the normal vectors of all faces of the microelement body, the sum of the normal vectors of the faces of the edge channel is 0, and the sum of the normal vectors of the faces of the corner channel is greater than 0, in the regular hexagonal assembly, the sum of the normal vectors of the faces of the corner channel is 0.732, which can avoid the error caused by marking the channel and improve the calculation efficiency;

[0068] Step 7-2: According to the velocity field obtained by the last thermal hydraulic calculation, update the information including the flow area A, the wetted perimeter P, and the hydraulic diameter D hReynolds number Re; these information will be used to calculate the resistance second order tensor K, turbulent viscosity coefficient μ t , turbulent thermal conductivity κ t , sub-channel grid shape factor A;

[0069] The calculation formula of flow area in different types of channels is as follows:

[0070]

[0071] A b =N1·A bI +N2·A bE +N3·A bC

[0072] In the formula, p pin is the fuel rod pitch, d pin is the fuel rod diameter, p edge is the sum of the fuel rod diameter and the edge channel width, A bI is the flow area of the inner channel, A bE is the flow area of the edge channel, A bC is the flow area of the corner channel, A b is the total flow area, N1 is the number of inner channels, N2 is the number of edge channels, and N3 is the number of corner channels;

[0073] The calculation formula of wetted perimeter in different types of channels is as follows:

[0074]

[0075] P b =N1×P bI +N2×P bE +N3×P bC

[0076] In the formula, P bI is the wetted perimeter of the inner channel, P bE is the wetted perimeter of the edge channel, P bC is the wetted perimeter of the corner channel, and P b is the total wetted perimeter;

[0077] The calculation formula of hydraulic diameter in different types of channels is as follows:

[0078]

[0079] In the formula, D hbI is the hydraulic diameter of the inner channel, D hbE is the hydraulic diameter of the edge channel, and D hbCis the hydraulic diameter of the corner channel, D hb is the average hydraulic diameter;

[0080] The calculation formula of the Reynolds number in different types of channels is as follows:

[0081]

[0082] In the formula, Re bI is the Reynolds number of the inner channel, Re bE is the Reynolds number of the side channel, Re bC is the Reynolds number of the corner channel, Re b is the average Reynolds number, and μ is the viscosity coefficient of the coolant;

[0083] The sub-channel grid topology can be realized by modifying the above characteristic quantity information in the grid, thereby ensuring the correctness of the calculation of the three-dimensional CFD control equation.

[0084] The step 8 is specifically as follows:

[0085] According to the principle of the SIMPLE algorithm, the problem of solving the mass conservation equation and the momentum conservation equation is converted into solving a velocity prediction equation, a pressure Poisson equation and a velocity correction equation;

[0086] Step 8-1: solving the energy conservation equation;

[0087] Step 8-2: solving the velocity prediction equation;

[0088] Step 8-3: solving the pressure Poisson equation;

[0089] Step 8-4: solving the velocity correction equation;

[0090] Step 8-5: if the temperature is less than 1e-5K and the velocity is less than 1e-8ms -1 and the pressure is less than 1e-8Pa, convergence is reached, and step 9 is entered; if not, return to step 8-1

[0091] The step 10 is specifically as follows:

[0092] Step 10-1: determining the necessity of physical thermal coupling research by analyzing the change of the maximum value of the temperature field and the effective multiplication coefficient k eff in the outer iteration process, optimizing the design of the liquid metal cooled fast reactor and improving the safety of the liquid metal cooled fast reactor;

[0093] Step 10-2: analyzing the fuel Doppler temperature feedback and the coolant density feedback phenomena of the liquid metal cooled fast reactor by analyzing the change of the neutron flux distribution and the change of the core fission power distribution before and after coupling;

[0094] Step 10-3: Analyze how the thermal parameters affect the three-dimensional neutron physics effect by affecting the cross sections of the core and structural materials, thereby affecting the neutron behavior, by coupling the changes in the pre-and post-core fission power distribution and the few-group constant distribution;

[0095] Step 10-4: Obtain the distribution of neutron physics affecting the thermal parameters by coupling the changes in the pre-and post-core fission power distribution and the coolant temperature field distribution and flow field distribution, and analyze the inlet effect and the distribution of flow through the core in the liquid metal cooled fast reactor under the physical feedback, including three-dimensional thermal-hydraulic phenomena such as cross flow, reverse flow, and backflow.

[0096] Compared with the prior art, the present application has the following advantages:

[0097] 1. High fidelity and high calculation accuracy. Both physics and thermal parameters use physical equations that can describe three-dimensional phenomena. The neutron physics part uses a three-dimensional steady-state neutron transport equation and solves it using the SP3 method, and the thermal-hydraulic part uses a three-dimensional CFD control equation and uses the SIMPLE algorithm to convert the mass conservation equation and the momentum conservation equation into a velocity prediction equation, a pressure Poisson equation, and a velocity correction equation, which are iteratively solved with the energy conservation equation, enabling high-fidelity analysis of three-dimensional physical phenomena. Compared with the traditional neutron diffusion method, it can analyze spatial scattering, compared with the traditional sub-channel method, it can analyze backflow and other three-dimensional physical phenomena, and compared with the porous medium method, it considers the relative position of the fuel rod cladding and the coolant, which is closer to the actual physical model and improves the calculation accuracy.

[0098] 2. Fast calculation speed. In the present application, the neutron physics part uses the SP3 method to solve the three-dimensional steady-state neutron transport equation, and in OpenFOAM, the linear interpolation algorithm is used to obtain the few-group constants under the thermal feedback effect according to different regions, temperatures, densities, and energy groups, which has a faster solving speed compared with the Monte Carlo method. The thermal-hydraulic part uses a sub-channel simplified grid, which significantly reduces the number of grids compared with fine modeling, improving the calculation speed. The physical and thermal parts are iterated separately, improving the convergence speed of the Picard outer iteration. In addition, the present application uses a unified platform, OpenFOAM, for data transfer, shared memory, and avoids the problem of slow transmission speed of read-write interface programs.

[0099] 3. Based on the open source platform, the code flexibility is high, and the application range is wide. The method is completely based on the open source platform development, including the OpenMC few-group constant generation part and the OpenFOAM core pin-by-pin level physical-thermal coupling solving part, wherein the physical-thermal coupling part is developed in a modular way, the code flexibility is high, and the method is convenient for subsequent improvement. In addition, in the method, the same grid is adopted for physical and thermal coupling, the complexity and error caused by cross-scale calculation are avoided, the boundary condition does not need to be specially processed, and the application range is wide. BRIEF DESCRIPTION OF DRAWINGS

[0100] Figure 1 It is a few-group constant library generation method schematic diagram.

[0101] Figure 2 It is a sub-channel grid division method schematic diagram.

[0102] Figure 3 It is a pin-by-pin level core steady-state physical-thermal coupling method schematic diagram. DETAILED DESCRIPTION

[0103] The application will be further described in detail in combination with the drawings and specific embodiments:

[0104] As shown in the drawing, Figure 3 The liquid metal cooled fast reactor pin-by-pin level steady-state physical-thermal coupling analysis method of the application comprises the following steps:

[0105] Step 1: OpenMC is used to carry out pin-by-pin fine modeling on the structure and material of the liquid metal cooled fast reactor, and the homogenized group constant extraction region is divided according to the actual material and structure of the core, and the few-group constant is extracted by Monte Carlo calculation under different global temperature and coolant density conditions, including neutron fraction χ g , energy release fission cross section (κΣ f ) g , total cross section, i.e. transport correction cross section effective fission neutron yield cross section (νΣ f ) g and scattering cross section Σ s,g'→g , and the post-processing is made into a few-group constant library readable by OpenFOAM;

[0106] Step 2: geometric modeling software and grid division software are used to carry out geometric modeling and grid division on the liquid metal cooled fast reactor, wherein the component region is divided into sub-channel grid, and the structure material region is freely divided into grid according to the demand; the homogenized group constant extraction region is named according to the OpenMC calculation;

[0107] Step 3: Set initial physical field and boundary conditions: the initial physical field includes neutron flux Field distribution and second moment of neutron flux density Field distribution, coolant temperature T 0 Field distribution, coolant density p 0 Field distribution and coolant velocity U 0 Field distribution; physical calculation outer boundary set vacuum boundary, thermal calculation inlet boundary set coolant mass flow, the rest are set according to default;

[0108] Step 4: Enter the physical-thermal coupling outer iteration, according to the region and grid of the few-group constant and the coolant temperature field distribution and the coolant density field distribution of the last iteration thermal calculation as global temperature, carry out interpolation calculation, update the few-group constant;

[0109] Step 5: Enter the neutron physics part inner iteration, solve the three-dimensional steady-state neutron transport equation by using the simplified spherical harmonic function SP3 method; the specific steps are as follows:

[0110] Step 5-1: Calculate the scattering source term S scatter and the fission source term S fission ;

[0111] Step 5-2: Solve the three-dimensional steady-state neutron transport equation after approximation by using the simplified spherical harmonic function SP3 method;

[0112] Step 5-3: Calculate the effective multiplication coefficient k eff ;

[0113] Step 5-4: Determine whether the neutron flux density and the effective multiplication coefficient k eff converge; if they converge, end the neutron physics inner iteration and enter step 6; if they do not converge, return to step 5-1;

[0114] Step 6: According to the neutron flux field distribution obtained by physical calculation and the design reference power of the liquid metal cooled fast reactor, calculate the actual power field distribution of the liquid metal cooled fast reactor, which is used as the heat source term for updating the thermal hydraulic calculation;

[0115] Step 7: Update the sub-channel grid information according to the coolant velocity field distribution and the coolant density field distribution obtained by the last calculation;

[0116] Step 7-1: First, distinguish the inner channel by the number of adjacent faces of the microelement, and then distinguish the edge channel and the corner channel by the calculation result value of the sum of the normal vectors of all faces of the microelement, so as to avoid the error caused by marking the channel and improve the calculation efficiency;

[0117] Step 7-2: Update the information including flow area A, wetted perimeter P and hydraulic diameter D according to the velocity field obtained by the last thermal hydraulic calculationh , Reynolds number Re; this information will be calculated with the resistance model, liquid metal heat transfer model, sub-channel shape factor transformation model and turbulent mixing model to obtain the second-order resistance tensor K, turbulent viscosity coefficient μ t , turbulent thermal conductivity κ t , subchannel grid shape factor A;

[0118] Step 8: Enter the thermal hydraulic part and iterate, using the SIMPLE algorithm to solve the three-dimensional CFD control equations;

[0119] The three-dimensional CFD control equation is as follows:

[0120] The mass conservation equation:

[0121] ▽·ρU=0

[0122] Momentum conservation equation:

[0123] ▽·(ρUU)+K·U-▽·(ρμ t ▽U)-▽·(μ▽U)=-▽p+ρg

[0124] Energy conservation equation:

[0125] ▽·(ρUh)-▽·(A▽h)-▽·(ρκ t ▽h)=Q;

[0126] Where: ρ is the coolant density, U is the coolant velocity, p is the coolant pressure, g is the gravity coefficient, h is the coolant specific enthalpy, and Q is the heat source term;

[0127] The SIMPLE algorithm decouples the velocity field from the pressure field and performs iterative solutions. After adding the energy conservation equation, the order of solving the inner iterative equations is: energy conservation equation, velocity prediction equation, pressure Poisson equation, and velocity correction equation. When the velocity field, pressure field, and temperature field converge, the thermal part of the inner iteration ends and the coolant temperature field distribution, coolant velocity field distribution, coolant density field distribution, and coolant pressure field distribution are obtained.

[0128] Step 9: Determine whether the physical thermal coupling external iteration has converged. If the maximum value of the temperature field calculated after two iterations is less than 0.01K, the effective multiplication coefficient k eff If it is less than 1e-5, the physical-thermal coupling calculation is terminated; if it does not converge, return to step 4;

[0129] Step 10: Perform physical thermal coupling characteristic analysis, including analyzing and determining the necessity of physical thermal coupling research, optimizing the design of liquid metal cooled fast reactor, improving the safety of liquid metal cooled fast reactor; analyzing the fuel Doppler temperature feedback and coolant density feedback phenomena of liquid metal cooled fast reactor; the three-dimensional neutron physical effect of how thermal parameters affect the cross section of the core and structural materials, thereby affecting the behavior of neutrons; how neutron physics affects the distribution of thermal parameters, and analyzing the three-dimensional thermal-hydraulic phenomena of cross flow, counter flow, and backflow in liquid metal cooled fast reactors under physical feedback.

[0130] As shown in Figure 1 , it is a specific implementation of generating the few-group constant library in step 1. The first step is to generate a reactor-independent multi-group constant library from the evaluation and database through the NJOY processing program, including ENDF / B-VII, ENDF / B-VIII, JEFF3.3 in the official database of OpenMC. The second step is to model the core material and geometry through the Monte Carlo program OpenMC, set the energy group division for the neutron energy spectrum of the liquid metal cooled fast reactor, set the thermal parameters and Monte Carlo calculation settings, perform Monte Carlo calculation and group constant extraction; through post-processing, including region division and arrangement according to structure and material, energy group structure arrangement, full-field temperature and coolant density arrangement, form an OpenFOAM readable few-group constant library. In subsequent physical thermal coupling calculations, the three-dimensional steady-state neutron transport equation is solved in the unified platform OpenFOAM, which directly calls interpolation, avoiding the need to call the Monte Carlo program each time and improving calculation speed.

[0131] As shown in Figure 2 , it is a specific implementation of generating the sub-channel grid in step 2. The center of the fuel rod is located at the grid node, so the solution can reflect the influence of the fuel rod position compared to the free grid of the traditional core porous medium method. When updating the grid in step 7, since the grid area contains fuel solid regions (as shown in the solid parts of the inner channel, side channel, and corner channel) and coolant fluid regions (as shown in the shaded parts of the inner channel, side channel, and corner channel), it is necessary to update the grid volume, flow area, and other information. By updating the sub-channel grid information in this way, the three-dimensional CFD control equations are solved in the sub-channel grid, which can analyze three-dimensional thermal-hydraulic phenomena such as cross flow and backflow compared to traditional sub-channel methods, and achieve three-dimensional pin-by-pin level physical thermal coupling.

Claims

1. A method for steady-state physical-thermal coupling analysis of a pin-by-pin liquid metal-cooled fast reactor, characterized by: The steps include: Step 1: Use OpenMC to perform pin-by-pin fine modeling of the structure and materials of the liquid metal cooled fast reactor, and divide the homogenized group constant extraction area according to the actual core materials and structure. Under different global temperature and coolant density conditions, Monte Carlo calculations are performed to extract the minority group constants, including the neutron fraction χ g , energy release fission cross section (κΣ f ) g , total cross section is transport correction cross section Effective fission neutron yield cross section (νΣ f ) g and scattering cross section ∑ s,g'→g , through post-processing, it is made into a minority group constant library that can be read by OpenFOAM; Step 2: Use geometric modeling software and meshing software to perform geometric modeling and meshing of the liquid metal-cooled fast reactor. The component area is meshed in the form of sub-channels, and the structural material area is meshed freely as needed. The regions are named according to the homogenized group constants extracted in the OpenMC calculation. Step 3: Set up the initial physics and boundary conditions: The initial physics includes the neutron flux Field distribution and second-order moment of neutron flux density Field distribution, coolant temperature T 0 Field distribution, coolant density ρ 0 Field distribution and coolant velocity U 0 Field distribution; set the vacuum boundary for the outer boundary of the physical calculation and the coolant mass flow rate for the inlet boundary of the thermal calculation, and follow the default settings for the rest; Step 4: Enter the external iteration of physical thermal coupling, perform interpolation calculations based on the coolant temperature field distribution and coolant density field distribution as the global temperature calculated by the initial or last iterative thermal engineering calculation in the region and grid of the few-group constants, and update the few-group constants; Step 5: Enter the neutron physics part and iterate, using the simplified spherical harmonic function SP3 method to solve the three-dimensional steady-state neutron transport equation; the details are as follows: Step 5-1: Calculate the scattering source term S of the previous iteration scatter and the fission source term S fission ; Step 5-2: Solve the three-dimensional steady-state neutron transport equation approximated by the simplified spherical harmonic function SP3 method; Step 5-3: Calculate the effective multiplication coefficient k eff ; Step 5-4: Determine the neutron flux density and effective multiplication factor k eff Whether it converges; if it converges, end the neutron physics inner iteration and go to step 6; if it does not converge, return to step 5-1; Step 6: Based on the neutron flux field distribution obtained from physical calculations and the design reference power of the liquid metal cooled fast reactor, the actual power field distribution of the liquid metal cooled fast reactor is calculated and updated as the heat source term in the thermal hydraulic calculation; Step 7: Update the sub-channel grid information based on the coolant velocity field distribution and coolant density field distribution calculated last time; Step 7-1: First, distinguish the inner channel by the number of adjacent faces of the microelement, and then distinguish the edge channel and corner channel by the calculated value of the sum of the normal vectors of all faces of the microelement, so as to avoid errors caused by marking channels and improve calculation efficiency; Step 7-2: Based on the velocity field obtained from the last thermal hydraulic calculation, update the information including the flow area A, wetted perimeter P, and hydraulic diameter D h , Reynolds number Re; this information will be calculated with the resistance model, liquid metal heat transfer model, sub-channel shape factor transformation model and turbulent mixing model to obtain the second-order resistance tensor K, turbulent viscosity coefficient μ t , turbulent thermal conductivity κ t , subchannel grid shape factor A; Step 8: Enter the thermal hydraulic part and iterate, using the SIMPLE algorithm to solve the three-dimensional CFD control equations; The three-dimensional CFD control equation is as follows: The mass conservation equation: Momentum conservation equation: Energy conservation equation: Where: ρ is the coolant density, U is the coolant velocity, p is the coolant pressure, g is the gravity coefficient, h is the coolant specific enthalpy, and Q is the heat source term; The SIMPLE algorithm decouples the velocity field and the pressure field and solves them iteratively. After adding the energy conservation equation, the order of solving the inner iterative equations is: energy conservation equation, velocity prediction equation, pressure Poisson equation, and velocity correction equation. When the velocity field, pressure field and temperature field converge, the internal iteration of the thermal part ends, and the coolant temperature field distribution, coolant velocity field distribution, coolant density field distribution and coolant pressure field distribution are obtained; Step 9: Determine whether the physical thermal coupling external iteration has converged. If the maximum value of the temperature field calculated after two iterations is less than 0.01K, the effective multiplication coefficient k eff If it is less than 1e-5, the physical-thermal coupling calculation is terminated; if it does not converge, return to step 4; Step 10: Conduct physical-thermal coupling characteristic analysis, including analyzing and determining the necessity of physical-thermal coupling research, optimizing the design of liquid metal-cooled fast reactors, and improving the safety of liquid metal-cooled fast reactors; analyzing the fuel Doppler temperature feedback and coolant density feedback phenomena of liquid metal-cooled fast reactors; how thermal parameters affect the three-dimensional neutron physics effect of neutron behavior by affecting the cross-section of the core and structural materials; how neutron physics affects the distribution of thermal parameters, and analyzing the three-dimensional thermal-hydraulic phenomena of crossflow, backflow, and reflux in liquid metal-cooled fast reactors under the action of physical feedback.

2. The method for pin-by-pin steady-state physical-thermal coupling analysis of a liquid metal-cooled fast reactor according to claim 1, characterized in that: The step 1 is specifically as follows: Step 1-1: First, perform pin-by-pin detailed modeling of the liquid metal cooled fast reactor geometry and materials in OpenMC, including fuel, air gap, cladding, coolant, reflector layer, and core structure support materials; Step 1-2: Divide and homogenize the low-enrichment fuel area, high-enrichment fuel area, axial reflector layer, radial reflector layer, control component area, external support structure material area, and coolant descending section area; Steps 1-3: Set the energy group structure division and number, set the number of Monte Carlo calculation particle iterations, set different global temperatures and coolant densities for Monte Carlo calculation and uniform small group constant extraction; Steps 1-4: After OpenMC generates the minority group constant library, it is post-processed into a minority group constant library readable by OpenFOAM based on different global temperatures and coolant densities.

3. The method for pin-by-pin steady-state physical-thermal coupling analysis of a liquid metal-cooled fast reactor according to claim 1, characterized in that: The step 2 is specifically as follows: Step 2-1: Use geometric modeling software to build a geometric model; Step 2-2: The assembly area adopts the sub-channel prism grid division method, and the center of the fuel rod is located at the intersection of the prism. The radial reflector layer, structural material and coolant drop section adopt the arbitrary free grid division method. Step 2-3: Extract the region with few group constants according to OpenMC homogenization, divide the geometry into the same region and name it. In step 4, the group constant interpolation calculation will be performed based on the region of group constants and the global temperature and coolant density in the grid.

4. The method for pin-by-pin steady-state physical-thermal coupling analysis of a liquid metal-cooled fast reactor according to claim 1, characterized in that: In step 4 The formula for the interpolation call is as follows: Among them, the influencing factors include global temperature and coolant density, the reference point is the global temperature T0 and coolant density ρ0, the temperature change point is the global temperature T1 and coolant density ρ0, and the density change point is the global temperature T0 and coolant density ρ1.

5. The method for pin-by-pin steady-state physical-thermal coupling analysis of a liquid metal-cooled fast reactor according to claim 1, characterized in that: The step 5 is specifically as follows: Step 5-1: Fission source term S fission and the scattering source term S scatter The calculation formula is: S scatter =∑ s,g'→g,i f 0,g',i Where χ g,i is the neutron fraction of energy group g in grid i, (νΣ f ) g',i is the effective fission neutron yield cross section for energy group g' in grid i, φ 0,g',i is the neutron flux density of energy group g' in grid i, ∑ s,g'→g,i is the scattering cross section from energy group g' to energy group g in grid i; Step 5-2: After simplification using the simplified spherical harmonic function SP3 method, the steady-state three-dimensional neutron transport equation is written as: Among them, D g,i =1 / (3Σ t,g,i ),∑ r,g,i =∑ t,g,i -∑ s,g→g,i Where D g,i is the diffusion coefficient of energy group g in grid i, Σ t,g,i is the transport cross section of energy group g in grid i, φ 0,g,i is the neutron flux density of energy group g in grid i, φ 2,g,i is the second-order moment of the neutron flux density of energy group g in grid i, ∑ r,g,i is the removal cross section of energy group g in grid i, Σ s,g→g,i is the scattering cross section from energy group g to energy group g in grid i; The boundary conditions adopt the approximate simplified MASHARK boundary conditions: where φ * i =φ 0,i +2φ 2,i ,φ * i is the adjoint neutron flux density in grid i, φ 0,i is the neutron flux density in grid i, φ 2,i is the second-order moment of the neutron flux density in grid i, α i is the ratio of the neutron flux entering the control volume to the neutron flux leaving the control volume. When α i = 0 is the vacuum boundary condition, when α i =1, it is the total reflection boundary condition; Step 5-3: Calculate the effective multiplication coefficient k eff , is the neutron flux density ratio between the current iteration step and the previous iteration step; Step 5-4: Determine the neutron flux density and effective multiplication factor k eff Convergence; if the neutron flux density of two iterations is less than 1e-4 and the effective multiplication coefficient k of two iterations is eff If the value is less than 1e-5, go to step 6; if it does not converge, return to step 5-1.

6. The method for pin-by-pin steady-state physical-thermal coupling analysis of a liquid metal-cooled fast reactor according to claim 1, characterized in that: The step 6 is specifically as follows: Step 6-1: The actual power distribution calculation formula of the liquid metal cooled fast reactor is as follows: Where, P i is the actual power distribution of the liquid metal cooled fast reactor, P0 is the design basis power of the liquid metal cooled fast reactor, (νΣ f ) g',i is the effective fission neutron yield cross section of energy group g in grid i, V i is the volume of grid i; Step 6-2: The heat source term calculation formula for thermal hydraulic calculation is as follows: Where q V is the heat source term, V i,fluid is the actual coolant fluid volume in grid i after removing the solid fuel.

7. The method for pin-by-pin steady-state physical-thermal coupling analysis of a liquid metal-cooled fast reactor according to claim 1, characterized in that: The step 7 is specifically as follows: Step 7-1: First, distinguish the inner channels by the number of adjacent faces of the microelement. The number of adjacent faces of the inner channel control volume is 5, while the number of adjacent faces of the side and corner channels control volumes is 6. Then, distinguish the side channels and corner channels by the calculated value of the sum of the normal vectors of all the faces of the microelement. The sum of the face normal vectors of the side channel is 0, while the sum of the face normal vectors of the corner channel is greater than 0. In the regular hexagonal component, the sum of the face normal vectors of the corner channel is 0.

732. This can avoid errors caused by labeling channels and improve calculation efficiency. Step 7-2: Based on the velocity field obtained from the last thermal hydraulic calculation, update the information including the flow area A, wetted perimeter P, and hydraulic diameter D h , Reynolds number Re; this information will be calculated with the resistance model, liquid metal heat transfer model, sub-channel shape factor transformation model and turbulent mixing model to obtain the second-order resistance tensor K, turbulent viscosity coefficient μ t , turbulent thermal conductivity κ t , subchannel grid shape factor A; The calculation formula for the flow area in different types of channels is calculated based on geometric relationships as follows: A b =N1·A bI +N2·A bE +N3·A bC Where p pin is the distance between the fuel rods, d pin is the fuel rod diameter, p edge is the sum of the fuel rod diameter and the side channel width, A bI is the flow area of ​​the inner channel, A bE is the flow area of ​​the side channel, A bC is the flow area of ​​the corner channel, A b is the total flow area, N1 is the number of internal channels, N2 is the number of side channels, and N3 is the number of corner channels; The calculation formula for the wetted perimeter in different types of channels is as follows: P b =N1×P bI +N2×P bE +N3×P bC Where, P bI is the wetted perimeter of the inner channel, P bE is the wetted perimeter of the side channel, P bC is the wetted perimeter of the angular channel, P b is the total wetted perimeter; The calculation formula for hydraulic diameter in different types of channels is as follows: Where D hbI is the hydraulic diameter of the inner channel, D hbE is the hydraulic diameter of the side channel, D hbC is the hydraulic diameter of the corner channel, D hb is the mean hydraulic diameter; The calculation formulas for the Reynolds number in different types of channels are as follows: Where, Re bI is the Reynolds number of the inner channel, Re bE is the Reynolds number of the side channel, Re bC is the Reynolds number of the angular channel, Re b is the average Reynolds number, μ is the viscosity coefficient of the coolant; By modifying the above characteristic quantity information in the grid, the sub-channel grid topology can be realized, thereby ensuring the correctness of the calculation of the three-dimensional CFD control equation.

8. The method for pin-by-pin steady-state physical-thermal coupling analysis of a liquid metal-cooled fast reactor according to claim 1, characterized in that: The step 8 is specifically as follows: According to the SIMPLE algorithm principle, the problem of solving the mass conservation equation and momentum conservation equation is transformed into solving the velocity prediction equation, pressure Poisson equation and velocity correction equation; Step 8-1: Solve the energy conservation equation; Step 8-2: Solve the velocity prediction equation; Step 8-3: Solve the pressure Poisson equation; Step 8-4: Solve the velocity correction equation; Step 8-5: If the temperature between two iterations is less than 1e-5K and the speed is less than 1e-8ms -1 If the pressure is less than 1e-8Pa ​​and convergence is achieved, proceed to step 9; if not, return to step 8-1.

9. The method for pin-by-pin steady-state physical-thermal coupling analysis of a liquid metal-cooled fast reactor according to claim 1, characterized in that: In step 10 Step 10-1: The maximum value of the temperature field and the effective multiplication coefficient k during the external iteration process eff The necessity of physical-thermal coupling research can be determined by analyzing the changes in the temperature, optimizing the design of liquid metal-cooled fast reactors, and improving the safety of liquid metal-cooled fast reactors; Step 10-2: Analyze the fuel Doppler temperature feedback and coolant density feedback phenomena of the liquid metal-cooled fast reactor by observing the changes in neutron flux distribution and core fission power distribution before and after coupling; Step 10-3: Analyze the three-dimensional neutron physics effects of how thermal parameters affect the cross-section of the core and structural materials, thereby affecting neutron behavior, by analyzing the changes in the core fission power distribution and the distribution of the minority group constant before and after coupling; Step 10-4: By coupling the changes in the core fission power distribution before and after the coupling, and the changes in the coolant temperature field distribution and flow field distribution, we can determine how neutron physics affects the distribution of thermal parameters. We can also analyze the inlet effect and the distribution of flow through the core in the liquid metal-cooled fast reactor under the action of physical feedback, including the three-dimensional thermal-hydraulic phenomena of crossflow, backflow, and reflux.

Images