Liquid metal cooling fast reactor pin-by-pin level steady state physical thermal coupling analysis method

Abstract

The invention discloses a pin-by-pin-level steady-state physical-thermal coupling analysis method for a liquid metal cooling fast reactor, which is based on an OpenFOAM platform, realizes fast and high-fidelity physical-thermal coupling calculation by coupling a three-dimensional steady-state neutron transport equation and a three-dimensional CFD control equation under a sub-channel grid, and realizes high-fidelity physical-thermal coupling analysis of the three-dimensional steady-state neutron transport equation and the three-dimensional CFD control equation. And three-dimensional physical and thermal phenomena such as neutron space scattering, cross flow and backflow are accurately analyzed. The method comprises the steps that firstly, OpenMC is adopted for fine modeling, and a few-group constant library is generated; secondly, geometric modeling and grid division are carried out; and then, setting an initial physical field and boundary conditions, respectively carrying out internal iteration solution of the physical thermal engineering, gradually converging the neutron flux and the coolant temperature field through external iteration between the physical thermal engineering, and finally obtaining a steady-state physical-thermal coupling characteristic. The efficient and accurate pin-by-pin level coupling analysis method is provided for the physical-thermal characteristics of the liquid metal cooled fast reactor, and effective guidance is provided for design optimization and safety analysis.

Description

Technical Field

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

[0002] In the operation of liquid metal-cooled fast reactors, core neutron physics and thermal hydraulics have a strong coupling relationship, which has a key impact on the design and safe operation of nuclear reactors. In liquid metal-cooled fast reactors, there are also characteristics different from pressurized water reactors. First, the fast neutron energy spectrum varies over a wide range, and the good thermal performance of liquid metal makes the temperature difference between the reactor inlet and outlet larger, resulting in more obvious changes in the nuclear reaction cross section. Second, fast neutrons have a long free path, which causes local power changes in the reactor to affect the power distribution of multiple components. Third, liquid metal-cooled fast reactors use a special structure of pool-type or hexagonal components, and there are three-dimensional thermal hydraulic phenomena such as thermal stratification, crossflow, reflux, and backflow. It is precisely because of the above reasons that the development of a high-fidelity three-dimensional pin-by-pin level physical thermal coupling analysis program is crucial to a deep understanding of physical thermal interactions, thereby optimizing reactor design and improving nuclear reactor safety.

[0003] In existing liquid metal cooled fast reactor pin-by-pin level physical and thermal coupling analysis programs, neutron physics is mostly described by point pile or diffusion equations, and thermal hydraulics mostly use sub-channel methods, which cannot reflect the scattering of neutrons in space, as well as key three-dimensional phenomena such as coolant crossflow, reflux, and backflow. With the development of computers, future research will need to integrate higher-precision 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 of the invention

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

[0005] In order to achieve the above object, the present invention adopts the following technical solution:

[0006] A method for analyzing the steady-state physical and thermal coupling of a pin-by-pin level liquid metal cooled fast reactor comprises the following steps:

[0007] 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 structures. Under different global temperature and coolant density conditions, Monte Carlo calculations are used to extract the minority group constants, including the neutron fraction χ g , energy release fission cross section (κΣ f ) g , total cross section is the 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 library of minority group constants readable by OpenFOAM;

[0008] 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 divided into sub-channel grids, and the structural material area is freely divided into grids according to needs. The area is named according to the homogenized group constant extraction area in the OpenMC calculation;

[0009] Step 3: Set up initial physics and boundary conditions: The initial physics includes the neutron flux Field distribution and second-order moment of neutron flux density φ2 0 Field distribution, coolant temperature T0 Field distribution, coolant density ρ 0 Field distribution and coolant velocity U 0 Field distribution; set the vacuum boundary for the outer boundary of physical calculations and the coolant mass flow rate for the inlet boundary of thermal calculations, and follow the default settings for the rest;

[0010] 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 in the area and grid of the few group constants, and update the few group constants;

[0011] 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:

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

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

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

[0015] Step 5-4: Determine the neutron flux density and effective multiplication factor k eff Whether it converges; if it converges, the neutron physics internal iteration ends and goes to step 6; if it does not converge, it returns to step 5-1;

[0016] 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, the actual power field distribution of the liquid metal cooled fast reactor is calculated and updated as the heat source item of the thermal hydraulic calculation;

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

[0018] Step 7-1: First, the inner channel is distinguished by the number of adjacent faces of the micro-element, and then the edge channel and the corner channel are distinguished by the calculation result value of the sum of the normal vectors of all faces of the micro-element, so as to avoid errors caused by marking channels and improve calculation efficiency;

[0019] 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, 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;

[0020] Step 8: Enter the internal iteration of the thermal hydraulic part and use the SIMPLE algorithm to solve the three-dimensional CFD control equations;

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

[0022] The mass conservation equation:

[0023]

[0024] Momentum conservation equation:

[0025]

[0026] Energy conservation equation:

[0027]

[0028] 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;

[0029] The SIMPLE algorithm decouples the velocity field and the pressure field for iterative solution. After adding the energy conservation equation, the order of solving the internal 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.

[0030] Step 9: Determine whether the external iteration of physical thermal coupling converges. If the maximum value of the temperature field calculated by 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;

[0031] 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 physical effects 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.

[0032] The step 1 is specifically as follows:

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

[0034] Step 1-2: Divide and homogenize the low-enrichment fuel area, the high-enrichment fuel area, the axial reflection layer, the radial reflection layer, the control component area, the external support structure material area, and the coolant descending section area;

[0035] Step 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 homogenized small group constant extraction;

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

[0037] The step 2 is specifically as follows:

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

[0039] 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 reflection layer, structural material and coolant descending section adopt the arbitrary free grid division method;

[0040] Step 2-3: Extract the area with few group constants according to OpenMC homogenization, divide the geometry into the same area and name it. In step 4, the group constant interpolation calculation will be performed based on the area of ​​group constants and the global temperature and coolant density in the grid.

[0041] In step 4

[0042] The formula for the interpolation call is as follows:

[0043]

[0044] 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.

[0045] The step 5 is specifically as follows:

[0046] Step 5-1: Fission source term S fission and the scattering source term S scatter The calculation formula is:

[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 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;

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

[0051]

[0052] Among them, 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 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;

[0054] The boundary conditions adopt the approximate simplified MASHARK boundary conditions:

[0055]

[0056] 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.i =1, it is the total reflection boundary condition;

[0057] 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;

[0058] 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 eff If the value is less than 1e-5, go to step 6; if it does not converge, return to step 5-1;

[0059] The step 6 is specifically as follows:

[0060] Step 6-1: The actual power distribution calculation formula of the liquid metal cooled fast reactor is as follows:

[0061]

[0062] 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;

[0063] Step 6-2: The heat source term calculation formula for thermal hydraulic calculation is as follows:

[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, the inner channel is distinguished by the number of adjacent faces of the microelement. The number of adjacent faces of the inner channel control body is 5, while the number of adjacent faces of the side channel and corner channel control body is 6; then the side channel and corner channel are distinguished by the calculation result value of the sum of the normal vectors of all the faces of the microelement. The sum of the normal vectors of the side channel is 0, and the sum of the normal vectors of the corner channel is greater than 0. In the regular hexagonal component, the sum of the normal vectors of the corner channel is 0.732, which can avoid the errors caused by the marked channels and improve the calculation efficiency;

[0068] 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, 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;

[0069] The calculation formula for the flow area in different types of channels is calculated based on geometric relationships as follows:

[0070]

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

[0072] In the formula, 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;

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

[0074]

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

[0076] 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;

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

[0078]

[0079] Where D hbI is the hydraulic diameter of the inner channel, D hbE is the hydraulic diameter of the side channel, D hbCis the hydraulic diameter of the corner channel, D hb is the mean hydraulic diameter;

[0080] The calculation formula of 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 angular channel, Re b is the average Reynolds number, μ is the viscosity coefficient of the coolant;

[0083] 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.

[0084] The step 8 is specifically as follows:

[0085] 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;

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

[0087] Step 8-2: Solve the speed prediction equation;

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

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

[0090] Step 8-5: If the temperature of two iterations is less than 1e-5K and the speed is less than 1e-8ms -1 If the pressure is less than 1e-8Pa, the system will converge and go to step 9. If the pressure is not converged, the system will return to step 8-1.

[0091] The step 10 is specifically as follows:

[0092] 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 is 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;

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

[0094] Step 10-3: Analyze the three-dimensional neutron physics effect of how thermal parameters affect the cross-section of the core and structural materials, through the changes in the core fission power distribution and the distribution of the minority group constant before and after coupling;

[0095] Step 10-4: By coupling the changes in the core fission power distribution before and after, the coolant temperature field distribution changes and the flow field distribution changes, we can obtain how neutron physics affects the distribution of thermal parameters, and analyze the inlet effect and flow distribution through the core in the liquid metal cooled fast reactor under the action of physical feedback, including the three-dimensional thermal-hydraulic phenomena of cross flow, backflow and reflux.

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

[0097] 1. High fidelity and high calculation accuracy. Both physics and thermal engineering use physical equations that can describe three-dimensional phenomena. The neutron physics part is described by the three-dimensional steady-state neutron transport equation and solved by the SP3 method. The thermal hydraulic part is described by the three-dimensional CFD control equation and the SIMPLE algorithm is used to transform the mass conservation equation and momentum conservation equation into the velocity prediction equation, the pressure Poisson equation and the velocity correction equation, and it is solved iteratively with the energy conservation equation, which can analyze three-dimensional physical phenomena with high fidelity. Compared with the traditional neutron diffusion method, it can analyze spatial scattering, compared with the traditional sub-channel method, it can analyze three-dimensional physical phenomena such as backflow and reverse flow, and compared with the porous medium method, it takes into account 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. The neutron physics part of the method of the present invention adopts the SP3 method to solve the three-dimensional steady-state neutron transport equation. In OpenFOAM, the linear interpolation algorithm is used to call according to different regions, temperatures, densities, and energy groups to obtain the few group constants under the action of thermal feedback, which has a faster solution speed than the Monte Carlo method. The thermal hydraulic part adopts sub-channel simplified grid, which greatly reduces the number of grids compared with fine modeling and improves the calculation speed. The physics and thermal parts perform internal iterations respectively to improve the Picard external iteration convergence speed. In addition, the present invention uses a unified platform OpenFOAM to transfer data and share memory during physical and thermal coupling calculations to avoid the problem of slow transmission speed of the read-write interface program.

[0099] 3. Based on the open source platform, the code is highly flexible and has a wide range of applications. The method of the present invention is completely developed based on the open source platform, including the OpenMC small group constant generation part and the OpenFOAM core pin-by-pin level physical-thermal coupling solution part. Among them, the physical thermal coupling part is developed in a modular way, with high code flexibility, which is convenient for subsequent improvements to the method. In addition, in the method of the present invention, the same grid is used for both physical and thermal engineering to avoid the complexity and errors caused by cross-scale calculations, and the boundary conditions do not need to be specially processed, which has a wide range of applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 Schematic diagram of the method for generating a small group constant library.

[0101] Figure 2 Schematic diagram of sub-channel grid division method.

[0102] Figure 3 Schematic diagram of the steady-state physical-thermal coupling method for pin-by-pin core. DETAILED DESCRIPTION

[0103] The present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments:

[0104] like Figure 3 As shown, the present invention provides a liquid metal cooled fast reactor pin-by-pin level steady-state physical thermal coupling analysis method, comprising the following steps:

[0105] 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 structures. Under different global temperature and coolant density conditions, Monte Carlo calculations are used to extract the minority group constants, including the neutron fraction χ g , energy release fission cross section (κΣ f ) g , total cross section is the 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 library of minority group constants readable by OpenFOAM;

[0106] 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 divided into sub-channel grids, and the structural material area is freely divided into grids according to needs. The area is named according to the homogenized group constant extraction area in the OpenMC calculation;

[0107] Step 3: Set up 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 physical calculations and the coolant mass flow rate for the inlet boundary of thermal calculations, and follow the default settings for the rest;

[0108] 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 in the area and grid of the few group constants, and update the few group constants;

[0109] 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:

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

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

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

[0113] Step 5-4: Determine the neutron flux density and effective multiplication factor k eff Whether it converges; if it converges, the neutron physics internal iteration ends and goes to step 6; if it does not converge, it returns 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, the actual power field distribution of the liquid metal cooled fast reactor is calculated and updated as the heat source item of the thermal hydraulic calculation;

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

[0116] Step 7-1: First, the inner channel is distinguished by the number of adjacent faces of the micro-element, and then the edge channel and the corner channel are distinguished by the calculation result value of the sum of the normal vectors of all faces of the micro-element, so as to avoid errors caused by marking channels and improve calculation efficiency;

[0117] 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, hydraulic diameter Dh , 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 internal iteration of the thermal hydraulic part and use 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 and the pressure field for iterative solution. After adding the energy conservation equation, the order of solving the internal 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.

[0128] Step 9: Determine whether the external iteration of physical thermal coupling converges. If the maximum value of the temperature field calculated by 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: 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 physical effects 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.

[0130] like Figure 1 As shown, it is a specific implementation method for generating a small group constant library in step 1. The first step is to generate a multi-group constant library that is independent of the reactor from the evaluation and database through the NJOY processing program. In the official database of OpenMC, it includes ENDF / B-VII, ENDF / B-VIII, and JEFF3.3. The second step is to use the Monte Carlo program OpenMC to model the core material and geometric structure, 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 calculations and extract group constants; through post-processing, including regional division and sorting according to structure and material, energy group structure sorting, and sorting the full-field temperature and coolant density, a small group constant library that can be read by OpenFOAM is formed. In the subsequent physical thermal coupling calculation, the three-dimensional steady-state neutron transport equation will be solved in the unified platform OpenFOAM. Interpolation calls will be made directly, avoiding the need to call the Monte Carlo program each time, thereby improving the calculation speed.

[0131] like Figure 2 As shown, it is a specific implementation method of the sub-channel grid generation 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 the fuel solid area (as shown in the solid part of the inner channel, side channel and corner channel) and the coolant fluid area (as shown in the shaded part of the inner channel, side channel and corner channel), it is necessary to update the grid volume, flow cross-sectional 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. Compared with the traditional sub-channel method, it can analyze three-dimensional thermal hydraulic phenomena such as cross flow and reflux, and realize 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 level liquid metal cooled fast reactor, characterized in that: 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 structures. Under different global temperature and coolant density conditions, Monte Carlo calculations are used to extract the minority group constants, including the neutron fraction χ g , energy release fission cross section (κΣ f ) g , total cross section is the 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 library of minority group constants readable 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 divided into sub-channel grids, and the structural material area is freely divided into grids according to needs. The area is named according to the homogenized group constant extraction area in the OpenMC calculation; Step 3: Set up 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; the outer boundary of physical calculation is set as vacuum boundary, the inlet boundary of thermal calculation is set as coolant mass flow rate, and the rest are set as default; 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 in the area 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 step 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 factor k eff ; Step 5-4: Determine the neutron flux density and effective multiplication factor k eff Whether it converges; if it converges, the neutron physics internal iteration ends and goes to step 6; if it does not converge, it returns to step 5-1; 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, the actual power field distribution of the liquid metal cooled fast reactor is calculated and updated as the heat source item of the thermal hydraulic calculation; Step 7: Update the sub-channel grid information according to the coolant velocity field distribution and coolant density field distribution calculated last time; Step 7-1: First, the inner channel is distinguished by the number of adjacent faces of the micro-element, and then the edge channel and the corner channel are distinguished by the calculation result value of the sum of the normal vectors of all faces of the micro-element, 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, 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 internal iteration of the thermal hydraulic part and use 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 distribution of coolant temperature field, coolant velocity field, coolant density field and coolant pressure field are obtained; Step 9: Determine whether the external iteration of physical thermal coupling converges. If the maximum value of the temperature field calculated by 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 physical effects 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 analyzing the steady-state physical thermal coupling of a pin-by-pin level 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 fine modeling of the geometric structure and materials of the liquid metal cooled fast reactor 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, the high-enrichment fuel area, the axial reflection layer, the radial reflection layer, the control component area, the external support structure material area, and the coolant descending section area; Step 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 homogenized small group constant extraction; Step 1-4: After OpenMC generates the few-group constant library, it is post-processed into a few-group constant library readable by OpenFOAM according to different global temperatures and coolant densities.

3. The method for steady-state physical thermal coupling analysis of a liquid metal cooled fast reactor pin-by-pin according to claim 1, characterized in that: The step 2 is specifically as follows: Step 2-1: Use geometric modeling software to establish 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 reflection layer, structural material and coolant descending section adopt the arbitrary free grid division method; Step 2-3: Extract the area with few group constants according to OpenMC homogenization, divide the geometry into the same area and name it. In step 4, the group constant interpolation calculation will be performed based on the area of ​​group constants and the global temperature and coolant density in the grid.

4. The method for steady-state physical thermal coupling analysis of a pin-by-pin level 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 steady-state physical thermal coupling analysis of a pin-by-pin level 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 fraction of neutrons 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. i =1, it is the total reflection boundary condition; Step 5-3: Calculate the effective multiplication factor 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 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 steady-state physical thermal coupling analysis of a liquid metal cooled fast reactor pin-by-pin 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: 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.

7. The method for steady-state physical thermal coupling analysis of a liquid metal cooled fast reactor pin-by-pin according to claim 1, characterized in that: The step 7 is specifically as follows: Step 7-1: First, the inner channel is distinguished by the number of adjacent faces of the microelement. The number of adjacent faces of the inner channel control body is 5, while the number of adjacent faces of the side channel and corner channel control body is 6; then the side channel and corner channel are distinguished by the calculation result value of the sum of the normal vectors of all faces of the microelement. The sum of the normal vectors of the side channel is 0, and the sum of the normal vectors of the corner channel is greater than 0. In the regular hexagonal component, the sum of the normal vectors of the corner channel is 0.732, which can avoid the errors caused by the marked channels and improve the 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, 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 In the formula, 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 formula of Reynolds number in different types of channels is as follows: 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 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 steady-state physical thermal coupling analysis of a pin-by-pin level 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 speed prediction equation; Step 8-3: Solve the pressure Poisson equation; Step 8-4: Solve the velocity correction equation; Step 8-5: If the temperature of two iterations is less than 1e-5K and the speed is less than 1e-8ms -1 If the pressure is less than 1e-8Pa, convergence is achieved and the process goes to step 9. If not, the process returns to step 8-1.

9. The method for steady-state physical thermal coupling analysis of a liquid metal cooled fast reactor pin-by-pin 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 is 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 analyzing the changes in neutron flux distribution and core fission power distribution before and after coupling; Step 10-3: Analyze the three-dimensional neutron physics effect of how thermal parameters affect the cross-section of the core and structural materials, through 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 coolant temperature field distribution changes and the flow field distribution changes, we can obtain how neutron physics affects the distribution of thermal parameters, and analyze the inlet effect and flow distribution through the core in the liquid metal cooled fast reactor under the action of physical feedback, including the three-dimensional thermal-hydraulic phenomena of cross flow, backflow and reflux.

Images