A Comparative Assessment on Different Aspects of the Non-Linear Instability Dynamics of Supercritical Fluid in Parallel Channel Systems

Abstract

The thermal-hydraulic behavior of supercritical water reactors with a parallel channel configuration was examined through a non-linear instability analysis. This analysis was performed under wide-ranging conditions and aspects, including different working supercritical fluids, varied heat-flux and flow-rate conditions, and channel inclinations. The supercritical fluid (SCFs) dynamics were captured using the density, enthalpy, and velocity analytical approximation functions. The major findings show that both SCFs (water and carbon dioxide) experienced density wave oscillations at a low pseudo-subcooling number. Static instability characteristics were observed for supercritical water, at a relatively high subcooling number. Further, it was found that at different heat flux, the hotter channel makes the overall system more unstable, whereas, at equal heat flux, parallel channels perform similar to a single-channel system. However, the effect of the inclination angle was found to be negligible owing to supercritical pressure conditions. Moreover, stable and unstable limit cycles along with out-of-phase oscillation characteristics were observed in dynamic stability regions. The present model was also compared with experimental and numerical data. Moreover, co-dimension and numerical simulations were performed to confirm the observed non-linear characteristics. This study helps to enhance the heat transfer characteristics during safe operation of heated channel systems, such as nuclear reactors and solar thermal systems.

1. Introduction

Over the last few decades, supercritical fluids (SCFs) have been considered as a valuable option for working fluids in several energy systems due to their high thermal efficiencies. It has been observed that close to the pseudo-critical temperature, SCFs’ properties experience significant changes [1]. Numerous analyses have been performed to study the impact of these drastic changes on the overall performance of the supercritical Brayton cycle [1,2,3]. In comparison to other two-phase flow energy systems, boiling and phase-change do not take the place. Since there is no clear understanding of the phase change effect on dynamical system behavior, it becomes quintessential to doubt the existence of flow instabilities in supercritical systems [4,5]. Flow instabilities are categorized into two parts: density wave, pressure drop, and thermal oscillations, which are characterized mainly by dynamic instability; and Ledinegg, which is commonly reported in static instability [6,7,8]. Several authors have adopted the two-phase stability analysis methodology to develop single and parallel channel models to study the stability phenomena in supercritical systems; however, these reported investigations have only analyzed linear stability [4,9,10]. Several authors have proposed a new dimensionless parameter and have compared it with two-phase flow parameters [4,10,11].

In several thermal energy systems, parallel channel configurations have been widely used to fulfill the system objective, where the channel-to-channel interaction, heat flux distribution, and flow rate have been found to be critical parameters. These parameters significantly influence the thermal-hydraulic behavior of the system. Numerous experimental [12,13,14] and numerical [4,15,16,17,18,19] research has been carried out on the flow instability in parallel channels under symmetric heating conditions. However, an absolute symmetric heating condition does not exist in actual operation since most industrial equipment always operates under dissymmetric heating conditions. Dissymmetric heating induces an uneven flow distribution among the channels and further affects the flow instability of the parallel channel system. A previous study indicated a considerable difference between the flow instability under symmetric and dissymmetric heating conditions [20,21,22]. However, a scarce amount of research has studied the flow instability in parallel channel systems under the dissymmetric heating condition.

Parallel channel instability is observed in most heat exchangers and nuclear reactor cores, in which the instability is significantly influenced by the channel-to-channel interactions. However, these studies concluded that in-phase and out-of-phase oscillations are expected for an identical two-channel system [4,5,23,24,25,26]. The different heat inputs between the channels and more channels in a system may cause complex channel-to-channel interaction modes, leading to system instability.

Ting Xiong experimentally [12] and numerically [15] studied the flow instability characteristics of a parallel channel system with SCF. In the experimental study, the combination of the heat flux and flow rate varied and was divided into different phases. It has been reported that the flow rate would experience more asymmetry at higher fluid temperatures, making it difficult to demonstrate flow instability experimentally.

Zhang, L. [13] also conducted experimental studies on a forced circulation loop under SCW with a parallel channel. Based on the oscillation characteristics, they reported that two types (I and II) of dynamic instabilities occurred. They observed in-phase and out-of-phase oscillations at low and high temperatures, respectively. Later, Xi [14,23] carried out out-of-phase oscillation using the CFX code for a three-dimensional (3-D) study, which, in comparison with experimental data, was able to better predict the flow instability than a one-dimensional model. However, improvement of the 3-D physical models for enhanced prediction of the period of oscillation has been recommended.

Su, Y. [18] adopted the time-domain method to theoretically study flow instability characteristics in the parallel channels for supercritical water. They compared the obtained characteristics with a two-phase flow system. They concluded that arc-shaped stability characteristics are significantly affected by the pressure drop, mass flow rate, and volume of supercritical water. Moreover, it has been reported that delayed properties and feedback effects should be considered to demonstrate instability phenomena.

Jin Der Lee [21,27] developed a time-domain model for heated channels based on a three-zone methodology. Earlier, only the polynomial approximation function for density and enthalpy had been included; however, later, they also added the impact of vertical seismic accelerations to investigate the non-linear transient characteristics of the supercritical water. Out-of-phase flow oscillations were observed in two channels, whereas three channels under different heat flux were used to demonstrate the in-phase oscillation. They only observed dynamic instability characteristics, whereas static and bifurcation analysis was missing in these studies.

Jialun Liu [19,20] also used the time-domain method to study flow instabilities due to dissymmetric heating conditions in two vertical upward channels. The impact of dissymmetric heating conditions was found to be dependent on the mass flow rate and inlet fluid temperatures. Based on this dependency, they evaluated the different positive and negative response characteristics of the system instability. The above-discussed literature review and concepts of bifurcation theory show that linear stability cannot demonstrate the entire system dynamics as it only predicts tiny or small perturbation impacts on the system dynamics. Thereby, this reveals the importance of non-linear studies in dynamical systems [6,7,8,24,25,26,27,28], such as two-phase flow, boiling water reactors, advanced heavy water reactors, heated channel systems, etc., for the normal and safe operation of such systems.

The present work exploits a previously developed methodology (nodalized reduced ordered model) for a single channel [29] to a parallel channel system. An ROM was developed to capture the non-linear characteristic of a parallel channel under different supercritical fluid (water and carbon dioxide) flow conditions. Apart from the in-house code, MATCONT software (bifurcation tool) was used [30], which works on a MATLAB platform. In order to demonstrate realistic engineering applications, the system was studied under different working conditions, such as an equal and un-equal heat flux distribution and varied flow rate. The major findings are the several bifurcation characteristics obtained along the stability boundary, depicting the system’s non-linear dynamical behavior. The model was validated using numerical and experimental literature data. The present study is instrumental in the design of the safe operation of nuclear reactors and solar thermal systems, in a wide parametric range.

2. Modelling Methodology

The objective of the present work is fulfilled using a system with two parallel channels. The physical schematic is shown in Figure 1, and the corresponding design and operating parameters are shown in

Table A1. The list of non-dimensional parameters used in the present study.

$\rho =\frac{{\rho }^{*}}{{\rho }_{pc}^{*}}$	$h=\frac{{\beta }_{pc}}{C_{p,pc}}\left(h^{*}-h_{pc}^{*}\right)$
$N_{tpc}=\frac{N_{tpc}^{\prime }}{{\rho }_{in}}$	$N_{spc}=\frac{{\beta }_{pc}}{C_{p,pc}}\left(h_{pc}^{*}-h_{in}^{*}\right)$
$w=\frac{w^{*}}{{w_{in}}^{*}}$	$N_{tpc}^{\prime }=\frac{q_o^{″}{}^{*}{\Pi }_h{}^{*}{L}^{*}}{{\rho }_{pc}^{*}{w}_o^{*}{A}^{*}}\frac{{\beta }_{pc}}{C_{p,pc}}$
$z=\frac{z^{*}}{L^{*}}$	$p=\frac{p^{*}}{{\rho }_{pc}^{*}{w_o^2}^{*}}$
$Fr=\frac{w_o^2{}^{*}}{g^{*}{L}^{*}}$	$t=\frac{t^{*}{w_o}^{*}}{L^{*}}$

. This physical system is replicated in the mathematical form using the 1-D Navier–Stokes equations [4,5,6,9,10]. These conservation equations are further simplified using non-dimensional parameters [11,29] under some appropriate assumptions as follows:

• Homogenous flow is considered in both channels.
• The inlet temperature remains constant.
• The system pressure remains constant to use the thermodynamics property of the fluids.

After using the above-mentioned assumptions and non-dimensional parameters, the non-dimensional 1-D Navier–Stokes equation is as follows:

Mass balance equation:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho w}{\partial z}=0$$

(1)

Momentum balance equation:

$$\frac{\partial \rho w}{\partial t}+\frac{\partial \rho w^2}{\partial z}+\frac{\partial p}{\partial z}=-\frac{\rho }{Fr}-\left[\Lambda +K_{in}{\delta }_d\left(z\right)+K_{exit}{\delta }_d\left(z-1\right)\right]\rho w^2$$

(2)

Energy balance equation:

$$\frac{\partial \rho h}{\partial t}+\frac{\partial \rho wh}{\partial z}=N_{tpc}^{\prime }{f}_q\left(z\right)$$

(3)

The equation of the state is:

$$\mathsf{\rho }=\mathsf{\rho }\left(h,P\right)$$

(4)

Since the above-mentioned supercritical fluids’ properties experience drastic changes [4,5,9,10,11,12,13,14,15,16,17], NIST steam table data (for steady-state solution) and appropriate approximation function for density, enthalpy, and velocity are used to capture the accurate variation as follows:

$${\rho }_{i,j}\left(t\right)=\frac{1}{\frac{1}{{\rho }_{i,j-1}\left(t\right)}+b_{i,j}\left(t\right) \left(z_{i,j}\left(t\right)-z_{i,j-1}\left(t\right)\right)}$$

(5)

$$h_{i,j}\left(t\right)=h_{i,j-1}\left(t\right)+a_{i,j}\left(t\right)\left(z_{i,j}\left(t\right)-z_{i,j-1}\left(t\right)\right)$$

(6)

$$w_{i,j}\left(t\right)=w_{i,j-1}\left(t\right)+D_{i,j}{N}_{tpc}\left(z_{i,j}\left(t\right)-z_{i,j-1}\left(t\right)\right)$$

(7)

where i = number of node; j = number of channels; a_i,j, b_i,j are phase variables; and D_i,j is a constant term (defined in Appendix A). These approximation functions are used in basic PDEs (Equations (1)–(3)) and applied to the weighted residual function, which transforms the PDEs to corresponding ordinary differential equations (ODEs). The detailed description is as follows:

• First ODE for the phase variable $a_{i,j}$ from the energy balance equation:

The energy balance equation is multiplied by the weighted function ${\Psi }_k$ and the weighted residual method is applied on the $i^{th}$ a node of a channel with limits z = $z_{i,j-1}\left(t\right)$ to z =$z_{i,j}\left(t\right)$ as:

$${{\displaystyle \int }}_{z_{i,j-1}\left(t\right)}^{z_{i,j}\left(t\right)}\left(M{h}_{i,j}\left(t\right)-R\right){\Psi }_k\left(z\right)dz=0$$

(8)

where M = $\left(\frac{\partial }{\partial t}+\frac{\partial w}{\partial z}\right){\rho }_{i,j}\left(t\right)$ and R = $N_{tpc}^{\prime }$. The weighted function is chosen as a unit and leads to the first ODE for the system as:

$$\frac{d{a}_{i,j}\left(t\right)}{dt}={\chi }_{i,j}\left(\overrightarrow{X},P\right)$$

(9)

• Second ODE for the phase variable $b_i$ from the energy balance equation:

Similarly, for the second phase variable, the mass balance equation is used and integrated over the $i^{th}$ node, and the equation is further simplified to obtain an ODE:

$$\frac{d{b}_{i,j}\left(t\right)}{dt}={{\rm Y}}_{i,j}\left(\overrightarrow{X},P\right)$$

(10)

• Final ODEs for the inlet velocity of each channel:

The momentum equation is integrated on each node to calculate all pressure drop components along the channels as follows:

$${{\displaystyle \int }}_{z_{i-1}\left(t\right)}^{z_i\left(t\right)}\Delta P_{i,j}dz={{\displaystyle \int }}_{z_{i-1}\left(t\right)}^{z_i\left(t\right)}{\left(\Delta P_{acc}+C\frac{d{w}_{in}\left(t\right)}{dt}+\Delta P_{fri}+\Delta P_{grav}\right)}_{i,j}dz$$

(11)

The inlet and outlet pressure drop of the parallel channels are as follows:

$$\Delta P_{k_{in},j}={\left(K_{in}{\rho }_{in}{w}_{in}^2\right)}_j$$

(12)

$$\Delta P_{k_{exit},j}={\left(K_{exit}{\rho }_{exit}{w}_{exit}^2\right)}_j$$

(13)

As discussed above, a previous single-channel methodology is adopted, and based on node convergence, the channels are nodalized into six nodes [29]. All pressure drop coefficients in each node are calculated (Equations (11)–(13)) and equated to the applied external pressure drop:

$${\displaystyle \sum }_{i=0,j=1}^{i=6,j=2}\Delta P_{i,j}+\Delta P_{k_{in}}+\Delta P_{k_{exit}}=\Delta P_{j,ext}$$

(14)

The final ODEs for the $w_{1in}\left(t\right)$ and $w_{2in}\left(t\right)$ are obtained using two boundary conditions for the parallel channel system as follows:

As the lower and upper plenum connect both channels, the applied pressure drop across both channels is as follows:

$$\Delta P_{1,ext}=\Delta P_{2,ext}\dots =\Delta P_{ext}$$

(15)

The total mass flow rate in can be calculated as follows:

$$w_{total}=w_{1,in}+w_{2,in}\dots +w_{N,in}$$

(16)

Using the above boundary conditions with Equations (14) and re-arranging the term leads to the final ODE for the $w_{1in}\left(t\right)$ as follows:

$$\frac{d{w}_{1in}\left(t\right)}{dt}=f_{i,j}(\overrightarrow{X},\overrightarrow{P})$$

(17)

$$\frac{d{w}_{2,in}\left(t\right)}{dt}=-\frac{d{w}_{1,in}\left(t\right)}{dt}$$

(18)

where ${\chi }_{i,j}(\overrightarrow{X},P),{{\rm Y}}_{i,j}(\overrightarrow{X},P),f_{i,j}(\overrightarrow{X},P),{\left(\Delta P_{acc},\Delta P_{grav},\Delta P_{fri}\right)}_{i,j},$ and $\Delta P_{exit}$ are defined in Appendix A. In the present analysis, a double channel is used; therefore, a total of 26 (2 equations from each node and 1 from momentum) coupled ODEs are used to examine the non-linear stability characteristics in a parallel channel system under supercritical fluid flow conditions.

3. Stability Analysis

Stability analysis predicts the dynamical system behavior when it experiences some disturbances in its initial operating conditions. To obtain the stability boundary, developed time-dependent coupled non-linear ODEs (Equations (9), (10), (17), and (18)) are linearized around the steady state by neglecting high-order terms in Taylor series expansion. The Jacobian matrix is created from these linearized equations, and the corresponding eigenvalues are obtained:

$$J\overrightarrow{x}=\lambda \overrightarrow{x}$$

(19)

$$J\overrightarrow{x}-\lambda \overrightarrow{x}=0$$

(20)

$$\left|J-\lambda I\right|=0$$

(21)

where $I$ and $\lambda$ represent the identity matrix and the eigenvalues of matrix J, respectively:

$${\lambda }_{}={\lambda }_R+i{\lambda }_I$$

(22)

After the application of a small perturbation, if all eigenvalues are ${\lambda }_R\le 0,$ the system is stable (decaying oscillations); however, if any ${\lambda }_R>0,$ the system displays unstable (growing oscillations) characteristics. The separating boundary between stable and unstable regions is known as a stability boundary, where all eigenvalues should be negative or purely imaginary (zero real part). Hence, a constant amplitude periodic oscillation would be observed.

The above-mentioned heated channel configuration is used in different energy systems, mainly in solar thermal applications, nuclear reactors, heat exchangers, etc., where the channel-to-channel interaction, varied heat flux (shading conditions), and flow rate affect the overall heat transfer performance of the system. Therefore, to study these effects on the non-linear dynamics of the system, two parallel heated channels are examined under different conditions.

3.1. Case I: Channels under Different Fluid (Supercritical $H_{2}0$ and $C{O}_2$) Flow Conditions

From the literature survey, it was observed that in several energy systems that are designed for supercritical cycles, supercritical water and carbon dioxide are commonly used as a working fluid. Both fluids work in different applications and operation conditions (pressure and temperature range, etc.). Therefore, NIST thermo-physical properties, as mentioned in

Table 1. The geometrical and operational parametric values that are assumed in the present study [9,29].

Property	${\mathbf{H}}_2\mathbf{O}$	$\mathbf{C}{\mathbf{O}}_2$	Units
System pressure	25	8	MPa
$P_c^{*}$	22.0640	7.3773	MPa
$L_H^{*}$	4.2672	4.2672	m
$D_H^{*}$	0.0034	0.0034	m
$T_{pc}^{*}$	373.95	307.82	$℃$
${\rho }_{pc}^{*}$	317.03	459.49	$\mathrm{kg}/{\mathrm{m}}^3$
$h_{pc}^{*}$	2152.54	3414.48	$\mathrm{kJ}/\mathrm{kg}$
${\beta }_{pc}^{*}$	0.129	0.299	$1/\mathrm{K}$
$C_{p,pc}^{*}$	76.445	35.267	$\mathrm{kJ}/\mathrm{kg}-\mathrm{K}$
$w_o^{*}$	1	1	m/s

, are used to capture the real fluid dynamics.

It can be evidently seen in Figure 2 that the obtained stability threshold is almost identical for both working fluids at lower subcooling numbers. In contrast, at high subcooling numbers, the dynamics change drastically. These static instability characteristics are observed only for supercritical water systems.

The Ledinegg instability characteristics are observed for the parallel channel similar to the single channel. Detailed analysis of the interaction between dynamic instability and static instability was reported in a previous study [29]. Co-dimension analysis was performed in a reported work to capture high non-linear stability phenomena via bifurcation analysis.

3.2. Case II: Channels under Equal and Unequal Heat Flux Conditions

The thermodynamic performance is highly dependent on the applied or generated heat to achieve the overall stable performance or higher efficiency. The heat distribution should be equally distributed. However, due to some internal (air voids or defects) or external disturbances (shading), unequal heat distribution and thermal fatigue occurs, and affects the thermal efficiency, which may lead to physical degradation. Therefore, this section examines parallel channel systems under equal and unequal heat flux conditions. To fulfil this objective, a new coefficient heat flux distribution factor ($h_{fd})$ is introduced as follows:

$$N{1}_{tpc}=h_{fd}N{2}_{tpc}$$

(23)

The above-mentioned equation is used with the developed ODEs (Equations (9), (10), (17) and (18)) to obtain the stability characteristics, drawn at different $h_{fd}$ values in the system operating parameter ($N_{tpc}-N_{spc}$) space as shown in Figure 3.

It is observed from Figure 3 that under an applied equal heat flux condition on both parallel channels, the obtained stability characteristics are almost similar to the single heated channel system [28]. Under different heat flux conditions, the obtained stability characteristics significantly influence the stability threshold as the hottest channel misbalances the heat transfer characteristics and makes the system unstable. Therefore, it can be observed form Figure 3 and Figure 4 that the stability boundary is shifted on the left side, making the system more unstable, as the heat flux ratio in both channels increases. Therefore, for a detailed non-linear analysis, an unequal heat flux condition model is used.

3.3. Case III: Channels under Different Inclination Conditions

Especially in solar energy applications, solar receivers or heated channels are primarily used in tilted configurations to capture maximum incoming solar irradiance. Therefore, this section is devoted to studying the inclination effect on the non-linear dynamics of parallel channels (Figure 4). The tilted angle, as the theta ($\theta$) term (from the horizontal) in the gravitational pressure drop in the momentum equation (Equation (3)), is included as follows:

$$\frac{\partial \rho w}{\partial t}+\frac{\partial \rho w^2}{\partial z}+\frac{\partial p}{\partial z}=-\frac{\rho sin\theta }{Fr}-\left[\Lambda +K_{in}{\delta }_d\left(z\right)+K_{exit}{\delta }_d\left(z-1\right)\right]\rho w^2$$

(24)

It is observed from Figure 5 that at a lower pseudo subcooling number (lighter fluid), the stability threshold remains almost the same. However, at a higher pseudo subcooling number (heavier fluid), it is significantly affected. In the supercritical cycles, the system pressure and temperature are observed to be significantly high; therefore, below the pseudo-critical region, the fluid is too dense whereas above the pseudo-critical region without phase change, the fluid becomes lighter, and can be treated similarly to ideal gas. Therefore, the gravitational pressure drop contribution to the overall pressure drop is much higher at a low subcooling number, in comparison to a high subcooling number, as shown in Figure 5.

3.4. Case IV: Different Flow Rates in a Channel

As aforementioned, due to some internal or external thermal incidents, the flow rate could be restricted or unequal. This can affect the overall thermal hydraulic performance of the system. This section examines the parallel channel system under different flow rate conditions in the heated channels. Both channels are simulated under the same flow rate conditions as the first condition whereas in the second condition, the flow rate is kept fixed for one channel while it is varied in the second channel with respect to heat. In Figure 6, it is evidently observed that the stability threshold is considerably affected as it shifts into the left-hand side (more unstable region) and makes the overall system more unstable. This reduces the stable parametric zone, which is the safe operation zone.

After examining the parallel heated channels under different aspect and working conditions, an unequal heat flux model is chosen to carry out detailed non-linear stability analysis in a widespread parametric space.

3.5. Validation and Comparison

The previously developed ROM is used to study the non-linear dynamics in parallel channel systems. In the previous study, spatial or temporal convergence was studied [29]; therefore, based on the results, the same six node nodalization scheme is executed in the parallel channel code. The present developed model is validated and compared with experimental studies [12,13] and numerical studies [18,31]. In Figure 7, line data are plotted at specific parametric values, whereas marker data are plotted without considering the inlet and exit loss coefficients.

Owing to the lack of experimental data, the present model is modified to establish a comparison with existing data. It can be seen from Figure 7 that the stability boundary observed in the present work follows the same trend observed in the literature data. It is also observed that the range of operating parameters is almost similar.

3.6. Non-Linear Analysis

As per the bifurcation theory, system behavior can be qualitatively and quantitatively changed, as some parameter values cross the critical threshold [7,32]. Therefore, linear stability analysis, which is only limited to a small perturbation, shows only the local behavior of the system. In case of any kind of internal or external perturbation of these critical parameters, the linear stability is lost, and a family of limit cycles bifurcates from this point. This is well known as a non-linear stability phenomenon. As mentioned above, coupled non-linear phenomenon are linearized as higher-order terms and are neglected in Taylor’s series expression. Therefore, some realistic characteristics are not captured through linear stability analysis. Henceforth, in non-linear analysis, these terms are considered and perturbed to correctly predict the system’s stability in a larger parametric space.

For this reason, detailed non-linear stability analysis is performed by execution co-dimensional analysis in a large parametric space using MATCONT bifurcation computational tools [30]. These tools used the in-built ODEs solver (ODEs45) to numerically simulate the developed ROM (Equations (9), (10), (17), and (18)).

Based on the eigenvalues’ nature, when the maximum critical pair of complex eigenvalues crossed the imaginary axis, it is termed Poincare–Andronov–Hopf or Hopf bifurcation. This bifurcation is categorized as co-dimension-1 bifurcation. Therefore, it is detected by changing only one free parameter at the moment. After detecting the Hopf bifurcation point and considering it as an initial solution, the non-linear stability threshold for both supercritical fluids is plotted in the ${\mathrm{N}}_{\mathrm{tpc}}-{\mathrm{N}}_{\mathrm{spc}}$ parametric space (Figure 8 and Figure 9) at unequal flux $h_{fd}=1.1$.

The Hopf bifurcation point associated with the dynamic instability characteristics appears in the system when the periodic solutions exist [28,32,33,34]. Therefore, numerical simulations are carried out on small and reasonably large perturbations around the Hopf bifurcation boundary at different parametric locations to observe the periodic characteristics. At the point on the stable side of the boundary, decaying oscillations with respect to time are observed as shown in Figure 10a, which indicates that the system has returned to its original equilibrium point. Whereas, on the unstable side, growing oscillations (shown in Figure 10c) are observed, confirming that the system is unstable. At the Hopf bifurcation point or stability boundary, the eigenvalues are purely imaginary; therefore, constant amplitude oscillations are observed in the system, as shown in Figure 10b. Almost similar oscillations characteristics are observed along the Hopf bifurcation boundary. Hence, only selected numerical simulations are presented to avoid repetition of the results.

3.7. Generalized Hopf Bifurcation

Apart from the eigenvalues, non-linear mathematical coefficients are also calculated. Herein, the first Lyapunov coefficient is one of the most critical parameters found in Hopf bifurcation analysis since it distinguishes the nature of stable and unstable limit cycles [7,32,34]. The first Lyapunov coefficient (FLC) is calculated, corresponding to Figure 8 and Figure 9, and plotted with the pseudo-subcooling number ($N_{spc}$) in Figure 11 and Figure 12. It is observed that the FLC value changes twice and thrice between positive to negative or vice versa for the SCW and SC-$C{O}_2$ system, respectively. The origin and disappearance of these limit cycles are known as generalized Hopf (GH) bifurcation points, where the FLC coefficient value is zero. Two and three generalized Hopf bifurcation points are observed along with bifurcation stability in the ($N_{tpc}$ − $N_{spc}$) plan, as shown in Figure 8 and Figure 9. It is observed that for different flow rate conditions, the GH point location is significantly changed for these cases.

3.8. Subcritical and Supercritical Hopf Bifurcation

As per the definition of the first Lyapunov coefficient, the nature if the limit cycles changes with its sign. The stable and unstable limit cycles are observed around the unstable (negative sign) and stable (positive sign) fixed points. Based on the trajectory behavior, the stable limit cycle works as an attractor whereas the unstable limit cycle works as a repeller. These non-linear dynamics are identified as supercritical and sub-critical Hopf bifurcation, respectively. The detailed dynamics of this phenomenon can be found in several studies [6,24,25,26,33,34].

As aforementioned, the bifurcation characteristic correspondingly changes its nature multiple times between supercritical and subcritical Hopf bifurcation. Therefore, to confirm this strange behavior, several numerical simulations are performed to demonstrate the appearance of limit cycles at different parameter values under similar flow rate conditions.

As shown in Figure 12, growing oscillations are observed and settle into a large constant amplitude, which confirms the existence of stable limit cycles on the unstable side, before $G{H}_1$ and between $G{H}_2-G{H}_3$. On the other hand, decay oscillations are observed initially between $G{H}_1-G{H}_2$ for small perturbations, which move into growing oscillations when the system quantitatively observes large perturbations. This confirms the existence of unstable limit cycles on the stable side of the boundary, as shown in Figure 13.

Similarly, numerical simulation is performed at different parametric locations to confirm both fluids’ other supercritical and subcritical Hopf regions in between other GH points. However, to avoid repetition, other numerical simulations are not presented. It should be noted that the supercritical region is identified as a safe bifurcation, whereas the subcritical region is considered more dangerous for normal system operations.

3.9. Bogdanov–Takens bifurcation

The Bogdanov–Takens (BT) bifurcation is categorized as co-dimension 2 bifurcation that appears in any system with ≥2 dimensions. The present system has 2N dimensions ($a_{i.j},b_{i,j};i=\mathrm{1..6},j=1,2$). The BT point represents the set of two operating parameters at which the linear stability boundary or Hopf bifurcation originates or terminates. Therefore, the observed BT point confirms the interaction of dynamic (DWOs) and static (Ledinegg) instability and vice versa along the stability boundary (Figure 8). For the equal flow-rate condition, the first to second and second to the third region of the stability boundary intersect by the BT point whereas for unequal flow rate conditions, a single BT point is observed, showing that the Hopf bifurcation characteristics are transformed into a saddle node (limit point) curve or vice versa.

In geometrical terms, BT bifurcation is a type where any autonomous ODEs have at-least two zero eigenvalues. Therefore, the system exhibits two equilibria: non-saddle and saddle. The Poincare–Andronov–Hopf bifurcation, which exists as a non-saddle equilibrium, continues to generate a limit cycle unidirectional to the BT point. In the other direction, the system becomes unstable. Detailed analysis of the Ledinegg instability region can be found in a previously reported work [28].

3.10. Wall Heat Effect

In this section, the temperature distribution profile is studied to investigate the wall heat effect, as shown in Figure 14. The linear temperature profile is considered without any heat loss (100% heat transfer rate). The wall temperature depends on the material properties and the wall heat transfer coefficient. Several authors conducted numerical and experimental studies to observe the wall heat effect on the system stability [19,20,21,22,35] and concluded that the wall heat effect plays a significant role in the system stability characteristics. Similar dynamics are observed in Figure 15, where stability boundary maps are obtained at different wall heat absorption rates.

$$T_{total}\left(q^{″}\right)=T_{wall}+T_{fuel}$$

(25)

4. Conclusions

The present work captured the thermal-hydraulic dynamics of a parallel channel system using a novel reduced order model (ROM) to convey different aspects of the non-linear characteristics observed in several energy systems, where the overall performance of the system depends on externalities. In this context, a parallel heated channel system was examined under different supercritical working fluids, different inclination conditions, varied heat flux, and flow rate conditions to examine the non-linear instability characteristics of the supercritical fluid flow system. The major observations and findings are detailed as follows:

• Different working fluid conditions, namely, supercritical water and carbon dioxide, were considered. Thermodynamic properties show significant changes near the pseudo-critical temperature; therefore, the NIST fluid properties and appropriate approximation functions (density, enthalpy, and velocity) were used to capture the real fluid dynamics. For both types of fluids, the stability characteristics were marginally different, especially at the high pseudo subcooling number ($N_{spc}$). The Ledinegg instability region was observed in between the dynamic instability region. This interaction was observed only for supercritical water through Bogdanov–Takens bifurcation (BT point). In comparison, only Hopf bifurcation was observed for the supercritical carbon dioxide. Similar stability characteristics were previously reported for a single heated channel system.
• The heat flux distribution plays a crucial role, and in order to apply the unequal heat condition on each channel, a heat distribution coefficient ($h_{fd}$ ) was introduced. It was observed that as the heat flux ratio of both channels increased, the stability threshold shifted towards making the overall system more unstable.
• The channels used in several energy systems are often oriented at an angle. Henceforth, the effect of the inclination on the stability characteristics was analyzed. It was observed that at a high subcooling number (heavy fluid), the major component of the pressure drop is due to gravitational force, leading to the dominant inclination effect.
• The parallel channel system was also examined under different flow rate conditions in respective channels at a fixed total mass flow rate where channel 1 counterbalances the flow rate in channel 2. Here, out-of-phase oscillations were observed. Additionally, when the flow rate was fixed in channel 1 and variable in channel 2, oscillation characteristics were observed only in the latter channel.

The present parallel channel model was also validated and compared with the experimental and numerical data. Furthermore, an unequal heat flux model was exploited to capture the non-linear stability dynamics of the parallel channel system. Herein, the observed stability threshold represents the dynamic and static instabilities and several bifurcation characteristics. In dynamic instability, Hopf bifurcation characteristics were observed, which confirmed the existence of oscillatory behavior and the presence of limit cycles. It was observed that the stable limit cycle acts as an attractor on the unstable side, whereas the unstable limit cycle acts as a repeller on the stable side. This unstable limit cycle is more dangerous for safe system operation. Henceforth, sub-critical and supercritical Hopf bifurcations were confirmed through numerical simulation and by calculating the first Lyapunov coefficient. In static instability, saddle node bifurcation was observed, which is associated with the Ledinegg excursive phenomena. Additionally, the Bogdanov–Takens bifurcation points were observed as an interaction of the saddle node with Hopf bifurcation. A co-dimension analysis was performed in order to capture the bifurcation characteristics along with several numerical simulations. The findings consist of different types of bifurcation phenomena, representing various non-linear dynamics features. This analysis helps to enhance the heat transfer characteristics in energy systems, such as nuclear reactors, solar thermal systems, boilers, heat exchangers, etc., which are instrumental in their safety and design.