- Email: [email protected]
Interaction of density wave oscillations and flow maldistribution for two-phase flow boiling parallel channels
International Journal of Thermal Sciences 145 (2019) 106026
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Interaction of density wave oscillations and flow maldistribution for twophase flow boiling parallel channels
T
Deepraj Paula,b, Suneet Singha,∗, Surendra Mishrab a b
Department of Energy Science and Engineering, Indian Institute of Technology-Bombay, Mumbai, India Nuclear Power Corp. of India Ltd, Mumbai, India
ARTICLE INFO
ABSTRACT
Keywords: Thermal hydraulics Stability Flow maldistribution Density wave oscillations Bifurcation
Density wave oscillations (DWO) and flow maldistribution (FMD) are two important instabilities among others in two-phase flow boiling systems. As both the phenomena can lead to undesirable consequences such as channel burnout, the design intent is to avoid the same. DWO is caused by the interaction between single-phase and twophase region as perturbation propagates slowly in the later region. Whereas in case of FMD, it leads to one channel receiving higher flowrate and other channel receiving lower flowrate for a twin parallel channels system under forced flow condition. Thus, channels subjected to similar thermal hydraulic and geometric condition receive non-identical flowrate. Bifurcation analysis is used to identify DWO and FMD instabilities in this work. DWO is identified by Hopf bifurcation, thus Hopf curve which constitutes Hopf bifurcations represents the DWO boundary. FMD has been observed earlier through numerical simulation but the region which depicts the same in a parametric space has not been demarcated in previous studies. In this work, FMD is identified by a pitchfork bifurcation. This leads to straight-forward identification of FMD boundary as a pitchfork curve. A combined analysis of multiple instabilities using a single stability map can demarcate regions with the dominance of particular instability. It can thus provide a better perspective for system design and system instabilities in terms of suitable operating region selection. A stability map which can capture and demarcate both the phenomena simultaneously i.e, DWO and FMD in a parametric space is the interest of present work. The boundary of DWO and FMD splits the parametric space into three distinct regions. These are DWO region, a region of symmetrical stable solutions and FMD region which contains non-identical (asymmetric) stable solutions. This breakdown of symmetry of the system is observed by non-identical flowrate in the channels. In some cases, boundary depicting DWO and FMD intersects. Hence, a fourth region is also observed which possess the characteristics of both the instabilities. The stability map shows the transition between these regions.
1. Introduction Flow instabilities can occur within a two-phase flow boiling heated channels. The best approach to deal with such phenomena is to design the system in such a way which eventually avoids the same. A thermal hydraulics model is based on a set of conservation equations and boundary conditions. In literature, several models have been proposed for the analysis of two-phase flow behaviour. The common approach to model the two-phase flow is by using Homogeneous equilibrium model (HEM). Other models such as drift-flux model by Rizwan-Uddin & Dorning [1]; separated flow model, slip flow model, stratified flow model etc. are also used. In the homogeneous flow model, it is assumed that velocity, temperature and pressure between the liquid and vapour phases are equal. The drift-flux model improves the homogeneous flow
∗
model, by considering slip velocities between the two phases. Whereas, a two-fluid approach considers different sets of conservation equations of each phases and their interactions. This signifies more details towards the mathematical modelling. Moreover from the system analysis perspective, optimization between model order and level of detailing needs to be analyzed first for reducing computational expense. The most popular classification, introduced in Boure et al. [2] divides two-phase flow instabilities in static and dynamic instabilities. For static instabilities, the threshold of unstable behaviour can be predicted from the steady-state conservation laws. On the other hand, to describe the behaviour of dynamic instabilities it is necessary to take into account different dynamic effects, such as the propagation time, inertia, compressibility, etc. In addition, the term compound instability is normally used when several of the basic mechanisms interact with each
Corresponding author. E-mail address: [email protected] (S. Singh).
https://doi.org/10.1016/j.ijthermalsci.2019.106026 Received 30 July 2018; Received in revised form 29 June 2019; Accepted 12 July 2019 1290-0729/ © 2019 Elsevier Masson SAS. All rights reserved.
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
other. A considerable amount of studies has been carried out in this field by Ruspini et al. [3]; Kakac & Bon [4] etc. A lot of 1D analyses has been carried by researchers for parallel channels viz. Lee & Pan [5]; Lu et al. [6]; Mishra & Singh [7] etc, using various numerical schemes. DWO for thermal hydraulics coupled with neutronics for a single channel was analyzed by Lee et al. [8]. The influence of void-reactivity feedback on the bifurcation phenomena and nonlinear characteristics of a single nuclear-coupled boiling channel in an advanced boiling water reactor was studied. As noted, DWO is the most widely studied dynamic instability which occurs due to the delayed feedback between single and two-phase regions. Since the perturbation propagates slowly along the two-phase region, a significant delay marks the onset of perturbations in two-phase region. Hence, the two-phase pressure drop and the single phase pressure drop oscillates out-of-phase. It is thus very important to identify the regions in the parametric space where such instabilities exists. Van Oevelen et al. [9] discussed that FMD can occur because of nonmonotonic nature of pressure drop vs flowrate curve in a heated channel, asymmetrical inlet header design, differences in channel geometries etc. His work proposes a new numerical approach to analyze the stability based on a generalized eigenvalue problem. This methodology is specifically designed for analyzing systems with hundreds of identical parallel channels. Due to the non-monotonicity of the pressure drop curve, there exist multiple solutions. For a single channel system, driven by a constant pressure boundary there are three solutions. For two channel system, as both the channels have the same pressure drop so the number of solutions are nine. Thus the number of solutions increases as the number of channels increase. Flow rate distribution among parallel evaporator tubes was studied by Akagawa et al. [10]. They experimentally obtained channel load curves for individual tubes experiencing flow boiling. Flow excursion was also observed, as the system may not be stable at all the intersection points of the internal and external pressure drop vs flowrate curve. Minzer et al. [11] experimentally observed the presence of maldistribution in heated channels for a forced flow system. Here, asymmetric but stable solutions exist, where one channel receives higher flowrate but another channel is subjected to lower flowrate. Minzer et al. [12] developed a model where each channel is represented by a single ordinary differential equation. Coupling between multiple channels of the thermal-hydraulic system was done using boundary condition. For 2 number of parallel channels, FMD was observed even though channels were subjected to the same heat flux. Baikin et al. [13] analyzed FMD in a system which comprises of 4 number of parallel channels. This lead to even more number of possible solutions. A combination of heated and unheated pipes in the system was analyzed in the work for FMD. Even for the case of all pipes being heated equally, there were ranges of inlet flowrate where severe FMD could occur. The model was further developed in the work of Taitel & Barnea [14] where the pipe was subdivided into multiple cells. The calculation of the pressure drop was done in each cell. The transient simulations were performed and was observed that initiation of the system at unstable fixed points leads to the divergence of trajectories which settles at an uneven flowrate being distributed among the channels. Thus, there is an asymmetric distribution of flowrate among channels. However, the model used in the above studies assumes the mass flowrate to be uniform along the channel (quasiequilibrium assumption). A perturbation propagates slowly in twophase region as compared to single-phase region, which means there is phase-difference between them. DWO occurs because of delayed feedbacks between the inertia of flow and compressibility of the two phase mixture. Thus with quasi-equilibrium assumption in the above approach, one cannot capture DWO dynamics. Pandey & Singh [15] studied the excursive instability and density wave oscillations (DWO) for two-phase flow in a natural circulation loop. The stability maps obtained through their work by using bifurcation analysis has three distinct regions, viz. the first region has only DWO, second region has Ledinegg as well as DWO, while the third
region has only Ledinegg instability. Using bifurcation analysis to identify various flow instabilities, Pandey & Singh [16] explored a new type of stability boundary for high sub-cooling number in a two-phase natural circulation loop. The region near to this boundary shows sustained oscillations for flow velocity; confirming that this is a dynamic instability like DWO. This stability boundary is near to Ledinegg instability region but does not show the characteristics of Ledinegg stability. Simanovskii et al. [17] studied the influence of a temperaturedependent inter-facial heat release/consumption on nonlinear oscillatory convective regimes in a two-layer system with rigid heat-insulated lateral walls under the action of an inclined temperature gradient. A wide range of parameter can lead to the change of the sequence of bifurcations and the development of new nonlinear regimes in the system. Luo et al. [18] studied the bifurcation phenomena and the existence of dual solutions in natural convection in eccentric annulus. The onset and evolution of pitchfork bifurcation were tracked by the lattice Boltzmann method (LBM) that was inherently transient and has a good flexibility in dealing with irregular geometries. Paruya et al. [19] analyzed DWO in several boiling channels with varying lengths (Froude number) by adopting a moving node scheme (MNS) and fixed node scheme (FNS). It was observed that MNS has a higher-order convergence and high computational efficiency compared to finite difference scheme and FNS. It was also demonstrated that the solutions bifurcate from a stable steady state to periodic solutions to chaotic attractor ones when the boiling channel was subjected to the periodic perturbation of pressure drop at different power level. Basu et al. [20] focused upon the development of an unified model for the analysis of single-phase natural circulation loops with constant heat flux heating and convective cooling. Characterizing parameters were defined such that they are applicable irrespective of any geometrical configuration. Karmakar & Paruya [21] reported the chaotic flow oscillations in a natural circulation boiling loop. It was observed that boiling regime changes from nucleate boiling to transition boiling as a result of decrease in inlet sub-cooling and increase in heat flux. The thermal oscillations were strongly coupled with pressure-drop oscillations. Nonlinear analysis of the time series of loop flow rate at various heater power and inlet sub-cooling were carried out using statistical analysis. Krishnani & Basu [22] focused on assessing the viability of Boussinesq approximation during nonlinear stability appraisal following transient simulation of single-phase natural circulation loops. Paul & Singh [23] studied the stability behaviour of various uniform and non-uniform axial heat flux. The study reveals that keeping the total heat rate same, the stability characteristics of the system significantly changes from uniform to non-uniform axial heat flux assumption. Xia et al. [24] studied the flow distribution and flow instability through numerical simulation at various parametric points in a stability map. Such an approach is a tedious one and is computationally expensive. Moreover, a clear demarcation of FMD in the parametric space was not shown. A combined analysis which captures DWO and FMD in a forced flow system has been lacking in the literature. Bifurcation analysis is used in this work for identification of various instabilities and their boundary in a parametric space. Unlike DWO which are identified by Hopf bifurcation, identification of FMD by the use of bifurcation analysis has not been done earlier. In this work FMD is identified by pitchfork bifurcation. Our focus also involves the identification of DWO and FMD boundary in the stability map. The two boundaries splits the parametric space in three regions. These are DWO region, a region of the symmetrical stable solutions and FMD region. For a different parametric combination, a fourth region can also be observed. A MATLAB based mathematical package named MATCONT developed by Dhooge et al. [25] is used in this work to detect and analyze DWO and FMD. MATCONT provides the means for continuing equilibria of systems of Ordinary Differential Equations (ODEs), and their bifurcations. Identification of DWO and FMD along with their boundaries are computed using the package. Numerical simulations are carried out to verify and understand the characteristics of the system 2
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
Fig. 2. Experimental data validation of stability boundary, Ns = 3,5 & 7 .
• Heat flux is uniform in the axial direction • Fixed inlet sub-cooling • Sub-cooled boiling is not considered here • Constant thermodynamic properties at system pressure as pressure
drop through the channel is much smaller than that of system pressure
Fig. 1. Single Channel System, Lee et al. [8].
The derivation of the mathematical model is provided in Appendix A. The thermal hydraulic model (A.30) is validated (Fig. 2) with the reevaluated experimental data of Freon-113 from Saha et al. [29] as given in Rizwan-Uddin & Dorning [30] & Paul & Singh [31] (Set I). Nsub and Npch (as described in the work of Ishii [32] and Saha [33]) are the two non-dimensional parameters used for validation. These are significant not only for stability analyses, but also for the description of steady state operational conditions. Nsub takes into the account of time lag effects in the liquid region due to sub-cooling of the fluid entering the heated duct. It scales the inlet sub-cooling and it is the dimensionless residence time in the single-phase region under the equilibrium assumption. Thus it is one of the important parameters for the stability analysis. Whereas, Npch scales the change of phase due to the heat transfer to the system. The single phase region is nodalized with 3, 5 & 7 nodes. It is observed that all the three different stability boundaries (corresponding to different nodalizations in the single-phase region) as calculated using the model almost overlap with each other. The calculated stability boundary matches reasonably with the boundary obtained from experimental data. For lower values of Nsub , there is some mismatch. Here, inlet enthalpy is closer to saturation enthalpy and may lead to superheated exit condition (following a transient). Such is beyond the scope of analysis of the present work. Our work focuses primarily on vapour quality, 0 < x < 1. A schematic representation of a parallel channel is shown in Fig. 3. For forced flow system, mathematical coupling between the channels is used as per the work of Lee & Pan [5];
around the stability boundaries. 2. Mathematical modelling and bifurcation 2.1. Model description A typical flow boiling channel is shown in Fig. 1. The liquid in subcooled condition enters the channel and heats up which enforces the change of phase of liquid. The point where the change of phase occurs i.e, the liquid starts to boil is called the boiling boundary. Beyond the boiling boundary, the temperature remains constant although an increase in enthalpy due to heat flux is accounted by positive vapour quality. The liquid thus leaves the channel as a two-phase mixture. 2.2. Model details Alobaid et al. [26] discussed various types of models with varied complexities viz. homogeneous equilibrium model (HEM), drift flux model and two-fluid separated flow model (which includes four-equation, five-equation, six-equation and seven-equation flow model) etc. HEM although simple to implement, it cannot capture the relative velocity between the two-phases. However, the uncertainty in modelling the interaction between phases at the inter-phase boundary poses considerable difficulties. These relations cannot be derived from fundamental physical laws and in most cases are based on empirical assumptions. Solving the resulting differential equations requires higher computational effort and entails parameters that may cause numerical instability, especially due to improper selection of interfacial terms. The difficulties associated with the two-fluid models can be significantly reduced by formulating the two-phase flow in terms of the mixture flow model. Increased computational expense, existence of uncertainty in modelling the interaction at the inter-phase boundary and mathematical complexities of multi equation two-fluid models drive us to the use Homogeneous equilibrium model (HEM). HEM is simple, computationally inexpensive and yet applicable to wide range of two-phase flow regimes. The details of the Homogenous Equilibrium Model used here can be found in the work of Lin et al. [27]; Lee et al. [8] and Clausse & Lahey Jr [28]. Assumptions for the flow boiling system are:
dui+, j dt +
= Aj
dui+,1 dt +
+ Bj , j = 2,3…M
here, j represents the channel number; sional inlet velocity of jth channel. and
dui+,1 dt +
+ dWTot
=
dt
1+
M j=2 M j=2
(1)
ui+, j
represents the non-dimen-
Ax+ s Bj
Ax+ s Aj
(2)
where, Ax+ s, j is the area ratio of jth channel to first channel. Thus, Ax+ s,1 = 1. When total mass flow rate is held constant for a forced flow system, 3
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
called co-dimension of a system. Codim-1 and Codim-2 are used to represent single and two parameters variation, respectively. For a dynamical system, the sign of real part of the largest eigenvalue characterizes the stability of the system. For a positive sign, a small perturbation leads to divergence of trajectory. For a negative sign, a small perturbation leads to convergence of trajectory towards a stable fixed point. Such analysis can provide us with a demarcation between the stable and unstable region in a parametric space. Nonlinear dynamics analysis or bifurcation analysis can identify the region in a parametric space which contains stable and unstable limit cycles. Hopf bifurcation, Generalized Hopf (GH) bifurcation and Pitchfork bifurcation are observed in this work and are discussed below. 2.3.1. Hopf bifurcation (density wave oscillations) Density Wave Oscillations are self-sustaining oscillations and are identified by a Hopf bifurcation. The Hopf bifurcation occurs when one pair of the complex eigenvalues become purely imaginary. The presence of complex eigenvalues signify oscillatory behaviour. The Hopf curve which constitutes Hopf bifurcation points represents the DWO boundary. The Hopf bifurcation is a codim-1 parametric bifurcation. Hopf curve is generated by codim-2 parametric bifurcation which constitutes the Hopf points. Evolution of limit cycles is the phenomena related to Hopf bifurcation. This type of Hopf bifurcation can be understood by the sign of the first Lyapunov coefficient (FLC). For a negative sign of FLC, it represents super-critical Hopf whereas a positive sign represents sub-critical Hopf. In sub-critical Hopf region, the presence of an unstable limit cycle in the linearly stable region is of prime concern. Here small perturbation to the system converges towards a stable fixed point and relatively large perturbation leads to diverging DWO. On the other hand, for super-critical Hopf bifurcation in a linearly unstable region, there is a presence of a stable limit cycle. It means, irrespective of the magnitude of perturbation, DWO converges towards a stable limit cycle.
Fig. 3. Parallel channel schematic.
+ WTot =
M
Wi+
(3)
j =1
where, m is the number of channels. So, + dWTot =0 + dt
dui+,1 dt +
=
(4)
1+
M A+ B j=2 x s j M A+ A j=2 x s j
2.3.2. Generalized Hopf (GH) bifurcation DWO boundary or Hopf boundary is obtained as codim-2 bifurcation. The curve splits the parametric space into two regions namely, stable DWO and unstable DWO regions. Generalized Hopf bifurcation or Bautin bifurcation is observed over the Hopf curve i.e, DWO boundary. Here, FLC is zero. GH point divides the DWO boundary in two regions namely sub-critical and super-critical Hopf regions characterized by positive and negative FLC respectively. It represents the transition from sub-critical Hopf to super-critical Hopf or vice-versa. As discussed earlier for sub-critical Hopf region, large perturbation leads to divergence of trajectories in linearly stable region of DWO boundary. So identification of GH point is very important.
(5)
where,
Aj = MH+,1/ MH+, j Bj =
(
PH+,1
(6)
PH+, j )
MH+, j
(7)
here, MH+, j is the non-dimensional channel mass of the jth channel (eqn (A.24)). PH+, j is the non-dimensional pressure drop of the jth channel (eqn (A.29)). Thus the model of 2 number of parallel channels can be represented as,
x = [Ln+,1
+ e,1
ui+,1 Ln+,2
+ e,2 ]
, n = 1…Ns
2.3.3. Pitchfork bifurcation (flow maldistribution) Flow excursion (in a single channel) occurs due to the presence of negative slope region in pressure drop vs flowrate curve. An operational point is stable when the slope of the internal characteristic curve is smaller than the slope of the external characteristic curve. In case of a multiple channels, connected by an inlet and exit plenum the pressure drop across them is same. So there exists multiple solutions as external pressure drop intersects the internal pressure drop multiple times. Naturally, due to the presence of negative slope region in pressure drop curve and system flowrate being constant (in a forced flow system) the operating point (if on the negative slope) shifts towards stable solutions. This leads to flow re-distribution in the system. Such multiple solutions can be identified by a pitchfork bifurcation. It is a codim-1 parametric bifurcation. In MATCONT, pitchfork bifurcation is identified using Branch Point (BP) bifurcation. BP marks the onset of FMD. As BP is observed in an equilibrium curve, it signifies the emancipation of other equilibrium curves from there. Thus BP is the influx/outflux of different equilibrium curves. Pitchfork bifurcation is of two types viz, sub-critical and super-critical. In the present work, sub-critical
(8)
The inlet velocity of the other channel is found using the boundary condition of total mass flowrate being constant. 2.3. Bifurcation analysis in context to flow instability The details related to bifurcation analysis and review of two-phase flow instabilities are available in Kuznetsov [34] and Ruspini et al. [3] respectively. However, a brief description of stability and bifurcation analysis for the identification flow instabilities and demarcation of flow stability boundaries are presented here. Bifurcation of a dynamical system refers to the quantitative changes in the dynamics of the system which is obtained by varying parameters. The parameter values at which they occur are called bifurcation points. The number of parameters which are varied for bifurcation to occur is 4
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
pitchfork bifurcation is observed and represents the presence of asymmetric solution but stable ones in the parallel channels. Symmetric solutions although present, are unstable. The stability of such equilibrium curves can be studied by eigenvalues. The locus of pitchfork bifurcation represents the FMD boundary. 3. Results and discussion This section is divided into three major parts. In the first part (section 3.1), the equilibrium points are generated using MATCONT by varying the non-dimensional inlet velocity us+ while remaining parameters are kept constant. These equilibrium points are represented by plotting boiling boundary ( +) as a function of us+ and various bifurcations are shown in these equilibrium curves. These are so-called codim-1 bifurcation curves. In the second part (section 3.2), the various type of stability boundaries are plotted in a parametric space spanned by Nsub and us+. Here, these two parameters are simultaneously varied while other parameters are kept constant. These are so-called codim-2 bifurcation curves. These curves divides the parameteric space into regions of different flow characteristics. In the third part (section 3.3), numerical simulations are carried out for different parametric values. The characteristics of the time evolution of state variables can be observed. Table 1 is used to generate the graphs for this study.
Fig. 4. Equilibrium curve corresponding to codim-1 parameter variation, boiling boundary ( +) vs us+ .
3.1. Equilibrium curves with variation in us+ Fig. 4 represents the variation of location of the boiling boundary ( +) with us+. As expected, the location of boiling boundary moves up with increase in us+ (represented by the solid lines). Moreover, it is observed that there is the presence of two bifurcation points namely, a BP (or pitchfork bifurcation) and a Hopf point. As us+ is increased, Hopf point is observed first and then a BP. From BP, additional equilibrium curves emanates (represented by dotted line) in addition to the curve represented by solid line. The channels have symmetric (solid lines) as well as asymmetric solutions (dotted lines). For symmetric part of the solution, both the channels have identical solution i.e, channels receive the same flowrate and have same boiling boundary. For asymmetric part of the solution, both the channels have different solutions viz. boiling boundaries, flowrate etc even though they are subjected to the same thermal hydraulic and geometric conditions. It is thus observed, that for a particular Npch value, there exist three different regions in an equilibrium curve. The first region is before the occurrence of Hopf point. The second region is in between the Hopf point and BP. The third region is beyond the occurrence of BP. For the green coloured equilibrium curve (Npch = 12 ), the regions are marked as X1 i.e, the first region, X2 i.e, the second region and X3 i.e, the third region. These three regions caters to the symmetrical solution of the system i.e, both the identical channels have same flowrate and boiling boundaries ( + ). The asymmetrical solution which occurs after BP, are marked as Y1/Y2. This marks the onset of FMD. (The marked points namely X1, X2, X3 and Y1/Y2 will be used later for numerical simulations, section 3.3.) To further analyze the solutions of these three regions marked as X1, X2 and X3 in Fig. 4, the real part of the largest eigenvalue along the equilibrium curve (of the solid line) are plotted in Fig. 5. For region X1, it is observed that as us+ is varied, for lower values before the observance of Hopf point, real part of the largest eigenvalue is positive.
Fig. 5. Variation of real part of the largest eigenvalue for symmetric part along the equilibrium solutions. (Truncated y-axis for proper representation of change of sign of largest eigenvalue).
This signifies unstable solutions. Numerical simulation is provided in section 3.3.1. As us+ is increased further, the sign of the eigenvalue changes from positive to negative as a Hopf point is encountered. This signifies stable solutions. This is marked as X2 region in Fig. 4. Numerical simulation is provided in section 3.3.2. Further increase in us+, BP is encountered on the equilibrium curve. The sign of the real part of the largest eigenvalue, Fig. 5 which caters the symmetric part of the solution changes to a positive value. This signifies the symmetric part of the solutions is unstable. This is marked as X3 region in Fig. 4. Numerical simulation is provided in section 3.3.3. Two set of equilibrium curves emanates (dotted lines) from the branch point (BP) in Fig. 4. These are marked as Y1/Y2. Each equilibrium curve corresponds to a channel which are complementary to each other. The real part of the largest eigenvalue vs non-dimensional boiling boundary (which itself is a function of velocity) are shown in Fig. 6. Negative sign of real part signifies stable solution. It can be observed from Fig. 4 that, as us+ is increased beyond BP, symmetry of the system breaks as these solutions have positive real parts of largest eigenvalue. Asymmetric solution exists and are stable, as observed by negative real part of eigenvalue. This confirms the presence of FMD. Such asymmetric solution leads to higher flowrate in one channel and lower flowrate in the other channel, which is undesirable. Numerical
Table 1 Parameters used for analysis. Pressure (MPa)
Ns
ui0 (m/s)
D (m)
L (m)
ki
ke
7
3
1.0
0.01
1
1
1
5
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
Fig. 6. Variation of real part of the largest eigenvalue for asymmetric part along the equilibrium solutions; (A)-Branch Y1 in Fig. 4; (B)-Branch Y2 in Fig. 4.
simulation is provided in section 3.3.4.
boundary is provided in section 3.3.5. Similarly, from the BP or pitchfork bifurcation, Nsub is varied simultaneously along with us+ as codim-2 parameters to generated a BP curve or the pitchfork curve i.e, FMD boundary. This is shown as black coloured curve in Fig. 7. The stability map thus generated is spanned by Nsub and us+. It is subdivided into three regions viz. DWO region (A), Symmetrical stable solution (B) and FMD region (C). Region A is unstable. Region B is the region of symmetrical stable solution. Region C is the region of FMD where asymmetric and symmetric solution co-exists. The symmetric solutions which exists in this region are not stable. Perturbation causes trajectories to diverge and settle to a stable solution. This causes one channel to receive higher flowrate and another channel receives lower flowrate. For different values of Npch , stability map for DWO and FMD are plotted in Fig. 8. It is observed that for higher values of Npch , the stability boundaries (i.e, the pair of DWO and FMD boundary) shifts right which demands higher value of us+ to attain the symmetrical stable solution region. An additional fourth region (D) can also be observed for a different parametric combination and is shown in Fig. 9. The presence of this additional region in the stability map is due to the overlap between the pitchfork curve (i.e, FMD) and the Hopf curve (i.e, DWO). Further discussion along with numerical simulation of this region is provided in section 3.3.6.
3.2. Stability boundary From the Hopf point which was found by varying us+ in Fig. 4, Nsub also is now varied simultaneously to generate the locus of Hopf boundary i.e, DWO boundary. This is codim-2 parametric stability boundary and is shown in Fig. 7. The blue curve represents the Hopf curve (or the DWO stability boundary), where each point on the curve is a Hopf point. This curve is the demarcation boundary between the stable and unstable region related to DWO. The characteristics of the Hopf curve is analyzed by computing FLC. Further analysis of this DWO
3.3. Numerical simulation In this section, the time domain simulations are described. Analysis for a sample case of Npch = 12 (in Fig. 4) is performed below. 3.3.1. Region X1 For lower value of us+ (marked as region X1 in Fig. 4), real part of largest eigenvalue is positive (Fig. 5). It signifies unstable system behaviour. If the system is initialized from an operating point in this region, a small perturbation results in the divergence of the trajectory. Time evolution of inlet velocity is shown in Fig. 10.
Fig. 7. Stability Map for Density Wave Oscillations and Flow Maldistribution (Npch = 10 ); Blue Curve - Density Wave Oscillations (DWO) boundary i.e, Hopf Curve, Region-A is unstable region; Black Curve - Flow Maldistribution (FMD) boundary i.e, Pitchfork Curve; Region-B is Symmetrical stable region; Region-C is the region of Flow Maldistribution i.e, Asymmetrical stable region. 6
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
Fig. 10. Numerical Simulation in region X1 of Fig. 4, perturbation leads to divergence of trajectory.
Fig. 8. Stability Map for Density Wave Oscillations (DWO) and Flow Maldistribution (FMD) for various Npch values; Each color represents the pair of DWO and FMD boundaries.
Fig. 11. Numerical Simulation in region X2 of Fig. 4, perturbation leads to convergence of trajectory.
Fig. 9. Stability Map for DWO and FMD with overlapping region(D) for Npch = 20 ; FMD boundary marked in blue color; DWO boundary marked in green color.
3.3.2. Region X2 As us+ is further increased in Fig. 4, Hopf point is observed. Beyond the Hopf point (marked as region X2), the equilibrium curve posses the real part of the largest eigenvalue to be negative (Fig. 5). It signifies stable system behaviour. If the system is initialized from an operating point in this region, a small perturbation results in the convergence of trajectory towards the fixed point. Time evolution of inlet velocity is shown in Fig. 11. 3.3.3. Region X3 As us+ is further increased in Fig. 4, BP(or pitchfork bifurcation) is observed. This marks the onset of FMD. From BP, two sets of equilibrium curve emanates. One corresponds to stable solution (marked as region Y1/Y2 in Fig. 4, asymmetric part, dotted lines) and other unstable solution (marked as region X3 in Fig. 4, symmetric part, solid line). The branches which emanates from the BP are the complementary solutions of the twin channel system. For symmetric solution, sign of the largest eigenvalue is positive (Fig. 5). If the system is initialized from an operating point in the region X3, trajectory diverges on account
Fig. 12. Numerical Simulation in region X3 of Fig. 4, perturbation causes trajectory to settle at an asymmetric stable solution. 7
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
Fig. 15. Variation of FLC vs Nsub (FLC is positive along the Hopf curve).
Fig. 13. Numerical Simulation in region Y1/Y2 of Fig. 4, perturbation leads to convergence of trajectories and are asymmetrically stable.
of trajectory. In region A, point 2 (Fig. 7) lies in the unstable region. Here, perturbation leads to divergence of the trajectory. Such is the nature of sub-critical Hopf region (section 2.3.1). This kind of identification of the sub-critical Hopf region around Hopf (DWO) boundary has been extensively studied in the past from the context of DWO. Therefore, the characterization of the Hopf stability boundary is not the interest of present work.
of small perturbation and settles on the solution of stable branches. Time evolution of inlet velocity is shown in Fig. 12. 3.3.4. Region Y1/Y2 For an asymmetric solution marked as Y1/Y2 in Fig. 4, sign of the largest eigenvalue is negative (Fig. 6). It signifies stable solution. If the system is initialized from an operating point in the region Y1/Y2, trajectory converges on account of small perturbation. Time evolution of inlet velocity is shown in Fig. 13.
3.3.6. Region D As us+ is increased we first observe FMD boundary marked in blue color in Fig. 9. If a system is initialized from an operating point marked in region A, we expect diverging DWO. Post occurrence of FMD boundary, we encounter region D. As we further increase us+ , DWO boundary marked in green color is observed. Such region D is observed for Nsub above 5. For lower Nsub values, symmetric stable solution region is observed. The same can be found in Fig. 7. Now, the region D as observed in Fig. 9, between FMD and DWO boundary has very interesting dynamics. Both the characteristics of FMD and DWO are observed here. If the system is initialized from an operating point in region D, system symmetry breaks and flowrate in the two channels become different. Flow maldistribution is thus observed first. After FMD boundary, it is the region of asymmetrical stable solutions. It is to be
3.3.5. DWO boundary characteristics The characteristics of the DWO curve (i.e, Hopf Curve) as shown in Fig. 7, is further analyzed by the plotting the FLC along the curve. It is shown in Fig. 14 and Fig. 15. The positive sign of FLC signifies the subcritical nature of the stability boundary. This implies the presence of an unstable limit cycle in the linearly stable region which emanates from the stability boundary. Here, small perturbation leads to convergence of DWO whereas large perturbation leads to divergence of DWO. In region B, point 1 (Fig. 7) lies in the stable region. Here, small perturbation leads to convergence of the trajectory. Also, when the system is initialized from the same point 1, large perturbation leads to divergence
Fig. 16. Numerical simulation from an operating point in region D of Fig. 9, presence of both the phenomenon i.e., FMD and DWO.
Fig. 14. Variation of FLC vs us+ (FLC is positive along the Hopf curve). 8
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
further noted from Fig. 7, that on the left of DWO curve, region A is unstable DWO region. So, post FMD in region D (Fig. 9), DWO in the system prevails. Thus dynamics of both flow maldistribution and density wave oscillations are observed in this region. It is to be noted that operating point of the system should be far away from this region. Here, symmetry break in the system on account of FMD is observed at the start of the transient. Furthermore, DWO is observed in the flow maldistributed system. Time evolution of inlet velocity is shown in Fig. 16.
pitchfork curve) demarcates symmetrical stable region and FMD region. Each point on the BP curve, there emanates multiple equilibrium curves. Stability of these curves are analyzed by the plotting the eigenvalues along the equilibrium curves. The region between the DWO curve (Hopf curve) and the FMD curve (pitchfork curve) represents the stable operating region and should be the region of candidature for operating point selection. For certain value of parameters, a fourth region is also observed. It arises because of the overlap between FMD boundary and DWO boundary. Interaction of DWO and FMD is observed here. Stability analyses using dedicated models which captures a particular instability only demarcates the parametric space into stable and unstable region for that instability. A combined analysis using a single model provides a demarcation for dominance of multiple instabilities in a single parametric space. Proper use of bifurcation analysis using a single model is shown here to analyze two different instabilities. Such unified approach provides a better perspective of system design and system instability to identify the region of stable operating points. The study also recommends a future scope to analyze multiple type of instabilities in a multi-channel flow boiling system.
4. Summary and conclusion This work analyzes the density wave oscillations (DWO) and flow maldistribution (FMD) in a combined fashion through the development of a stability map. Bifurcation analysis is used to identify various instabilities and their boundaries. For the first time, FMD is identified by pitchfork bifurcation in this work. Three regions in the stability map spanned by non-dimensional numbers Nsub and us+ are obtained by DWO and FMD boundaries. They are demarcated as DWO region, a region of the symmetrical stable solutions and FMD region. The DWO (Hopf) boundary demarcates the density wave oscillations and symmetrical stable region. Branch Point (BP) curve (i.e, FMD boundary or Nomenclature Acronyms
BP Branch Point DWO Density Wave Oscillations FLC First Lyapunov Coefficient FMD Flow Maldistribution GH Generalized Hopf Greek Symbols δ λ
Dirac delta function Boiling boundary (m) Single phase friction number
( ) Two phase friction number( )
1
f1 LH 2DH
f2 LH
2
2DH
Sign of the real part of largest eigenvalue ρ Density (kg / m3) Subscripts 0 Steady state e Channel exit f Liquid g Vapour H Heated channel i Channel inlet nth node in single phase region n Total Tot Other Symbols
P A Ax+ D f g h k L M Ns
s
Npch
Pressure drop (Pa) Area (m2 ) Non-dimensional cross sectional area of the jth heated channel Diameter (m) Friction factor Gravity acceleration (m2 /s ) Enthalpy (kJ / kg ) Resistance coefficient Length (m) Mass (kg ) No. of nodes in single phase region Phase change number
(
Q vfg W hfg vf
) 9
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
P Q q+ q t u us+ W W+ z Nsub x
Pressure (Pa) Heating power (W ) q Non-dimensional heat flux q0 2 Heat flux (W / m ) Time (s) Velocity (m /s ) Non-dimensional velocity, (ui,1 + ui,2) Mass flowrate Non-dimensional mass flowrate Axial coordinate (m) vfg hi ) ) Sub-cooling number ( h v (hf fg f Vapour quality
Appendix A. Derivation of mathematical model The detailed mathematical model is given in the work of Lin et al. [27]; Lee et al. [8] and Clausse & Lahey Jr [28]. However, a brief description is presented here. Mass Conservation Equation,
t
+
z
( u) = 0
(A.1)
Energy Conservation Equation,
t
( h) +
z
( hu ) =
q pH
(A.2)
AH
Momentum Conservation Equation, t
( u) +
z
i=1
f 2De
P z
( u2) =
N k (z i
z i)
u2 2
u2 g
(A.3)
here, ρ is density (kg / m3) ; t, time (sec ) ; u, velocity (m /s ); z, axial coordinate (m) ; f, friction factor; De , equivalent Diameter (m) ; k, loss Coefficient; δ, dirac delta function; g, gravity constant (m /s 2) ; P, Pressure (Pa) ; h, enthalpy (kJ / kg ) ; q", heat flux (W / m2) ; pH , heated perimeter (m) ; AH , crosssectional area (m2 ) . The density ρ can be written as,
=
(
= vf +
vfg
f,
(h
hfg
for h hf )
)
hf (Single 1
Phase region)
, for h > hf (Two
Phase region)
(A.4)
vf , specific volume of saturated liquid hfg , here, subscript f is fluid; vfg is difference in specific volume of saturated liquid and vapour latent heat of evaporation ( J / kg ). As shown in Fig. 1, the channel is divided into two regions, single and two-phase region. Entire channel is considered to be heated. Here, single phase region is divided into Ns number of nodes having variable length. Linear enthalpy profile is considered in single phase region between the two neighbouring nodes. So, the fluid enthalpy at the node can be defined as, (m3/ kg );
hn = h i +
n (h f
(m3/ kg );
hi ) (A.5)
Ns
here, subscript n is node in single-phase region; subscript i, inlet of channel; Ns , no. of nodes in single-phase region. A set of differential equations can be obtained by integration of energy balance equation (eqn. (A.2)) between two successive nodes (in singlephase region),
nth
Ln Ln 1 dLn dt
t
( h) dz +
= 2u i
Ln Ln 1
z
( hu) dz =
q p 2Ns A (hH h ) (Ln H f f i
Ln q pH Ln 1 AH
L n 1)
dz
dLn 1 dt
(A.6)
here, L is the length (m) ; L Ns , boiling boundary of the channel (m) . So, eqn (A.6) refers to the dynamic boiling boundary of the channel. The exit velocity ue can be calculated from eqn (A.1), eqn. (A.2) and density ρ (eqn. (A.4)),
q pH vfg u = = z AH hfg
(A.7)
Integration of the above equation from the boiling boundary to exit of the channel gives,
ue = ui +
(LH
L Ns )
(A.8)
here, subscript H is heated section of the channel; e, channel exit. Integration of the mass equation (eqn. (A.1)) from the inlet to the exit of the channel yields the following equation for the channel mass,
10
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
dMH = ( f ui dt
MH =
f
e ue ) AH
(A.9)
AH L Ns + M2
(A.10)
where, M is the coolant mass of the channel and subscript 2 is two phase. Also,
M2 = AH (LH
ln( f / e )
L Ns )
f
f/ e
1
(A.11)
Now, combining eqn. (A.9) and eqn. (A.10) we get ln
d e dt
=
f
ui
e ue
f / e
1
f
f / e
dL Ns dt
1
×
2 e 2 f (LH
L Ns )
f
e/ f
1
ln
e/ f
(A.12)
The dynamics of inlet velocity can be obtained by balancing the dynamic pressure drop in the channel and the external pressure head corresponding to the steady state. Pressure drop ( P ) in the channel is calculated by integrating the momentum balance equation (eqn. (A.3)), (A.13)
Pexternal = PI + Pg + Pf + Pa where, Inertial Pressure Drop, LH
PI = =
d dt
( u) dz t
0
(ue
MH ui AH
=
LH
d dt
ui) LH f
+
u dz
0
MH / AH
f / e
1
(A.14)
Gravitation Pressure Drop, LH
Pg =
g dz = gMH /AH
0
(A.15)
Frictional Pressure Drop, LH
Pf = =
f DH
0
L Ns f u2 dz 2DH
0
=
2 f ui L Ns
f1
2DH
2ui (ue ui) [ f f / e 1
(
f
3
e
)
+
+ +
i=1
N k (z i
LH f u2 dz L Ns 2DH
+
f2
L Ns
LH
2DH
1/ e
ki u2 2 i i
1/ f
MH / AH ] +
LH
ui2 ln
L Ns 2
ki
2 f ui
2
(LH
+
+ ke
L Ns )ln f / e
ke u2 2 e e
+
( )+ f
e
ui )2
(u e f / e
LH
u2 dz 2
zi)
2 f×
1
f / e
1
2 e ue
(A.16)
2
Acceleration Pressure Drop,
Pa = =
0 2 e ue
LH
u2 dz z 2 f ui
(A.17)
All the above equations are non-dimensionalized using the following definitions, A+ = A/ AH , non-dimensional cross-section area of the channel P + = P / f us2 , non-dimensional pressure drop t + = t /(LH /us ), non-dimensional time + = / , non-dimensional density f + u = u/ us , non-dimensional inlet velocity us = ui , characteristic velocity L+ = L / LH , non-dimensional channel length + = /L H , non-dimensional channel length z+ = z / LH , non-dimensional length
11
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al. h
h+ =
hf
Q0 /
f Aus
, non-dimensional enthalpy
M+ = M /( f LH AH ) , non-dimensional channel mass W + = W /( f Aus ) , non-dimensional mass flowrate Eu = Pexternal /( f us2 ) , Euler Number Fr = us2/(gLH ) , Froude number Subscript s refers to standard value(characteristic value). Hence, the conservation equations (refer eqn. (A.1), refer eqn. (A.2), refer eqn. (A.3)) can be written in non-dimensional form as, +
+
t+
t+ t+
z+
( +u+) = 0
( +h+) +
z+
( +u+) +
z+
i=1
(A.18)
( +h+u+) = 1
( +u+2) =
N k (z+ i
(A.19)
Eu
zi+)
P+ z
+u+2
+u+2
+
2
(A.20)
Fr
where, + = 1, for single phase region = (1 + Npch h+) 1, for two phase region
and = 1 , single-phase region; = 2 , two-phase region Fluid enthalpy at the node in the single phase region (refer eqn. (A.5)) can be written as,
Ns
hn+ =
n Ns
hi+
(A.21)
The differential equation of boiling boundary (refer eqn. (A.6)) can be written as,
dLn+ dt +
= 2ui+
2Ns
Npch Nsub
where, n = 1…Ns and
(Ln+
L N+s
Ln+ 1) +)
(or
dLn+ 1 dt +
(A.22)
is the dynamic boiling boundary. Npch =
Q vfg W hfg vf
, phase change number; Nsub =
vfg (hfg vf )
Integrating the mass balance equation (eqn. (A.18)) from inlet to exit of the channel yields (refer eqn. (A.9)),
dMH+ = ui+ dt +
+ + e ue
(hf
hi ) , sub-cooling number.
(A.23)
where (refer eqn. (A.10), eqn. (A.11)),
MH+ =
+
+ ln( e+) +) e + 1 e
+ (1
(A.24)
and (refer eqn. (A.8))
ue+ = ui+ + Npch (1
(A.25)
+)
The equation for dynamic fluid density (refer eqn. (A.12)) at the channel exit can be obtained from eqn. (A.23), eqn. (A.24) and eqn. (A.25), d e+ dt +
=
1+
+ + e ln( e ) + 1 e
(1 (1
+)[1
d + dt +
+
+ + e ue
ui+ ×
+ 2 e ) + + e + ln( e )]
(A.26)
Dynamics of inlet velocity can be obtained using (refer eqn. (A.13)), + Eu = PI+ + Pg+ + P + f + Pa
(A.27)
Using non-dimensionless form of inertial pressure drop (refer eqn. (A.14)), gravitation pressure drop (refer eqn. (A.15)), frictional pressure drop (refer eqn. (A.16)) and acceleration pressure drop (refer eqn. (A.17)), rearrangement yields,
dui+ (Eu PH+ ) = + + dt MH
(A.28)
where,
12
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al. +) (1 MH (1 / e+ 1)
PH+ =
+)
Npch (1 (1 / e+)
+
+ MH
Fr
(
+ 1+
+
1 ke 2
)
(
d +
d e+
Npch
dt +
+ 2 e ue
(
+
2ui+Npch (1
Npch (1 (1 / e+)
+) 2
1
1
+ e
+) MH
+
1 2
+) Mch
2 e)
)
1)
3
(1
ln(1/
ui+
dt +
1 ui2 +
+)(1
(1 / e+
+ dMH
+)(1
(1
ki 2
+ 1 u2 (1 / e+) 1 i
2
)
Npch dt + +
+u ++ i
1
+ e )+
+
+ MH+
+
(A.29)
Thus, the model can be represented as,
x = [Ln+
+ e
ui+] , n = 1…Ns
(A.30)
Appendix B. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.ijthermalsci.2019.106026.
References [15]
[1] Rizwan-Uddin, J. Dorning, Some nonlinear dynamics of a heated channel, Nucl. Eng. Des. 93 (1986) 1–14 http://www.sciencedirect.com/science/article/pii/ 0029549386901901 https://doi.org/10.1016/0029-5493(86)90190-1. [2] J. Boure, A. Bergles, L. Tong, Review of two-phase flow instability, Nucl. Eng. Des. 25 (1973) 165–192 http://www.sciencedirect.com/science/article/pii/ 0029549373900435 https://doi.org/10.1016/0029-5493(73)90043-5. [3] L.C. Ruspini, C.P. Marcel, A. Clausse, Two-phase flow instabilities: a review, Int. J. Heat Mass Transf. 71 (2014) 521–548 http://www.sciencedirect.com/science/ article/pii/S0017931013010922 https://doi.org/10.1016/j.ijheatmasstransfer. 2013.12.047. [4] S. Kakac, B. Bon, A review of two-phase flow dynamic instabilities in tube boiling systems, Int. J. Heat Mass Transf. 51 (2008) 399–433 http://www.sciencedirect. com/science/article/pii/S0017931007005959 https://doi.org/10.1016/j. ijheatmasstransfer.2007.09.026. [5] J. Lee, C. Pan, Dynamics of multiple parallel boiling channel systems with forced flows, Nucl. Eng. Des. 192 (1999) 31–44 http://www.sciencedirect.com/science/ article/pii/S0029549399000850 https://doi.org/10.1016/S0029-5493(99) 00085-0. [6] X. Lu, Y. Wu, L. Zhou, W. Tian, G. Su, S. Qiu, H. Zhang, Theoretical investigations on two-phase flow instability in parallel channels under axial non-uniform heating, Ann. Nucl. Energy 63 (2014) 75–82 http://www.sciencedirect.com/science/ article/pii/S030645491300385X https://doi.org/10.1016/j.anucene.2013.07.030. [7] A.M. Mishra, S. Singh, Non-linear stability analysis of uniformly heated parallel channels for different inclinations, Appl. Therm. Eng. 98 (2016) 1189–1200 http:// www.sciencedirect.com/science/article/pii/S1359431115013666 https://doi.org/ 10.1016/j.applthermaleng.2015.11.118. [8] J.-D. Lee, Y.-G. Lin, S.-W. Chen, C. Pan, The influence of void-reactivity feedback on the bifurcation phenomena and nonlinear characteristics of a single nuclear-coupled boiling channel, Ann. Nucl. Energy 94 (2016) 814–825 http://www. sciencedirect.com/science/article/pii/S0306454916302110 https://doi.org/10. 1016/j.anucene.2016.04.044. [9] T. Van Oevelen, J.A. Weibel, S.V. Garimella, Predicting two-phase flow distribution and stability in systems with many parallel heated channels, Int. J. Heat Mass Transf. 107 (2017) 557–571 http://www.sciencedirect.com/science/article/pii/ S0017931016324395 https://doi.org/10.1016/j.ijheatmasstransfer.2016.11.050. [10] K. Akagawa, M. Kono, T. Sakaguchi, M. Nishimura, Study on distribution of flow rates and flow stabilities in parallel long evaporators, Bull. JSME 14 (1971) 837–848, https://doi.org/10.1299/jsme1958.14.837. [11] U. Minzer, D. Barnea, Y. Taitel, Evaporation in parallel pipes splitting characteristics, Int. J. Multiph. Flow 30 (2004) 763–777 http://www.sciencedirect.com/ science/article/pii/S0301932204000527 https://doi.org/10.1016/j. ijmultiphaseflow.2004.04.006. [12] U. Minzer, D. Barnea, Y. Taitel, Flow rate distribution in evaporating parallel pipes―modeling and experimental, Chem. Eng. Sci. 61 (2006) 7249–7259 http:// www.sciencedirect.com/science/article/pii/S0009250906005173 https://doi.org/ 10.1016/j.ces.2006.08.026. [13] M. Baikin, Y. Taitel, D. Barnea, Flow rate distribution in parallel heated pipes, Int. J. Heat Mass Transf. 54 (2011) 4448–4457 http://www.sciencedirect.com/science/ article/pii/S0017931011002535 https://doi.org/10.1016/j.ijheatmasstransfer. 2011.04.034. [14] Y. Taitel, D. Barnea, Transient solution for flow of evaporating fluid in parallel pipes using analysis based on flow patterns, Int. J. Multiph. Flow 37 (2011) 469–474
[16]
[17]
[18]
[19]
[20]
[21] [22]
[23]
[24]
[25]
[26]
[27]
[28]
13
http://www.sciencedirect.com/science/article/pii/S030193221100022X https:// doi.org/10.1016/j.ijmultiphaseflow.2011.01.002. V. Pandey, S. Singh, Characterization of stability limits of ledinegg instability and density wave oscillations for two-phase flow in natural circulation loops, Chem. Eng. Sci. 168 (2017) 204–224 http://www.sciencedirect.com/science/article/pii/ S0009250917302865 https://doi.org/10.1016/j.ces.2017.04.041. V. Pandey, S. Singh, Bifurcation analysis of density wave oscillations in natural circulation loop, Int. J. Therm. Sci. 120 (2017) 446–458 http://www.sciencedirect. com/science/article/pii/S1290072916312248 https://doi.org/10.1016/j. ijthermalsci.2017.06.029. I. Simanovskii, A. Viviani, F. Dubois, Nonlinear dynamics of two-layer systems with a temperature-dependent heat release, Int. J. Therm. Sci. 119 (2017) 51–67 http:// www.sciencedirect.com/science/article/pii/S1290072916317975 https://doi.org/ 10.1016/j.ijthermalsci.2017.05.022. K. Luo, H.-L. Yi, H.-P. Tan, Eccentricity effect on bifurcation and dual solutions in transient natural convection in a horizontal annulus, Int. J. Therm. Sci. 89 (2015) 283–293 http://www.sciencedirect.com/science/article/pii/S1290072914003329 https://doi.org/10.1016/j.ijthermalsci.2014.11.020. S. Paruya, S. Maiti, A. Karmakar, P. Gupta, J. Sarkar, Lumped parameterization of boiling channel―bifurcations during density wave oscillations, Chem. Eng. Sci. 74 (2012) 310–326 http://www.sciencedirect.com/science/article/pii/ S0009250912001352 https://doi.org/10.1016/j.ces.2012.02.039. D.N. Basu, S. Bhattacharyya, P. Das, Performance comparison of rectangular and toroidal natural circulation loops under steady and transient conditions, Int. J. Therm. Sci. 57 (2012) 142–151 http://www.sciencedirect.com/science/article/pii/ S1290072912000567 https://doi.org/10.1016/j.ijthermalsci.2012.02.011. A. Karmakar, S. Paruya, Nonlinear analysis of chaotic time series in a natural circulation boiling loop, Ind. Eng. Chem. Res. 51 (2012) 16467–16481 doi:10.1021/ ie300788t https://doi.org/10.1021/ie300788t. M. Krishnani, D.N. Basu, On the validity of boussinesq approximation in transient simulation of single-phase natural circulation loops, Int. J. Therm. Sci. 105 (2016) 224–232 http://www.sciencedirect.com/science/article/pii/S1290072915300065 https://doi.org/10.1016/j.ijthermalsci.2016.03.004. S. Paul, S. Singh, On nonlinear dynamics of density wave oscillations in a channel with non-uniform axial heating, Int. J. Therm. Sci. 116 (2017) 172–198 http:// www.sciencedirect.com/science/article/pii/S1290072916307761 https://doi.org/ 10.1016/j.ijthermalsci.2017.02.008. G. Xia, G. Su, M. Peng, Analysis of flow distribution instability in parallel thin rectangular multi-channel system, Nucl. Eng. Des. 305 (2016) 604–611 http:// www.sciencedirect.com/science/article/pii/S0029549316300449 https://doi.org/ 10.1016/j.nucengdes.2016.04.016. A. Dhooge, W. Govaerts, Y.A. Kuznetsov, MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs, ACM Trans. Math Software 29 (2003) 141–164 doi:10.1145/779359.779362 http://doi.acm.org/10.1145/779359. 779362. F. Alobaid, N. Mertens, R. Starkloff, T. Lanz, C. Heinze, B. Epple, Progress in dynamic simulation of thermal power plants, Prog. Energy Combust. Sci. 59 (2017) 79–162 http://www.sciencedirect.com/science/article/pii/S0360128516300375 https://doi.org/10.1016/j.pecs.2016.11.001. Y. Lin, J. Lee, C. Pan, Nonlinear dynamics of a nuclear-coupled boiling channel with forced flows, Nucl. Eng. Des. 179 (1998) 31–49 http://www.sciencedirect.com/ science/article/pii/S0029549397002422 https://doi.org/10.1016/S00295493(97)00242-2. A. Clausse, R. Lahey Jr., An investigation of periodic and strange attractors in boiling flows using chaos theory, Heat Transfer 1990. Proceedings of the Ninth
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al. International Heat Transfer Conference, 1990. [29] P. Saha, M. Ishii, N. Zuber, An experimental investigation of the thermally induced flow oscillations in two-phase systems, ASME J. Heat Transf. 98 (1976) 616–622 http://heattransfer.asmedigitalcollection.asme.org/article.aspx?articleid= 1436574 https://doi.org/10.1115/1.3450609. [30] Rizwan-Uddin, J. Dorning, Some nonlinear dynamics of a heated channel, Nucl. Eng. Des. 93 (1986) 1–14 http://www.sciencedirect.com/science/article/pii/ 0029549386901901 https://doi.org/10.1016/0029-5493(86)90190-1. [31] S. Paul, S. Singh, A density variant drift flux model for density wave oscillations, Int. J. Heat Mass Transf. 69 (2014) 151–163 http://www.sciencedirect.com/ science/article/pii/S0017931013008727 https://doi.org/10.1016/j.
ijheatmasstransfer.2013.10.012. [32] M. Ishii, Thermally Induced Flow Instabilities in Two-phase Mixtures in Thermal Equilibrium, Georgia Institute of Technology, 1971, https://books.google.co.in/ books?id=eWurtgAACAAJ. [33] P. Saha, Thermally Induced Two-phase Flow Instabilities, Including the Effect of Thermal Non-equilibrium between the Phases, Georgia Institute of Technology, 1974, https://books.google.co.in/books?id=oc2RbwAACAAJ. [34] Y. Kuznetsov, Elements of Applied Bifurcation Theory. Applied Mathematical Sciences, Springer, New York, 2004https://books.google.co.in/books?id= B6Z2TQSRJuMC.
14
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Interaction of density wave oscillations and flow maldistribution for twophase flow boiling parallel channels
T
Deepraj Paula,b, Suneet Singha,∗, Surendra Mishrab a b
Department of Energy Science and Engineering, Indian Institute of Technology-Bombay, Mumbai, India Nuclear Power Corp. of India Ltd, Mumbai, India
ARTICLE INFO
ABSTRACT
Keywords: Thermal hydraulics Stability Flow maldistribution Density wave oscillations Bifurcation
Density wave oscillations (DWO) and flow maldistribution (FMD) are two important instabilities among others in two-phase flow boiling systems. As both the phenomena can lead to undesirable consequences such as channel burnout, the design intent is to avoid the same. DWO is caused by the interaction between single-phase and twophase region as perturbation propagates slowly in the later region. Whereas in case of FMD, it leads to one channel receiving higher flowrate and other channel receiving lower flowrate for a twin parallel channels system under forced flow condition. Thus, channels subjected to similar thermal hydraulic and geometric condition receive non-identical flowrate. Bifurcation analysis is used to identify DWO and FMD instabilities in this work. DWO is identified by Hopf bifurcation, thus Hopf curve which constitutes Hopf bifurcations represents the DWO boundary. FMD has been observed earlier through numerical simulation but the region which depicts the same in a parametric space has not been demarcated in previous studies. In this work, FMD is identified by a pitchfork bifurcation. This leads to straight-forward identification of FMD boundary as a pitchfork curve. A combined analysis of multiple instabilities using a single stability map can demarcate regions with the dominance of particular instability. It can thus provide a better perspective for system design and system instabilities in terms of suitable operating region selection. A stability map which can capture and demarcate both the phenomena simultaneously i.e, DWO and FMD in a parametric space is the interest of present work. The boundary of DWO and FMD splits the parametric space into three distinct regions. These are DWO region, a region of symmetrical stable solutions and FMD region which contains non-identical (asymmetric) stable solutions. This breakdown of symmetry of the system is observed by non-identical flowrate in the channels. In some cases, boundary depicting DWO and FMD intersects. Hence, a fourth region is also observed which possess the characteristics of both the instabilities. The stability map shows the transition between these regions.
1. Introduction Flow instabilities can occur within a two-phase flow boiling heated channels. The best approach to deal with such phenomena is to design the system in such a way which eventually avoids the same. A thermal hydraulics model is based on a set of conservation equations and boundary conditions. In literature, several models have been proposed for the analysis of two-phase flow behaviour. The common approach to model the two-phase flow is by using Homogeneous equilibrium model (HEM). Other models such as drift-flux model by Rizwan-Uddin & Dorning [1]; separated flow model, slip flow model, stratified flow model etc. are also used. In the homogeneous flow model, it is assumed that velocity, temperature and pressure between the liquid and vapour phases are equal. The drift-flux model improves the homogeneous flow
∗
model, by considering slip velocities between the two phases. Whereas, a two-fluid approach considers different sets of conservation equations of each phases and their interactions. This signifies more details towards the mathematical modelling. Moreover from the system analysis perspective, optimization between model order and level of detailing needs to be analyzed first for reducing computational expense. The most popular classification, introduced in Boure et al. [2] divides two-phase flow instabilities in static and dynamic instabilities. For static instabilities, the threshold of unstable behaviour can be predicted from the steady-state conservation laws. On the other hand, to describe the behaviour of dynamic instabilities it is necessary to take into account different dynamic effects, such as the propagation time, inertia, compressibility, etc. In addition, the term compound instability is normally used when several of the basic mechanisms interact with each
Corresponding author. E-mail address: [email protected] (S. Singh).
https://doi.org/10.1016/j.ijthermalsci.2019.106026 Received 30 July 2018; Received in revised form 29 June 2019; Accepted 12 July 2019 1290-0729/ © 2019 Elsevier Masson SAS. All rights reserved.
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
other. A considerable amount of studies has been carried out in this field by Ruspini et al. [3]; Kakac & Bon [4] etc. A lot of 1D analyses has been carried by researchers for parallel channels viz. Lee & Pan [5]; Lu et al. [6]; Mishra & Singh [7] etc, using various numerical schemes. DWO for thermal hydraulics coupled with neutronics for a single channel was analyzed by Lee et al. [8]. The influence of void-reactivity feedback on the bifurcation phenomena and nonlinear characteristics of a single nuclear-coupled boiling channel in an advanced boiling water reactor was studied. As noted, DWO is the most widely studied dynamic instability which occurs due to the delayed feedback between single and two-phase regions. Since the perturbation propagates slowly along the two-phase region, a significant delay marks the onset of perturbations in two-phase region. Hence, the two-phase pressure drop and the single phase pressure drop oscillates out-of-phase. It is thus very important to identify the regions in the parametric space where such instabilities exists. Van Oevelen et al. [9] discussed that FMD can occur because of nonmonotonic nature of pressure drop vs flowrate curve in a heated channel, asymmetrical inlet header design, differences in channel geometries etc. His work proposes a new numerical approach to analyze the stability based on a generalized eigenvalue problem. This methodology is specifically designed for analyzing systems with hundreds of identical parallel channels. Due to the non-monotonicity of the pressure drop curve, there exist multiple solutions. For a single channel system, driven by a constant pressure boundary there are three solutions. For two channel system, as both the channels have the same pressure drop so the number of solutions are nine. Thus the number of solutions increases as the number of channels increase. Flow rate distribution among parallel evaporator tubes was studied by Akagawa et al. [10]. They experimentally obtained channel load curves for individual tubes experiencing flow boiling. Flow excursion was also observed, as the system may not be stable at all the intersection points of the internal and external pressure drop vs flowrate curve. Minzer et al. [11] experimentally observed the presence of maldistribution in heated channels for a forced flow system. Here, asymmetric but stable solutions exist, where one channel receives higher flowrate but another channel is subjected to lower flowrate. Minzer et al. [12] developed a model where each channel is represented by a single ordinary differential equation. Coupling between multiple channels of the thermal-hydraulic system was done using boundary condition. For 2 number of parallel channels, FMD was observed even though channels were subjected to the same heat flux. Baikin et al. [13] analyzed FMD in a system which comprises of 4 number of parallel channels. This lead to even more number of possible solutions. A combination of heated and unheated pipes in the system was analyzed in the work for FMD. Even for the case of all pipes being heated equally, there were ranges of inlet flowrate where severe FMD could occur. The model was further developed in the work of Taitel & Barnea [14] where the pipe was subdivided into multiple cells. The calculation of the pressure drop was done in each cell. The transient simulations were performed and was observed that initiation of the system at unstable fixed points leads to the divergence of trajectories which settles at an uneven flowrate being distributed among the channels. Thus, there is an asymmetric distribution of flowrate among channels. However, the model used in the above studies assumes the mass flowrate to be uniform along the channel (quasiequilibrium assumption). A perturbation propagates slowly in twophase region as compared to single-phase region, which means there is phase-difference between them. DWO occurs because of delayed feedbacks between the inertia of flow and compressibility of the two phase mixture. Thus with quasi-equilibrium assumption in the above approach, one cannot capture DWO dynamics. Pandey & Singh [15] studied the excursive instability and density wave oscillations (DWO) for two-phase flow in a natural circulation loop. The stability maps obtained through their work by using bifurcation analysis has three distinct regions, viz. the first region has only DWO, second region has Ledinegg as well as DWO, while the third
region has only Ledinegg instability. Using bifurcation analysis to identify various flow instabilities, Pandey & Singh [16] explored a new type of stability boundary for high sub-cooling number in a two-phase natural circulation loop. The region near to this boundary shows sustained oscillations for flow velocity; confirming that this is a dynamic instability like DWO. This stability boundary is near to Ledinegg instability region but does not show the characteristics of Ledinegg stability. Simanovskii et al. [17] studied the influence of a temperaturedependent inter-facial heat release/consumption on nonlinear oscillatory convective regimes in a two-layer system with rigid heat-insulated lateral walls under the action of an inclined temperature gradient. A wide range of parameter can lead to the change of the sequence of bifurcations and the development of new nonlinear regimes in the system. Luo et al. [18] studied the bifurcation phenomena and the existence of dual solutions in natural convection in eccentric annulus. The onset and evolution of pitchfork bifurcation were tracked by the lattice Boltzmann method (LBM) that was inherently transient and has a good flexibility in dealing with irregular geometries. Paruya et al. [19] analyzed DWO in several boiling channels with varying lengths (Froude number) by adopting a moving node scheme (MNS) and fixed node scheme (FNS). It was observed that MNS has a higher-order convergence and high computational efficiency compared to finite difference scheme and FNS. It was also demonstrated that the solutions bifurcate from a stable steady state to periodic solutions to chaotic attractor ones when the boiling channel was subjected to the periodic perturbation of pressure drop at different power level. Basu et al. [20] focused upon the development of an unified model for the analysis of single-phase natural circulation loops with constant heat flux heating and convective cooling. Characterizing parameters were defined such that they are applicable irrespective of any geometrical configuration. Karmakar & Paruya [21] reported the chaotic flow oscillations in a natural circulation boiling loop. It was observed that boiling regime changes from nucleate boiling to transition boiling as a result of decrease in inlet sub-cooling and increase in heat flux. The thermal oscillations were strongly coupled with pressure-drop oscillations. Nonlinear analysis of the time series of loop flow rate at various heater power and inlet sub-cooling were carried out using statistical analysis. Krishnani & Basu [22] focused on assessing the viability of Boussinesq approximation during nonlinear stability appraisal following transient simulation of single-phase natural circulation loops. Paul & Singh [23] studied the stability behaviour of various uniform and non-uniform axial heat flux. The study reveals that keeping the total heat rate same, the stability characteristics of the system significantly changes from uniform to non-uniform axial heat flux assumption. Xia et al. [24] studied the flow distribution and flow instability through numerical simulation at various parametric points in a stability map. Such an approach is a tedious one and is computationally expensive. Moreover, a clear demarcation of FMD in the parametric space was not shown. A combined analysis which captures DWO and FMD in a forced flow system has been lacking in the literature. Bifurcation analysis is used in this work for identification of various instabilities and their boundary in a parametric space. Unlike DWO which are identified by Hopf bifurcation, identification of FMD by the use of bifurcation analysis has not been done earlier. In this work FMD is identified by pitchfork bifurcation. Our focus also involves the identification of DWO and FMD boundary in the stability map. The two boundaries splits the parametric space in three regions. These are DWO region, a region of the symmetrical stable solutions and FMD region. For a different parametric combination, a fourth region can also be observed. A MATLAB based mathematical package named MATCONT developed by Dhooge et al. [25] is used in this work to detect and analyze DWO and FMD. MATCONT provides the means for continuing equilibria of systems of Ordinary Differential Equations (ODEs), and their bifurcations. Identification of DWO and FMD along with their boundaries are computed using the package. Numerical simulations are carried out to verify and understand the characteristics of the system 2
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
Fig. 2. Experimental data validation of stability boundary, Ns = 3,5 & 7 .
• Heat flux is uniform in the axial direction • Fixed inlet sub-cooling • Sub-cooled boiling is not considered here • Constant thermodynamic properties at system pressure as pressure
drop through the channel is much smaller than that of system pressure
Fig. 1. Single Channel System, Lee et al. [8].
The derivation of the mathematical model is provided in Appendix A. The thermal hydraulic model (A.30) is validated (Fig. 2) with the reevaluated experimental data of Freon-113 from Saha et al. [29] as given in Rizwan-Uddin & Dorning [30] & Paul & Singh [31] (Set I). Nsub and Npch (as described in the work of Ishii [32] and Saha [33]) are the two non-dimensional parameters used for validation. These are significant not only for stability analyses, but also for the description of steady state operational conditions. Nsub takes into the account of time lag effects in the liquid region due to sub-cooling of the fluid entering the heated duct. It scales the inlet sub-cooling and it is the dimensionless residence time in the single-phase region under the equilibrium assumption. Thus it is one of the important parameters for the stability analysis. Whereas, Npch scales the change of phase due to the heat transfer to the system. The single phase region is nodalized with 3, 5 & 7 nodes. It is observed that all the three different stability boundaries (corresponding to different nodalizations in the single-phase region) as calculated using the model almost overlap with each other. The calculated stability boundary matches reasonably with the boundary obtained from experimental data. For lower values of Nsub , there is some mismatch. Here, inlet enthalpy is closer to saturation enthalpy and may lead to superheated exit condition (following a transient). Such is beyond the scope of analysis of the present work. Our work focuses primarily on vapour quality, 0 < x < 1. A schematic representation of a parallel channel is shown in Fig. 3. For forced flow system, mathematical coupling between the channels is used as per the work of Lee & Pan [5];
around the stability boundaries. 2. Mathematical modelling and bifurcation 2.1. Model description A typical flow boiling channel is shown in Fig. 1. The liquid in subcooled condition enters the channel and heats up which enforces the change of phase of liquid. The point where the change of phase occurs i.e, the liquid starts to boil is called the boiling boundary. Beyond the boiling boundary, the temperature remains constant although an increase in enthalpy due to heat flux is accounted by positive vapour quality. The liquid thus leaves the channel as a two-phase mixture. 2.2. Model details Alobaid et al. [26] discussed various types of models with varied complexities viz. homogeneous equilibrium model (HEM), drift flux model and two-fluid separated flow model (which includes four-equation, five-equation, six-equation and seven-equation flow model) etc. HEM although simple to implement, it cannot capture the relative velocity between the two-phases. However, the uncertainty in modelling the interaction between phases at the inter-phase boundary poses considerable difficulties. These relations cannot be derived from fundamental physical laws and in most cases are based on empirical assumptions. Solving the resulting differential equations requires higher computational effort and entails parameters that may cause numerical instability, especially due to improper selection of interfacial terms. The difficulties associated with the two-fluid models can be significantly reduced by formulating the two-phase flow in terms of the mixture flow model. Increased computational expense, existence of uncertainty in modelling the interaction at the inter-phase boundary and mathematical complexities of multi equation two-fluid models drive us to the use Homogeneous equilibrium model (HEM). HEM is simple, computationally inexpensive and yet applicable to wide range of two-phase flow regimes. The details of the Homogenous Equilibrium Model used here can be found in the work of Lin et al. [27]; Lee et al. [8] and Clausse & Lahey Jr [28]. Assumptions for the flow boiling system are:
dui+, j dt +
= Aj
dui+,1 dt +
+ Bj , j = 2,3…M
here, j represents the channel number; sional inlet velocity of jth channel. and
dui+,1 dt +
+ dWTot
=
dt
1+
M j=2 M j=2
(1)
ui+, j
represents the non-dimen-
Ax+ s Bj
Ax+ s Aj
(2)
where, Ax+ s, j is the area ratio of jth channel to first channel. Thus, Ax+ s,1 = 1. When total mass flow rate is held constant for a forced flow system, 3
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
called co-dimension of a system. Codim-1 and Codim-2 are used to represent single and two parameters variation, respectively. For a dynamical system, the sign of real part of the largest eigenvalue characterizes the stability of the system. For a positive sign, a small perturbation leads to divergence of trajectory. For a negative sign, a small perturbation leads to convergence of trajectory towards a stable fixed point. Such analysis can provide us with a demarcation between the stable and unstable region in a parametric space. Nonlinear dynamics analysis or bifurcation analysis can identify the region in a parametric space which contains stable and unstable limit cycles. Hopf bifurcation, Generalized Hopf (GH) bifurcation and Pitchfork bifurcation are observed in this work and are discussed below. 2.3.1. Hopf bifurcation (density wave oscillations) Density Wave Oscillations are self-sustaining oscillations and are identified by a Hopf bifurcation. The Hopf bifurcation occurs when one pair of the complex eigenvalues become purely imaginary. The presence of complex eigenvalues signify oscillatory behaviour. The Hopf curve which constitutes Hopf bifurcation points represents the DWO boundary. The Hopf bifurcation is a codim-1 parametric bifurcation. Hopf curve is generated by codim-2 parametric bifurcation which constitutes the Hopf points. Evolution of limit cycles is the phenomena related to Hopf bifurcation. This type of Hopf bifurcation can be understood by the sign of the first Lyapunov coefficient (FLC). For a negative sign of FLC, it represents super-critical Hopf whereas a positive sign represents sub-critical Hopf. In sub-critical Hopf region, the presence of an unstable limit cycle in the linearly stable region is of prime concern. Here small perturbation to the system converges towards a stable fixed point and relatively large perturbation leads to diverging DWO. On the other hand, for super-critical Hopf bifurcation in a linearly unstable region, there is a presence of a stable limit cycle. It means, irrespective of the magnitude of perturbation, DWO converges towards a stable limit cycle.
Fig. 3. Parallel channel schematic.
+ WTot =
M
Wi+
(3)
j =1
where, m is the number of channels. So, + dWTot =0 + dt
dui+,1 dt +
=
(4)
1+
M A+ B j=2 x s j M A+ A j=2 x s j
2.3.2. Generalized Hopf (GH) bifurcation DWO boundary or Hopf boundary is obtained as codim-2 bifurcation. The curve splits the parametric space into two regions namely, stable DWO and unstable DWO regions. Generalized Hopf bifurcation or Bautin bifurcation is observed over the Hopf curve i.e, DWO boundary. Here, FLC is zero. GH point divides the DWO boundary in two regions namely sub-critical and super-critical Hopf regions characterized by positive and negative FLC respectively. It represents the transition from sub-critical Hopf to super-critical Hopf or vice-versa. As discussed earlier for sub-critical Hopf region, large perturbation leads to divergence of trajectories in linearly stable region of DWO boundary. So identification of GH point is very important.
(5)
where,
Aj = MH+,1/ MH+, j Bj =
(
PH+,1
(6)
PH+, j )
MH+, j
(7)
here, MH+, j is the non-dimensional channel mass of the jth channel (eqn (A.24)). PH+, j is the non-dimensional pressure drop of the jth channel (eqn (A.29)). Thus the model of 2 number of parallel channels can be represented as,
x = [Ln+,1
+ e,1
ui+,1 Ln+,2
+ e,2 ]
, n = 1…Ns
2.3.3. Pitchfork bifurcation (flow maldistribution) Flow excursion (in a single channel) occurs due to the presence of negative slope region in pressure drop vs flowrate curve. An operational point is stable when the slope of the internal characteristic curve is smaller than the slope of the external characteristic curve. In case of a multiple channels, connected by an inlet and exit plenum the pressure drop across them is same. So there exists multiple solutions as external pressure drop intersects the internal pressure drop multiple times. Naturally, due to the presence of negative slope region in pressure drop curve and system flowrate being constant (in a forced flow system) the operating point (if on the negative slope) shifts towards stable solutions. This leads to flow re-distribution in the system. Such multiple solutions can be identified by a pitchfork bifurcation. It is a codim-1 parametric bifurcation. In MATCONT, pitchfork bifurcation is identified using Branch Point (BP) bifurcation. BP marks the onset of FMD. As BP is observed in an equilibrium curve, it signifies the emancipation of other equilibrium curves from there. Thus BP is the influx/outflux of different equilibrium curves. Pitchfork bifurcation is of two types viz, sub-critical and super-critical. In the present work, sub-critical
(8)
The inlet velocity of the other channel is found using the boundary condition of total mass flowrate being constant. 2.3. Bifurcation analysis in context to flow instability The details related to bifurcation analysis and review of two-phase flow instabilities are available in Kuznetsov [34] and Ruspini et al. [3] respectively. However, a brief description of stability and bifurcation analysis for the identification flow instabilities and demarcation of flow stability boundaries are presented here. Bifurcation of a dynamical system refers to the quantitative changes in the dynamics of the system which is obtained by varying parameters. The parameter values at which they occur are called bifurcation points. The number of parameters which are varied for bifurcation to occur is 4
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
pitchfork bifurcation is observed and represents the presence of asymmetric solution but stable ones in the parallel channels. Symmetric solutions although present, are unstable. The stability of such equilibrium curves can be studied by eigenvalues. The locus of pitchfork bifurcation represents the FMD boundary. 3. Results and discussion This section is divided into three major parts. In the first part (section 3.1), the equilibrium points are generated using MATCONT by varying the non-dimensional inlet velocity us+ while remaining parameters are kept constant. These equilibrium points are represented by plotting boiling boundary ( +) as a function of us+ and various bifurcations are shown in these equilibrium curves. These are so-called codim-1 bifurcation curves. In the second part (section 3.2), the various type of stability boundaries are plotted in a parametric space spanned by Nsub and us+. Here, these two parameters are simultaneously varied while other parameters are kept constant. These are so-called codim-2 bifurcation curves. These curves divides the parameteric space into regions of different flow characteristics. In the third part (section 3.3), numerical simulations are carried out for different parametric values. The characteristics of the time evolution of state variables can be observed. Table 1 is used to generate the graphs for this study.
Fig. 4. Equilibrium curve corresponding to codim-1 parameter variation, boiling boundary ( +) vs us+ .
3.1. Equilibrium curves with variation in us+ Fig. 4 represents the variation of location of the boiling boundary ( +) with us+. As expected, the location of boiling boundary moves up with increase in us+ (represented by the solid lines). Moreover, it is observed that there is the presence of two bifurcation points namely, a BP (or pitchfork bifurcation) and a Hopf point. As us+ is increased, Hopf point is observed first and then a BP. From BP, additional equilibrium curves emanates (represented by dotted line) in addition to the curve represented by solid line. The channels have symmetric (solid lines) as well as asymmetric solutions (dotted lines). For symmetric part of the solution, both the channels have identical solution i.e, channels receive the same flowrate and have same boiling boundary. For asymmetric part of the solution, both the channels have different solutions viz. boiling boundaries, flowrate etc even though they are subjected to the same thermal hydraulic and geometric conditions. It is thus observed, that for a particular Npch value, there exist three different regions in an equilibrium curve. The first region is before the occurrence of Hopf point. The second region is in between the Hopf point and BP. The third region is beyond the occurrence of BP. For the green coloured equilibrium curve (Npch = 12 ), the regions are marked as X1 i.e, the first region, X2 i.e, the second region and X3 i.e, the third region. These three regions caters to the symmetrical solution of the system i.e, both the identical channels have same flowrate and boiling boundaries ( + ). The asymmetrical solution which occurs after BP, are marked as Y1/Y2. This marks the onset of FMD. (The marked points namely X1, X2, X3 and Y1/Y2 will be used later for numerical simulations, section 3.3.) To further analyze the solutions of these three regions marked as X1, X2 and X3 in Fig. 4, the real part of the largest eigenvalue along the equilibrium curve (of the solid line) are plotted in Fig. 5. For region X1, it is observed that as us+ is varied, for lower values before the observance of Hopf point, real part of the largest eigenvalue is positive.
Fig. 5. Variation of real part of the largest eigenvalue for symmetric part along the equilibrium solutions. (Truncated y-axis for proper representation of change of sign of largest eigenvalue).
This signifies unstable solutions. Numerical simulation is provided in section 3.3.1. As us+ is increased further, the sign of the eigenvalue changes from positive to negative as a Hopf point is encountered. This signifies stable solutions. This is marked as X2 region in Fig. 4. Numerical simulation is provided in section 3.3.2. Further increase in us+, BP is encountered on the equilibrium curve. The sign of the real part of the largest eigenvalue, Fig. 5 which caters the symmetric part of the solution changes to a positive value. This signifies the symmetric part of the solutions is unstable. This is marked as X3 region in Fig. 4. Numerical simulation is provided in section 3.3.3. Two set of equilibrium curves emanates (dotted lines) from the branch point (BP) in Fig. 4. These are marked as Y1/Y2. Each equilibrium curve corresponds to a channel which are complementary to each other. The real part of the largest eigenvalue vs non-dimensional boiling boundary (which itself is a function of velocity) are shown in Fig. 6. Negative sign of real part signifies stable solution. It can be observed from Fig. 4 that, as us+ is increased beyond BP, symmetry of the system breaks as these solutions have positive real parts of largest eigenvalue. Asymmetric solution exists and are stable, as observed by negative real part of eigenvalue. This confirms the presence of FMD. Such asymmetric solution leads to higher flowrate in one channel and lower flowrate in the other channel, which is undesirable. Numerical
Table 1 Parameters used for analysis. Pressure (MPa)
Ns
ui0 (m/s)
D (m)
L (m)
ki
ke
7
3
1.0
0.01
1
1
1
5
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
Fig. 6. Variation of real part of the largest eigenvalue for asymmetric part along the equilibrium solutions; (A)-Branch Y1 in Fig. 4; (B)-Branch Y2 in Fig. 4.
simulation is provided in section 3.3.4.
boundary is provided in section 3.3.5. Similarly, from the BP or pitchfork bifurcation, Nsub is varied simultaneously along with us+ as codim-2 parameters to generated a BP curve or the pitchfork curve i.e, FMD boundary. This is shown as black coloured curve in Fig. 7. The stability map thus generated is spanned by Nsub and us+. It is subdivided into three regions viz. DWO region (A), Symmetrical stable solution (B) and FMD region (C). Region A is unstable. Region B is the region of symmetrical stable solution. Region C is the region of FMD where asymmetric and symmetric solution co-exists. The symmetric solutions which exists in this region are not stable. Perturbation causes trajectories to diverge and settle to a stable solution. This causes one channel to receive higher flowrate and another channel receives lower flowrate. For different values of Npch , stability map for DWO and FMD are plotted in Fig. 8. It is observed that for higher values of Npch , the stability boundaries (i.e, the pair of DWO and FMD boundary) shifts right which demands higher value of us+ to attain the symmetrical stable solution region. An additional fourth region (D) can also be observed for a different parametric combination and is shown in Fig. 9. The presence of this additional region in the stability map is due to the overlap between the pitchfork curve (i.e, FMD) and the Hopf curve (i.e, DWO). Further discussion along with numerical simulation of this region is provided in section 3.3.6.
3.2. Stability boundary From the Hopf point which was found by varying us+ in Fig. 4, Nsub also is now varied simultaneously to generate the locus of Hopf boundary i.e, DWO boundary. This is codim-2 parametric stability boundary and is shown in Fig. 7. The blue curve represents the Hopf curve (or the DWO stability boundary), where each point on the curve is a Hopf point. This curve is the demarcation boundary between the stable and unstable region related to DWO. The characteristics of the Hopf curve is analyzed by computing FLC. Further analysis of this DWO
3.3. Numerical simulation In this section, the time domain simulations are described. Analysis for a sample case of Npch = 12 (in Fig. 4) is performed below. 3.3.1. Region X1 For lower value of us+ (marked as region X1 in Fig. 4), real part of largest eigenvalue is positive (Fig. 5). It signifies unstable system behaviour. If the system is initialized from an operating point in this region, a small perturbation results in the divergence of the trajectory. Time evolution of inlet velocity is shown in Fig. 10.
Fig. 7. Stability Map for Density Wave Oscillations and Flow Maldistribution (Npch = 10 ); Blue Curve - Density Wave Oscillations (DWO) boundary i.e, Hopf Curve, Region-A is unstable region; Black Curve - Flow Maldistribution (FMD) boundary i.e, Pitchfork Curve; Region-B is Symmetrical stable region; Region-C is the region of Flow Maldistribution i.e, Asymmetrical stable region. 6
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
Fig. 10. Numerical Simulation in region X1 of Fig. 4, perturbation leads to divergence of trajectory.
Fig. 8. Stability Map for Density Wave Oscillations (DWO) and Flow Maldistribution (FMD) for various Npch values; Each color represents the pair of DWO and FMD boundaries.
Fig. 11. Numerical Simulation in region X2 of Fig. 4, perturbation leads to convergence of trajectory.
Fig. 9. Stability Map for DWO and FMD with overlapping region(D) for Npch = 20 ; FMD boundary marked in blue color; DWO boundary marked in green color.
3.3.2. Region X2 As us+ is further increased in Fig. 4, Hopf point is observed. Beyond the Hopf point (marked as region X2), the equilibrium curve posses the real part of the largest eigenvalue to be negative (Fig. 5). It signifies stable system behaviour. If the system is initialized from an operating point in this region, a small perturbation results in the convergence of trajectory towards the fixed point. Time evolution of inlet velocity is shown in Fig. 11. 3.3.3. Region X3 As us+ is further increased in Fig. 4, BP(or pitchfork bifurcation) is observed. This marks the onset of FMD. From BP, two sets of equilibrium curve emanates. One corresponds to stable solution (marked as region Y1/Y2 in Fig. 4, asymmetric part, dotted lines) and other unstable solution (marked as region X3 in Fig. 4, symmetric part, solid line). The branches which emanates from the BP are the complementary solutions of the twin channel system. For symmetric solution, sign of the largest eigenvalue is positive (Fig. 5). If the system is initialized from an operating point in the region X3, trajectory diverges on account
Fig. 12. Numerical Simulation in region X3 of Fig. 4, perturbation causes trajectory to settle at an asymmetric stable solution. 7
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
Fig. 15. Variation of FLC vs Nsub (FLC is positive along the Hopf curve).
Fig. 13. Numerical Simulation in region Y1/Y2 of Fig. 4, perturbation leads to convergence of trajectories and are asymmetrically stable.
of trajectory. In region A, point 2 (Fig. 7) lies in the unstable region. Here, perturbation leads to divergence of the trajectory. Such is the nature of sub-critical Hopf region (section 2.3.1). This kind of identification of the sub-critical Hopf region around Hopf (DWO) boundary has been extensively studied in the past from the context of DWO. Therefore, the characterization of the Hopf stability boundary is not the interest of present work.
of small perturbation and settles on the solution of stable branches. Time evolution of inlet velocity is shown in Fig. 12. 3.3.4. Region Y1/Y2 For an asymmetric solution marked as Y1/Y2 in Fig. 4, sign of the largest eigenvalue is negative (Fig. 6). It signifies stable solution. If the system is initialized from an operating point in the region Y1/Y2, trajectory converges on account of small perturbation. Time evolution of inlet velocity is shown in Fig. 13.
3.3.6. Region D As us+ is increased we first observe FMD boundary marked in blue color in Fig. 9. If a system is initialized from an operating point marked in region A, we expect diverging DWO. Post occurrence of FMD boundary, we encounter region D. As we further increase us+ , DWO boundary marked in green color is observed. Such region D is observed for Nsub above 5. For lower Nsub values, symmetric stable solution region is observed. The same can be found in Fig. 7. Now, the region D as observed in Fig. 9, between FMD and DWO boundary has very interesting dynamics. Both the characteristics of FMD and DWO are observed here. If the system is initialized from an operating point in region D, system symmetry breaks and flowrate in the two channels become different. Flow maldistribution is thus observed first. After FMD boundary, it is the region of asymmetrical stable solutions. It is to be
3.3.5. DWO boundary characteristics The characteristics of the DWO curve (i.e, Hopf Curve) as shown in Fig. 7, is further analyzed by the plotting the FLC along the curve. It is shown in Fig. 14 and Fig. 15. The positive sign of FLC signifies the subcritical nature of the stability boundary. This implies the presence of an unstable limit cycle in the linearly stable region which emanates from the stability boundary. Here, small perturbation leads to convergence of DWO whereas large perturbation leads to divergence of DWO. In region B, point 1 (Fig. 7) lies in the stable region. Here, small perturbation leads to convergence of the trajectory. Also, when the system is initialized from the same point 1, large perturbation leads to divergence
Fig. 16. Numerical simulation from an operating point in region D of Fig. 9, presence of both the phenomenon i.e., FMD and DWO.
Fig. 14. Variation of FLC vs us+ (FLC is positive along the Hopf curve). 8
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
further noted from Fig. 7, that on the left of DWO curve, region A is unstable DWO region. So, post FMD in region D (Fig. 9), DWO in the system prevails. Thus dynamics of both flow maldistribution and density wave oscillations are observed in this region. It is to be noted that operating point of the system should be far away from this region. Here, symmetry break in the system on account of FMD is observed at the start of the transient. Furthermore, DWO is observed in the flow maldistributed system. Time evolution of inlet velocity is shown in Fig. 16.
pitchfork curve) demarcates symmetrical stable region and FMD region. Each point on the BP curve, there emanates multiple equilibrium curves. Stability of these curves are analyzed by the plotting the eigenvalues along the equilibrium curves. The region between the DWO curve (Hopf curve) and the FMD curve (pitchfork curve) represents the stable operating region and should be the region of candidature for operating point selection. For certain value of parameters, a fourth region is also observed. It arises because of the overlap between FMD boundary and DWO boundary. Interaction of DWO and FMD is observed here. Stability analyses using dedicated models which captures a particular instability only demarcates the parametric space into stable and unstable region for that instability. A combined analysis using a single model provides a demarcation for dominance of multiple instabilities in a single parametric space. Proper use of bifurcation analysis using a single model is shown here to analyze two different instabilities. Such unified approach provides a better perspective of system design and system instability to identify the region of stable operating points. The study also recommends a future scope to analyze multiple type of instabilities in a multi-channel flow boiling system.
4. Summary and conclusion This work analyzes the density wave oscillations (DWO) and flow maldistribution (FMD) in a combined fashion through the development of a stability map. Bifurcation analysis is used to identify various instabilities and their boundaries. For the first time, FMD is identified by pitchfork bifurcation in this work. Three regions in the stability map spanned by non-dimensional numbers Nsub and us+ are obtained by DWO and FMD boundaries. They are demarcated as DWO region, a region of the symmetrical stable solutions and FMD region. The DWO (Hopf) boundary demarcates the density wave oscillations and symmetrical stable region. Branch Point (BP) curve (i.e, FMD boundary or Nomenclature Acronyms
BP Branch Point DWO Density Wave Oscillations FLC First Lyapunov Coefficient FMD Flow Maldistribution GH Generalized Hopf Greek Symbols δ λ
Dirac delta function Boiling boundary (m) Single phase friction number
( ) Two phase friction number( )
1
f1 LH 2DH
f2 LH
2
2DH
Sign of the real part of largest eigenvalue ρ Density (kg / m3) Subscripts 0 Steady state e Channel exit f Liquid g Vapour H Heated channel i Channel inlet nth node in single phase region n Total Tot Other Symbols
P A Ax+ D f g h k L M Ns
s
Npch
Pressure drop (Pa) Area (m2 ) Non-dimensional cross sectional area of the jth heated channel Diameter (m) Friction factor Gravity acceleration (m2 /s ) Enthalpy (kJ / kg ) Resistance coefficient Length (m) Mass (kg ) No. of nodes in single phase region Phase change number
(
Q vfg W hfg vf
) 9
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
P Q q+ q t u us+ W W+ z Nsub x
Pressure (Pa) Heating power (W ) q Non-dimensional heat flux q0 2 Heat flux (W / m ) Time (s) Velocity (m /s ) Non-dimensional velocity, (ui,1 + ui,2) Mass flowrate Non-dimensional mass flowrate Axial coordinate (m) vfg hi ) ) Sub-cooling number ( h v (hf fg f Vapour quality
Appendix A. Derivation of mathematical model The detailed mathematical model is given in the work of Lin et al. [27]; Lee et al. [8] and Clausse & Lahey Jr [28]. However, a brief description is presented here. Mass Conservation Equation,
t
+
z
( u) = 0
(A.1)
Energy Conservation Equation,
t
( h) +
z
( hu ) =
q pH
(A.2)
AH
Momentum Conservation Equation, t
( u) +
z
i=1
f 2De
P z
( u2) =
N k (z i
z i)
u2 2
u2 g
(A.3)
here, ρ is density (kg / m3) ; t, time (sec ) ; u, velocity (m /s ); z, axial coordinate (m) ; f, friction factor; De , equivalent Diameter (m) ; k, loss Coefficient; δ, dirac delta function; g, gravity constant (m /s 2) ; P, Pressure (Pa) ; h, enthalpy (kJ / kg ) ; q", heat flux (W / m2) ; pH , heated perimeter (m) ; AH , crosssectional area (m2 ) . The density ρ can be written as,
=
(
= vf +
vfg
f,
(h
hfg
for h hf )
)
hf (Single 1
Phase region)
, for h > hf (Two
Phase region)
(A.4)
vf , specific volume of saturated liquid hfg , here, subscript f is fluid; vfg is difference in specific volume of saturated liquid and vapour latent heat of evaporation ( J / kg ). As shown in Fig. 1, the channel is divided into two regions, single and two-phase region. Entire channel is considered to be heated. Here, single phase region is divided into Ns number of nodes having variable length. Linear enthalpy profile is considered in single phase region between the two neighbouring nodes. So, the fluid enthalpy at the node can be defined as, (m3/ kg );
hn = h i +
n (h f
(m3/ kg );
hi ) (A.5)
Ns
here, subscript n is node in single-phase region; subscript i, inlet of channel; Ns , no. of nodes in single-phase region. A set of differential equations can be obtained by integration of energy balance equation (eqn. (A.2)) between two successive nodes (in singlephase region),
nth
Ln Ln 1 dLn dt
t
( h) dz +
= 2u i
Ln Ln 1
z
( hu) dz =
q p 2Ns A (hH h ) (Ln H f f i
Ln q pH Ln 1 AH
L n 1)
dz
dLn 1 dt
(A.6)
here, L is the length (m) ; L Ns , boiling boundary of the channel (m) . So, eqn (A.6) refers to the dynamic boiling boundary of the channel. The exit velocity ue can be calculated from eqn (A.1), eqn. (A.2) and density ρ (eqn. (A.4)),
q pH vfg u = = z AH hfg
(A.7)
Integration of the above equation from the boiling boundary to exit of the channel gives,
ue = ui +
(LH
L Ns )
(A.8)
here, subscript H is heated section of the channel; e, channel exit. Integration of the mass equation (eqn. (A.1)) from the inlet to the exit of the channel yields the following equation for the channel mass,
10
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al.
dMH = ( f ui dt
MH =
f
e ue ) AH
(A.9)
AH L Ns + M2
(A.10)
where, M is the coolant mass of the channel and subscript 2 is two phase. Also,
M2 = AH (LH
ln( f / e )
L Ns )
f
f/ e
1
(A.11)
Now, combining eqn. (A.9) and eqn. (A.10) we get ln
d e dt
=
f
ui
e ue
f / e
1
f
f / e
dL Ns dt
1
×
2 e 2 f (LH
L Ns )
f
e/ f
1
ln
e/ f
(A.12)
The dynamics of inlet velocity can be obtained by balancing the dynamic pressure drop in the channel and the external pressure head corresponding to the steady state. Pressure drop ( P ) in the channel is calculated by integrating the momentum balance equation (eqn. (A.3)), (A.13)
Pexternal = PI + Pg + Pf + Pa where, Inertial Pressure Drop, LH
PI = =
d dt
( u) dz t
0
(ue
MH ui AH
=
LH
d dt
ui) LH f
+
u dz
0
MH / AH
f / e
1
(A.14)
Gravitation Pressure Drop, LH
Pg =
g dz = gMH /AH
0
(A.15)
Frictional Pressure Drop, LH
Pf = =
f DH
0
L Ns f u2 dz 2DH
0
=
2 f ui L Ns
f1
2DH
2ui (ue ui) [ f f / e 1
(
f
3
e
)
+
+ +
i=1
N k (z i
LH f u2 dz L Ns 2DH
+
f2
L Ns
LH
2DH
1/ e
ki u2 2 i i
1/ f
MH / AH ] +
LH
ui2 ln
L Ns 2
ki
2 f ui
2
(LH
+
+ ke
L Ns )ln f / e
ke u2 2 e e
+
( )+ f
e
ui )2
(u e f / e
LH
u2 dz 2
zi)
2 f×
1
f / e
1
2 e ue
(A.16)
2
Acceleration Pressure Drop,
Pa = =
0 2 e ue
LH
u2 dz z 2 f ui
(A.17)
All the above equations are non-dimensionalized using the following definitions, A+ = A/ AH , non-dimensional cross-section area of the channel P + = P / f us2 , non-dimensional pressure drop t + = t /(LH /us ), non-dimensional time + = / , non-dimensional density f + u = u/ us , non-dimensional inlet velocity us = ui , characteristic velocity L+ = L / LH , non-dimensional channel length + = /L H , non-dimensional channel length z+ = z / LH , non-dimensional length
11
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al. h
h+ =
hf
Q0 /
f Aus
, non-dimensional enthalpy
M+ = M /( f LH AH ) , non-dimensional channel mass W + = W /( f Aus ) , non-dimensional mass flowrate Eu = Pexternal /( f us2 ) , Euler Number Fr = us2/(gLH ) , Froude number Subscript s refers to standard value(characteristic value). Hence, the conservation equations (refer eqn. (A.1), refer eqn. (A.2), refer eqn. (A.3)) can be written in non-dimensional form as, +
+
t+
t+ t+
z+
( +u+) = 0
( +h+) +
z+
( +u+) +
z+
i=1
(A.18)
( +h+u+) = 1
( +u+2) =
N k (z+ i
(A.19)
Eu
zi+)
P+ z
+u+2
+u+2
+
2
(A.20)
Fr
where, + = 1, for single phase region = (1 + Npch h+) 1, for two phase region
and = 1 , single-phase region; = 2 , two-phase region Fluid enthalpy at the node in the single phase region (refer eqn. (A.5)) can be written as,
Ns
hn+ =
n Ns
hi+
(A.21)
The differential equation of boiling boundary (refer eqn. (A.6)) can be written as,
dLn+ dt +
= 2ui+
2Ns
Npch Nsub
where, n = 1…Ns and
(Ln+
L N+s
Ln+ 1) +)
(or
dLn+ 1 dt +
(A.22)
is the dynamic boiling boundary. Npch =
Q vfg W hfg vf
, phase change number; Nsub =
vfg (hfg vf )
Integrating the mass balance equation (eqn. (A.18)) from inlet to exit of the channel yields (refer eqn. (A.9)),
dMH+ = ui+ dt +
+ + e ue
(hf
hi ) , sub-cooling number.
(A.23)
where (refer eqn. (A.10), eqn. (A.11)),
MH+ =
+
+ ln( e+) +) e + 1 e
+ (1
(A.24)
and (refer eqn. (A.8))
ue+ = ui+ + Npch (1
(A.25)
+)
The equation for dynamic fluid density (refer eqn. (A.12)) at the channel exit can be obtained from eqn. (A.23), eqn. (A.24) and eqn. (A.25), d e+ dt +
=
1+
+ + e ln( e ) + 1 e
(1 (1
+)[1
d + dt +
+
+ + e ue
ui+ ×
+ 2 e ) + + e + ln( e )]
(A.26)
Dynamics of inlet velocity can be obtained using (refer eqn. (A.13)), + Eu = PI+ + Pg+ + P + f + Pa
(A.27)
Using non-dimensionless form of inertial pressure drop (refer eqn. (A.14)), gravitation pressure drop (refer eqn. (A.15)), frictional pressure drop (refer eqn. (A.16)) and acceleration pressure drop (refer eqn. (A.17)), rearrangement yields,
dui+ (Eu PH+ ) = + + dt MH
(A.28)
where,
12
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al. +) (1 MH (1 / e+ 1)
PH+ =
+)
Npch (1 (1 / e+)
+
+ MH
Fr
(
+ 1+
+
1 ke 2
)
(
d +
d e+
Npch
dt +
+ 2 e ue
(
+
2ui+Npch (1
Npch (1 (1 / e+)
+) 2
1
1
+ e
+) MH
+
1 2
+) Mch
2 e)
)
1)
3
(1
ln(1/
ui+
dt +
1 ui2 +
+)(1
(1 / e+
+ dMH
+)(1
(1
ki 2
+ 1 u2 (1 / e+) 1 i
2
)
Npch dt + +
+u ++ i
1
+ e )+
+
+ MH+
+
(A.29)
Thus, the model can be represented as,
x = [Ln+
+ e
ui+] , n = 1…Ns
(A.30)
Appendix B. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.ijthermalsci.2019.106026.
References [15]
[1] Rizwan-Uddin, J. Dorning, Some nonlinear dynamics of a heated channel, Nucl. Eng. Des. 93 (1986) 1–14 http://www.sciencedirect.com/science/article/pii/ 0029549386901901 https://doi.org/10.1016/0029-5493(86)90190-1. [2] J. Boure, A. Bergles, L. Tong, Review of two-phase flow instability, Nucl. Eng. Des. 25 (1973) 165–192 http://www.sciencedirect.com/science/article/pii/ 0029549373900435 https://doi.org/10.1016/0029-5493(73)90043-5. [3] L.C. Ruspini, C.P. Marcel, A. Clausse, Two-phase flow instabilities: a review, Int. J. Heat Mass Transf. 71 (2014) 521–548 http://www.sciencedirect.com/science/ article/pii/S0017931013010922 https://doi.org/10.1016/j.ijheatmasstransfer. 2013.12.047. [4] S. Kakac, B. Bon, A review of two-phase flow dynamic instabilities in tube boiling systems, Int. J. Heat Mass Transf. 51 (2008) 399–433 http://www.sciencedirect. com/science/article/pii/S0017931007005959 https://doi.org/10.1016/j. ijheatmasstransfer.2007.09.026. [5] J. Lee, C. Pan, Dynamics of multiple parallel boiling channel systems with forced flows, Nucl. Eng. Des. 192 (1999) 31–44 http://www.sciencedirect.com/science/ article/pii/S0029549399000850 https://doi.org/10.1016/S0029-5493(99) 00085-0. [6] X. Lu, Y. Wu, L. Zhou, W. Tian, G. Su, S. Qiu, H. Zhang, Theoretical investigations on two-phase flow instability in parallel channels under axial non-uniform heating, Ann. Nucl. Energy 63 (2014) 75–82 http://www.sciencedirect.com/science/ article/pii/S030645491300385X https://doi.org/10.1016/j.anucene.2013.07.030. [7] A.M. Mishra, S. Singh, Non-linear stability analysis of uniformly heated parallel channels for different inclinations, Appl. Therm. Eng. 98 (2016) 1189–1200 http:// www.sciencedirect.com/science/article/pii/S1359431115013666 https://doi.org/ 10.1016/j.applthermaleng.2015.11.118. [8] J.-D. Lee, Y.-G. Lin, S.-W. Chen, C. Pan, The influence of void-reactivity feedback on the bifurcation phenomena and nonlinear characteristics of a single nuclear-coupled boiling channel, Ann. Nucl. Energy 94 (2016) 814–825 http://www. sciencedirect.com/science/article/pii/S0306454916302110 https://doi.org/10. 1016/j.anucene.2016.04.044. [9] T. Van Oevelen, J.A. Weibel, S.V. Garimella, Predicting two-phase flow distribution and stability in systems with many parallel heated channels, Int. J. Heat Mass Transf. 107 (2017) 557–571 http://www.sciencedirect.com/science/article/pii/ S0017931016324395 https://doi.org/10.1016/j.ijheatmasstransfer.2016.11.050. [10] K. Akagawa, M. Kono, T. Sakaguchi, M. Nishimura, Study on distribution of flow rates and flow stabilities in parallel long evaporators, Bull. JSME 14 (1971) 837–848, https://doi.org/10.1299/jsme1958.14.837. [11] U. Minzer, D. Barnea, Y. Taitel, Evaporation in parallel pipes splitting characteristics, Int. J. Multiph. Flow 30 (2004) 763–777 http://www.sciencedirect.com/ science/article/pii/S0301932204000527 https://doi.org/10.1016/j. ijmultiphaseflow.2004.04.006. [12] U. Minzer, D. Barnea, Y. Taitel, Flow rate distribution in evaporating parallel pipes―modeling and experimental, Chem. Eng. Sci. 61 (2006) 7249–7259 http:// www.sciencedirect.com/science/article/pii/S0009250906005173 https://doi.org/ 10.1016/j.ces.2006.08.026. [13] M. Baikin, Y. Taitel, D. Barnea, Flow rate distribution in parallel heated pipes, Int. J. Heat Mass Transf. 54 (2011) 4448–4457 http://www.sciencedirect.com/science/ article/pii/S0017931011002535 https://doi.org/10.1016/j.ijheatmasstransfer. 2011.04.034. [14] Y. Taitel, D. Barnea, Transient solution for flow of evaporating fluid in parallel pipes using analysis based on flow patterns, Int. J. Multiph. Flow 37 (2011) 469–474
[16]
[17]
[18]
[19]
[20]
[21] [22]
[23]
[24]
[25]
[26]
[27]
[28]
13
http://www.sciencedirect.com/science/article/pii/S030193221100022X https:// doi.org/10.1016/j.ijmultiphaseflow.2011.01.002. V. Pandey, S. Singh, Characterization of stability limits of ledinegg instability and density wave oscillations for two-phase flow in natural circulation loops, Chem. Eng. Sci. 168 (2017) 204–224 http://www.sciencedirect.com/science/article/pii/ S0009250917302865 https://doi.org/10.1016/j.ces.2017.04.041. V. Pandey, S. Singh, Bifurcation analysis of density wave oscillations in natural circulation loop, Int. J. Therm. Sci. 120 (2017) 446–458 http://www.sciencedirect. com/science/article/pii/S1290072916312248 https://doi.org/10.1016/j. ijthermalsci.2017.06.029. I. Simanovskii, A. Viviani, F. Dubois, Nonlinear dynamics of two-layer systems with a temperature-dependent heat release, Int. J. Therm. Sci. 119 (2017) 51–67 http:// www.sciencedirect.com/science/article/pii/S1290072916317975 https://doi.org/ 10.1016/j.ijthermalsci.2017.05.022. K. Luo, H.-L. Yi, H.-P. Tan, Eccentricity effect on bifurcation and dual solutions in transient natural convection in a horizontal annulus, Int. J. Therm. Sci. 89 (2015) 283–293 http://www.sciencedirect.com/science/article/pii/S1290072914003329 https://doi.org/10.1016/j.ijthermalsci.2014.11.020. S. Paruya, S. Maiti, A. Karmakar, P. Gupta, J. Sarkar, Lumped parameterization of boiling channel―bifurcations during density wave oscillations, Chem. Eng. Sci. 74 (2012) 310–326 http://www.sciencedirect.com/science/article/pii/ S0009250912001352 https://doi.org/10.1016/j.ces.2012.02.039. D.N. Basu, S. Bhattacharyya, P. Das, Performance comparison of rectangular and toroidal natural circulation loops under steady and transient conditions, Int. J. Therm. Sci. 57 (2012) 142–151 http://www.sciencedirect.com/science/article/pii/ S1290072912000567 https://doi.org/10.1016/j.ijthermalsci.2012.02.011. A. Karmakar, S. Paruya, Nonlinear analysis of chaotic time series in a natural circulation boiling loop, Ind. Eng. Chem. Res. 51 (2012) 16467–16481 doi:10.1021/ ie300788t https://doi.org/10.1021/ie300788t. M. Krishnani, D.N. Basu, On the validity of boussinesq approximation in transient simulation of single-phase natural circulation loops, Int. J. Therm. Sci. 105 (2016) 224–232 http://www.sciencedirect.com/science/article/pii/S1290072915300065 https://doi.org/10.1016/j.ijthermalsci.2016.03.004. S. Paul, S. Singh, On nonlinear dynamics of density wave oscillations in a channel with non-uniform axial heating, Int. J. Therm. Sci. 116 (2017) 172–198 http:// www.sciencedirect.com/science/article/pii/S1290072916307761 https://doi.org/ 10.1016/j.ijthermalsci.2017.02.008. G. Xia, G. Su, M. Peng, Analysis of flow distribution instability in parallel thin rectangular multi-channel system, Nucl. Eng. Des. 305 (2016) 604–611 http:// www.sciencedirect.com/science/article/pii/S0029549316300449 https://doi.org/ 10.1016/j.nucengdes.2016.04.016. A. Dhooge, W. Govaerts, Y.A. Kuznetsov, MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs, ACM Trans. Math Software 29 (2003) 141–164 doi:10.1145/779359.779362 http://doi.acm.org/10.1145/779359. 779362. F. Alobaid, N. Mertens, R. Starkloff, T. Lanz, C. Heinze, B. Epple, Progress in dynamic simulation of thermal power plants, Prog. Energy Combust. Sci. 59 (2017) 79–162 http://www.sciencedirect.com/science/article/pii/S0360128516300375 https://doi.org/10.1016/j.pecs.2016.11.001. Y. Lin, J. Lee, C. Pan, Nonlinear dynamics of a nuclear-coupled boiling channel with forced flows, Nucl. Eng. Des. 179 (1998) 31–49 http://www.sciencedirect.com/ science/article/pii/S0029549397002422 https://doi.org/10.1016/S00295493(97)00242-2. A. Clausse, R. Lahey Jr., An investigation of periodic and strange attractors in boiling flows using chaos theory, Heat Transfer 1990. Proceedings of the Ninth
International Journal of Thermal Sciences 145 (2019) 106026
D. Paul, et al. International Heat Transfer Conference, 1990. [29] P. Saha, M. Ishii, N. Zuber, An experimental investigation of the thermally induced flow oscillations in two-phase systems, ASME J. Heat Transf. 98 (1976) 616–622 http://heattransfer.asmedigitalcollection.asme.org/article.aspx?articleid= 1436574 https://doi.org/10.1115/1.3450609. [30] Rizwan-Uddin, J. Dorning, Some nonlinear dynamics of a heated channel, Nucl. Eng. Des. 93 (1986) 1–14 http://www.sciencedirect.com/science/article/pii/ 0029549386901901 https://doi.org/10.1016/0029-5493(86)90190-1. [31] S. Paul, S. Singh, A density variant drift flux model for density wave oscillations, Int. J. Heat Mass Transf. 69 (2014) 151–163 http://www.sciencedirect.com/ science/article/pii/S0017931013008727 https://doi.org/10.1016/j.
ijheatmasstransfer.2013.10.012. [32] M. Ishii, Thermally Induced Flow Instabilities in Two-phase Mixtures in Thermal Equilibrium, Georgia Institute of Technology, 1971, https://books.google.co.in/ books?id=eWurtgAACAAJ. [33] P. Saha, Thermally Induced Two-phase Flow Instabilities, Including the Effect of Thermal Non-equilibrium between the Phases, Georgia Institute of Technology, 1974, https://books.google.co.in/books?id=oc2RbwAACAAJ. [34] Y. Kuznetsov, Elements of Applied Bifurcation Theory. Applied Mathematical Sciences, Springer, New York, 2004https://books.google.co.in/books?id= B6Z2TQSRJuMC.
14