Enhancing Heat Transport in Multiphase Rayleigh-B\'enard Turbulence by Changing the Plate-Liquid Contact Angles

This numerical study presents a simple but extremely effective way to considerably enhance heat transport in turbulent multiphase flows, namely by using oleophilic walls. As a model system, we pick the Rayleigh-B\'enard setup, filled with an oil-water mixture. For oleophilic walls, e.g. using only $10\%$ volume fraction of oil in water, we observe a remarkable heat transport enhancement of more than $100\%$ as compared to the pure water case. In contrast, for oleophobic walls, the enhancement is then only about $20\%$ as compared to pure water. The physical explanation of the highly-efficient heat transport for oleophilic walls is that thermal plumes detach from the oil-rich boundary layer and are transported together with the oil phase. In the bulk, the oil-water interface prevents the plumes to mix with the turbulent water bulk. To confirm this physical picture, we show that the minimum amount of oil to achieve the maximum heat transport is set by the volume fraction of the thermal plumes. Our findings provide guidelines of how to optimize heat transport in thermal turbulence. Moreover, the physical insight of how coherent structures are coupled with one phase of a two-phase system has very general applicability for controlling transport properties in other turbulent multiphase flows.


Introduction
The key property of a turbulent flow is its ability to efficiently transport heat, mass, and/or momentum. Understanding this property allows to control the global transport, enhance or reduce it, which is not only of fundamental interest, but also highly relevant in various industrial applications. For example one can think of heat transfer enhancement for cooling applications (Dhir 1998) or drag reduction in pipe or channel flow (Bushnell 2015;Chung et al. 2021). A typical way to control integral properties is by manipulating coherent turbulent structures, which play a crucial role in the global transport (Smits et al. 2011;Holmes et al. 2012;Haller 2015;Graham & Floryan 2021). An example is in thermally driven turbulence where the global heat transport can be enhanced by manipulating the coherent thermal plumes -a major heat carrier in thermal convection Lohse & Xia 2010;Chilla & Schumacher 2012;Shishkina 2021).
Indeed, previous studies have proposed many different methods to enhance heat transport related to thermal plume manipulation: One method is to promote the detachment of plumes from the boundary layers by adding surface roughness (Shen et al. 1996;Ciliberto & Laroche 1999;Du & Tong 2000;Ahlers et al. 2009;Salort et al. 2014;Wagner & Shishkina 2015;Xie & Xia 2017;Zhu et al. 2017;Jiang et al. 2018) or adding shear (Bergman et al. 2011;Pirozzoli et al. 2017;Blass et al. 2020;Wang et al. 2020a). Another method is to confine the system in span-wise direction in order to increase the coherence of the thermal plumes and thus the heat transport (Huang et al. 2013;Chong et al. 2015Chong et al. , 2017.
However, also the addition of a new phase can enhance heat transport. For example, the mechanical injection of air bubbles (Gvozdić et al. 2018;Ng et al. 2020) or nucleating vapor bubbles by boiling (Dhir 1998;Zhong et al. 2009;Biferale et al. 2012;Lakkaraju et al. 2013;Guzman et al. 2016) are both quite efficient heat transfer enhancement strategies. While in the former case this is purely due to the extra mechanical stirring by the rising buoyant bubbles, in the latter case the vapor bubbles also act as direct heat carriers. But as shown by Wang et al. (Wang et al. 2019, 2020b, a considerable heat transfer can also emerge when the second phase has comparable density, even for small volume fraction of the second phase. In this study, we present another, novel, and simple way to substantially enhance heat transfer in wall-bounded thermally driven two-phase turbulence, namely by manipulating the wettability of the walls. We then reveal the underlying physics of the enhancement and in particular the crucial role of the interaction between the coherent structures (plumes) and the second phase. As paradigmatic example to study turbulent heat transfer in wall-bounded systems, we pick Rayleigh-Bénard (RB) convection Lohse & Xia 2010;Chilla & Schumacher 2012;Shishkina 2021). The system consists of a two-phase fluid in a box heated from below and cooled from above. Next to water, we use oil as a second phase, with its higher heat-conduction properties as compared to water. In our setup the top and the bottom walls are either both oleophilic or both oleophobic. With oleophilic walls, we reveal a novel flow dynamics, namely an "expressway" of oil for the thermal plumes. When thermal plumes detach from the oil-rich thermal boundary layers, they can be transported together with the oil phase. These oil-rich thermal plumes act as the dominant heat transport phase, bypassing the water phase. Thus the global heat transfer is significantly enhanced. We further reveal a strong correlation between the oil phase and the thermal plumes in the oleophilic cases, and show that the minimum amount of oil to achieve maximum heat transport is set by the volume fraction of thermal plumes.
The organization of this paper is as follows. The numerical methodology is introduced in Section 2. The main results on turbulent RB convection with multiphase is presented in Section 3. The mechanism of heat transfer enhancement is revealed in Section 4. The paper ends with conclusions and an outlook.

Methodology
The numerical method used here combines the phase-field model (Jacqmin 1999;Ding et al. 2007;Liu & Ding 2015;Soligo et al. 2021) and a direct numerical simulations solver for the Navier-Stokes equations namely AFiD, which is a second-order finitedifference open-source solver (Verzicco & Orlandi 1996;van der Poel et al. 2015). This method is developed to simulate the turbulence with two immiscible fluids and well validated in Liu et al. (2021a,b) by comparing the results to previous experimental, theoretical, and numerical results. The mass conservation is ensured by using the method of Wang et al. (2015). Various test cases in multiphase turbulence are extensively discussed in Liu et al. (2021b).

Governing equations
In the phase-field model, the evolution of the volume fraction of water, C, is governed by the Cahn-Hilliard equation, where u is the flow velocity, and ψ = C 3 −1.5C 2 +0.5C −Cn 2 ∇ 2 C the chemical potential. We set the Péclet number Pe = 0.9/Cn and the Cahn number Cn = 0.75h/H according to the criteria proposed in (Ding et al. 2007;Liu & Ding 2015), where h is the mesh size and H the domain height. The flow is governed by the Navier-Stokes equations, the heat transfer equation and the incompressibility condition, where θ is the dimensionless temperature. We define all dimensionless fluid properties (indicated by q), including the density ρ, the dynamic viscosity µ, the kinematic viscosity ν = µ/ρ, the thermal conductivity k, the thermal expansion coefficient β, the thermal diffusivity κ, and the specific heat capacity c p = k/(κρ), in a uniform way as follows, where λ q = q o /q w is the ratio of the material properties of oil and water. Here the surface tension force is defined as F st = 6 √ 2ψ∇C/(CnWe), and the gravity G = [C + λ β λ ρ (1 − C)] θj with j being the unit vector in vertical direction. Besides the property ratios, the global dimensionless numbers controlling the flow are Ra = β w gH 3 ∆/(ν w κ w ), Pr = ν w /κ w , and We = ρ w U 2 H/σ, where ∆ is the temperature difference between the bottom and top plates and U = √ β w gH∆ the free-fall velocity, and σ the surface tension coefficient. The response parameter is the Nusselt number Nu = QH/(k w ∆) with Q being the dimensional heat transfer. Note that the governing equations are normalized by the properties of water. Therefore, only global dimensionless numbers (e.g. Ra and Pr , without subscripts) and the ratios λ q of the material properties show up in eqs. (2.2) (2.4). More numerical details can be found in (Liu et al. 2021b,a).

Configurations of turbulent Rayleigh-Bénard convection with two phases
In this study, the control parameters are the volume fraction of the oil, α, the wettability of the wall surface (oleophobic and oleophilic), and the ratio of the Rayleigh number Ra o /Ra w . Here the subscripts "o" and "w" respectively indicate oil and water. The ratio Ra o /Ra w is varied by changing the ratio of thermal expansion coefficients, λ β , that of the thermal conductivity coefficients, λ k , and that of the dynamic viscosity, λ µ . Here Ra o > Ra w is used since the oil used has higher heat-conduction properties. Similarly to in the single-phase RB convection, we found that Nu always depends on the Rayleigh number, which depends on combinations of the fluid properties. Therefore, Ra o /Ra w is used here to characterize the ratio of fluid properties of oil and water. We take the density ratio λ ρ as 1 in an attempt to separate the various effects. Here we restrict the density variations only to temperature (i.e. oil and water are assumed to have the same density at the same temperature) and focus on the effects of β, k and µ. Note that the viscosity is assumed to be independent on temperature for simplicity (Oberbeck-Boussinesq approximation (Oberbeck 1879;Boussinesq 1903)). Other parameters are kept fixed, including the Weber number We (the ratio of inertia to surface tension), Ra w (the dimensionless temperature difference between the plates) and the Prandtl numbers Pr w and Pr o (a material property), at values allowing considerable turbulence and based on the properties of water and oil. We therefore fixed Ra w = 10 8 , Pr o = Pr w = 4.38, and We = 800. This value is e.g. obtained by choosing the realistic values of ρ w = 1000 kg/m 3 , σ = 0.021 N/m, U ≈ 0.4 m/s, and H ≈ 0.1 m. We varied 1 Ra o /Ra w 12. The latter value is motivated by Ra o /Ra w = 11.8 as obtained in a system consisting of KF-96L-1cs silicone oil and water. Next we varied the volume fraction 0% α 100% and the wettability (oleophilic or oleophobic), where we only took the two extreme contact angles 0 • and 180 • .
The simulations were performed in a three-dimensional cubic domain of size H 3 . Periodic boundary conditions were used in the horizontal directions, and no-slip walls at the top and bottom. The oleophobic and oleophilic conditions of the wall surface were realized by the Dirichlet boundary condition used in the phase-field model. We use a multiple resolution strategy (Liu et al. 2021b) in the simulation: the uniform mesh with 576 3 gridpoints to capture the two phases and the stretched mesh with 288 3 gridpoints for the velocity and temperature fields. The mesh is sufficiently fine and is comparable to corresponding single-phase flow studies ( , the black symbols denote the cases with Nu in agreement with the prediction of the Grossmann-Lohse theory with 100% oil (dash-dotted line corresponding to data with and dotted line to ♦), and the blue ones with Nu more than 2% smaller than the prediction.

Heat transfer in the RB convection with oleophilic and oleophobic walls
Initially, one big oil drop is placed at the center of the domain full of water. Then the oil drop breaks up into many smaller ones due to the thermally driven turbulent convection. As shown in figure 1(a), the oil phase is repelled by the oleophobic walls, whereas in figure  1(b), the oil phase is attracted by the oleophilic boundaries. This preferential wetting causes the boundary layer region to be occupied by water in the former and by oil in the latter case. These different near-wall flow structures significantly affect the heat transfer through the flow. In order to further elucidate this phenomenon, we have performed simulations with oleophobic and oleophilic walls, respectively, at various ratio Ra o /Ra w from 1 to 12 for fixed Ra w = 10 8 and α = 10%, as shown in figure 2(a). Since we define Ra o > Ra w , for all cases, both with the oleophobic walls and with the oleophilic walls, the heat transport is enhanced as compared to the case with pure water, but the amount of enhancements for the oleophilic cases is significantly larger than for the oleophobic cases.
In the oleophobic case, the amount of heat transfer enhancement linearly depends on Ra o /Ra w , and an enhancement of 21% is achieved at Ra o /Ra w = 12. Based on the observation, we define an effective Rayleigh number as Ra ef f = αRa o + (1 − α)Ra w . With this definition, the values of Nu for the oleophobic cases well agree with the Grossmann-Lohse (GL) (Grossmann & Lohse 2000, 2001Stevens et al. 2013) prediction for N u(Ra ef f ).
Next, we examine the more intriguing case with oleophilic walls. Surprisingly, the heat transfer enhancement in the oleophilic cases dramatically increases up to 103% at Ra o /Ra w = 12 (see figure 2a), which is 4 times larger than the enhancement in the oleophobic cases for the same parameters. It is remarkable to observe that the heat transfer efficiency can be as large as the cases with 100% oil, despite only 10% volume fraction of oil has been used. One question naturally arises: what is the minimum amount of oil required for the system to transfer the maximum heat, assumed as the value for 100% oil? To answer this question, we performed simulations with oleophilic walls for various α from 0% to 100% at fixed Ra o /Ra w = 4 and Ra w = 10 8 . For α < 8%, Nu in the oleophilic cases increases with α, while for α larger than this critical value (α c ≈ 8%), the heat transfer increase saturates to a value close to Nu o (the value for 100% oil) within ±2%, as shown in figure  2(b). The same result is also achieved at Ra o /Ra w = 8 and Ra w = 10 7 , but then with α c ≈ 10%.

Mechanism of heat transfer enhancement
What causes the significant enhancement of heat transfer in the cases with oleophilic walls despite using only a small amount of oil? The major reason is related to the thermal plumes, which are the primary heat carrier in thermal turbulence Lohse & Xia 2010;Chilla & Schumacher 2012;Shishkina 2021). To show the effect of the thermal plumes, we visualize them in figure 3(a) for the oleophobic cases, and in figure 3(b) for the oleophilic cases. In the cases with oleophobic walls, the thermal plumes mainly stay outside the oil phase after travelling into the bulk. The reason is that surface tension prevents the entrainment of external fluid into the oil phase. The only way for the heat carried by plumes to be transported across the interface of the two immiscible fluids is via thermal diffusion, the effects of which however are sufficiently small at such high Ra.
In the oleophobic cases the thermal plumes are mainly exposed to the turbulent bulk. In contrast, in the oleophilic cases, the thermal plumes detach from the oil-rich boundary layer (figure 3b). At the same time, the oil-water interface not only keeps thermal plumes inside the oil phase, but also prevents the heat carried by plumes to be mixed into the turbulent bulk. As a result, most hot (cold) thermal plumes can travel to the opposite plate with little heat loss (gain) to the turbulent bulk as they are thermally strongly coupled to the oil phase. As shown in figure 3(b), the oil carrying the hot plumes rises up in the shape of column, then breaks up into drops, and eventually coalesces with the oil in the upper boundary layer. During this process, the oil builds up an "expressway" to transfer the heat efficiently. Thus, the heat transfer in the oleophilic cases is always significantly enhanced up to the same value as in the system with 100% oil, and only a small amount of oil (larger than α c ) is required due to the strong coupling of the oil and the thermal plumes.
To quantitatively confirm this physical picture, we measure the amount of thermal plumes, using as thermal plumes definition condition |T − T | > (T − T ) 2 with T (x, t) being the local temperature and the spatial and temporal average. In figure  4, we show the volume fraction of the total thermal plumes in the domain, α plume , the ratio of thermal plumes in oil over the total thermal plumes, φ, and the correlation coefficient r xy = (x −x)(y −ȳ)/ (x −x) 2 (y −ȳ) 2 between the distribution of thermal plumes (x) and oil (y) in the whole domain, where the bar represents the spatial average and r xy ranges from 0 (no correlation) to 1 (perfect correlation).
For the oleophilic cases, for various Ra o /Ra w and fixed α = 10% (figure 4a), we observe φ = 80% ± 5% and r xy = 0.75 ± 0.05, whereas for the oleophobic cases we always have φ < 6% and r xy < 0.09. This finding shows that, in the oleophilic cases, most thermal plumes indeed are strongly coupled to the oil, contrasting to the weak coupling in the oleophobic cases.
Finally how to determine the critical volume fraction of oil α c to achieve maximum heat transfer efficiency? In figure 4(b), we show the oleophilic cases at various α and fixed (Ra w , Ra o ) = (10 7 , 8 × 10 7 ) and (10 8 , 4 × 10 8 ), and the critical condition for cases with Nu smaller than or close to Nu o , α c = α plume , (4.1) which well agrees with our numerical results. For α < α plume (blue symbols in figure  4b), although oil and the thermal plumes are still well coupled (r xy ≈ 0.75), there is not enough oil to generate all the plumes (φ ≈ 50%), leading to Nu smaller than Nu o (figure 2b). On the other hand, for α > α plume (black symbols in figure 4b), the excess amount of oil does not contribute to generating additional thermal plumes, and thus the heat transfer enhancement levels off and Nu remains close to Nu o as shown in figure 2(b). Therefore, α c is clearly determined by the volume of thermal plumes in the domain.

Conclusions and outlook
We have shown a novel way to significantly enhance the heat flux in wall-bounded thermally driven turbulence, namely by adding an oil phase (even of relatively small volume fraction) to the water and at the same time using oleophilic walls. The heat transport enhancement is brought about by the newly-found flow structure, where the thermal plumes are strongly coupled with the oil phase for cases with oleophilic walls. With the thermal plumes fully encapsulated inside the oil phase, the maximum enhancement can be obtained.
Our method and findings can easily and directly be applied and extended to other wall-bounded turbulent multiphase transport systems relevant in science and technology, way beyond thermal transport, namely to the turbulent multiphase mass or momentum transfer, such as drag in pipe or channel flow of oil-water mixtures (Pal 1993;Roccon et al. 2021), or of bubbly flow (Lu et al. 2005;Sanders et al. 2006;Ceccio 2010;Elbing et al. 2013;Murai 2014), and angular momentum transfer in turbulent multiphase Taylor-Couette flow (van den Berg et al. 2007;Murai et al. 2008;van Gils et al. 2013;Spandan et al. 2018;Srinivasan et al. 2015;Hu et al. 2017;Yi et al. 2021;Bakhuis et al. 2021). The reason is that what we manipulate by the wall coating and the qualitative and quantitative choice of the two phases are the coherent structures of the turbulent multiphase flow, and these determine not only the heat transfer in turbulent flow, but also mass and momentum transfer and therefore the overall drag.

Supplementary movies
Supplementary movies are available at URL