## Abstract

The accuracy of computational fluid dynamic (CFD)-based heat transfer predictions have been examined of relevance to liquid cooling of IC engines at high engine loads where some nucleate boiling occurs. Predictions based on (i) the Reynolds Averaged Navier-Stokes (RANS) solution and (ii) large eddy simulation (LES) have been generated. The purpose of these simulations is to establish the role of turbulence modeling on the accuracy and efficiency of heat transfer predictions for engine-like thermal conditions where published experimental data are available. A multiphase mixture modeling approach, with a volume-of-fluid interface-capturing method, has been employed. To predict heat transfer in the boiling regime, the empirical boiling correlation of Rohsenow is used for both RANS and LES. The rate of vapor-mass generation at the wall surface is determined from the heat flux associated with the evaporation phase change. Predictions via CFD are compared with published experimental data showing that LES gives only slightly more accurate temperature predictions compared to RANS but at substantially higher computational cost.

## 1 Introduction

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 [13]. 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 [46].

Despite extensive research [711], 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 [812] particularly in automotive applications [1317].

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+.

## 2 The Numerical Modeling Framework

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].

The convective part of the heat transfer is computed from [23]
$qc=hc(Tw−Tc)$
(1)
and
$hc=ρf|ycCp.f|ycuτT+y+|yc$
(2)

where qc is the convective heat flux, hc is the convective heat transfer coefficient, Tw is the wall temperature, Tc is the local near-wall cell temperature, uτ is the friction velocity, T+ is the dimensionless temperature, yc 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.

### 2.1 Phase-Change Model.

The phase change model used accounts for the onset of boiling through a number of submodels. The heat transfer at the wall-to-fluid boundary is used to calculate the phase change mass transfer rate. The vapor-phase temperature is assumed constant at the saturation temperature Tsat, and the liquid temperature Tf is approximated by the mixture temperature T. The total heat interchange between the liquid and vapor phases is then used to specify the mass transfer between phases, i.e.,
$m˙ec=CHTC×Area(T−Tsat)hfg$
(3)
where $CHTC×Area$ is the product of the heat transfer coefficient between the vapor bubbles and the adjacent liquid, with the contact area separating the two phases, where $hfg$ is the phase change enthalpy. Boiling occurs at the liquid–solid surface interface where the wall temperature is higher than Tsat. The surface heat flux qbw, as a result of boiling, is calculated using the empirical correlation of Rohsenow [8] given as follows:
$qbw=μfhfgg(ρf−ρg)σ(cp,f(Tw−Tsat)CqwhfgPrf1.7)3.03$
(4)

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 Tw is the wall temperature.

## 3 Model Parameters and Experimental Verification Geometry

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 Cqw in Eq. (4) is a function of the particular liquid-surface combination. A chosen value of Cqw = 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/m3, μ = 0.00085 kg/ms, k = 0.424 W/mK, and Cp = 3620 J/kgK. The aluminum alloy was that used in the rig from which the experimental data were obtained in Refs. [14] and [25].

### 3.1 Geometry and Boundary Conditions.

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.

Fig. 1
Fig. 1
Close modal
Fig. 2
Fig. 2
Close modal

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/m2.

## 4 Results

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/m2 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 Tsat 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.

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal

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.

Fig. 5
Fig. 5
Close modal

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.

Fig. 6
Fig. 6
Close modal

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/m2. 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.

Fig. 7
Fig. 7
Close modal

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.

Fig. 8
Fig. 8
Close modal

## 5 Conclusions

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.

## Acknowledgment

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

## Funding Data

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

## References

1.
Kandlikar
,
S. G.
,
2002
, “
Fundamental Issues Related to Flow Boiling in Minichannels and Microchannels
,”
Exp. Therm. Fluid Sci.
,
26
(
2–4
), pp.
389
407
.10.1016/S0894-1777(02)00150-4
2.
Mesquita
,
A. Z.
, and
Rodrigues
,
R. R.
,
2013
, “
Detection of the Departure From Nucleate Boiling in Nuclear Fuel Rod Simulators
,”
Int. J. Nucl. Energy
,
2013
, p. 950129.
3.
Suzuki
,
K.
, and
Kawamura
,
H.
,
2004
, “
Microgravity Experiments on Boiling and Applications: Research Activity of Advanced High Heat Flux Cooling Technology for Electronic Devices in Japan
,”
,
1027
(
1
), pp.
182
195
.10.1196/annals.1324.017
4.
Torregrosa
,
A. J.
,
Broatch
,
A.
,
Olmeda
,
P.
, and
Cornejo
,
O.
,
2014
, “
Experiments on Subcooled Flow Boiling in IC Engine-Like Conditions at Low Flow Velocities
,”
Exp. Therm. Fluid Sci.
,
52
, pp.
347
354
.10.1016/j.expthermflusci.2013.10.004
5.
Robinson
,
K.
,
Campbell
,
N. A. F.
,
Hawley
,
J. G.
, and
Tilley
,
D. G.
,
1999
, “
A Review of Precision Engine Cooling
,”
SAE
Paper No. 1999-01–0578.10.4271/1999-01-0578
6.
Ap
,
N.
, and
Tarquis
,
M.
,
2005
, “
Innovative Engine Cooling Systems Comparison
,”
SAE
Paper No. 2005-01–1378.10.4271/2005-01-1378
7.
Sanna
,
A.
,
Hutter
,
C.
,
Kenning
,
D. B. R.
,
Karayiannis
,
T. G.
,
Sefiane
,
K.
, and
Nelson
,
R. A.
,
2014
, “
Numerical Investigation of Nucleate Boiling Heat Transfer on Thin Substrates
,”
Int. J. Heat Mass Transfer
,
76
, pp.
45
64
.10.1016/j.ijheatmasstransfer.2014.04.026
8.
Rohsenow
,
W. M.
,
1952
, “
A Method of Correlating Heat Transfer Data for Surface Boiling of Liquids
,”
Trans. ASME
,
74
, pp.
969
976
.
9.
Le Martelot
,
S.
,
Saurel
,
R.
, and
Nkonga
,
B.
,
2014
, “
Towards the Direct Numerical Simulation of Nucleate Boiling Flows
,”
Int. J. Multiph. Flow
,
66
, pp.
62
78
.10.1016/j.ijmultiphaseflow.2014.06.010
10.
Mohanty
,
R. L.
, and
Das
,
M. K.
,
2017
, “
A Critical Review on Bubble Dynamics Parameters Influencing Boiling Heat Transfer
,”
Renewable Sustainable Energy Rev.
,
78
, pp.
466
494
.10.1016/j.rser.2017.04.092
11.
Sun
,
T.
,
Li
,
W.
, and
Yang
,
S.
,
2013
, “
Numerical Simulation of Bubble Growth and Departure During Flow Boiling Period by Lattice Boltzmann Method
,”
Int. J. Heat Fluid Flow
,
44
, pp.
120
129
.10.1016/j.ijheatfluidflow.2013.05.003
12.
Chen
,
J. C.
,
1966
, “
Correlation for Boiling Heat Transfer to Saturated Fluids in Convective Flow
,”
Ind. Eng. Chem. Process Des. Dev.
,
5
(
3
), pp.
322
329
.10.1021/i260019a023
13.
Steiner
,
H.
,
Brenn
,
G.
,
Ramstorfer
,
F.
, and
,
B.
,
2011
, “
Increased Cooling Power With Nucleate Boiling Flow in Automotive Engine Applications
,”
New Trends and Developments in Automotive System Engineering
,
InTech
, Rijeka, Croatia.
14.
Robinson
,
K.
,
Hawley
,
J. G.
, and
Campbell
,
N. A. F.
,
2003
, “
Experimental and Modelling Aspects of Flow Boiling Heat Transfer for Application to Internal Combustion Engines
,”
Proc. Inst. Mech. Eng. Part D
,
217
(
10
), pp.
877
889
.10.1243/095440703769683289
15.
Cardone
,
M.
,
Senatore
,
A.
,
Buono
,
D.
,
Polcino
,
M.
,
De Angelis
,
G.
, and
Gaudino
,
P.
,
2008
, “
A Model for Application of Chen's Boiling Correlation to a Standard Engine Cooling System
,”
SAE
Paper No. 2008–01–1817.10.4271/2008-01-1817
16.
Fontanesi
,
S.
, and
Giacopini
,
M.
,
2013
, “
Multiphase CFD-CHT Optimization of the Cooling Jacket and FEM Analysis of the Engine Head of a V6 Diesel Engine
,”
Appl. Therm. Eng.
,
2
(
2
), pp.
293
303
.10.1016/j.applthermaleng.2012.12.005
17.
Carpentiero
,
D.
,
Fontanesi
,
S.
,
Gagliardi
,
V.
,
Malaguti
,
S.
,
Margini
,
S.
,
Giacopini
,
M.
,
Strozzi
,
A.
,
Arnone
,
L.
,
Bonanni
,
M.
, and
Franceschini
,
D.
,
2007
, “
Thermo-Mechanical Analysis of an Engine Head by Means of Integrated CFD and FEM
,”
SAE
Paper No. 2007–24–0067.10.4271/2007-24-0067
18.
Fontanesi
,
S.
,
Carpentiero
,
D.
,
Malaguti
,
S.
,
Giacopini
,
M.
,
Margini
,
S.
, and
Arnone
,
L.
,
2008
, “
A New Decoupled CFD and FEM Methodology for the Fatigue Strength Assessment of an Engine Head
,”
SAE
Paper No. 2008–01–0972.10.4271/2008-01-0972
19.
Shala
,
M.
,
2012
, “
Simulation of Nucleate Boiling Flow Using a Multiphase Mixture Modelling Approach
,”
IMA J. Appl. Math
,
77
(
1
), pp.
47
58
.10.1093/imamat/hxr074
20.
Yun
,
B.-J.
,
Splawski
,
A.
,
Lo
,
S.
, and
Song
,
C.-H.
,
2012
, “
Prediction of a Subcooled Boiling Flow With Advanced Two-Phase Flow Models
,”
Nucl. Eng. Des.
,
253
, pp.
351
359
.10.1016/j.nucengdes.2011.08.067
21.
Deen
,
N. G.
,
Solberg
,
T.
, and
Hjertager
,
B. H.
,
2001
, “
Large Eddy Simulation of the Gas–Liquid Flow in a Square Cross-Sectioned Bubble Column
,”
Chem. Eng. Sci.
,
56
(
21–22
), pp.
6341
6349
.10.1016/S0009-2509(01)00249-4
22.
Hirt
,
C. W.
, and
Nichols
,
B. D.
,
1981
, “
Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries
,”
J. Comput. Phys.
,
39
(
1
), pp.
201
225
.10.1016/0021-9991(81)90145-5
23.
2016
, “STAR-CCM+ User Guide, Version 11,” CD-ADAPCO, Melville, NY.
24.
Nicoud
,
F.
, and
Ducros
,
F.
,
1999
, “
Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor
,”
Flow, Turbul. Combust
,
62
(
3
), pp.
183
200
.10.1023/A:1009995426001
25.
Robinson
,
K.
,
2001
, “
IC Engine Coolant Heat Transfer Studies
,”
Ph.D. thesis
, University of Bath, Bath, UK.http://ethos.bl.uk/OrderDetails.do?uin=uk.bl.ethos.275444