# Study on the Density Wave Instability in Natural Circulation System under Rolling Conditions

Authors: Houjun Gong; Xiaoxi Zhang; Nan Gui; Yanping Huang; Xingtuan Yang; Shengyao Jiang
Journal: International Journal of Energy Research
 (2023)

---

## Abstract

In the ocean environment, the thermal fluid system installed on the ship moves with the waves. This unsteady transient motion imposes additional body forces on the liquid. Thereby, the thermohydraulic properties of the thermal fluid system are changed, especially for natural circulation systems with a low driving force. Flow instability related to system security is a common two-phase flow system phenomenon. It is an essential subject in the study of thermal-hydraulic characteristics. In this paper, the density wave instability in a low-pressure natural circulation system with several parallel heating channels in the circumstances of rolling motion is studied by PNCMC (Program for Natural Circulation under Motion Condition). PNCMC is a newly developed program by adding extra body force induced by ship motions in two-phase’s momentum equations. The mechanism of rolling motion and its impact on flow oscillation and unstable boundary power are investigated. The supercooled boiling caused by the higher surface heat flux firstly promotes the occurrence of density wave oscillations between parallel heating channels. When enough steam enters the riser, the system density wave oscillation occurs under the gravity pressure drop-flow rate-vapor fraction feedback. A more complicated compound oscillation is created by superimposing the interchannel oscillation and the system density wave oscillation. If the rolling motion has little change to the natural circulation flow, when density wave instability arises, the flux oscillation rule and the instability boundary power are basically consistent with the vertical conditions. When the rolling angle is large, the fairly large amplitude flow oscillation generated by rolling will prevent density wave instability from occurring and force the flow oscillation to obey the law of rolling effects.


---

## Body

## 1. Introduction

In the ocean environment, ships move with sea wave, and we generally deconstruct ship motion into a collection of simple movements, including surge, sway, yaw, pitch, roll, and heave, as shown in Figure1. Unlike a stationary system on land, the additional motion in an ocean environment can affect the thermohydraulics of the fluid systems.Figure 1
Decompositions of additional motions.Numerous studies have been conducted on thermal hydraulics for reactors in the ocean environment [1–4]. Recently, there have been several review papers on the thermal-hydraulic properties of reactors in oceanic environments [3–5]. Especially the flow heat transfer and instability issues, since these are important for system safety and dependability. The heat transfer characteristic of ocean conditions is much more complex than that of static states. Murata et al. [6] found that the rolling motion improves the efficiency of heat transmission inside the reactor core due to the increased internal flow. When the intensity of the rolling motion is great, as a result of the rolling motion, inertial force predominates over other factors in the heat transmission process. Natural convection regulates heat transfer when the rolling motion is minimal. The influence of ocean conditions on the thermal-hydraulic characteristics of a passive residual heat removal system is theoretically studied by Xi et al. [1]. It is observed that the amplitudes and periods of oscillations in the system parameters are controlled by the movements of the ocean. Flow heat transfer experiments can be put into two categories: those that utilize natural circulation and those that use forced circulation. For forced circulation, the effect of ocean motions on heat transfer is not as clear as it is for natural circulation [4]. The heat transfer coefficient for natural circulation is impacted by flow oscillations generated by ocean movements and oscillates with the flow together.Yu et al. [7] performed an experiment in a minirectangular channel under rolling circumstances to evaluate the flow instability of forced circulation for flow instabilities, such as trough-type flow oscillation, single-phase or two-phase flow fluctuation, and coupled flow oscillation. The flow instability zone may be separated into three sections based on the effect difference between rolling motion and thermal-hydraulic characteristics: thermal-hydraulic dominant region, rolling dominant region, and coupling resonance region. Natural circulation flow is intrinsically less stable and encounters flow instabilities than forced circulation due to the comparatively small hydraulic driving head [8].Tan et al. [9] carried out an experiment to investigate the fluid instability of natural circulation in a rolling motion. According to their findings, the rolling motion could initiate flow instability early and change the types of instability. The stability of a natural circulation system reduces as the amplitude and frequency of rolling motions increase. Zhang et al. [10] analyzed the nonlinear properties of natural circulation fluid instabilities during rolling movements, which used the same apparatus as Tan et al. [9]. As the frequency or amplitude of a rolling process increases, the threshold for chaotic oscillation rises, in accordance with Tan’s results. Wang and Yang [11] claimed that under the resonance condition, the system stability reduces considerably for type I as well as type II density wave instability. The system destabilizes slightly where the rolling period is rather close to the natural period. The influence of the rolling period on system stability is negligible in the other scenarios.The following serves as a brief introduction to the related work of simulating the problems of rolling and other dynamic motions with a variety of thermal-hydraulic programs in ocean Conditions [12]. Japan Atomic Energy Research Institution developed RETRAN-02/GRAV to simulate natural circulation when the MUTSU nuclear ship is in motion [13]. By upgrading the body force component and process of coordinate transformation related to ocean conditions, its improved version, RETRAN-03, was changed into RETRAN03/MOV [14]. The rolling motion effect, flow resistance, and heat transfer models were then included to RELAP5/MOD3.3 for further refinement [15]. A reactor system analysis code called MARS also has a dynamic motion model that can simulate thermal-hydraulic interactions in three dimensions by using the momentum equation to calculate the body force term. Using conceptual problems, Beom et al. [2] demonstrate that the dynamic motion model is accurate enough to simulate flow. MARS-KS [16] includes a model of dynamic motion for simulating a marine or floating reactor developed by the Korea Atomic Energy Research Institute. Recently, we studied the density wave instability in natural circulation with parallel channels under inclined and heaving conditions [17]. And the stability boundaries under different thermal-hydraulic parameters are investigated to support the safe operation of the system.This research studies density wave instability in a rolling natural circulation system. The PNCMC program [18] is also employed to study the density wave instability of parallel three-channels under simple harmonic rolling, which is described in detail in Sec. 2 [17]. Sec. 3 introduces our research object. The mechanism of the impact of rolling conditions on natural circulation is given in Sec. 4. Sec. 5 discusses the two-phase natural circulation characteristics under rolling conditions. Density wave instability under static and rolling conditions is analyzed in Sec. 6 and Sec. 7, respectively. The flow oscillation law and stability boundary are concluded in Sec. 8.


## 2. PNCMC Program

### 2.1. Governing Equations

The PNCMC program was developed to study thermal hydraulics under marine scenarios using a one-dimensional two-fluid flow model [19]. The two-phase flow equations are modified by adding extra body force terms.Gas momentum equation:(1)αgρg∂Vg∂t+αgρgVg∂Vg∂z=−αg∂P∂z+ΓgVi−Vg−pwαgwτgw+Mgd+αgρggz+fine,z.Liquid momentum equation:(2)αfρf∂Vf∂t+αfρfVf∂Vf∂z=−αf∂P∂z+ΓfVi−Vf−pwαfwτfw−Mgd+αfρfgz+fine,z.


### 2.2. Extra Body Force Terms

This formula gives fluid particle inertial forces,(3)fine=a0+ω̇×r+ω×ω×r+2ω×V.The first terma0 in the RHS accounts for the acceleration of the noninertial frame with respect to a stationary frame. The second term ω̇×r is the acceleration resulting from the noninertial frame’s angular acceleration. The third term ω×ω×r is the radial acceleration, which is highly sensitive to both rotational velocity and the fluid particle’s displacement from the rotational axis. The last element 2ω×V represents the Coriolis acceleration, which is independent of the fluid particle’s location but rather its velocity.Figure2 shows that a noninertial coordinate system OXYZ with unit vectors ijk and a stationary frame OXsYsZs are established. The origin O is positioned on the axis of rolling. In the initial stage, stationary frame OXsYsZs and the noninertial frame OXYZ are overlapping; the unit vector is straight up, which is in the opposite direction of gravity. Noninertial frame OXYZ has a rotational motion around an axis S=Sx,Sy,Sz and acceleration a0=a0,x,a0,y,a0,z. The unit vector n=cosθ,cosβ,cosγ symbolizes the positive mainstream of the channel in the OXYZ. Σ is an arbitrary channel cross-section whose points correspond to the domain x0±Δx,y0±Δy,z0±Δz.Figure 2
Flow channel and noninertial frame geometry.A ship’s rolling is analogous to a sine wave, which may be written as(4)ϕ=ϕmsinϖdt,(5)ω=ϕmϖdcosϖdt,(6)ω̇=−ϕmϖd2sinϖdt,where ϕ, ϕm, ϖd are the rolling angle, amplitude, and frequency, respectively.The reason why the pipeline cross-sectional scaleΔx,Δy,Δz is substantially less than the spatial scale x0,y0,z0 of the system is because the mainstream velocity is much higher than the cross-sectional velocity, fine,n can be interpreted as
(7)fine,n=a0,x−ω2SxSyy0−Sy2x0−Sz2x0+SxSzz0−ω̇Syz0+ω̇Szy0cosθ+a0,y−ω2SySzz0−Sz2y0−Sx2y0+SxSyx0−ω̇Szx0+ω̇Sxz0cosβ+a0,z−ω2SxSzx0−Sx2z0−Sy2z0+SySzy0−ω̇Sxy0+ω̇Syx0cosγ.fine,n is changed to fine,z in two-phase momentum equations. More detailed explanations can be found in Reference [18].During the rolling process,gz shifts depending on the angle between the channel and the vertical direction in a static frame OXsYsZs. Rolling is a variable-speed rotation around a fixed axis, the main stream’s new unit vector n′=cosθ′,cosβ′,cosγ′ in stationary frame OXsYsZs at the angle ϕ of rolling appears to result from the following operations:
(8)Aˆ=SxSxSxSySxSzSySxSySySySzSzSxSzSySzSz,(9)A∗=0−SzSySz0−Sx−SySx0,(10)I=100010001,(11)M=Aˆ+cosϕ·I−Aˆ+sinϕ·Aˆ,(12)n′=n·MT,(13)gz=g·cosγ′.where MT is the transpose matrix of M, and g is the gravity.


### 2.3. Solution and Validation

The two-phase six governing partial differential equations were discretized on a staggered grid by finite volume method with the semi-implicit, first-order upwind scheme and backward time difference scheme. The successive over relaxation (SOR) method is used to solve all control volume pressure equations, and then the new time pressures are substituted in momentum equations to obtain the new time velocities. At last, pressures and velocities are substituted in mass and energy equations to obtain the new time gas volume fraction and temperatures. The constitutive relations include models for flow regimes, interfacial shear force, wall friction, wall heat transfer, and interphase heat transfer are in accordance with that used in RELAP5/MOD3.3.The PNCMC program uses the first-order upwind format to discretize the governing equations, which has good numerical stability. It has been verified by single-phase and two-phase natural circulation experiments, and the results show that the PNCMC program can be applied to thermal-hydraulic analysis under ocean conditions.For more information on PNCMC program development and validation, the reader is referred to the literature [18].


## 2.1. Governing Equations

The PNCMC program was developed to study thermal hydraulics under marine scenarios using a one-dimensional two-fluid flow model [19]. The two-phase flow equations are modified by adding extra body force terms.Gas momentum equation:(1)αgρg∂Vg∂t+αgρgVg∂Vg∂z=−αg∂P∂z+ΓgVi−Vg−pwαgwτgw+Mgd+αgρggz+fine,z.Liquid momentum equation:(2)αfρf∂Vf∂t+αfρfVf∂Vf∂z=−αf∂P∂z+ΓfVi−Vf−pwαfwτfw−Mgd+αfρfgz+fine,z.


## 2.2. Extra Body Force Terms

This formula gives fluid particle inertial forces,(3)fine=a0+ω̇×r+ω×ω×r+2ω×V.The first terma0 in the RHS accounts for the acceleration of the noninertial frame with respect to a stationary frame. The second term ω̇×r is the acceleration resulting from the noninertial frame’s angular acceleration. The third term ω×ω×r is the radial acceleration, which is highly sensitive to both rotational velocity and the fluid particle’s displacement from the rotational axis. The last element 2ω×V represents the Coriolis acceleration, which is independent of the fluid particle’s location but rather its velocity.Figure2 shows that a noninertial coordinate system OXYZ with unit vectors ijk and a stationary frame OXsYsZs are established. The origin O is positioned on the axis of rolling. In the initial stage, stationary frame OXsYsZs and the noninertial frame OXYZ are overlapping; the unit vector is straight up, which is in the opposite direction of gravity. Noninertial frame OXYZ has a rotational motion around an axis S=Sx,Sy,Sz and acceleration a0=a0,x,a0,y,a0,z. The unit vector n=cosθ,cosβ,cosγ symbolizes the positive mainstream of the channel in the OXYZ. Σ is an arbitrary channel cross-section whose points correspond to the domain x0±Δx,y0±Δy,z0±Δz.Figure 2
Flow channel and noninertial frame geometry.A ship’s rolling is analogous to a sine wave, which may be written as(4)ϕ=ϕmsinϖdt,(5)ω=ϕmϖdcosϖdt,(6)ω̇=−ϕmϖd2sinϖdt,where ϕ, ϕm, ϖd are the rolling angle, amplitude, and frequency, respectively.The reason why the pipeline cross-sectional scaleΔx,Δy,Δz is substantially less than the spatial scale x0,y0,z0 of the system is because the mainstream velocity is much higher than the cross-sectional velocity, fine,n can be interpreted as
(7)fine,n=a0,x−ω2SxSyy0−Sy2x0−Sz2x0+SxSzz0−ω̇Syz0+ω̇Szy0cosθ+a0,y−ω2SySzz0−Sz2y0−Sx2y0+SxSyx0−ω̇Szx0+ω̇Sxz0cosβ+a0,z−ω2SxSzx0−Sx2z0−Sy2z0+SySzy0−ω̇Sxy0+ω̇Syx0cosγ.fine,n is changed to fine,z in two-phase momentum equations. More detailed explanations can be found in Reference [18].During the rolling process,gz shifts depending on the angle between the channel and the vertical direction in a static frame OXsYsZs. Rolling is a variable-speed rotation around a fixed axis, the main stream’s new unit vector n′=cosθ′,cosβ′,cosγ′ in stationary frame OXsYsZs at the angle ϕ of rolling appears to result from the following operations:
(8)Aˆ=SxSxSxSySxSzSySxSySySySzSzSxSzSySzSz,(9)A∗=0−SzSySz0−Sx−SySx0,(10)I=100010001,(11)M=Aˆ+cosϕ·I−Aˆ+sinϕ·Aˆ,(12)n′=n·MT,(13)gz=g·cosγ′.where MT is the transpose matrix of M, and g is the gravity.


## 2.3. Solution and Validation

The two-phase six governing partial differential equations were discretized on a staggered grid by finite volume method with the semi-implicit, first-order upwind scheme and backward time difference scheme. The successive over relaxation (SOR) method is used to solve all control volume pressure equations, and then the new time pressures are substituted in momentum equations to obtain the new time velocities. At last, pressures and velocities are substituted in mass and energy equations to obtain the new time gas volume fraction and temperatures. The constitutive relations include models for flow regimes, interfacial shear force, wall friction, wall heat transfer, and interphase heat transfer are in accordance with that used in RELAP5/MOD3.3.The PNCMC program uses the first-order upwind format to discretize the governing equations, which has good numerical stability. It has been verified by single-phase and two-phase natural circulation experiments, and the results show that the PNCMC program can be applied to thermal-hydraulic analysis under ocean conditions.For more information on PNCMC program development and validation, the reader is referred to the literature [18].


## 3. Research Objects

In order to study the effect of ocean conditions on the flow characteristics and two-phase instability of the natural circulation system, a natural circulation system shown in Figure3 was established. The primary loop is a pipeline symmetrical double-loop structure, including three electric heating channels (HC1, HC2, and HC3), the double-pipe heat interchanger (HEX1 and HEX2), one nitrogen gas pressurizer, and connecting pipes. The height between the upper and lower plenum is 3.6 m, the horizontal distance between HEX1 and HEX 2 is 1.6 m, and the vertical distance between the centers of heating channels and heat exchangers is 1.6 m. The rolling axis is 3.0 m higher than the lower plenum. There are 16 stainless steel heating tubes arranged in 4×4 in the electric heating channel, the length of the heating tubes is 1000 mm, the diameter is 10 mm, the wall thickness is 1 mm, and the grid spacing is 13 mm. This structure can better simulate the symmetrical structure of the fluid system, and at the same time, it can also investigate the difference in the circulation characteristics of different branches under ocean conditions. The primary circuit system is installed on the motion platform, and the motion platform can perform inclination and simple harmonic rolling motions from front to back, left to right.Figure 3
Schematic diagram of the natural circulation system.As shown in Figure4, the constant pressure control body TDP0 is used to set the primary loop operating pressure, TDP1 is used to set the cooling water loop operating pressure, and TDV1 is used to set the temperature of the cooling water loop, and the alterable flow rate junction TDJ is used to set the flow rate of the cooling water. The flow rate of the cooling water changes with the power of the heating channels HC1, HC2, and HC3, and the goal is to maintain the inlet temperature of the heating channels unchanged. According to the measured data, the inlet resistance coefficient of the heating channel is kin=10.0, and the outlet resistance coefficient is kout=3.0. In this paper, the number of system nodes is 169, the maximum length of the control body is 0.4 meters, and the maximum time step is 0.01 s, which is the result of node independence verification.Figure 4
Node diagram of natural circulation for PNCMC.


## 4. The Effect Principle of Rolling Conditions on Natural Circulation

### 4.1. Extra Body Force Analysis

The additional force supplied by rolling can be used to identify the effect mechanism of the rolling conditions on the natural circulation. Figure5 is a schematic diagram of the added inertial force caused by the rolling. Tangential force, normal force, and Coriolis force are the most important extra inertial forces of fluid in a rolling motion. Because the Coriolis force is perpendicular to the velocity vector, its effect is neglected in the one-dimensional model.Figure 5
Schematic illustrating rolling’s induced inertia.According to Equation (14) of the tangential force Fe, The tangential force’s changing period is one rolling interval. The tangential force’s component forces along the natural circulation loop can produce two closed loops depending on the effective direction of the force. The inner loop is made up of HC1, HC3, and some of the lower connecting boxes, and the outer loop is made up of the heat exchanger, the downcomer, and the upper and lower connecting boxes.
(14)Fe=ρ×ω̇×r.According to Equation (15) of the normal force Fr, the normal force’s change period is half the rolling period. The component of normal acceleration along a left-right symmetrical flow channel is completely symmetrical.
(15)Fr=ρ×ω×ω×r.


### 4.2. Effect Principle of Rolling Conditions

In theory, for a closed loop, the rolling motion introduces an inertial force, and this force has the potential to cause the internal fluid to flow. Figure6 shows the oscillation of the flow behind the state that the loop fluid is at room temperature, the rolling period is 8 s, and the rolling amplitude is 30°. The fluid in the HC1 and HC3 changed from static to flowing, indicating that the flow in the inner loop oscillates sinusoidally, the oscillation period is the same as the rolling period. The oscillation amplitudes of HC1 and HC3 are equal, but the trends are opposite, and the phase difference is 180°; at the same time, the fluids of HEX1 and HEX2 heat exchangers also changed from static to flowing, indicating that the flow of the outer loop oscillates sinusoidally, the oscillation period and rolling period are the same, and their amplitudes are identical, the trends are opposite, and the phase difference is 180°; the fluid in the HC2 and the riser always remain stationary. This result shows that in the absence of heating, the rolling condition makes the fluid flow on both the inner and outer loops, but no fluid flows from the heating section to the heat exchanger. According to Equation (14) of the additional tangential force, for a system with a constant geometric shape, the magnitude of the tangential force can be represented by the magnitude of the angular acceleration. The greater the angular acceleration, the greater the tangential force, resulting in a greater oscillation of the inner and outer loop flow rate. For unheated conditions, the loop density distribution is uniform, and neither gravity nor normal force will make the fluid flow. It is theoretically and computationally explained that when fluids are impacted by a rolling motion, one of the mechanisms at effect is the introduced additional tangential force; it can make the fluids in the HC1 and HC3 heating sections and HEX1 and HEX2 heat exchangers oscillate sinusoidally, and the oscillation period is one rolling period.Figure 6
Flow rate oscillation under the unheated rolling condition.For the natural circulation prompted by density differences in the fluids, the flow is not only affected by the tangential force and the normal force, but also through the changes in the driving height difference generated by the rolling process’ inclination. Therefore, the flow rate oscillation form is determined by magnitudes of the driving height difference, the added normal force, and the extra tangential force. The frequencies of changes in these three quantities are shown in Table1.Table 1
The frequencies of changes in additional inertia force and driving head.
Additional tangential forceAdditional normal forceDriving height difference1 times the rolling frequency2 times the rolling frequency2 times the rolling frequency


## 4.1. Extra Body Force Analysis

The additional force supplied by rolling can be used to identify the effect mechanism of the rolling conditions on the natural circulation. Figure5 is a schematic diagram of the added inertial force caused by the rolling. Tangential force, normal force, and Coriolis force are the most important extra inertial forces of fluid in a rolling motion. Because the Coriolis force is perpendicular to the velocity vector, its effect is neglected in the one-dimensional model.Figure 5
Schematic illustrating rolling’s induced inertia.According to Equation (14) of the tangential force Fe, The tangential force’s changing period is one rolling interval. The tangential force’s component forces along the natural circulation loop can produce two closed loops depending on the effective direction of the force. The inner loop is made up of HC1, HC3, and some of the lower connecting boxes, and the outer loop is made up of the heat exchanger, the downcomer, and the upper and lower connecting boxes.
(14)Fe=ρ×ω̇×r.According to Equation (15) of the normal force Fr, the normal force’s change period is half the rolling period. The component of normal acceleration along a left-right symmetrical flow channel is completely symmetrical.
(15)Fr=ρ×ω×ω×r.


## 4.2. Effect Principle of Rolling Conditions

In theory, for a closed loop, the rolling motion introduces an inertial force, and this force has the potential to cause the internal fluid to flow. Figure6 shows the oscillation of the flow behind the state that the loop fluid is at room temperature, the rolling period is 8 s, and the rolling amplitude is 30°. The fluid in the HC1 and HC3 changed from static to flowing, indicating that the flow in the inner loop oscillates sinusoidally, the oscillation period is the same as the rolling period. The oscillation amplitudes of HC1 and HC3 are equal, but the trends are opposite, and the phase difference is 180°; at the same time, the fluids of HEX1 and HEX2 heat exchangers also changed from static to flowing, indicating that the flow of the outer loop oscillates sinusoidally, the oscillation period and rolling period are the same, and their amplitudes are identical, the trends are opposite, and the phase difference is 180°; the fluid in the HC2 and the riser always remain stationary. This result shows that in the absence of heating, the rolling condition makes the fluid flow on both the inner and outer loops, but no fluid flows from the heating section to the heat exchanger. According to Equation (14) of the additional tangential force, for a system with a constant geometric shape, the magnitude of the tangential force can be represented by the magnitude of the angular acceleration. The greater the angular acceleration, the greater the tangential force, resulting in a greater oscillation of the inner and outer loop flow rate. For unheated conditions, the loop density distribution is uniform, and neither gravity nor normal force will make the fluid flow. It is theoretically and computationally explained that when fluids are impacted by a rolling motion, one of the mechanisms at effect is the introduced additional tangential force; it can make the fluids in the HC1 and HC3 heating sections and HEX1 and HEX2 heat exchangers oscillate sinusoidally, and the oscillation period is one rolling period.Figure 6
Flow rate oscillation under the unheated rolling condition.For the natural circulation prompted by density differences in the fluids, the flow is not only affected by the tangential force and the normal force, but also through the changes in the driving height difference generated by the rolling process’ inclination. Therefore, the flow rate oscillation form is determined by magnitudes of the driving height difference, the added normal force, and the extra tangential force. The frequencies of changes in these three quantities are shown in Table1.Table 1
The frequencies of changes in additional inertia force and driving head.
Additional tangential forceAdditional normal forceDriving height difference1 times the rolling frequency2 times the rolling frequency2 times the rolling frequency


## 5. Characteristics of Two-Phase Natural Circulation in a Rolling Environment

For the parallel three heating channels structure, taking the steady state point of 1% of the outlet dryness of a single heating channel as the starting point of transient calculation, working conditions such as the rolling periods of 3 s, 8 s, 13 s, 18 s, and 23 s, and the rolling amplitudes of 5°, 15°, 22.5°, 30°, 45°, etc. were carried out.Figure7 shows the variation of outlet dryness and mass flow rate under different amplitude conditions with a rolling period of 13 s. Since the system has formed stable oscillation, only the flow oscillation of the HC1 heating section and HC1 heat exchanger is given here. It is clear that the outlet dryness and flow rate oscillation amplitudes are positively associated with the rolling amplitude. The flow rate on the primary side of the heat exchanger deviates farther from sinusoidal waveform as the rolling amplitude rises. If the rolling amplitude is 45°, the heat exchanger flow rate has a secondary disturbance in the process of its increase, which is no longer a simple harmonic sinusoidal fluctuation.Figure 7
Flow and dryness fluctuations with a 13 s period rolling condition.Figure8 demonstrates the results of the Fourier spectrum analysis performed on the flow rate at the primary side of the heating section and the heat exchanger. There are two frequency oscillations in the flow rates on both the primary side of the heating section and the heat exchanger, and the corresponding periods are 6.5 s and 13 s, respectively, namely, 1/2 rolling period and 1 rolling period. The frequency of 0.154 Hz controls the oscillation of flow rate inside the heat section, the additional normal force and the change of the driving pressure head play a leading role, and the tangential force plays a secondary role. The frequency of 0.077 Hz is dominant in the flow rate oscillation of the primary side of the heat exchanger, and the tangential force plays a dominant role.Figure 8
FFT spectrum analysis of mass flow.Figure9 shows the variation of outlet dryness and mass flow rate under different period conditions with a rolling amplitude of 22.5°. It can be found from the figure that the oscillation amplitude of the outlet dryness and mass flow rate is inversely proportional to the rolling period. When the period is 3 s, the larger rolling angular acceleration introduces a larger tangential force resulting in the flow oscillation of the heating section dominated by the tangential force. Hence, the flow oscillation period of the HC1 heating section is one rolling period.Figure 9
Flow rate and dryness oscillations with an amplitude of 22.5° rolling condition.


## 6. Density Wave Instability in Static Case

### 6.1. Instability Types in Natural Circulation Systems

Numerous experimental and theoretical studies have proved that flashing induced instability [20–31] and density-wave instability [32–36] usually occur in natural circulation systems.Instabilities caused by flashing are common in natural circulation systems that feature a long adiabatic chimney or riser. A drop in static head along the chimney causes a rapid boiling at a specific place where the liquid temperature exceeds the local saturation temperature. And then, because of the increased vapor volume, the natural circulation flow rate has risen. Its temperature just at the chimney’s input falls as a result of the increased flow. Until after chimney has been loaded with subcooled liquid, the flow rate reduces till the liquid boils again suddenly.Fukuda and Kobori’s [33] natural circulation experiment revealed two fundamental types of oscillation. One mode is the so-called U-tube oscillation, which is marked by channel flow rate fluctuations with a 180° phase difference. The fluctuation of channel flow rates within phase with one another is another mode. The latter scenario describes that the overall flow rate varies in phase with the fluctuations in the channel flow rate. Also, the two extreme phases of flow oscillation—the “U-tube” oscillation and the total flow oscillation—sometimes coexist in their natural circulation tests. According to the experimental results, Fukuda proposed two types of density wave instabilities. Type-I, caused by gravity, is seen in upward-vertical systems that feature a lengthy adiabatic chimney at low pressure and low quality. Any disturbance can significantly alter the gravity head, void fraction, and flow conditions at low pressure and quality. Consequently, oscillating situations might result from the feedback between flow, void fraction, and gravitational head. This process is crucial for natural circulation loops. Figure 10 shows Fukuda’s experimental facility and type I oscillation under natural circulation. The study object in this paper is analogous to the experimental setup shown in Figure 11. The difference lies in the three heating channels and the symmetrically arranged heat exchanger.Figure 10
Experimental device and type-I density wave instability oscillations [33].Figure 11
Density wave oscillation under static conditions.As for type-II, due to friction, the primary reason for this is because flow disturbances travel at different speeds in the single-phase and two-phase regions. Pressure-drop changes occur when the flow or void fraction in the two-phase area changes. Due to the slow propagation of the perturbation along the two-phase region, the initial perturbations in the two-phase region are marked by a large delay. And type-II density wave instability usually occurs in high pressure conditions.


### 6.2. Thermohydraulic Oscillations of Type-I Density Wave Oscillation Instability

Figure12 illustrates the variation in flow rate with heating power in system pressure of 1.5 MPa and inlet subcooling of 40 K.Figure 12
Vapor void fraction and liquid temperature vary with time.
(a)(b)(c)(d)Figure12 shows the void fraction of vapor and temperature of the liquid at the outlet of the heating section, this demonstrates that, firstly, subcooled boiling appears in the heating channel. The liquid subcooling degree at the outlet of the heating section at the beginning of boiling is 4.11 K, and the vapor entering the riser is condensed by the liquid. The period of flow oscillation given in Figure 1 is short, which is different from the mechanism and the phenomenon of flashing-induced instability, so it can be judged that the flow rate oscillation is caused by instability of density wave.The process for flow rate oscillation changing is described in detail as follows. When the single-channel heating power reaches 72 kW, the density wave oscillation occurs between the heating channels. As shown in Figure1(c), the out-of-phase oscillations between the channels have a phase difference of 120° and an oscillation period of 3.6 s. The void fraction of vapor at the outlet of the heating section reaches 0.122. Due to the asymmetry of the oscillating waveform, the out-of-phase oscillation of the heating channels caused a small amplitude oscillation of the flow in the riser. When the heating power is further increased, as shown in Figures 1(b) and 2(c), if the void fraction at the riser exit is larger than zero, the flow oscillation amplitude of the riser increases significantly. As shown in Figure 1(d), the flow rate in the riser is approximately sinusoidal, with an oscillation period of about 5.8 s; the flow oscillation of the heating channel becomes more complex, and a compound oscillation of the system density wave oscillation and the interchannel density wave oscillation occurs, while the interchannel density wave oscillation plays a dominant role. The flow oscillation between the heating channels is in a competitive relationship with the system flow oscillation, which further increases the heating power. As stated in Figure 1(e), the feedback between the pressure difference-flow rate-void in the riser is much larger than that between the three heating channels. The system flow oscillation is enhanced, and the out-of-phase oscillation between the heating channels is suppressed, so that the flow oscillation in heating channel obeys the density wave oscillation of system, and its oscillation amplitude is much smaller than the interchannel out-of-phase oscillation shown in Figure 1(c). As shown in Figure 1(f), the compound oscillation of the system density wave oscillation and the interchannel density wave oscillation occurs again. Figure 1(h) shows the FFT analysis results of the HC1 flow rate from 6500 s to 6700 s. The flow oscillation of HC1 is dominated by two frequencies, the system density wave oscillation (0.176 Hz, the period is 5.68 s, and the amplitude is 0.11 kg/s); the interchannel density wave oscillation (0.317 Hz, the period is 3.15 s, and the amplitude is 0.1 kg/s), respectively. When the power of the heating channel reaches 110 kW, the density wave oscillation of the system disappears, and pure out-of-phase oscillation with a phase difference of 120° left in the heating channel. When the power is 117 kW, the system becomes a steady state and the type-I density wave oscillation completely disappears.


## 6.1. Instability Types in Natural Circulation Systems

Numerous experimental and theoretical studies have proved that flashing induced instability [20–31] and density-wave instability [32–36] usually occur in natural circulation systems.Instabilities caused by flashing are common in natural circulation systems that feature a long adiabatic chimney or riser. A drop in static head along the chimney causes a rapid boiling at a specific place where the liquid temperature exceeds the local saturation temperature. And then, because of the increased vapor volume, the natural circulation flow rate has risen. Its temperature just at the chimney’s input falls as a result of the increased flow. Until after chimney has been loaded with subcooled liquid, the flow rate reduces till the liquid boils again suddenly.Fukuda and Kobori’s [33] natural circulation experiment revealed two fundamental types of oscillation. One mode is the so-called U-tube oscillation, which is marked by channel flow rate fluctuations with a 180° phase difference. The fluctuation of channel flow rates within phase with one another is another mode. The latter scenario describes that the overall flow rate varies in phase with the fluctuations in the channel flow rate. Also, the two extreme phases of flow oscillation—the “U-tube” oscillation and the total flow oscillation—sometimes coexist in their natural circulation tests. According to the experimental results, Fukuda proposed two types of density wave instabilities. Type-I, caused by gravity, is seen in upward-vertical systems that feature a lengthy adiabatic chimney at low pressure and low quality. Any disturbance can significantly alter the gravity head, void fraction, and flow conditions at low pressure and quality. Consequently, oscillating situations might result from the feedback between flow, void fraction, and gravitational head. This process is crucial for natural circulation loops. Figure 10 shows Fukuda’s experimental facility and type I oscillation under natural circulation. The study object in this paper is analogous to the experimental setup shown in Figure 11. The difference lies in the three heating channels and the symmetrically arranged heat exchanger.Figure 10
Experimental device and type-I density wave instability oscillations [33].Figure 11
Density wave oscillation under static conditions.As for type-II, due to friction, the primary reason for this is because flow disturbances travel at different speeds in the single-phase and two-phase regions. Pressure-drop changes occur when the flow or void fraction in the two-phase area changes. Due to the slow propagation of the perturbation along the two-phase region, the initial perturbations in the two-phase region are marked by a large delay. And type-II density wave instability usually occurs in high pressure conditions.


## 6.2. Thermohydraulic Oscillations of Type-I Density Wave Oscillation Instability

Figure12 illustrates the variation in flow rate with heating power in system pressure of 1.5 MPa and inlet subcooling of 40 K.Figure 12
Vapor void fraction and liquid temperature vary with time.
(a)(b)(c)(d)Figure12 shows the void fraction of vapor and temperature of the liquid at the outlet of the heating section, this demonstrates that, firstly, subcooled boiling appears in the heating channel. The liquid subcooling degree at the outlet of the heating section at the beginning of boiling is 4.11 K, and the vapor entering the riser is condensed by the liquid. The period of flow oscillation given in Figure 1 is short, which is different from the mechanism and the phenomenon of flashing-induced instability, so it can be judged that the flow rate oscillation is caused by instability of density wave.The process for flow rate oscillation changing is described in detail as follows. When the single-channel heating power reaches 72 kW, the density wave oscillation occurs between the heating channels. As shown in Figure1(c), the out-of-phase oscillations between the channels have a phase difference of 120° and an oscillation period of 3.6 s. The void fraction of vapor at the outlet of the heating section reaches 0.122. Due to the asymmetry of the oscillating waveform, the out-of-phase oscillation of the heating channels caused a small amplitude oscillation of the flow in the riser. When the heating power is further increased, as shown in Figures 1(b) and 2(c), if the void fraction at the riser exit is larger than zero, the flow oscillation amplitude of the riser increases significantly. As shown in Figure 1(d), the flow rate in the riser is approximately sinusoidal, with an oscillation period of about 5.8 s; the flow oscillation of the heating channel becomes more complex, and a compound oscillation of the system density wave oscillation and the interchannel density wave oscillation occurs, while the interchannel density wave oscillation plays a dominant role. The flow oscillation between the heating channels is in a competitive relationship with the system flow oscillation, which further increases the heating power. As stated in Figure 1(e), the feedback between the pressure difference-flow rate-void in the riser is much larger than that between the three heating channels. The system flow oscillation is enhanced, and the out-of-phase oscillation between the heating channels is suppressed, so that the flow oscillation in heating channel obeys the density wave oscillation of system, and its oscillation amplitude is much smaller than the interchannel out-of-phase oscillation shown in Figure 1(c). As shown in Figure 1(f), the compound oscillation of the system density wave oscillation and the interchannel density wave oscillation occurs again. Figure 1(h) shows the FFT analysis results of the HC1 flow rate from 6500 s to 6700 s. The flow oscillation of HC1 is dominated by two frequencies, the system density wave oscillation (0.176 Hz, the period is 5.68 s, and the amplitude is 0.11 kg/s); the interchannel density wave oscillation (0.317 Hz, the period is 3.15 s, and the amplitude is 0.1 kg/s), respectively. When the power of the heating channel reaches 110 kW, the density wave oscillation of the system disappears, and pure out-of-phase oscillation with a phase difference of 120° left in the heating channel. When the power is 117 kW, the system becomes a steady state and the type-I density wave oscillation completely disappears.


## 7. Instability of Density Waves in Rolling Situations

### 7.1. Compound Oscillations of Thermal Hydraulics in Rolling Situations

Rolling motion is equivalent to an external excitation force, forcing the natural circulation flow to oscillate periodically and approximately sinusoidally. The void fraction-pressure drop-flow rate feedback is equivalent to the self-excitation effect inside the system, forcing the flow rate to oscillate according to the feedback transfer period. The oscillation of the flow generated by rolling motion will likely increase or decrease the transfer rate of the internal void fraction-pressure drop-flow rate feedback, while the reduction of the gravitational pressure drop caused by the rolling will also weaken the feedback strength of the void fraction-gravity pressure drop-flow rate. So the final oscillation law of natural circulation flow rate depends on the competitive relationship between density wave instability oscillation and rolling flow oscillation. When density wave oscillation and rolling oscillation coexist, the flow oscillation law will be very complex.Figure13 shows the oscillation of flow during the 13 s rolling period, the amplitude of 22.5°, the inlet subcooling degree of 40.0 K, increasing the heating power continuously and slowly. Since rolling also causes flow rate oscillations, we can determine what kind of oscillations occurs from the magnitude and waveform of the flow oscillations. Obviously, during the whole heating process, oscillations caused by rolling (Figure 3(c)), out-of-phase interchannel oscillations (Figure 3(d)), compound oscillations between the system and channels (Figure 3(e)), and interchannel oscillations, out-of-phase oscillations (Figure 3(f)), as well as oscillations caused by simply rolling and transitions of different oscillations appeared successively.Figure 13
Flow rate oscillation under rolling conditions (period 13 s, amplitude 22.5°).


### 7.2. Decomposition of Compound Oscillations

As shown in Figure14, the changes in flow rate during compound oscillations are highly complex and become unpredictable. Therefore, we performed an FFT analysis on the flow rate (such as Figure 14). From the FFT results, it can be found that the flow oscillation of the heating section is mainly dominated by the interchannel oscillation of 0.294 Hz, superimposed the rolling oscillation of 0.154 Hz and the system density wave oscillation of 0.171 Hz as well as high-frequency small-amplitude oscillations. The flow oscillation in the riser is mainly composed of the superposition of the system density wave oscillation of 0.171 Hz and the rolling oscillation of 0.154 Hz, and their oscillation amplitudes are comparable in size. As we know, when two simple harmonics with similar frequencies and amplitudes are superimposed, the phenomenon of “beat” will occur, that is, the amplitude of Figure 3(a) between 5500 and 7000 shows a crescendo-decrescendo phenomenon.Figure 14
FFT analysis of heating channel flow rate and riser flow rate.It is seen from Figure7 that the change of the flow rate by 13 s and 22.5° rolling motion is very small. Therefore, compared with the vertical static condition (Figure 11), although there are differences in the details of the oscillation waveforms between them, the overall oscillation law and the initial power of interchannel oscillation, the initial power of system oscillation, the terminated power of system oscillation, and the terminated power of interchannel oscillation are consistent.


### 7.3. Effect of Rolling Amplitude on Compound Oscillation

As shown in Figure15, keeping the rolling period unchanged, if the rolling angle is increased to 30°, the flow oscillation amplitude caused by rolling increases compared with the 22.5° case, and the flow oscillation law changes accordingly. In the range of 80~94 kW, the compound oscillation between the heating interchannels and the system occurred. Compared with the 22.5° rolling motion, the initial instability power is significantly increased. But after 94 kW, the oscillation between the heating channels disappears, and flow rates in the heating channel and the riser are synchronized to approximate simple harmonic oscillation, and the period of oscillation is equal to one-half of the rolling period. In Figures 5(a) and 5(b), it can be found that the oscillation amplitude is not continuous. In the A region, that is, if the heating power is 125 kW, the amplitude of flow oscillation changes suddenly. It is shown that although the oscillation period is exactly the same as the rolling period, in the range of 94 kW-125 kW, the void fraction-gravity pressure drop-flow rate feedback inside the system plays a strengthening role in the flow rate oscillation and increases the flow rate oscillation amplitude.Figure 15
Flow rate oscillation under rolling conditions of 13 s period and 30° angle.As shown in Figure16, keeping the rolling period unchanged and increasing the rolling angle to 45°, in the whole process of increasing the power, there is neither interchannel density wave oscillation nor system density wave oscillation, and the flow rate oscillation is synchronized with the rolling motion. Similarly, we found a sudden change in the oscillation amplitude in regions B and C, that is, in the range of 65 kW-130 kW, the void fraction-gravity pressure drop-flow rate feedback inside the system played a strengthening role in the flow rate oscillation.Figure 16
Flow rate oscillation under rolling conditions with a period of 13 s and an angle of 45°.In general, from the perspective of the flow oscillation law, whether the density wave’s instability arises under the rolling situation in comparison to the static state, the flow oscillation amplitude does not change significantly.


## 7.1. Compound Oscillations of Thermal Hydraulics in Rolling Situations

Rolling motion is equivalent to an external excitation force, forcing the natural circulation flow to oscillate periodically and approximately sinusoidally. The void fraction-pressure drop-flow rate feedback is equivalent to the self-excitation effect inside the system, forcing the flow rate to oscillate according to the feedback transfer period. The oscillation of the flow generated by rolling motion will likely increase or decrease the transfer rate of the internal void fraction-pressure drop-flow rate feedback, while the reduction of the gravitational pressure drop caused by the rolling will also weaken the feedback strength of the void fraction-gravity pressure drop-flow rate. So the final oscillation law of natural circulation flow rate depends on the competitive relationship between density wave instability oscillation and rolling flow oscillation. When density wave oscillation and rolling oscillation coexist, the flow oscillation law will be very complex.Figure13 shows the oscillation of flow during the 13 s rolling period, the amplitude of 22.5°, the inlet subcooling degree of 40.0 K, increasing the heating power continuously and slowly. Since rolling also causes flow rate oscillations, we can determine what kind of oscillations occurs from the magnitude and waveform of the flow oscillations. Obviously, during the whole heating process, oscillations caused by rolling (Figure 3(c)), out-of-phase interchannel oscillations (Figure 3(d)), compound oscillations between the system and channels (Figure 3(e)), and interchannel oscillations, out-of-phase oscillations (Figure 3(f)), as well as oscillations caused by simply rolling and transitions of different oscillations appeared successively.Figure 13
Flow rate oscillation under rolling conditions (period 13 s, amplitude 22.5°).


## 7.2. Decomposition of Compound Oscillations

As shown in Figure14, the changes in flow rate during compound oscillations are highly complex and become unpredictable. Therefore, we performed an FFT analysis on the flow rate (such as Figure 14). From the FFT results, it can be found that the flow oscillation of the heating section is mainly dominated by the interchannel oscillation of 0.294 Hz, superimposed the rolling oscillation of 0.154 Hz and the system density wave oscillation of 0.171 Hz as well as high-frequency small-amplitude oscillations. The flow oscillation in the riser is mainly composed of the superposition of the system density wave oscillation of 0.171 Hz and the rolling oscillation of 0.154 Hz, and their oscillation amplitudes are comparable in size. As we know, when two simple harmonics with similar frequencies and amplitudes are superimposed, the phenomenon of “beat” will occur, that is, the amplitude of Figure 3(a) between 5500 and 7000 shows a crescendo-decrescendo phenomenon.Figure 14
FFT analysis of heating channel flow rate and riser flow rate.It is seen from Figure7 that the change of the flow rate by 13 s and 22.5° rolling motion is very small. Therefore, compared with the vertical static condition (Figure 11), although there are differences in the details of the oscillation waveforms between them, the overall oscillation law and the initial power of interchannel oscillation, the initial power of system oscillation, the terminated power of system oscillation, and the terminated power of interchannel oscillation are consistent.


## 7.3. Effect of Rolling Amplitude on Compound Oscillation

As shown in Figure15, keeping the rolling period unchanged, if the rolling angle is increased to 30°, the flow oscillation amplitude caused by rolling increases compared with the 22.5° case, and the flow oscillation law changes accordingly. In the range of 80~94 kW, the compound oscillation between the heating interchannels and the system occurred. Compared with the 22.5° rolling motion, the initial instability power is significantly increased. But after 94 kW, the oscillation between the heating channels disappears, and flow rates in the heating channel and the riser are synchronized to approximate simple harmonic oscillation, and the period of oscillation is equal to one-half of the rolling period. In Figures 5(a) and 5(b), it can be found that the oscillation amplitude is not continuous. In the A region, that is, if the heating power is 125 kW, the amplitude of flow oscillation changes suddenly. It is shown that although the oscillation period is exactly the same as the rolling period, in the range of 94 kW-125 kW, the void fraction-gravity pressure drop-flow rate feedback inside the system plays a strengthening role in the flow rate oscillation and increases the flow rate oscillation amplitude.Figure 15
Flow rate oscillation under rolling conditions of 13 s period and 30° angle.As shown in Figure16, keeping the rolling period unchanged and increasing the rolling angle to 45°, in the whole process of increasing the power, there is neither interchannel density wave oscillation nor system density wave oscillation, and the flow rate oscillation is synchronized with the rolling motion. Similarly, we found a sudden change in the oscillation amplitude in regions B and C, that is, in the range of 65 kW-130 kW, the void fraction-gravity pressure drop-flow rate feedback inside the system played a strengthening role in the flow rate oscillation.Figure 16
Flow rate oscillation under rolling conditions with a period of 13 s and an angle of 45°.In general, from the perspective of the flow oscillation law, whether the density wave’s instability arises under the rolling situation in comparison to the static state, the flow oscillation amplitude does not change significantly.


## 8. Conclusion

In the ocean environment, sea waves can change the thermal hydraulic characteristics of natural circulation thermal fluid systems by introducing extra body force on liquid. This is an important issue for the system safety. In order to study the thermal hydraulic characteristics under ocean conditions, PNCMC was developed by adding extra body force in one-dimensional two-phase’s momentum equations. In this paper, flow rate variation and type-I density wave instability of natural circulation system with three parallel heating channels under rolling conditions was studied by PNCMC. The mechanism of rolling motion and its impact on flow oscillation and unstable boundary power are investigated.The effect of rolling motion on the natural circulation system can be mainly described as three aspects, namely the additional tangential force, the additional normal force and the gravity pressure drop’s variation caused by the loop’s inclination. The integral of additional tangential force along the loop is simple harmonic. This additional tangential force directly drives the fluid flow. The additional normal force is equivalent to exerting a harmonically varying volume force on the fluid along the flow direction, while the gravitational height difference is directly affected by inclination, thereby reducing gravity-driven head. The magnitude and direction of the additional tangential force, the additional normal force, and the gravity head, depending on the spatial arrangement of the system and the position of the rolling axis.The high surface heat flux density of the heating section generates supercooled boiling in a low-pressure natural circulation system with several parallel heating channels, and when the void fraction increases to a certain value, the density wave oscillation between parallel channels occurs. The vapor-liquid mixture entering the riser from the heating section is gradually condensed in the riser due to the supercooled liquid phase. When enough vapor enters the riser, and the vapor occupies a long enough distance in the riser, the gravitational pressure drop is significantly changed, so the system density wave oscillation occurs under the feedback between the pressure drop-flow rate-vapor fraction. The inter-channel oscillation and the system density wave oscillation are overlapped to produce a much more complicated compound oscillation.If rolling motion has little change to the natural circulation flow rate, when instability of density wave occurs, the flow rate oscillation rule, the initial power of inter-channel oscillation, the initial power of system oscillation, and the terminated power of system oscillation as well as the terminated power of inter-channel oscillation, basically are the same as the vertical condition. When the rolling angle is large, the large amplitude flow oscillation resulting from the rolling will suppress the occurrence of density wave instability and force the flow oscillation to obey the law caused by rolling. Within a certain heating power range, the void-gravity pressure drop-flow rate self-feedback inside system will enhance the flow rate’s oscillation amplitude.



---

*Source: 3249806-2023-02-03.xml*