Annals of Nuclear Energy 163 (2021) 108552

Contents lists available at ScienceDirect

Annals of Nuclear Energy

journal homepage: [Link]/locate/anucene

Stability analysis of a parallel channel-loop system with single-phase

natural circulation under asymmetric conditions
Enping Zhu a, Ting Wang a, Pengcheng Zhao a,b,⇑, Shengwen Deng a, Shenquan Zeng a, Jie Li a
a
School of Nuclear Science and Technology, University of South China, Hengyang, Hunan 421001, China
b
Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu, Sichuan 610213, China

a r t i c l e i n f o a b s t r a c t

Article history: Flow instability can cause mechanical oscillations in the equipment as well as periodic changes in the
Received 5 March 2021 thermal stress, both of which have a great impact on the safety of the reactor. In addition, flow instability
Received in revised form 4 July 2021 can interfere with the control system, making it difficult to capture stable parameters. In this study, to
Accepted 6 July 2021
examine the change in the stability of a parallel channel-loop system with natural circulation under
Available online 23 July 2021
asymmetric conditions, the single-loop natural circulation system is taken as the starting point, and
the control equations are obtained by using a dimensionless analysis method. The Fourier expansion of
Keywords:
displacement term in the control equations is carried out to obtain the Jacobian matrix representing
Single-phase natural circulation
Jacobian matrix
the single-loop natural circulation system. Based on this matrix, the Jacobian matrix model of parallel
Asymmetric condition channel-loop system is constructed. According to the model, the natural circulation stability of parallel
Parallel channel loop Stability boundary channel-loop under different load and resistance difference is analyzed, and the effects of geometric char-
acteristics, coolant and power on the stability boundary under asymmetric condition are examined. The
results show that there are two critical Reynolds numbers, and the system becomes unstable when the
load difference introduced by the left and right loops is greater than the critical Reynolds numbers. The
area of stable region can be increased by increasing the aspect ratio and heating zone length, and reduc-
ing the pipe diameter and cooling zone length. Furthermore, the stability boundary of natural circulation
is sensitive to the aspect ratio and the length of heating zone. In addition, the stability can be improved
within the allowable range of the natural circulation by choosing liquid metal coolant and increasing the
pressure drop, and this model cannot accurately describe the changes in the stability boundary at low
power. Overall, the above results can serve as useful reference for the design of the parallel channel-
loop systems and can be used to improve the stability of natural circulation. The stability boundary
can be used to preliminarily judge the stability distribution and the variation trend of the parallel
channel-loop system with a certain degree of reliability.
Ó 2021 Elsevier Ltd. All rights reserved.

1. Introduction of natural circulation is disturbed, which can lead to the instability

of the natural circulation system and a series of safety problems.
The natural circulation can effectively remove the heat gener- According to the fluid state, natural circulation flow instability
ated by a reactor without requiring any external power. It is useful can be divided into single-phase flow instability and two-phase
for improving the inherent safety of reactor, removing the decay flow instability. Two-phase flow instability can be divided into sta-
heat, and simplifying the system. The reactor operating in natural tic instability and dynamic instability (Hou et al., 2017). At present,
circulation is a typical multi-physical coupling heat transport sys- there are more researches on the instability of two-phase flow, but
tem, therefore its flow rate and temperature are not independent less research on the instability of single-phase flow. However,
variables but are related to parameters such as power, operating there are big differences between the two, and their control equa-
pressure, geometrical parameters of the loop, etc. When a small tions and flow patterns are not all the same. In addition, liquid
disturbance occurs in the power or loop resistance, the steady state metal has a high boiling point, so only single-phase flow instability
occurs in the primary loop of liquid metal reactors. Therefore, it is
urgent to carry out research on the instability of single-phase nat-
ural circulation flow.
⇑ Corresponding author.
E-mail address: zpc1030@[Link] (P. Zhao).

[Link]
0306-4549/Ó 2021 Elsevier Ltd. All rights reserved.
E. Zhu, T. Wang, P. Zhao et al. Annals of Nuclear Energy 163 (2021) 108552

To elucidate the significance of natural circulation in enhancing

the inherent safety performance of reactors, several theoretical and
experimental studies have focused on the instability of natural cir-
culation loop. Lee and Pan, 2005 proposed a non-linear numerical
model to study the stability and nonlinear dynamics of a
nuclear-thermal coupled two-phase natural circulation loop, and
they investigated the stability diagram, parametric effects, and
transient characteristics of this system. Satoh et al. (1998) used
the finite difference method to solve the one-dimensional (1D)
energy and momentum equations for a parallel channel-loop nat-
ural circulation system with one heating tube and two cooling
tubes and obtained the variation in the flow pattern as a function
of heating power under different flow differences. Demazière
et al. (2008) studied two-phase-flow oscillations at the natural-
circulation CIRCUS test facilit, and adopted the multifractal
approach to study characterizing the chaos. Research results show
that this chaos is deterministic, which is in good agreement with
the experimental results (Marcel et al., 2010). Nadella et al.
(2018) obtained a semi-analytical model for the stability analysis
of a single-phase natural circulation loop with conventional local
surface heating and cooling based on linear stability analysis, and
they found that the local pressure losses and the wall thermal iner-
tia had a considerable influence on the stability behavior of the
studied loop. Zhu et al. (2013) used a marine working condition
simulation test facility to investigate the characteristics of a paral-
lel channel-loop system with single-phase natural circulation
under inclined conditions through theoretical analysis and
experiments. Fig. 1. Simplified single-loop model of natural circulation.
By analyzing the research status, it can be inferred that nonlin-
ear methods are generally used to solve the governing equations of
4) The effect of the viscosity term is ignored in the energy
asymmetric multi-loop system, and the linear methods are mostly
equation.
suitable for analyzing 1D single-loop system. In the rectangular
multi-loop natural circulation system, each loop is asymmetric
Based on the above assumptions, the continuity equation,
due to the variations in the heat exchange capacity or the resis-
momentum equation, and energy equation of the 1D single-
tance of the flow path. In this case, it is obviously uneconomical
phase flow natural circulation loop are obtained as follows (Lu
to use a nonlinear method to calculate each possibility, so there
et al., 2014):
is an urgent need to develop more efficient and facile methods.
In this study, an effective linear analysis method is used to derive @ qu @ðquÞ
þ ¼0 ð1Þ
dimensionless governing equations for the natural circulation in a @t @v
single-loop system. Subsequently, the displacement term in these
equations is Fourier expanded, and the Jacobian matrix, which rep- @ðquÞ @ðqu2 Þ @P 1 f
resents the stability of the natural circulation in the single-loop þ ¼  qg  qu2 ð2Þ
@t @x @x 2D
system, is obtained and verified. Finally, the parallel channel-loop
Jacobian matrix model is constructed based on the single-loop @ðTÞ @ðTÞ k @2T 1 4q
þu ¼ þ aðxÞ  ð3Þ
Jacobian matrix, which is used to obtain the stability boundary of @t @x qC p @x2 qC p D
the rectangular loop under different load and resistance differ-
8
ences, and the influence of the loop geometry, coolant and power >
< 1; heating zone
characteristics on the stability boundary is analyzed.
aðxÞ ¼ 1; colding zone
>
:
0; adiabatic zone
2. Construction and verification of linear analysis method
Here, Eqs. (1), (2), and (3) are the continuity, momentum, and
2.1. Nondimensionalization of governing equations energy equations, respectively. q is the fluid density, kg/m3; t is
the time, s; u is the flow velocity, m/s; g is the gravitational accel-
To simplify the calculation, the single-loop system with single- eration, m/s2; P is the pressure, Pa; h is the convective heat transfer
phase natural circulation should be simplified, as shown in Fig. 1. coefficient, W/(m2K); D is the hydraulic diameter, m; f is the fric-
Here, the following hypotheses are considered: tion coefficient; T is the temperature, K; k is the thermal conductiv-
ity, W/(mK); Cp is the isobaric specific heat capacity, J/(kgK); q is
1) A 1D circuit model is considered, thus the radial distribution the heat flux density, W/m2; x is the fluid position in 1D coordinate
of the coolant temperature as well as the drop in local resis- system of single loop, m.
tance is ignored. If the fluid is regarded as an incompressible fluid, the flow rate
2) It is believed that the buoyancy force satisfies the Boussi- is only related to time, not position. According to the hypothesis 2,
nesq assumption, so all the density terms in the governing the following expression can be obtained by integrating Eq. (2)
equation are constant except for the volume force term. along the entire loop after substituting into the equation:
3) Except heating zone and cooling zone, the other zones are Z
du 1 f 2 gb
considered to be adiabatic. þ l ¼ ðT  T 0 Þdx ð4Þ
dt 2 D Ltot

2
E. Zhu, T. Wang, P. Zhao et al. Annals of Nuclear Energy 163 (2021) 108552

where bis thermal expansion coefficient, K1, and Ltot is the total @h @h @2h
þ Re ¼ Fo 2 þ hðsÞ ð9Þ
length of loop, m. To obtain the stability boundary of loop, the @r @s @s
equation should be dimensionless (Cammarata et al. 2003), i.e.
Z
dRe Lt p 2d Gr 0
@h @h @2h þ Re ¼ gðsÞ  h  ds ð10Þ
þ Re ¼ Fo 2 þ aðxÞ ð5Þ ds 4D 2
@s @s @s
Z 8
dRe Lt Gr0 >
>
> 1 ; LHtot1 < s < H1LþH h
þ f ðReÞRejRej ¼ h  ds ð6Þ >
>
tot
ds 4D 2 >
>
>
>
< - 1 ; 1=2 þ H2 < s < 1=2 þ H2 þHc
ah hðsÞ ¼ Ltot Ltot ð11Þ
Fo0 ¼ ð7Þ >
>
ðLtot =2Þ2 >
>
>
>
>
>
>
:
q2 gbD2 L2tot  q 0 ; otherwise
Gr0 ¼ ð8Þ
l3 c p A 8
>
< 1; 0 < s < H
Ltot
where h is the dimensionless temperature; s is the dimensionless
gðsÞ ¼ - 1; 0:5 < s < 0:5 þ H ð12Þ
time; Re is the Reynolds number; Fo is the Fourier number; s is >
:
Ltot

the dimensionless displacement; Gr0 is the modified Graschev num- 0; otherwise

ber; A is the fluid flow area, m2; l is the dynamic viscosity, Pa  s; a
Here, the non-dimensional abscissa s ¼ x=Ltot , and the value of p
is the thermal diffusivity, m2/s;
and d depends on the fluid flow state. For example, if it is a turbu-
lent fluid, p is 0.316 and d is 0.25; If it is laminar flow, p is 64, d is 1.
2.2. Construction of single-loop Jacobian matrix If Fourier series expansion is performed on the variable dis-
placement terms in Eqs. (9) and (10), then a finite system of ordi-
There are two typical methods for describing natural circulation nary differential equations can be obtained, i.e., the functions h(s),
flow: nonlinear method and linear method. In the non-linear g (s), and h(s,s)can be expressed as follows:
method, the conservation equations of the system are discretized
8 P
1
and solved numerically. This method is effective for analyzing >
>
> hðsÞ ¼ hk eikps
the transient response of natural circulation. However, when the >
>
>
>
k¼1
geometric parameters, resistance, power, and other factors change, < P1

it is obviously uneconomical and resource-consuming to use non- gðsÞ ¼ g k eikps ð13Þ

>
>
linear methods for calculating each steady state. In this case, the >
>
k¼1
>
> P1
linear method is preferred to make the governing equations >
: hðs; sÞ ¼ ðak ðsÞ þ ibk ðsÞÞeikps
k¼1
dimensionless so that their analytical solutions can be obtained
to better examine the various factors affecting natural circulation. Substituting Eq. (13) into Eqs. (9) and (10), the kth order equa-
There are two commonly used linear methods. The first one is tion can be obtained as follows:
the Nyquist stability criterion method. In this method, the roots (
dak ¼ 4p2  k Foak þ 2kpRebk þ 2k1 p sinðkrÞ
2
of the characteristic equations are derived to obtain the stability
and instability regions of the system according to the positive and ð14Þ
¼ 4p2  k Fobk þ 2kpReak þ 2k1 p ð1 þ cosðkrÞÞ
dbk 2
ds
negative changes of the roots. However, this method cannot obtain
the circuit temperature distribution, and it is extremely difficult to
solve the loop with more segments (Goyal et al., 2020). The second dRe L1 p 2d
¼ Re
ds 4D
linear method is to solve the Jacobian matrix of partial differential 0
Gr
equations (Cammarata et al. 2003). In this method, the eigenvalues þ fðg 1 þ g 2 þ . . . þ g k Þ½ða1 þ a2 þ ::: þ ak Þ þ iðb1 þ b2 þ ::: þ bk Þg
2
of the Jacobian matrix are obtained, and the stability region of the ð15Þ
system is established according to the variation in the real part of
the eigenvalues, i.e., when the real part of the eigenvalue is greater
pH
than 0, the system is in an unstable region. In this paper, the second r¼ ð16Þ
method is used to calculate the stability of natural circulation. HþL
To obtain the Jacobian matrix, the lower left corner of Fig. 1 is If the partial derivatives on the left side of Eqs. (14) and (15) are
considered as the starting point of the 1D coordinate system, and 0, then the steady state can be obtained. Consequently, the (2 k-
the clockwise flow direction is the positive direction. Therefore, 1)(2 k-1) Jacobian matrix of the partial differential equations
Eqs. (5) – (6) can be expressed as follows (Wu and Sienicki, 2003): about Re, ak , and bk can be derived.

2 3
Gr0 Gr 0 Gr 0 Gr 0
ðd  2Þ LDtot pRe1d p  sinðrÞ p  ð1  cosðrÞÞ    kp
 sinðkr Þ kp
 cosð1 þ krÞ
6 7
6 2pb1 4p2 Fo 2pRe  0 0 7
6 7
6 2pa 2pRe 4p2 Fo  7
6 1 0 0 7
A¼6
6 .. .. .. .. .. 7
7 ð17Þ
6. . .  . . 7
6 7
6 2ðk  1Þpb p2 Fo
   ½2ðk  1Þ 2ðk  1ÞpR
2 7
4 n 0 0 5
2ðk  1Þpan 2ðk  1ÞpR ½2ðk  1Þ p2 Fo
2
0 0  ð2k1Þð2k1Þ

3
E. Zhu, T. Wang, P. Zhao et al. Annals of Nuclear Energy 163 (2021) 108552

The Jacobian matrix is composed of the first k terms of the Four-

ier series, and the linear equations of a loop can be described by
Jacobian matrix. When k is large enough, the stability boundary
of the loop is almost constant. In addition, if k is greater than
100, a more accurate non-dimensional temperature distribution
can be calculated by Eq. (18).

hðs; sÞ ¼ 2a1 cosðpsÞ  2b1 sinðpsÞ þ 2a2 cosð2psÞ  :::

þ 2ak cosðkpsÞ  2bk sinðkpsÞ ð18Þ

2.3. Verification of the single-loop Jacobian matrix

In this section, to verify the reliability of the linear analysis

method, the stability boundary of loop is derived using the Jaco-
bian matrix based on the Argonne Lead Loop Facility (ALLF) loop
(Zhang et al., 2020), which is then compared with the already
known calculation results. Table 1 shows the main geometric
parameters and thermal hydraulic parameters of the ALLF loop.
The critical Jacobian matrix of ALLF is obtained by substituting
the main parameters in Table 1 into the Jacobian matrix and using
C++ programming (Nayak et al., 1995). This matrix is shown in
Fig. 2.
Fig. 2. Critical Jacobian matrix of ALLF loop.
It can be seen from Fig. 2 that the maximum real part of the
Jacobian matrix is close to 0, and the entire system lies in the
boundary between stability and instability, i.e., the system is in a
critical state. The dimensionless numbers Re and Gr’ can reflect
the characteristics of the buoyancy and viscous force, which can
be used to judge whether the loop is in a stable state. Fig. 3 shows
the Re-Gr’ stability diagram of the ALLF loop in the critical state.
In Fig. 3, the solid line is the stability boundary, and its left and
right regions are the unstable and stable regions. When the ther-
mal power is increased, the steady-state operating point moves
from right to left, and the loop gradually enters the unstable region
from the stable region, i.e., the instability of natural circulation
appears. Fig. 3 shows the Jacobian calculation results of this paper
and Lu et al., 2014, and Fig. 4 shows the Re-Gr’ stability diagram of
the ALLF loop, which is obtained by Wu and Sienicki, 2003 based
on the Nyquist stability criterion.
According to the assumption 1, the Jacobian matrix is obtained
based on a 1D model without considering the pressure drop in
local resistance, therefore K = 0, where K is the shape resistance
coefficient. If the local resistance is considered, the drop in local
resistance can be converted to that in frictional resistance as Fig. 3. Re-Gr 0 stability diagram of ALLF loop.
follows:
 
fL
fLt
þ K ¼ t;eff (19) the heat generated by the core is absorbed by the heat exchanger.
D D
Therefore, the obtained results are conservative.
where Lt;eff is the equivalent total length of loop, m.
Comparing Fig. 3 and in Fig. 4 (where K = 0), it can be found that
the obtained stability boundary of the ALLF loop is consistent with 3. Stability boundary analysis under asymmetric conditions
that reported by Zhu et al., 2013. Numerically, it is closer to the sta-
bility boundary of the ALLF loop obtained by Wu and Sienicki, 2003 3.1. Construction of parallel channel-loop Jacobian matrix
using the Nyquist stability criterion. The reason for this difference
is that the Jacobian matrix constructed in this paper does not con- The linear method can be used to obtain the stable and unstable
sider the insufficient heat exchange of the heat exchanger, i.e., all regions of the single-loop natural circulation system. However, it
cannot be directly applied to a two-dimensional (2D) model
Table 1
Main parameters of ALLF loop.
because the Jacobian matrix is only applicable to a 1D model, thus
some assumptions need to be made for the parallel channel-loop
Parameters Value Parameters Value
system. Assuming that the fluids in each loop are completely sep-
Height H 5m Power 0.206 MW arated in the parallel channel-loop system, then there is no mass
Width L 1.1 m Hot leg temperature 627 K exchange between the loops, and only the momentum and energy
Height of heating 1.52 m Cold leg temperature 542 K
zoneHh
is exchanged. Under this assumption, the parallel channel-loop
Height of cooling 1.52 m Temperature rise 85 K natural circulation system can be visualized as two parallel chan-
zoneHc nels, as shown in Fig. 5.
Pipe diameterD 0.063 m Flow rate of lead–bismuth 0.501 ms- Suppose that the temperature at the inlet and outlet of the core
eutectic (LBE) 1
and the flow rate remain constant under slight disturbances. Then,
4
E. Zhu, T. Wang, P. Zhao et al. Annals of Nuclear Energy 163 (2021) 108552

Fig. 6. Re-Gr’ stability diagram of the left and right loops under asymmetric load
difference.

Fig. 4. Re-Gr’ stability diagram of ALLF loop that is obtained based on Nyquist It can be seen from Eq. (20) that the heat transfer capacity of
stability criterion. one loop can decrease or increase when disturbance occurs, and
the heat transfer capacity of the other loop changes accordingly.
However, the total heat transfer capacity remains unchanged, so
the temperature variation of the two loops can be obtained. From
the variation of temperature distribution, the changes in Re1 of
the left loop and Re2 of the right loop can be obtained, and then
the Re-Gr’ variation curve of the two loops can be derived, as
shown in Fig. 6.
It is clear from Fig. 6 that if the loop is initially symmetrical,
there are two critical Reynolds numbers, so that when a small dis-
turbance is generated in a loop, both the loops are still in the stable
region. However, the loop changes from the stable region to the
unstable region when the disturbance is greater than the critical
value, and the greater the disturbance, the more it deviates from
the stability boundary.
According to the assumptions, the temperature and flow veloc-
ity of the core inlet and outlet of the loop do not change under
slight disturbances. Therefore, this model cannot accurately
Fig. 5. Analysis model for the Jacobian matrix of parallel channel-loop system. describe the unstable conditions under large disturbances. Regard-
ing the flow and heat transfer characteristics of the core, the sys-
the total heat generated by the core remains unchanged, and it is tem experiences periodic or non-periodic oscillations of flow and
completely absorbed by the left and right heat exchangers, i.e., temperature when it operates in an unstable region under sym-
metrical conditions. The changes in the thermal parameters, such
Q ¼ Q1 þ Q2 ð20Þ as flow and temperature, can lead to changes in the reactivity,
When Eq. (20) is satisfied, the Jacobian matrix of the left and which in turn affects the power produced by the core, and power
right loops can be constructed separately. Further, the influence changes affect the thermal parameters because there is feedback
of the left and right loops can be analyzed by this equation. Then, between the reactor power and thermal parameters.
the variations in the stability boundary under disturbance can be The large disturbance can be divided into two types according
analyzed to determine whether the loop fluid can finally remain to whether the loop can finally maintain a stable state. In the first
stable. In addition, the Jacobian matrix analysis model with parallel one, the loop can return to stability after a period of time when the
channel-loop is composed of two single-loop models in parallel, so disturbance is greater than the critical value. In addition, the
the result is calculated based on the single-loop model. The verifi- curves of the two loops should exhibit periodic or non-periodic
cation of the single-loop Jacobian matrix has been given, thereby oscillations with gradually decreasing amplitude above and below
the calculation results of parallel channel-loop model have a cer- the stability boundary (Zhou et al., 2017). In the second one, the
tain degree of reliability. loop cannot return to a stable state. In this case, the left and right
loop curves progressively deviate from the stability boundary, and
the entire system becomes increasingly unstable over time.
3.2. Analysis of stability boundary under asymmetric thermal load

In the parallel channel-loop system, due to the influence of var- 3.3. Analysis of stability boundary under asymmetric resistance
ious factors, such as changes in the steam load of each loop in the difference
secondary loop, changes in feed water flow, changes in feed water
temperature, or blockage of a loop flow path, the entire system is in If a slight resistance disturbance is generated in a loop, the pres-
a state of asymmetrical load difference or resistance difference, sure drop along the entire loop increases according to Eq. (21),
which leads to the instability of natural circulation. which may make the flow of natural circulation unstable.
5
E. Zhu, T. Wang, P. Zhao et al. Annals of Nuclear Energy 163 (2021) 108552

X qu2 f  l X qu2 f  Lt;eff Table 2 shows the geometric parameters of different natural cir-
DP 1 ¼ ð þ KÞ ¼ ð Þ ð21Þ
2 D 2 D culation loops, where Hl is the height of the natural circulation loop
and Loop1 is the control loop. The influence of loop geometry on
To analyze the variation in the stability boundary with the pres-
the stability boundary can be analyzed by changing the geometric
sure drop, four different Re-Gr’ working conditions are analyzed.
parameters of Loop1, for example, Hh and Hh , and keeping the total
The pressure drop DP 0 is used as the initial value, and the pressure
length of the loop unchanged. According to the nature of geometric
drop along the entire loop after the disturbance is generated is
variables, Table 2 can be divided into three groups for comparative
indicated by DP. The variations in the stability boundary of the loop
analysis:
are analyzed by using as a variable, and the results are shown in
Fig. 7.
1) First group: Loop1, Loop2, and Loop6 are used to study the
It is evident from Fig. 7 that as the pressure drop of the entire
influence of heating and cooling zone length on the stability.
loop increases, the stability boundary of the system gradually
2) Second group: Loop1, Loop3, and Loop4 are used to study
shifts upward and the stability region increases, i.e., the system
the influence of different aspect ratios on the stability.
becomes increasingly stable. However, from the perspective of nat-
3) Third group : Loop1 and Loop5 are used to study the effect of
ural circulation capability, the driving head of the natural circula-
pipe diameter on the stability of natural circulation.
tion continuously decreases due to the increase in the pressure
drop along the entire loop, which is harmful for the formation of
The area of the stability region under small disturbances is
natural circulation.
defined as the stability boundary margin. Fig. 8 shows the variation
By analyzing the different Re-Gr’ curves in Fig. 7, it can be seen
in the stability boundary margin of the loop under different geo-
that the larger the Reynolds number Re in the initial state of the
metric parameters when the reactor loop is in symmetrical
loop, the greater the changes in the loop stability boundary with
operation.
the pressure drop. Therefore, the stability boundary is more sensi-
It is clear from Fig. 8 that a change in the geometric parameters
tive to the change in pressure drop. In a large pressurized water
has a substantial impact on the stability boundary. It can be seen
reactor, even a small disturbance can cause a large change in Re
from Loop1, 2, and 6 that increasing the heating zone and reducing
because the temperature and mass flow in the primary circuit of
the length of the cooling zone can improve the loop’s ability to
the reactor are very large. Therefore, if the resistance distribution
maintain stability under small disturbances, and changing the
of each loop has a large deviation during the design of the reactor,
heating zone length has a stronger impact on the stability of the
a small disturbance can cause a significant change in the stability
loop. Further, it can be seen from Loop1, 3, and 4 that increasing
boundary, which may increase the system instability and affect
the aspect ratio can improve the loop’s ability to maintain stability
the safety of the reactor.
under small disturbances. The results of Loop1 and 5 indicate that
reducing the pipe diameter can improve the loop’s ability to main-
3.4. Sensitivity analysis of the system under asymmetric conditions
tain stability under small disturbances. The analysis shows that
geometric characteristics are the key factors that determine the
3.4.1. Influence of geometric parameters on the stability under
natural circulation driving force, as well as the frictional pressure
asymmetric conditions
drop and local pressure drop.
When the geometric parameters of the loop are different, the
Increasing the length of the heating zone, reducing the length of
transient characteristics and stability of its natural circulation
the cooling zone etc. can increase the driving force, and reducing
may be very different (Karami and Aghaie, 2017). Here, the influ-
the pipe diameter can reduce the friction pressure drop, both of
ence of the geometric parameters on the stability boundary is ana-
which can improve the stability of the loop.
lyzed using the controlled variable method.

Fig. 7. Re-Gr’ stability diagram of the left and right loops under asymmetric resistance difference.

6
E. Zhu, T. Wang, P. Zhao et al. Annals of Nuclear Energy 163 (2021) 108552

Table 2
Geometric parameters of natural circulation loop with different sizes (unit: m).

Loop Hl Hh Hc D L Lt
1 1.25 0.75 0.75 0.2 1.8 10
2 1.25 0.4 0.75 0.2 1.8 10
3 1.50 0.75 0.75 0.2 1.8 10
4 1.25 0.75 0.75 0.2 1.5 10
5 1.25 0.75 0.75 0.15 1.8 10
6 1.25 0.75 0.6 0.2 1.8 10

ent thermophysical properties. For example, the density of lead–

bismuth metal coolant is relatively large, and the specific heat
capacity of water is relatively large. Therefore, the natural circula-
tion stability may also be very different. Here, liquid metal lead–
bismuth eutectic (LBE), water, and liquid sodium are considered
as the research objects to analyze the influence of different cool-
ants on the stability. The Re-Gr’ stability diagrams of each coolant
under Loop1, 2, and 3 are shown in Fig. 10.
It can be seen from Fig. 10 that different coolants have a minor
effect on the stability boundary, and the stability boundary is
mainly related to the geometric structure and heating method.
However, when different coolants are in the same thermal–hy-
draulic environment, the physical parameters such as density, vis-
cosity, and flow rate are highly different under different coolants,
which can cause a significant change in the Re-Gr’ operating range
of the system. In order to analyze the variation of the stability
Fig. 8. Stability margin of circuits with different geometric parameters under
symmetrical condition.
operation interval under the same loop and power, and different
coolants, the relevant numerical simulation analysis is carried
out based on the FLUENT software. Fig. 11 shows the Re contour
of different coolants based on Loop1, where the power is 4 MW.
It is obvious from Fig. 11 that in the same loop, the power are
equal, but the Re range of natural circulation formed by different
coolants is not the same. The Re of liquid metal coolant is much lar-
ger than that of liquid water. Further, it can be seen from Fig. 10
that as Re increases, Gr’ also increases, and the stability boundary
moves upward. In addition, the slope also increases. This indicates
that if the natural circulation loop is formed in the high Re area, its
stability is higher than that in the low Re area. In the Loop1, the Re-
Gr’ stability boundaries of water, LBE, and sodium are almost the
same, which verifies that the coolant has a negligible effect on
the Re-Gr’ stability boundaries. However, the natural circulation
loop formed by liquid metal is in the high Re area. Therefore, com-
pared with water, liquid metal provides higher safety in the same
geometric loop and power level.

Fig. 9. Re-Gr’ stability boundary of loops with different geometric parameters.

Fig. 9 shows the changes in stability boundary of Loop1, 2, 3

under asymmetric thermal load disturbance.
It is clear from the Fig. 9 that increasing the aspect ratio can
increase the range of the stable region when a slight disturbance
is generated, while reducing the length of the heating zone reduces
the range of the stable region, and the stable region is highly sen-
sitive to both the factors. If there is a large deviation in the design
of each loop of the reactor, the stability boundary of each loop
deviates greatly, which can increase the instability of the entire
system.

3.4.2. Influence of coolant on the stability under asymmetric

conditions
Even in the same loop, the heat transfer and flow of different
coolants are quite different because different coolants have differ- Fig. 10. Re-Gr’ stability diagram under different coolants.

7
E. Zhu, T. Wang, P. Zhao et al. Annals of Nuclear Energy 163 (2021) 108552

Fig. 11. Contour of Reynolds number (left: LBE, middle: sodium, right: water).

3.4.3. Influence of power on the stability under asymmetric conditions occurrence of flow instability. Therefore, this model cannot accu-
There are big differences when the system is operating in differ- rately describe the changes in the stability boundary at low power.
ent power ranges, and the natural circulation are not completely
the same. In order to clarify the influence of power on stability, this
4. Conclusion
study is based on Loop1, using lead–bismuth as coolant, inlet tem-
perature of 483 K, flow rate of 1 m/s, and analyzing the changes of
In this study, a dimensionless analysis method was employed to
stability boundary by changing the system operating power. The
derive the Jacobian matrix of a single-loop heat transport system
result is shown in Fig. 12:
with natural circulation. Based on this matrix, the Jacobian matrix
It can be seen from Fig. 12 that the stability boundary continues
model for a parallel channel-loop system was constructed. Accord-
to rise with the power increases, and the system operating state
ing to the model, the natural circulation stability of parallel
also continues to rise. The operating state will exceed the stability
channel-loop system under different load and resistance was ana-
boundary when the power exceeds the limit, and as the power
lyzed, and the effect of loop geometry and coolant on the stability
increases, the more it deviates from the stability boundary. There-
of the system under asymmetric condition was investigated. The
fore, there is a maximum power limit and an optimal power for a
main results of the study are summarized as follows:
particular loop, which can be used for evaluation during the design
of the reactor.
(1) The non-dimensional analysis method was used to analyze
But in the low power range, the system operating state slowly
the stability boundary of the single-loop natural circulation
approaches 0 when the power gradually approaches 0, the curve
system, and the Jacobian matrix, which characterized this
is smoother, and the stability margin is larger. Fig. 12 shows that
system, was obtained and compared with the results of
the system has better stability at low power. However, the system
other researchers. The obtained results were in good agree-
is a steady-state ideal model without external interference,
ment with the previous reports, which proved the accuracy
thereby it does not consider the interference from external factors,
and reliability of the proposed model.
nor does it consider the transient effects. Since the natural circula-
(2) A parallel channel-loop Jacobian matrix model was con-
tion state of the system is weak when operating in the low power
structed based on the Jacobian matrix theory of single-
range, the anti-disturbance ability is poor, and the stability bound-
loop, and the changes in the stability of parallel channel-
ary is low. At this time, the introduction of an external disturbance
loop system under different load and resistance were ana-
will have a greater impact on the operating state, which lead to the
lyzed. The results showed that for parallel channel-loop sys-
tem, there were two critical Reynolds numbers. When the
load difference generated in the left and right loops was
greater than these two critical Reynolds numbers, the sys-
tem became unstable. The range of the stability region could
be increased by increasing the aspect ratio and the length of
the heating zone and reducing the pipe diameter and the
length of the cooling zone. In addition, the stability bound-
ary was found to be more sensitive to the aspect ratio and
the length of the heating zone.
(3) The change in the stability of the left and right loops under
the asymmetric resistance difference was analyzed based
on the parallel channel-loop Jacobian matrix model. In addi-
tion, the influence of different geometric characteristics and
coolants on the stability boundary of the circuit was ana-
lyzed. The results showed that the system stability can be
improved within the allowable range of natural circulation
by choosing liquid metal coolant and increasing the pressure
drop, and this model cannot accurately describe the changes
Fig. 12. Curve of stability boundary versus power variation. in the stability boundary at low power.
8
E. Zhu, T. Wang, P. Zhao et al. Annals of Nuclear Energy 163 (2021) 108552

Overall, this study provides useful insights on the design of the References
parallel channel-loop system, and the obtained results can be used
Hou, X., Sun, Z., Lei, W., 2017. Capability of RELAP5 code to simulate the thermal-
to improve the stability of natural circulation. The stability bound- hydraulic characteristics of open natural circulation. Ann. Nucl. Energy 109,
ary can be effectively used to preliminarily assess the stability dis- 612–625.
tribution and the variation trend of the system. Lee, J.D., Pan, C., 2005. Nonlinear analysis for a nuclear-coupled two-phase natural
circulation loop. Nucl. Eng. Des. 235 (5), 613–626.
Satoh, A., Okamoto, K., Madarame, H., 1998. Instability of single-phase natural
circulation under double loop system. Chaos, Solitons Fractals 9 (9), 1575–1585.
CRediT authorship contribution statement Demazière, C., Marcel, C., Rohde, M., van der Hagen, T., 2008. Multifractal analysis of
chaotic flashing-induced instabilities in boiling channels in the natural-
circulation CIRCUS facility. Nucl. Sci. Eng. 158 (2), 164–193.
Enping Zhu: Data curation, Writing - original draft. Ting Wang:
Marcel, C.P., Rohde, M., Van Der Hagen, T.H.J.J., 2010. Experimental investigations
Conceptualization, Methodology, Software. Pengcheng Zhao: on flashing-induced instabilities in one and two-parallel channels: A
Supervision. Shengwen Deng: Visualization, Investigation. Shen- comparative study. Exp. Therm Fluid Sci. 34 (7), 879–892.
Nadella, S., Srivastava, A.K., Maheshwari, N.K., 2018. A semi-analytical model for
quan Zeng: Software, Validation. Jie Li: Writing - review & editing.
linear stability analysis of rectangular natural circulation loops. Chem. Eng. Sci.
192, 892–905.
Zhu, H., Yang, X., Gong, H., Jiang, S., 2013. Theoretical and experimental study on
Declaration of Competing Interest single-phase natural circulation under inclined conditions. J. Nucl. Sci. Technol.
50 (3), 304–313.
Lu, D., Zhang, X., Guo, C., 2014. Stability analysis for single-phase liquid metal
The authors declare that they have no known competing finan- rectangular natural circulation loops. Ann. Nucl. Energy 73, 189–199.
cial interests or personal relationships that could have appeared Cammarata, L., Fichera, A., Pagano, A., 2003. Stability maps for rectangular
to influence the work reported in this paper. circulation loops. Appl. Therm. Eng. 23 (8), 965–977.
Goyal, V., Hassija, V., Pandey, V., Singh, S., 2020. Non-linear dynamics of single
phase rectangular natural circulation loop. Prog. Nucl. Energy 130, 103530.
[Link]
Acknowledgements Wu, Q., Sienicki, J.J., 2003. Stability analysis on single-phase natural circulation in
Argonne lead loop facility. Nucl. Eng. Des. 224 (1), 23–32.
This work is supported by the ‘‘Coupled Response Mechanism of Zhang, Y., Wang, C., Lan, Z., Wei, S., Chen, R., Tian, W., Su, G., 2020. Review of
thermal-hydraulic issues and studies of lead-based fast reactors. Renew.
Lead-based Fast Reactor Hot Pool and Cold Pool Thermal stratifica- Sustain. Energy Rev. 120, 109625. [Link]
tion under Asymmetric Thermal Load Condition and its Influence Nayak, A.K., Vijayan, P.K., Saha, D., Venkat Raj, V., 1995. Mathematical modelling of
on Natural Circulation Performance” of National Natural Science the stability characteristics of a natural circulation loop. Math. Comput. Modell.
22 (9), 77–87.
Foundation of China (Grant No. 11905101) and ‘‘Characteristics Zhou, X., Xu, Y., Yao, J., Zheng, K., Zhao, Y., 2017. Stiffness modelling and comparison
Analysis and Influence Research of Core Power Redistribution and of the 5-UPS/PRPU parallel machine tool with its non-redundant counterpart.
Pool Complex Thermal Stratification for Pump-driven Lead- Proc. Inst. Mech. Eng. B: J. Eng. Manuf. 231 (9), 1646–1657.
Karami, I., Aghaie, M., 2017. Sensitivity analysis of numerical schemes in natural
cooled Fast Reactor under Asymmetric Operation Conditions” of
cooling flows for low power research reactors. Adv. Energy Res. 5 (3), 255.
Natural Science Foundation of Hunan Province, China (Grant No.
2019JJ40239).