A two-fluid model based nuclear reactor core sub-channel analysis method

Abstract

The application discloses a nuclear reactor core sub-channel analysis method based on a two-fluid model, comprising the following steps: 1, establishing a fluid domain calculation grid according to the structure of a nuclear reactor core or an assembly; 2, establishing a control equation set and a fuel rod heat conduction model for continuous gas-liquid phases and discrete droplet phases respectively; 3, considering flow heat exchange, resistance pressure drop and turbulent flow mixing and other physical phenomena in the core, implanting corresponding constitutive models; 4, discretizing the control equation based on a second-order discrete format to obtain a discrete equation set; 5, dividing parallel tasks and parallel regions through an MPI parallel strategy to form a large coefficient matrix of a pressure equation; and 6, solving the discrete equation set and the large coefficient matrix through a pressure coupling semi-implicit algorithm and a stable double-conjugate gradient method. Through the method, key physical parameter distribution of a pressurized water reactor can be more accurately and reasonably calculated, thereby providing technical guidance for reactor safety analysis and engineering design.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of method invention, and particularly relates to a nuclear reactor core sub-channel analysis method based on a two-fluid model, so as to realize stable and high-precision calculation of the distribution of key parameters such as coolant flow, temperature, pressure and DNBR of a pressurized water reactor core under steady-state and transient-state conditions, and provide technical guidance for safety analysis and engineering design of the pressurized water reactor. BACKGROUND

[0002] Nuclear power is an important part of the energy structure in China, and in recent years, it has been attracting more and more attention. In the design of a nuclear reactor system, the thermal-hydraulic analysis of the core is very critical. In order to better design the reactor and improve the safety and economy of the reactor operation, the fine distribution of physical fields such as temperature, flow rate and pressure in the core needs to be obtained to guide the overall design of the reactor. The sub-channel analysis method divides the core coolant region into different sub-channels, and can maintain high calculation precision while performing rapid calculation. At present, it is still widely used in the thermal-hydraulic calculation of the reactor core.

[0003] Design / base design basis accident analysis is an important part of the thermal-hydraulic safety analysis of a nuclear reactor. However, the traditional sub-channel analysis method does not consider the influence of the droplet phase, and uses a first-order discrete format, which is difficult to capture the trend of changes in key parameters under transient conditions. It is very important to further optimize the sub-channel analysis method for transient condition calculation and promote the safety analysis and engineering construction of China's advanced pressurized water reactor. SUMMARY

[0004] In order to solve the problems existing in the prior art, the purpose of the present application is to provide a nuclear reactor core sub-channel analysis method based on a two-fluid model, which has higher calculation stability and precision, so as to better calculate the accurate distribution of key parameters of a pressurized water reactor core under various transient conditions, and thus provide technical guidance for the thermal-hydraulic safety analysis of China's advanced pressurized water reactor.

[0005] A nuclear reactor core sub-channel analysis method based on a two-fluid model, the steps are as follows:

[0006] Step 1: According to the shape of the reactor core or component, the fluid domain calculation grid for sub-channel analysis is divided based on the center of the coolant. For each sub-channel, the number of axial grids should be consistent;

[0007] Step 2: Based on the two-fluid three-flow field model, the continuous gas phase, the continuous liquid phase and the discrete droplet phase are considered respectively, and the control equation set for sub-channel analysis is established, including the mass conservation equation, the axial momentum conservation equation, the energy conservation equation and the transverse momentum conservation equation. The specific model is:

[0008]

[0009] where:

[0010] i - coolant state, representing continuous liquid phase, continuous gas phase, or discrete droplet phase;

[0011] a i - coolant volume fraction;

[0012] p i - coolant density, kg / m 3 ;

[0013] x - channel axial length, m;

[0014] U i - coolant axial flow rate, m / s;

[0015] k - control volume adjacent to the current coolant channel;

[0016] m - set of control volumes adjacent to the current coolant channel;

[0017] e k - cross flow indicator;

[0018] V i - coolant cross flow rate, m / s;

[0019] G i - coolant phase change flow per unit volume, kg / m 3 · s;

[0020] C i - coolant cross mass exchange per unit volume, kg / m 3 · s;

[0021] p - coolant pressure, Pa;

[0022] f iw - frictional resistance coefficient;

[0023] D h - hydraulic diameter, m;

[0024] K i - local resistance coefficient;

[0025] K I - phase interface resistance coefficient;

[0026] i' - another phase corresponding to i;

[0027] U i' - coolant axial flow rate in another state corresponding to i, m / s;V i'- transverse coolant velocity of another state corresponding to i, m / s; I - indicates a phase change process;

[0028] - axial coolant velocity participating in phase change, m / s;

[0029] M i - transverse momentum exchange per unit volume of coolant, kg / m 2 ·s 2 ;

[0030] g - gravitational acceleration, m / s 2 ;

[0031] S - gap width, m;

[0032] l - turbulent flow length, defined as the center distance of the coolant channel, m; K Gi - gap friction resistance coefficient;

[0033] - transverse coolant velocity participating in phase change, m / s;

[0034] H i - specific enthalpy of coolant, J / kg;

[0035] P w - wet perimeter of the coolant channel, m;

[0036] Φ mn - heating fraction corresponding to channel m;

[0037] η wv - direct wall heating fraction;

[0038] q" - wall heat flux, W / m 2 ;

[0039] q Ii - heat transferred per unit volume of coolant and phase interface, W / m 3 ;

[0040] E i - transverse momentum exchange per unit volume of coolant, W / m 3 ;

[0041] A fuel rod heat conduction model is established, considering radial, axial and circumferential heat conduction, and the specific equation is:

[0042]

[0043] In the formula:

[0044] r - radial grid length, m;

[0045] λ - thermal conductivity of fuel rod, W / m·K;

[0046] T - temperature of fuel rod, K;

[0047] - circumferential angle, deg;

[0048] q v - volumetric heat source, W / m 3 ;

[0049] c p - specific heat capacity of fuel rod, J / kg·K;

[0050] t - calculation time, s.

[0051] Step 3: considering the flow heat transfer, resistance pressure drop, turbulent mixing and phase interface heat and mass transfer physical phenomena in the core, implanting the corresponding constitutive model, closing the control equation group;

[0052] Step 4: discretizing the control equation by staggered semi-implicit difference format, using the central difference format and the second-order discrete method composed of numerical flux limiting function to discretize the convection term in the control equation, obtaining the second-order discrete format of the control equation, reducing the numerical oscillation phenomenon in transient calculation;

[0053] Step 5: dividing the whole core into several parallel regions by expanding and dividing the calculation region through the introduction of imaginary fuel rods and subchannels, forming a large coefficient matrix of pressure equation through MPI parallel strategy;

[0054] Step 6: controlling the solution process of the control equation group composed of mass conservation, momentum conservation and energy conservation by pressure coupling semi-implicit algorithm, and iteratively solving the large coefficient matrix of pressure equation by stable double conjugate gradient method, finally obtaining the accurate distribution of coolant pressure, temperature, flow and deviation saturation boiling ratio physical parameters in the reactor core.

[0055] In step 3, the flow heat transfer, resistance pressure drop, turbulent mixing and phase interface heat and mass transfer physical phenomena in the core are considered, and the corresponding constitutive model is implanted. The specific model is as follows:

[0056] The flow heat transfer model is Gorenflo relationship:

[0057]

[0058] In the formula:

[0059] q nb - nucleate boiling heat flux, W / m 2 ;

[0060] q fc— Single-phase water forced convection heat transfer region heat flux, W / m 2 ;

[0061] q pb — Pool boiling heat transfer region heat flux, W / m 2 ;

[0062] The resistance pressure drop model is the Churchill correlation:

[0063]

[0064] In the formula:

[0065] f f — Fanning friction factor;

[0066] Re — Reynolds number;

[0067] a and b — coefficients related to the hydraulic diameter and wall roughness;

[0068] The turbulent mixing model is the Sadatomi correlation:

[0069]

[0070] In the formula:

[0071] D h — Hydraulic diameter, m;

[0072] β — Turbulent mixing coefficient;

[0073] F i — Channel shape factor;

[0074] S — Gap width, mD rod — Fuel rod diameter, m.

[0075] The second-order discrete method of the second-order discrete format of the convection term is specifically:

[0076]

[0077] In the formula:

[0078] n — Current calculation time;

[0079] n+1 — Next calculation time;

[0080] j — Mass / energy equation control volume;

[0081] j+1 — Next mass / energy equation control volume;

[0082] J — Axial momentum equation control volume;

[0083] Control the axial flow rate of the control body J next calculation time, m / s;

[0084] The numerical flux limiting function of the control body J current calculation time;

[0085] Delta x J Control body length, m;

[0086] Delta t - time step, s;

[0087] Numerical flux limiting function in formula (6) Specifically, the OSPRE format is adopted:

[0088]

[0089] In step 5, the whole core is divided into several parallel regions by introducing the imaginary fuel rods and subchannels, and the calculation region is expanded and divided, and the specific steps are as follows:

[0090] Step 1: Introduce imaginary fuel rods and imaginary subchannels, and expand the entire solution domain to a rectangle;

[0091] Step 2: The number of processes set is divided into two closest integers multiplied;

[0092] Step 3: The entire solution domain is divided according to the row and column, and the number of sub-regions is obtained, if the region has real fuel rods and subchannels, the sub-region is effective, otherwise, the sub-region is invalid;

[0093] Step 4: Remove the invalid region, and number the effective sub-regions in turn.

[0094] Compared with the prior art, the present application has the following advantages and beneficial effects:

[0095] 1. The method uses a two-fluid three-flow field model, and independently considers the continuous gas phase, continuous liquid phase and discrete droplet phase in the control equation, compared with the commonly used two-fluid model, the distribution of each physical field in the core under high gas rate working condition can be calculated more accurately;

[0096] 2. The method disperses the control equation set by the second-order dispersion method composed of the central difference format and the OSPRE second-order numerical flux limiting function, compared with the traditional first-order upwind dispersion format, the numerical oscillation phenomenon in transient calculation can be reduced;

[0097] 3. The method expands the solution domain by introducing imaginary fuel rods and subchannels, and then removes the invalid sub-regions to flexibly divide the parallel calculation regions, compared with the traditional parallel division method of one component corresponding to one thread, the calculation resources can be more flexibly allocated. BRIEF DESCRIPTION OF DRAWINGS

[0098] Figure 1 Flow chart of two-fluid three-flow subchannel analysis method based on second-order discrete format of the present application.

[0099] Figure 2 Schematic diagram of the method for dividing parallel regions. DETAILED DESCRIPTION

[0100] The present application will be further described in detail below in combination with the drawings and specific embodiments.

[0101] The present application provides a full-coupling nuclear reactor core subchannel analysis method based on finite elements, and the specific method is as follows:

[0102] The method first divides the control volume of the coolant region according to the structure of the nuclear reactor core or assembly, and establishes the fluid domain calculation grid; secondly, control equations are established for continuous gas-liquid phase and discrete droplet phase, including mass conservation equation, axial momentum conservation equation, transverse momentum conservation equation and energy conservation equation, and a heat conduction model of the fuel rod is established; meanwhile, the physical phenomena such as flow heat exchange, resistance pressure drop and turbulent mixing in the core are considered, and the corresponding constitutive model is implanted; then the control equations are discretized based on the second-order discrete format, the discrete equation set is obtained, and the numerical oscillation phenomenon in the calculation is reduced; the parallel tasks and parallel regions are divided through the MPI parallel strategy, and a large coefficient matrix of the pressure equation is formed; finally, the discrete equation set and the large pressure matrix are solved through the pressure-coupled semi-implicit (SIMPLE) algorithm and the stable bi-conjugate gradient (BiCGSTAB) method. Through the method of the present application, the high gas content transient condition of the pressurized water reactor can be more accurately and reasonably calculated, and the accurate distribution of the physical parameters such as coolant pressure, temperature, flow and DNBR is obtained, thereby providing technical guidance for reactor safety analysis and engineering design.

[0103] As shown in Figure 1 The full-coupling nuclear reactor core subchannel analysis method based on finite elements includes the following six steps:

[0104] Step 1: According to the shape of the reactor core or assembly, the fluid domain calculation grid for subchannel analysis is divided based on the coolant center, and for each subchannel, the number of axial grids should be consistent;

[0105] Step 2: Based on the two-fluid three-flow field model, the continuous gas phase, continuous liquid phase and discrete droplet phase are considered respectively, and the control equation set for subchannel analysis is established, including mass conservation equation, axial momentum conservation equation, energy conservation equation and transverse momentum conservation equation, and the specific model is:

[0106]

[0107] where:

[0108] i - coolant state, indicating continuous liquid phase, continuous gas phase, or discrete droplet phase;

[0109] a i - coolant volume fraction;

[0110] p i - coolant density, kg / m 3 ;

[0111] x - channel axial length, m;

[0112] U i - coolant axial flow rate, m / s;

[0113] k - control volume adjacent to the current coolant channel;

[0114] m - set of control volumes adjacent to the current coolant channel;

[0115] e k - cross-flow indicator;

[0116] V i - coolant cross-flow rate, m / s;

[0117] G i - coolant phase change flow per unit volume, kg / m 3 · s;

[0118] C i - coolant cross-flow mass exchange per unit volume, kg / m 3 · s;

[0119] p - coolant pressure, Pa;

[0120] f iw - frictional resistance coefficient;

[0121] D h - hydraulic diameter, m;

[0122] K i - local resistance coefficient;

[0123] K I - phase interface resistance coefficient;

[0124] i' - another phase corresponding to i;

[0125] U i' - coolant axial flow rate in another state corresponding to i, m / s; V i'- the transverse flow velocity of the coolant in another state corresponding to i, m / s; I - indicates a phase change process;

[0126] - the axial flow velocity of the coolant participating in the phase change, m / s;

[0127] M i - the transverse momentum exchange per unit volume of the coolant, kg / m 2 ·s 2 ;

[0128] g - the acceleration due to gravity, m / s 2 ;

[0129] S - the gap width, m;

[0130] l - the turbulent flow length, defined as the center distance of the coolant channel, m; K Gi - the gap friction resistance coefficient;

[0131] - the transverse flow velocity of the coolant participating in the phase change, m / s;

[0132] H i - the specific enthalpy of the coolant, J / kg;

[0133] P w - the wet perimeter of the coolant channel, m;

[0134] Φ mn - the heating share corresponding to the channel m;

[0135] η wv - the wall direct heating share;

[0136] q" - the wall heat flux density, W / m 2 ;

[0137] q Ii - the heat transferred per unit volume of the coolant and the phase interface, W / m 3 ;

[0138] E i - the transverse momentum exchange per unit volume of the coolant, W / m 3 .

[0139] A fuel rod heat conduction model is established, which considers radial, axial and circumferential heat conduction, and the specific equation is:

[0140]

[0141] In the formula:

[0142] r - the radial grid length, m;

[0143] λ - thermal conductivity of fuel rod, W / m·K;

[0144] T - temperature of fuel rod, K;

[0145] - circumferential angle, deg;

[0146] q v - volumetric heat source, W / m 3 ;

[0147] c p - specific heat capacity of fuel rod, J / kg·K;

[0148] t - calculation time, s.

[0149] Step 3: Considering the physical phenomena such as flow heat exchange, resistance pressure drop, turbulent mixing and phase interface heat and mass transfer in the core, implanting the corresponding constitutive model, and closing the control equation group;

[0150] Step 4: Discretize the control equation by staggered semi-implicit difference format, for example, the specific discrete form of the mass conservation equation is:

[0151]

[0152] In the formula:

[0153] J - center position of momentum grid;

[0154] j - center position of mass / energy grid;

[0155] n+1 - next calculation time;

[0156] n - current calculation time;

[0157] * - difference format of convection term.

[0158] The second-order discrete method composed of central difference format and numerical flux limiting function is used to discretize the convection term in the control equation, to obtain the second-order discrete format of the control equation, and to reduce the numerical oscillation phenomenon in transient calculation;

[0159] Step 5: The whole core is divided into several parallel regions by introducing imaginary fuel rods and subchannels, expanding and dividing the calculation region, and forming a large coefficient matrix of pressure equation by MPI parallel strategy;

[0160] Step 6: The solution process of the control equations of mass conservation, momentum conservation and energy conservation is controlled by the pressure coupling semi-implicit (SIMPLE) algorithm, and the large coefficient matrix of the pressure equation is iteratively solved by the stable bi-conjugate gradient (BiCGSTAB) method, and finally the accurate distribution of the physical parameters such as coolant pressure, temperature, flow rate and DNBR in the reactor core is obtained.

[0161] In step 3, the flow heat exchange, resistance pressure drop, turbulent mixing and phase interface heat and mass transfer physical phenomena in the core are considered, and the corresponding constitutive models are implanted, and the specific models are as follows:

[0162] The flow heat exchange model is the Gorenflo relationship:

[0163]

[0164] In the formula:

[0165] q nb — Nucleate boiling heat flux, W / m 2 ;

[0166] q fc — Single-phase water forced convection heat transfer area heat flux, W / m 2 ;

[0167] q pb — Pool boiling heat transfer area heat flux, W / m 2 .

[0168] The resistance pressure drop model is the Churchill relationship:

[0169]

[0170] In the formula:

[0171] f f — Fanning friction resistance coefficient;

[0172] Re — Reynolds number;

[0173] a and b — Coefficients related to hydraulic diameter and wall roughness.

[0174] The turbulent mixing model is the Sadatomi relationship:

[0175]

[0176] In the formula:

[0177] D h — Hydraulic diameter, m;

[0178] β — Turbulent mixing coefficient;

[0179] F i * —Channel shape factor;

[0180] S—Gap width, in meters;

[0181] D rod — Fuel rod diameter, in meters (m).

[0182] When processing the convection term using the second-order discretization method, the specific second-order discretization scheme is as follows:

[0183]

[0184] In the formula:

[0185] n—the current calculation time;

[0186] n+1 — the next calculation time;

[0187] j—Mass / Energy Equation Control Volume;

[0188] j+1 — the next mass / energy equation control volume;

[0189] J—Control volume of the axial momentum equation;

[0190] —Axial flow velocity of control volume J at the next calculation moment, m / s;

[0191] —The numerical flux constraint function of control volume J at the current calculation time;

[0192] Δx J — Length of the control volume, in meters;

[0193] Δt — time step, seconds;

[0194] Numerical flux constraint function in formula (6) Specifically, the OSPRE format is used:

[0195]

[0196] like Figure 2 As shown, taking the cross-shaped coolant region as an example, the specific steps for dividing parallel tasks and parallel regions using the MPI strategy are as follows:

[0197] Step 1: Introduce hypothetical fuel rods and hypothetical subchannels to expand the entire solution domain into a rectangle;

[0198] Step 2: Decompose the set number of processes into the product of the two closest integers;

[0199] Step 3: divide the whole solution domain by row and column to get sub-regions, if the region contains real fuel rods and sub-channels, the sub-region is valid, otherwise, the sub-region is invalid;

[0200] Step 4: remove the invalid region, and number the valid sub-regions in sequence.

[0201] The above is a further detailed description of the present application in combination with specific preferred embodiments, and cannot be deemed to limit the specific embodiments of the present application to this. For ordinary skilled persons in the technical field to which the present application belongs, several simple deductions or substitutions can be made without departing from the concept of the present application, and all of them shall be deemed to belong to the scope of patent protection determined by the submitted claims.

Claims

1. A method for analyzing core channels of a nuclear reactor based on a two-fluid model, characterized in that: The steps are as follows: Step 1: Based on the shape of the reactor core or assembly, divide the fluid domain computational grid for subchannel analysis in a way that is based on the coolant center. The number of axial grids should be consistent for each subchannel. Step 2: Based on the two-fluid, three-field model, considering the continuous gas phase, continuous liquid phase, and discrete droplet phase, establish the governing equations for sub-channel analysis, including the mass conservation equation, axial momentum conservation equation, energy conservation equation, and transverse momentum conservation equation. The specific model is as follows: In the formula: i—coolant state, indicating continuous liquid phase, continuous gas phase, or discrete droplet phase; α i —Coolant volume fraction; ρ i —Coolant density, kg / m³ 3 ; x — Axial length of the channel, in meters; U i —Axial velocity of coolant, m / s; k — the control unit adjacent to the current coolant channel; m — the set of control bodies adjacent to the current coolant channel; ε k — Crossflow sign; V i —Transverse flow velocity of coolant, m / s; Γ i —Coolant phase change flow rate per unit volume, kg / m³ 3 ·s; C i — Lateral mass exchange of coolant per unit volume, kg / m 3 ·s; p—coolant pressure, Pa; f iw — Friction coefficient; D h —Hydraulic diameter, m; K i —Local drag coefficient; K I —Interfacial drag coefficient; i' — the other phase corresponding to i; U i' —The axial velocity of the coolant in the other state corresponding to i, in m / s; V i' —The lateral flow velocity of the coolant in the other state corresponding to i, in m / s; I — Indicates a phase transition process; —Axial velocity of the coolant participating in the phase change, m / s; M i —Transverse momentum exchange per unit volume of coolant, kg / m 2 ·s 2 ; g — acceleration due to gravity, m / s² 2 ; S—Gap width, in meters; l — Turbulent length, defined as the center distance of the coolant channels, in meters; K Gi — Gap friction resistance coefficient; — Lateral flow velocity of the coolant participating in the phase change, m / s; H i — Specific enthalpy of coolant, J / kg; P w —Coolant passage wetted perimeter, m; Φ mn —The heating share corresponding to channel m; η wv —The proportion of direct wall heating; q"——wall heat flux density, W / m 2 ; q Ii —Heat transferred per unit volume of coolant to the phase interface, W / m 3 ; E i —Transverse momentum exchange per unit volume of coolant, W / m 3 ; A heat conduction model for the fuel rods was established, considering radial, axial, and circumferential heat conduction. The specific equations are as follows: In the formula: r—radial grid length, m; λ—thermal conductivity of the fuel rod, W / m·K; T—Fuel rod temperature, K; —Circumferential angle, deg; q v —Volume heat source, W / m 3 ; c p —Specific heat capacity of fuel rods, J / kg·K; t — Calculation time, s. Step 3: Consider the physical phenomena of flow heat transfer, drag pressure drop, turbulent mixing and phase interface heat and mass transfer in the reactor core, and implant the corresponding constitutive model to close the control equation set. Step 4: Discretize the control equations using an interleaved semi-implicit difference scheme. Use a second-order discretization method consisting of a central difference scheme and a numerical flux constraint function to discretize the convection terms in the control equations, thereby obtaining a second-order discretization scheme for the control equations and reducing numerical oscillations in transient calculations. Step 5: By introducing hypothetical fuel rods and sub-channels, the entire reactor core is divided into several parallel regions by expanding and dividing the computational region. A large coefficient matrix of the pressure equation is formed through the MPI parallel strategy. Step 6: The solution process of the control equations consisting of mass conservation, momentum conservation, and energy conservation is controlled by a pressure-coupled semi-implicit algorithm. The large coefficient matrix of the pressure equations is solved iteratively by the stable double conjugate gradient method, and finally the accurate distribution of coolant pressure, temperature, flow rate, and physical parameters deviating from saturated boiling ratio in the reactor core is obtained.

2. The method for analyzing nuclear reactor core channels based on a two-fluid model according to claim 1, characterized in that: Step 3 considers the physical phenomena of flow heat transfer, drag pressure drop, turbulent mixing, and phase interface heat and mass transfer within the reactor core, and incorporates the corresponding constitutive model, as follows: The flow heat transfer model uses the Gorenflo relation: In the formula: q nb —Nuclear boiling heat flux density, W / m 2 ; q fc —Heat flux density in the forced convection heat transfer region of single-phase water, W / m 2 ; q pb —Heat flux density of the pool-type boiling heat exchange zone, W / m 2 ; The drag-pressure drop model uses the Churchill relation: In the formula: f f —Fanning coefficient of frictional resistance; Re – Reynolds number; a and b – coefficients related to hydraulic diameter and wall roughness; The turbulent mixing model is based on the Sadatomi relation: In the formula: D h —Hydraulic diameter, m; β—Turbulent mixing coefficient; F i * —Channel shape factor; S—gap width, m D rod — Fuel rod diameter, in meters (m).

3. The method for analyzing nuclear reactor core channels based on a two-fluid model according to claim 1, characterized in that: The specific second-order discretization scheme for the convection term in the second-order discretization method is as follows: In the formula: n—the current calculation time; n+1 — the next calculation time; j—Mass / Energy Equation Control Volume; j+1 — the next mass / energy equation control volume; J—Control volume of the axial momentum equation; —Axial flow velocity of control volume J at the next calculation moment, m / s; —The numerical flux constraint function of control volume J at the current calculation time; Δx J — Length of the control volume, in meters; Δt — time step, seconds; Numerical flux constraint function in formula (6) Specifically, the OSPRE format is used:

4. The method for analyzing nuclear reactor core channels based on a two-fluid model according to claim 1, characterized in that: In step 5, the entire reactor core is divided into several parallel regions by introducing hypothetical fuel rods and sub-channels to expand and partition the computational domain. The specific steps are as follows: Step 1: Introduce hypothetical fuel rods and hypothetical subchannels to expand the entire solution domain into a rectangle; Step 2: Decompose the set number of processes into the product of the two closest integers; Step 3: Divide the entire solution domain into rows and columns to obtain several sub-regions of the process. If a sub-region contains real fuel rods and sub-channels, the sub-region is valid; otherwise, the sub-region is invalid. Step 4: Remove invalid regions and number the valid sub-regions sequentially.

Images