Simulation Study of Flow Instability in Parallel Multi-Channel Systems Based on Modelica

In parallel channels of a nuclear reactor core, flow instability can cause a significant decrease in critical heat flux (CHF) or mechanical oscillation of the fuel components, endangering the normal operation of the reactor. Nuclear reactor Modeling and Analysis Platform ， NUMAP, developed based on the two-fluid six-equation theory and using the Modelica language, is a multi-domain unified modeling and simulation platform for nuclear power plants. In this paper, a parallel dual-channel system model was constructed based on the NUMAP, referencing a high-temperature and high-pressure steam-water two-phase thermohydraulic experimental device, to simulate flow instability phenomena. The comparison with experimental data validated the transient analysis ability of the NUMAP for flow instability phenomena. Based on this, the flow instability boundary of a parallel multi-channel system was calculated under the same operating conditions. When the number of parallel channels was 2, 3, and 4, the calculated flow instability boundary error did not exceed ±5%, verifying that a parallel dual-channel structure can be used to obtain the flow instability boundary when there are multiple parallel heating channels.


Introduction
Flow instability refers to the phenomenon where a system deviates from its stable state and experiences nonperiodic drift or periodic oscillation in thermal parameters after being subjected to instantaneous disturbances (Xu 2001).For systems with multiple parallel heating channels, flow instability may occur between the parallel channels, manifested as periodic flow pulsation between the heating channels.When flow instability occurs between parallel channels, it can significantly reduce the critical heat flux density or cause mechanical vibration of the equipment, thereby endangering equipment safety.
The mechanism of flow instability is shown in Figure1 which shows the pressure drop and flow rate curve of water in a straight pipe with constant heating power.Within the entire range of mass flow rate ( M ), when the flow rate decreases to the point E in a flow pipeline with constant total heating power, two-phase flow occurs due to the generation of a large amount of steam, leading to an increase in pressure drop ( P  ) in the pipeline as the flow rate decreases.At point D, the flow rate is small enough to form single-phase steam flow in the pipeline, and the D-O region is the pressure drop characteristic curve of single-phase steam.In Figure 1, we can see that when the differential pressure ( P  ) between the inlet and outlet is between point D and point E, for the same differential pressure ( P  ), there may be three flow rates in the channel, which may lead to flow instability.For parallel channels, periodic flow pulsations may occur in the flow rates in each channel in this situation.
Many scholars have conducted experimental research on flow instability in parallel channels.Cheng et al. (Cheng 2018) studied the flow instability of natural circulation rod bundle parallel channels through experiments.The study mainly focused on the experimental phenomenon and mechanism of pressure drop flow instability when the stabilizer was connected to the upstream of the heating channel, and the effect of inlet subcooling and heat flux density on flow instability was also studied.Tang Yu et al. (Tang 2014) conducted experimental research on flow instability in rectangular parallel channels under forced circulation conditions, considering the effects of system pressure, mass flow rate, and inlet subcooling on the flow instability boundary, and obtained a non-dimensional relationship for flow instability through experimental data fitting.
Currently, experimental methods are mainly used to study flow instability in parallel channels (Peng 2021;Zang 2014).However, experimental research has the disadvantages of long cycle length, limited experimental conditions, and inability to capture the internal mechanism of flow instability phenomena.Therefore, it is necessary to use simulation methods to simulate flow instability phenomena, study the flow instability mechanism under different operating conditions of parallel heating channels, and determine the flow instability boundary to provide a reference for experimental research.
In engineering practice, one-dimensional system simulation programs are commonly used to simulate and analyze flow instability phenomena in parallel heating channels.Since the fluid state in parallel heating channels under forced circulation changes dramatically, the flow path is generally divided into several control volumes.This method reduces the complexity of the calculation compared to the three-dimensional numerical simulation method for studying the flow and temperature fields in the flow path, and achieves the expression of the thermalhydraulic characteristics of flow instability, which is widely used in modeling and simulation analysis of flow instability (Qian 2014;Xiong 2015).Commonly used simulation programs are developed using Fortran, C language, Matlab/Simulink, or directly using nuclear engineering professional simulation software such as RELAP5, THEATRE (Xia 2011;Xing 2010).These programs and software use card-based modeling, and the modeling process lacks a graphical interface, making it difficult for users to obtain the composition of the system equipment, equipment connection relationships, and equipment parameter information intuitively.Graphical simulation software such as MMS in the United States, APROS in Finland, and STAR-90 in China have been widely used in the field of dynamic system simulation.Their drag-and-drop modeling and visualization of equipment parameters greatly reduce the workload of modeling, but the code of these commercial software is relatively closed and not conducive to users for secondary development according to their needs.
Modelica is a physical modeling language that is object-oriented, equation-based, and uses reusable hierarchical component models (Zhao 2006).Its numerical model is open-source and visible, and is written in equation form, which has strong comparability with empirical relationships obtained from experiments and theoretical equations.Users can efficiently and accurately construct and optimize numerical models using the advantages of the Modelica language.The Modelica language supports graphical drag-and-drop modeling, which can correspond one-to-one between the topological structure design of the numerical model and the system schematic diagram, greatly improving the modeling efficiency and model readability (Zhang 2022).Based on the advantages of the Modelica language, Suzhou Tongyuan Soft & Control Technology Co., Ltd. and China Nuclear Power Research and Design Institute have constructed a two-phase thermal-hydraulic characteristic model architecture, established a hierarchical structure from the basic control equations and closed equations to the system model, and developed a nuclear reactor unified modeling and analysis platform NUMAP based on Modelica, which is multi-disciplinary, multi-level, and multi-system (Huang 2021).

Two-Fluid Six-Equation Modeling Based on Modelica
The thermal hydraulic model library of NUMAP is described using the Modelica language by two-fluid sixequation formulation.The two-fluid six-equation model is a commonly used mathematical equation model for thermal-hydraulic system analysis of nuclear reactors, which can accurately simulate the mass and heat transfer behavior of the two-phase flow in normal and accident operating conditions of the nuclear reactor system.It can be used for nuclear reactor system design verification, operation simulation, safety analysis, and other scenarios.Based on the Modelica language, the two-fluid six-equation model considers the different properties, flow rates, temperatures, and mass, energy, and momentum exchange between the two-phase fluids in the actual twophase fluid flow process。Mass, momentum, and energy conservation equations are established separately for the vapor and liquid phases.In order to make the control equation group closed, constitutive equations such as interfacial friction, interfacial heat transfer, interfacial mass exchange, wall friction, and wall heat transfer are added.

Conservation equation model
The mass conservation equation of two phase： The momentum conservation equation of two phase： The thermal energy conservation equation of two phase： (3) where the following quantities and symbols appear:

Wall-to-Fluid Heat Transfer Model
the total wall heat flux ( q ) is the heat flux to the vapor plus the heat flux to the liquid.( ) ( ) where the following quantities and symbols appear:

Interfacial Heat Transfer Model
In NUMAP, the interfacial heat transfer between the gas liquid phase actually involves both heat and mass transfer.The form used in defining the heat transfer correlations for superheated liquid (SHL), subcooled liquid (SCL), superheated gas (SHG), and subcooled gas (SCG) is that for a volumetric heat transfer coefficient (W/(m 3. K)).Since heat transfer coefficients are often given in the form of a dimensionless parameter (usually Nusselt number, Nu), the volumetric heat transfer coefficients are coded as follows: where the following quantities and symbols appear:

Interphase Friction Model
The interface friction in the phasic momentum equations ( 2) is expressed in terms of phasic interfacial friction coefficients as： ( -)

Wall Friction Model
The pressure loss of fluid flow includes wall frictional pressure drop, gravity pressure drop, and accelerated pressure drop.The wall frictional pressure drop is calculated by The Wall Friction Model which include Frictional resistance of single-phase flow and two-phase flow.
(1) single-phase flow Darcy's formula is used for the frictional pressure drop of single-phase flow where f is the frictional resistance factor.(2) two-phase flow The HTFS correlation is used to calculate the twophase friction multipliers.This correlation was chosen because it is correlated to empirical data over very broad ranges of phasic volume fractions, phasic flow rates and The HTFS correlation for the two-phase friction multiplier is expressed as: where C is the correlation coefficient and 2 X is the Lockhart-Martinelli ratio.
In order to calculate the frictional pressure drop of two-phase flow ， NUMAP represents the two-phase pressure drop as the product of the single-phase frictional pressure drop and the HTFS correlation.

Water and Steam Property Calculation Models
In thermal hydraulic simulation software, it is necessary to accurately and quickly calculate all thermodynamic parameters of water and steam.When calculating the thermodynamic parameters of water and steam through fitting formulas, high-order polynomials are necessary to improve calculation accuracy, but highorder polynomials can slow down the calculation speed and consume time when the number of control volume is large.
In order to improve calculation accuracy and speed, NUMAP uses the lookup table method to solve the water and steam parameters, using the 1967 ASME Steam Table .The water vapor table consists of three tables, namely the steam property table (two-dimensional), the water property table (two-dimensional), and the saturated property table (one-dimensional).

Discrete solution of equations
The conservation equations and closed constitutive equations form a nonlinear equation system for solving thermal-hydraulic problems.In the process of solving the equation system, numerical methods are used to solve the coupled "pressure-velocity" relationship in the model.In the discretization process, a semi-implicit method is used, which was proposed by Patankar and Spalding in 1972 and is commonly known as the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) method, which is a semi-implicit method for solving the pressure coupling equation.The semi-implicit solution method is based on the coupling of velocity and pressure on a staggered grid (Tao 2001).
For one-dimensional fluids, a staggered grid refers to storing the velocity and pressure separately on two different grids.The pressure is stored at the control volume center (node), and the velocity is stored at the interface between two adjacent control volumes.When calculating the velocity correction value, the influence of the velocity correction value on the adjacent interface is not considered, that is, the velocity correction value on the interface is only determined by the pressure correction value between adjacent control volume nodes.The parallel heating dual-channel experiment was conducted on a high-temperature and high-pressure twophase thermal-hydraulic experimental device, which has a design pressure of 17.2 MPa, a design temperature of 350℃, and an experimental mass flow rate range of 100-5000 kg/(m 2 .s)(Wang 2021).Referring to the experimental platform, a parallel heating dual-channel system model was built based on NUMAP, and the system model is shown in Figure 5.The way to carry out simulation calculation conditions is to gradually increase the heating power in a step-by-step manner while keeping the system pressure and inlet parameters unchanged.After each power increase, the heating power is kept constant for 2 minutes.When the flow instability occurs, the heating power is kept constant, and the quality gas rate at the channel outlet is recorded.Experiments and simulation calculations were conducted using 10 groups of uniform heating conditions and 16 groups of non-uniform heating conditions.The comparison between the simulation results and the experimental results is shown in Figures 6 and 7.   Based on the reference experimental conditions, simulations were carried out for the flow instability of uniform and non-uniform parallel heating dual-channel systems.The simulation results of the quality gas rate at the channel outlet when the flow instability occurs were selected and compared with the experimental values.In terms of calculation error, the maximum errors between the simulation results and the experimental results for uniform and non-uniform heating conditions were 18.87% and 15.7% respectively, which met the error requirement of 20% for this experimental condition; In terms of calculation speed, the simulation setting time for both uniform and non-uniform heating test cases is 400 seconds, and the actual average simulation time is 356 seconds, meeting the real-time requirements of simulation calculation.This indicates that the NUMAP can be used for simulation analysis of the transient characteristics of flow instability in parallel channels.

The effect of different numbers of heating channels on the unstable boundary of parallel multi-channel flow.
Parallel multi-channel systems exist in various nuclear power equipment, and their flow instability can lead to a significant decrease in critical heat flux (CHF) or mechanical oscillation of the equipment, affecting the normal operation of nuclear power plants.Therefore, how to avoid the occurrence of flow instability in parallel multi-channel systems during the operation of nuclear power plants has always been a hot research topic.After verifying the transient calculation capability of NUMAP for the flow instability of parallel heating dual-channel flow, this section will simulate the flow instability boundary of parallel multi-channel flow.

Calculation method for flow instability boundary.
Ishii (Ishii 1942) showed that the subcooling number Nsub and the dimensionless number of phase change Npch are the two most important dimensionless numbers for characterizing flow instability under forced circulation conditions.Most researchers prefer to use the Npch -Nsub diagram to plot the flow instability boundary and region, and this method is also used in this study.
The dimensionless subcooling number Nsub and the dimensionless phase change number Npch and are defined as follows:     (2) Simulation study of a parallel multi-channel heating model was conducted using the NUMAP.The values of subcooling degree number Nsub and phase change number Npch at the onset of flow instability were calculated for 2, 3, and 4 heating channels.The difference of the flow instability boundary was within ±5%.This indicates that when there are multiple parallel heating channels, a parallel dual-channel structure can be used to obtain the flow instability boundary.Compared to traditional nuclear thermal hydraulic software, the NUMAP developed based on Modelica language has the following advantages: in terms of user experience, the drag and drop modeling method has greatly improved interactivity and intuitiveness in the modeling process; In terms of post processing of results, NUMAP provides real-time monitoring of variable parameter processes based on data points and XY curves, which can clearly and intuitively obtain the variation patterns of parameters and the linkage relationships of various parameters during the simulation process; In terms of application, the characteristics of Modelica multi domain unified modeling endow NUMAP with the ability of joint simulation in different fields of Mechanical, electrical, hydraulic and control systems.It provides not only thermal hydraulic system model library for nuclear power direction simulation, but also mechanical, electrical, control and other multi domain model libraries.

prospects
In the next stage, we will use Modelica language to continuously develop nuclear domain model libraries on the NUMAP, such as nuclear Nuclear reactor core model, pump model, steam generators model ， Steam turbine model and other thermal hydraulic models, as well as nuclear instrumentation and control, mechanical and electrical system models.We will build NUMAP into a large multi discipline, multi-level and multi system unified modeling and simulation platform based on Modelica that promote high-quality research and development of nuclear power systems.

Figure 1 .
Figure 1.Characteristic curves of pressure drop and flow rate in the heating channel the thermal-hydraulic conservation equations used by NUMAP are presented.To solve these equations in a mathematically closed manner, sufficient closed constitutive equation models need to be supplemented, including Flow Regime Maps Model, Wall-to-Fluid Heat Transfer Model, Interfacial Heat Transfer Model, Interphase Friction Model, Wall Drag Model, water and steam property calculation Models, etc.(Rockville 2001).
Two flow-regime maps for two-phase flow are used in the NUMAP: a horizontal map and a vertical map for flow regime in pipes.(1)The horizontal flow regime mapThe horizontal flow regime map is for the flow pipeline whose elevation angle  is such that 0 45 Simulation Study of Flow Instability in Parallel Multi-Channel Systems Based on Modelica 570 Proceedings of the Modelica Conference 2023 October 9-11, 2023, Aachen, Germany DOI 10.3384/ecp204569 degrees.The map contain Bubbly Flow, Slug Flow, Annular Mist Flow, Mist Flow, Horizontal Stratification Flow and Transition Flow pattern.

Figure 2 .
Figure 2. Schematic of horizontal flow regime map with hatchings, indicating transition regions.(2) The vertical flow regime map The vertical flow regime map is for the flow pipeline whose elevation angle  is such that 45 90   degrees.The map contain Inverted Annular Flow, Inverted Slug Flow, Mist Flow, Bubbly Flow, Slug Flow, Annular mist Flow, Vertically Stratification Flow and Transition Flow pattern.

Figure 3 .
Figure 3. Schematic of vertical flow regime map with hatchings, indicating transition regions.

F
where ig is the magnitude of the interfacial friction force per unit volume on the vapor and if F is the magnitude of the interfacial friction force per unit volume on the liquid.FIG and FIF are the phasic interfacial friction coefficients.
phasic flow regimes.The correlation has also been shown to give good agreement with empirical data.

Figure 4 .
Figure 4. Storage of pressure and velocity.3 Development and engineering validation of a parallel channel system model 3.1 Development and validation of a parallel dual-channel system model.

Figure 5 .
Figure 5. Model of parallel heating dual-channel system.The way to carry out simulation calculation conditions is to gradually increase the heating power in a step-by-step manner while keeping the system pressure and inlet parameters unchanged.After each power increase, the heating power is kept constant for 2 minutes.When the flow instability occurs, the heating power is kept constant, and the quality gas rate at the channel outlet is recorded.Experiments and simulation calculations were conducted using 10 groups of uniform heating conditions and 16 groups of non-uniform heating conditions.The comparison between the simulation results and the experimental results is shown in Figures6 and 7.

Figure 6 .
Figure 6.Calculation error of uniform heating.

Figure 7 .
Figure 7. Calculation error of non-uniform heating.Based on the reference experimental conditions, simulations were carried out for the flow instability of uniform and non-uniform parallel heating dual-channel systems.The simulation results of the quality gas rate at the channel outlet when the flow instability occurs were selected and compared with the experimental values.In terms of calculation error, the maximum errors between the simulation results and the experimental results for uniform and non-uniform heating conditions were 18.87% and 15.7% respectively, which met the error requirement of 20% for this experimental condition; In terms of calculation speed, the simulation setting time for both uniform and non-uniform heating test cases is 400 seconds, and the actual average simulation time is 356 seconds, meeting the real-time requirements of simulation calculation.This indicates that the NUMAP can be used for simulation analysis of the transient characteristics of flow instability in parallel channels.
quantities and symbols appear: Q the total heating power, W.W the total mass flow rate at the inlet of the channels, kgThe flow instability boundary of parallel dual, triple, and quadruple channel systems was studied using the experimental conditions selected in Section 3.1.The typical simulation curves obtained are shown in Figures 6

Figure 8 .
Figure 8.Typical flow rate curve for parallel dual-channel system with symmetric uniform heating.

Figure 9 .
Figure 9.Typical flow rate curve for parallel triplechannel system with symmetric uniform heating.

Figure 10 .
Figure 10.Typical flow rate curve for parallel quadruplechannel system with symmetric uniform heating.In the simulation of parallel heating multi-channel system, as the heating power increases, two-phase flow of gas and liquid begins to appear when the outlet flow of the heating channel reaches saturation.As the outlet gas content of two-phase flow gradually increases, the flow rate of each channel in the parallel multi-channel system shows periodic flow pulsation when the outlet gas content reaches a certain level, which indicates the occurrence of flow instability, as shown in the local enlarged figures in Figures8-10.Under different operating conditions, the values of the subcooling degree number Nsub and phase change

Figure 11 .
Figure 11.Flow instability boundaries under different heating channel numbers 4 Conclusion (1) Flow instability is one of the important research topics in nuclear reactor thermal safety, and simulation analysis can serve as an important means for conducting flow instability research.In this paper, a parallel dualchannel heating model was built using the NUMAP developed based on Modelica language.Simulations and analysis of flow instability by this model were conducted.The simulation results met the error requirements of the experimental conditions, verifying the transient calculation and analysis ability of NUMAP for flow instability phenomena.(2)Simulation study of a parallel multi-channel heating model was conducted using the NUMAP.The values of subcooling degree number Nsub and phase change number Npch at the onset of flow instability were calculated for 2, 3, and 4 heating channels.The difference of the flow instability boundary was within ±5%.This indicates that when there are multiple parallel heating channels, a parallel dual-channel structure can be used to obtain the flow instability boundary.

Table 1 .
Wall convection heat transfer models Simulation Study of Flow Instability in Parallel Multi-Channel Systems Based on Modelica Simulation analysis software is one of the indispensable core technologies for Nuclear Reactor System research and development, and also an important cornerstone of Nuclear Reactor System innovation.The NUMAP is developed by giving full play to the characteristics of Modelica physical Modeling language's Simulation Study of Flow Instability in Parallel Multi-Channel Systems Based on Modelica 574 Proceedings of the Modelica Conference 2023 October 9-11, 2023, Aachen, Germany DOI 10.3384/ecp204569 object-oriented, equation based, reusable and reconfigurable hierarchical component model technology and its advantages in multi physical and multi-scale collaborative simulation。