Nuclear reactor core channel analysis method based on two fluid models

Abstract

The invention discloses a nuclear reactor core channel analysis method based on two fluid models, which comprises the following steps of: 1, establishing a fluid domain computational grid according to the structure of a nuclear reactor core or an assembly; 2, respectively establishing a control equation set and a heat conduction model of the fuel rod for the continuous gas-liquid phase and the discrete liquid drop phase; 3, considering physical phenomena such as flow heat exchange, resistance pressure drop and turbulence mixing in a reactor core, and implanting a corresponding constitutive model; 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 a discrete equation set and a large coefficient matrix through a pressure coupling semi-implicit algorithm and a stable biconjugate gradient method. According to the method, the key physical parameter distribution of the pressurized water reactor can be more accurately and reasonably calculated, so that technical guidance is provided for safety analysis and engineering design of the reactor.

Description

Technical Field

[0001] The present invention belongs to the technical field of method inventions, and specifically relates to a nuclear reactor core sub-channel analysis method based on a two-fluid model, so as to achieve 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 / transient conditions, and provide technical guidance for the safety analysis and engineering design of pressurized water reactors. Background Art

[0002] As an important part of my country's energy structure, nuclear power has received more and more attention in recent years. In the design of nuclear reactor systems, core thermal-hydraulic analysis is very critical. In order to better design the reactor and improve the safety and economy of reactor operation, it is necessary to obtain the fine distribution of physical fields such as temperature, flow rate and pressure in the core to guide the overall design of the reactor. The subchannel analysis method divides the core coolant area into different subchannels, which can maintain high calculation accuracy while performing fast calculations. It is still widely used in thermal-hydraulic calculations of reactor cores.

[0003] Design / beyond-design basis accident analysis is an important part of thermal-hydraulic safety analysis of nuclear reactors. However, the traditional sub-channel analysis method does not consider the influence of the droplet phase and uses a first-order discrete format, which makes it difficult to capture the changing trend of key parameters under transient conditions. It is crucial to further optimize the sub-channel analysis method for transient condition calculations and promote the safety analysis and engineering construction of my country's advanced pressurized water reactors. Summary of the invention

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

[0005] A method for analyzing a nuclear reactor core subchannel based on a two-fluid model, the steps are as follows:

[0006] Step 1: According to the shape of the reactor core or component, divide the fluid domain calculation grid for sub-channel analysis based on the coolant center. 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, continuous liquid phase and discrete droplet phase are considered respectively, and the control equations for sub-channel analysis are established, including the mass conservation equation, axial momentum conservation equation, energy conservation equation and lateral momentum conservation equation. The specific model is:

[0008]

[0009] Where:

[0010] i – coolant state, indicating continuous liquid phase, continuous gas phase or discrete droplet phase;

[0011] α i - coolant volume share;

[0012] ρ i ——Coolant density, kg / m 3 ;

[0013] x——axial length of the channel, m;

[0014] U i ——Coolant axial velocity, m / s;

[0015] k——the control body adjacent to the current coolant channel;

[0016] m——the set of control bodies adjacent to the current coolant channel;

[0017] ε k - Cross-flow identification;

[0018] V i ——Coolant lateral flow velocity, m / s;

[0019] Γ i ——Coolant phase change flow rate per unit volume, kg / m 3 ·s;

[0020] C i ——Transverse mass exchange per unit volume of coolant, kg / m 3 ·s;

[0021] p——coolant pressure, Pa;

[0022] f iw —— resistance coefficient along the way;

[0023] D h ——hydraulic diameter, m;

[0024] K i ——local resistance coefficient;

[0025] K I ——interface resistance coefficient;

[0026] i'——the other phase corresponding to i;

[0027] U i' ——Axial flow velocity of coolant 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——represents the phase change process;

[0028] ——Axial flow velocity of the coolant involved in phase change, m / s;

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

[0030] g——acceleration due to gravity, m / s 2 ;

[0031] S——gap width, m;

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

[0033] ——Transverse flow velocity of the coolant involved in phase change, m / s;

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

[0035] P w ——Wetted perimeter of coolant channel, m;

[0036] Φ mn ——Heating share corresponding to channel m;

[0037] η wv ——Direct heating share of the wall;

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

[0039] q Ii ——The heat transferred between the coolant unit volume and the phase interface, W / m 3 ;

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

[0041] The fuel rod thermal conduction model is established, taking into account radial, axial and circumferential heat conduction. The specific equation is:

[0042]

[0043] Where:

[0044] r——radial grid length, m;

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

[0046] T——fuel rod temperature, K;

[0047] ——circumferential angle, deg;

[0048] q v ——Volume 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: Consider the flow heat transfer, resistance pressure drop, turbulent mixing and phase interface heat and mass transfer physical phenomena in the core, implant the corresponding constitutive model, and close the control equation group;

[0052] Step 4: Discretize the control equations through the staggered semi-implicit difference format, and use the second-order discretization method consisting of the central difference format and the numerical flux limit function to discretize the convection terms in the control equations to obtain the second-order discretization format of the control equations, thereby reducing the numerical oscillation phenomenon in transient calculations;

[0053] Step 5: By introducing imaginary fuel rods and sub-channels, the whole core is divided into several parallel regions by expanding and dividing the calculation area, and a large coefficient matrix of the pressure equation is formed through the MPI parallel strategy;

[0054] Step 6: The solution process of the control equations consisting of conservation of mass, momentum and energy is controlled by the pressure-coupled semi-implicit algorithm. The large coefficient matrix of the pressure equation is iteratively solved by the stable biconjugate gradient method to finally obtain the accurate distribution of the coolant pressure, temperature, flow rate and deviation from the saturated boiling ratio 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 the Gorenflo relationship:

[0057]

[0058] Where:

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

[0060] q fc——Heat flux density in single-phase water forced convection heat transfer area, W / m 2 ;

[0061] q pb ——Heat flux density in pool boiling heat transfer area, W / m 2 ;

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

[0063]

[0064] Where:

[0065] f f ——Fanning friction coefficient;

[0066] Re——Reynolds number;

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

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

[0069]

[0070] Where:

[0071] D h ——hydraulic diameter, m;

[0072] β——turbulent mixing coefficient;

[0073] F i *--channel shape factor;

[0074] S——gap width, mD rod ——Diameter of fuel rod, m.

[0075] The second-order discretization format of the convection term in the second-order discretization method is as follows:

[0076]

[0077] Where:

[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——control volume of axial momentum equation;

[0083] ——The axial velocity of the control body J at the next calculation time, m / s;

[0084] ——Numerical flux limit function of control volume J at the current calculation time;

[0085] Δx J ——control body length, m;

[0086] Δt——time step, s;

[0087] The numerical flux limiting function in formula (6) is The specific format is OSPRE:

[0088]

[0089] In step 5, the entire core is divided into several parallel regions by introducing imaginary fuel rods and sub-channels, expanding and dividing the calculation area. The specific steps are as follows:

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

[0091] Step 2: Decompose the set number of processes into two nearest integers and multiply them;

[0092] Step 3: Divide the entire solution domain equally by rows and columns to obtain a number of sub-regions. If there are real fuel rods and sub-channels in the region, the sub-region is valid, otherwise, the sub-region is invalid.

[0093] Step 4: Remove the invalid area and number the valid sub-areas in sequence.

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

[0095] 1. This method uses a two-fluid three-flow field model, which independently considers the continuous gas phase, continuous liquid phase and discrete droplet phase in the control equation. Compared with the currently common two-fluid model, it can more accurately calculate the distribution of various physical fields in the core under high gas content conditions;

[0096] 2. This method discretizes the control equations through a second-order discretization method consisting of a central difference format and an OSPRE second-order numerical flux limit function. Compared with the traditional first-order upwind discretization format, it can reduce the numerical oscillation phenomenon in transient calculations.

[0097] 3. This method introduces imaginary fuel rods and sub-channels, expands the solution domain, and then removes invalid sub-regions to flexibly divide the parallel computing area. Compared with the traditional parallel division method of one component corresponding to one thread, it can allocate computing resources more flexibly. BRIEF DESCRIPTION OF THE DRAWINGS

[0098] Figure 1 The present invention is a flow chart of the two-fluid three-flow field sub-channel analysis method based on the second-order discrete format.

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

[0100] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0101] The present invention provides a finite element-based method for analyzing a fully coupled nuclear reactor core subchannel. The specific method is as follows:

[0102] The method first divides the coolant area into control bodies according to the structure of the nuclear reactor core or component, and establishes a fluid domain calculation grid; secondly, establishes control equations for the continuous gas-liquid phase and the discrete droplet phase, including the mass conservation equation, the axial momentum conservation equation, the lateral momentum conservation equation and the energy conservation equation, and establishes a thermal conductivity model of the fuel rod; at the same time, considers the physical phenomena such as flow heat transfer, resistance pressure drop and turbulent mixing in the core, and implants the corresponding constitutive model; then discretizes the control equation based on the second-order discrete format, obtains the discrete equation group and reduces the numerical oscillation phenomenon in the calculation; divides the parallel tasks and parallel regions through the MPI parallel strategy to form a large coefficient matrix of the pressure equation; finally, solves the discrete equation group and the large pressure matrix through the pressure coupling semi-implicit (SIMPLE) algorithm and the stable double conjugate gradient (BiCGSTAB) method. The method of the present invention can more accurately and reasonably calculate the transient working conditions of the high gas content of the pressurized water reactor, obtain the accurate distribution of physical parameters such as coolant pressure, temperature, flow rate and DNBR, thereby providing technical guidance for reactor safety analysis and engineering design.

[0103] like Figure 1 As shown in the figure, the finite element-based fully coupled nuclear reactor core subchannel analysis method includes the following six steps:

[0104] Step 1: According to the shape of the reactor core or component, divide the fluid domain calculation grid for sub-channel analysis based on the coolant center. For each sub-channel, 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 equations for sub-channel analysis are established, including the mass conservation equation, axial momentum conservation equation, energy conservation equation and lateral momentum conservation equation. The specific model is:

[0106]

[0107] Where:

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

[0109] α i - coolant volume share;

[0110] ρ i ——Coolant density, kg / m 3 ;

[0111] x——axial length of the channel, m;

[0112] U i ——Coolant axial velocity, m / s;

[0113] k——the control body adjacent to the current coolant channel;

[0114] m——the set of control bodies adjacent to the current coolant channel;

[0115] ε k - Cross-flow identification;

[0116] V i ——Coolant lateral flow velocity, m / s;

[0117] Γ i ——Coolant phase change flow rate per unit volume, kg / m 3 ·s;

[0118] C i ——Transverse mass exchange per unit volume of coolant, kg / m 3 ·s;

[0119] p——coolant pressure, Pa;

[0120] f iw —— resistance coefficient along the way;

[0121] D h ——hydraulic diameter, m;

[0122] K i ——local resistance coefficient;

[0123] K I ——interface resistance coefficient;

[0124] i'——the other phase corresponding to i;

[0125] U i' ——Axial flow velocity of coolant 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——represents the phase change process;

[0126] ——Axial flow velocity of the coolant involved in phase change, m / s;

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

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

[0129] S——gap width, m;

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

[0131] ——Transverse flow velocity of the coolant involved in phase change, m / s;

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

[0133] P w ——Wetted perimeter of coolant channel, m;

[0134] Φ mn ——Heating share corresponding to channel m;

[0135] η wv ——Direct heating share of the wall;

[0136] q"——wall heat flux, W / m 2 ;

[0137] q Ii ——The heat transferred between the coolant unit volume and the phase interface, W / m 3 ;

[0138] E i ——Transverse momentum exchange per unit volume of coolant, W / m 3 .

[0139] The fuel rod thermal conduction model is established, taking into account radial, axial and circumferential heat conduction. The specific equation is:

[0140]

[0141] Where:

[0142] r——radial grid length, m;

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

[0144] T——fuel rod temperature, K;

[0145] ——circumferential angle, deg;

[0146] q v ——Volume 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: Consider the physical phenomena such as flow heat transfer, resistance pressure drop, turbulent mixing and heat and mass transfer at the phase interface in the core, implant the corresponding constitutive model, and close the control equation group;

[0150] Step 4: Discretize the control equations using the staggered semi-implicit difference format. Taking the mass conservation equation as an example, the specific discretization form is:

[0151]

[0152] Where:

[0153] J——center position of momentum grid;

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

[0155] n+1——next calculation time;

[0156] n——current calculation time;

[0157] *——Difference format of convection terms.

[0158] The convection term in the control equation is discretized by using the second-order discretization method consisting of the central difference format and the numerical flux limit function, and the second-order discretization format of the control equation is obtained to reduce the numerical oscillation phenomenon in the transient calculation.

[0159] Step 5: By introducing imaginary fuel rods and sub-channels, the whole core is divided into several parallel regions by expanding and dividing the calculation area, and a large coefficient matrix of the pressure equation is formed through the MPI parallel strategy;

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

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

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

[0163]

[0164] Where:

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

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

[0167] q pb ——Heat flux density in pool boiling heat transfer area, W / m 2 .

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

[0169]

[0170] Where:

[0171] f f ——Fanning friction coefficient;

[0172] Re——Reynolds number;

[0173] a and b are coefficients related to hydraulic diameter and wall roughness.

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

[0175]

[0176] Where:

[0177] D h ——hydraulic diameter, m;

[0178] β——turbulent mixing coefficient;

[0179] F i * — channel shape factor;

[0180] S——gap width, m;

[0181] D rod ——Diameter of fuel rod, m.

[0182] When the convection term is processed by the second-order discretization method, its second-order discretization format is as follows:

[0183]

[0184] Where:

[0185] n——current calculation time;

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

[0187] j——mass / energy equation control volume;

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

[0189] J——control volume of axial momentum equation;

[0190] ——The axial velocity of the control body J at the next calculation time, m / s;

[0191] ——Numerical flux limit function of control volume J at the current calculation time;

[0192] Δx J ——control body length, m;

[0193] Δt——time step, s;

[0194] The numerical flux limiting function in formula (6) is The specific format is OSPRE:

[0195]

[0196] like Figure 2 As shown in the figure, taking the "cross" shaped coolant area as an example, the steps of dividing the parallel tasks and parallel areas through the MPI strategy are as follows:

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

[0198] Step 2: Decompose the set number of processes into two nearest integers and multiply them;

[0199] Step 3: Divide the entire solution domain equally by rows and columns to obtain a number of sub-regions. If there are real fuel rods and sub-channels in the region, the sub-region is valid, otherwise, the sub-region is invalid.

[0200] Step 4: Remove the invalid area and number the valid sub-areas in sequence.

[0201] The above content is a further detailed description of the present invention in combination with a specific preferred embodiment. It cannot be determined that the specific embodiments of the present invention are limited to this. For ordinary technicians in the technical field to which the present invention belongs, several simple deductions or substitutions can be made without departing from the concept of the present invention, which should be regarded as belonging to the present invention and the scope of patent protection determined by the submitted claims.

Claims

1. A nuclear reactor core subchannel analysis method based on a two-fluid model, characterized in that: Here are the steps: Step 1: According to the shape of the reactor core or component, divide the fluid domain calculation grid for sub-channel analysis based on the coolant center. For each sub-channel, the number of axial grids should be consistent; 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 equations for sub-channel analysis are established, including the mass conservation equation, axial momentum conservation equation, energy conservation equation and lateral momentum conservation equation. The specific model is: Where: i – coolant state, indicating continuous liquid phase, continuous gas phase or discrete droplet phase; α i - coolant volume share; ρ i ——Coolant density, kg / m 3 ; x——axial length of the channel, m; U i ——Axial flow velocity of coolant, m / s; k——the control body adjacent to the current coolant channel; m——the set of control bodies adjacent to the current coolant channel; ε k - Cross-flow identification; V i ——Coolant lateral flow velocity, m / s; Γ i ——Coolant phase change flow rate per unit volume, kg / m 3 ·s; C i ——Transverse mass exchange per unit volume of coolant, kg / m 3 ·s; p——coolant pressure, Pa; f iw —— resistance coefficient along the way; D h ——hydraulic diameter, m; K i ——local resistance coefficient; K I ——interface resistance coefficient; i'——the other phase corresponding to i; U i' ——Axial flow velocity of coolant 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 the phase change process; ——Axial flow velocity of the coolant involved in 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, m; l – turbulence length, defined as the center distance of the coolant channel, m; K Gi ——gap friction resistance coefficient; ——Transverse flow velocity of the coolant involved in phase change, m / s; H i ——specific enthalpy of coolant, J / kg; P w ——Wetted perimeter of coolant channel, m; Φ mn ——Heating share corresponding to channel m; η wv ——Direct heating share of the wall; q"——wall heat flux, W / m 2 ; q Ii ——The heat transferred between the coolant unit volume and the phase interface, W / m 3 ; E i ——Transverse momentum exchange per unit volume of coolant, W / m 3 ; The fuel rod thermal conduction model is established, taking into account radial, axial and circumferential heat conduction. The specific equation is: Where: r——radial grid length, m; λ——thermal conductivity of 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 rod, J / kg·K; t——calculation time, s. Step 3: Consider the flow heat transfer, resistance pressure drop, turbulent mixing and phase interface heat and mass transfer physical phenomena in the core, implant the corresponding constitutive model, and close the control equation group; Step 4: Discretize the control equations through the staggered semi-implicit difference format, and use the second-order discretization method consisting of the central difference format and the numerical flux limit function to discretize the convection terms in the control equations to obtain the second-order discretization format of the control equations, thereby reducing the numerical oscillation phenomenon in transient calculations; Step 5: By introducing imaginary fuel rods and sub-channels, the whole core is divided into several parallel regions by expanding and dividing the calculation area, and 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 conservation of mass, momentum and energy is controlled by the pressure-coupled semi-implicit algorithm. The large coefficient matrix of the pressure equation is iteratively solved by the stable biconjugate gradient method to finally obtain the accurate distribution of the coolant pressure, temperature, flow rate and deviation from the saturated boiling ratio in the reactor core.

2. The method for analyzing a nuclear reactor core subchannel based on a two-fluid model according to claim 1, characterized in that: 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: The flow heat transfer model is the Gorenflo relationship: Where: q nb ——Nucleate boiling heat flux, W / m 2 ; q fc ——Heat flux density in single-phase water forced convection heat transfer area, W / m 2 ; q pb ——Heat flux density in pool boiling heat transfer area, W / m 2 ; The resistance pressure drop model is the Churchill relationship: Where: f f ——Fanning friction coefficient; Re——Reynolds number; a and b – coefficients related to hydraulic diameter and wall roughness; The turbulent mixing model is the Sadatomi relationship: Where: D h ——hydraulic diameter, m; β——turbulent mixing coefficient; F i * — channel shape factor; S——gap width, m D rod ——Diameter of fuel rod, m.

3. The method for analyzing a nuclear reactor core subchannel based on a two-fluid model according to claim 1, characterized in that: The second-order discretization format of the convection term in the second-order discretization method is as follows: Where: n——current calculation time; n+1——next calculation time; j——mass / energy equation control volume; j+1——next mass / energy equation control volume; J——control volume of axial momentum equation; ——The axial velocity of the control body J at the next calculation time, m / s; ——Numerical flux limit function of control volume J at the current calculation time; Δx J ——control body length, m; Δt——time step, s; The numerical flux limiting function in formula (6) is The specific format is OSPRE:

4. The method for analyzing a nuclear reactor core subchannel based on a two-fluid model according to claim 1, characterized in that: In step 5, the entire core is divided into several parallel regions by introducing imaginary fuel rods and sub-channels, expanding and dividing the calculation area. The specific steps are as follows: Step 1: Introduce imaginary fuel rods and imaginary subchannels, and expand the entire solution domain into a rectangle; Step 2: Decompose the set number of processes into two nearest integers and multiply them; Step 3: Divide the entire solution domain equally by rows and columns to obtain a number of sub-regions. If there are real fuel rods and sub-channels in the region, the sub-region is valid, otherwise, the sub-region is invalid. Step 4: Remove the invalid area and number the valid sub-areas in sequence.

Images