Numerical Study on Pool Boiling of Hybrid Nanofluids Using RPI Model

Abstract

The performance of deionized (DI) water and hybrid nanofluids for pool boiling from a horizontal copper heater under atmospheric pressure conditions is numerically examined in the current study. The Eulerian–Eulerian scheme is adopted with a Rensselaer Polytechnic Institute (RPI) sub-boiling model to simulate the boiling phenomena and predict the heat and mass transfer in the interior of the pool boiling vessel. This paper attempts to correct the coefficient of the bubble waiting time (BWTC) in the quenching heat flux partition as a proportion of the total heat flux and then correlate this coefficient to the superheat temperature. The pool boiling curve and pool boiling heat transfer coefficient (PBHTC) obtained for the present model are verified against experimental data from the literature and show good agreement. In addition, this work comprehensively discusses the transient analysis of the vapor volume fraction contours, the vapor velocity vectors, and the streamlines of water velocity at different superheat temperatures. Finally, for BWTC, new proposed correlations with high coefficients of determination of 0.999, 0.932, and 0.923 are introduced for DI water and 0.05 vol.% and 0.1 vol.% hybrid nanofluids, respectively.

1. Introduction

The continuous demand for efficient heat exchange systems with reduced size and high heat flux removal has encouraged researchers to propose and develop new techniques. One of the most essential passive techniques is the use of novel fluids with some efficient thermal properties such as thermal conductivity, viscosity, density, heat capacity, and surface tension [1,2,3,4]. On the other hand, using some surface modification techniques by the means of micro- or/nanostructured geometries is considered another proven method to increase the ability of heat removal from thermal systems [5]. The idea of using new thermofluids with superior thermophysical properties by combining solid nanomaterial with a liquid was firstly introduced by Choi and his team [6], who called them nanofluids, which are colloidal suspensions of nanoscale materials within conventional fluids to enhance the thermal conductivity of those fluids. In recent years, nanofluids have been used as a new class of thermal fluids for cooling applications by applying two-phase flow systems [7,8,9].

As an efficient mode, the boiling heat transfer process could dissipate a lot of heat via latent heat of vaporization, especially during the nucleate boiling regime. Utilizing this phenomenon together with nanofluids may enhance the performance of heat exchange systems and, thus, increase efficiency and save energy. Recently, there has been a clear development by researchers of thermofluids through applying new fluids, so-called hybrid nanofluids, which were achieved by mixing two or more nanomaterials based on conventional liquids [10,11]. Hybrid nanofluids are a new form of mono nanofluids. They have shown some good thermal properties, especially thermal conductivity. This gives us an excellent insight to use them rather than single nanofluids in heattransfer applications as stated in previous studies [3,10,12]. In the present work, the nucleate PBHT process from a horizontal, typical copper heater for DI water and a new type of hybrid nanofluid with a mixing ratio of 50:50 was predicted using an Eulerian–Eulerian model. The coefficient of a bubble waiting time in the quenching heat flux section was adjusted and then correlated to the superheat temperature for purposes of validation for both DI water and hybrid nanofluids.

In the next section, the authors present the relevant studies that deal with the boiling of nanofluids, employing the numerical works that deal with the pool boiling of nanofluids.

2. Literature Overview

In recent years, many studies have been reported that deal with experimental pool boiling using different passive and active improvement methods such as surface modification or addition of solid nanoscale materials [13,14,15,16,17,18,19]. On the other side, the knowledge of the PBHT modeling utilizing nanofluids is still weak. Mortezazadeh et al. [20] studied the pool boiling of ferrofluids with the presence of a moving electric charge within which the force of magnetism acts using the Eulerian–Eulerian scheme. They adopted the heat flux partitioning sub-boiling model to predict the boiling performance, and their results depicted that using this type of nanofluids resulted in degradation in the BPHTC and vapor volume fraction. In addition, there was an enhancement in BPHTC when using the magnetic field. Aminfar et al. [21] also studied the PBHT performance of mono nanofluids on flat surfaces using a multiphase model (mixture model). Their results showed that the two-phase model is more accurate than the three-phase model. Niknam et al. [22] studied the effect of particle size on pool boiling behavior using the RPI sub-boiling model. They introduced a nucleation site density ratio, and a correlation was derived based on experimental pool boiling data from the literature to include this correlation in the RPI model.

Mahdavi et al. [23] included the influence of the bubble departure diameter and nucleation site density on PBHT of nanofluids from a set of two horizontal heaters immersed in a pool boiling chamber. They implanted the closure correlations related to the nucleation parameters’ as user-defined functions in their Eulerian–Lagrangian approach by adopting a discrete model. Kamel et al. [24,25] introduced new correlations related to the bubble waiting time coefficient by correcting this coefficient and correlate it to the superheat temperature to modify the quenching heat flux as an important item inside the RPI model when using mono nanofluids. Mousavi et al. [26] studied the effect of surface roughness on the PBHT behavior represented by PBHTC of water and mono nanofluids using the RPI model. They used various nanofluid concentrations with different applied heat fluxes at atmospheric pressure. Their results demonstrated that when the surface roughness increased, the PBHTC was enhanced for nanofluids.

Zaboli et al. [27] also used the Eulerian–Eulerian approach to simulate the PBHT of different concentrations of SiO₂ nanofluids using the RPI model on a flat plate heater. A single-phase mixture model was adopted to predict this phenomenon. They proposed new correlations for PBHTC, bubble departure diameter, and nucleation site density based on their numerical data. The numerical results were improved when they considered the modification that happen using nanofluids, and the PBHTC was enhanced for nanofluids with 0.1 vol.%, especially in high heat flux region. Majdi et al. [28] investigated two types of nanomaterials (Al₂O₃ and CuO) based on water mono nanofluids for their pool boiling performance. They used various shapes of fins to increase the surface area of the boiling heater. The results indicated that, with the increase in the concentration of nanomaterials, the velocity of the vapor increased and affected the pressure and velocity of flow during the boiling process. Mao et al. [29] also numerically studied the PBHT of R134a from a horizontal heater under atmospheric pressure conditions. The Eulerian–Eulerian multiphase model was adopted with an RPI nucleate sub-boiling model. They used various heat fluxes and different locations for the heater inside the pool boiling chamber. The results showed good agreement with the experimental data when using the modified correlation of nucleation site density of the pool boiling heat transfer. According to the simulation of azimuthal variations of temperature around the tube, it was found that the region between 90 and 130 degrees has the highest heat transfer behavior. The maximum temperature happens at the top of the heater due to the velocity and turbulence intensity in this region being low. Their study helps to explain the phenomena and understand the temperature distributions in the experiment.

3. The Aim of This Work

This work aims to predict the pool boiling of (Al₂O₃ + MgO) hybrid nanofluids from a typical horizontal heater inserted inside a boiling chamber at atmospheric pressure conditions. Various hybrid nanoparticle concentrations with a mixing ratio (50:50) based on DI water hybrid nanofluids and heat fluxes are used to simulate this phenomenon. New correlations for the bubble waiting time coefficient (BWTC) are proposed to modify the quenching heat flux section that is included inside the total heat flux of the wall boiling submodel.

4. Theoretical Formulation and Numerical Simulation

4.1. Governing Equations

As shown in previous works, the conservation equations (Equations (1)–(3)) control the overall balance of mass, momentum, and energy at each phase [30,31,32,33].

$$\mathrm{Continuity} \mathrm{equation}: \frac{\partial \left({\rho }_h{\alpha }_h\right)}{\partial t}+𝛻.\left({\rho }_h{\alpha }_h{\overrightarrow{v}}_h\right)={\dot{m}}_{hp}-{\dot{m}}_{ph}$$

(1)

$$\mathrm{Momentum} \mathrm{equation}: \frac{\partial \left({\rho }_h{\alpha }_h{v}_h\right)}{\partial t}+𝛻.\left({\rho }_h{\alpha }_h{\overrightarrow{v}}_h{\overrightarrow{v}}_h\right)=-{\alpha }_h𝛻P+{\rho }_h{\alpha }_h\overrightarrow{g}+𝛻{\tau }_h+\left({\dot{m}}_{hp}{\overrightarrow{v}}_h-{\dot{m}}_{ph}{\overrightarrow{v}}_h\right)+S_{1h}$$

(2)

$$\mathrm{Energy} \mathrm{equation}: \frac{\partial \left({\rho }_h{\alpha }_h{E}_h\right)}{\partial t}+𝛻.\left({\rho }_h{\alpha }_h{\overrightarrow{v}}_h{E}_h\right)={\alpha }_h\frac{\partial P}{\partial t}+𝛻{\dot{q}}_h+q_{exchange,ph}+\left({\dot{m}}_{hp}{E}_h-{\dot{m}}_{ph}{E}_p\right)+S_{2h}$$

(3)

where the subscript h denotes $h$-th phase (h =$l$ for the water or nanofluid phase and h = $g$ for the vapor phase). $P$ refers to the pressure. ${\rho }_h, {\alpha }_h,{\overrightarrow{v}}_h$ are the density, friction of the volume, and velocity of the $h$-th phase, respectively. ${\dot{m}}_{hp}$ is mass transfer across the contacting surfaces’ interface; in the liquid phase (water), this term is equal to zero. This is attributable to the fact that the boiling process begins at the saturation temperature. The term $\left({\dot{m}}_{hp}{\overrightarrow{v}}_h+{\dot{m}}_{ph}{\overrightarrow{v}}_h\right)$ in the equation of momentum explains that water evaporation or vapor condensation causes a momentum transfer. Furthermore, the term $\left({\dot{m}}_{hp}{E}_h+{\dot{m}}_{ph}{E}_p\right)$ in the equation of energy indicates the energy transfer that arises as a result of phase change. Additionally, the terms $\overrightarrow{g}$, ${\dot{q}}_h$, $S_{1h}$, $S_{2h},$ and $q_{exchange,ph}$ are the gravitational acceleration, heat flux, the interphase momentum transfer term, the interphase energy transfer term, and the direct heat transfer in phase “$h$”, which can be calculated as follows:

$$q_{exchange,ph}=HT{C}_{if}{A}_{if}\left(T_p-T_h\right)$$

(4)

where $HT{C}_{if}$, $A_{if},$ and $T$ are the interface heat transfer coefficient, interfacial area, and fluid temperature, respectively. In this work, to find out the solution to the mathematical model that formulates the exchange of phase interactions linked to interfacial momentum, heat and mass transfer were introduced. There are two types of flow in the viscous model, i.e., the laminar and the turbulent regimes, and because of the vapor bubbles’ nature and their dynamics, it is considered a chaotic region. Moreover, the rough model with two well-known, realizable equations, the $k-epsilon$ model, was adopted as it is recommended by previous studies [24,30,32] due to its appropriate performance for boiling heat transfer. Equations (5)–(7) show the formulation of the turbulence model used in this study.

$$\frac{\partial \left({\alpha }_h{\rho }_h{K}_h\right)}{\partial t}+𝛻.\left({\alpha }_h{\rho }_h{K}_h{\overrightarrow{v}}_h\right)=𝛻.\left({\alpha }_h\frac{{\mu }_{T,h}}{{\sigma }_K}\right)+{\alpha }_h{\tau }_{T,h}𝛻v_h-{\alpha }_h{\rho }_h{\epsilon }_h+S_{K,h}^{int}$$

(5)

$$\frac{\partial \left({\alpha }_h{\rho }_h{\epsilon }_h\right)}{\partial t}+𝛻.\left({\alpha }_h{\rho }_h{\epsilon }_h{\overrightarrow{v}}_h\right)=𝛻.\left({\alpha }_h\frac{{\mu }_{T,h}}{{\sigma }_h}\right)+{\alpha }_h\frac{{\epsilon }_h}{K_h}(C_{\epsilon 1}{\tau }_{T,h}𝛻v_h-C_{\epsilon 2}{\rho }_h{\epsilon }_h)+S_{\epsilon ,h}^{int}$$

(6)

$${\tau }_{T,h}={\mu }_{T,h}\left[𝛻{\overrightarrow{v}}_h+\left(𝛻{\overrightarrow{v}}_h\right)-\frac{2}{3}\left(𝛻.{\overrightarrow{v}}_h\right)I\right]-\frac{2}{3}{\rho }_h{K}_{h}I$$

(7)

In Equations (5) and (6), $K$ and $\epsilon$ denote the turbulence kinetic energy and dissipation rate of the turbulent model, respectively. ${\mu }_{T,h}=C_{\mu }{K}^2/\epsilon$ is the turbulent viscosity coefficient, $S_{K,h}^{int}$, $S_{\epsilon ,h}^{int}$ are the source terms of the interaction of the nucleate bubbles and the turbulent fluid regime for water. The constants in Equations (5) and (6) are set similar to the ones presented in previous studies [30,32], which were set as follows: $C_{\mu }=0.09, {\sigma }_K=1.0, {\sigma }_{\epsilon }=1.3, C_{\epsilon 1}=1.44, C_{\epsilon 2}=1.92$.

4.2. Phase Interaction (Interfacial Exchange)

The interaction between the two-phase flow (liquid–vapor) could lead to creating some important forces such as drag force, viscous force, lift force, virtual mass force, wall lubrication force as well as turbulent dispersion force, especially the nature of the boiling including vapor bubbles, which, in turn, make this flow regime typically turbulent due to the chaos of the bubbles inside the water phase. The term $S_{1k}$ is the interphase momentum transfer term, which is written on the left side of the momentum-governing equation (Equation (2)) refers to all the above interaction forces [31]. Equation (8) shows the forces during the interphase interaction.

$$S_{1h}=F_{h, drag}+F_{h, lift}+F_{h, VM}+F_{h, dispersion}+F_{h, L wall}$$

(8)

In the current simulation, detailed information on all the interfacial forces is presented in

Table 1. Models for interfacial forces implemented in this simulation.

Forces	Model [Ref.]	Formulation
Virtual mass forces	Explicit source term [35]	$F_{h, VM}=C_{VM}{\alpha }_{vapor}{\rho }_{water}\left(\frac{d_{water}{\overrightarrow{v}}_{water}}{dt}-\frac{d_{vapor}{\overrightarrow{v}}_{vapor}}{dt}\right)$ $C_{VM}$ is the virtual mass constant, which is equal by default (0.5), and the term $\frac{d_{water}}{dt}$ denotes the liquid-phase time derivative.
Drag force	Ishii [36]	$C_D=Min\left(C_{D, viscous}, C_{D,distorted}\right)$ $C_D$ is the drag coefficient and is achieved by choosing the minimum of $C_{D, viscous}$ and $C_{D,distorted},$ which are the coefficients of the viscous and distorted regimes, respectively. $C_{D, viscous}=\frac{24}{R{e}_{vapor}}\left(1+0.15R{e}_{vapor}{}^{0.75}\right)$ $C_{D,distorted}=\frac{2}{3} \frac{d_{bubble}}{\sqrt{\frac{\sigma }{G\left({\rho }_{water}-{\rho }_{vapor}\right)}}}$
Turbulent dispersion force	Lopez-de-bertodano [37]	$F_{water, dispersion}=-F_{vapor, dispersion}=C_{TD}{\rho }_l{K}_{water}𝛻{\alpha }_{vapor}$ $C_{TD}$ is a user-modifiable constant that is set to 1, $K_{water}$ is the turbulent kinetic energy in the water phase. $𝛻{\alpha }_{vapor}$ is the gradient of the vapor-phase volume fraction.
Lift force	Tomiyama [38]	$F_{h, lift}=-C_L{\alpha }_{vapor}{\rho }_{water}\left({\overrightarrow{v}}_{water}-{\overrightarrow{v}}_{vapor}\right)\left(𝛻{\overrightarrow{v}}_{water}\right)$ $C_L$ is the lift coefficient $C_l=\left\{\begin{array}{c}min\left[0.288tanh\left(0.121R{e}_{vapor}\right),f\left(E{\acute{o}}_{modif}\right)\right] E{\acute{o}}_{modif}\le 4\\ f\left(E{\acute{o}}_{modif}\right) 4<E{\acute{o}}_{modif}\le 10\\ -0.27 10<E{\acute{o}}_{modif}\end{array}\right.$ $f\left(E{\acute{o}}_{modif}\right)=0.00105E{\acute{o}}_{modif}{}^3-0.0159E{\acute{o}}_{modif}{}^2-0.020E{\acute{o}}_{modif}+0.474$ $Eo$ is a modified Eotvos number that is defined as follows: $E{\acute{o}}_{modif}=\frac{g\left({\rho }_{vapor}-{\rho }_{water}\right)d_h^2}{\sigma }$ $d_h=d_{bubble} {\left(1+0.163E{o}^{0.757}\right)}^{\frac{1}{3}}$ $d_h,$ $d_{bubble}$ are deformable bubbles and bubble diameter. $Eo=\frac{g\left({\rho }_{vapor}-{\rho }_{water}\right)d_{b.d}^2}{\sigma }$
Wall lubrication force	Antal et al. [39]	$F_{h, L wall}=C_{L wall}{\rho }_{water}{\alpha }_{vapor}{\left|\left({\overrightarrow{v}}_{water}-{\overrightarrow{v}}_{vapor}\right)\right|}^2{\overrightarrow{n}}_W$ $C_{L wall}$$, {\overrightarrow{n}}_W$ are the wall lubrication coefficient and the normal unit pointing away from the wall. $C_{L wall}=max\left(0, \frac{C_{W1}}{d_{bubble}}+\frac{C_{W2}}{y_W}\right)$; $C_{W1}=-0.01$, $C_{W2}=0.05$ are non-dimensional coefficients.

. The heat transfer from the liquid to the vapor phase was also included in this work, and due to the non-thermal equilibrium across the interphase, this transmitted energy might develop [31]. To calculate the interfacial HTC between the liquid and gas phases $\left(HT{C}_{if}=k_{l}N{u}_g/d_{bubble}\right)$, the Ranz–Marshall model [33] was adopted in this simulation in order to simulate heat transfer between phase interactions, Equation (9) introduces this model.

$$N{u}_{vapor}=2.0+0.6P{r}_{water }^{0.333}R{e}_{vapor}^{0.5}$$

(9)

where $R{e}_{vapor}$ is the relative Reynolds number based on the diameter of the bubble, relative velocity, and $P{r}_{water}$ is the Prandtl number of the water phase, as shown in Equation (10).

$$P{r}_l=\frac{C_{p,water}{\mu }_{water}}{k_{water}}$$

(10)

where $C_{p,water}$, ${\mu }_{water}$, and $k_{water}$ are the water-phase specific heat, viscosity, and thermal conductivity, respectively. For turbulent interaction (mixture turbulence model), the Troshko–Hassan [34] model is used as shown in Equations (11)–(13).

$$S_{K,h}^{int}=C_{ke}{\displaystyle \sum }_{p=1}^M{K}_{vaporl water}{\left|{\overrightarrow{v}}_{water}-{\overrightarrow{v}}_{vapor}\right|}^2$$

(11)

$$S_{\epsilon ,h}^{int}=C_{td}\frac{1}{{\tau }_p}{S}_{K,h}^{int}$$

(12)

where $C_{ke}=0.75$ and $C_{td}=0.45$,${\tau }_p$ is the characteristic time of the induced turbulence, which is defined as follows:

$${\tau }_p=\frac{2C_{VM}{d}_{bubble}}{3C_D\left|{\overrightarrow{v}}_{water}-{\overrightarrow{v}}_{vapor}\right|}$$

(13)

4.3. Heat Flux Partitioning Boiling Submodel

In this simulation, the Rensselaer Polytechnic Institute (RPI) boiling sub-model was utilized to predict the heat and mass transfer from the heating surface through the nucleate pool boiling of pure water and hybrid nanofluid two-phase flow [40,41]. In the RPI classical boiling submodel, the total heat flux from the heater wall consists of three heat flux mechanisms: first, the heat flux caused by the energy transfer from the generation of bubbles represented by natural convection to the bulk water, and this called ${\dot{q}}_{convection};$ second, the heat flux, caused by the latent heat of vaporization ${\dot{q}}_{evaporation}$; and last, the heat flux caused by quenching process ${\dot{q}}_{quenching},$ which is the repeated averaged transient energy transfer process associated with water filling the wall region following bubble separation from the surface. Equations (14)–(17) present the three terms of heat fluxes that used in this boiling submodel.

$${\dot{q}}_{total}={\dot{q}}_{convection}+{\dot{q}}_{evaporation}+{\dot{q}}_{quenching}$$

(14)

$${\dot{q}}_{convection}=A_{convection}\times h_{convection}\times \left(T_{wall}-T_{water}\right)$$

(15)

$${\dot{q}}_{evaporation}=V_{bubble}{N}_a{\rho }_{vapor}{h}_{lg}f=\frac{\pi }{6}{d}_{bd}^3{\rho }_{vapor}f{N}_a{h}_{lg}$$

(16)

$${\dot{q}}_{quenching}=\mathrm{B}\mathrm{W}\mathrm{T}\mathrm{C}\frac{2k_{water}}{\sqrt{\pi {\alpha }_l{t}_{waiting}}}\left(T_{wall}-T_{water}\right)A_{bubble}$$

(17)

where ${\dot{q}}_{convection}$, ${\dot{q}}_{evaporation},$ and ${\dot{q}}_{quenching}$ indicate the elements of the total heat flux induced by natural convection, evaporation, and quenching, respectively. Furthermore, $d_{bubble}$, $f$, $N_a$, $A_{convection}$, $A_{quenching,}$ and $t_{waiting}$ are the diameter of the bubble, the frequency of the bubble, the sites of active nucleation, the convection and quenching area fractions of the heated wall, and the time being spent waiting for a bubble, respectively. $C_w$ is the bubble waiting time coefficient, BWTC; the default value of this coefficient in Ansys Fluent solver was assumed to be 1, and this value can be modifered according to the user’s needs.

Equations (14)–(17) represent the heat flux partitioning boiling submodel required for a closure correlation to predict the nucleation of boiling parameters. There is a wide range of correlations presented in the literature to predict the boiling parameters of pure liquids that are related to bubble dynamics. According to the authors’ best knowledge, the parameters of nucleate boing, especially the nucleation sites’ density, the bubbles’ departure diameters, and the frequency of bubbles could be affected during the boiling of nanofluids since nanoparticles are deposited on the heated surface during this process. As a result, efforts must be made to account for the heater wall alteration to reveal new closure correlations associated with the characteristics of bubble dynamics through the boiling process. In the current simulation, the nucleation site density of water, which depends on the superheat temperature, was predicted using the Lemmert and Chawla model [42] as shown in Equation (18).

$$N_a=C^n{\left(T_{wall}-T_{sat}\right)}^n$$

(18)

where $n=1.805 \mathrm{and} C=210$ are the experimental variables that this model relies on. It can be revealed from the literature that the deposition of nanoparticles during the boiling process on the heater wall affects the sites of nucleation of the bubble by altering the surface structure, thereby changing the active nucleation sites as well as improving the wettability of the heater surface. Ganapathy and Sajith [43] proposed a semi-analytical model for pool boiling of nanofluids by considering the change of the nucleation sites’ density. In their model, the influence of the deposition of nanomaterials on nucleation sites was included. Moreover, the nanoparticle size and the wettability augmentation parameter related to the surface roughness were adopted, as shown in Equation (19).

$$N_a=218.8 \frac{1}{\gamma }P{r}_l^{1.63}{\left\{\left[14.5-4.5\left(\frac{P{R}_a}{\sigma }\right)+0.4{\left(\frac{P{R}_a}{\sigma }\right)}^2\right] {\beta }^{-3}{\left(\frac{R_a}{d_p}\right)}^{-0.5}\right\}}^{-0.4}∆T_{sup}$$

(19)

where $P$, $d_p,$ and $R_a$ represent the pressure, nanoparticle diameter, and average surface roughness, respectively. Furthermore, $\gamma$ is a surface–liquid interaction criterion that the surface and liquid materials describe, and $\beta$ is a parameter that enhances surface wettability as described in Equation (20).

$$\beta =\frac{1-cos\theta }{1-cos{\theta }^{*}}$$

(20)

where $\theta \mathrm{and} {\theta }^{*}$ are the contact angles on the nanoporous and clean surfaces, respectively.

Kocamustafaogullari and Ishii [44] proposed a correlation to predict the bubble departure diameter during the boiling heat transfer, and due to the equilibrium between the surface tension and the buoyancy force on the heater wall as stated in [26], this correlation was used in the current simulation as shown in Equation (21).

$$d_{bd}=0.0012{\rho }_{water}^{-0.9}\times 0.0208\theta \sqrt{\frac{\sigma }{\left({\rho }_{water}-{\rho }_{vapor}\right)g}}$$

(21)

The frequency of bubble departure is linked to the bubble departure diameter which has been widely proved to have a direct effect on this factor, and it was shown to decrease as the diameter of bubble departure increases, for both pure liquids and nanofluids. Equation (22) shows the main relationship between the frequency and diameter of bubbles in the boiling process. In this work, the equation introduced by Cole was used as stated in previous studies [26,30,45] based on Equation (22). The index $n$ has a different values in the majority of common correlations, for example, for Cole correlationit has the value 0.5 as in Equation (23).

$$f~\frac{1}{d_{bd}^n}$$

(22)

$$f=\sqrt{\frac{4g\left({\rho }_{water}-{\rho }_{vapor}\right)}{3d_{bd}{\rho }_{water}}}$$

(23)

4.4. Thermophysical Properties

The thermophysical properties of working fluids and nanoparticles are listed in

Table 2. Thermophysical properties of working fluids and nanoparticles used in the current simulation.

Property	DI Water [46]	Vapor [46]	Al2O3 [26]	MgO [47]	Hybrid (50:50)
Density $\left[{\mathrm{kg} \mathrm{m}}^{-3}\right]$	958.35	0.59817	3490	3580	3535
Specific heat $\left[{\mathrm{J} \mathrm{kg}}^{-1}{ \mathrm{K}}^{-1}\right]$	4215.7	2080	451	874	662.5
Thermal conductivity $\left[{\mathrm{W} \mathrm{m}}^{-1}{ \mathrm{K}}^{-1}\right]$	0.67909	0.02509	25	55	40
Dynamic viscosity $\left[\mathrm{Pa}·\mathrm{s}\right]$	0.000281	0.0000122	-	-	-
Surface tension $\left[{\mathrm{N} \mathrm{m}}^{-1}\right]$	0.0589	-	-	-	-

. In this study, due to the dilute concentrations of hybrid nanofluid (less than 1 vol.%), the nanofluids are assumed to be in a single-phase flow (homogeneous flow) as stated in previous works [26,30]. The mixing ratio of hybrid nanoparticles was 50:50, so the averaging of thermophysical properties for hybrid nanoparticles was used in this simulation. In addition, the properties of water and vapor included in this simulation were those at a saturation temperature of 100 °C according to NIST Chemistry WebBook [46]. Moreover, the surface tension of nanofluid was supposed to be the same as that of water due to the diluted volumetric concentration. The following correlations were utilized to calculate the properties of thermophysical quantities of effective hybrid nanofluids in this study with acceptable accuracy and cover the range of conditions for the concentration, shape, and nanoparticle types used in the current work. Equations (24)–(27) show the correlations of thermophysical properties, which were used also in previous studies [30,47,48].

$${\rho }_{HNFs}=\phi {\rho }_{HNPs}+\left(1-\phi \right){\rho }_{water}$$

(24)

$$C_{p,HNFs}=\frac{\phi {\left(\rho C_p\right)}_{HNPs}+\left(1-\phi \right){\left(\rho C_p\right)}_{water}}{\phi {\rho }_{HNPs}+\left(1-\phi \right){\rho }_{water}}$$

(25)

$$K_{HNFs}=K_{water}\left[\frac{K_{HNPs}+2K_{water}-2\phi \left(K_{water}-K_{HNPs}\right)}{K_{HNPs}+2K_{water}+\phi \left(K_{water}-K_{HNPs}\right)}\right]$$

(26)

$${\mu }_{HNFs}=\left(1+7.3\phi +123{\phi }^2\right){\mu }_{water}$$

(27)

4.5. Domain Description

To investigate the pool boiling heat transfer behavior of alumina/magnesium oxide hybrid nanofluids at atmospheric pressure conditions, a 2-D square pool boiling chamber was selected with a whole domain of 120 mm length and 120 mm width, containing a typical copper tube with an outer diameter of 22 mm. To make a validation of the recent simulation results, our previous experimental study [31] was selected. A schema diagram of the current domain is provided in Figure 1.

4.6. Assumptions

In this simulation, an appropriate strategy is implemented to obtain an approximate solution due to the high density ratio of the water and vapor phases. A phase-coupled SIMPLE scheme (PC-SIMPLE) was selected to handle the pressure–velocity coupling and give stability in predicting the pool boiling process [26,32]. The following assumptions and simplifications are considered in this simulation.

• The problem under consideration is the transient and turbulent flow due to bubble formation in the nucleate pool boiling regime.
• As a result of molecular mixing of the low concentration of hybrid nanoparticles and water, the hydrodynamic behavior of the hybrid nanofluids would be similar to that of the single-phase nanofluid. Therefore, a single-phase model is considered in this study.
• Under the specified operating temperature and pressure, the thermophysical characteristics of water and vapor phases are assumed to be constant.
• The operating pressure of the chamber is controlled to be the same as the atmospheric pressure condition, then the pressure is $P=101.325 \mathrm{kPa}$.
• Due to the low volume fraction employed in this analysis, the surface tension parameter of hybrid nanofluids is presumed to be the same as that of water.
• The temperature of the DI water inside the chamber is the saturation temperature.
• Due to the high density ratio between the water and vapor phases, the vapor phase is also believed to be quite light to take the nanoparticles within. Therefore, it is considered that the stable nanoparticles do not affect the thermal properties of the vapor phase.
• For this simulation, a time interval size of (0.001 s) is used. Moreover, following a trial-and-error strategy to ensure that the solution is convergence at each time step, the maximum number of iterations per time step was adjusted to 100.

4.7. Numerical Methods and Boundary Conditions

In this work, the numerical simulation of pool boiling of hybrid nanofluids from the horizontal copper heater is performed for purpose of validating and predicting pool boiling performance. The geometry of this simulation was created according to previous experimental work [31], and it was built as a 2-D chamber by using Ansys Design Modeler Toolbox. The Finite Volume Method (FVM) was utilized to conserve the computational model spatially and to transform the partial differential equations into linear algebraic equations. Hence, Ansys Fluent Software was used to solve the conservation equations for mass, momentum, and energy for water and vapor phases using the FVM. As mentioned earlier, a phase-coupled SIMPLE algorithm was selected to treat the pressure–velocity coupling. In this method, the velocities are solved coupled by two phases in a segregated technique, whilst the pressure correction is handled based on total continuity. Therefore, the coefficients of the pressure correction were obtained from the coupled per phase momentum equations in this technique. In addition, a second-order upwind method was selected to converge the momentum, turbulent kinetic energy, and turbulent dissipation rate. Moreover, the first-order upwind method was used for energy equation, and the gradient of all flow variables was determined utilizing the least-square cell-based technique. The current simulation was continued until the scaled residuals reached the convergence setting criteria, which was ${10}^{-4}$. Regarding the boundary conditions of this model, it was subjected to the following conditions: A constant temperature was assumed to be at the heater surface. In addition, heat flux was assumed to be zero at the adiabatic walls of the boiling vessel. Moreover, the top of the boiling vessel was assumed to be at atmospheric pressure as shown in the following Equations (28)–(30).

$$T=T_{wall}$$

(28)

$${\dot{q}}_{total}=-K_{wall}\frac{\partial T}{\partial x}=0$$

(29)

$$P=P_{atmospheric}$$

(30)

4.8. Grid Test and Validation

The structural grids were created through the mesh tool, which is available in the Ansys software toolbox as shown in Figure 2. To test the grid dependency of this simulation, a different number of elements were checked (29,276, 46,256, and 68,345) to see the sensitivity of the results of two variables, which are the total heat flux and vapor volume friction with constant superheat temperature 6 K for water pool boiling process. The maximum relative errors between the three elements for heat flux and vapor volume friction were found to be 1.83% and 2.99%, respectively. Hence, mesh refinement did not result in any significant changes to the mean values of heat flux and the vapor volume fraction. Therefore, the number of moderate mesh elements (46,256) was considered in this simulation to obtain a balance between the simulation time and the accuracy of the obtained data. Figure 3 shows the grid dependency test for the present simulation.

The details of the present simulation related to the geometry structure, working fluids, boundary conditions, and the investigated parameters are shown in

Table 3. The details of the geometry structures, working fluids, boundary conditions, and investigated parameters for the present simulation.

Domain Structure	Working Fluids	Model/Submodel	Boundary Conditions (BCs)	Purpose of Simulation	Studied Parameters
2-D square chamber with a horizontal copper heater as in the experimental work of [31]	Water–vapor two-phase flow Alumina/magnesium oxide hybrid nanofluid–vapor two-phase flow	Eulerian–Eulerian multiphase model/RPI boiling submodel	$∆T=4.5–13.3 \mathrm{K}$ ${\dot{q}}_{total}=15–116 \mathrm{kW}/{\mathrm{m}}^2$ $\phi =0.05 \mathrm{vol}.\%$ $\phi =0.1 \mathrm{vol}.\%$ $R_a=0.382 \mathsf{\mu }\mathrm{m}$ ${\theta }_{HNFs}=38°$	Validation Pool boiling heat transfer performance prediction	Pool boiling curve PBHTC Heat flux partitioning portions V-V-V * V-V-F ** W.S.V ***

V-V-V * vapor volume friction, V-V-F ** vapor velocity vector, W.S.V *** water streamline velocity.

. The purpose of this work was the validation and determination of PBHT performance for both DI water and hybrid nanofluids. The PBHT of DI water and (Al₂O₃ + MgO) hybrid nanoparticles based on deionized water hybrid nanofluids from a typical horizontal copper heater at atmospheric conditions was numerically studied. The current simulation was validated with an extended heat flux partitioning RPI sub-boiling model for deionized water and hybrid nanofluids with volume concentrations (0.05 vol.% and 0.1 vol.%). To validate the accuracy of the present simulation results, the previous experimental results of pool boiling heat transfer for DI water and hybrid nanofluids [31] at atmospheric pressure conditions were compared to the current numerical data. Figure 4 and Figure 5 depict the data of the pool boiling curve and BPHTC for DI water and hybrid nanofluids for 0.05 vol.% and 0.1 vol.% volume concentrations. Comparison was undertaken between the two types of RPI models, i.e., the classical RPI model and the extended RPI model, for both the DI water and hybrid nanofluids. It was observed from both figures that the numerical results for DI water and hybrid nanofluids were in good agreement with the experimental results of [31] when using the properties of water and vapor phases together with modifying the bubble waiting time coefficient (BWTC) and correlating this coefficient to superheat temperature. Hence, this indicates that the extended RPI model using the modified BWTC to adjust the timing between the subsequent bubbles’ escape during the nucleate regime gives a good prediction for boiling of DI water as well the hybrid nanofluids in parallel with using the closure correlations of boiling parameters during boiling of hybrid nanofluids. In addition, the validity of the model was examined through the relative root mean square error (rRMSE). The value of rRMSE was estimated as 3.1%, 1.5%, and 2.5% for the cases of DI water, 0.05 vol.% hybrid nanofluid, and 0.1 vol.% hybrid nanofluid, respectively.

5. Results and Discussion

In this section, important parameters such as the proposed correlations of the bubble waiting time coefficient (BWTC), the heat flux portions of the RPI model, vapor volume fractions, vapor velocity vectors, and water streamlines velocity are investigated for both DI water and hybrid nanofluids.

5.1. Bubble Waiting Time Coefficient Correlations

In this simulation, new polynomial correlations for the BWTC are proposed to correct the quenching heat flux item under the RPI model by modifying this coefficient and correlating it to superheat temperature as an important parameter during the pool boiling process for DI water as well as hybrid nanofluids. The default setting for the BWTC is equal to 1 during the quenching boiling correction model, which is inside the RPI boiling model, and there is a possibility in Ansys Fluent to modify this coefficient as the user needs. However, the BWTC is a coefficient that is introduced in the quenching heat flux part to adjust the timing between the subsequent bubbles’ departures, and this coefficient can only be a constant value. To validate the present numerical data with experimental results as mentioned before, this coefficient was corrected through a trial-and-error method to modify the total heat flux through the quenching heat flux part due to its importance during the nucleate boiling regime, and then it correlated to the superheat temperature as shown in Figure 6. A simple fit was performed using polynomial fitting according to the Equation (31).

$$\mathrm{B}\mathrm{W}\mathrm{T}\mathrm{C}=Interecpt+\sum B_i\times ∆T_{sup}^i\left(i=1,2,3\dots .\right)$$

(31)

where BWTC and $∆T_{sup}$ are the waiting time coefficient of the bubble and the superheat temperature, respectively.

Table 4. The statistics and parameters of the proposed correlations.

Statistics and Parameters	DI Water	0.05 Vol.% HNFs	0.1 Vol.% HNFs
Number of points	5	5	5
Degrees of freedom	2	2	2
Residual sum of residual	6.497 × 10−5	0.0077	0.0069
R-square (COD)	0.999	0.966	0.962
Adj. R-square	0.998	0.932	0.923
Intercept value	−0.1499	−0.4313	0.2757
Intercept standard error	0.0235	0.1955	0.1909
B1 value	0.0594	0.2058	0.1312
B1 standard error	0.0059	0.0605	0.0545
B2 value	−8.797 × 10−4	−0.0095	−0.0046
B2 standard error	3.309 × 10−4	0.0041	0.0034

shows the statistics and parameters of the proposed correlations for DI water and two volume fractions of hybrid nanofluids.

5.2. Portions of the RPI Model

The portions of the total heat flux portioning sub-boiling model, which are the convection heat flux, evaporation heat flux, and quenching heat flux for both sub-boiling models (the classical model and the extended model), are plotted in Figure 7A,B for two concentrations of hybrid nanofluids. Both RPI sub-boiling models demonstrate that the element quenching heat flux plays an important role among other heat fluxes, particularly at the nucleate boiling regime as stated in previous investigations [24,25,31]. It can be noticed that the dissipation of heat flux by the quenching process compared to evaporation and convection mechanisms has a greater impact on the range of total heat flux (50–116 kW/m²), and this could be attributed to the dominance of the nucleate pool boiling regime and the mechanism of bubble formation. In the classical RPI model (before correction), the conditions for those portions were quite different after the improvement of the BWTC (extended model). When using the classical RPI model, the heat dissipation of the quenching part was near to that of the convective part, which means that the convective heat dissipation part dominated in this area and was similar with that of the quenching process at the conducted heat flux range. However, when correcting the BWTC by modifying the quenching heat flux value, the condition became quite different, with the convection part dominiating in the case of using both hybrid nanofluids, and this was due to the fact that, in this range of heat flux, there is a possibility for the domination of the natural convection process before the nucleation mechanism starts to dominate in large portion. It can be concluded that the quenching and convective parts dominate in the regions of boiling, as studied in the previous experimental investigations [17,18], and this was confirmed in this simulation.

5.3. Contours of Vapor Volume Fraction

Figure 8A depicts the vapor volume friction patterns for DI water for two different superheat temperatures and various time steps. It can be seen from the pictures that the bubbles start to nucleate from the heater’s surface at a time of 100 ms. Over time, bubbles continue to form from all surface directions (the bottom, sides, and top) in the case of the horizontal heater by a sliding process, and then the bubbles try to coalesce and separate from the top of the heater. At the time of step 400 ms, the vapor bubbles try to rise along the vertical direction, heading to the top of the chamber. By comparison with the higher superheat temperature for the same time steps, it is noticed that the formation and the movement of bubbles was faster at the high superheat temperature, and this could be attributed to the advanced nucleation regime in the case of high applied heat flux during the boiling process. Figure 8B shows the formation and growth of the bubbles at different time steps for approximately the same superheat temperatures for DI water and 0.05 vol.% hybrid nanofluids. It can be noted that the formation of the bubbles begins, and over time, the vapor volume fraction increases in both cases due to the nucleation process during the boiling phenomenon. In the comparison between the DI water and hybrid nanofluids, it can be seen that the void fraction for hybrid nanofluids was lower than that of DI water, and this is could be attributed to the deposition of nanoparticles during boiling on the heater surface, which enhances the wettability of the surface and then improves the heat transfer rate by delaying the growth of the bubbles and their coalescence from the heater. The vapor volume fraction of hybrid nanofluids, on the other hand, was smaller than that of DI water, which was owing to the combination of hybrid nanoparticles with bulk fluid bubbles, and this tendency was recognized in [24,25].

5.4. Vectors of Vapor Velocity

Figure 9A depicts the vectors of vapor velocity over time and superheat temperatures for DI water. It can be noted that the velocity of vapor bubbles through the formation and growth from the heater for both superheat temperatures increased over time, and this was due to the bubbles’ sliding process from the horizontal heater. The bottom side of the vapor bubbles had a velocity of around (0.2 m/s), and this was high enough to slide on both sides of the circular heater to head away to the top of the boiling vessel for a time step of 100 ms and higher temperature. It indicates that velocity vectors’ distribution during vapor bubble formation at the heater begins to build up as the applied heat flux (superheat temperature) increases at the same time steps. This could have been attributed to the intensified vapor columns from the sides heater along the vertical direction. However, the arrows showed that the velocity of the vapor bubbles increased over time as well as the superheat temperatures. It is noteworthy that the bubbles’ vapor velocity in the horizontal heater started to get bigger and faster at the sides of the heater, and then a column of bubbles was formed, vertically heading to the top of boiling chamber. This mechanism was visualized and captured in other investigations in literature [17,31], and therefore, it can be considered that the results of the present simulation physically match the experimental studies. Figure 9B illustrates the vectors of vapor bubbles velocity for DI water and 0.05 vol.% hybrid nanofluids at different time steps and approximately the same superheat temperatures. It can be seen from the compared pictures, that the velocity of vapor bubbles for hybrid nanofluids was faster than that for DI water, and this could be attributed to the deposition of nanoparticles on the heater surface, which results in a decreased bubble diameter and then an increase in the frequency of vapor bubbles.

5.5. Velocity of Water Streamlines

Figure 10A depicts the water velocity directions over time and different superheat temperatures for DI water. The results indicate that the agitation of vapor bubbles’ formation from the heater surface toward the top of the boiling chamber makes the water move in the vertical direction and these eddies will arrive back downward and replace the agitated water with fresh stationary water through liquid circulation, and this mechanism will increase as the vapor bubbles improve due to the large superheat temperature (more applied heat flux). A plausible explanation for the concept of quenching heat flux on the heater surface and water recalling near the departing bubbles is introduced during this process of the so-called transient quenching mechanism, and this fact is presented in the RPI sub-boiling model. Figure 10B presents the streamlines of water velocity for both DI water as well as the 0.05 vol.% hybrid nanofluids at different time steps and approximately the same superheat temperatures. It can be seen from the pictures that the velocity of hybrid nanofluids represented by vorticities helps to circulate the nanofluids faster than the DI water over time, and this is due to the frequency of vapor bubbles toward the top of the boiling chamber. This mechanism helps to push the nanofluids upward, and this results in the vortices and substituted the nanofluids with stagnant nanofluids coming from both sides of the chamber. However, the circulation of nanofluids enhanced the stability of nanofluids during the test, and this was confirmed by experimental results in the literature [17,19,31]; the circulation prevented the deposition of nanoparticles, especially at the top and sides of the heater, which makes the heat transfer better in case of hybrid nanofluids, and the results give a good indication the mechanism of bubble formation from the horizontal heater.

6. Conclusions and Future Direction

In the present work, the numerical simulation of pool boiling heat transfer for DI water and (Al₂O₃ + MgO) hybrid nanofluids from a horizontal heater at atmospheric pressure conditions was studied. The Eulerian–Eulerian approach by adopting the RPI sub-boiling model is utilized to predict the heat transfer behavior for various volume concentrations of hybrid nanofluids with a range of applied heat flux (15–120 kW/m²). The Rensselaer Polytechnic Institute (RPI) model was extended by modifying the quenching heat flux part to simulate the DI water and hybrid nanofluids from a horizontal copper heater. The obtained numerical data were validated with experimental studies in the literature. In this simulation, new BWTC correlations are presented by correlating BWTC to superheat temperature in the range of experimental conditions. The results demonstrated that, to predict boiling heat transfer behavior for hybrid nanofluids, efforts should take into account the bulk properties effect of hybrid nanofluids and the surface modification effect during the deposition of hybrid nanoparticles on the heater surface. In addition, the contours of vapor volume fraction, vectors of vapor velocity, and the streamlines of water velocity were introduced in this work. It was concluded that the superheat temperature has a significant effect on boiling parameters, and these parameters could give a considerable insight for predicting temperature measurements from a horizontal heater through the pool boiling process. Moreover, the obtained results show that the quenching heat flux item plays a vital role among other heat fluxes items. Finally, to introduce a comprehensive model of pool boiling using mono and hybrid nanofluids, more data from experimental investigations should be collected in the future to obtain more closure correlations related to the bubble dynamics parameters and surface improvement during pool boiling of hybrid nanofluids.