Mathematical modeling and numerical investigation of density wave instability of supercritical natural circulation loop

In present paper, a mathematical model based on the one dimensional nonlinear mass, momentum and energy conservation equations has been developed to study the density wave instability (DWI) in horizontal heater and horizontal cooler supercritical water natural circulation loop (HHHCSCWNCL). The one dimensional nonlinear mass, momentum and energy conservation equations are discretized by using finite difference method (FDM). The numerical model is validated with the benchmark results (NOLSTA model). Numerical simulations are performed to find the threshold stability zone (TSZ) and draw the stability map for natural circulation loop. Further, effect of change in diameter and riser height on the density wave instability of SCWNCL has been investigated.


Introduction
Supercritical water reactor (SCWR) is one of the most important technologies among the Generation-IV forum in the nuclear reactor. SCWR is the very potential. Their unique features and design after various advantages over the boiling water reactor. For example, the thermal efficiency of the SCWR is close to 44-45%, which is much higher than the thermal efficiency of 30-32 % of BWR [1][2][3][4]. In the SCWRs, the size of the turbine is compact due to high steam enthalpy which reduced the capital cost of the load. Due to their exceptional properties such as heat removing capacity, better fuel utilization, and minimum capital cost to generate electricity [5]. It is produced clean and economical energy compared to the current water reactor. SCWRs will experience drastic change in the density along the core, therefore there is lot of chances to occur instability in the SCWRs, which is similar to the BWRs insatiability namely static, Ledinegg, dynamic, density wave oscillation and pressure wave oscillation instability etc [6][7][8][9][10][11][12]. These instabilities are not desirable in the reactor because failure of control and damage to the wall of the reactor. So, one should be required to show their skill to predict the instability in the reactor. Accordingly, it is very important to understand the mechanism of instability in the reactor and recognize the variables which affect the instability. Before going to next, first we want to explain here the flow instability in the reactor. When providing the perturbation in the mass flow rate in the SCWNCL , if the mass flow rate come back to the first state or next stable point, where the oscillation of mass flow rate is not growing with respect to time, they system is known as stable. if the oscillation of mass flow rate is continuous growing with respect to time, such system is known as unstable, whereas the oscillation of mass flow rate neither growing nor shrinking, but it oscillate with constant amplitude with respect to time, such system is known as neutral stable [13][14][15][16][17][18][19][20][21][22].
In this article, the primary work is validated the present model with the benchmark problem and search the existence of dynamic wave instability of a SCWNCL. The density wave instability (DWI) can be affected by the various parameters of the loop such as diameter and riser heights, therefore, the main works are draw the threshold stability zone (TSZ) of the loop for DWI and incorporate the effect of them.

Initial and boundary conditions
The initial and boundary circumstances influence the solution of balancing coupled equations. For transient conservation equations, the steady state solution is utilised as the starting condition. The applied boundary conditions at the inlet and exit for the steady state and transient study are known as constant pressure and inlet specific enthalpy. To maintain a zero pressure drop in the loop, the inlet and outlet of the loop are linked in a huge tank with constant intake pressure and inlet temperature.

development of the mathematical model
The present mathematical model is a simple one dimension nonlinear coupled balance mass, momentum and energy equations, which properties are varying along the axial direction. These balance equations are solved numerically [27][28][29][30][31][32][33].

Equation of state
To avoid loop discontinuity, the connection between fraction factor and dimensionless Reynolds number is provided below to estimate the friction factor for laminar and turbulent flow in the loop. The essential Reynolds number is compared to the calculated value to determine if the flow is laminar or turbulent, and then the friction factor for each node is calculated.

Steady state numerical solution methodology
In the governing balance coupled nonlinear steady state equation (1), the temporal terms are first eliminated and then space terms are discretizated by using finite difference technique. The obtained discretizated equations are given below:

Transient numerical solution approach
For solving the transient governing equations, first converted the above equation (1) into a primitive form, which is given below as equation (8) > @ > @ > @ The characteristics of equation (8) depend on the eigenvalues of the matrix A . The eigenvalues of A are found to be real and distinct, and the equation (8) identified as a hyperbolic in nature. Further, the equation (8) is modified in a characteristic form which is given in equation (9).
Where, / is diagonal element of eigenvalues of matrix A .
The governing equation (9) is discretized by using finite difference technique. The discretized equations are given below. All the set of equations are consisting spatial and temporal terms, the spatial derivative terms of first two equations are discretized by backward difference and third equation is discretized by forward difference technique. All discretized equations are combined together and apply Bc's and solve it till convergence is found [30].

Comparison of present model against literature result
The present numerical mathematical model has been compared against the literature result, which is given by Sharma [20]. The obtained result by code is compared with the NOLSTA code for HHHC loop reported by M Sharma [20].

Results and discussions
The unsteady steady instability of SCWNCL has been investigated. Also, the effect of geometry (such as diameter and riser height) and operating parameters like (pressure and frication factor) on unsteady steady instability of SCWNCL have been examined.

Density wave instability
The density wave instability (DWI) has been determined by using the Von Neumann stability criterion. The amplitude of mass flow rate has been calculated for a particular input such as applied heater power, inlet pressure and inlet specific enthalpy with respect to time. The ratio of mass flow rate of new time step to the previous time step shows the stability of the SCWNCL. The stability criterion is given below: If the value of G is equal to the one for all time step, this phenomena is known as neutral stable and when the value of G increased and go beyond the one and continuously increase its value, this phenomena is known as unstable, whereas the value of G is less than one and continuously decrease its value, this phenomena is known as stable system. The DWI criterion is further given below: A simulation has been carried out to find the stable, unstable and neutral point for three different applied heater powers 21, 24.08 and 26 kW, inlet specific enthalpy (1500 kJ/kg) and inlet system pressure (25 MPa). The geometry parameter has been kept constant for all three simulations such as diameter 14 mm and riser height 4.1m. From the figure 3, it can be observed that the amplitude of mass flow rate is growing against the time at heater power 26.0 kW, this system becomes unstable, while at heater power 21.0 kW, and its value reduced against the time and become the stable system. At applied heater power 24.08 kW, the amplitude of mass flow rate is neither growing against the time nor reduced against the time but amplitude of mass flow rate is constant against the time, this system is known as neutral stable. Next, all the neutral point has been joined together and draw the stability map for loop.

Effect of loop diameter
The SCWNCL flow instabilities highly depend on the loop diameter. In fact, at a larger diameter, the mass flow rate was observed as high compared to the smaller diameter. It can be observed that from figure 4, the MSB's of DWI for SCWNCL is increased with an increase in the loop diameter.

Effect of riser height
The effect of riser height on DWI of the SCWNCL is depicted in figures 5 for three different riser height of 4.1m, 52.1m and 6.1m for keeping other parameters constant such as diameter, inlet pressure and friction factor. From the figure 9, it can be observed that the marginal stability boundary of DWI of SCWNCL is increased with increase in the riser height. The MSB of DWI of SCWNCL has been significantly increased with increase in the riser height 4.1m to 5.1 m, whereas the rate of MSB of DWI at 6.1 m is increased less compared to the 5.1m.

Conculsions
To accurately identify the zone of density wave instability at supercritical temperature and pressure, a one dimensional nonlinear thermal-hydraulic (TH) model was developed for a natural circulation loop in FORTRAN 90. The TH model first verified with NOLSTA code and next various simulations carry out to analyze the DWI of the SCWNCL. The present work identified the occurrence of DWI in the loop. The nonlinear TH model investigation reveals that smaller diameter of the NCLs are less stable than bigger diameter of the NCLs with respect to power for SCWNCL. Based on the above results, the following conclusions have been drawn which is given below: x Increasing in diameter of the loop, the stable zone of DWI is increased. x Increasing in riser height of the loop, the stable zone of DWI is increased.