Conjugate Heat Transfer Predictions for Subcooled Boiling Flow in a Horizontal Channel Using a Volume-of-Fluid Framework

Predicting the heat transfer levels that occur in nucleate boiling is important for many engineering applications such as the cooling of micro-electronic devices and nuclear reactors [1–3]. This is of increasing importance for highly boosted IC engines which, at full load, relies heavily on a degree of nucleate boiling to operate at the limits of conventional liquid cooling systems [4–6].

Despite extensive research [7–11], the mechanism of nucleate boiling heat transfer is still not fully understood. As a result, mechanistic models, based on bubble formation and the bubble departure phenomenon, are not yet able to predict nucleate boiling heat transfer. To provide a solution to this difficulty, empirical correlations are generally used in practical situations [8–12] particularly in automotive applications [13–17].

Appropriate thermal boundary layer and vapor formation modeling by contrast is known to play a significant role in accurate predictions of heat transfer in the boiling regime [18,19]. Shala [19] used a “mixture” two-phase model with a mechanistic nucleate boiling model to predict the wall heat flux and its partitioning into separate heat transfer mechanisms (i.e., single-phase convection, wall quenching, and evaporation). The results agree reasonably well with experimental data for both a horizontal channel flow and flow in a vertical annulus, but significant under-prediction of heat flux occurs at higher velocities in both the convective and boiling regions. This shortcoming was associated with the turbulence model and the wall function used. Elsewhere, Yun et al. [20] used the k–ɛ model with a special logarithmic velocity wall function for both liquid and vapor phases, and a mechanistic bubble-size model in their two-fluid Eulerian computational fluid dynamic (CFD) simulation for subcooled boiling flow prediction. Although improved prediction of phase velocities was achieved with this wall function, they highlighted the need for improved turbulence models for two-phase boiling flows. And motivated by the possibility that an large eddy simulation (LES) approach to turbulence modeling may provide better thermal predictions in flow boiling situations, Deen et al. [21] investigated the performance of Eulerian-Eulerian Reynolds Averaged Navier-Stokes (RANS) and LES approaches for a liquid–gas flow (i.e., with no heat transfer). They showed that LES gave better agreement with experimental data than the RANS k–ε model simulations.

Despite the past decade seeing massive growth in CFD modeling with complex geometries, accurate prediction of pointwise heat transfer between a solid surface and a two-phase flow still remains a major challenge. In this paper, a two-phase LES approach for heat transfer prediction in flow boiling conditions is examined. This is accomplished by comparing RANS and LES simulations on a horizontal channel flow with heating from the lower surface. Simulations are undertaken by means of the finite volume CFD solver STAR_CCM+.

Liquid and vapor flows are simulated using the VOF Multiphase modeling approach within the Eulerian framework (originally proposed in Ref. [22]). One set of governing equations is solved for the mixture flow with shared velocity, pressure, and temperature fields, assumed for the two phases. An additional phase fraction equation is used to calculate the spatial phase concentrations within the domain. The equations are solved for an equivalent fluid whose physical properties are calculated as functions of the physical properties of its constituent phases and their volume fractions. The governing equations and the conservation equation that describe the transport of volume fraction of phase are fairly standard, and hence, will not be presented here, but can be found in Ref. [23].

where q_c is the convective heat flux, h_c is the convective heat transfer coefficient, T_w is the wall temperature, T_c is the local near-wall cell temperature, u_τ is the friction velocity, T⁺ is the dimensionless temperature, y_c is the normal distance of the near wall cell, and y⁺ is the dimensionless wall distance $uτyc/vf$ with $vf$ being the kinematic viscosity of the fluid.

A grid independency test has been undertaken for RANS simulations where a final grid of 434,300 cells is chosen with particular attention paid to the discretization of near-wall regions. Boundary layers are resolved by appropriate near-wall prism layers, with the near-wall thickness resolution in the region 2.0 < y⁺ < 25. Three turbulence models, the realizable k–ɛ model, the SST k–ω model, and the k–ɛ v2-f model, are employed with a hybrid wall function called the “All y⁺” wall approach in STAR_CCM+ [23] for near-wall turbulence quantities. In the RANS simulations, the mesh is not fine enough to resolve the interface between the liquid and the vapor phase. The interface, however, is captured directly in the LES simulations as a much finer mesh is used.

Large eddy simulation computations have been achieved on a 950 × 40 × 60 structured grid (i.e., ∼2,280,000 cells) with a refined near-wall grid to ensure the nearest wall cell y⁺ <1. The wall adapting local eddy subgrid scale model is employed [24]. The LES inlet flow is generated by mapping instantaneous velocity and turbulent field data tables extracted from a precursor channel flow simulation. The temporal and spatial fidelity of LES is ensured by applying bounded-central-differencing and second-order temporal discretization. The time-step size is chosen to maintain a convective Courant Number smaller than 1. For each thermal load, simulations are continued until enough samples have been collected to obtain the mean wall temperature.

where $μf$⁠, $cp,f$⁠, and $Prf$ are, respectively, the dynamic viscosity, the specific heat, and the Prandtl Number of the liquid phase; g is the gravitational acceleration, ρ_g is the density of the vapor-phase, $σ$ is the surface tension at the liquid-vapor interface, and T_w is the wall temperature.

The geometry and flow conditions of an experimental study [25] are used for the current RANS and LES simulations including conjugate heat transfer (CHT). The value of the empirical coefficient C_qw in Eq. (4) is a function of the particular liquid-surface combination. A chosen value of C_qw = 0.0029 is considered to represent the water-organic acid technology coolant mixture. At 90 °C inlet flow temperature, used in the present study, the properties of the mixture are ρ = 1038 kg/m³, μ = 0.00085 kg/ms, k = 0.424 W/mK, and C_p = 3620 J/kgK. The aluminum alloy was that used in the rig from which the experimental data were obtained in Refs. [14] and [25].

Figure 1 shows the geometry of the three-dimensional CFD/CHT model. Figure 2 shows a view of the grid for the fluid (top) and solid domains. The fluid is a mixture of water-organic acid technology coolant with 50% volumetric ratios and with the same physical properties as defined in Refs. [14] and [25] where a set of test conditions is available at different pressures, inlet velocities, and temperatures.

Simulations have been obtained at the highest boiling potential condition, i.e., the lowest coolant velocity of 0.25 m/s with an inlet temperature of 90 °C, a coolant pressure of 1.0 bar (absolute), and a corresponding saturation temperature of 108 °C. The RANS simulation inlet boundary is defined with a uniform inlet velocity profile. The inlet turbulence intensity has been assumed to be 5%, typical for low Re number pipe-flows (where Re = 5800 in the current study based on channel height and inlet condition). A convective outflow boundary condition is applied at the outlet. The heating power applied to the heating section (to match the experimental conditions) gives heat fluxes in the range 83–1300 kW/m².

Figure 3 shows the predicted wall temperatures compared with the experimental data from Ref. [25] at different heat flux levels. It can be seen that the temperature predictions obtained by different turbulence models of the RANS approach are compared very well with the experimental data, where the maximum difference is smaller than 3 °C except for the predictions obtained by the v2-f model at higher thermal loads. The LES predictions are a little closer to the experimental data, especially at the heat flux around 271 kW/m² where there is almost perfect agreement between the LES prediction and the experimental data. At lowest thermal loads, the measured wall temperature of 107.4 °C is just below the saturation temperature of the liquid T_sat at 108 °C. For this case, the heat transfer is essentially the result of pure convection (i.e., no boiling contribution), whereas the predicted wall temperature of 109.6 °C obtained by both RANS and LES is slightly above the saturation temperature. This means that the boiling model (Eq. (4)) is activated in the prediction to work out the heat transfer contribution from boiling. This can be confirmed in Fig. 4 which shows contours of the predicted vapor volume fraction obtained by LES and the realizable k–e model at the lowest heat flux, where it can be clearly seen that the predicted vapor volume fraction is not zero.

Although the LES simulations provide much better predictions of the turbulent flow field, in the present study, the complex process involved in the subcooled nucleate boiling flow, such as bubble nucleation, growth, detachment, coalescence, and collapse, cannot be directly captured by LES. Therefore, submodels for boiling using empirical correlations have been employed in both the LES and the RANS. As a consequence, the LES predictions of wall temperature show only marginal improvement over the RANS predictions as the temperature prediction depends strongly on the boiling submodels. Nevertheless, LES predictions can provide more information such as the instantaneous wall surface temperature distribution as shown in Fig. 5. Another distinguishing feature that can be observed is the instantaneous temperature distribution that is very different from uniform in the spanwise direction—there are a few hot spots scattered in the downstream region. This may have serious implications for real engineering application as the averaged wall surface temperature could be well below a critical value, whereas the instantaneous value might well not be. This kind of information can only be obtained from LES predictions.

It is evident from the current study that, in terms of the averaged wall temperature prediction, it is not essential to resolve the interface between the liquid and the vapor phase because the RANS predictions (without resolving the interface) are very close to the LES predictions (resolving the interface). Figure 6 shows that a clear interface can be seen in the LES prediction but not in the RANS prediction.

With regard to the v2-f model, the wall temperature is overpredicted at higher thermal loads. This may stem from the fact that boiling is underpredicted so there is less heat transfer contribution from boiling, leading to slightly lower total heat transfer predictions compared with other turbulence models. This results in slightly higher surface temperature. Figure 7 shows the predicted vapor volume fraction obtained by different turbulence models at the highest thermal loads of 1300 kW/m². It can be seen that the vapor volume fraction obtained by the v2-f model is slightly lower than the predictions obtained by the other two turbulence models.

Figure 8 shows contours of the predicted turbulent intensity obtained by the realizable k–e model. Predictions by other turbulence models are actually very similar. It is evident that turbulence is mainly generated in the heated surface area as a result of heat addition and boiling.

Numerical simulations have been undertaken for conjugate heat transfer predictions in thermal and flow conditions representative of coolant flow within IC engine cooling jackets. A multiphase mixture modeling approach with a volume-of-fluid interface capturing method has been employed for two phase flow combined with empirical correlations to predict heat transfer and phase change in the boiling regime. Both LES and RANS approaches have been employed in the study. The performance of LES and a number of different turbulence models in RANS are assessed by comparing the predicted wall temperatures with published experimental data. The main findings are as follows:

• The predicted wall temperatures obtained by both the LES and RANS approaches are in good agreement with measured data, apart from the predictions using the v2-f model at high values of heat flux. This suggests that the multiphase mixture modeling approach, combined with the empirical correlations for boiling, employed in the present study, is an appropriate choice for simulating such flows.

• Compared to the RANS models, LES predictions of the wall temperature are only marginally closer to the experimental data. This stems from the fact that the same boiling submodels are used in both the LES and RANS approaches. For subcooled nucleate boiling flow conditions (as in the present study, i.e., without flow separation and swirl), this finding suggests that it may not be advantageous to employ LES as the computational cost is significantly higher.

• It is not essential to resolve the interface between the liquid and the vapor phases. This is because the RANS predictions (without resolving interface) are very close to the LES predictions where the interface is resolved.

• Since Reynolds number in the current study is quite low, turbulence is mainly generated in the heated surface region as a consequence of heat addition and the action of boiling. This may be another reason for the relatively good agreement between RANS and LES predictions. However, this may not be the case for higher Reynolds number, where further study is needed.

The technical support is acknowledged of colleagues at Ford Dunton UK and Dearborn USA, Denso Italy, and the Ricardo Technical Centre Shoreham, UK.

• Engineering and Physical Sciences Research Council (Contract No. EP/M005755/1).