# Single Phase Natural Circulation Behaviors of the Integral Type Marine Reactor Simulator under Rolling Motion Condition

Authors: Hou-jun Gong; Xing-tuan Yang; Yan-ping Huang; Sheng-yao Jiang
Journal: Science and Technology of Nuclear Installations
 (2015)

---

## Abstract

During operation in the sea the reactor natural circulation behaviors are affected by ship rolling motion. The development of an analysis code and the natural circulation behaviors of a reactor simulator under rolling motion are described in this paper. In the case of rolling motion, the primary coolant flow rates in the hot legs and heating channels oscillated periodically, and the amplitude of flow rate oscillation was in direct proportion to rolling amplitude, but in inverse proportion to rolling period. The total mass flow rate also oscillated with half the rolling period, and the average total mass flow rate was less than that in steady state. In the natural circulation under a rolling motion, the flow rate oscillations in the hot legs were controlled by the tangential force; however, the mass flow rate oscillations in the total natural circulation and the heating channels were a result of the combined action of the change of inclination angle, flow resistance, and the extra force arising from the rolling motion. The extra tangential force brought about intense flow rate oscillations in the hot legs, which resulted in increasing total flow resistance; however the extra centrifugal force played a role in increasing thermal driving head.


---

## Body

## 1. Introduction

Apart from designed as land-based nuclear power plants, another significant application of small and medium sized LWRs is to be mounted on barges as floating nuclear power plants FNPPs (FNPPs) or as propulsion power of commercial cargo ship or icebreaker, such as SMART and KLT-40s. However, during operation in the sea, the ship inclination or rolling motion will affect the reactor thermal hydraulics, especially for natural circulation behaviors. The effects of ship motion focus on two principles. Firstly, inclination decreases the height difference between the source (core) and sink (SG) and thus decreases the thermal driving head of natural circulation. Secondly, inertial forces imparted by the ship rolling motion in addition to the gravity act upon the primary coolant, which control the natural circulation flow rate with the loop resistances together. The direction and magnitude of inertial forces depend on the position with respect to rolling axis, rolling amplitude, and period. So the natural circulation behaviors under ship motion conditions are much more complicated than those operated on the land.Experimental studies [1–4] showed that the loop flow rates varied periodically, total mass flow rate decreased, and the heat transfer was enhanced in the case of rolling. Beside experiments, some thermal-hydraulic codes were developed to study the influence of ocean conditions. JAERI developed RETRAN-02/GRAV code, which is the modified version of RETRAN-02 code, to simulate the experimental voyages of Mutsu [5] and analyse the effects of ship motions on natural circulation of deep sea research reactor DRX [6]. Seoul National University [7] developed RETRAN-03/MOV to perform thermal-hydraulic analysis of ship reactor under multidimensional motions.To understand the natural circulation behaviors of the integral type marine reactor under rolling motion conditions, a scaled test facility was built in Institute of Nuclear and New Energy Technology of Tsinghua University (INET) and a one-dimensional time-domain simulation program development was synchronized with test facility design. This paper describes the physical models and solution method in the program and the natural circulation behaviors of the reactor simulator under rolling motion condition based on the analysis results.


## 2. Test Facility

The test facility is a reactor simulator (natural circulation test loop) mounted on rolling equipment shown schematically in Figure1(a) and the geometry parameters are listed in Table 1. To investigate the symmetric characteristics of the real reactor system as shown in Figure 1(b) with a comparable simple structure, the primary loop is designed to have two symmetrical circuits. Each circuit contains a heat exchanger and a downcomer connecting to hot legs and lower plenum. Three electrical heating channels are connected to lower plenum, playing the role of the reactor core and each of them consists of 16 Φ
10
×
1 mm and 1080 mm long stainless steel tubes with the arrangement shown in Figure 1(c). During circulation, fluid flows from the lower plenum to the electrical heater and is heated to be with higher temperature and lower density and then flows up and separates in hot legs; in the heat exchanger of each circuit, the fluid is cooled down and flows back into the lower plenum through the downcomer finishing a complete circulation. The primary loop is installed on rolling equipment and the rolling axis is perpendicular to primary loop.Table 1
Geometry parameters of test facility.
Items
Number
Length (m)
Flow area (m2)
Inlet section
3
0.306
0.00237
Heating channel
3
1.08
0.00178
Mid riser
3
0.575
0.00126
Riser
1
1.59
0.0104
Hot leg
2
0.67
0.0113
Heat exchanger
2
1.6
0.00246
Downcomer
2
1.951
0.00594
Lower plenum
1
1.34
0.0113(a) Configuration of the test facility. (b) Relation between test section and real reactor system. (c) Connection of electrical heating tubes.
(a)
(b)
(c)


## 3. Code Development

### 3.1. Physical Model

Be similar to study the fluid mechanics of rotating machinery that is often best analyzed in a rotating frame of reference, a noninertial frame of reference with the same rolling motion of test facility was chosen in our study. In a noninertial frame of reference, the continuity and energy equations are unchanged but the momentum equation must be modified. The extra body force term that arises from the motion of the noninertial frame must be added to the right of momentum equation to take account of the motion of the frame.One-dimensional single phase conservation equations of hydrodynamic component can be written as follows:mass continuum equation:(1)
∂
ρ
∂
t
+
∂
ρ
u
∂
z
=
0
,
momentum continuum equation:(2)
∂
ρ
u
∂
t
+
∂
ρ
u
2
∂
z
=
-
∂
p
∂
z
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
e
ρ
u
2
2
,
energy continuum equation:(3)
∂
ρ
h
∂
t
+
∂
ρ
u
h
∂
z
=
q
′′
P
H
A
s
.The extra accelerationf
ex arising from rolling was obtained in this general expression,
(4)
f
ex
=
ω
˙
×
r
+
ω
×
ω
×
r
+
2
ω
×
U
,
where ω is angular velocity of noninertial frame, ω
˙ is angular acceleration of noninertial frame, r is radial vector of fluid particle in noninertial frame, and U is velocity of fluid particle in noninertial frame.The three terms in the right side of (4) may each be significant. The first term ω
˙
×
r is the acceleration caused by angular acceleration of the noninertial frame, so it is of little importance for geophysical flows or for flows in machinery that rotate at a constant rate about a fixed axis. However, it does play an important role in rolling motion. The second term ω
×
ω
×
r is the centripetal acceleration, depending strongly on the rotation rate and the distance of the fluid particle from the axis of rotation. The final term 2
ω
×
U is the Coriolis acceleration, depending on the fluid particle’s velocity, not on its position. Because Coriolis acceleration is perpendicular to main stream, it is omitted in the code.Figure2 shows a control volume with three-dimensional position in noninertial frame of reference; the center of control volume is p
=
x
,
y
,
z and unit vector of main flow direction is n
=
cos
⁡
α
,
cos
⁡
β
,
cos
⁡
γ. In the simulation of rolling motion, the X-axis overlapped rolling axis, and motion parameters of the noninertial reference system were the same as rolling, so coordinate of control volume was unchanged and the one-dimensional flow was always along the main direction of control volume in the noninertial reference system.Figure 2
Three-dimensional spatial information of a control volume in noninertial frame of reference.The parameters of rolling motion were defined as follows:rolling angle:(5)
φ
=
φ
m
sin
ϖ
d
t
,
rolling angle speed:(6)
ω
=
φ
m
ϖ
d
cos
⁡
ϖ
d
t
,
rolling angel acceleration:(7)
ω
˙
=
-
φ
m
ϖ
d
2
sin
ϖ
d
t
.The component of extra acceleration along main flow direction in our test facility was described by the following equation:(8)
f
ex
,
z
=
ω
2
y
+
ω
˙
z
cos
⁡
β
+
ω
2
z
-
ω
˙
y
cos
⁡
γ
.During rolling,g
z was changed with the angle between flow channel and vertical direction. Rolling is in some sense a rotation with varying speed around a fixed axis; g
z at the rolling angle φ can be obtained as
(9)
g
z
=
g
cos
⁡
β
sin
φ
-
cos
⁡
φ
+
cos
⁡
γ
cos
⁡
φ
.


### 3.2. Constitute Correlations

The friction factor model used in the code was simply an interpolation scheme linking the laminar, laminar-turbulent transition, and turbulent flow regimes. The laminar friction factor was calculated as(10)
λ
L
=
64
Re
,
for
50
≤
Re
≤
2200
.The friction factor in the transition region between laminar and turbulent flows was computed by reciprocal interpolation as(11)
λ
L
,
T
=
3.75
-
8250
Re
λ
3000
-
λ
2200
+
λ
2200
,
for 2200
<
Re
<
3000.The turbulent friction factor is given by the Zigrang-Sylvester [8]:
(12)
λ
T
=
-
2
lo
g
10
ε
3.7
D
e
+
2.86
Re
-
5.02
Re
lo
g
10
ε
D
e
-
21.25
R
e
0.9
-
0.5
for Re
≥
3000.The wall heat transfer coefficienth
wl was obtained:
(13)
h
wl
=
k
l
D
h
Nu
.In the code,Nu was taken as the maximum of the values for laminar and turbulent forced convection and natural convection. That was,
(14)
Nu
=
max
⁡
N
u
lam
,
N
u
turb
,
N
u
nc
,
where N
u
lam is Nusselt number for laminar forced convection, N
u
turb is Nusselt number for turbulent forced convection, and N
u
nc is Nusselt number for natural convection.The Nusselt number for laminar flow was set to the analytical value for fully developed flow with a constant heat flux boundary condition. Thus,(15)
N
u
lam
=
4.36
.The Nusselt number for fully developed turbulent flow was that of Weisman [9]; it was given by
(16)
N
u
turb
=
0.023
R
e
0.8
P
r
n
C
w
,
where Pr
⁡ is liquid Prandtl number and n equals 0.4 for heating and 0.3 for cooling. Heat transfer correction factor C
w for rods bundle is a function of bundle pitch P
b and rod diameter D
R as
(17)
C
w
=
0.042
P
b
D
R
-
0.024
.The Nusselt number for natural circulation was given by(18)
N
u
nc
=
max
⁡
N
u
nc
,
lam
,
N
u
nc
,
turb
,
where
(19)
N
u
nc
,
lam
=
0.59
Gr
·
Pr
⁡
1
/
4
was used for laminar natural convection, and
(20)
N
u
nc
,
turb
=
0.1
Gr
·
Pr
⁡
1
/
3
was used for turbulent natural convection.The constitute correlations discussed above were widely used in analysis codes. Although these accuracy for real cases of the test loop should be verified, it was the best choice in our studies before deep experimental studying.


### 3.3. Solution Method

To solve (1)–(3) and calculate the accurate density distribution along the test loop, a flow channel with abrupt area change was divided into a group of control volumes shown in Figure 3.Figure 3
Nodalization of a flow channel.For single phase flow, during a small time interval, we can assume(21)
∂
ρ
∂
t
=
0
.Then (1) transforms into
(22)
d
W
d
z
=
0
.Integrating (3) over control volume j, we can obtain
(23)
A
j
Δ
z
j
ρ
∂
h
∂
t
+
h
∂
ρ
∂
t
j
+
W
h
˙
n
j
+
1
/
2
-
h
˙
n
j
-
1
/
2
=
q
j
′′
n
+
1
P
H
,
j
Δ
z
j
.Fluid enthalpy at the joint of control volumes was obtained by using upwind scheme as(24)
h
˙
j
+
1
/
2
=
0.5
h
j
+
h
j
+
1
+
0.5
W
j
+
1
/
2
W
j
+
1
/
2
h
j
-
h
j
+
1
.In the code, temperatureT was chosen as independent state variable. Omitting the effect of pressure on density and enthalpy, we can obtain
(25)
∂
ρ
∂
t
=
∂
ρ
∂
T
n
∂
T
∂
t
,
∂
h
∂
t
=
∂
h
∂
T
n
∂
T
∂
t
.For simplicity, the long expressions were replaced by single variables:(26)
X
j
,
1
=
A
j
Δ
z
j
ρ
j
n
∂
h
∂
T
j
n
+
h
j
n
∂
ρ
∂
T
j
n
,
X
j
,
2
=
q
j
′′
n
+
1
P
H
,
j
Δ
z
j
.Combining (23) and (26), we can obtain
(27)
T
j
n
+
1
-
T
j
n
Δ
t
=
X
j
,
2
-
W
h
˙
j
+
1
/
2
n
-
h
˙
j
-
1
/
2
n
X
j
,
1
.Integrating (2) along flow channel, we can obtain
(28)
∫
L
1
A
∂
W
∂
t
d
l
+
W
out
2
A
out
ρ
out
-
W
in
2
A
in
ρ
in
=
-
P
out
-
P
in
-
∫
L
ρ
g
z
d
l
+
∫
L
ρ
f
ex
,
z
d
l
-
∫
L
λ
D
h
W
W
2
ρ
A
2
d
l
-
H
L
o
s
s
.And the integrating can be solved approximately as follows:(29)
∫
L
1
A
∂
W
∂
t
d
l
=
∑
j
=
1
M
∂
W
n
+
1
∂
t
Δ
z
j
A
j
,
∫
L
ρ
g
z
d
l
=
∑
j
=
1
M
ρ
j
n
+
1
g
z
,
j
n
+
1
Δ
z
j
,
∫
L
ρ
f
ex
,
z
d
l
=
∑
j
=
1
M
ρ
j
n
+
1
f
ex
,
z
,
j
n
+
1
Δ
z
j
,
∫
L
λ
D
h
W
W
2
ρ
A
2
d
l
=
∑
j
=
1
M
Δ
z
j
λ
j
n
D
h
,
j
W
j
W
j
n
2
ρ
j
n
A
j
2
,
H
L
o
s
s
=
k
W
W
2
ρ
A
j
2
.
Integrating (2) along four close loops shown in Figure 4, we can obtain the following:Figure 4
Nodalization of test facility.loop 1 composed of left heating channel, middle risers, right heating channel, and part of lower plenum,(30)
-
Δ
W
1
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
+
L
9
A
3
+
Δ
W
2
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
A
3
+
Δ
W
4
-
Δ
W
1
Δ
t
L
9
A
9
=
∫
loop
1
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
h
W
W
2
ρ
A
2
d
l
-
∑
loop
1
H
L
o
s
s
,loop 2 composed of right heating channel, middle risers, heating channel 2, and part of lower plenum,(31)
-
Δ
W
3
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
+
L
9
A
3
+
Δ
W
2
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
A
3
+
Δ
W
5
-
Δ
W
3
Δ
t
L
9
A
9
=
∫
loop
2
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
h
W
W
2
ρ
A
2
d
l
-
∑
loop
2
H
L
o
s
s
,loop 3 composed of left heating channel, riser, left hot leg, left heat exchanger, downcomer, and part of lower plenum,(32)
Δ
W
1
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
+
L
9
A
3
+
Δ
W
1
+
Δ
W
2
+
Δ
W
3
Δ
t
L
4
A
4
+
Δ
W
4
Δ
t
L
5
A
5
+
L
6
A
6
+
L
7
A
7
+
L
8
A
8
=
∫
loop
3
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
h
W
W
2
ρ
A
2
d
l
-
∑
loop
3
H
L
o
s
s
,loop 4 composed of right heating channel, riser, right hot leg, right heat exchanger, downcomer, and part of lower plenum,(33)
Δ
W
3
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
+
L
9
A
3
+
Δ
W
1
+
Δ
W
2
+
Δ
W
3
Δ
t
L
4
A
4
+
Δ
W
5
Δ
t
L
5
A
5
+
L
6
A
6
+
L
7
A
7
+
L
8
A
8
=
∫
loop
4
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
h
W
W
2
ρ
A
2
d
l
-
∑
loop
4
H
L
o
s
s
.W
1, W
2, W
3, W
4, and W
5 are the mass flow rates of left heating channel, mid heating channel, right heating channel, left hot leg, and right hot leg, respectively, and their positive directions are identical with those in steady natural circulation.Applying the mass continuity, we can obtain(34)
Δ
W
1
+
Δ
W
2
+
Δ
W
3
=
Δ
W
4
+
Δ
W
5
.In the code, Gaussian elimination is used to solve (30)–(34), and final W
n
+
1 is obtained by using
(35)
W
n
+
1
=
W
n
+
Δ
W
.


## 3.1. Physical Model

Be similar to study the fluid mechanics of rotating machinery that is often best analyzed in a rotating frame of reference, a noninertial frame of reference with the same rolling motion of test facility was chosen in our study. In a noninertial frame of reference, the continuity and energy equations are unchanged but the momentum equation must be modified. The extra body force term that arises from the motion of the noninertial frame must be added to the right of momentum equation to take account of the motion of the frame.One-dimensional single phase conservation equations of hydrodynamic component can be written as follows:mass continuum equation:(1)
∂
ρ
∂
t
+
∂
ρ
u
∂
z
=
0
,
momentum continuum equation:(2)
∂
ρ
u
∂
t
+
∂
ρ
u
2
∂
z
=
-
∂
p
∂
z
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
e
ρ
u
2
2
,
energy continuum equation:(3)
∂
ρ
h
∂
t
+
∂
ρ
u
h
∂
z
=
q
′′
P
H
A
s
.The extra accelerationf
ex arising from rolling was obtained in this general expression,
(4)
f
ex
=
ω
˙
×
r
+
ω
×
ω
×
r
+
2
ω
×
U
,
where ω is angular velocity of noninertial frame, ω
˙ is angular acceleration of noninertial frame, r is radial vector of fluid particle in noninertial frame, and U is velocity of fluid particle in noninertial frame.The three terms in the right side of (4) may each be significant. The first term ω
˙
×
r is the acceleration caused by angular acceleration of the noninertial frame, so it is of little importance for geophysical flows or for flows in machinery that rotate at a constant rate about a fixed axis. However, it does play an important role in rolling motion. The second term ω
×
ω
×
r is the centripetal acceleration, depending strongly on the rotation rate and the distance of the fluid particle from the axis of rotation. The final term 2
ω
×
U is the Coriolis acceleration, depending on the fluid particle’s velocity, not on its position. Because Coriolis acceleration is perpendicular to main stream, it is omitted in the code.Figure2 shows a control volume with three-dimensional position in noninertial frame of reference; the center of control volume is p
=
x
,
y
,
z and unit vector of main flow direction is n
=
cos
⁡
α
,
cos
⁡
β
,
cos
⁡
γ. In the simulation of rolling motion, the X-axis overlapped rolling axis, and motion parameters of the noninertial reference system were the same as rolling, so coordinate of control volume was unchanged and the one-dimensional flow was always along the main direction of control volume in the noninertial reference system.Figure 2
Three-dimensional spatial information of a control volume in noninertial frame of reference.The parameters of rolling motion were defined as follows:rolling angle:(5)
φ
=
φ
m
sin
ϖ
d
t
,
rolling angle speed:(6)
ω
=
φ
m
ϖ
d
cos
⁡
ϖ
d
t
,
rolling angel acceleration:(7)
ω
˙
=
-
φ
m
ϖ
d
2
sin
ϖ
d
t
.The component of extra acceleration along main flow direction in our test facility was described by the following equation:(8)
f
ex
,
z
=
ω
2
y
+
ω
˙
z
cos
⁡
β
+
ω
2
z
-
ω
˙
y
cos
⁡
γ
.During rolling,g
z was changed with the angle between flow channel and vertical direction. Rolling is in some sense a rotation with varying speed around a fixed axis; g
z at the rolling angle φ can be obtained as
(9)
g
z
=
g
cos
⁡
β
sin
φ
-
cos
⁡
φ
+
cos
⁡
γ
cos
⁡
φ
.


## 3.2. Constitute Correlations

The friction factor model used in the code was simply an interpolation scheme linking the laminar, laminar-turbulent transition, and turbulent flow regimes. The laminar friction factor was calculated as(10)
λ
L
=
64
Re
,
for
50
≤
Re
≤
2200
.The friction factor in the transition region between laminar and turbulent flows was computed by reciprocal interpolation as(11)
λ
L
,
T
=
3.75
-
8250
Re
λ
3000
-
λ
2200
+
λ
2200
,
for 2200
<
Re
<
3000.The turbulent friction factor is given by the Zigrang-Sylvester [8]:
(12)
λ
T
=
-
2
lo
g
10
ε
3.7
D
e
+
2.86
Re
-
5.02
Re
lo
g
10
ε
D
e
-
21.25
R
e
0.9
-
0.5
for Re
≥
3000.The wall heat transfer coefficienth
wl was obtained:
(13)
h
wl
=
k
l
D
h
Nu
.In the code,Nu was taken as the maximum of the values for laminar and turbulent forced convection and natural convection. That was,
(14)
Nu
=
max
⁡
N
u
lam
,
N
u
turb
,
N
u
nc
,
where N
u
lam is Nusselt number for laminar forced convection, N
u
turb is Nusselt number for turbulent forced convection, and N
u
nc is Nusselt number for natural convection.The Nusselt number for laminar flow was set to the analytical value for fully developed flow with a constant heat flux boundary condition. Thus,(15)
N
u
lam
=
4.36
.The Nusselt number for fully developed turbulent flow was that of Weisman [9]; it was given by
(16)
N
u
turb
=
0.023
R
e
0.8
P
r
n
C
w
,
where Pr
⁡ is liquid Prandtl number and n equals 0.4 for heating and 0.3 for cooling. Heat transfer correction factor C
w for rods bundle is a function of bundle pitch P
b and rod diameter D
R as
(17)
C
w
=
0.042
P
b
D
R
-
0.024
.The Nusselt number for natural circulation was given by(18)
N
u
nc
=
max
⁡
N
u
nc
,
lam
,
N
u
nc
,
turb
,
where
(19)
N
u
nc
,
lam
=
0.59
Gr
·
Pr
⁡
1
/
4
was used for laminar natural convection, and
(20)
N
u
nc
,
turb
=
0.1
Gr
·
Pr
⁡
1
/
3
was used for turbulent natural convection.The constitute correlations discussed above were widely used in analysis codes. Although these accuracy for real cases of the test loop should be verified, it was the best choice in our studies before deep experimental studying.


## 3.3. Solution Method

To solve (1)–(3) and calculate the accurate density distribution along the test loop, a flow channel with abrupt area change was divided into a group of control volumes shown in Figure 3.Figure 3
Nodalization of a flow channel.For single phase flow, during a small time interval, we can assume(21)
∂
ρ
∂
t
=
0
.Then (1) transforms into
(22)
d
W
d
z
=
0
.Integrating (3) over control volume j, we can obtain
(23)
A
j
Δ
z
j
ρ
∂
h
∂
t
+
h
∂
ρ
∂
t
j
+
W
h
˙
n
j
+
1
/
2
-
h
˙
n
j
-
1
/
2
=
q
j
′′
n
+
1
P
H
,
j
Δ
z
j
.Fluid enthalpy at the joint of control volumes was obtained by using upwind scheme as(24)
h
˙
j
+
1
/
2
=
0.5
h
j
+
h
j
+
1
+
0.5
W
j
+
1
/
2
W
j
+
1
/
2
h
j
-
h
j
+
1
.In the code, temperatureT was chosen as independent state variable. Omitting the effect of pressure on density and enthalpy, we can obtain
(25)
∂
ρ
∂
t
=
∂
ρ
∂
T
n
∂
T
∂
t
,
∂
h
∂
t
=
∂
h
∂
T
n
∂
T
∂
t
.For simplicity, the long expressions were replaced by single variables:(26)
X
j
,
1
=
A
j
Δ
z
j
ρ
j
n
∂
h
∂
T
j
n
+
h
j
n
∂
ρ
∂
T
j
n
,
X
j
,
2
=
q
j
′′
n
+
1
P
H
,
j
Δ
z
j
.Combining (23) and (26), we can obtain
(27)
T
j
n
+
1
-
T
j
n
Δ
t
=
X
j
,
2
-
W
h
˙
j
+
1
/
2
n
-
h
˙
j
-
1
/
2
n
X
j
,
1
.Integrating (2) along flow channel, we can obtain
(28)
∫
L
1
A
∂
W
∂
t
d
l
+
W
out
2
A
out
ρ
out
-
W
in
2
A
in
ρ
in
=
-
P
out
-
P
in
-
∫
L
ρ
g
z
d
l
+
∫
L
ρ
f
ex
,
z
d
l
-
∫
L
λ
D
h
W
W
2
ρ
A
2
d
l
-
H
L
o
s
s
.And the integrating can be solved approximately as follows:(29)
∫
L
1
A
∂
W
∂
t
d
l
=
∑
j
=
1
M
∂
W
n
+
1
∂
t
Δ
z
j
A
j
,
∫
L
ρ
g
z
d
l
=
∑
j
=
1
M
ρ
j
n
+
1
g
z
,
j
n
+
1
Δ
z
j
,
∫
L
ρ
f
ex
,
z
d
l
=
∑
j
=
1
M
ρ
j
n
+
1
f
ex
,
z
,
j
n
+
1
Δ
z
j
,
∫
L
λ
D
h
W
W
2
ρ
A
2
d
l
=
∑
j
=
1
M
Δ
z
j
λ
j
n
D
h
,
j
W
j
W
j
n
2
ρ
j
n
A
j
2
,
H
L
o
s
s
=
k
W
W
2
ρ
A
j
2
.
Integrating (2) along four close loops shown in Figure 4, we can obtain the following:Figure 4
Nodalization of test facility.loop 1 composed of left heating channel, middle risers, right heating channel, and part of lower plenum,(30)
-
Δ
W
1
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
+
L
9
A
3
+
Δ
W
2
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
A
3
+
Δ
W
4
-
Δ
W
1
Δ
t
L
9
A
9
=
∫
loop
1
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
h
W
W
2
ρ
A
2
d
l
-
∑
loop
1
H
L
o
s
s
,loop 2 composed of right heating channel, middle risers, heating channel 2, and part of lower plenum,(31)
-
Δ
W
3
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
+
L
9
A
3
+
Δ
W
2
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
A
3
+
Δ
W
5
-
Δ
W
3
Δ
t
L
9
A
9
=
∫
loop
2
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
h
W
W
2
ρ
A
2
d
l
-
∑
loop
2
H
L
o
s
s
,loop 3 composed of left heating channel, riser, left hot leg, left heat exchanger, downcomer, and part of lower plenum,(32)
Δ
W
1
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
+
L
9
A
3
+
Δ
W
1
+
Δ
W
2
+
Δ
W
3
Δ
t
L
4
A
4
+
Δ
W
4
Δ
t
L
5
A
5
+
L
6
A
6
+
L
7
A
7
+
L
8
A
8
=
∫
loop
3
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
h
W
W
2
ρ
A
2
d
l
-
∑
loop
3
H
L
o
s
s
,loop 4 composed of right heating channel, riser, right hot leg, right heat exchanger, downcomer, and part of lower plenum,(33)
Δ
W
3
Δ
t
L
1
A
1
+
L
2
A
2
+
L
3
+
L
9
A
3
+
Δ
W
1
+
Δ
W
2
+
Δ
W
3
Δ
t
L
4
A
4
+
Δ
W
5
Δ
t
L
5
A
5
+
L
6
A
6
+
L
7
A
7
+
L
8
A
8
=
∫
loop
4
-
ρ
g
z
+
ρ
f
ex
,
z
-
λ
D
h
W
W
2
ρ
A
2
d
l
-
∑
loop
4
H
L
o
s
s
.W
1, W
2, W
3, W
4, and W
5 are the mass flow rates of left heating channel, mid heating channel, right heating channel, left hot leg, and right hot leg, respectively, and their positive directions are identical with those in steady natural circulation.Applying the mass continuity, we can obtain(34)
Δ
W
1
+
Δ
W
2
+
Δ
W
3
=
Δ
W
4
+
Δ
W
5
.In the code, Gaussian elimination is used to solve (30)–(34), and final W
n
+
1 is obtained by using
(35)
W
n
+
1
=
W
n
+
Δ
W
.


## 4. Results and Discussion

The node model for test facility is shown in Figure4, and the distance from rolling axis to bottom of test loop was 3.2 m. The coordinates of hydrodynamic volumes’ center and unit vector of main stream in the noninertial frame were added when establishing the input card. Boundary condition of the analysis was that the feed water properties (temperature, pressure, and enthalpy) and the pressure of feed water outlet were constants. The initial condition of transient simulation corresponds to vertical steady state, with the pressure of gas space being 5 MPa and the subcooling of heating channel inlet being 44 K. The mass flow rate W
0 of feed water was determined by a control system in the code. The function of the control system was to satisfy the balance of energy during steady state. If the total heat of heat exchanger transfer exceeded the total heating power, the mass flow rate W
0 was decreased and vice versa.The transient simulation began with the variation of body force due to rolling motion. During the transient simulation, the feed water mass flows and total heating power kept being identical with those in steady state. Analyses were carried out for 3 s, 8 s, 13 s, 18 s, and 23 s rolling periods and 5°, 10°, 22.5°, and 45° rolling amplitude with heating power 175 kw, operating pressure 5 MPa.A case of natural circulation mass flow rates and subcooling of heating channels outlet under rolling motion was shown in Figure5. As was observed, the heat and flow symmetry between heating channels and hot legs were broken. The mass flow rates of heating channels began to oscillate when starting the transient simulation, the oscillation periods of left and right heating channels were the same as rolling motion, but the oscillation period of middle heating channel was half the rolling motion. The same oscillations were also observed in mass flow rates of hot legs with rolling period, but oscillation amplitudes were much larger than those of heating channels. Total natural circulation mass flow rate also oscillated with half the rolling period, and the average within a period was smaller than that of steady state. The oscillation of subcooling of heating channels outlets corresponded to the oscillation of mass flow rates, but the subcooling varied lag related to mass flow rate because of the time needed for water to flow through heating channel. The phenomenon was similar to that found by Murata et al. [2].Figure 5
Natural circulation behavior under rolling condition.Figure6 shows the mass flow rate under same rolling period 8 s, different rolling amplitude. Obviously, the oscillation amplitude of mass flow rate increased as the rolling amplitude increased. When the rolling amplitude equaled 45°, the mass flow rate left heating channel was only 10% of that in vertical steady state. This situation was very harmful and CHF may occur at high power level for a real reactor. When the rolling amplitude was greater than 22.5°, the fluid of hot leg changed the velocity direction in a small range near the trough, which can degrade the quality of steam generator. So it was a great challenge to use natural circulation in primary circuit under severe ocean conditions.Mass flow rate under the same rolling period, different rolling amplitude, (a) left heating channel, (b) left hot leg.
(a)
(b)Figure7 shows the mass flow rate under the same rolling amplitude 22.5°, different rolling period. As was observed, the oscillation amplitude decreased as the period increased. An interesting phenomenon was also observed that the mass flow rate oscillation waveform of left heating channel deviated from sine wave more and more as rolling period increased. However the oscillation waveform of hot legs was very close to sine wave. The clearer phenomenon can be seen in Figure 8. Although the oscillation period of left and right heating channels was equal to that of rolling motion, there were two maximum and minimum in one period and the oscillation period had the trend to be half of rolling period. A discussion of the causes of these forms is given later in this paper.Mass flow rate under the same rolling amplitude, different rolling period, (a) left heating channel, (b) left hot leg.
(a)
(b)Figure 8
Mass flow rate of heating channels under rolling period 23 s, rolling amplitude 45°.Natural circulation flow rates under rolling motion were affected by extra forces, thermal driving head, and flow resistance. The extra force varied with the rolling angle. On the other hand, the thermal driving head also varied with the inclination angle, and flow resistance varied with velocity. Therefore, it was difficult to judge how much do these three factors contribute to overall thermal hydraulic behaviors ostensibly. If the natural circulation was only affected by thermal driving head decrease due to inclination, the normalized total natural circulation mass flow was always less than unit. However, as shown in Figure5, the peak value of normalized total mass flow was greater than unit, so there was influence of rolling motion on natural circulation without doubt.The mechanisms of influence of rolling motion can be interpreted by the analysis of extra forces. Figure9 shows the extra forces induced by rolling motion; the vectors only denote the direction of extra forces; they do not denote the modulus. As was observed, the extra tangential forces F
e along the pipes were in two circles, which were thought to generate theoretical flows along outer and inner loops marked in Figure 9. The tangential force varied periodically with the same period of rolling motion; thus the theoretical flows oscillations were approximate to sine wave. Therefore, the oscillations of W
1, W
3, W
4, and W
5 were the superposition of theoretical flows and flows driven by density difference. The integration of tangential force along outer loop was greater than that of inner loop, and that was why flow rate oscillation amplitude of hot legs was larger than that of heating channels. Due to the intense flow rate oscillation of hot legs, the total flow resistance was larger than that in steady state. The extra centrifugal force F
r was parallel to radial vector, because the rolling axis was higher than the geometric center; thus the effect of centrifugal force was similar to gravity, which caused the peak value of total mass flow to be greater than unit. The centrifugal force and thermal driving head varied periodically with half the rolling period, thus the total mass flow oscillation with half rolling period. Under large rolling period and large rolling amplitude condition, ratio the variation of thermal driving head and tangential force along inner loop become large, that was why the oscillation waveform deviated from sine wave.Figure 9
Analysis of extra forces arising from rolling.From above discussion, it was concluded that in the natural circulation under rolling motion, the flow rates oscillations in hot legs were controlled by tangential force; however, the mass flow rate oscillations of total natural circulation and heating channels were a result of the combined action of the change of inclination angle, flow resistance, and extra force due to rolling motion. The centrifugal force and thermal driving head varied periodically with half the rolling period, and the tangential force varied periodically with the same rolling period, so the oscillation waveforms of left and right heating channels depended on the quantity of the centrifugal force, the tangential force, and thermal driving head, that was why there were two maximum and minimum in one period shown in Figure8.Figure10 shows the average of normalized total mass flow rates of different rolling cases. The average total mass flow rate was less than unit in rolling condition and was in direct proportion to rolling period but was in inverse proportion to rolling amplitude.Figure 10
Average normalized total mass flow rates of different rolling cases.


## 5. Conclusions

A transient system analysis code was developed to study natural circulation under motion condition. The momentum equation was modified and it adopted the fluid mechanics in noninertial frame of reference. The natural circulation behavior of integral type marine reactor simulator under rolling condition can be concluded as follows.(1)
The heat and flow symmetry between heating channels and hot legs were broken. The flow rates of hot legs and left and right heating channels oscillated periodically with the same rolling period; however, the flow rates of middle heating channel and total mass flow rate oscillated with half rolling period. The average total mass flow rate in rolling motion was less than that in steady state.
(2)
The tangential forces along the pipes were in two circles, which generate theoretical flows in outer and inner loops. The superposition of theoretical flows and flows driven by density difference resulted in the oscillation of hot legs and left and right heating channels.
(3)
Natural circulation behaviors in rolling motion were results of the combined action of the variation of inclination angle, flow resistance (extra tangential force), and extra centrifugal force. Extra tangential force caused intense flow rates oscillation of hot legs, which resulted in increasing total flow resistance; however extra centrifugal force played a role in increasing the thermal driving head.



---

*Source: 815603-2015-01-26.xml*