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

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.


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].
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.

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. 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: 1.
Homogenous flow is considered in both channels.
The system pressure remains constant to use the thermodynamics property of the fluids.
the entire system dynamics as it only predicts tiny or small perturbation impacts system dynamics. Thereby, this reveals the importance of non-linear studies in dyn systems [6][7][8][24][25][26][27][28], such as two-phase flow, boiling water reactors, advanced heav reactors, heated channel systems, etc., for the normal and safe operation of such syste The present work exploits a previously developed methodology (nodalized r ordered model) for a single channel [29] to a parallel channel system. An ROM wa oped to capture the non-linear characteristic of a parallel channel under different critical fluid (water and carbon dioxide) flow conditions. Apart from the in-hous MATCONT software (bifurcation tool) was used [30], which works on a MATLA form. In order to demonstrate realistic engineering applications, the system was under different working conditions, such as an equal and un-equal heat flux distr and varied flow rate. The major findings are the several bifurcation characteris tained along the stability boundary, depicting the system's non-linear dynamical ior. The model was validated using numerical and experimental literature data. T sent study is instrumental in the design of the safe operation of nuclear reactors an thermal systems, in a wide parametric range.

Modelling Methodology
The objective of the present work is fulfilled using a system with two paralle nels. The physical schematic is shown in Figure 1, and the corresponding design a erating parameters are shown in Table A1. This physical system is replicated in th ematical form using the 1-D Navier-Stokes equations [4][5][6]9,10]. These conservatio tions are further simplified using non-dimensional parameters [11,29] under some priate assumptions as follows: 1. Homogenous flow is considered in both channels. 2. The inlet temperature remains constant. 3. The system pressure remains constant to use the thermodynamics property fluids.
Energy balance equation: ∂ρh ∂t The equation of the state is: 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: 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 Ψ k and the weighted residual method is applied on the i th a node of a channel with limits z = z i,j−1 (t) to z = z i,j (t) as: where M = ∂ ∂t + ∂w ∂z ρ i,j (t) and R = N tpc . The weighted function is chosen as a unit and leads to the first ODE for the system as: • 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: • 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: The inlet and outlet pressure drop of the parallel channels are as follows: 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: The final ODEs for the w 1in (t) and w 2in (t) 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: The total mass flow rate in can be calculated as follows: Using the above boundary conditions with Equations (14) and re-arranging the term leads to the final ODE for the w 1in (t) as follows: where , and ∆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.

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: where I and λ represent the identity matrix and the eigenvalues of matrix J, respectively: After the application of a small perturbation, if all eigenvalues are λ R ≤ 0, the system is stable (decaying oscillations); however, if any λ 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.

Case I: Channels under Different Fluid (Supercritical H 2 0 and CO 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, 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. The Ledinegg instability characteristics are observed for the parallel chan to the single channel. Detailed analysis of the interaction between dynamic inst static instability was reported in a previous study [29]. Co-dimension analys formed in a reported work to capture high non-linear stability phenomena via analysis.

Case II: Channels under Equal and Unequal Heat Flux Conditions
The thermodynamic performance is highly dependent on the applied or heat to achieve the overall stable performance or higher efficiency. The heat d should be equally distributed. However, due to some internal (air voids or external disturbances (shading), unequal heat distribution and thermal fatig and affects the thermal efficiency, which may lead to physical degradation. The section examines parallel channel systems under equal and unequal heat flux To fulfil this objective, a new coefficient heat flux distribution factor (ℎ ) is as follows: The above-mentioned equation is used with the developed ODEs (Equatio (17) and (18)) to obtain the stability characteristics, drawn at different ℎ va system operating parameter ( − ) space as shown in Figure 3.

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 f d ) is introduced as follows: 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 f d 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 Figures 3 and 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.  It is observed from Figure 3 that under an applied equal heat flux condition allel channels, the obtained stability characteristics are almost similar to the s channel system [28]. Under different heat flux conditions, the obtained stability ch significantly influence the stability threshold as the hottest channel misbalances th fer characteristics and makes the system unstable. Therefore, it can be observed  Energies 2022, 15, x FOR PEER REVIEW as the heat flux ratio in both channels increases. Therefore, for a detailed nonan unequal heat flux condition model is used.

Case III: Channels under Different Inclination Conditions
Especially in solar energy applications, solar receivers or heated chann ily used in tilted configurations to capture maximum incoming solar irra fore, this section is devoted to studying the inclination effect on the non-li of parallel channels ( Figure 4). The tilted angle, as the theta ( ) term (from in the gravitational pressure drop in the momentum equation (Equation (3 as follows: It is observed from Figure 5 that at a lower pseudo subcooling number the stability threshold remains almost the same. However, at a higher pseu number (heavier fluid), it is significantly affected. In the supercritical cyc pressure and temperature are observed to be significantly high; therefo pseudo-critical region, the fluid is too dense whereas above the pseudo without phase change, the fluid becomes lighter, and can be treated similar Therefore, the gravitational pressure drop contribution to the overall pr much higher at a low subcooling number, in comparison to a high subcool shown in Figure 5.

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 (θ) term (from the horizontal) in the gravitational pressure drop in the momentum equation (Equation (3)), is included as follows: ∂ρw ∂t 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 pseudocritical 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.

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

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 w conditions, an unequal heat flux model is chosen to carry out detailed non-linear s analysis in a widespread parametric space.

Validation and Comparison
The previously developed ROM is used to study the non-linear dynamics in p channel systems. In the previous study, spatial or temporal convergence was studi

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.

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 codimensional 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 N tpc − N spc parametric space (Figures 8 and 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.

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 codimensional 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 N tpc − N spc parametric space (Figures 8 and 9) at unequal flux h f d = 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.   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. 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.

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 Figures 8 and 9, and plotted with the pseudo-subcooling number (N spc ) in Figures 11 and 12. It is observed that the FLC value changes twice and thrice between positive to negative or vice versa for the SCW and SC-CO 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 Figures 8 and 9. It is observed that for different flow rate conditions, the GH point location is significantly changed for these cases.

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 Figures 8  and 9, and plotted with the pseudo-subcooling number ( ) in Figures 11 and 12. It is observed that the FLC value changes twice and thrice between positive to negative or vice versa for the SCW and SC-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 ( − ) plan, as shown in Figures 8 and 9. It is observed that for different flow rate conditions, the GH point location is significantly changed for these cases.

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 1 and between 2 − 3 . On the other hand, decay oscillations are observed initially between 1 − 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.

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 GH 1 and between GH 2 − GH 3 . On the other hand, decay oscillations are observed initially between GH 1 − GH 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.

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 ( . , , ; = 1. .6, = 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 flowrate 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 atleast 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].

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.

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 = 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].

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.

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 nonlinear 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

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 nonlinear 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

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 f d ) 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.

Appendix A
In the present work, the complete study is performed in a dimensionless form. In order to non-dimensionalize the set of conservation Equations (1)-(4), the pseudo-critical temperature point is considered as a threshold value, as follows [14]: Table A1. The list of non-dimensional parameters used in the present study.