Flow boiling heat transfer characteristics using the modified eulerian and wall heat balance model

Flow boiling heat transfer is routinely encountered in nuclear reactors, steam engines and other engineering applications. Although several researchers have carried out different numerical and experimental investigations on flow boiling, the underlying physics of the interfacial interaction is still a complex phenomenon to understand in detail. Hence, the numerical simulation and optimisation regarding the adoption of engineering flow boiling parameters have been conducted in this study using the Modified Eulerian-Eulerian Model (MEEM) and Wall Heat Balance Model (WHBM). To predict interpenetrating flow fields and to provide detailed relevant information on the flow behaviour, this study considered a uniform axial heating profile for a cylindrical flow channel. The Raynolds Average Navier Stokes (RANS) equation with an appropriate turbulence model are used to predict the effect of turbulence on the mean flow field, while the MEEM multiphase sub-models are employed to predict the temperature distribution along the wall, the average void fraction, tracking of the single bubble detachment diameter, heat balance at the wall, effect of surface roughness on heat transfer, the effect of aspect ratio, and the critical heat flux. The results obtained from this study are compared with the selected numerical investigation and experimental data presented in the open literature. The present study shows a better approximate prediction (with minimal uncertainties) of both the subcooled boiling heat transfer and the saturated boiling heat transfer. In summary, this study agrees with extant theories and experimental predictions. Thus, it has provided more profound insights into flow boiling heat transfer particularly for flow in a vertical pipe.


Introduction
Flow boiling occurs when a fluid circulates over a heated surface by external means such as a pump or due to the natural buoyancy effect [1]. Several flow patterns characterise such a flow field, and the transition between flow regimes dependent upon the channel geometry and thermophysical properties of the fluid [2]. A continuous flow of liquid and vapour in a two-phase flow boiling and heat transfer encounters some flow resistance, that causes the fluid pressure to drop at a certain height (i.e., the non-boiling and the boiling height) which is one of the main factors of characterising the flow boiling in a channel. The non-boiling height is the singlephase height in which the subcooled liquid entering the channel receives a quantity of sensible heat. This sensible heat causes a change in subcooled temperature, resulting in subcooled boiling without any significant effects on the pressure drops and fluid density. In the boiling height, the pressure drop is higher than that in a non-boiling height for the same flow rate [2]. This higher pressure drop is a result of an increase in the flow velocity along the channel and reduction in cross-sectional area due to the presence of both liquid and gas phases.
Flow boiling heat transfer is primarily affected by factors such as the inertia forces, viscous forces, pressure forces, interfacial tension forces, liquid-surface contact angle, exchange of mass, momentum, and energy between the interacting interface between the liquid and vapour phases. These interfaces are bubbles in a continuous liquid flow, a droplet in continuous vapour flow, a vapour film in continuous liquid flow and a liquid film in continuous vapour flow [3]. Due to a large number of dependencies, it is not possible to obtain a single heat transfer correlation for the entire flow regimes. The nature of the flow regimes depends on the operating conditions such as pressure, temperature, mass flux, channel orientation, and fluid properties [2,4].
The onset of sub-cooled nucleate boiling initiates when the wall temperature at the single-phase forced convection region reaches to the saturation temperature at a given pressure, and also the wall superheats necessary to cause nucleation. The amount of superheat necessary for nucleate boiling to occur is dependent on the inner surface heat flux. Hence, no boiling can occur when the inner surface temperature is less than the saturation temperature at a given location along the channel [2]. Therefore, subcooled boiling begins and develops when the bulk temperature or the mean enthalpy of the liquid phase is below the saturation temperature. In the highly subcooled region, bubbles grow and collapse while attached to the heated surface and do not penetrate far into the bulk subcooled flow. Bubbles grow and detach in the low subcooled region because the void fraction in this region rises with its length. The surface tension and inertia forces are responsible for the bubble formation in a highly subcooled region, which also holds the bubble to the heated surface. Buoyancy and the frictional (drag) forces are responsible for the bubble detachment from a heated surface in a low subcooled region [2].

Background literature
In order to understand the state-of-the-art of flow boiling, a comprehensive review of the background literature has been conducted. Flow boiling research has seen significant advances in recent years with particular attention given by the nuclear industry. Over the years, numerous numerical and experimental techniques were employed to understand flowboiling phenomena. For instance, experimental investigations have used a combination of high-speed cameras and image processing techniques to visualise the different flow regimes and analyse the different bubble dynamic parameters. Numerical simulations were applied to study and visualise the bubble dynamics and different flow parameters for different working fluids and other operating conditions. Analytical approaches based on theoretical correlations have been widely used for studying the interfacial transport phenomenon. This approach is particularly vital in the subcooled boiling region where bubble nucleation at the interfacial boundary is the essential heat transfer mechanism. This method also explores the use of interfacial area equation to study particle interaction at the interface and phase change.

Experimental framework
Prodanovic et al., [5] performed experimental studies using a vertical annular test section to observe the effect of parameters such as subcooling, pressure, and flow rate on bubble behaviour in subcooled flow boiling. High-speed photography at rates between 6-8 kHz revealed at low heat flux, the bubble shape is spherical and does not change in size with the flow distance. The final region termed the bubble coalescence region was characterised by large bubbles of different sizes and shapes. Work by Michel et al. [6] experimently investigated the effect of gravity on heat transfer of upward flow boiling. This experiment was performed by varying the gravity levels while making local temperature measurements and observing the different flow transitions. The effect of gravity on heat transfer decreases as the mass flux and the heat flux increases. Forced convection was experienced in high mass flux while film evaporation was seen dominant during applied high heat flux. Burak et al. [7], carried out an experimental investigation of flow boiling using a single 1.8 mm × 0.9 mm microchannel at low mass velocities to ascertain corresponding effects on the bubble dynamics. Their experimental observations indicated that the rate at which heat transferred from the heated surface to the main flow stream is dependent on the heat flux and vapour quality, and it increases with increasing mass flux. Recent work on flow boiling by Colgan et al., [8] studied the Critical Heat Flux (CFH) under subatomic atmospheric conditions. Based on their findings, the mass flux and the inlet subcooling do not have any substantial effect on the critical heat flux, and also the wall superheat is prolonged at the subatmospheric conditions compared to the atmospheric conditions.
Experimental studies of Zhu et al. [9] on natural circulation of two-phase flow resistance in a vertical 3 × 3 rod bundle showed a rapid occurrence of flow instability under natural circulation conditions. This instability results in a rapid increase in acceleration pressure drop and frictional pressure drop under a low mass quality condition. With the gravity, channel pressure drop decreases under the same condition. Sira et al., [10] investigated the effect of channel orientations on flow boiling behaviours. As expected, their investigation was in agreement with the expected results as achieved by previous researches over the decades. They specifically examined the effect of change in channel orientation on pressure drop, flow pattern, and heat transfer coefficient. Cheng et al. [11] studied the effect of diameter on two-phase gas-liquid flow using a small tube of 28.9 mm in diameter and a larger tube of diameter 150 mm. They reported that slug flow does appear only in the pipe with small diameter but did not appear in large diameter. Chen et al. [12] studied the two-phase flow of R134 refrigerant in vertical mini pipes and found out that the transition regions of slug-churn and churn-annular, strongly depend on tube the diameter. However, the transition boundaries for bubbly to churn and bubbly to slug are slightly affected by a change in the magnitude of diameter. Furukawa et al. [13] studied the effect of viscosity on transition regions between flow patterns in vertical tubes with inner diameter of 19.2 mm and height of 54 mm. They used water and aqueous glycerol solutions with viscosities of 1 × 10-6 up to 14.7 × 10-6 m2/s and reported that an increase in liquid viscosity shifts the bubbly-slug transition boundary to lower superficial air velocities and the froth-annular transition boundary to regions of higher superficial viscosity.
Taylor performed visualisation tests on a vertical annulus to investigate the bubble departure frequency [14]. Conclusions from this study established that heat flux, mass flux, degree of subcooling, and system pressure are the four major thermalhydraulic parameters which control the bubble departure frequency. Their investigation shows that the bubble departure frequency increases with an increase in heat flux and pressure, but it decreases with the increase of mass flux and subcooling. Thorncroft et al. [15], measured the waiting time, bubble growth, and bubble departure size in an upward vertical subcooled flow-boiling channel. They correlated the waiting time against the wall superheat and the growth time with the bulk subcooling.

Computational framework
A validation and assessment study to predict the critical heat flux in flow boiling inside a 4 × 4 rod bundle was carried out by Harish et al. [16]. They tested the effect of axial power distribution by the applied heat flux, which affects the temperature distribution that, in turn, affects the coolant phase change and the critical heat flux. From detailed and systematic validations, mean relative errors of 15% and 7% for tube and annulus cases respectively were recorded. These simulations identify the critical heat flux (CHF) location in an annulus and a rod bundle. The wall temperature variations predicted for the concentric and eccentric annulus as well as a 4 × 4 rod bundle geometry. Eccentricity found to promote the early occurrence of CHF in the narrow gap and delay in the broad gap regions. For the 4 × 4 rod bundle considered, the corner rod facing the pressure tube was found to have an early shoot up in the wall temperature. The effect of different turbulence model was also carried out and they concluded that RNG k-e accurately predicts both the CHF and the post-dryout wall temperature; while standard k-e and SST k-w models underpredicted and overpredicted the peak wall temperature, respectively.
In order to understand the nature of bubble behaviour at the reactor's operating conditions, Eyitayo et al. [17], carried out a numerical study using the Large Eddy Simulation (LES)-micro layer model to investigate the growth rate, the microlayer thickness, and distortion of a single bubble in a highly turbulent flow boiling. They investigated the effect of system pressure, contact angle, initial bubble size, and bulk velocity on bubble growth. For the analysis of bubble microlayer thickness, the results showed that the microlayer thickness increases with the increase of both the radius of the micro-region and the system pressure. Furthermore, they studied the bubble growth rate for different operating pressures. In particular, a slow growth rate was experienced at system pressures between 1-21 MPa while the growth was rapid at higher pressures due to the decrease in liquid-vapour density ratio and surface tension. This report shows that using Large Eddy Simulation (LES) to analyse the bubble microlayer and growth in subcooled flow boiling ensured a faster computational speed with an excellent accuracy than carrying out a full Direct Numerical Simulation (DNS) analysis. Overall, this study covers the operating conditions for a boiling water reactor (BWR) and pressurised water reactor (PWR).
Adrian et al. [18], studied the two-phase flow of air/water mixture using OpenFOAM and compared their findings with the mechanistic model developed by Petalas et al. [19]. Their findings showed a good correlation with the mechanistic approach as all the flow regimes observed were consistent with the mechanistic prediction. Fangfang et al. [20] studied the instability mechanisms leading to slug formation by using Direct Numerical Simulation (DNS) analysis of two-phase flow in an inclined pipe. They identified different flow regimes (i.e., stratified, slug, dispersed, and annular) and the findings were compared with the theoretical models. Nemitallah et al., [21], carried out an intensive numerical study to analyse the current status of two-phase flow boiling characteristics using the Eulerian Multiphase flow boiling model. In their research, they considered the influence of a non-uniform axial heating profile for a highly pressurised system of 2 m long with inner and outer diameters of 15.4 mm and 25.4 mm respectively. The mass flow rate and the inlet temperature of the water was 0.161 kg/s and 200 o C, respectively. Francisco et al. [22], also carried out a numerical investigation to validate the subcooled flow boiling in a vertical channel using the model developed by Kurul et al. [23]. They considered several factors (such as the mesh aspect ratio, heat flux, system pressure, and thermodynamic quantities.) which could enhance the active nucleation along the heated wall. Their numerical investigation was validated using the benchmark of the experimental work of Chenturiya [24]. Adrian et al. [25], modelled two-phase flow boiling in a Boiling Water Reactor Fuel Assembly using a CFD code STAR-CD which has been seen to have the capability of tracking the mass, momentum, and energy of the liquid and vapour phases in each computational cell. Huiying Li et al. [26] modelled the boiling flow within the framework of the Eulerian multifluid model, using the RPI model for nucleate boiling and its extension to the non-equilibrium boing and critical heat flux regime. Mukesh et al. [27], of the Bhabha Atomic Research Centre, Mumbai, India, carried out a CFD simulation to predict the void fraction and temperature distribution inside the rod bundle to identify the growth of any hot spots during the removal of decay heat in an Advanced Heavy-Water Reactor (AHWR). This boiling water reactor is a natural circulation based reactor with many passive safety features and utilizes an isolation condenser system for decay heat removal. It is believed that, during the removal of the decay heat, there is a possibility of a hot spot occurring inside the fuel rod cluster. Hence, a boiling Eulerian model was developed to obtain the flow field during the decay heat removal by first considering the effect of different models on the flow field.
Multiphase flow boiling modelling using CFD is very challenging due to the inherent complications such as the instabilities arising from the turbulent mixing of the different phases that have different characteristic properties, the wallfluid interaction encountered during the process, and the interfacial heat and mass transfer in the two-phase flow-mixing region. Despite the numerous research, the two-phase flow boiling phenomenon still remains a challenging case study in the nuclear industry and other engineering applications. Notwithstanding, the importance of CFD cannot be overemphasised because it is an indispensable tool used in the absence of an experimental setup to study and understand several complex structures in fluid flows. This present study is aimed at exploring the capabilities of Computational Fluid Dynamics (CFD) using the modified eulerian and wall heat balance model [28] to optimise the accuracy of some predicted flow boiling parameters. The modified Eulerian two-fluid coupled model is more robust and can predict flow boiling parameters more accurately in an upward flow through a heated vertical channel. The results in this study can be used as a benchmark for the validation of future study.

Modelling description
The pictorial representation, as shown in Fig. 1, explains the concept established during vertical flow boiling. In practical, the thermodynamic nonequilibrium experienced in two-phase flow, causes the actual vapour quality to be less than the equilibrium quality. The present study has been implemented using ANSYS Fluent. Some key information about modelling with ANSYS Fluent cand be found in the user guide manual [52].

Turbulence model
The Raynolds Average Navier-Stokes (RANS) method targets resolving turbulence quantities, and the gradients of the Renolds stress tensor components entering the equation of motion represent the total stress tensor in a fluid obtained from the averaging operation over Navier-Stokes equations to account for turbulent fluctuations in fluid momentum. This component is of decisive importance for the correct capturing of the flow physics. Equation (1), is a modified standard k − ϵ turbulence model with an additional dissipation energy term (R ϵ ) which account for the effective viscosity. This model provides information for low Reynolds number effects and improves the accuracy for rapidly strained flows [1].

Multiphase model
The baseline wall boiling model used for this simulation is the Two-Phase Eulerian model developed by Kurul et al. [28]. This wall heat balance model, as shown in Eq. (2), is used to quantify the wall-to-fluid heat transfer. It does not automatically compute the vapour temperature when solving the phase energy equation; instead, the model partitioned the total heat flux H B (W) from the wall to the fluid free stream into three components namely: the liquid phase convective heat fluxq, the quenching heat fluxq Q ð Þ and the evaporation heat fluẋ To account for the non-equilibrium boiling and critical heat flux, the value for the vapour temperature, which corresponds to the saturation temperature at a given operating system, the pressure is inputted using the user defined function (UDF). Thus, Eq. (3), shows the modified wall heat-partitioning model which was used in this study.
where; α_(v,2)=0.95 and α_(v,1)=0.90 Table 1 shows different models available in the published literature to predict the bubble departure diameter (d b ) which is one of the most significant dynamic parameters observed during flow boiling heat transfer. Most of the existing bubble departure diameter correlations (such as Fritz [32], Cole and Rohsenow [33], Tolubinski and Kostanchuk [34], Kocamustafaogullari-ishii [35], Klauser et al. [36] and Situ et al. [29]) were formulated based on the balance of forces acting on the bubble at the departure. Several researchers, such as Huiying Li et al. [26], used the Tolubinski & Kostanchurch model for their research. Unal correlation [37] is considered in this present study because it depends solely on the degree of local subcooling and the wall superheat, and also it shows the best agreement with the experimental data reported by Situ et al. [29]. The general formulation for bubble departure diameter is as shown in Eq. (4).
During the transition from a bubbly to a mist flow, the d r o pl e t d i am et er ( D d ) w a s e s t im a t e d us i n g t h e Kocamustafaogullari-Ishii correlation provided in Eq. (5) [38].
Where C ds is a coefficient, which has a value of about 0.028, and j v is the superficial velocity.
The active nucleation site density, defined as the number of nucleation sites per unit area, has been resolved in the current study using the Lammet correlation [39]. Kocamustafaogullari-Ishii correlation [40] was developed to  Table 1 Bubble Departure Diameter Models [29][30][31] Model Cole and Rohsenow C fs S: Cl ur p Pr f 2 account for cases of forced convective boiling at high pressure, the Lammert-Chawla model as used in this current research gives a more precise, accurate and stable solution as compared to the Hibiki and Ishii model [41] used by Colombo and Fairweather [42] in their research. All the published models have provided the nucleation site density (N w ) as a function of the wall superheat, as shown in Eq. 6. Where C and m are dimensionless constants with values 1.805 and 15,545.54, respectively.
The drag coefficient ( f) was calculated using the Universal Drag Function [43], computed based on the relative Reynolds number as shown in Eq. (7) Also, the disperse bubbly, and mist flows were also calculated using Eq. (8); In the bubbly flow regime, the lift coefficient (C L lv ) was calculated using the formulations of Moraga et al. [26]. However, this coefficient is assumed to be zero for droplet flows. The general formulation, as used for the simulation, is shown in Eq. (9), The wall lubrication force ( F ! wl ) acts as a repulsive force along the wall by pushing the bubbles away from the wall. Previous studies as reported in the literature used the Antal et al. model [25] [26] [22] whereas this present study used the model proposed by Hosokawa et al. [44] and the wall lubrication coefficient (C wl ) value of 0.0217 was calculated using Eq. (10). Also, this model employs the Tomiyama formulation, which is defined based on the channel hydraulic diameter.
Lopez-de-Bertadano model [45] was considered for the turbulence drift force (see Eq. 11), although several authors used Burns-et-al Model [46] as recommended in the literature.
Where ke c is the turbulent kinetic energy. During flow boiling, there is an interfacial heat transfer (i.e., the interface to liquid transfer and the interface to vapour transfer). At the onset of the nucleate boiling, the bubble departure from the wall transport both heat and mass from to bulk liquid. This phenomenon is known as the interface-liquid transfer and numerically by Eq. (12); Where h lit is the heat transfer coefficient of the liquid interface. This study used two resistance Phase-Interphase heat transfer models switching from Ranz-Marshal to Hughmark formulation [47], shown in Eq. (13).
Similarly, the interface to vapour heat transfer is calculated based on the Lavieville et al. correlation [48], which proposed that vapour retains the saturation temperature by rapid evaporation/condensation. Equation (14) is a mathematical expression for vapour phase heat transfer.
Where δt is the time scale. Also, the wall-vapour mass transfer (evaporation mass flow,ṁ E ) is calculated based on the evaporation heat fluẋ q E defined in Eq. (15); Furthermore, the interfacial mass transfer relationship used in the present study is given by Eq. (16);

Boundary conditions
Appropriate thermal boundary conditions specified for all flow boundaries and all walls. Setting up the correct boundary conditions is very important due to the relevance in the drag and lift exerted on bodies as well as heat and mass transfer across boundaries. Specifically, when considering continuous vapour and liquid phases where evaporation and condensation occur, their boundary conditions becomes very complicated as they require comprehensive information about the molecular phenomena at the interface [49].
The following initial conditions used for our numerical analysis; At the Inlet; the velocity and the temperature of the working fluid as it enters the domain were specified to be 1.04 m/s and 474.15 K respectively. While at the outlet, a pressure outlet was specified and the value of 530.15 K was used as the backflow temperature to account for any reverse flow occurs during part of the iterative process. This value can only use it if there is reverse flow, but if there is no reverse flow, the backflow value will not have any effect on the solution.
The heat flux of 0.57 MW/m2 was specified at the wall boundary. By defining the heat flux, the CFD algorithm will calculate what the wall temperature will be in other to achieve that heat flux.
Turbulence settings for near-wall boundary layer are dependent upon the turbulence approach used. The boundary layer is a thin layer of fluid adjacent to the solid surface in whuch the velocity develops from zero at the wall to the free stream value at soe distance away from the wall. Hence, it is essential to use a very fine mesh to resolve the steep profile within the boudary layer because the frictional drag on the surface and the heat transfer takes place within the boundary layer. The mesh resolution requirements are normally expressed in terms of the dimensional wall distace, which should be maintained within a acceptable range for the model used to achieve accurate results. The enhanced wall treatment for the k-epsilon realisable model used because this near-wall treatment will automatically switch between the wall function and viscous sublayer resolution based on the value of y + . Also, it will blend between the two if the first grid point happens to fall above the buffer region. With the RANS model, no turbulence was resolved in the simulation; instead, the mean flow field was calculated. The turbulence specification method used for this study is the intensity and hydraulic diameter. Turbulence intensity is a measure of the strength of turbulence fluctuations. The value of 0.05 used in this study means that the velocity of the turbulence flunctuations is 5% of the mean velocity.
The time step size Δt was set to be small enough to resolve time-dependent behaviour in the flow and was calculated using Eq. (17).
L=the characteristic length scale, β,ΔT, is the thermal properties of the fluid Also, the y + value estimated using Eq. (18);

Results and discussion
This section presents the current case study analysis for the numerical framework and results validated with the most re-

Steady state analysis
Sensitivity analysis has been conducted showing the effect of mesh resolution and turbulence on the flow field.
Mesh dependence A 2D and 3D Hexahedral mesh presented in Fig. 2 was created in such a way that the grid points are close enough together to predict the regions where there are high gradients in the flow or temperature field. Also, such that all the geometric features of interest resolved. These were performed for a different number of grid points, as shown in Fig. 3. The results show that the number of grid points has an overriding on the entry length and the void fraction at the outlet. The results obtained with the refined mesh gave a shorter entry length, (hence longer computational times and large memory consumption), and a high vapour volume   5 shows the output void fraction for different test cases and compared with the experiment and the tracking of the single bubble departure diameter for three test cases. This study explains the flow dynamics when liquid moves along the heated pipe. The fluid near the heated wall gains thermal energy, thereby resulting in the temperature distribution from the heated surface along the fluid stream. Though the temperature of the fluid increased along the pipe, the primary fluid stream will not boil unless it reaches the saturation temperature, which in this case is 530.15 K. The length of transition to nucleate boiling is where the fluid gains thermal energy albeit not enough to cause bubble nucleation. This phenomenon has been investigated for different test cases as tabulated in Table 2. Hence, the onset of nucleate boiling begins when the wall temperature and the fluid near the wall are the same as the saturation temperature. Bubble growth also begins at this stage until detachment occurs at the onset of fully developed subcooled nucleate boiling (ON-FDSCNB). Spontaneous bubble growth and detachment continues as the vapour nucleus attains a size greater than that for unstable equilibrium, as observed in Fig. 5. Figure (6, 7 and 8) show a comparative analysis of the subcooled boiling length and the average volume fraction reported by Adrian et al. [18], Fransisco el al. [22], Huiying Li et al. [26], compared with the current numerical study and experiment. The values, as reported in  from the experiment with the non-boiling length of 0.4 m, which is perfect for a stable, steady-state condition. The previous numerical investigation and experimental data used as a reference for the current study introduced the wall heat flux at some distance (0.4 m) away from the channel inlet. This adiabatic entry length accounts for a single-phase flow region where there is no transfer of heat between the thermodynamic system and the environment, neither there is a change in the temperature of the subcooled liquid flowing along this region. The onset of the single-phase convective heat transfer began at the point where the uniform wall heat flux is introduced up to about 1 m from the channel inlet. The set criteria results to an increase in the ratio of the entry length to the channel length, which inturns underpredicted the average void fraction at the channel outlet. Whereas, the current simulation assumed that the entire channel length from the inlet is uniformly heated and the wall temperature at the inlet is assumed to be the same as the liquid inlet temperature. It is thereby resulting in a smooth transition of the heating effect up to the onset of incipient boiling and better prediction of the average void fraction.  Figure 9 and Fig. 10 show the centreline temperature, and the degree of wall superheat, respectively. The flow channel shows that the wall temperature increases exponentially along the vertical tube. At the channel inlet, the temperature of the primary phase is 474.15 K. As the liquid moves along the heated pipe, it gains thermal energy from the heated wall thereby resulting in the temperature distribution from the heated surface along the fluid stream starting from the liquid near the heated wall. Though the temperature of the fluid increases along the pipe, the primary fluid stream will not boil unless it reaches the saturation temperature, which in this case is 530.15 K. The entry length where the fluid gains thermal energy, but it is not enough to cause bubble nucleation known as the non-boiling height. Hence, the onset of nucleate boiling begins (if and only if) the wall temperature and the fluid near the wall are the same as the saturation temperature. Spontaneous bubble growth continues when the vapour nucleus attains a size greater than that for unstable equilibrium. Previous studies show that for a normal operating condition of a water reactor heat transport channel, the wall superheat should not be higher than 10 K [2]. The result obtained in this study falls within the benchmark reported in the literature, and it is in good agreement with the experiment.  Figures 11 and 12 illustrate the significant effects of surface roughness on the total pressure drop and heat transfer coefficient which determines the amount of vapour volume fraction at the channel outlet. In Fig. 11; P1, P2, W1 and W2 represent the average pressure drops along the mainstream liquid phase, in the mainstream vapour phase, along the liquid phase pipe wall, and along the vapour phase, respectively. When a fluid is flowing along a heated pipe, the heat transfer and flow behaviour vary with the heat flux and the conditions of the fluid. Where the fluid temperature reaches saturation point, there is an increase in volume and velocity of the mixture. The ressure loss encountered due to the increase in momentum of the mixture as it flows through the channel and vapourises, and increase in frictional pressure drop due to the effective roughness of the inner channel. The pressure loss is also caused by the buoyancy effect and can lead to an increase in gravitational potential energy as the boiling mixture rises in the channel. Table 3 gives a comprehensive summary of data predicted for regions of both subcooled nucleate boiling As the wall temperature increased further to the point of ONFDSUB, the degree of subcooling, bulk liquid temperature, and wall superheats are 21.29 K, 510.64 K, and 11.21 K respectively. The heat transfer until incipient boiling is a single-phase flow heat transfer of the liquid. Thus, the heat transfer coefficient is almost a constant along the channel, but when the wall temperature gradually increases, there is a corresponding increase in the liquid temperature near the wall thereby causing a transition from one regime to another. In order to control a smooth transition from continuous liquid bubbly flow to a continuous vapour droplet, it was assumed that the interfacial surface topology contains a multi-connected interface and the transition from one flow regime to another controlled by the vapour volume fraction. There exists a point (PDSAT) where the bulk liquid enthalpy equals the saturation enthalpy, and then boiling is said to occur. After boiling starts, the heat transfer coefficient begins to increase and then become almost constant after the saturated boiling initiation. Hence, the interfacial length from the inlet to this point is about 0.9 m, while the degree of subcooling, the bulk liquid temperature, and the degree of the wall superheat are approximately, 16.68 K, 531.94 K, and 21.29 K respectively. This situation continues until the fully developed saturated nucleate boiling by transition from nucleate boiling to forced convection evaporation heat transfer.
In summary, within the bubbly flow region, the vapour phase is dispersed in the continuous liquid in the form of bubble and the vapour volume fraction is said to <=0.28 in this region. While in the intermediate regime, the vapour   Figure 13, shows the vapour volume fractions (for different heat fluxes) at the onset of the partially developed saturated boiling (VVF-ONPDSAT), fully developed saturated boiling (VVF-ONFDSAT) and vapour volume fraction at the pipe exit (VVF-EXIT). This scenario can be observed when the bulk liquid enthalpy is higher than the saturation enthalpy. The amount of heat transfer in the saturated boiling region controlled by the mass quality or the vapour volume fraction. As reported in Table 3, the vapour volume fraction is about 28% at the onset of the saturated boiling, and as the heat flux increased along the pipe, the vapour volume fraction increases up to 46%. Further increase in the surface heat flux could result to burnout scenario. Burnout is characterised by a noticeable increase in surface temperature in response to a change in the wall heat flux. This change in temperature can be very drastic, as in the case of a Pressurized Water Reactor (PWR) or the change can The condition for critical heat flux approached for a given system by varying any of the independent variables. The four independent variables that influence critical heat flux, as investigated in this study are; the mass velocity, the channel inner diameter, the channel length, and the degree of inlet subcooling. For uniformly heated channels where the flow is assumed to be stable, the onset of critical heat flux condition occurs at the exit end of the channel. At a high mass velocity, the critical heat flux increases as the steam mass fraction increases along the flow channel, whereas, at a low mass flux, the steam mass fraction approaches a critical point due to a drastic decrease in the heat transfer coefficient, thereby causing a high localised heating effect at the channel exit. For predicting the critical heat flux, it is assumed that the fluid fed to the channel is subcooled and contains no entrained vapour, the flow is stableinand oscillations do not occur. Figure 14 shows the plot of the wall superheats at different aspect ratios for the onset of partially developed subcooled nucleate boiling and the onset of fully developed subcooled nucleate boiling, and the pressure drop (DP) against enthalpy rise for different heat flux and system pressure (45 bar) and mass flux (900 kg/m2 s) respectively. This chat is a lookup chat (with minimal error) for design modelling, and the results agree with any experimental data. Figure 15 shows how the wall heat flux partitioned along the heated channel, as explained earlier in section 3.2. The total heat flux balance can be given by splitting the heat flux into three different wall regions; the evaporation heat flux, quenching heat flux, and convective heat flux. All the heat balance accurately predicted and compared with the work of Adrian et al., as tabulated in Table 4. Figure 15 shows a smooth spontaneous exponential decay pattern of the convective heat flux up to a channel length of 1 m where the bubble nucleation is dominant. In addition, both the quenching and evaporation heat flux values show a better trend, which is a reflection of an accurate prediction as compared to the work of Adrian et al.
Most importantly, the results reported by Adrian is contradictory with the present work. Here, their results do not reflect the theory of the physics of the flow. According to Collier [2], the effect of both evaporation heat flux and the quenching heat flux can only be experienced when the bulk liquid gain enough thermal energy from the heated wall to cause bubble nucleation. In essence, both the quenching and evaporation heat flux values should be zero at the channel inlet. The report of Adrian et al., contradict this norm as it records the values of 63.48 W/m2 and 1745.755 W/m2 for hese two types of heat fluxes, respectively. Figure 16 shows the effect of aspect ratio on the singlephase entry length and the maximum void fraction recorded along the axial length. The ratio of the length of the channel to the channel diameter (L/D) investigated for three different cases using a 2D hexahedral refined mesh with grid points of 13,000. The first case which represents the diameter (D: 0.00385 m) with an aspect ratio of 519. 48 Figure 17 and Fig. 18 show the transient behaviour of bubble growth and detachment for flow boiling along the vertical pipe using the Volume of Fluid (VOF) multiphase submodels. As shown in Fig. 17 (Step 2), the fluid near the wall had gained enough thermal energy to initiate bubble nucleation at the wall cavities. Hence, forces responsible for bubble growth and bubble detachment plays a vital role in this stage. The formation, growth, and detachment of the vapour bubble depend on the influence of the forces acting on the bubble [51]. The surface tension and inertia forces are responsible for the bubble formation in the highly subcooled region and also holds the bubble to the heated surface, while the buoyancy and the frictional (drag) force are responsible for the bubble detachment from the heated surface in the low subcooled region [2]. For a spherical bubble to remain in the bulk fluid after detachment from the heated surface, the pressure in it must be higher than the pressure of the surrounding fluid, and its size must be larger than the critical size of a bubble. This phenomenon is practically observed by the bubble contours in the current case study, as seen in Fig. 17 and Fig. 18.

Conclusions
This study considers the 2D and 3D numerical approachs for simulating a boiling phenomenon for upward vertical flow boiling in a pipe using the Modified Eulerian Model and Wall Heat Balance Model (WHBM). The model and sub-models used in this study were able to provide more accurate results than the few reported in the open literature. Parameters investigated were the axial distribution of bulk liquid temperature, axial distribution of the wall adjacent temperature, axial bulk average vapour volume fraction, superheat along the wall surface, the bulk liquid subcooling with axial length, tracking of the single bubble departure diameter, heat balance at the wall, the effect of aspect ratio, the effect of inlet subcooling and the critical heat flux. All the results obtained are in good agreement with the experimental data used for the validation and agree with the physics of the theory. The two critical features established in this study were the effect of grid points and turbulence on the non-boiling length and the output void fraction. The results obtained for different test cases show how the number of grid points and the choice of turbulence model could influence the accuracy in resolving the different flow fields. Runs1-Runs4, as shown in Table 2 and Fig. 3 have different numbers of grid points, and the results obtained show reasonable percentage deviation of the entry length and void fraction, using either the k-epsilon and k-omega models. Longer entry length is experienced for low mesh cells, whereas, high mesh cells results in a shorter entry length, thereby enhancing effective/efficient heat transfer out of the channel.
Other key observations deduced from this study can be summarised as below; & The Modified Eulerian Two-phase Coupled Model is more robust and can predict more accurately flow boiling parameters in an upward flow through a heated vertical channel. & The amount of superheat necessary for nucleate boiling to occur is dependent on the inner surface heat flux. Therefore, no boiling can occur when the inner surface temperature is less than the saturation temperature at a given location along the channel. & The adjacent wall temperature close to the inlet is the same as the fluid inlet temperature because the heating process of the wall undergoes a gradual and sustained thermal conduction effect. As such, one should not expect an abrupt increase in heat flux distribution along the wall. The Onset of Sub-cooled Nucleate Boiling initiates when the temperature of the wall at the single-phase region equals the sum of the saturation temperature at the given system pressure, and the wall superheats necessary to cause nucleation. & The onset of the critical heat flux condition occurs at the exit end of the channel In conclusion, this current study employed alternative numerical phase coupling models as compared to ones used by previous researchers. The results obtained in this study is comprehensive, explicit, and more accurate than those reported in the literature. It is thus, recommended as a validation baseline for future study.