Nuclear Engineering and Design 342 (2019) 10–19

Contents lists available at ScienceDirect

Nuclear Engineering and Design

journal homepage: www.elsevier.com/locate/nucengdes

Numerical study on the density wave oscillation of supercritical water in T

parallel multichannel system
⁎
Jialun Liu, Huixiong Li , Qian Zhang, Xiangfei Kong, Xianliang Lei
State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China

A R T I C LE I N FO A B S T R A C T

Keywords: A time-domain model is proposed to study the density wave instability of supercritical water in parallel mul-
Parallel multichannel system tichannel system. The present model splits the flow domain along each channel based on the collocated mesh
Density wave oscillation method, and adopts an iterative solution to solve the coupling among multiple channels. The model is verified
Supercritical water firstly by the experimental data obtained in open literature. After that, the model is used to study the effect of the
Time-domain method
key parameters on the density wave instability of supercritical water in parallel multichannel system (such as,
Critical inlet fluid temperature
adding different perturbations, inlet fluid temperature, and so on). Moreover, the influence mechanisms of these
key parameters on the density wave oscillation of supercritical water are also explained. It is found that out-of-
phase oscillation occurs between the perturbation channels (which are added with the perturbation) and the rest
other channels, and the amount of perturbation channels has little effect on oscillation period and system sta-
bility. The system stability first decreases and then increases with the inlet fluid temperature, and the system has
the worst stability when the inlet fluid temperature is equal to a critical value, denoted by Tin-cr. This is because
that the variation of the fluid density along the system will reach the maximum when the inlet fluid temperature
is near Tin-cr, which makes the system most unstable. Tin-cr is always lower than the pseudo-critical temperature
of supercritical water at the corresponding pressure. Tin-cr has no obvious change with system pressure.

1. Introduction 2011; Xiong et al., 2012; Guo et al., 2010; Xi et al., 2014a; Yuan et al.,
2012; Papini et al., 2014, and so on), the theoretical research by the
Parallel channels are widely applied to transfer heat to working analysis models (Xiong et al., 2013, Zhang et al., 2009; Lu et al., 2014;
fluid in such industrial equipment as boiler heater surfaces, steam Su et al., 2013; Hu et al., 2015; Papini et al., 2012; Muñoz-Cobo et al.,
evaporators, water reactors, refrigerating systems and various heat 2002, and so on), and the numerical research by RELAP5 code or CFD
exchangers, due to their high cycle performance and small equipment simulation (Colombo et al., 2012; Papini et al., 2012; Xi et al., 2014b; Li
size (Matsuyama et al., 1995; Ablanque et al., 2010). However, many et al., 2015, Li et al., 2018, and so on). These above researches focused
kinds of flow instability phenomena always occur and threaten the on the characteristics of the DWO in two parallel channels. However,
normal running of the parallel channel system. As one of the most most industrial equipment (such as the reactor cores, water-wall of
common flow instabilities, density wave oscillation (DWO) is char- boiler, etc.) is composed of multiple channels and the amount of
acterized by continuous oscillation of such operation parameters as the channels is always huge. Moreover, the relevant study had indicated
mass flux, the fluid temperature, tube wall temperature and so on, that the amount of channels and the difference in working conditions
which may further induce the mechanical vibrations and the thermal among channels may produce great effects on the DWO and system
fatigue of the channel structure, and also results in problems in con- stability due to strong interactions existing among channels (Aritomi
trolling of the system (Kakaç and Cao, 2009). Therefore, it is necessary et al., 2012). Unfortunately, scarce models and experimental research
to study and prevent the DWO in parallel channel system. were found for the DWO in parallel multichannel system, and this may
In the past, many researchers had carried out a lot of research work be because that more channels bring some difficulties in solving the
including the numerical studies and experimental studies on the DWO complex coupling among channels, and also increase the difficulty of
in parallel channel system. Among these past researches, it was found experimental work. In the existing studies on the DWO in multiple
that most of them focused on the DWO between two channels, in- channels, Li et al. (2005) conducted an experiment about the DWO of
cluding the experimental research (Gao et al., 2005a,b; Huang et al., high-pressure steam-water two-phase flow in parallel three channels,

⁎
Corresponding author.
E-mail address: huixiong@mail.xjtu.edu.cn (H. Li).

https://doi.org/10.1016/j.nucengdes.2018.11.014
Received 25 April 2018; Accepted 19 November 2018
Available online 03 December 2018
0029-5493/ © 2018 Elsevier B.V. All rights reserved.
J. Liu et al. Nuclear Engineering and Design 342 (2019) 10–19

Nomenclature QT threshold heat power for the occurrence of DWO, kW

q heat power, kW
A cross-sectional area, m2 Re Reynolds number
Cf resistance coefficient T fluid temperature, °C
D pipe inner diameter, m t time, s
f frictional resistance coefficient z the axial distance from the pipe inlet, m
G mass flow velocity, kg/(m2·s) ΔP pressure drop, MPa
g gravitational acceleration, m/s2
h fluid enthalpy, kJ/kg Greek letters
i sequence number of channels
j sequence number of calculated sections ρ density, kg/m3
K local resistance coefficient θ inclination angle with respect to horizontal direction
k time layer
L pipe length, m Subscripts
M mass flow, kg/s
NB amount of total parallel channels in inlet
NC amount of total sections contained in one channel max maximum
P pressure, MPa min minimum
Q heat flux per length transferred to the fluid, kW/m out outlet

and researched the influences of the system pressure, inlet sub-cooling, study:
and that of volumetric compressibility on instability. Yun et al. (2008)
built a calculation model by the method of nonlinear analysis to study (1) The flow in each channel is regarded as one-dimensional and
the two-phase flow instability of a multichannel system. Xia et al. homogeneous.
(2012, 2016) studied the two-phase flow instabilities of parallel narrow (2) There only exists radial heat transfer between the working fluid and
multi-channel system using RELAP5/MOD3.4 codes. Lee and Pan the pipe wall.
(1999) extended the nonlinear analysis model for a single channel de- (3) The heat load is uniformly distributed along the axial direction.
veloped by Clausse and Lahey (1990) to a multichannel system by using (4) The effect of the viscous dissipation, kinetic energy, and potential
the Galerkin nodal approximation method, and then used the model to energy on the energy conversation equation is ignored.
study the two-phase flow instability of a boiling water nuclear reactor.
Dutta et al. (2015) also extended the 1-D thermal-hydraulic model Due to the fact that the distinct difference between the liquid phase
(THRUST) for single channel to analyze DWO in parallel channels in the and the gas phase of the fluid disappears at supercritical pressures, the
reactor core. Even so, the past researches on the DWO in parallel supercritical pressure fluid can be considered as a single-phase fluid
multichannel system are still very limited, especially the situation at with no obvious phase transition, but the thermo-physical properties of
supercritical pressure conditions. Moreover, the influence mechanisms the fluid may change remarkably along the heated channels in the su-
of the factors (such as inlet fluid temperature, system pressure, etc.) on percritical pressure cases. The flow and heat transfer in each channel
the DWO of supercritical fluid in parallel multichannel system are also are based on the following one-dimensional homogeneous unsteady N-S
needed to be further discussed. equations.
In view of this situation, a time-domain model suitable for the DWO
of the supercritical water in parallel multichannel system (the amount ∂ρ ∂M
A + =0
of channels is more than 2) is established. The present model splits the ∂t ∂z (1)
flow domain along each channel into a series of collocated grids, and
adopts an iterative solution to solve the coupling among multiple
channels. The model is then verified by the experimental data obtained
in open literature. After that, the effects of the factors such as adding
the perturbation in different ways, system pressure, and inlet fluid
temperature on the DWO in the parallel multichannel system are in-
vestigated under supercritical pressure conditions. The results in this
paper can provide a theoretical guidance for the design and safe op-
eration of the parallel channel system in the boilers, reactors, etc.

2. Mathematical model

Fig. 1 shows the schematic diagram of a parallel multichannel

system studied in this paper. As shown in Fig. 1, supercritical pressure
water first entries the system from a lower plenum, and then flows
upward into the heated parallel channels where the fluid property
changes along the flow direction. The total number of channels is set as
NB. Finally, the fluid exits the system from the outlet of an upper
plenum.

2.1. Basic equations and assumptions

For ease of analysis, the following assumptions are made in this Fig 1. The schematic diagram of a parallel multichannel system.

11
J. Liu et al. Nuclear Engineering and Design 342 (2019) 10–19

∂M ∂ ⎛ M2 ⎞ ∂P M2 A ⎛ [ρ ·h](i, j, k + 1) + [ρ ·h](i, j − 1, k + 1)
+ ⎜ ⎟ + A + Cf · + ρgA sin(θ) = 0
∂t ∂z ⎝ ρA ⎠ ∂z ρA (2) Δt ⎝ 2
[ρ ·h](i, j, k ) + [ρ ·h](i, j − 1, k ) ⎞
−
∂ (ρh) ∂ (Mh) 2 ⎠
A + = Ql
∂t ∂z (3) [M·h](i, j, k + 1) − [M·h](i, j − 1, k + 1)
+ = Q (i. j. k + 1)
where Cf is the pressure drop resistance coefficient, which considers the Δz (9)
frictional resistance loss and local resistance loss. Cf is determined by:
2.3. Boundary conditions
f K δ (z ) K δ (z − L)
Cf = + in d + out d
2D 2 2 (4)
The boundary conditions used in this paper are given below:
where δd is dimensional Dirac delta function, and f is frictional re-
sistance coefficient. For supercritical flow, f can be calculated by Eq. (5) (1) The inlet pressure of the system is given, i.e., Pin = constant;
according to the literature (Filonenko, 1954). (2) The total mass flow of all channels is given, i.e., Mtotal = constant;
(3) The heat load on each channel is given, i.e., Q(i) = constant, i = 1,
1
f= 2, 3, …, NB;
(1.82lg Re−1.64)2 (5)
(4) The inlet fluid temperature of the system is given, i.e.,
where Re is the local Reynolds number. Tin = constant;
The thermo-physical properties of supercritical water are calculated (5) Outlet pressure of all channels are equal.
by the following state equation based on the IAPWS Industrial
Formulation 1997. 2.4. Solution algorithm

ρ = f (P , h) (6) The coupling among channels become more complex in the system
with a larger amount of channels, and an iterative solution is designed
for the present model. An overview of the solution algorithm is given
2.2. Discretized governing equations below.

Based on the collocated mesh method, the flow domain along each (1) The initial values of fluid parameters (such as M, P, h, etc.) at time
channel is split into a series of sections, and the fluid parameter (such layer (k) for all channels are obtained by solving the steady-state
as: the fluid pressure, density, temperature, mass flow, etc.) of each equations;
section is stored in the center of each section. Fig. 2 shows the dis- (2) A perturbation is imposed on the heat load of the channels, and the
cretization of the flow domain along a channel numbered as channel i calculation for the next time layer (k + 1) is started as below.
(i = 1, 2, 3,…, NB). (3) The inlet mass flow of channel i (i = 1, 2, 3,…, NB-1) is assigned by
In Fig. 2, the flow domain along the channel i is split into a series of a guessed value M(i, 0, k + 1), and the inlet mass flow of the last
sections (the total number of sections contained in each channel is NC).
To facilitate the modelling process, X(i, j, k) is defined to denote the
fluid parameters at the outlet interface of section j in channel i, and X
may denote the fluid pressure P, or the fluid temperature T, or the fluid
enthalpy h, or the mass flow M, or the fluid density ρ, etc. Here, k de-
notes the time layer. Q(i, j, k) denotes the heat load in section j of
channel i at time layer k. By integration of the conservation equations
(Eq. (1), Eq. (2), and Eq. (3)) from the inlet interface to the out interface
of section j and forward-difference approximation for the time deriva-
tive (Tao, 2010), the discretized equations of section j in channel i are
given below.
Discretized mass equation:

A ⎛ ρ (i, j, k + 1) + ρ (i, j − 1, k + 1) ρ (i, j, k ) + ρ (i, j − 1, k ) ⎞

−
Δt ⎝ 2 2 ⎠
M (i, j, k + 1) − M (i, j − 1, k + 1)
+ =0
Δz (7)

Discretized momentum equation:

1 ⎛ M (i, j, k + 1) + M (i, j − 1, k + 1) M (i, j, k ) + M (i, j − 1, k ) ⎞

−
Δt ⎝ 2 2 ⎠
1 ⎡ ⎛ M (i, j, k + 1)2 ⎞ ⎛ M (i, j − 1, k + 1)2 ⎞ ⎤
+ ⎜ ⎟ − ⎜ ⎟
AΔz ⎢ ⎣ ⎝ ρ (i, j, k + 1) ⎠ ⎝ ρ (i, j − 1, k + 1) ⎠ ⎦
⎥

P (i, j, k + 1) − P (i, j − 1, k + 1) Cf (i, j, k + 1)

+A + ·
Δz 2A
2 2
⎡ ⎛⎜ M (i, j, k + 1) ⎞⎟ + ⎛⎜ M (i, j − 1, k + 1) ⎞⎟ ⎤
⎢ ρ (i, j, k + 1) ⎥
⎣⎝ ⎠ ⎝ ρ (i, j − 1, k + 1) ⎠ ⎦
ρ (i, j, k + 1) + ρ (i, j − 1, k + 1)
+ ·g ·A·sin θ = 0
2 (8)

Discretized energy equation: Fig 2. Discretization of the flow domain along each channel.

12
J. Liu et al. Nuclear Engineering and Design 342 (2019) 10–19

channel can be obtained by:

NB − 1
M (NB, 0, k + 1) = Mtotal − ∑ M (i, 0, k + 1)
i=1 (10)

(4) The discretized equations of each section are calculated from

channel 1 to channel NB, and then the fluid parameters for all
sections of all channels (X(i, j, k + 1), i = 1, 2, 3,…, NB, j = 1, 2, 3,
…, NC) are obtained.
(5) According to the boundary condition, the outlet pressure of each
channel (denoted by Pout (i, k + 1), i = 1, 2, 3,…, NB) should be
equal to each other. To check that whether this boundary condition
is satisfied, the relative deviation of the outlet pressure of each
channel from the average outlet pressure (denoted by δP(i, k + 1))
is calculated by the following equation:
NB
∑i = 1 Pout (i, k + 1)
δP (i, k + 1) = Pout (i, k + 1) − , i = 1, 2, 3, ...,NB
NB
(11)
If δP(i, k + 1) for each channel is less than the specified tolerance
limit, the guessed values in step (3) are regarded as correct, and the
calculation is moved to the next time layer. Otherwise, the guessed
values in step (3) are modified and the steps (3)∼(5) are repeated until
the boundary condition is satisfied.

3. Model verifications

3.1. Grid refinement study

Fig. 4. Schematic diagram of the experimental system (Xiong et al., 2012).
To ensure independence of numerical results from temporal and
spatial grid, a grid refinement study is performed beforehand. During variation in the threshold heat power of the system with the time step is
the calculation, the amount of channels is set as for example of three, shown in Fig. 3(a). In the calculation, the density wave oscillations of
the length of each channel is 5.0 m, and the inner diameter of each the system with different heat powers are calculated, and the heat
channel is 0.015 m. The inlet total mass flow is 0.3 kg/s, the inlet fluid power where the system exhibits a sustained oscillation with constant
temperature is 280 °C, the fluid pressure is 25 MPa, and the inlet and amplitude is regarded as the threshold heat power under this inlet
outlet local resistance coefficients are both 1.0. condition. When the difference between two adjacent oscillation
Firstly, a temporal grid refinement study is discussed, and the

Fig 3. Threshold heat power variation with the time step and spatial step.

13
J. Liu et al. Nuclear Engineering and Design 342 (2019) 10–19

amplitudes is less than a given threshold (for example of 0.1%), the mass flow is 0.3 kg/s, and the heat power on each channel is equal to
oscillation process is considered to be a sustained oscillation with 24.285 kW.
constant amplitude. Here, threshold heat power (denoted by QT) is the It can be seen from Fig. 5(a) that out-of-phase oscillation occurs
sum of all channels’ heat power when the system presents a sustained between one perturbation channel (channel 1) and six rest other
oscillation with unchanged amplitude. It can be found in Fig. 3(a) that channels (channel 2, 3, 4, 5, 6, 7) after the perturbation is added on
threshold heat power almost has no change with the time step when the channel 1. It should be pointed that the rest other channels have the
time step is less than 0.02 s, and the time step can be set as 0.02 s in this same oscillation due to the same working condition, and the oscillation
study. Similar to the temporal grid refinement study, Fig. 3(b) shows curves of these channels are overlapped with each other in Fig. 5(a).
the variation of threshold heat power with the increase in the spatial The ratio of the oscillation amplitude of channel 1 (9.0E−5 kg/s) to
step, and the spatial step is set as 0.2 m in this study. that of each other channel (1.5E−5 kg/s) is about 6:1. This is mainly
because that the total mass flow is fixed in the out-of-phase oscillation,
3.2. Validation of the model and the sum of all perturbation channels’ oscillation amplitude should
be equal to the sum of all rest other channels’ oscillation amplitude. The
In this section, the present model is verified by the experimental similar oscillation characteristic can also be found in Fig. 5(b), and (c).
data collected by Xiong et al. (2012). Fig. 4 shows the corresponding Moreover, the system presents the same sustained oscillation in
channel system used in the experiment. Fig. 5(a), (b), and (c), indicating that the system stability has no change
In Fig. 4, the experimental system is consisting of the upper plenum, when different channels are added with perturbation. Besides this, the
lower plenum, and two vertical upward channels. The inner diameter of oscillation period also has no obvious change in different cases (7.22 s,
each channel is 0.006 m. Each channel is composed of an entrance 7.28 s, and 7.34 s). Therefore, just the way of adding a perturbation in
section (0.2 m), a heated section (3.0 m), and a riser section (0.35 m). Case (1-a) is selected in the following study, because the system sta-
The heated section is heated uniformly along the channel. An orifice is bility is independent of the above ways of adding perturbations.
set at the inlet and outlet of each channel. The inlet and outlet orifice
pressure drop coefficients of channel 1 is 5.4 and 4.93, respectively. The 4.2. Effect of inlet pressure on DWO in parallel multichannel system
inlet and outlet orifice pressure drop coefficients of channel 2 is 5.5 and
6.46, respectively. The total mass flux in the experimental cases is kept Fig. 6 shows the density wave oscillation curves of the inlet mass
as 0.033 kg/s. The threshold heat fluxes (Qc) with different inlet fluid flow with the inlet fluid pressure being 23, 25, and 27 MPa in three
temperatures and different pressures are calculated by the present parallel channels. The inlet total mass flow is kept as 0.3 kg/s, the inlet
model, and Table 1 shows the calculated results and the corresponding fluid temperature is kept as 280 °C, the heat power on each channel is
experimental results in all cases. kept as 58 kW.
From Table 1, it can be seen that the relative errors between the From Fig. 6, it can be obviously seen that, the oscillation curve of
calculated results and experimental results are less than 4% in all cases, the inlet mass flow presents a divergent trend at 23 MPa, a sustained
showing good agreement of the calculated results with experimental trend at 25 MPa, and a convergent trend at 27 MPa, which indicates
data, and confirming the capability of the present model in studying that the system stability increases with the pressure. This is mainly
DWO of the supercritical flow in the parallel channel system. because the fluid density varies more slowly with the fluid temperature
at a higher pressure, as shown in Fig. 7.
4. Results and discussion
4.3. Effect of inlet fluid temperature on DWO in parallel multichannel
In this section, three cases are designed to discuss the effect of system
adding different perturbations, system pressure, and inlet fluid tem-
perature on the DWO of supercritical flow in parallel multichannel Fig. 8 shows the threshold heat power with the inlet fluid tem-
system. In the calculation, the length of each channel is 5.0 m, and the peratures from 240 °C to 350 °C at 23 MPa, 25 MPa, and 28 MPa in three
inner diameter of each channel is 0.015 m. parallel channels, the inlet total mass flow is kept as 0.3 kg/s.
In Fig. 8, the system has a larger threshold heat power at a higher
4.1. Effect of adding different perturbations on DWO in parallel pressure, indicating that the system stability increases with the fluid
multichannel system pressure. Moreover, it is found that, similar to the situation in single
channel (Li et al., 2014) and two parallel channels (Xiong et al., 2013),
The perturbations may be added in various ways for the parallel the system stability first decreases and then increases with the inlet
multichannel system, and the effect of adding the perturbation in dif- fluid temperature, and there exists a critical value about the inlet fluid
ferent channels on the DWO is discussed firstly in the parallel multi- temperature, denoted by Tin-cr (called ‘Critical inlet fluid temperature’).
channel system (for example of seven channels). In this case, three ways When the inlet fluid temperature is equal to Tin-cr, the threshold heat
of adding perturbation are studied as below. power is reduced to the minimum, denoted by QTmin, and at this time

Case (1-a): the perturbation is added on only one channel (the heat Table 1
load on channel 1 is suddenly increased by 1% and then back to the The comparison between the calculated results and the experimental results.
original value at t = 10 s);
Case System Inlet fluid Calculated Experimental Relative
Case (1-b): the perturbation is added on four channels (the heat load pressure temperature Qc (kW/m2) Qc (kW/m2) error
on channel 1, 2, 3, 4 is suddenly increased by 1% and then back to (MPa) (°C)
the original value at t = 10 s);
1 23 180 620 600.07 0.0321
Case (1-c): the perturbation is added on six channels (the heat load
2 23 200 590 583.8 0.0106
on channel 1, 2, 3, 4, 5, 6 is suddenly increased by 1% and then back 3 23 220 570 579.6 0.0166
to the original value at t = 10 s); 4 24 200 598 592.2 0.0098
5 24 220 580 583.7 0.0063
Fig. 5 shows the oscillation curves of the inlet mass flow in parallel 6 24 240 550 571.1 0.0370
7 25 200 610 612.4 0.0039
seven channels calculated in Case (1-a), Case (1-b) and Case (1-c). In
8 25 220 600 609 0.0148
different cases, the other conditions are kept as the same, in which the 9 25 240 590 600.6 0.0177
pressure is 25 MPa, the inlet fluid temperature is 200 °C, the inlet total

14
J. Liu et al. Nuclear Engineering and Design 342 (2019) 10–19

Fig. 5. Oscillation curves of the inlet mass flow in seven parallel channels calculated in Case (1-a), (1-b) and (1-c).

the system has the worst stability. Besides, Tin-cr has almost no obvious stability of the system. However, the mechanism of Tin-cr in DWO of
change at pressures of 23 MPa, 25 MPa, and 28 MPa under this condi- supercritical water was rarely discussed before, and almost no criterion
tion. was proposed to predict Tin-cr under different working conditions. In
Clearly, it is needed to pay special attention to the working condi- this section, the mechanism of Tin-cr in DWO of supercritical water in
tions where the inlet fluid temperature is close to Tin-cr, due to the worst parallel multichannel system is explained below from the relation

Fig. 6. The oscillation curves of three parallel channels with the inlet fluid pressure being 23, 25, and 27 MPa.

15
J. Liu et al. Nuclear Engineering and Design 342 (2019) 10–19

enthalpy is near the pseudo-critical point, and this is mainly because

that the fluid density varies most sharply with the fluid enthalpy near
the pseudo-critical point. Accordingly, it can be expected that there
exists a value about the inlet fluid temperature where Δρ will achieve
the maximum. Fig. 11 further shows the variation of Δρ with the inlet
fluid temperature (Tin) at different supercritical pressures when Δh is
kept as a constant of 500 kJ/kg.
Fig. 11 shows that Δρ first increases and then drops down with the
increase in Tin under different pressures. Tin-max is used to denote the
corresponding inlet fluid temperature when Δρ reaches the maximum.
It can be found that Tin-max slightly increases with the pressure, and Tin-
max at different pressures is near 350 °C when Δh is equal to 500 kJ/kg.
Fig. 12 further shows the variation of Tin-max with different Δh at dif-
ferent pressures.
It can be seen from Fig. 12 that the Tin-max is always lower than the
pseudo-critical temperature under the corresponding system pressure,
and Tin-max is much closer to the pseudo-critical temperature with the
decrease in Δh of the system. It can be deduced that Tin-max will be equal
Fig. 7. The variation of water density with fluid temperature at different to the pseudo-critical temperature under the corresponding system
pressures. pressure when Δh approaches zero, and this also means that the fluid
density of supercritical water varies most sharply near the pseudo-cri-
tical temperature. The above inference is fully consistent with the
property variation of supercritical water, and this indirectly proves the
rationality of the analysis in this paper. Besides, Tin-max slightly in-
creases with the pressure.
Then, let’s analyze the relation between Tin-cr and Tin-max. When the
inlet fluid temperature reaches Tin-max, Δρ reaches the max value, the
fluid density has the most drastic variation along the system, and the
system should be in the most unstable condition. As discussed before,
the system has the worst stability when the inlet fluid temperature is
equal to Tin-cr. Therefore, it can be concluded that Tin-cr should be ap-
proximately equal to Tin-max, as shown in Fig. 13.
Then, the above theories are used to analyze the variation rule of
Tin-cr with the number of parallel channels. Fig. 14(a) and (b) further
shows the critical inlet temperature (Tin-cr) and the corresponding
minimum threshold heat power (QTmin) with the number of channels
being two, three, four, and five, under the pressures of 23 MPa, 25 MPa,
28 MPa, respectively. In the calculation, the inlet fluid mass flux is kept
the same as 0.3 kg/s. In Fig. 14, Tin-cr increases with the number of
Fig. 8. Threshold heat power of three parallel channels with the inlet fluid channels, while QTmin decreases with the number of channels. This can
temperature at different pressures. be explained below. Because the interactions among channels are more
complex with the increase in the number of channels, the stability of
parallel multichannel system generally decreases with the number of
between the fluid property variation and system stability.
channels, accordingly, QTmin decreases with the number of channels.
For a simple heated channel system as shown in Fig. 9, the occur-
Considering the total mass flow is kept as 0.3 kg/s in all calculation
rence of the DWO is mainly induced by the density variation of the
cases, so Δh of the system is proportional to the heat load on the system.
heated fluid along the system. The difference in density between the
Hence, with the increase in the number of channels, QTmin decreases, so
fluid entering the heated channel system and the fluid exiting (denoted
Δh of the system under the corresponding QTmin also decreases. As
by Δρ) is the key factor in the density wave oscillation, and the system
discussed in Fig. 12 and Fig. 13, Tin-cr is approximately equal to Tin-max,
always becomes more unstable with a larger Δρ (Yadigaroglu and
and the corresponding Tin-max is higher at a lower Δh of the system.
Bergles, 1972; Kakac and Bon, 2008). In general, Δρ is not only asso-
Therefore, Tin-cr increases with the number of channels. Besides, it is
ciated with the fluid enthalpy difference between the inlet and outlet of
also found that Tin-cr has no obvious change with the system pressure.
the system (denoted by Δh), but also associated with the inlet fluid
Fig. 15 further shows the comparison of Tin-cr and the corresponding
temperature of the system (denoted by Tin), as shown in Eq. (12).

Δρ = f (Δh, Tin ) (12)

Firstly, Δρ increases with Δh of a heated system, due to that the fluid

density generally decreases with the fluid enthalpy. Secondly, the
variation of fluid density with the fluid enthalpy is non-linear at the
supercritical pressures, and this causes that Δρ also changes with the
inlet fluid temperature Tin, even if Δh is fixed. Fig. 10 shows the fluid
density variation of supercritical water with the fluid enthalpy at
25 MPa.
It can be seen from Fig. 10 that Δρ is obviously different with the
inlet fluid enthalpy (or temperature) within the same Δh (for example of
500 kJ/kg). Moreover, Δρ is relatively larger when the inlet fluid Fig. 9. The schematic diagram of a simple heated channel system.

16
J. Liu et al. Nuclear Engineering and Design 342 (2019) 10–19

Fig 13. The relation between Tin-cr and Tin-max.

difference between Tin-max and Tin-cr is within 6% under all twelve

calculated conditions, showing good agreement with each other. This
confirms that the variation rule of Tin-cr is determined by the density
variation along the system, and the system has the worst stability when
the inlet fluid temperature is near Tin-cr is just because the fluid density
has the most drastic variation along the system at this working condi-
tion. Moreover, the variation rule of Tin-max (as shown in Fig. 12) can be
Fig. 10. The density variation with the fluid enthalpy at 25 MPa.
suitable for the prediction of Tin-cr in the parallel multichannel system.
This paper mainly discussed the Tin-cr phenomenon for supercritical
water. In fact, the Tin-cr phenomena also exists in subcritical water. The
mechanism of DWO is similar between the supercritical water and
subcritical water, and the occurrence of the DWO is mainly induced by
the density variation of the heated fluid along the system. Moreover,
the fluid density variation with fluid temperature is similar between the
supercritical water and subcritical water, as shown in Fig. 7. Therefore,
the rule of Tin-cr for subcritical water should be similar to that for su-
percritical water, and it can be obtained by the analysis on Δρ of the
system as discussed before in this paper.

5. Conclusions

In this paper, a time-domain model for studying the density wave

instability of the supercritical water in parallel multichannel system is
established. The present model is verified firstly with the experimental
data in the open literature, and then is used to study the characteristics
of the density wave oscillation in the parallel multichannel system at
Fig 11. The variation of Δρ with the inlet fluid temperature at different pres- supercritical pressure conditions. The following conclusions are ob-
sures. tained:

(1) Out-of-phase oscillation occurs between the perturbation channels

and the rest other channels, and the number of perturbation
channels has little effect on the oscillation period and the system
stability.
(2) The system stability increases with the system pressure.
(3) The system stability first decreases and then increases with the inlet
fluid temperature, and there exists a critical value about the inlet
fluid temperature, denoted by Tin-cr (called ‘critical inlet fluid
temperature’). When the inlet fluid temperature is equal to Tin-cr,
the system has the worst stability.
(4) The mechanism of Tin-cr in DWO of supercritical water in parallel
multichannel system is explained. With the increase of inlet fluid
temperature, the variation of the fluid density along the system will
reach the maximum when the inlet fluid temperature is near Tin-cr,
which makes the system most unstable.
(5) Tin-cr is always lower than the pseudo-critical temperature of the
supercritical water under the corresponding system pressure. Tin-cr
is much closer to the pseudo-critical temperature with the decrease
in the enthalpy difference between the inlet and outlet of the
Fig 12. The variation of Tin-max with Δh at different pressures.
system. Tin-cr has no obvious change with the system pressure.

Tin-max under those total twelve conditions given in Fig. 14. In Fig. 15, it
can be seen that both Tin-cr and the corresponding Tin-max present a Acknowledgements
consistent upward trend with the increase in the number of channels,
and Tin-max is slightly greater than Tin-cr. Moreover, the relative This work is supported by the National Key Research and
Development Program of China (Grant No. 2018YFB0604400).

17
J. Liu et al. Nuclear Engineering and Design 342 (2019) 10–19

Fig 14. The variation of Tin-cr (a) and the corresponding threshold heat power (b) with the amount of channels under the pressures of 23 MPa, 25 MPa, 28 Mpa.

phase flow instabilities in low-mass-flux water wall tubes of supercritical circular

fluidized bed boilers. Heat Transfer Eng. 32, 928–935.
Hu, L., Chen, D., Huang, Y., Yuan, D., Wang, Y., Pan, L., 2015. Numerical investigation of
the mechanism of two-phase flow instability in parallel narrow channels. Nucl. Eng.
Des. 287, 78–89.
Kakaç, S., Cao, L., 2009. Analysis of convective two-phase flow instabilities in vertical and
horizontal in-tube boiling systems. Int. J. Heat Mass Transf. 52, 3984–3993.
Kakac, S., Bon, B., 2008. A review of two-phase flow dynamic instabilities in tube boiling
systems. Int. J. Heat Mass Transf. 51, 399–433.
Lee, J.D., Pan, C., 1999. Dynamics of multiple parallel boiling channel systems with
forced flows. Nucl. Eng. Des. 192, 31–44.
Li, H.X., Wang, B., Xi, A., 2005. Experimental research on pulsating density instability of
high-temperature and high-pressure steam-water two-phase flow. Power Eng. 25,
55–59.
Li, J., Zhou, T., Song, M., Huo, Q., Huang, Y., Xiao, Z., 2015. CFD analysis of supercritical
water flow instability in parallel channels. Int. J. Heat Mass Transf. 86, 923–929.
Li, S., Chatoorgoon, V., Ormiston, S.J., 2018. Numerical study of oscillatory flow in-
stability in upward flow of supercritical water in two heated parallel channels. Int. J.
Heat Mass Transf. 116, 16–29.
Li, Y., Lu, D., Liu, Y., Zeng, X., 2014. Numerical analysis on supercritical flow instabilities
using APROS. Nucl. Eng. Des. 277, 154–162.
Lu, X., Wu, Y., Zhou, L., Tian, W., Su, G., Qiu, S., Zhang, H., 2014. Theoretical in-
vestigations on two-phase flow instability in parallel channels under axial non-uni-
form heating. Ann. Nucl. Energy 63, 75–82.
Fig 15. The comparison of critical inlet temperature (Tin-cr) and the corre- Aritomi, Masanori, Aoki, Shigebumi, Inoue, Akira, 2012. Instabilities in parallel channel
sponding Tin-max. of forced-convection boiling upflow system, (II). J. Nucl. Sci. Technol. 14, 88–96.
Matsuyama, K., Chikamatsu, K., Asada, H., 1995. General characteristics of two-phase
flow distribution in a multipass tube. IEEE Trans. Magnetics 24, 32–44.
Muñoz-Cobo, J.L., Podowski, M.Z., Chiva, S., 2002. Parallel channel instabilities in
References
boiling water reactor systems: boundary conditions for out of phase oscillations. Ann.
Nucl. Energy 29, 1891–1917.
Ablanque, N., Oliet, C., Rigola, J., Pérez-Segarra, C.D., Oliva, A., 2010. Two-phase flow Papini, D., Cammi, A., Colombo, M., Ricotti, M.E., 2012. Time-domain linear and non-
distribution in multiple parallel tubes. Int. J. Therm. Sci. 49, 909–921. linear studies on density wave oscillations. Chem. Eng. Sci. 81, 118–139.
Clausse, A., Lahey, R.T., Jr., 1990. An investigation of periodic and strange attractors in Papini, D., Colombo, M., Cammi, A., Ricotti, M.E., 2014. Experimental and theoretical
boiling flows using chaos theory. In: Proceedings of the Ninth International Heat studies on density wave instabilities in helically coiled tubes. Int. J. Heat Mass Transf.
Transfer Conference, Jerusalem, vol. 2, pp. 3–8. 68, 343–356.
Colombo, M., Cammi, A., Papini, D., Ricotti, M.E., 2012. RELAP5/MOD3.3 study on Su, Y., Feng, J., Zhao, H., Tian, W., Su, G., Qiu, S., 2013. Theoretical study on the flow
density wave instabilities in single channel and two parallel channels. Prog. Nucl. instability of supercritical water in the parallel channels. Prog. Nucl. Energy 68,
Energy 56, 15–23. 169–176.
Dutta, G., Zhang, C., Jiang, J., 2015. Analysis of parallel channel instabilities in the Tao, W.Q., 2010. Numerical Heat Transfer, second ed. Xi’an Jiaotong University, Xi’an.
CANDU supercritical water reactor. Ann. Nucl. Energy 83, 264–273. Xi, X., Xiao, Z., Yan, X., Li, Y., Huang, Y., 2014a. An experimental investigation of flow
Filonenko, G.K., 1954. Hydraulic resistance of pipelines. Therm. Eng. 40–44. instability between two heated parallel channels with supercritical water. Nucl. Eng.
Gao, F., Luo, Y.S., Chen, T.K., Lu, D.H., Li, H.X., 2005a. Comparative research on steam- Des. 278, 171–181.
water two-phase flow instabilities between internally ribbed tubes and smooth tubes. Xi, X., Xiao, Z., Yan, X., Xiong, T., Huang, Y., 2014b. Numerical simulation of the flow
Nucl. Power Eng. 26, 334–339. instability between two heated parallel channels with supercritical water. Ann. Nucl.
Gao, F., Chen, T.K., Luo, Y.S., Fei, Y., Liu, W.M., 2005b. Experimental research on density Energy 64, 57–66.
wave oscillation of steam-water two-phase flow in parallel inclined internally ribbed Xia, G., Peng, M., Guo, Y., 2012. Research of two-phase flow instability in parallel narrow
pipes. Nucl. Power Eng. 26, 11–14. multi-channel system. Ann. Nucl. Energy 48, 1–16.
Guo, Y., Huang, J., Xia, G., Zeng, H., 2010. Experiment investigation on two-phase flow Xia, G.L., Su, G.H., Peng, M.J., 2016. Analysis of flow distribution instability in parallel
instability in a parallel twin-channel system. Ann. Nucl. Energy 37, 1281–1289. thin rectangular multi-channel system. Nucl. Eng. Des. 305, 604–611.
Huang, F., Luo, Y., Wang, H., Chen, T., Wu, Y., 2011. Experimental investigation of two- Xiong, T., Yan, X., Xiao, Z., Li, Y., Huang, Y., Yu, J., 2012. Experimental study on flow

18
J. Liu et al. Nuclear Engineering and Design 342 (2019) 10–19

instability in parallel channels with supercritical water. Ann. Nucl. Energy 48, 60–67. phase flow instability in parallel narrow rectangular channels. Ann. Nucl. Energy 50,
Xiong, T., Yan, X., Huang, S., Yu, J., Huang, Y., 2013. Modeling and analysis of super- 103–110.
critical flow instability in parallel channels. Int. J. Heat Mass Transf. 57, 549–557. Yun, G., Qiu, S.Z., Su, G.H., Jia, D.N., 2008. Theoretical investigations on two-phase flow
Yadigaroglu, G., Bergles, A.E., 1972. Fundamental and higher-mode density-wave oscil- instability in parallel multichannel system. Ann. Nucl. Energy 35, 665–676.
lations in two-phase flow. J. Heat Transfer 94, 189–195. Zhang, Y.J., Su, G.H., Yang, X.B., Qiu, S.Z., 2009. Theoretical research on two-phase flow
Yuan, Z., Xiao, Y., Yanlin, W., Yanjun, L., Yanping, H., 2012. Experimental study of two instability in parallel channels. Nucl. Eng. Des. 239, 1294–1303.

19