The analysis of density wave oscillation in ocean motions with a density variant drift-flux model

International Journal of Heat and Mass Transfer 115 (2017) 138–147

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

The analysis of density wave oscillation in ocean motions with a density variant drift-flux model B.H. Yan ⇑, R. Li, L. Wang Sino-French Institute of Nuclear Engineering and Technology, Sun Yat-sen University, Zhuhai, Guangdong, China

a r t i c l e

i n f o

Article history: Received 12 April 2017 Received in revised form 26 July 2017 Accepted 9 August 2017

Keywords: Ocean motion Density wave oscillation Density variation Drift-flux model Bifurcation

a b s t r a c t The two-phase flow instability in heated channels is an important phenomenon in industrial facilities. The effect of ocean motions upon the flow instability should be predicted accurately for the heat exchangers in floating platform. The reduced order model with drift-flux approach is developed to analyze the density wave oscillation in a uniformly heated channel in ocean motions. The density variation of coolant flow is considered. The stability maps can be predicted precisely on experimental data. The effect of ocean motions on the stability boundary is investigated quantitatively by the defining of the degree of instability. The subcritical Hopf bifurcation in ocean motions is also studied. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction The stability of two-phase flow in a uniformly heated channel is important in several industrial domains like refrigeration systems, turbo-machinery, boiling water reactors, solar steam generating systems and two-phase flow heat exchangers [1]. In the case of a reactor or a heated channel in a floating platform, the ocean motions could introduce an additional force to the system and induce the periodical fluctuation of primary coolant in the reactor core and destroy system stability. Therefore, it is important to clarify the effect of ocean motions upon the two-phase flow instability. In nuclear reactors and thermal systems operated in high pressure, density wave oscillation (DWO) is observed most frequently. It is found that the marginal stability boundary (MSB) of regular complex flow oscillation is similar to that of the DWO in stationary state. The influence of rolling parameters on the MSB is not significant [2,3]. The models of two-phase flow instability between multi-channels (FIBM) in ocean motions were developed and used to obtain the instability oscillating trajectories of multi-channel systems. The instability zone of a nine-channel system in rolling motion was also analyzed and it was found that some of the trajectories showed chaotic characteristics [4]. Then these scholars analyzed the complex curve of mass flow with Fast Fourier Transformation (FFT) method and analyzed the onset of inherent parallel-channel instability [5]. A parallel nine- channel model in ⇑ Corresponding author. E-mail address: [email protected] (B.H. Yan). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.08.022 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

rolling motion was established based on homogeneous flow model [2]. The marginal stability boundary in rolling motion in that model was obtained numerically. It was observed that the unstable regions occurred in both low and high equilibrium-quality regions. In high equilibrium-quality region, the multiplied period phenomenon was observed and the chaotic phenomenon appeared on the right of marginal stability boundary. The flow instability of forced circulation in a mini-rectangular channel in rolling motion was also investigated experimentally [6]. Some researchers found that the flow oscillation in ocean motions was the superposition of thermal-induced oscillation and motion-induced oscillation [7]. It was found that the difference of boundary heat flux induced by ocean motions was not significant by comparing experimental results in stationary states and ocean motions. The instability was mainly affected by thermal parameters rather than ocean motions. An empirical formula which was capable of predicting the instability boundary under both static and motion conditions was also summarized. Shen et al developed a lumped model of flow instability in parallelrectangular channels in heaving motion [8]. Their results indicated that heaving amplitude contributed slightly to the flow instability, which was in agreement with the results of [3,7]. It was also found that the threshold power was decreased as the heaving frequency which was next to the system characteristic frequency. Until now, the numerical investigations for the two-phase flow instability in ocean motions were mainly carried out with time domain method and homogeneous equilibrium model (HEM) [9]. A one-dimensional HEM was used as the basis for the development

B.H. Yan et al. / International Journal of Heat and Mass Transfer 115 (2017) 138–147

Nomenclature a C0 D g Fr h Dhfg k L Nf Npch Nsub DP qw s u

vgj

x

phase variable for single-phase enthalpy void distribution parameter channel diameter gravitational acceleration Froude number coolant enthalpy latent heat loss coefficient channel length Friction number phase change number Subcooling number pressure drop heat flux phase variable for two-phase quality velocity drift velocity vapor quality

Greek letters perimeter onset of boiling density

e l q

Superscripts – dimensionless parameter Subscripts 1/ single-phase 2/ two-phase acc acceleration exit channel exit ext external f liquid g vapor m mixture r reference

of a nodal model of a boiling channel. Then it was found that the simple model yielded a rich variety of nonlinear phenomenon [10]. This moving nodes model was also introduced by Bertodano et al in the analysis of the occurrence of sustained oscillations under operating conditions close to the linear stability margins [11]. However, it is known that HEM may not be useful to low flow rate conditions. Therefore, this approach was extended by Dokhane using drift flux model (DFM) [12]. Although these models predict the stability characteristics of the system having good agreement with experimental data sets [13], the analysis is difficult due to the large and tedious form of the set of coupled ordinary difference equations [14]. Besides that, the single-phase density of the coolant flow is considered to be constant (saturation density) in many aforementioned studies. Due to the fact that the density of a subcooled liquid is a function of temperature and enthalpy in these kinds of isobaric systems, it arises a discrepancy between the accounted density and actual density, which results in an unsatisfactory agreement between experimental data and numerical results. Therefore, the density variation should be included to improve the accuracy of prediction. According to the Liapunov theorem [15], the stability of a linearized system corresponds to the stability of the nonlinear system operating under quasi-equilibrium conditions. Therefore, the analysis with reduced order model with linearization analysis in quasi-equilibrium conditions is satisfactory. In this work, the DWO in a uniformly heated channel in ocean motions is investigated by using reduced order model. This model does allow a deep insight into the complex processes determining two-phase flow instability, and provides a valuable tool for fast and detailed semi-analytical bifurcation analysis [12,14,16]. This is helpful for the analysis of two-phase flow instability in ocean motions in the future. The density variant drift flux equations were introduced to correlate the two-phase properties. The simulation results were validated with experimental data and some valuable results were obtained and analyzed. 2. Theoretical models 2.1. Flow channel and assumptions The analyzed system is a uniformly heated channel, as shown in Fig. 1, which is identical to the test section designed by Saha [13], for the convenience of validation with experimental data. The

Fig. 1. Schematic of the uniformly heated channel [11].

139

140

B.H. Yan et al. / International Journal of Heat and Mass Transfer 115 (2017) 138–147

analyzed system is a 9 feet long circular tube made of stainless steel with outer diameter of 0.5 inch and inner diameter of 0.402 inch. The axial heat flux was essentially uniform because of the constant wall thickness. The channel length is L. The singlephase coolant enters from the bottom with a velocity of uin and starts boiling at a certain point inside the channel. The coolant continues as a two-phase mixture of liquid and vapor beyond the onset of boiling (ONB). The inlet temperature is T in and the coolant temperature increases linearly in the single-phase region before ONB. The coolant temperature remains as the saturation temperature after the ONB point, despite the quality increases. The channel size is obtained from [13] in order to validate the calculation results with experimental data. In ocean motions, the rolling and pitching axis is located at the center of the channel bottom. The heated channel is analyzed by a 1-D axial flow model using three basic conservation partial difference equations (PDE) of mass, energy and momentum along the channel length. These PDEs are reduced to corresponding ODEs, using weighted residual procedure [14]. These subsequent assumptions are made in these subsequent models: 1. The system pressure remains constant. 2. The heat flux is uniformed along with the channel length. 3. Boussinesq approximation is valid for density variation in twophase region. 4. Subcooled boiling is neglected. 2.2. Energy balance equation

2.3. Mass balance equation The void propagation equation in non-dimensional form in ocean motions is written as:

@ aðz; tÞ @ aðz; tÞ þ ½C 0 jðz; tÞ þ V gj  ¼ Npch ½Nr  C 0 aðz; tÞ @z @t

ð7Þ

jðz; tÞ ¼ qin uin ðtÞ þ Npch ½z  lðtÞ

ð8Þ

aðz; tÞ ¼

jðz; tÞ xðz; tÞNr  ½C 0 jðz; tÞ þ V gj xðz; tÞ þ Nq N r

ð9Þ

where a is void fraction. x is vapor quality. C 0 is void distribution parameter. V gj is vapor drift velocity. The vapor quality is written as:

xðz; tÞ ¼ sðtÞ½z  lðtÞ

ð10Þ

Integrating Eq. (7) from z ¼ lðtÞ to z ¼ 1, one can obtain dsðtÞ=dt as:

ds F 1 ðtÞ  F 3 ðtÞ ¼ dt F 2 ðtÞ

ð11Þ

where F 1 ðtÞ, F 2 ðtÞ and F 3 ðtÞ are functions of phase variables and parameters as defined in Appendix A. 2.4. Momentum conservation equation The dimensionless momentum conservation equation in singlephase region in ocean motions could be written as:

With these aforementioned assumptions, the dimensionless energy conservation equation for incompressible liquid in ocean motions can be written as:

q

@ qh @ quh þ ¼ Npch Nq Nr @z @t

ð1Þ

! ! ! ! ! ! ! ! F a;1U ¼ q  a ¼ q  ½ b  r þ x  ð x  r Þ þ 2 x  u 

q ¼ N1  N2 h

ð2Þ

where P1U is the pressure drop in single-phase region. Nf ;1U is the

means dimensionless parameter. q is fluid

single-phase friction number. Fr is the Froude number. F a;1U is the additional force. The subscript 1/ and 2/ denote single-phase and ! ! two-phase, respectively. x and b stand for rolling angular velocity ! and acceleration, respectively. r is the vector radius. It is usually assumed that the rolling motion is in a cosine order [17]. The rolling angular velocity and acceleration can be written as:

where the superscript



density, h is fluid enthalpy, u is flowing velocity, t is time, z is coordinate in flowing direction, N pch is phase change number. The expressions of Nq and Nr are listed in the Appendix A. N 1 and N 2 are dimensionless constants. The reason that the density is expressed in such form could be found in Ref. [14]. It should be noted that these two constants should be changed in cases with different pressures and coolant flows. The enthalpy is expressed as:

h ¼ hi þ aðtÞz

0

# @ quh þ  Npch Nq Nr dz ¼ 0 @z @t

ð4Þ

where

ð5Þ

b¼

4p2 hm

sin

T2

ð6Þ

where N sub is sub-cooling number with the expression shown in the Appendix A.

2pt T

ð12Þ

ð13Þ

ð14Þ

ð15Þ

where hm and T are rolling amplitude and period, respectively. Integrating Eq. (12) from z ¼ 0 to z ¼ lðtÞ, the equation is changed to:

duin þ F 4 ðtÞ dt

ð16Þ

The dimensionless momentum conservation equation in twophase region can be given as:

DP2U ¼ F 5 ðtÞ

lðtÞ is the location of ONB and could be expressed as:

aðtÞlðtÞ ¼ N q Nr Nsub

2phm 2p t cos T T

DP1U ¼ qin l

Then, one can obtain daðtÞ=dt as:

Npch Nq Nr  qin uin a daðtÞ 2 ¼ l ðN1  2N2 hin Þ  4N2 al=3 dt

x¼

ð3Þ

aðtÞ is phase variable. Substituting Eq. (3) into Eq. (1) and integrating from z ¼ 0 to z ¼ lðtÞ, the equation is changed to:

Z lðtÞ " @ qh

@u @u @P1U þ qu ¼   Nf ;1U qu2  qFr 1 sinðh0 þ hr Þ  F a;1U @z @z @t

duin þ F 6 ðtÞ dt

ð17Þ

DPin ðtÞ ¼ K in qin uin 2 ðtÞ

ð18Þ

DPexit ðtÞ ¼ K exit qm ð1; tÞum 2 ð1; tÞ

ð19Þ

141

B.H. Yan et al. / International Journal of Heat and Mass Transfer 115 (2017) 138–147

where DPin and DP exit are the inlet and exit channel pressure drop, respectively. K in and K exit are inlet and exit local loss coefficients, respectively. Then, one can obtain duin ðtÞ=dt as:

duin ðtÞ=dt ¼ F 7 ðtÞ

ð20Þ

The expressions of intermediate variables F 4 ðtÞ, F 5 ðtÞ, F 6 ðtÞ and F 7 ðtÞ are given in Appendix A. 2.5. Stability analysis The dynamical system consists of three ODEs (ordinary differential equations), viz. Eqs. (5), (11) and (20), which can be written in a compact form as:

dXðtÞ ¼ FðXðtÞ; gÞ dt

ð21Þ

here XðtÞ is the vector of phase variables, and g is the vector of parameters that includes both operating and design parameters as defined in Appendix A.

XðtÞ ¼ ðaðtÞ; sðtÞ; uin ðtÞÞ

ð22Þ

g ¼ ðN1 ; N2 ; Nq ; Nr ; Npch ; Nsub ; K in ; K exit ; Fr; Nf ;1U ; Nf ;2U ; F a;1U ; F a;2U Þ ð23Þ For a given set of parameters an equilibrium (or steady state) point can be expressed as:

FðX 0 ; g0 Þ ¼ 0

ð24Þ

Then, the stability of the system is obtained by analyzing the eigenvalues of the Jacobian matrix Dx FðX 0 ; g0 Þ obtained from the ODEs. If all the eigenvalues of Jacobian matrix have negative real part, then the system is stable. However, if the real part of at least one eigenvalue is positive, the system is unstable for the given set of parametric values. Whereas if the real parts of the eigenvalue are zero, the system is said to have an oscillatory behavior [16]. At the stability boundary, the real part of the eigenvalue is zero. The stability boundary is obtained by solving the following equations along with the ODEs. In ocean motions, the forces contributing to the fluid are different from that in stationary state. As a result, the flow and heat transfer should be revaluated due to the effect of additional force. However, some scholars suggested that the traditional empirical correlations are applicable in ocean motions [17–19]. Their results were also validated with experiments. Therefore, the traditional correlations and models are introduced but the additional forces are included in the momentum equations.

In this work, a parameter R is defined to describe the degree of instability, which could be expressed as:

R ¼ Aun =Adia

ð25Þ

R is the ratio of the area of unstable region and the total area in the two-phase operating region which is also called as the diamond region. If the stability boundary moves to the stable region, the stable region decreases and the parameter R increases. On the contrary, if the stability boundary moves to the unstable region, the stable region increases and the parameter R decreases. As for a thermal hydraulic system, R is the smaller the better. It should be noted that the value of R is also affected by the maximum N sub . Therefore, the comparison of different R should be conducted in cases with identical maximum N sub . In this work, the definition of degree of instability aims at the quantitative analyzing of the stability boundary. 3.2. Validation with experimental data Saha obtained the DWO instability map experimentally, using Freon-113 as the operating fluid in a uniformly heated channel in stationary platform [13]. Since the experimental data in ocean motions in open literature is mainly about the identification of the types of instability, this instability map is used to validate the mathematical model. Fig. 3 shows the comparison with experimental data of sets I, III and IV, respectively. The results show a very close match with that in [13]. The qualitative behavior is similar to that of Mishra and Singh [16]. Compared with the models used by Karve [20] and Uddin and Dorning [21], the density variant drift-flux model captures the stability boundary much better. Although the current model and the Uddin and Dorning model are based on the same drift-flux approach, they are constructed differently. The Uddin and Dorning model is an exact model but ignored the density variation of coolant flow. The present DFM is based on an assumption implying an approximate spatial treatment, viz. that the single-phase enthalpy and the two-phase quality have a linear dependence on the spatial z direction. Besides that, the density variation of coolant flow along the channel is simulated. This partly explains why the current model fits the data better than the Uddin and Dorning model for low values of subcooling number. In summary, the current thermal hydraulic model has been found to be in satisfactory agreement with the experimental data

10 9

Maximum Nsub

8

3. Results and discussions

7

3.1. Degree of instability

X eq =0

Nsub

6

It is evident from past studies that a thermal hydraulic system exhibit two-types of instabilities, classified as Type-I and Type-II [14]. The Type-I instability arises mainly at low pressure conditions where the steam quality is low. In high pressure conditions, the Type-II instability is more important. As for the nuclear power system, its operating pressure is usually very high. Therefore, only the Type-II instability is investigated in this work. The two-phase operating region of a thermal hydraulic system could be expressed by Npch-Nsub plane. It is enclosed by X eq ¼ 0, X eq ¼ 1, N sub ¼ 0 and N sub ¼ Max, as shown in Fig. 2. The unstable region is the area enclosed by X eq ¼ 1, N sub ¼ Max and the stability boundary. The other operating conditions belong to the stable region.

5

Xeq =1

Unstable

4 3

Stable

2 1 0 0

2

4

6

8

10

12

14

16

Npch Fig. 2. The stable region and unstable region.

18

20

22

142

B.H. Yan et al. / International Journal of Heat and Mass Transfer 115 (2017) 138–147

10

10 Cal results Exp data Mishra and Singh Karve Uddin and Dorning

9 8 7

8 7 6 Nsub

Nsub

6 5

5

4

4

3

3

2

2

1

1

0

Cal results Exp data Karve Saha

9

4

6

8

10

12

0

14

4

6

8

10

Npch

Npch

(a) Set I

(b) Set III

12

14

10 Cal results Exp data Paul and Singh

9 8

Nsub

7 6 5 4 3 2 1

7

8

9

10 Npch

11

12

(c) Set Fig. 3. Validation with experimental data.

Table 1 Margin of ocean motions. Parameters

Rolling

Pitching

Heaving

Heeling

Amplitude period

45° 3–28 s

20° 3–14 s

1.0 g 3–14 s

45° –

for large sub-cooling number, as is the case for several earlier developed models. For lower values of sub-cooling number, the current model agrees better than most of others, so that its validation against the experimental data can be considered satisfied. 3.3. Effect of ocean motions The magnitudes of ocean motions of a floating platform have been defined by Lewis [22] and Hiroshi and Kazuo [23]. However, some amplitudes introduced by Hiroshi and Kazuo are too serious to be observed in engineering application. Therefore, the magni-

tudes of ocean motions in [9,17,23] are introduced, as shown in Table 1. Since the length of a floating platform is usually several times greater than its width, the rolling amplitude is higher than the pitching amplitude. In fact, the rolling amplitude of a floating platform is usually less than 30° or 35° since the floating platform with a nuclear reactor is very heavy. The rolling motions with amplitudes greater than 35° could only been observed in small ships. Heeling motion could change the inclination angle of the heated channel and contributes to the gravitational pressure drop. Fig. 4 shows that heeling motion and the variation of inclination angle contribute slightly to the stability boundary. The variation of stability boundary is very limited. The variation of the degree of instability is less than 1% as the heeling angle increases from 15° to 45°. Therefore, the effect of gravitational pressure drop induced by heeling motion is not significant. In heaving motion, the coolant flow in the heated channel is affected by periodical gravitational acceleration. Therefore, the stability boundary oscillates periodically as shown in Fig. 5. The max-

143

B.H. Yan et al. / International Journal of Heat and Mass Transfer 115 (2017) 138–147

8

8

t=0 t=T/4 t=T/2 t=3T/4

o

θr=0

7

7

o

θr=15

o

θr=30

6

6

o

5

Nsub

Nsub

θr=45

5

4

4

3

3

2

2

1

6

7

8

9

10

11

12

13

14

1

15

6

7

8

9

10

Npch

11

12

13

14

15

Npch Fig. 7. Stability boundary in pitching motion.

Fig. 4. Stability boundaries in heeling motions.

8 t=0 t=T/4

7

t=T/2 t=3T/4

Nsub

6 5 4 3 2 1 5

6

7

8

9

10

11

12

13

14

15

Npch Fig. 5. Stability boundary in heaving motion.

Fig. 8. Parameter R in pitching motion.

8 t=0 t=T/4

66

7

t=T/2 t=3T/4

65

6

Nsub

R(%)

64 63

4

62

3

61

2

60 59 0

5

0.1

0.2

0.3

0.4

0.5

0.6

0.7

t/T Fig. 6. Parameter R in heaving motion.

0.8

0.9

1

1 5

6

7

8

9

10

11

12

Npch Fig. 9. Stability boundary in rolling motion.

13

14

15

144

B.H. Yan et al. / International Journal of Heat and Mass Transfer 115 (2017) 138–147

Fig. 10. Parameter R in rolling motion.

Fig. 11. Location of five typical points.

imum heaving acceleration is 1.0g. Fig. 6 describes the variation of the degree of instability. The parameter R oscillates in a sinusoidal form with the period identical to heaving period. Both the heaving motion and heeling motion contribute to the stability boundary by changing the gravitational pressure drop. The effect of heaving motion on the stability boundary is more significant than that of heeling motion. The variation of stability boundary is next to 5% within a heaving period. The pitching motion is similar to the rolling motion. The pitching amplitude is several times less than the rolling amplitude because of the structure characteristics of floating platforms. Therefore, the effect of pitching motion shown in Figs. 7 and 8 is less than that of rolling motion shown in Figs. 9 and 10. The amplitudes and periods shown in Figs. 7 and 8 and Figs. 9 and 10 are 20° & 3 s and 35° & 3 s, respectively. Figs. 7 and 8 and Figs. 9 and 10 reveal that the effect of rolling and pitching motions on the stability boundary is limited, which is in agreement with the experimental results in Refs. [3,7]. The oscillating periods of the degree of instability in rolling and pitching motions are half of the rolling and pitching periods. That is because the rolling and pitching motions contribute to the stability bound-

ary by the additional pressure drop which is determined by the square of angular velocity shown in Eq. (14). The oscillating period of the square of angular velocity is just equal to the half of the rolling and pitching period. In order to find the system behavior beyond the local bifurcation, numerical integrations have been carried out at different regions in the predicted stability map using 4th order RungeKutta method to the aforementioned set of ODEs. The time series of inlet velocity surrounding the stability boundary in heaving motion is investigated to analyze the effect of ocean motion on the local bifurcation phenomenon. The locations of five typical points are listed in Fig. 11. Fig. 12a–d represent the time evolution of inlet velocity of the coolant flow for different sets of parameter values denoted by point A1, A2 and A3 respectively. Point A1 being in the stable side of the SB, the time series (Fig. 12a) shows a damping nature of oscillation and correspondingly the applied perturbation dies out soon. Point A2 being in the stable region of stability map near the stability boundary, the perturbation applied grows with time (Fig. 12b). This indicates the presence of an unstable limit cycle, as shown in Fig. 12c, in the vicinity of stability boundary on the stable side. Point A3 being present away from the SB in the unstable region shows a diverging nature of oscillation (Fig. 12d). This growth led the system variable into a non-realistic value where the system collapses. Point B1 is at the stable side of the stability boundary. The time series in Fig. 12e shows a damping nature of oscillation for small perturbation. But it shows unstable characteristics and the oscillation becomes diverging in nature, as shown in Fig. 12f, for a comparably large perturbation. This indicates the existence of an unstable limit cycle (Fig. 12g). Point B2 being present in the unstable region of the stability map shows a diverging nature of oscillation (Fig. 12h). The results of pioneer researchers revealed that a supercritical Hopf bifurcation is characterized by the appearance of stable limit cycle solutions inside the linear unstable region close to the stability boundary, while a subcritical Hopf bifurcation is characterized by a stable fixed point and an unstable limit cycle solution inside the linear stable region close to the stability boundary [12]. Thus the presence of an unstable limit cycle in the stable region of the stability boundary indicates the occurrence of Subcritical Hopf bifurcation in that region. The aforementioned analysis reveals that the subcritical Hopf bifurcation could be observed in ocean motions. But the supercritical Hopf bifurcation is not observed despite that both subcritical and supercritical Hopf bifurcations have been observed in stationary state. 4. Conclusions The objective of current research is to develop a model for DWO in ocean motions including the effect of density variation and drift. The reduced order model using drift-flux correlations was developed to analyze the DWO in heated channels in ocean motions. In the present model, the density variation of coolant flow is considered. The stability maps can be predicted closely to experimental data. The effect of ocean motions on the stability boundary is investigated quantitatively with the degree of instability. The occurrence of subcritical Hopf bifurcation in heaving motion is studied using numerical simulation of the system dynamics across the stability boundary. The detailed conclusions could be listed as followings: 1. The current thermal hydraulic models are in satisfactory agreement with the experimental data for both large and small sub-cooling numbers.

B.H. Yan et al. / International Journal of Heat and Mass Transfer 115 (2017) 138–147

(a) Point A1

(b) Point A2

(c) unstable limit cycle for Point A2

(d) Point A3

(e) Point B1 for small perturbation

145

(f) Point B1 for large perturbation

Fig. 12. Time evolution of inlet velocity in heaving motion.

2. The ocean motions do not contribute significantly to the DWO instability. Compared with the rolling, pitching and heeling motions, the effect of heaving motion is the most significant.

3. The oscillating periods of the degree of instability in rolling and pitching motions are half of the rolling and pitching periods due to the effect of angular velocity.

146

B.H. Yan et al. / International Journal of Heat and Mass Transfer 115 (2017) 138–147

(g) unstable limit cycle for Point B1

(h) Point B2 Fig. 12 (continued)

4. The subcritical Hopf bifurcation is observed in ocean motions but the supercritical Hopf bifurcation is not observed.

Acknowledgement

 um ðz; tÞ ¼ jðz; tÞ þ V gj 1 

ðA6Þ

V gj ¼ V gj þ ðC 0  1Þjðz; tÞ

ðA7Þ

F 1 ðtÞ ¼ Npch Nr ð1  lÞ

ðA8Þ

The authors are grateful to Dr. Ashish Mani Mishra in Indian Institute of Technology for providing valuable advices for this work.

F 2 ðtÞ ¼

Conflict of interest statement

F 3 ðtÞ ¼ T 2 þ

We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.

 1 qm ðz; tÞ

Z

Nq N2r ðz  lÞ j dz 2 l C 0 j þ V gj ðx þ N q N r Þ 1

dl duin T3 þ T4 dt dt  xNr j  T2 ¼ z¼1 x þ Nq Nr  Z T3 ¼ 

1

ðA9Þ

ðA10Þ ðA11Þ

(

l

) Nq N2r s Npch V gj xNr j dz þ x þ Nq N r ðC 0 j þ V gj Þ2 C 0 j þ V gj ðx þ Nq Nr Þ2 ðA12Þ

Appendix A

Z

1

xNr qin V gj dz x þ Nq Nr ðC 0 j þ V gj Þ2

A.1. Dimensionless variables and expressions

T4 ¼

The dimensionless variables and expressions used in this work are listed below.

The acceleration pressure drop, frictional pressure drop, gravitational pressure drop and additional pressure drop in single-phase region are given by

t¼ h¼

tur q q z u um uin ; uin ¼ ; ; q ¼ ; qm ¼ m ; z ¼ ; u ¼ ; um ¼ qf L ur L qf ur ur

qg qf h ; Nq ¼ ; Nr ¼ hfg qf Dq

N1 ¼

P 1U

P 2U

N1

qf

; N2 ¼

N2 hfg

qf

; Npch ¼

DP1U0 ¼

ðA1Þ qw DqeL ðhsat  hin ÞDq ; Nsub ¼ ðA2Þ Ahfg ur qf hfg qg

P 1U k1U L u2 F a;1U L ; Fr ¼ r ; F a;1U ¼ ¼ ; Nf ;1U ¼ 2 2D qf ur gL qf u2r P 2U k2U L F a;2U L ; F a;2U ¼ ¼ ; Nf ;2U ¼ 2D qf u2r qf u2r

qm ðz; tÞ ¼ 1  aðz; tÞ=Nr

l

DP1U1 ¼ DP1U2 ¼

ðA3Þ

DP1U3 ¼ ðA4Þ

DPa;1U ¼ ðA5Þ

Z l @q dz u @t 0 Z l 0

Z l 0

Z l

ðA13Þ

ðA14Þ

@u dz @z

ðA15Þ

Nf ;1U qu2 dz

ðA16Þ

qFr1 sin ðh0 þ hr Þdz

ðA17Þ

F a;1U dz

ðA18Þ

qu

0

Z l 0

B.H. Yan et al. / International Journal of Heat and Mass Transfer 115 (2017) 138–147

F 4 ðtÞ ¼ DP1U0 þ DP1U1 þ DP1U2 þ DP1U3 þ DPa;1U

ðA19Þ

The acceleration pressure drop, frictional pressure drop, gravitational pressure drop and additional pressure drop in two-phase region are given by

Z

DP2U1 ¼

1

l

Z

DP2U2 ¼ Nf ;2U Z

DP2U4 ¼

1

Nq

l

Z

DPa;2U ¼

@um dz @z

ðA20Þ

qm um 2 dz

ðA21Þ

qm um 1

l

! @ a V gj 2 dz @z 1  a qm

ðA22Þ

1

ðA23Þ

F a;2U dz

l

! ! ! ! ! ! ! ! F a;2U ¼ qm  a ¼ qm  ½ b  r þ x  ð x  r Þ þ 2 x  u  Z F 5 ðtÞ ¼

1

(

qin ½1 þ C 0 ðqm  1Þ 

l

ðA24Þ )

V gj qin V gj x dz qm x þ Nq Nr ðC 0 j þ V gj Þ2 ðA25Þ

ds dl þ F 9 ðtÞ þ DP2U1 þ DP2U2 þ DP 2U3 dt dt þ DP 2U4 þ DPa;2U

F 6 ðtÞ ¼ F 8 ðtÞ

F 7 ðtÞ ¼

DP ex  DPin ðtÞ  DPexit ðtÞ  F 4 ðtÞ  F 6 ðtÞ qin l þ F 5 ðtÞ Z

F 8 ðtÞ ¼

1

( 

l

) N q Nr ðz  lÞ j dz qm C 0 j þ V gj ðx þ Nq Nr Þ2 V gj

ðA26Þ ðA27Þ

ðA28Þ

( " # Nq Nr s Npch V gj V gj x j F 9 ðtÞ ¼ þ qm x þ Nq Nr ðC 0 j þ V gj Þ2 C 0 j þ V gj ðx þ N q Nr Þ2 l ) Z

1

 Npch ½1 þ C 0 ðqm  1Þ dz

ðA29Þ

Appendix B. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.08.022.

147

References [1] L.C. Ruspini, C.P. Marcel, A. Clausse, Two-phase flow instabilities: a review, Int. J. Heat Mass Transfer 71 (2014) 521–548. [2] Y.J. Zhang, G.H. Su, S.Z. Qiu, W.X. Tian, H. Li, L.B. Qian, Y. Li, X. Yan, Y.P. Huang, Numerical research on oscillation of two-phase flow in multi-channels under rolling motion, Nucl. Eng. Des. 241 (2011) 4704–4713. [3] S.C. Tan, G.H. Su, P.Z. Gao, Experimental study on two-phase flow instability of natural circulation under rolling motion condition, Ann. Nucl. Energy 36 (2009) 103–113. [4] Y. Guo, S.Z. Qiu, G.H. Su, D.N. Jia, The influence of ocean conditions on twophase flow instability in a parallel multi-channel system, Ann. Nucl. Energy 35 (2008) 1598–1605. [5] Y. Guo, C. Chen, H.Y. Zeng, The application of Fast Fourier Transform (FFT) method in the twin-channel system instability under ocean conditions, Ann. Nucl. Energy 37 (2010) 1048–1055. [6] Z.T. Yu, S.C. Tan, H.S. Yuan, C. Chen, X.B. Chen, Experimental investigation on flow instability of forced circulation in a mini-rectangular channel under rolling motion, Int. J. Heat. Mass. Transfer 92 (2016) 732–743. [7] Y. Tang, B.D. Chen, W.Y. Xiong, X.Z. Liu, Comparison of flow instabilities under static condition and marine motion conditions based on experiments, Ann. Nucl. Energy 70 (2014) 11–20. [8] Y.O. Shen, W. Chen, S.H. Ding, L.B. Qian, Theoretical research on effect of heaving condition on flow instability in parallel channels based on Galerkin nodal method, Nucl. Pow. Eng. 36 (2015) 102–107. [9] B.H. Yan, Review of the nuclear reactor thermal hydraulic research in ocean motions, Nucl. Eng. Des. 313 (2017) 370–385. [10] A. Clausse, R.T. Lahey, The analysis of periodic and strange attractors during density wave oscillations in boiling flows, Chaos. Solitons. Frac. 1 (1991) 167– 178. [11] M.L. Bertodano, W. Fullmer, A. Clausse, V.H. Ransom, Two-Fluid Model Stability, Simulation and Chaos, Springer, New York, 2016. [12] A. Dokhane, BWR Stability and Bifurcation Analysis Using a Novel Reduced Order Model and the System Code RAMONA (Ph.D. Thesis), Institute of Physics for Energy and particles, Switzerland, 2004. [13] P. Saha, Thermally Induced Two-Phase Flow Instabilities Including the Effect of Thermal Non-equilibrium Between Phases (Ph.D. Thesis), Georgia Institute of Technology, Georgia, Atlanta, 1974. [14] S. Paul, S. Singh, Analysis of sub- and supercritical Hopf bifurcation with a reduced order model in natural circulation loop, Int. J. Heat Mass Transfer 77 (2014) 344–358. [15] B. Porter, Stability Criteria for Linear Dynamical Systems, Academic Press, New York, 1968. [16] A.M. Mishra, S. Singh, Subcritical and supercritical bifurcations for two-phase flow in a uniformly heated channel with different inclinations, Int. J. Heat Mass Transfer 93 (2016) 235–249. [17] B.H. Yan, Q.L. Wen, Scaling analysis for the ocean motions in single-phase natural circulation, Ann. Nucl. Eng. 75 (2015) 521–526. [18] X.Z. Liu, Y.P. Huang, J. Ma, J. Huang, Study on friction resistance characteristics of narrow rectangular channel single-phase isothermal flow in heaving motion, Nucl. Power Eng. 32 (2011) 71–75. [19] J. Ma, Y.P. Huang, X.Z. Liu, R.J. Li, Fully developed laminar flow in a rolling narrow rectangular duct, Nucl. Power Eng. 34 (2013) 44–50. [20] A. Karve, Nuclear-Coupled Thermal-Hydraulic Stability Analysis of Boiling Water Reactors (Ph.D. Thesis), University of Virginia, 1998. [21] R. Uddin, J.J. Dorning, Some nonlinear dynamics of a heated channel, Nucl. Eng. Des. 93 (1986) 1–14. [22] E.V. Lewis, The Motion of Ships in Waves, Principles of Naval Architecture, New york, 1967. [23] O. Hiroshi, T. Kazuo, The ship design of the first nuclear ship in Japan, Nucl. Eng. Des. 10 (1969) 211–219.