Non-Linear Qualitative Dynamic Analysis of Supercritical Water-Heated Channels under External Vertical Accelerations

Abstract

Featured Application

The dimensionless nonlinear model can be modified to the qualitative dynamic analysis of the heated systems with different supercritical fluids under external vertical accelerations.

Abstract

Employing the external force method to regard seismic impact and the three-region methodology to analyze the supercritical heated channel, a non-linear dynamic model was developed to investigate the transient characteristics of single channel or parallel channels under the impacts of vertical sinusoidal and seismic accelerations. The present model was validated against the experimental data, which could suitably estimate the additional pressure drop caused by the vertical vibrations. The influences of parameters on the seismic-induced oscillation conducted in a supercritical heated channel indicated that a longer heated length, uprating operation power and a larger outlet loss coefficient all exhibit unstable effects, while the increase of inlet loss coefficient, a larger tube diameter and a lower inlet fluid temperature would tend to stabilize the system. Moreover, the supercritical fluid would present a high natural frequency in the very small N_P-SUB region. The parametric effects on the parallel channel system are related to the inherent stability nature of initial state and the interactions among channels. The more uneven heat flux distribution among channels would cause a larger vibration-induced oscillation. In particular, when it is combined with the resonance effect, the system may exhibit much larger oscillations than in the case of non-resonance.

1. Introduction

The applications of supercritical fluids with their peculiar characteristics have been considered in many technical fields, such as renewable energy system [1], refrigeration system [2], next generation nuclear technology [3] and commercial ultra-supercritical coal power plant. In addition, a supercritical water nuclear reactor (SCWR) with the merits of comparatively high thermal performance, compact and inherent safety design, and economic competitiveness is listed as a promising type in the international development program of the Generation-IV advanced nuclear reactor [4]. At a supercritical pressure, the heating fluids exhibit dramatic changes in fluid properties near their critical or pseudocritical points. These features can substantially affect the thermal characteristics and stabilities of the systems with supercritical fluids, which should be addressed in the development and applications of such systems.

At supercritical pressures, flow boiling phenomenon will not appear during the heating process. For a supercritical heated channel, the fluid behaves in a liquid-like manner below the corresponding pseudocritical point; however, it presents gas-like behavior when the fluid temperature is beyond the pseudocritical value [5]. Therefore, some researchers suggested three-region methodology to separate the supercritical flow into heavy fluid region, heavy and light fluid mixture region, and light fluid region, while considering the investigations in the supercritical heating systems [6,7].

Upon the drastic variation of flow properties in the pseudocritical region, the supercritical fluid systems could be subject to various kinds of instability, particularly the typical type of density-wave instability [8,9]. Since the instability problems may lead to the fluctuations of working fluids and system components, and thus affect the operation safety of a supercritical heated system, much more attentions need to be paid on the stability issues. In the literature, a linear model in the frequency domain and a non-linear model in the time domain are two major categories for carrying out theoretical studies of stability. Linear analysis methodology is a familiar tool to explore the stability threshold of a supercritical heated system due to its simplicity [8,9]. However, such systems are inherently non-linear, which cannot be explored by the linear methodology. This suggests the need of a proper non-linear model to clarify the instability phenomena of supercritical heated systems.

Several researchers had predicted non-linear characteristics of supercritical heated systems by the tools of computational fluid dynamics (CFD) [10] or the modifications of the existing large system codes from two-phase flows [11,12]. However, because their complexity and time-consumption, these tools are, in general, not favorable for the detailed investigation of non-linear phenomena. Thus, our recent studies developed simple non-linear dynamic models to analyze the stability problems of supercritical water uniformly heated in a single channel [13] and parallel channels [14], respectively, particularly for that after the occurrence of instability. The results indicated that complex and abundant non-linear phenomena, such as supercritical Hopf and period-doubled bifurcations, and their evolutions to chaotic oscillations, could appear in a uniformly-heated channel with supercritical water [13]. On the other hand, parallel supercritical heated channels essentially presented out-of-phase mode of oscillations if total flow rates were kept constant. Channel-to-channel interactions due to the more difference among radial power distributions would have unstable effects on these systems [14]. Such out-of-phase instabilities have also been identified by experimental studies in parallel two channel systems with supercritical water [15,16]. Moreover, the parallel supercritical heated channel systems could be stabilized by the increases in inlet flow resistance or channel diameter, while the unstable influences on the system stability could be attributed to the increases in the outlet flow resistance or channel length [14].

The seismic issues are very crucial for the safe operations of various power generation systems. However, the relevant studies of seismic impacts on the flow stabilities were very limited in the literatures, which majorly conducted in nuclear boiling systems [17,18,19,20,21]. Among these studies, two major approaches, i.e., external force method [17,18,19,20,21] and moving grid method [21], were employed to evaluate the seismic impact on two phase flows. The external force method is usually considered in the rigid systems, where no deformation occurs in the system structure during the seismic motions; otherwise, it would require the moving grid method instead. The previous study indicated that the impact of seismic vibration on two-phase boiling flows may depend on its intensity and the inherent stability natures of operating states. The resonance oscillation can be triggered by the seismic wave if its oscillation frequency is equivalent to the natural frequency of working flow [17]. In a parallel nuclear-coupled boiling channel system, the external vibration acceleration could result in the resonance phenomena with either regional or core-wide instabilities. The impact of vertical acceleration would generate a much stronger effect on the core power oscillations with respect to those in the other directions [18]. Our recent studies also illustrated that complicated seismic-induced oscillations could exist in a nuclear-coupled boiling channel system [19]. For a natural circulation loop involving parallel nuclear-coupled boiling channels, the resonance effects in different natural circulation stable states would dominate seismic-induced flow oscillations. Under the impact of vertical seismic motion, the parallel channels tended to exhibit in-phase mode of oscillations. A more significant resonance effect could be triggered by the natural circulation state with a higher core inlet subcooling [20].

To the authors’ best knowledge, at present, the research regarding the influence of external vibration, such as seismic acceleration, on either single channel or parallel channels of the supercritical heated system is very sparse. The unique nature of parallel supercritical heated channels exhibits the complex interactions distributing among them, in which the unstable phenomenon is commonly the compound type of density-wave instability merged with parallel channel instability. The seismic impact interacting with the distinctive stability characteristics of such systems should need to be explored. Thus, the major objective of the present study is aimed at developing a non-linear dynamic model to deal with the effects of external vertical acceleration and seismic acceleration on the transient behaviors of parallel uniformly-heated channels with supercritical water. The non-linear model can also be applied to examine the seismic effects on the single supercritical heated channel after the proper simplifications.

2. The Model and Solution Method

Based on our preliminary work [22] employing the external force method to evaluate seismic impact, the non-linear dynamic model of parallel supercritical water uniformly-heated channels in the static system developed previously by the authors [14] is extended to the vertical vibration system. The modelling structure involving M parallel uniformly-heated channels subject to vertical seismic acceleration at a supercritical pressure is displayed in Figure 1. Adopting the three-region methodology suggested by Zhang et al. [7], each supercritical water-heated channel is separated into region 1 of heavy fluid, region 2 of heavy and light fluid mixture and region 3 of light fluid, as shown in Figure 1. At a supercritical pressure of 25 MPa, the dividing point “A” between region 1 and 2 is about 349 °C and “B” between region 2 and 3 is nearly 449 °C, respectively.

Considering the j-th channel in the parallel supercritical heated channel system in Figure 1, the analytical problem can be simplified on the basis of following assumptions:

• The whole system is assumed as a rigid body;

• The system is retained at a supercritical pressure of 25 MPa under both static and external vibration conditions, while the pressure drop through each channel is usually far less than the system pressure;

• Each channel is supposed to have a uniform heat flux distribution along the axial direction;

• The same constant inlet enthalpy is distributed over all channels under both static and external vibration conditions;

• The parallel channels possess the same pressure drop under both static and external vibration conditions, as they are connected by the common upper plenum and lower plenum.

For future applications to different supercritical fluids besides supercritical water, the state variables for each channel are non-dimensionalized based on the fluid properties corresponding to the pseudocritical point at a supercritical pressure, as defined in the nomenclature. When considering the effect of seismic motion on the present system, the entire system is supposed as a rigid body, which infers no deformation occurring in system components during seismic motions. This demonstrates that all the parts are subject to the same but variant vibration acceleration, inferring they face the external force induced by the seismic motions. Thus, based on the assumptions and the external force method, with the non-dimensional methodology [8], the one-dimensional dimensionless conservation equations for the j-th supercritical uniformly-heated channel subject to vertical seismic acceleration can be written as:

(1)$\frac{\partial {\rho }_j^{+}}{\partial t^{+}}+\frac{\partial {\rho }_j^{+}{u}_j^{+}}{\partial z^{+}}=0$

(2)$\frac{\partial {\rho }_j^{+}{h}_j^{+}}{\partial t^{+}}+\frac{\partial {\rho }_j^{+}{h}_j^{+}{u}_j^{+}}{\partial z^{+}}=N_{P-PCH,j}$

(3)$\frac{\partial {\rho }_j^{+}{u}_j^{+}}{\partial t^{+}}+\frac{\partial {\rho }_j^{+}{u}_j^{+}{}^2}{\partial z^{+}}=-{\mathsf{\Lambda }}_j{\rho }_j^{+}{u}_j^{+}{}^2-{\displaystyle \sum _{m=1}^N{k}_{m,j}\delta \left(z_j^{+}-z_{m,j}^{+}\right)\frac{{\rho }_j^{+}{u}_j^{+}{}^2}{2}}-\frac{{\rho }_j^{+}{a}_V^{+}}{Fr}-\frac{\partial P_j^{+}}{\partial z^{+}}$

where N_P-PCH,j is pseudo phase-change number for the j-th channel that was originally suggested by Ambrosini and Sharabi [8], as defined in the nomenclature. The parameter $a_V^{+}$ represents the dimensionless vertical acceleration including both the gravity and vertical seismic ($a_{ses}^{+}$) accelerations. Thus,

(4)$a_V^{+}=\{\begin{array}{c}1+a_{ses}^{+},underverticalvibrations\\ 1,understaticconditions\end{array}$

The parallel heated channels shown in Figure 1 exhibit the same dynamic pressure drop under seismic conditions as they are linked with the common upper plenum and lower plenum.

(5)$\Delta P_{ch,1}^{+}=\Delta P_{ch,j}^{+}=\cdots =\Delta P_{ch,M}^{+}$

The dynamic pressure drop for each channel subject to vertical seismic acceleration could be estimated by the summation of the dynamic pressure drop in the static condition plus the additional pressure drop induced by the seismic vibration:

(6)$\Delta P_{ch,j}^{+}=\Delta P_{A,j}^{+}+\Delta P_{B,j}^{+}+\Delta P_{C,j}^{+}+\Delta P_{ses,j}^{+}=Coe{f}_j\frac{d{u}_{i,j}^{+}}{d{t}^{+}}+\Delta P_{ch0,j}^{+}+\Delta P_{ses,j}^{+},j=1,2,\cdots ,M$

(7)$Coe{f}_j=C_{A,j}{\lambda }_{A,j}^{+}+C_{B,j}({\lambda }_{B,j}^{+}-{\lambda }_{A,j}^{+})+C_{C,j}(1-{\lambda }_{B,j}^{+})$

where $\Delta P_{A,j}^{+}$, $\Delta P_{B,j}^{+}$ and $\Delta P_{C,j}^{+}$ indicate the dynamic pressure drop for region 1, region 2 and region 3 in the static condition, respectively. All these terms, including $\Delta P_{ch0,j}^{+}$, $C_{A,j}$, $C_{B,j}$ and $C_{C,j}$, are defined in Lee et al. [14]. Moreover, $\Delta P_{ses,j}^{+}$ is seismic-induced pressure drop under vertical vibration condition, as given below:

(8)$\Delta P_{ses,j}^{+}=\left[C_{A,j}{\lambda }_{A,j}^{+}+C_{B,j}({\lambda }_{B,j}^{+}-{\lambda }_{A,j}^{+})+C_{C,j}(1-{\lambda }_{B,j}^{+})\right]a_{ses}^{+}/Fr$

Through Equations (5) and (6), it can be deduced that the dynamic equations of inlet flow velocity, $u_{i,j}^{+}$, among parallel supercritical heated channels under vertical seismic acceleration are:

(9)$\frac{d{u}_{i,j}^{+}}{d{t}^{+}}=Coe{1}_j\frac{d{u}_{i,1}^{+}}{d{t}^{+}}+Coe{2}_j,j=2,3,\cdots ,M$

(10)$Coe{1}_j=Coe{f}_1/Coe{f}_j$

(11)$Coe{2}_j=(\Delta P_{ch0,1}^{+}+\Delta P_{ses,1}^{+}-\Delta P_{ch0,j}^{+}-\Delta P_{ses,j}^{+})/Coe{f}_j$

In addition, for this parallel channel system, the total flow rate is the summation of every channel flow rate. The variation in total flow rate ($W_{tot}^{+}$) subject to seismic motion can be considered by the time-derivative term:

(12)$\frac{d{W}_{tot}^{+}}{d{t}^{+}}={\displaystyle \sum _j{\rho }_{i,j}^{+}\frac{d{u}_{i,j}^{+}}{d{t}^{+}}}$

Substituting Equation (9) into Equation (12) can lead to the dynamic equation of inlet flow velocity,$u_{i,1}^{+}$, for the first channel under vertical seismic acceleration:

(13)$\frac{d{u}_{i,1}^{+}}{d{t}^{+}}=\left(\frac{d{W}_{tot}^{+}}{d{t}^{+}}-{\displaystyle \sum _j{\rho }_{i,j}^{+}Coe{2}_j}\right)/{\displaystyle \sum _j{\rho }_{i,j}^{+}Coe{1}_j}$

Through Equation (13) together with Equation (9) and adopting all the other dimensionless dynamic equations given in Lee et al. [14], the present non-linear dynamic model can be employed to explore the influences of vertical seismic accelerations on the system with parallel supercritical water uniformly-heated channels. If the analytical system has a constant total flow rate during seismic motion, the term of $d{W}_{tot}^{+}/d{t}^{+}$ is set to zero.

On the other hand, if the system only comprises single supercritical heated channel, the boundary condition can be simplified as: the dynamic pressure drop of the supercritical uniformly-heated channel under vertical seismic accelerations ($\Delta P_{ch}^{+}$) is equal to the summation of pressure drop for the static system and that additionally contributed by the external vertical acceleration ($\Delta P_{ses}^{+}$). Thus,

(14)$\Delta P_{ch}^{+}=\Delta P_A^{+}+\Delta P_B^{+}+\Delta P_C^{+}+\Delta P_{ses}^{+}=Coef\frac{d{u}_i^{+}}{d{t}^{+}}+\Delta P_{ch,0}^{+}+\Delta P_{ses}^{+}$

where $\Delta P_A^{+}$, $\Delta P_B^{+}$ and $\Delta P_C^{+}$ represent the dynamic pressure drop for region 1, region 2 and region 3, respectively, inside single supercritical heated channel in the static condition. The corresponding formulas in a parallel channel system [14] can be applied to a single channel case.

Through Equation (14), the dynamic equation of inlet flow velocity, $u_i^{+}$, for the single supercritical heated channel under vertical seismic acceleration can be acquired:

(15)$\frac{d{u}_i^{+}}{d{t}^{+}}=(\Delta P_{ch}^{+}-\Delta P_{ch,0}^{+}-\Delta P_{ses}^{+})/Coef$

The procedure of the solution method is described as following. In the static condition, by eliminating the time-derivative terms in all the dynamic equations, it can result in a set of nonlinear algebraic equations. The subroutine SNSQE of Kahaner et al. [23] that takes the Powell hybrid scheme is applied to solve the steady states of parallel channel system or single channel system. By imposing the vertical seismic acceleration on the initial steady state of the system, the corresponding transient responses of parallel heated channels or single heated channel are solved by the set of non-linear, ordinary differential equations. The subroutine SDRIV2 of Kahaner et al. [23] employing the gear multi-value method is applied to obtain the numerical solutions. The dimensionless time step is set to $\Delta t^{+}={10}^{-4}$ through numerical calculations. The requested relative accuracy of 10⁻¹⁰ is set in all solution components.

3. Validation of the Model

Our previous study [14] has validated the present nonlinear dynamic model with supercritical water in the static condition ($a_{ses}^{+}=0$), i.e., without seismic motion, against the experimental data of a supercritical water two uniformly-heated channel system [15]. The results indicated that the present non-linear model can predict the stability thresholds well against the experimental data, within a relative error of 10% [14].

To evaluate the applicability of the external force method in estimating the vibration impact that was not well validated in the past studies [17,18,19,20,21], the authors established a vertical vibration experimental system in Figure 2a to investigate the thermal-hydraulic characteristics of boiling two-phase flow under vertical vibrations [24]. The test section comprises a single-phase stand pipe (D_SP = 0.012 m, L_SP = 0.17 m), a heated section (D_H = 0.012 m, L_H =1.27 m) and an adiabatic riser (D_R = 0.012 m, L_R = 0.38 m). The additional pressure drop caused by the vertical vibration ($\Delta P_{VIB}$) can be estimated by the measured pressure drop ($\Delta P(Exp.)$) under vertical vibration subtracting the average pressure drop ($\Delta P_{ave}$). Thus, the non-dimensional form is:

(16)$\Delta P_{VIB}^{+}(Exp.)=[\Delta P(Exp.)-\Delta P_{ave}]/{\rho }_f{u}_s^2$

In terms of external force method, the calculated value of $\Delta P_{VIB}^{+}$ can be determined by Equation (17):

(17)$\Delta P_{VIB}^{+}(Predict.)=[({\rho }_i^{+}+{\rho }_e^{+})/2+{\rho }_e^{+}{L}_R^{+}+{\rho }_i^{+}{L}_{SP}^{+}]\times a_{\max }\sin (2\pi f^{+}{t}^{+})/Fr$

For a single heated channel vertical vibration system at atmosphere pressure, Figure 2b compares the peak magnitude of the calculated $\Delta P_{VIB}^{+}(Predict.)$ against the peak magnitude of the experimental $\Delta P_{VIB}^{+}(Exp.)$ among nine test cases with heat fluxes of q” = 8444, 11,029 and 13,959 W/m² and vibration frequencies of f = 0.5, 1.0 and 1.5 Hz, respectively.

The results exhibit that the vibration-induced pressure drop increases as the increase in the vibration frequency, i.e., vibration acceleration. The theoretical predictions using Equation (17) can predict the experimental values based on Equation (16) within a relative error of about 20% among all the cases. It can demonstrate that external force method is appropriate to explore the impact of external vertical vibration on the rigid system. It should note that the discrepancy between the experimental data and the predicted value is smaller in the relatively small vibrations with f = 0.5 and 1.0 Hz than that in the relatively large vibration with f = 1.5 Hz. This may suggest that the external force method is more applicable for the relatively small vibration acceleration, non-resonance and insignificant resonance conditions under the vibration impacts, as in most of the cases in the present study, in which the probability of the deformation occurring in the structure parts is very limited during such vertical vibrations. In addition, at present, the experiment data of the supercritical fluid system under external vibrations is deficient in the literature. The external force method may need a further validation against the vibration experimental data of supercritical fluid system if it is available in the future. On the other hand, the actual seismic impact on a power plant is very complicated, which could be coupled with the fluctuations of coolant flow and system structures. If either two of their vibration frequencies are approaching, it may trigger the resonance oscillations to deteriorate the system operation and cause the deformation of structure parts. The assumption of a rigid body would be an oversimplification under such conditions. If we consider the possible deformation of structure components facing the extreme large intensity of seismic wave or very strong resonance oscillations, the other methodologies, such as the moving grid method, need to be developed in evaluating the seismic impact. However, under such extreme conditions, the power plants with the auto scram function will shut down the system immediately to prevent it from the progressive event. Therefore, much more attention should be paid to the seismic wave below the scram limit of the power plant, which the external force method could be suitably applied to.

4. Results and Discussion

4.1. Stability Map of the Static System

In SCWR, the reactor core involves parallel supercritical heated channels. Table 1 lists the reference condition for three representative supercritical heating channels, which are extracted from the design data for U.S. reference SCWR [25]. In the present model, the fluid properties of supercritical water are obtained from the database of Nation Institute of Standards and Technology (NIST) [26], which have been widely used in the literature [5,6,7,8,9,10,11,12,13,14,15,16]. Based on the parameter values given in Table 1, this non-linear model can explore the stability maps of the static system by setting $a_{ses}^{+}=0$ through time domain evolutions. With a perturbation in the inlet velocity of the highest heating channel (Ch.1) at a constant pseudo-subcooling number of N_P-SUB = 1.55, Figure 3a illustrates the examples of transient responses corresponding to a stable state, a boundary state and an unstable state of a three supercritical water-heated channel system with a radial power distribution of 1.2:1.0:0.8. Under the constraint by constant total flow rate, all the results in Figure 3a show that the highest amplitude oscillation in the highest heating channel (Ch.1) is out-of-phase with those of the other less-heated channels (Ch.2 and Ch.3). The stability boundary could be judged by the transient responses among three asymmetric heating channels. For a stable state of N_P-SUB = 1.55 and N_TPC = 5.645, the perturbation would decay to its original steady state, while for an unstable state of N_P-SUB = 1.55 and N_TPC = 5.704, the time domain oscillations among channels are growing and eventually trigger the reversed flow. In addition, the results in Figure 3a display the fact that transient responses among channels finally evolve to stable limit cycle oscillations at a boundary state of N_P-SUB = 1.55 and N_TPC = 5.669.

Through non-linear analyses as the examples illustrated in Figure 3a, the results in Figure 3b reveal the comparison of stability boundary between parallel three heated channels and single heated channel with supercritical water under the system pressure of 25 MPa. These stability maps are presented based on the pseudo-subcooling number (N_P-SUB) and the average true trans-pseudo-critical number (N_TPC), which were originally introduced by Ambrosini and Sharabi [8].

(18)$N_{TPC}=\frac{Q_{tot}}{W_{tot}}\frac{{\beta }_{PC}}{C_{P,PC}},N_{P-SUB}=\frac{{\beta }_{PC}}{C_{P,PC}}\left(h_{PC}-h_i\right)$

The three-region methodology for the supercritical water-heated system at 25 MPa can be applied when the inlet flow temperature is between 0 °C and 349 °C while two-region model, excluding region 1, is employed instead when inlet flow temperature is beyond 349 °C (N_P-SUB < 0.904). Moreover, the case with the inlet fluid temperature below 0 °C (N_P-SUB > 3.59) implies a restricted area that prohibits the supercritical water flow. The results in Figure 3b suggest that the system of parallel three channels with a heat flux ratio of 1.2:1.0:0.8 is more stable than the single channel system in this study. Our previous study [27] demonstrated that the stable effect due to system mass and the unstable effect by channel-to-channel interaction could prevail in a multi-channel system. In this case, it seems that the stable effect would dominate over the unstable effect. However, when a more asymmetric heat flux distribution of 1.4:1.0:0.6 is set, the parallel three channel system would become more unstable than the single channel one. This resulted from the more dramatic unstable channel-to-channel interaction being triggered by the more difference among channels. In addition, the normal operating state of N_P-SUB = 1.55 and N_TPC = 3.232, as denoted as S1 in Figure 3b, is substantially distant from the stability boundaries of all systems. This infers that it is a very stable operating point.

4.2. Non-Linear Seismic-Induced Oscillations

4.2.1. Seismic-Induced Oscillation of Single Channel System

An actual vertical seismic acceleration wave in Figure 4a is considered in the present study, which is extracted from the data recorded by the Linkou Station during the Chi-Chi earthquake in Taiwan that took place on 21 September 1999 [28]. The data were collected in a time step of 0.005 s during a time interval of 88 s upon this earthquake, which corresponds to a dimensionless time interval of nearly t⁺ = 35 and a dimensionless time step of almost Δt⁺ = 0.002. Supposing this actual vertical seismic acceleration wave is directly imposed on the S1 state, i.e., normal operating state as denoted in Figure 3b, of the single channel system, the inlet flow velocity oscillation in Figure 4b shows that the seismic wave dominates the channel oscillation. The flow oscillates chaotically following the vibration of seismic acceleration with the same frequency; however, the oscillation amplitudes remain bound in a certain range, inferring that seismic-induced oscillation of this supercritical uniformly-heated channel at S1 state is still small. The seismic-induced oscillation could result from the compound results of the inherent stability characteristics interacting with the effect of vertical seismic acceleration. The S1 state is very stable, as revealed in Figure 3b, such that the seismic acceleration just triggers a small oscillation.

On the other hand, at the same inlet fluid temperature with N_P-SUB = 1.55, the operating power is increased to another S2 stable state with N_TPC = 5.645, which is rather near the stability boundary of a single channel system, as revealed in Figure 3b. Assuming the same vertical seismic acceleration in Figure 4a is imposed on the S2 stable state, the results in Figure 5a illustrate that the seismic-induced oscillation in the S2 state is much larger than that in S1 state. This demonstrates that the strength of seismic-induced oscillation is related to the inherent stability characteristics of the initial state. For the S1 stable state, Figure 5b explores the effect of inlet fluid temperature on the seismic-induced oscillation of the single heated channel system. Under the same heating power (N_TPC = 3.232), the inlet fluid temperature is step decreased from N_P-SUB = 1.2, N_P-SUB = 1.55 (S1) to N_P-SUB = 2.02. The results in Figure 5b indicate that the reduction in inlet fluid temperature tends to stabilize the seismic-induced oscillation. It will increase the inlet fluid density to extend the heavy fluid region (region 1) inside the channel, thereby increasing the pressure drop contributed by the heavy fluid area; in addition, it will also drive the system to have a heavier inertial mass. Our previous study [14] showed that increasing the pressure drop of the heavy fluid region in the channel would have a stable effect on the supercritical forced convection heating channel system. Therefore, Figure 5b shows that due to the effects of both heavy fluid region and inertial mass, lowering the inlet fluid temperature will have a stable effect on the transient oscillation of the system caused by the seismic acceleration wave.

In a supercritical water-heated channel system, such as the SCWR, the geometrical parameters may affect the transient oscillations of the system induced by the seismic acceleration wave. Under the same vertical seismic acceleration wave shown in Figure 4a and based on the S1 stable state, Figure 6a investigates the effect of inlet flow resistance (k_i) on the seismic-induced oscillation of a single supercritical heated channel. The results illustrate that the increase of inlet loss coefficient will drive the system to show a smaller non-linear oscillation. Increasing the flow resistance at channel inlet will increase the pressure drop of heavy fluid region. Moreover, the increase of inlet flow resistance will also attenuate the seismic-induced oscillation of the system. Thus, the increase of the inlet loss coefficient would generate a stable effect on the seismic-induced oscillation.

Figure 6b further evaluates the influence of the outlet loss coefficient (k_e) on the seismic-induced oscillation of a single supercritical heated channel. The analytical results show that a larger outlet loss coefficient will result in a more prominent seismic-induced oscillation. The increase of outlet loss coefficient will increase the pressure drop of light fluid region inside the channel. Our previous study [14] indicated that increasing the pressure drop in the light fluid area could destabilize the system. On the other hand, the increase of outlet flow resistance would also produce an attenuated influence on the seismic-induced oscillation. It seems that the unstable effect is dominant such that increasing the outlet loss coefficient (k_e) tends to destabilize the seismic-induced oscillation of this single heated channel system.

Figure 7 explores the influence of the heated length (L_H) and channel diameter (D_H), respectively, on the transient oscillations of a single supercritical heated channel caused by vertical seismic acceleration wave in Figure 4a. At a fixed uniform heat flux, a longer heating channel implies a higher input power, which will extend the length of the light fluid region. Moreover, lengthening the channel will enlarge the seismic-induced pressure drop estimated by Equation (8). Thus, the analytical results in Figure 7a show that a longer channel tends to destabilize the seismic-induced oscillation of single supercritical heated channel. On the other hand, at a fixed total flow rate, increasing the channel diameter will reduce the flow velocity; thereby a smaller channel pressure drop would stabilize the system. Therefore, Figure 7b illustrates that a larger channel diameter could generate a stable effect on the seismic-induced oscillation of a single supercritical heated channel.

4.2.2. Seismic-Induced Oscillation of Parallel Three Channel System

The seismic impact on the parallel channel system is related to the inherent natures of the initial operating state, particularly the resonance oscillation could be triggered while the seismic frequency approaches the natural frequency of the working flow. Figure 8 illustrates the variation in the dimensionless natural frequency (f⁺) along with the stability boundary of three supercritical water heated channels with a heat flux ratio of 1.2:1.0:0.8 at a supercritical pressure of 25 MPa. This figure shows that the distribution of the dimensionless natural frequency on the stability boundary is ranged between 0.93 and 4.16. The parameter of N_P-SUB indicates the effect of inlet fluid temperature. A larger value of N_P-SUB infers a lower inlet fluid temperature, and vice versa. In the region with N_P-SUB ≥ 1.67, where the heavy fluid region would contribute a substantial role, the dimensionless natural frequency on the stability boundary is mildly increased as N_P-SUB is decreased. In addition, at a fixed value of N_TPC = 6.15, the dimensionless natural frequencies at the stable states with N_P-SUB = 2.882, 2.811 and 2.727 are 0.98, 1.1 and 1.14, respectively, as labeled in Figure 8. It also suggests the similar slight increasing trend in this region if the inlet fluid temperature is increased. However, when N_P-SUB is less than 1.67, the dimensionless natural frequency on the stability boundary would drastically increase as N_P-SUB is further decreased due to the expansion of the light fluid region dominating the channels.

On the other hand, the parameter of N_TPC may include the effects of heat flux, heated length and channel diameter, which is generally proportional to heat flux (q”) and heated length (L_H) but inversely proportional to channel diameter (D_H). The increase of N_TPC can be reached by the increase of heating power (Q) or the reduction of total flow rate (W_tot), and vice versa. At a fixed value of N_P-SUB = 2.882, the dimensionless natural frequencies at the stable states with N_TPC = 6.15, 6.234 and 6.298 are 0.98, 1.05 and 1.09, respectively, as labeled in Figure 8. This implies that the dimensionless natural frequency would increase as the increase of N_TPC due to the expansion of light fluid region.

For the supercritical water three-heated channels with an asymmetric radial power ratio of 1.2:1.0:0.8 considered in this study, the seismic-induced effect on the system transients is investigated in S1 stable state, i.e., normal operating state as denoted in Figure 3b. The inherent natural frequency of S1 state is approximately f⁺ = 1.1. Figure 9a assumes a set of external vertical sinusoidal acceleration waves, including two vertical acceleration waves with the same peak acceleration (a_peak = 0.1 g) but with different oscillation frequencies of f⁺ = 0.5 (non-resonance condition) and 1.1 (resonance condition), respectively. This study supposes that the total flow rate of parallel supercritical heated channels keeps constant under the impact of the external vertical acceleration wave, as shown in Figure 9b. As the vertical sinusoidal acceleration wave in Figure 9a is imposed on the system, Figure 9c shows the vibration-induced amplitude oscillations of the inlet flow velocity ($\delta u_i^{+}$) among these three supercritical heated channels. The channel oscillations triggered by the vertical acceleration wave can be regarded as the result of the inherent stability characteristics of the initial state combined with the influence of external vertical acceleration wave. Since the normal operating point (S1) is very stable, the channel-to-channel interaction among parallel heated channels is fairly weak. In the case of non-resonance condition (f⁺ = 0.5), Figure 9c indicates that the oscillation of each channel caused by the external vertical acceleration wave appeared to be very small due to the inherent stability nature of S1 state. In addition, under the external vertical acceleration, out-of-phase mode of inlet velocity oscillations is presented among channels to maintain a fixed total flow arte. In the case of resonance condition (f⁺ = 1.1), the vibration-induced oscillation of each channel will become larger due to the resonance effect. However, as the normal operating point (S1) is very stable, the resonance effect seems to be not significant, which does not trigger the uncontrollable oscillations under the impact of the external vertical acceleration wave in Figure 9a.

For a parallel supercritical heated channel system, such as SCWR, the inlet loss coefficient in each channel is set, as listed in Table 1, such that the similar power-to-flow ratio among channels is reached in the normal operating point (S1) of the forced circulation system. In general, the most heated channel, i.e., channel 1, is with the smallest inlet loss coefficient and thus the highest mass flow rate, while the least heated channel, i.e., channel 3, has the largest inlet loss coefficient that makes it have the lowest mass flow rate. Based on the normal operating point (S1) in Figure 3b as a reference case and under the same external vertical acceleration wave in Figure 9a, Figure 10a investigates the effect of inlet loss coefficient at the most heated channel (channel 1) on the vibration-induced oscillations among parallel supercritical heated channels. Assuming the group of inlet loss coefficient in Table 1 is step set to channel 1 from k_i,1 = 22.7 (S1), 50.0, 93.0 to 241.1, the results in Figure 10a reveal that with the increase of inlet flow resistance set at the highest heating channel, the vibration-induced oscillation of the most heated channel firstly become smaller while the inlet loss coefficient is set to k_i,1 = 50.0, due to both the stable effect and the attenuated effect of inlet flow resistance. However, if the inlet flow resistance of the hottest channel is further increased to k_i,1 = 93.0 and 241.1, the most heated channel presents a larger vibration-induced oscillation, which infers an unstable influence on the system dynamics. This is caused by the more uneven power-to-flow distribution among channels resulted by the increase in the inlet loss coefficient of channel 1, particularly for the case with k_i,1 = 241.1. The more difference among channels will induce a stronger channel-to-channel interaction to destabilize the vibration-induced oscillation when its effect is greater than the stable and attenuated effects of inlet flow resistance.

Under the same external vertical acceleration wave in Figure 9a, Figure 10b further explores the influence of outlet flow resistance (k_e) on the transient oscillations of the parallel three supercritical heated channel system caused by external vertical sinusoidal acceleration wave. Based on the case of normal operating state (S1), the results illustrate that as the outlet loss coefficient is increased from k_e = 0.5, 1.0 (S1), to 2.0, the most heated channel will show a smaller transient oscillation, which implies that vibration-induced oscillation is attenuated by the outlet flow resistance. However, if the outlet loss coefficient is further increased to an extraordinary condition with k_e = 20.0, the most heated channel will exhibit a dramatically larger vibration-induced oscillation in the case of resonance condition (f⁺ = 1.1), while a very small one is presented in the non-resonance condition (f⁺ = 0.5). Our previous study [14] revealed that increasing the outlet loss coefficient would destabilize the system of parallel supercritical heated channels in the static system, which strengthen the instability degree of the system. Thus, the unstable effect of outlet flow resistance combined with the resonance effect will cause the system to exhibit much larger oscillation than in the case of non-resonance. Moreover, the resonance oscillation with k_e = 20.0 is much stronger than that with k_e = 1.0 in Figure 9c. This suggests that the strength of vibration-induced resonance oscillation is closely related to the inherent stability characteristics of initial state.

The feature of a multi-channel system would exhibit an asymmetric radial heat flux ratio distributed among channels, which may affect the seismic-induced oscillation. The power range relative to the average channel for three channel groups in SCWR is above 1.2 for group 1, 0.9 to 1.19 for group 2 and below 0.9 for group 3, as displayed in Table 1 [25]. Therefore, supposing the same external vertical acceleration wave in Figure 9a is imposed on the normal operation point S1 of the parallel three supercritical heated channel system, as denoted in Figure 3b, the heat flux ratios of 1.1:1.0:0.9, 1.2:1.0:0.8 and 1.3:1.0:0.7 are selected to explore their effects on the vibration-induced oscillations, as indicated in Figure 11. The results suggest that the more difference in the heat flux among channels would drive a stronger channel-to-channel interaction, and thus cause a much larger vibration-induced oscillation as in the case of 1.3:1.0:0.7 shown in Figure 11.

Based on the normal operation point S1 of the parallel three supercritical-heated channel system, as marked in Figure 3b, and assuming vertical seismic acceleration wave in Figure 4a imposed on the system, Figure 12 reveals the seismic-induced inlet velocity oscillations of parallel three supercritical heated channels with a heat flux ratio of 1.2:1.0:0.8 while a constant total flow rate is maintained at S1 state. It is found that since the normal operating point S1 is very stable, the interactions among parallel three channels are rather weak, thereby a vertical seismic acceleration wave would dominate the non-linear behaviors of the system oscillations. The inlet flow velocities among channels oscillate chaotically following the motion of vertical seismic acceleration wave with the same oscillation frequency. However, their oscillation amplitudes are still limited to a small range, indicating that the normal operating point S1 is still safe subject to this seismic wave. In addition, since the total flow rate is supposed to keep constant during the seismic motion, the inlet velocity oscillations among parallel three channels are out-of-phase, as illustrated in Figure 12.

With the same pseudo subcooling number (N_P-SUB = 1.55), the seismic-induced effect is explored in another stable state of S3 with a high operating power of N_TPC = 5.505, which is more approximate to the stability boundary of a parallel three supercritical heated channel system with a heat flux ratio of 1.2:1.0:0.8, as denoted in Figure 3b. Supposing the same seismic acceleration in Figure 4a is imposed on the S3 state, Figure 13 exhibits the transient oscillations of parallel three supercritical heated channels caused by this vertical seismic acceleration wave. The results show that the parallel three channels will present larger oscillations at S3 stable state as compared to those at S1 state in Figure 12. This has resulted from the inherent stability characteristics of the initial state. Because S3 state rather than S1 state is closer to the stability boundary of the system with a heat flux ratio of 1.2:1.0:0.8, the instability degree and interaction among channels in the S3 state are stronger than those in the S1 state. Thus, the seismic-induced oscillation is more significant in the S3 state versus those in the S1 state.

5. Conclusions

Incorporating the external force method with the three-region methodology in analyzing the supercritical heated channel, this study develops a non-linear dynamic model to investigate the transient behaviors of single channel or parallel channels with supercritical water under the impacts of vertical sinusoidal and seismic accelerations. The present model is validated reasonably based on the experimental data of a single heated channel vertical vibration system. Selecting an actual vertical seismic acceleration wave and the normal operating state as a reference case, the conclusions in the present study can be summarized as:

• The strength of seismic-induced oscillation is relevant to the inherent stability natures of the initial state.

• The effects of parameters on the seismic-induced oscillation conducted in single heated channel indicate that unstable influences could be caused by a longer heated length, uprating operation power and the increase of outlet loss coefficient, while stable effects could be illustrated by the increase of inlet loss coefficient, a larger channel diameter and a lower inlet fluid temperature.

• For a three supercritical heated channel system with a heat flux ratio of 1.2:1.0:0.8, in the region with N_P-SUB greater than 1.67, the dimensionless natural frequency on the stability boundary is mildly increased as N_P-SUB is decreased, where the heavy fluid region contributes a substantial role in the channels. However, in the area with N_P-SUB less than 1.67, it would significantly increase as N_P-SUB is further decreased, which is due to the effect of light fluid region dominating over the channels.

• The seismic-induced effects on parallel heated channels are related to the inherent stability characteristics of initial state and the interactions among channels. The effects of inlet loss coefficient and outlet loss coefficient on the vibration-induced oscillations among channels may depend on the operating conditions, which are affected by the power-to-flow difference among channels and stability nature of the system. The more uneven heat flux distribution among channels would cause a larger vibration-induced oscillation.

• The strength of vibration-induced resonance oscillation among parallel channels is closely related to the inherent stability characteristics of initial state. It may trigger much larger resonance oscillations than in the case of non-resonance.

Author Contributions

Conceptualization, J.D.L.; methodology, J.D.L.; software, J.D.L.; validation, J.D.L. and S.W.C.; investigation, J.D.L.; data curation, S.W.C.; writing—original draft preparation, J.D.L.; writing—review and editing, S.W.C.; visualization, J.D.L.; supervision, J.D.L. and S.W.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Ministry of Science and Technology of Taiwan, grant number MOST 107-2221-E-007-096-MY2 and MOST 108-2221-E-007-046-MY3.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

a	acceleration (ms−2)	NP-PCH,j	pseudo-phase change number for j-th channel, $=\frac{Q_j}{{\rho }_{PC}{A}_H{u}_s}\frac{{\beta }_{PC}}{C_{P,PC}}$
apeak	peak acceleration (ms−2)
amax	non-dimensional peak magnitude
$a_{ses}^{+}$	non-dimensional seismic acceleration, = ases/g	NTPC,j	true trans-pseudo-critical number for j-th channel, $=\frac{Q_j}{{\rho }_{i,j}{A}_H{u}_s}\frac{{\beta }_{PC}}{C_{P,PC}}$
$a_V^{+}$	non-dimensional vertical acceleration, = aV/g
A	cross sectional area (m2)	NTPC	Average true trans-pseudo-critical number, $=\frac{Q_{tot}}{W_{tot}}\frac{{\beta }_{PC}}{C_{P,PC}}$
AH	cross sectional area of heated channel (m2)
Cp	constant pressure specific heat (Jkg−1K−1)
D	diameter (m)	NP-SUB,j	pseudo-subcooling number for j-th channel, $=\frac{{\beta }_{PC}}{C_{P,PC}}\left(h_{PC}-h_{i,j}\right)$
DH	diameter of heated channel (m)
f	friction factor or frequency (Hz)
$f^{+}$	non-dimensional frequency, $=f{L}_H/u_s$	P	pressure (Pa)
		Q	heating power (W)
Fr	Froude number, = $u_s^2$/gLH	q”	heat flux (Wm−2)
		T	temperature (K)
g	gravity acceleration (ms−2)	TPC	pseudo-critical temperature (K)
		t	time (s)
h	enthalpy (Jkg−1)	tref	time scale, = LH/us
hPC	fluid enthalpy at pseudo-critical point (Jkg−1)	t +	non-dimensional time, = t/tref
hi	inlet fluid enthalpy (Jkg−1)	u	velocity (ms−1)
h +	non-dimensional fluid enthalpy, $=\left(h-h_{PC}\right){\beta }_{PC}/C_{P,PC}$	ui0	steady state inlet velocity (ms−1)
		us	average steady state inlet velocity (ms−1), = $W_{tot}/{\rho }_i{A}_{tot}$
		u+	non-dimensional velocity, $=u/u_s$
k	loss coefficient	W	mass flow rate (kgs−1)
L	length (m)	W +	non-dimensional mass flow rate, $=W/{\rho }_{PC}{A}_H{u}_s$
LH	channel length (m)	z	axial coordinate (m)
L+	non-dimensional length, = L/LH	z+	non-dimensional axial coordinate, = z/LH
M	mass (kg)
M+	non-dimensional mass, $=M/{\rho }_{PC}{A}_H{L}_H$
Greek Symbols		Subscripts
β	isobaric thermal expansion coefficient (K−1)	ave	average
		A	region 1 or boundary point of region 1
		B	region 2 or boundary point of region 2
ΔP	pressure drop (Pa)	ch	channel
		C	region 3
ΔP+	non-dimensional pressure drop, $=\Delta P/{\rho }_{PC}{u}_s^2$	e	exit of heated channel
		f	saturated liquid
		H	heated section
δx	(x-x0) for variable x, x0 represents the steady-state value	i	inlet of heated channel
		j	j-th channel
		P	constant pressure
ρ	density (kgm−3)	PC	pseudo-critical
ρ+	non-dimensional density, $=\rho /{\rho }_{PC}$	R	riser
		ses	seismic
		SP	stand pipe
${\rho }_{PC}$	fluid density at pseudo-critical point (kgm−3)	tot	total
		VIB	vibration
Λ	friction number, $=f{L}_H/2D_H$	0	steady state
λ	the boundary of region (m)
λ+	non-dimensional boundary of region, = λ/LH

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Modelling structure of parallel supercritical water uniformly-heated channels analyzed by three-region methodology [7] under vertical seismic acceleration. (“A” and “B” indicate the separation points between regions).

Figure 2. The validation of external force method: (a) single heated channel vertical vibration experimental system [24]; (b) comparisons of the vibration-induced pressure drop using Equation (17) against the experimental values through Equation (16).

Figure 3. Investigations of supercritical water heated systems through non-linear analysis: (a) examples of stable state, boundary state and unstable state; (b) stability map of the static system with parallel three heated channels or single heated channel.

Figure 4. The influence of seismic acceleration on the single channel system at S1 state: (a) an actual vertical seismic acceleration recorded during the Chi-Chi earthquake in Taiwan [28]; (b) seismic-induced inlet velocity oscillation.

Figure 5. The influence of parameter on the seismic-induced oscillation: (a) the effect of operating state with uprating power; (b) the effect of inlet flow temperature.

Figure 6. The influence of parameter on the seismic-induced oscillation: (a) the effect of inlet loss coefficient; (b) the effect of outlet loss coefficient.

Figure 7. The influence of parameter on the seismic-induced oscillation: (a) the effect of heated length; (b) the effect of channel diameter.

Figure 8. The variation in the dimensionless natural frequency along with the stability boundary of three supercritical heated channels with a heat flux ratio of 1.2:1.0:0.8 at 25 MPa.

Figure 9. The impact of external vertical acceleration on the parallel three channel system at S1 state: (a) external vertical sinusoidal acceleration with apeak= 0.1 g; (b) constant total flow rate under vertical vibrations; (c) vibration-induced inlet velocity oscillations among channels.

Figure 10. The parametric effects on vibration-induced oscillation of the most heated channel: (a) the effect of inlet loss coefficient at channel 1; (b) the effect of outlet loss coefficient.

Figure 11. The effects of heat flux distribution among channels on vibration-induced oscillation of the most heated channel.

Figure 12. Seismic-induced inlet velocity oscillations of parallel three channel system having a constant total flow rate at S1 state subject to a real vertical seismic acceleration.

Figure 13. Seismic-induced inlet velocity oscillations of parallel three channel system having a constant total flow rate at S3 state subject to a real vertical seismic acceleration.

Table 1

The reference case of three supercritical heated channels [25].

Parameter	Channel 1	Channel 2	Channel 3
Heat flux ratio	above 1.2	0.9 to 1.19	below 0.9
Length (LH)	4.2672 m	4.2672 m	4.2672 m
Diameter (DH)	3.4 mm	3.4 mm	3.4 mm
Outlet loss coefficient (ke)	1.0	1.0	1.0
Inlet loss coefficient (ki)	22.7	93.0	241.1