The dynamic analysis of a uniformly heated channel at a supercritical pressure using a simple nonlinear model

Annals of Nuclear Energy 126 (2019) 95–109

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

The dynamic analysis of a uniformly heated channel at a supercritical pressure using a simple nonlinear model Jin Der Lee a,⇑, Shao Wen Chen b a

Nuclear Science and Technology Development Center, National Tsing Hua University, 101, Sect.2, Kuang Fu Road, Hsinchu 30013, Taiwan, ROC Institute of Nuclear Engineering and Science and Department of Engineering and System Science, National Tsing Hua University, 101, Sect.2, Kuang Fu Road, Hsinchu 30013, Taiwan, ROC b

a r t i c l e

i n f o

Article history: Received 1 March 2018 Received in revised form 31 October 2018 Accepted 2 November 2018

Keywords: Supercritical water Stability Nonlinear analysis Period-doubled bifurcation Chaos

a b s t r a c t This study develops a simple nonlinear dynamic model of a supercritical uniformly heated channel by dividing the channel into three regions coupled with linear profile approximations between flow enthalpy and flow density. Stability map, nonlinear characteristics and parametric effects of a uniformly heated channel with supercritical water are investigated. The results indicate that the nonlinear oscillations between inlet flow velocity and outlet flow velocity tend to present out-of-phase in the supercritical heated system. The parametric studies on system stability suggest that increasing inlet flow resistance or enlarging the channel diameter would stabilize the system, while increasing outlet flow resistance or lengthening the channel length would destabilize the system. In addition, complex nonlinear phenomena, i.e. supercritical Hopf bifurcations and period-doubled bifurcations, could exist in this supercritical uniformly heated channel system. In the deep unstable region of high pseudo-subcooling number, various periodic oscillations through a series of period-doubled bifurcations would eventually lead to chaotic oscillations. The appearance of period three (P-3) oscillation deduces that unimaginable periodic types of nonlinear oscillations could occur in the limited unstable space of this system. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction The advanced technology of a heated system with supercritical fluids has been applied to the current new type of thermal power plants and adopted in next-generation nuclear power plants. In the Generation-IV international advanced nuclear reactor development program, supercritical water nuclear reactor (SCWR) is selected as one of the promising types, which has attracted a lot of attentions, even in Japan after the Fukushima Daiichi accident (Miwa et al., 2018). In boiling water reactors (BWRs), flow boiling occurs inside the reactor core, but boiling never occurs in SCWRs due to a supercritical pressure (approximately 25 MPa) and the coolant flow remains in a single phase throughout the system. Fluid properties and transport phenomena are dramatically changed near a pseudo-critical point. Therefore, a heated system with supercritical fluids may be susceptible to different types of instability, particularly density-wave instability. The consideration of

⇑ Corresponding author at: Nuclear Science and Technology Development Center, National Tsing Hua University, 101, Sect.2, Kuang Fu Road, Hsinchu 30013, Taiwan, ROC. E-mail address: [email protected] (J.D. Lee). https://doi.org/10.1016/j.anucene.2018.11.004 0306-4549/Ó 2018 Elsevier Ltd. All rights reserved.

instability problems is crucial for the design, operation, and safety of a supercritical fluid system. Because the phase change process never takes place in the supercritical heated system, Ambrosini and Sharabi (2008), Ortega G’omez et al. (2008) and T’Joen et al. (2011) all suggested new definitions of dimensionless parameters, i.e. pseudosubcooling number (sub-pseudo-critical number) and pseudophase change number (apparent trans-pseudo-critical number), to illustrate the stability maps of the supercritical heated systems. Ambrosini and Sharabi (2008), Ambrosini (2009, 2011) and Ampomah-Amoako and Ambrosini (2013) conducted a series of studies on stability concerns of a forced convection heated channel with supercritical fluids. On the basis of fluid properties at a pseudo-critical point, they suggested the formulation of dimensionless density and dimensionless enthalpy. Therefore, the curves of dimensionless density versus dimensionless enthalpy among different supercritical fluids (H2O, CO2, NH3 and R23) were almost consistent in a reasonably wide range of system pressures. In addition, Ambrosini (2011) and Ampomah-Amoako and Ambrosini (2013) illustrated that stability boundaries among these supercritical fluids could meet quite well based on the dimensionless parameters, namely sub-pseudo-critical number and true trans-

96

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

Nomenclature AH Cp DH f Fr g h hPC hi h+ k L LH L+ M M+ NA NB NP
PCH

cross sectional area of heated channels (m2) constant pressure specific heat (Jkg
1K
1) diameter (m) friction factor Froude number, =u2i0/gLH gravity acceleration (ms
2) enthalpy (Jkg
1) fluid enthalpy at pseudo-critical point (Jkg
1) inlet fluid enthalpy (Jkg
1) non-dimensional fluid enthalpy,¼ ðh
hPC ÞbPC =C P;PC loss coefficient length (m) channel length (m) non-dimensional length, =L/ LH mass (kg) non-dimensional mass,¼ M=qPC AH LH number of nodes in the region 1 number of nodes in the region 2 pseudo-phase change number (apparent trans-pseudoPC critical number),¼ q AQH ui0 CbP;PC PC

NTPC

bPC Q true trans-pseudo-critical number,¼ NP
PCH qþ ¼ q AH ui0 C P;PC

NP
SUB

pseudo-subcooling number (sub-pseudo-critical numPC ber),¼ CbP;PC ðhPC
hi Þ system pressure (Pa) heating power (W) linear power (Wm
1) heat flux (Wm
2) Reynolds number, =uD/m temperature (K) pseudo-critical temperature (K) time (s) time scale, = LH/ui0 non-dimensional time, =t/tref velocity (ms
1)

P Q q0 00 q Re T TPC t tref t+ u

i

i

pseudo-critical number, for the stability map of the supercritical heated channel system. In previous studies, density-wave oscillations (DWOs) have been identified in some forced convection heated channel systems with supercritical fluids, while the inconsistent results on other types of instabilities were reported in the relevant literatures. Ambrosini and Sharabi (2008), Ambrosini (2009, 2011) and Ampomah-Amoako and Ambrosini (2013) demonstrated that Ledinegg instability can be predicted in their heated systems with different supercritical fluids through numerical studies, whereas Zhang et al. (2015) did not predict the Ledinegg instability in the similar analytical system. Ortega G’omez et al. (2008) inferred that a supercritical heated system might not experience the instabilities of Ledinegg excursion and pressure drop oscillations (PDOs). Swapnalee et al. (2012) suggested that Ledinegg instability could appear in a heated system with supercritical water; however, it was not found for a heated system with supercritical CO2. Zhao (2005) deduced that there would be no Ledinegg type static instability in the U.S. reference SCWR design at steady-state operation conditions. In addition, Yu et al. (2011) studied Ledinegg instability in a natural circulation system at supercritical pressure by using a numerical analysis code verified through experimental results, and demonstrated that there would be no Ledinegg instabilities occurring at supercritical pressure in the loop. Concerning DWOs in the supercritical heated system, parametric effects on system stability were found to be similar to those in

ui0 u+ W z z+

steady state inlet velocity (ms
1) non-dimensional velocity,¼ u=ui0 mass flow rate (kgs
1) axial coordinate (m) non-dimensional axial coordinate, =z/LH

Greek symbols b isobaric thermal expansion coefficient (K
1) m kinematic viscosity (m2/s) DP pressure drop (Pa) DP+ non-dimensional pressure drop,¼ DP=qPC u2i0 q density (kgm
3) + q non-dimensional density,¼ q=qPC qPC fluid density at pseudo-critical point (kgm
3) K friction number,¼ fLH =2DH k the boundary of region (m) k+ non-dimensional boundary of region, =k/LH Subscripts a accelerational pressure drop A the point of the boundary of region 1 B the point of the boundary of region 2 ch channel e exit of heated channel f frictional pressure drop g gravitational pressure drop H heated channel I inertial pressure drop i inlet of heated channel j j-th node in the region 2 n n-th node in the region 1 P constant pressure PC pseudo-critical 0 steady state

two-phase flow systems. Ambrosini and Sharabi (2008), Ortega G’omez et al. (2008) and Shankar et al. (2018) all indicated that the inlet flow resistance had a stable effect, while outlet flow resistance generated an unstable influence, on the stability of the forced convection supercritical heated system. Moreover, Shankar et al. (2018) also revealed that the length effect would cause the system more unstable, whereas the diameter effect would drive the system more stable, in the analytical cases of a supercritical nuclear heating system. Linear (frequency domain) and nonlinear (time domain) analysis methods can be employed to conduct theoretical stability analyses. The linear analysis is prevalently adopted to investigate the stability of a heated system with supercritical fluids because of its simplicity (Ambrosini and Sharabi, 2008, Ambrosini, 2009, 2011; Ortega G’omez et al., 2008; T’Joen et al., 2011; Shankar et al., 2018; Zhao, 2005, Zhao et al., 2007. Zhao et al. (2007) proposed a three-region methodology, which included three regions namely heavy fluid (H-F) region, heavy and light fluid mixture (H-L-M) region, and light fluid (L-F) region, to treat supercritical flows. Here, two boundary points, i.e. the pseudo-saturated water point between the H-F and H-L-M regions and the pseudosaturated steam point between the H-L-M and L-F regions, divide supercritical water flow into three regions for a linear stability analysis. Most previously published studies regarding the stability and nonlinear phenomena of the supercritical heated system deter-

97

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

mined their stability maps through linear methodologies and then predicted nonlinear characteristics by using the large system codes and CFD tools, such as COMSOL (Ortega G’omez et al., 2008), RELAP5 (Ambrosini and Sharabi, 2008; Ambrosini, 2009, 2011), and STAR-CCM+ (Ampomah-Amoako and Ambrosini, 2013). Ortega G’omez et al. (2008) suggested that with increasing operating power, Hopf bifurcation and period-doubled bifurcation might occur in their single supercritical uniformly heated channel system. These nonlinear phenomena are carried out through the nonlinear analysis in the present study. The nonlinear phenomena of supercritical fluid heated systems are ambiguous and must be clarified, while large system codes are generally not suitable for the detailed analysis of these phenomena because of their complexity and time-consumption. Therefore, in this study, the methodology employed by Ambrosini and Sharabi (2008), the three-region methodology adopted by Zhao et al. (2007), and the linear profile approximation between flow enthalpy and flow density proposed by Swapnalee et al. (2012) are integrated to develop a simple nonlinear dynamic model of a uniformly heated channel at a supercritical pressure. The stability map, nonlinear dynamics and parametric effects of a supercritical uniformly heated channel can be completely analyzed using the present nonlinear dynamic model. This study primarily aims to survey the nonlinear nature and bifurcation phenomena occurring in the density-wave unstable region of the uniformly heated channel system with supercritical water. 2. The models 2.1. Three region model with linear profile approximations Fig. 1 illustrates the modeling structure of a uniformly heated channel at a supercritical pressure, which was the similar system investigated by Ambrosini and Sharabi (2008). The three-region methodology suggested by Zhao et al. (2007) is employed to ana-

lyze the heated channel with supercritical water. Thus, the channel is discretized into three regions (Fig. 1): (i) the heavy fluid region (region 1) that is like an incompressible liquid, (ii) the heavy and light fluid mixture region (region 2), which is similar to a homogeneous equilibrium two-phase mixture, and (iii) the light fluid region (region 3) that behaves like an ideal gas or superheated steam. The boundary point between regions 1 and 2, denoted as ‘‘A”, represents a fluid temperature of approximately 350 °C þ (hA ¼
0:892). The boundary point between regions 2 and 3, denoted as ‘‘B”, can be determined using the following relation (Zhao et al., 2007):

T B ¼ 572:545 þ ð

P
13:9189
3

1:0193  10

0:5

Þ ðKÞ

ð1Þ

where the system pressure (P) is in the unit of MPa. Therefore, the þ value of TB is nearly 676.81 K  404 °C (hB ¼ 0:8) when the system pressure is 25 MPa. The following assumptions are made to simplify the problem: (1) The system pressure is kept at 25 MPa under both steady and dynamic conditions, as the pressure drop through the channel is much smaller than the system pressure, (2) The heat flux is assumed to be uniform in the axial direction, (3) The channel has a constant inlet enthalpy under both steady and dynamic conditions, (4) The channel is supposed to have a fixed pressure drop, which is equivalent to that of the corresponding steady state in initial operation condition. As defined in the Nomenclature section, state variables can be non-dimensionalized based on thermal physical properties at pseudo-critical points listed in Table 1 for different system pressures used in this study. On the basis of the assumptions and the methodology proposed by Ambrosini and Sharabi (2008), onedimensional dimensionless conservation equations for a supercritical uniformly heated channel can be expressed as follows:

@ qþ @ qþ uþ þ ¼0 @t þ @zþ

qþ

þ

ð2Þ þ

@h @h þ qþ uþ þ ¼ NP
PCH @t þ @z

ð3Þ

N X @ qþ uþ @ qþ uþ2 qþ uþ2 qþ @Pþ 2 þ ¼
Kqþ uþ
km dðzþ
zþm Þ

@tþ @zþ 2 Fr @zþ m¼1

ð4Þ where NP-PCH is pseudo phase-change number (apparent transpseudo-critical number) originally given by Ambrosini and Sharabi (2008), as defined in Nomenclature. Ambrosini and Sharabi (2008) revealed that estimating the dimensionless density independent of pressure and being the function of enthalpy is appropriate at a supercritical pressure. Thus,

qþ ¼ qþ ðhþ Þ

ð5Þ

Through Eqs. (2) and (3) together with Eq. (5), it can get the velocity gradient along the axial direction as follows:

@uþ dqþ NP
PCH ¼
þ @zþ qþ2 dh

Fig. 1. Schematic of a uniformly heated channel with the three-region methodology (Zhao et al., 2007) at a supercritical pressure.

ð6Þ

For a system pressure of 25 MPa, Fig. 2 illustrates the relation between dimensionless density (qþ ) and dimensionless enthalpy þ (h ) extracted from the REFPROP database of NIST (2010). On the basis of three-region model (Fig. 1) and regarding the continuities of flow properties at the boundary points between regions, i.e.

98

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

Table 1 The thermal physical properties at pseudo critical points corresponding to different supercritical pressures (NIST, 2010). P(MPa)

T PC (K)

qPC (kg/m3)

hPC (kJ/kg)

bPC (1/K)

C P;PC (kJ/kgK)

23 23.5 25

650.62 652.5 658.04

316.60 315.95 317.23

2115.5 2126.9 2152.2

0.51395 0.30467 0.12852

286.05 172.54 76.44

Assuming the linear enthalpy profile and adopting the linear þ relation between qþ and h by using Eq. (7), integrating the momentum conservation equation, Eq. (4), from the inlet of the channel to the boundary of region 1 (kþ A ) can get the dynamic pres: sure drop of region 1, DPþ H;A

DPþH;A ¼ DPþI;A þ DPþa;A þ DP þf ;A þ DP þg;A

ð12Þ

where DPþ I;A is the inertial pressure drop of region 1 as:

ðqþA þ qþi Þ þ duþi kA þ þ DPþI;A0 2 dt

DPþI;A ¼

ð13Þ

ðqþ þ qþi Þ þ NP
PCH kþA dkþ ui
b1 DPþI;A0 ¼ ½ A
qþA uþA  þA þ 2 qi dt

ð14Þ

DP þa;A is the accelerational pressure drop of region 1 as:

DPþa;A ¼ qþA uþA
qþi uþi 2

DP þf ;A þ

qþ ðhþA ð
ÞÞ ¼ qþ ðhþA ðþÞÞ at point ‘‘A” and qþ ðhþB ð
ÞÞ ¼ qþ ðhþB ðþÞÞ at point ‘‘B”, the linear profile approximation suggested by Swapnalee et al. (2012) is employed to approximate the relation between þ dimensionless density (qþ ) and dimensionless enthalpy (h ) in each region, as shown in Fig. 2. Thus, þ

q ¼ a1 þ b1 h

þ

þ

; h 

þ hA

for region 1

ð7Þ

qþ ¼ a2 þ b2 hþ ; hþA  hþ  hþB for region 2

ð8Þ

qþ ¼ a3 þ b3 hþ ; hþ  hþB for region 3

ð9Þ

where the constants in each region at 25 MPa are as follows: þ þ hA ¼
0:892, hB ¼ 0:8, a1 = 1.53924, b1 =
0.4861, a2 = 1.157285, b2 =
0.9143, a3 = 0.477717, and b3 =
0.06484. þ

2.2. Region 1 (heavy fluid region, h+hA ) Region 1 in Fig. 1 is divided into NA nodes. For the uniformly heated channel, one can assume the linear enthalpy profile between two neighboring nodes, employ the linear relation þ between qþ and h by using Eq. (7), and integrate the energy conservation equation, Eq. (3), to obtain the dynamics of the boundary of region 1 (kþ A ) and of other nodes in region 1: þ

þ

qþn hn uþn
qþn
1 hn
1 uþn
1
NP
PCH ðLþn
Lþn
1 Þ dLþn NA þ þ þ ¼ dt ð2qþn =3 þ qþn
1 =3
a1 =2Þ ðhA
hi Þ


ð2qþn
1 =3 þ qþn =3
a1 =2Þ dLþn
1 ; n ¼ 1; 2; . . . ; NA ð2qþn =3 þ qþn
1 =3
a1 =2Þ dtþ

NP
PCH kþA 1 1 þ þ ð þ
þ Þ; n ¼ 1; 2; . . . ; N A ðhA
hi Þ qþ ðhn Þ qi

 þ ðqA þ qþi Þ þ2 b1 uþi NP
PCH kþA ui
2 qþi !2 " #9 ‘nðqþA =qþi Þ 2 ðqþA þ qþi Þ = 1 NP
PCH kþA 2 þ k qþ uþ þ
2 þ þ ; 2 i i i ðqþA
qþi Þ qþi hA
hi 2qþi

DPþf ;A ¼ KA kþA þ

ð16Þ DP þg;A

DPþg;A ¼

is the gravitational pressure drop of region 1 as:

ðqþA þ qþi Þ kþA 2 Fr

ð17Þ þ

þ

2.3. Region 2 (heavy and light fluid mixture region, hA  h+ hB ) Region 2 in Fig. 1 is divided into NB nodes. Under a uniform heating condition, assuming the linear enthalpy profile between two neighboring nodes and employing the linear relation between qþ and hþ by Eq. (8), it can integrate the energy conservation equation, Eq. (3), to yield the dynamics of the boundary (kþ B ) of region 2 and of other nodes in the region 2: þ

þ

qþj hj uþj
qþj
1 hj
1 uþj
1
NP
PCH ðLþj
Lþj
1 Þ dLþj NB þ þ þ ¼ ð2qþj =3 þ qþj
1 =3
a2 =2Þ dt ðhB
hA Þ


ð2qþj
1 =3 þ qþj =3
a2 =2Þ dLþj
1 ; j ¼ 1; 2; . . . ; NB ð2qþj =3 þ qþj
1 =3
a2 =2Þ dt þ

ð18Þ

where non-dimensional node velocity (uþ j ) in region 2 can be obtained by integrating Eq. (6) from the boundary of region 1 (kþ A) to each node:

ð10Þ

where non-dimensional node velocity (uþ n ) in region 1 can be obtained by integrating Eq. (6) from the inlet of the channel to each node, which is as follows:

uþn ¼ uþi þ

ð15Þ

is the frictional pressure drop of region 1 as:

þ

Fig. 2. The relation between q and h extracted from the REFPROP database of NIST (2010), and the three-region model with linear profile approximations in this study.

2

ð11Þ

uþj ¼ uþA þ

NP
PCH ðkþB
kþA Þ þ þ ðhB
hA Þ



 1 1 ; j ¼ 1; 2; . . . ; NB
qþ ðhþj Þ qþA

ð19Þ

Assuming the linear enthalpy profile and taking the linear relaþ tion between qþ and h by Eq. (8), the dynamic pressure drop of þ region 2, DP H;B , can be got by integrating the momentum conservation equation, Eq. (4), from the boundary of region 1 (kþ A ) to the boundary of region 2 (kþ B ):

99

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

DPþH;B ¼ DPþI;B þ DPþa;B þ DPþf ;B þ DPþg;B

ð20Þ

ð1
kþB Þðqþe þ qþB Þ NP
PCH b2 b1 dkþ ð þ
þ Þ þA þ 2 qA qB qi dt  b3 NP
PCH ð1
kþB Þ ðqþB
qþe Þ þ uB þ þ 2 qþB  ð1
kþB Þðqþe þ qþB Þ b2 NP
PCH dkþB ð1
kþB ÞuþB dqþe þ
þ þ 2 qB qA dtþ 2 dt þ

DPþI;C0 ¼

where DPþ I;B is the inertial pressure drop of region 2 as:

duþ ðqþA þ qþB Þ þ ðkB
kþA Þ þi þ DPþI;B0 dt 2

DPþI;B ¼ DPþI;B0

ð21Þ

þ dkþA b1 ðqþA þ qþB Þ dkþ þ þ dkB ¼q NP
PCH ðkþB
kþA Þ þA þ
qB uB þ
þ þ dt dt 2 qi qA dt  þ   þ ðqA þ qþB Þ þ b2 NP
PCH ðkþB
kþA Þ dkB dkþA ð22Þ þ
uA
qþA dt þ dtþ 2 þ þ A uA

DPþa;B ¼ qþB uB
qþA uA þ2

DP þf ;B

DP þa;C is the accelerational pressure drop of region 3 as:

DPþa;C ¼ qþe uþe
qþB uþB 2

DP þf ;C

DP þa;B is the accelerational pressure drop of region 2 as: þ2

ð23Þ

is the frictional pressure drop of region 2 as:

 þ ðqA þ qþB Þ þ2 b2 uþA NP
PCH ðkþB
kþA Þ uA
qþA 2     2 ‘nðqþB =qþA Þ ðqþA þ qþB Þ NP
PCH ðkþB
kþA Þ 2 ð24Þ þ
þ 2 þ þ qþB
qþA qþA hB
hA 2qþA

DPþg;B ¼

ðqþA þ qþB Þ ðkþB
kþA Þ Fr 2

ð25Þ þ

2.4. Region 3 (light fluid region, h+hB ) Integrating the continuity equation, Eq. (2), from the inlet to the outlet of the supercritical heated channel can result in the dynamic equation of channel mass:

dM þch þ dt

¼q

þ þ i ui


q

þ þ e ue

ð26Þ

where the channel mass, Mþ ch , is the summation of the region mass, which can be expressed as follows:

M þch

ðqþ þ qþi Þ ðqþ þ qþB Þ ðqþ þ qþB Þ ¼ kþA A þ ð1
kþB Þ e þ ðkþB
kþA Þ A 2 2 2 ð27Þ uþ e,

The outlet velocity of the supercritical heated channel, can be obtained by integrating Eq. (6) from the boundary of region 2 þ (kþ B ) to the outlet of the channel (LH ¼ 1):

uþe ¼ uþB þ

NP
PCH ð1
kþB Þ 1 1 ð þ
þÞ þ þ q q ðhe
hB Þ e B

ð28Þ

Substituting Eq. (27) into Eq. (26) can lead to the dynamic equation of outlet density, qþ e: dkþ

dkþ

2ðqþi uþi
qþe uþe Þ þ ðqþB
qþi Þ dtþA þ ðqþe
qþA Þ dtþB dqþe þ ¼ dt ð1
kþB Þ

ð29Þ

Assuming the linear enthalpy profile and employing the linear þ relation between qþ and h by Eq. (9), integrating the momentum conservation equation, Eq. (4), from the boundary of region 2 (kþ B) to the outlet of the supercritical heated channel (Lþ H ¼ 1) can obtain the dynamic pressure drop of region 3, DPþ H;C :

DPþH;C

¼

where

DPþI;C

DPþI;C

DPþI;C ¼

þ

DPþa;C

þ

DPþf ;C

þ

DPþg;C

ð30Þ

2

ð33Þ

is the frictional pressure drop of region 3 as:

 þ ðqe þ qþB Þ þ2 b3 uþB NP
PCH ð1
kþB Þ : uB
2 qþB 2    NP
PCH ð1
kþB Þ ‘nðqþe =qþB Þ 2 ðqþe þ qþB Þ
þ þ 2 þ þ þ þ qþe
qB qB he
hB 2qþB

DPþf ;C ¼ KC ð1
kþB Þ

1 2 þ ke qþe uþe 2

DPþf ;B ¼ KB ðkþB
kþA Þ

DP þg;B is the gravitational pressure drop of region 2 as:

ð32Þ

ð34Þ

DP þg;C is the gravitational pressure drop of region 3 as:

DPþg;C ¼

ðqþe þ qþB Þ ð1
kþB Þ 2 Fr

ð35Þ

2.5. Inlet flow velocity dynamics of supercritical uniformly heated channel The supercritical heated channel is supposed to have a fixed dynamic pressure drop, which is equivalent to that of the corresponding steady state in an initial operation condition (DP þ ch;0 ). Thus,

DPþch ¼ DPþH;A þ DPþH;B þ DPþH;C ¼ DPþch;0

ð36Þ

Finally, the inlet flow velocity dynamics of supercritical uniformly heated channel, uþ i , can be obtained by substituting Eqs. (12), (20), and (30) into Eq. (36):

duþi ¼ ðDPþch;0
DPþI;A0
DP þa;A
DP þf ;A
DPþg;A
DPþI;B0
DPþa;B dtþ
DP þf ;B
DP þg;B
DPþI;C0
DPþa;C
DPþf ;C
DP þg;C Þ=½kþA ðqþA þ qi þ Þ=2 þ ðkþB
kþA ÞðqþB þ qAþ Þ=2 þ ð1
kþB ÞðqþB þ qeþ Þ=2

ð37Þ

3. Solution method The transient responses of this supercritical uniformly heated channel system following a perturbation in a state variable, i.e. inlet flow velocity (uþ i ), at a given initial steady state can be determined by solving the set of nonlinear, ordinary differential equations. This numerical analysis can be completed using the subroutine SDRIV2 of Kahaner et al. (1989), which employs the Gear multi-value method. The non-dimensional time step is ensured through numerical tests and the value used for the results in this study is set to Dt+=10
4. The requested relative accuracy in all solution components is set to 10
10. 4. Results and discussion 4.1. The validation of the present model

is the inertial pressure drop of region 3 as:

duþ ðqþe þ qþB Þ ð1
kþB Þ þi þ DPþI;C0 2 dt

ð31Þ

Table 2 lists the geometrical and operating conditions of a typical supercritical heated channel in SCWR used in the present analysis that was introduced in Ambrosini and Sharabi (2008).

100

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109 Table 2 The condition of a typical supercritical heated channel in SCWR adopted in the present analysis that was given in Ambrosini and Sharabi (2008). Parameter

Value

P LH DH ki ke Wi AH f

25 MPa 4.2672 m 3.4 mm 20 1 0.055 kg/s 5.4972  10
5 m2 0.0352

Pseudo-subcooling number (N P
SUB ) and true trans-pseudo-critical number (N TPC ), which were originally proposed by Ambrosini and Sharabi (2008), can be applied to present the stability map of the heated system with supercritical fluids. These two parameters are defined as follows:

NTPC ¼

Q

bPC

qi AH ui0 C P;PC

; NP
SUB ¼

bPC ðhPC
hi Þ C P;PC

ð38Þ

On the basis of the parameter values presented in Table 2, the present model can determine the stability map of the supercritical heated channel system through nonlinear analysis methodology. By setting NA = NB = 3, the examples illustrated in Fig. 3 have a fixed pseudo-subcooling number of N P
SUB = 2.468. The corresponding state is stable if the magnitude of oscillation dampens after perturbation, as shown in the example (Fig. 3(a)) of a stable state with

Fig. 3. The transient response of the corresponding state at a fixed pseudo-subcooling number of N P
SUB ¼ 2:468 following a perturbation, (a) stable state, (b) boundary state and (c) unstable state.

101

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

N P
SUB = 2.468 and N TPC = 4.482, whereas it is considered unstable if the magnitude of oscillation grows continuously, which is revealed in Fig. 3(c) for an unstable state with N P
SUB = 2.468 and N TPC = 4.651. Moreover, the results in Fig. 3(c) demonstrate that the unstable oscillation finally goes to a limit cycle; that implies a kind of supercritical Hopf bifurcation appearing in the unstable region of the present system. This is similar to that reported by Ortega G’omez et al. (2008). In addition, the boundary state generally presents a limit cycle oscillation with the lowest power one, as illustrated in Fig. 3(b) with N P
SUB = 2.468 and N TPC = 4.585. Fig. 4 evaluates the effect of node numbers on the stability boundary of the uniformly heated channel with coolant water at a supercritical pressure of 25 MPa. The present three-region model is valid for supercritical water heating systems with an inlet flow temperature ranging from 0 °C to 350 °C. If the inlet flow temperature is beyond 350 °C, the supercritical water heating system must be analyzed using the two-region model (region 2 plus region 3). If the inlet flow temperature is lower than 0 °C, the system may be in a frozen state, which represents an unrealistic condition for a supercritical water heating system. The results in Fig. 4 indicate that stability boundaries predicted using the present nonlinear dynamic model change slightly when the node number is greater than three for both of regions 1 (NA) and 2 (NB). Therefore, a set of node numbers with NA = NB = 3 is selected in the present numerical analysis for modeling accuracy and computational simplicity. Fig. 4 reveals that the non-monotonic effects of inlet temperature on the system stability of the uniformly heated channel with

supercritical water. A particular effect of inlet temperature on the system stability is presented in a relatively low pseudosubcooling number region. The result in Fig. 4 with NA = NB = 3 illustrates that the effect of inlet temperature firstly makes the system more stable as increasing the inlet temperature, i.e. the decrease in pseudo-subcooling number, when the system is with a pseudo-subcooling number (NP-SUB) less than 1.87. However, if the inlet temperature is continuously increased such that the value of NP-SUB is less than 1.4, the nonlinear effect of inlet temperature would generate an unstable influence on the system stability. This may be resulted from the contribution of the light fluid region, which is suggested to have a destabilizing effect as in the later discussion, and the limitations of the present model with three region linear approximations that needs a further discussion. The light fluid region (region 3) would be expanded more and more if the inlet temperature is consecutively increased, and thus drive the system more unstable in the small pseudo-subcooling number region. In addition, the similar phenomenon, i.e. the inversion of the stability boundary at low subcooling numbers, was also predicted by Zhao et al. (2012) in two-phase flow systems. It may need more attention on the stability subject of the heated system with near boiling or near pseudo-critical inlet fluids. Considering the density-wave instability, the present nonlinear dynamic model is validated against the experimental data in Khabensky and Gerliga (1995). In the aforementioned study, six cases were experimentally evaluated in a long heated coil (a height of 12.6 m and a diameter of 0.01 m) by using coolant water at

Fig. 4. The effect of node numbers on the stability boundary of the uniformly heated channel at a supercritical pressure of 25 MPa.

Fig. 5. The validation of the present nonlinear model against the experimental data reported in Khabensky and Gerliga (1995).

Table 3 The validation results of the present model against experimental data reported in Khabensky and Gerliga (1995). Case

P (MPa)

hi (kJ/kg)

U (kg/m2s)

q00 (Exp.) (kW/m2)

NTPC (Exp.)

NTPC (Predict.)

Relative error (%)

1 2 3 4 5 6

23 23 23 23 23.5 23.5

1200 1200 1590 840 1000 1000

180 150 300 135 150 100

81.3 69.5 116 69.5 116 58

4.09 4.196 3.501 4.662 6.882 5.162

3.587 3.734 2.726 3.958 6.526 6.052

12.30 11.01 22.14 15.10 5.17 17.24

102

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

pressures of 23 and 23.5 MPa, as displayed in Table 3. Therefore, the present model is non-dimensionlized based on thermal physical properties at pseudo critical points, as listed in Table 1, corresponding to the system pressure of 23 and 23.5 MPa, respectively. By using the relation of Eq. (1), the calculated value of TB is approximately 666.93 K at a system pressure of 23 MPa, while it is nearly

669.50 K at a system pressure of 23.5 MPa. Moreover, the sets of coefficients used in linear approximations at supercritical pressures of 23 and 23.5 MPa can be determined through a similar procedure described in Section 2.1. Zhao (2005) suggested that the friction factor for a supercritical fluid in a smooth tube can be calculated through the relation proposed by Filonenko (1954):

Fig. 6. Detail time domain oscillations of state variables in the density-wave unstable state of N P
SUB ¼ 3:313 and N TPC ¼ 5:581, (a) inlet flow velocity (uþ i ), (b) the boundary of þ þ þ þ region 1 (kþ A ), (c) the boundary of region 2 (kB ), (d) outlet flow velocity (ue ), (e) the outlet flow enthalpy (he ) and (f) the outlet flow density (qe ).

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

f ¼ ½ð1:82log10 ðReÞ
1:64Þ


2

ð39Þ

Because the information regarding inlet and outlet restrictions was not available in Khabensky and Gerliga (1995), the effects of inlet and outlet restrictions are neglected during model validation. Fig. 5 and Table 3 compare the predicted stability boundary with the experimental data for the six cases based on the true transpseudo-critical number (N TPC ). The results indicate that the stability boundaries predicted using the present nonlinear model are within a relative error of 25%, which should be reasonable considering the unavailability of the detailed experimental information.

Fig. 7. The effect of inlet loss coefficient (ki) on the stability boundary of the supercritical heated channel.

Fig. 8. The effect of outlet loss coefficient (ke) on the stability boundary of the supercritical heated channel.

103

4.2. Parametric effects on system stability Fig. 6 illustrates the detailed time-domain oscillations of state variables in another density-wave unstable state with N P
SUB ¼ 3:313 and N TPC ¼ 5:581. If the inlet flow velocity (uþ i ) in Fig. 6(a) is increased, it would enlarge heavy fluid region (region 1) and thus lengthen the boundary of region 1 (kþ A ) with a time delay, as indicated in Fig. 6(b). With the increasing length of kþ A, the boundary of region 2 (kþ B ) would move downward to increase the value of kþ B with a time delay, as illustrated in Fig. 6(c). This

Fig. 9. The effect of channel length (LH) on the stability boundary of the supercritical heated channel.

Fig. 10. The effect of channel diameter (DH) on the stability boundary of the supercritical heated channel.

104

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

change implies that the light fluid region (region 3) would get shorter with increasing inlet flow velocity. Thus, the acceleration effects become weaker such that outlet flow velocity (uþ e ) is decreased with a certain time delay, as revealed in Fig. 6(d). It further represents that the nonlinear oscillations between inlet flow velocity and outlet flow velocity are nearly out-of-phase. This

phenomenon is similar to that occurring in the two-phase flow system. Moreover, the increase in inlet velocity indicates an increase in the total flow rate. At the same heating power, the difference in the flow enthalpy between the channel inlet and outlet would reduce to cause the decrease of flow enthalpy at the channel outlet þ (he ) in Fig. 6(e) with a time delay. Consequently, it further leads to

þ Fig. 11. The detail nonlinear oscillations with a periodic cycle of two (P-2) and the phase diagram onto the kþ A
ui plane for this supercritical uniformly heated channel system at NP-SUB = 3.313, NTPC = 5.842.

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

the increase in the fluid density at the channel outlet (qþ e ), as revealed in Fig. 6(f). And, vice versa. Fig. 7 illustrates the effect of inlet loss coefficient (ki ) on the stability boundary of the present system based on the stability map of

105

a supercritical uniformly heated channel in Fig. 4 with ki ¼ 20. The pressure drop of region 1 (heavy fluid region) is in phase with the inlet flow velocity. The results indicate that the increase in the inlet loss coefficient would enlarge the pressure drop of region 1, thus

þ Fig. 12. Various periodic types of inlet velocity oscillations and their corresponding phase diagrams onto the kþ A
qe plane appearing in this supercritical uniformly heated channel system, (a) periodic cycle of four (P-4), (b) P-4 phase diagram, (c) periodic cycle of eight (P-8), (d) P-8 phase diagram, (e) periodic cycle of six (P-6) and (f) P-6 phase diagram.

106

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

stabilizing the system. This result is consistent with the results of Ambrosini and Sharabi (2008), Ortega G’omez et al. (2008) and Shankar et al. (2018). Fig. 8 reveals the effect of outlet loss coefficient (ke ) on the stability boundary of the present system based on the stability map of a supercritical uniformly heated channel in Fig. 4 with ke ¼ 1. The

result in Fig. 6 indicates that outlet flow velocity is generally outof-phase with the inlet flow velocity in this supercritical heated channel system. This represents that the pressure drop of region 3 (light fluid region) is out-of-phase with the inlet velocity. The results in Fig. 8 show that the increase of outlet loss coefficient would enlarge the pressure drop of region 3, thus destabilizing

þ Fig. 13. The detail nonlinear oscillations with a periodic cycle of three (P-3) and the phase diagram onto the kþ A
ui plane for this supercritical uniformly heated channel system at NP-SUB = 3.313, NTPC = 5.848.

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

the system. Furthermore, this result is consistent with the results of Ambrosini and Sharabi (2008), Ortega G’omez et al. (2008) and Shankar et al. (2018). On the basis of the stability map of a supercritical uniformly heated channel in Fig. 4 with LH , Fig. 9 explores the effect of channel length on the stability boundary of the uniformly heated

107

channel with supercritical water. At the same heat flux, lengthening the channel length would extend the light fluid portion and thus increase pressure drop of region 3. Therefore, increasing the channel length would generally result in an unstable effect on the system, as shown in Fig. 9. This is consistent with the trend reported by Shankar et al. (2018).

þ Fig. 14. The detail nonlinear chaotic oscillations and the phase diagram onto the kþ A
ui plane for this supercritical uniformly heated channel system at NP-SUB = 3.313, NTPC = 5.846.

108

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

Fig. 10 presents the diameter effect on the stability of this system compared with the stability map of a supercritical uniformly heated channel with DH in Fig. 4. For the system with the same flow rate, the enlargement of channel diameter would result in the channel with a lower flow velocity. It implies a smaller pressure drop across the channel, particularly the part of frictional pressure drop. Therefore, in general, the enlargement of channel diameter would reduce the channel pressure drop to stabilize the system, as indicated in Fig. 10. Moreover, this result is consistent with the trend suggested by Shankar et al. (2018). 4.3. Nonlinear periodic oscillations and chaotic oscillations The results in Fig. 6 demonstrate the limit cycle oscillation, i.e. a supercritical Hopf bifurcation, occurring in the unstable region of the present system at N P
SUB = 3.313 and N TPC = 5.581. If the true trans-pseudo-critical number (N TPC ), i.e. the operating power, is continually increased, period-doubled bifurcations could take place in the deep unstable region of this system depending on the operating conditions. A series of interesting nonlinear bifurcation phenomena could be further identified in the unstable region at a fixed pseudo-subcooling number of N P
SUB = 3.313. Fig. 11 presents an example of period-doubled bifurcation exhibiting in an unstable state with N P
SUB = 3.313 and N TPC = 5.842 in this uniformly heated channel with supercritical water. The detailed time-domain oscillations among state variables indicate a nonlinear oscillation with a periodic cycle of two (P-2), which evolves from limit cycle oscillation through period-doubled bifurcation. þ The corresponding phase trajectory onto the kþ A
ui plane confirms that it is indeed a type of P-2 oscillation, as illustrated in Fig. 11(f). By increasing the value of true trans-pseudo-critical number (N TPC ) at a fixed pseudo-subcooling number of N P
SUB = 3.313, various periodic types of nonlinear oscillation could appear in the unstable region of this system. The inlet velocity oscillation in Fig. 12(a) and the corresponding phase diagram onto þ the kþ A
qe plane in Fig. 12(b) reveal the nonlinear oscillation with a periodic cycle of four (P-4) at N P
SUB = 3.313 and N TPC = 5.843. This oscillation implies the occurrence of another period-doubled bifurcation. If the true trans-pseudo-critical number is increased to N TPC = 5.844, inlet velocity oscillation in Fig. 12(c) and the correþ sponding phase diagram onto the kþ A
qe plane in Fig. 12(d) indicate the nonlinear oscillation with a periodic cycle of eight (P-8) through period-doubled bifurcation. If the true trans-pseudocritical number is further increased to N TPC = 5.845, another interesting type of nonlinear oscillation with a periodic cycle of six (P-6) could be recognized by the inlet velocity oscillation in þ Fig. 12(e) and the corresponding phase diagram onto the kþ A
qe plane in Fig. 12(f). At an unstable state of N P
SUB = 3.313 and N TPC = 5.848, a particular type of periodic oscillation with a cycle of three (P-3) is identified in the present supercritical uniformly heated channel system, as illustrated in Fig. 13. The detailed time-domain oscillations among state variables for this system reveal that the transient responses eventually evolve to P-3 oscillations. The corresponding þ phase diagram onto the kþ A
ui plane, as shown in Fig. 13(f), further confirms that it is indeed a type of P-3 oscillation. A wellknown paper on ‘‘Period Three Implies Chaos” (Li and Yorke, 1975) suggests that if a system has a periodic point of period three, then it must have periodic points of every period. Consequently, it could deduce that there are infinite periodic types of nonlinear oscillation existing in the limited unstable space of this uniformly heated channel with supercritical water. At an unstable state with N P
SUB = 3.313 and N TPC = 5.846, a distinct type of chaotic oscillations could appear in this uniformly heated channel with supercritical water, as indicated in Fig. 14.

The results illustrate that time evolutions among state variables oscillate chaotically. The corresponding phase diagram onto the þ kþ A
ui plane shows that it is a kind of strange attractor.

5. Conclusions This study develops a simple nonlinear dynamic model of a uniformly heated channel with supercritical water by dividing the channel into three regions coupled with linear profile approximations between flow density and flow enthalpy. For a representative uniformly heated channel with supercritical water introduced in Ambrosini and Sharabi (2008), the stability map, nonlinear characteristics and parametric effects of a supercritical uniformly heated channel are investigated. The conclusions can be summarized as follows: (1) This study indicates that the nonlinear oscillations between inlet flow velocity and outlet flow velocity tend to exhibit out-of-phase. This phenomenon is similar to that occurring in the two-phase flow system. It implies that increasing the pressure drop of heavy fluid region (region 1) would stabilize the system, whereas increasing the pressure drop of light fluid region (region 3) would destabilize the system. (2) The results of the parametric study on system stability suggest that increasing inlet flow resistance or enlarging the channel diameter would have a stable effect, while increasing outlet flow resistance or lengthening the channel length would generate an unstable influence, on the system stability. (3) The results demonstrate that complex nonlinear phenomena, i.e. supercritical Hopf bifurcations and period-doubled bifurcations, could exist in this uniformly heated channel system with supercritical water. In the deep unstable region of the high pseudo-subcooling number, various periodic oscillations could evolve through a series of perioddoubled bifurcations and subsequently lead to chaotic oscillations. A distinct type of period three (P-3) oscillation identified in the present system infers that unimaginable periodic types of nonlinear oscillations could exist in the limited unstable space of the present system. The present nonlinear dynamic model is in a dimensionless form. Under the thermal physical properties at the pseudo critical points corresponding to different supercritical fluids, it could be extended to investigate the stability problems of a uniformly heated channel with different supercritical fluids. The advanced works are undergoing to overcome the limitations of the present nonlinear model, including the partition method of three regions, the more accurate approximation between flow enthalpy and flow density, etc., and regarding the power shape and heated wall dynamics as well.

Acknowledgement This work was supported by the Ministry of Science and Technology, Taiwan of the Republic of China through contract MOST 106-2221-E-007-091-MY3.

Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2018.11.004.

J.D. Lee, S.W. Chen / Annals of Nuclear Energy 126 (2019) 95–109

References Ambrosini, W., Sharabi, M., 2008. Dimensionless parameters in stability analysis of heated channels with fluids at supercritical pressure. Nucl. Eng. Design 238, 1917–1929. Ambrosini, W., 2009. Discussion on the stability of heated channels with different fluids at supercritical pressures. Nucl. Eng. Design 239, 2952–2963. Ambrosini, W., 2011. Assessment of flow stability boundaries in a heated channel with different fluids at supercritical pressure. Ann. Nucl. Energy 38, 615–627. Ampomah-Amoako, E., Ambrosini, W., 2013. Developing a CFD methodology for the analysis of flow stability in heated channels with fluids at supercritical pressures. Ann. Nucl. Energy 54, 251–262. Filonenko, G.K., 1954. Hydraulic resistance of the pipelines. Therm. Eng. 4, 40–44. Kahaner, D., Moler, C., Nash, S., 1989. Numerical Methods and Software. Prentice Hall. Khabensky, V.B., Gerlig, a,V.A., 1995. Hydrodynamic flow instabilities. Power Equipment Components. CBS Publisher & Distributors, New Delhi, India. Li, T.Y., Yorke, J.A., 1975. Period three implies chaos. Am. Mathemat. Monthly 82, 985–992. Miwa, S., Yamamoto, Y., Chiba, G., 2018. Research activities on nuclear reactor physics and thermal-hydraulics in Japan after Fukushima-Daiichi accident. J. Nucl. Sci. Technol. 55, 575–598. NIST, 2010. REFPROP, NIST Standard Reference Database 23, Version 9.0. Ortega G’omez, T., Class, A., Lahey, R.T., Schulenberg, T., 2008. Stability analysis of a uniformly heated channel with supercritical water. Nucl. Eng. Design 238, 1930–1939.

109

Shankar, D., Pandey, M., Basu, D.N., 2018. Parametric effects on coupled neutronicthermohydraulic stability characteristics of supercritical water cooled reactor. Ann. Nucl. Energy 112, 120–131. Swapnalee, B.T., Vijayan, P.K., Sharma, M., Pilkhwal, D.S., 2012. Steady state flow and static instability of supercritical natural circulation loops. Nucl. Eng. Design 245, 99–112. T’Joen, C., Gilli, L., Rohde, M., 2011. Sensitivity analysis of numerically determined linear stability boundaries of a supercritical heated channel. Nucl. Eng. Design 241, 3879–3889. Yu, J., Che, S., Li, R., Qi, B., 2011. Analysis of Ledinegg flow instability in natural circulation at supercritical pressure. Progr. Nucl. Energy 53 (6), 775–779. Zhang, Y., Li, H., Li, L., Wang, T., Zhang, Q., Lei, X., 2015. A new model for studying the density wave instabilities of supercritical water flows in tubes. Appl. Therm. Eng. 75, 397–409. Zhao, J., 2005. Stability analysis of supercritical water cooled reactors. In: Department of Nuclear Science and Engineering, Doctor of Philosophy in Nuclear Science and Engineer. Massachusetts Institute of Technology, pp. 51– 97. Zhao, J., Saha, P., Kazimi, M.S., 2007. Hot channel stability of supercritical watercooled reactors-I: steady state and sliding pressure startup. Nucl. Technol. 158, 158–173. Zhao, J., Tso, C.P., Tseng, K.J., 2012. Nonhomogeneous-nonequilibrium two-phaseflow model for nuclear reactor single-channel stability analysis. Nucl. Technol. 180, 78–88.