Investigation of Forced Convective and Subcooled Flow Boiling Heat Transfer Coefficients of Water-Ethanol Mixture: Numerical Study

The current technology in the nuclear industries, petroleum refineries, chemical industries, process industries, automobile industries, refrigeration systems etc. endeavor to meet the demand and supply in daily life of mankind. In these industries, heat is dissipated from the catalytic reactors, batteries, electronic devices, burners, radiators etc. The forced convective and subcooled flow boiling of water-ethanol mixture is pertinent to the operation of heat dissipative devices like small catalytic reactors, electronic devices and HEV battery module. The heat transfer coefficient of water-ethanol mixture is required to design cooling equipments for the heat dissipative devices [1-3]. During the past decades, only experimental studies to determine the heat transfer coefficient had been performed. But performing experiment will be tedious and therefore the alternative path is to determine by numerical methods. The experimental studies to date on the subcooled flow boiling have yielded only limited information on the basic two-phase hydrodynamic characteristics of water-ethanol mixture. Subcooled boiling is a kind of boiling that reaches the saturation temperature in the vicinity of a heated wall with the bulk temperature being below the wall temperature [4, 5]. On the other hand, benefit from the numerical method can be considered as alternative tool to provide useful physical information with limited cost. The numerical schemes for the prediction of two-phase flows can be classified into surface and volume methods [6-8]. Surface methods maintain sharp interfaces throughout the calculation. A disadvantage of these methods is that special treatment needs to be implemented to solve for interfaces which are exposed to large deformations. In volume methods, the different fluids are marked by an indicator function known as volume fraction or a level set [9-12]. The advantage of these methods is their ability to deal with arbitrarily shaped interfaces and to handle large deformations. The convective scalar transport is discretized for the bubble void fraction equations with a differencing scheme that guarantees physical values and thus preventing streaking of the transitional area over the grids [13, 14]. Bubble void fraction is the presence of bubble in the subcooled fluid and its presence affects the thermo physical and thermodynamic properties of the fluid.

Detailed investigation on subcooled flow boiling of water-ethanol mixture is limited in literature which is essential to design the cooling devices. In view of this, the understanding of the bubble behavior and its effect on heat transfer significantly contributes to a better understanding of physical phenomena in subcooled flow boiling [15, 16]. Hence flow visualization is essential to study the bubble dynamics of water-ethanol mixture during subcooled flow boiling. The experimental studies on the flow boiling have shown limited information on the basic two-phase hydrodynamic characteristics of binary mixtures. The benefit from the recent studies on numerical analysis can be considered as an alternative path to provide useful physical information with limited cost. However, the detailed literature shows that the numerical analysis for determining the forced convective and subcooled flow boiling heat transfer coefficients for water-ethanol mixtures are not available. The two-fluid model essentially demands momentum equation for each of the phases. The interfacial transport of momentum and heat depends on the relations derived from the empirical equations and are not accurate till the date [17]. In view of the great significance of the reliable data for water-ethanol mixture as an alternative coolant, an independent study is undertaken here to determine the subcooled flow boiling characteristics of this mixture.

The present work consists of determining the forced convective and subcooled flow boiling heat transfer coefficients of water ethanol mixture. Heat transfer related to forced convection and subcooled flow boiling of water-ethanol mixture is determined from mathematical modelling and numerical simulation. The Thermophysical and thermodynamic properties of the subcooled boiling fluid of water-ethanol mixture are determined from mathematical modelling. The bubble void fraction is determined and substituted in the mixing rule equation of thermodynamic and Thermophysical properties. The thermodynamic and Thermophysical properties are substituted in the x-momentum and energy equation to determine the values of pressure drop, velocity and temperature of the fluid. From the temperature values, the subcooled flow boiling heat transfer coefficient is obtained and compared with that of experimental results.

3.1 Forced convective heat transfer coefficient: Mathematical modelling  

The continuity equation is represented by Eq. (1)

$\frac{\partial u}{\partial x}=0$    (1)

The x-momentum-Navier Stokes equation is represented by Eq. (2).

$\frac{\partial u}{\partial t}=-\frac{\partial p}{\partial x}+\mu\left(\frac{\partial^{2} u}{\partial y^{2}}\right)$      (2)

The energy equation is represented by Eq. (3).

$C_{p} \frac{\partial T}{\partial t}=-u \frac{d p}{d x}+k \frac{\partial^{2} T}{\partial x^{2}}+k \frac{\partial^{2} T}{\partial y^{2}}+\mu\left(\frac{d u}{d y}\right)^{2} v$     (3)

The wall temperature determined from the experiment is considered as the Dirichlet boundary condition. The pressure drop and velocity in the x-momentum equation are determined by pressure correction method. The determined values of pressure drop and velocity are substituted in the energy equation and solved by the Lax Wandroff method to determine the temperature of the fluid. The forced convective heat transfer coefficient is determined from the known values of wall temperature, heat flux and determined values of fluid temperature as shown in Eq. (4).

$h=\frac{q^{\prime \prime}}{\left(T_{W}-T_{f}\right)}$     (4)

Specific heat of the mixture is calculated using simple mixing rule. Thermal conductivity and liquid viscosity are calculated by Flippov and McLaughlin Equation [21, 22] represented in the Eq. (5) and Eq. (6).

$\frac{k_{m}-k_{i}}{k_{j}-k_{i}}=C m_{f j}^{2}-m_{f i}(1-C)$     (5)

$\ln \left(\mu_{m}\right)=x_{i} \ln \mu_{i}+x_{j} \ln \mu_{j}$    (6)

Volume of fluid method (VOF) is adopted in the present analysis. In this analysis, the thermophysical and thermodynamic properties of the subcooled boiling fluid are determined by incorporating bubble void fraction (α) by simple mixture rule. The thermodynamic and thermophysical properties are substituted in the x-momentum and energy equation to determine the values of pressure drop, velocity and temperature of the fluid. From the temperature values, the subcooled flow boiling heat transfer coefficient is determined and compared with that of experimental results. The code is developed by using Matlab R 2018 programming to solve for heat transfer coefficient. The sample image of bubble formation of water at heat flux=90.4 kW/m², mass flux=76.67 kg/m²-s and inlet temperature=303 K is shown in Figure 7. This sample is selected from the experiment. The channel containing the bubbles are divided into suitable grids to calculate the bubble void fraction. The conservative form of scalar convection equation for the bubble void fraction is given by Eq. (7).

$\frac{\partial \alpha}{\partial t}=u \frac{\partial \alpha}{\partial x}$    (7)

The Crank Nicolson implicit scheme is used as technique for discretization as shown in Eq. (8).

$\frac{\alpha_{i}^{t+1}-\alpha_{i}^{t}}{\Delta t}=u_{i n}\left\{\frac{\left(\frac{\alpha_{i+1}^{t+1}+\alpha_{i+1}^{t}}{2}\right)-\left(\frac{\alpha_{i}^{t+1}+\alpha_{i}^{t}}{2}\right)}{\left(\frac{\Delta x}{2}\right)}\right\}$     (8)

7.png

Figure 7. Bubble formation of water at Heat flux=90.4 kW/m² and mass flux=76.67 kg/m²-s

Algebraic equations are obtained from Eq. (6) is solved to determine bubble void fraction by using Tridiagonal Matrix Algorithm (TDMA). The bubble void fraction is solved in the control surface as upwind, donor and accepter cells with the face values in between the interfaces of grid as shown in Figure 8. Initially the values are assumed suitably varying from 0 to 1, depending upon occupation of bubbles in the channel surface. The bubble void fraction with values more than one and less than zero are obtained. Hence the corrector-predictor steps are involved to solve these void fractions. The new amount of fluid to be convected over the face is determined by subtracting the unboundedness error from the original amount of fluid convected over the face [23].

8.jpg

Figure 8. Control area to solve bubble void fraction

3.1.1 Prediction of $\alpha$ in grid centre by corrector predictor method

Following steps are followed when $\alpha_{D}$ values obtained exceed 1 or negative.

$\beta_{f}=\left\{\min \left\{\frac{\cos \left(2 \theta_{f}+1\right)}{2}\right\}, 1.0\right\}$     (9)

$\beta_{f}$ is called weighing factor. The weighing factor incorporates the unboundedness error while calculating the bubble void fraction at the grid center and the face.

$\Delta \alpha=\frac{\alpha_{A}{ }^{t}+\alpha_{A}{ }^{t+1}}{2}-\frac{\alpha_{D}{ }^{t}+\alpha_{D}{ }^{t+1}}{2}$    (10)

$\beta_{f}^{\prime}=\min \frac{E^{-}\left(2+C_{f}-2 C_{f} \beta_{f}\right)}{2 C_{f}\left(\Delta \alpha^{*}-E^{-}\right)}, \beta_{f} \quad$ When $\Delta \alpha^{*}>E^{-}$    (11)

$\beta_{f}^{\prime}=0$, When $\Delta \alpha^{*}<E^{-}$    (12)

$E^{-}=\max \left\{-\alpha_{\mathrm{D}}^{\mathrm{t}+\delta \mathrm{t}}, 0\right\}$ When $\alpha_{D}^{t+1}<0$     (13)

$E^{+}=\max \left\{\alpha_{\mathrm{D}}^{\mathrm{t}+\delta \mathrm{t}}-1,0\right\}$ When $\alpha_{D}^{t+1}>1$    (14)

$\beta_{f}^{\prime}=\min \frac{E^{-}\left(2+C_{f}-2 C_{f} \beta_{f}\right)}{2 C_{f}\left(-\Delta \alpha^{*}-E^{+}\right)}, \beta_{f}$ When $-\Delta \alpha^{*}>E^{+}$     (15)

$\beta_{f}^{\prime}=0$, When $-\Delta \alpha^{*}<E^{+}$    (16)

$\beta_{f}^{*}=\beta_{f} \sim \beta_{f}^{\prime}$    (17)

The corrected weighting factor $\beta_{f}^{*}$ should be always less than or equal to the previous weighting factor. Otherwise, the contribution of the downwind cell starts to increase and so also the degree of unboundedness increases. The lower limit on $\beta_{f}^{*}$ remains zero and this is applied to Eq. (18) to obtain bounds for $\beta_{f}^{\prime}$.

$\alpha^{* *}{ }_{f}=\left(1-\beta_{f}^{*}\right)\left[\frac{\alpha_{D}^{t}+\alpha_{D}^{t+\delta t}+E^{-}}{2}\right]+\left[\frac{\alpha_{A}^{t}+\alpha_{A}^{t+\delta t}+E^{-}}{2}\right]$    (18)

Or

$\alpha_{f}^{* *}=\left(1-\beta_{f}^{*}\right)\left[\frac{\alpha_{D}^{t}+\alpha_{D}^{t+\delta t}+E^{+}}{2}\right]+\left[\frac{\alpha_{A}^{t}+\alpha_{A}^{t+\delta t}+E^{+}}{2}\right]$   (19)

where, $\alpha_{f}^{* *}$ is the new face value and and $E^{-}$ and $E^{+}$ the magnitude of the unbounded void value.

$\alpha_{f}^{*}=\alpha_{f}^{* *} \pm \frac{E^{-}}{c}$     (20)

Or

$\alpha_{f}^{*}=\alpha_{f}^{* *} \pm \frac{E^{+}}{c}$    (21)

$\alpha_{f}^{\prime \prime}=\left(1-\beta_{f}\right) \frac{\alpha_{D}^{t}+\alpha_{D}^{t+\delta t}}{2}+\beta_{f} \frac{\alpha_{A}^{t}+\alpha_{A}^{t+\delta t}}{2}$    (22)

Eq. (23) is used for calculating $\alpha_{D}$ for next time step

$\alpha_{D}^{c}=\frac{2 \alpha_{f}^{*}+\beta_{f} \alpha_{D}^{t}-\alpha_{D}^{t}-\beta_{f} \alpha_{A}^{t}-\beta_{f} \alpha_{A}^{t+\Delta t}}{\left(1-\beta_{f}\right)}$     (23)

3.1.2 Prediction of $\alpha$ in face centre by corrector predictor method

Following steps are followed when $\alpha_{f}$ (face values) obtained exceed 1 or negative.

$\widetilde{\alpha_{D}}=\frac{\alpha_{D-\alpha_{U}}}{\alpha_{A}-\alpha_{I I}}$      (24)

$\widetilde{\alpha_{f}}=\frac{\alpha_{f}-\alpha_{U}}{\alpha_{A}-\alpha_{U}}$     (25)

$\tilde{\alpha}_{f C B C}=\min \left\{\frac{\widetilde{\alpha_{D}}}{c}, 1.0\right\}$ When $0 \leq \tilde{\alpha}_{D} \leq 1$     (26)

$\tilde{\alpha}_{f C B C}=\widetilde{\alpha_{D}}$ When $\tilde{\alpha}_{D}<0, \widetilde{\alpha_{D}}>1$    (27)

$\widetilde{\alpha_{f U Q}}=\min \left\{\frac{8 c{\widetilde{\alpha}}_{D}+(1-c)\left(6 \widetilde{\alpha_{D}{+3}}\right)}{8}\right\}$ When $0 \leq \tilde{\alpha}_{D} \leq 1$    (28)

$\widetilde{\alpha_{f U Q}}=\tilde{\alpha}_{D}$ When $\tilde{\alpha}_{D}<0, \widetilde{\alpha_{D}}>1$    (29)

The above derivation for the volume of fluid method is a resolution scheme to be carried out only in one-dimension.

$\alpha_{f}^{c}=\beta_{f} \tilde{\alpha}_{f C B C}+\left(1-\beta_{f}\right) \tilde{\alpha}_{f U Q}$    (30)

where, $\tilde{\alpha}_{f B C}$ is bubble void fraction in face center to satisfy boundedness criteria. $\tilde{\alpha}_{\text {fQuick }}$ is bubble void fraction in face center to satisfy conservative criteria by QUICK scheme.  When  is more than 1 or less than 0 in the grid face and grid centre, the average values obtained from Eq. (28) and Eq. (36) are considered for the new value at the grid face. 

3.1.3 Mixture rule

The viscosity, specific heat and thermal conductivity are determined by simple mixture rule as given by Eqns. (31), (32) and (33).

$\mu=\alpha \mu_{l m}+(1-\alpha) \mu_{v m}$     (31)

$k=\alpha k_{l m}+(1-\alpha) k_{v m}$      (32)

$C_{p}=\alpha C_{p l m}+(1-\alpha) C_{p v m}$    (33)

The pressure drop and velocity in the Eq. (2) are determined by pressure correction method. The values of pressure drop and velocity are substituted in the Eq. (3) and solved by the Lax Wandroff method to determine the temperature of the fluid. The subcooled flow boiling heat transfer coefficient is determined from the known values of wall temperature, heat flux and determined values of fluid temperature by Eq. (4).

3.1.4 Pressure correction technique and Lax Wandroff Explicit scheme

The procedure to solve the momentum and energy equations by using Pressure correction method and Lax Wandroff Explicit method is discussed in detail in the literature [24].

3.2 Grid independent study

The pressure drop, velocity and temperature are calculated by numerical techniques such as pressure correction method and the Lax Wandroff explicit method for mathematical modelling. The grid independence for the pressure drop of water at heat flux=90.4 kW/m² and mass flux=76.67 kg/m²-s and inlet temperature=303 K is checked. The change in pressure drop and temperature of water are not very significant. The grid size of 15 X 8 are chosen in the present study. During the numerical study, it is found that there is no significant variation in the single phase heat transfer coefficient with grid numbers.

The photographic images of bubbles obtained using high speed camera are used to determine the bubble void fractions. The interaction between liquid and local vapour is analysed by solving the bubble volume of fraction in the numerical study. This bubble volume of fraction is considered in the momentum and energy equations to determine the heat transfer coefficient of the binary mixture. The values obtained are compared with that of experiment. Following conclusions are arrived from the present work.

• It is observed that the experimentally determined subcooled flow boiling heat transfer coefficient increases with the addition of ethanol to water initially upto 25% ethanol volume fraction, but at 50% and 75% ethanol volume fractions the heat transfer coefficient reduces. The pure ethanol has marginally higher value than the mixture of ethanol volume fraction 75%.
• From the numerical analysis it is concluded that the addition of ethanol to water decreases the forced convective and subcooled flow boiling heat transfer coefficient of the water-ethanol mixture.
• It is found that the forced convective heat transfer coefficient of water obtained from numerical simulation deviated by 17.1% from that of the experiment and 21.3 % from that of mathematical modelling.
• The average deviation between the experimentally determined and numerically determined subcooled flow boiling heat transfer coefficient of water ethanol-mixture is 24.13%.

The possible future research can be carried out experimentally and numerically as follows:

• The present experiment may be extended upto critical heat flux to determine the heat transfer coefficient in the different regimes of boiling.
• The experiment can be conducted with different binary mixtures. Repeating these experiments with different binary mixtures for a wide range of aspect ratios would further expand understanding the channel geometries which influences the bubble departure.
• The forces acting on the bubbles during subcooled boiling liquid in low-aspect ratio and microchannels may be identified.
• The recently discovered confinement pressure effects deserve extensive experimentation to determine the degree to which these effects influence the flows. Such studies should attempt to demonstrate the additional nucleation due to bubble induced water hammer propagation and add to the dataset of bubble growth rates for a broader range of channels cross-sections and bubble pressures.
• The development of fiber optic sensors will further extend the experimental parameter set to include slug velocity, size, and growth rates as well as allow the first maps in time and space of liquid temperature and void fraction. Such maps would be extremely valuable in understanding the flow instabilities during the boiling.