Heat transfer in drop-laden turbulence

Abstract Heat transfer by large deformable drops in a turbulent flow is a complex and rich-in-physics system, in which drop deformation, breakage and coalescence influence the transport of heat. We study this problem by coupling direct numerical simulation (DNS) of turbulence with a phase-field method for the interface description. Simulations are run at fixed-shear Reynolds and Weber numbers. To evaluate the influence of microscopic flow properties, like momentum/thermal diffusivity, on macroscopic flow properties, like mean temperature or heat transfer rates, we consider four different values of the Prandtl number, which is the momentum to thermal diffusivity ratio: $Pr=1$, $Pr=2$, $Pr=4$ and $Pr=8$. The drop volume fraction is $\varPhi \simeq 5.4\,\%$ for all cases. Drops are initially warmer than the turbulent carrier fluid and release heat at different rates depending on the value of $Pr$, but also on their size and on their own dynamics (topology, breakage, drop–drop interaction). Computing the time behaviour of the drops and carrier fluid average temperatures, we clearly show that an increase of $Pr$ slows down the heat transfer process. We explain our results by a simplified phenomenological model: we show that the time behaviour of the drop average temperature is self-similar, and a universal behaviour can be found upon rescaling by $t/Pr^{2/3}$. Accordingly, the heat transfer coefficient $\mathcal {H}$ (respectively its dimensionless counterpart, the Nusselt number $Nu$) scales as $\mathcal {H}\sim Pr^{-2/3}$ (respectively $Nu\sim Pr^{1/3}$) at the beginning of the simulation, and tends to $\mathcal {H}\sim Pr^{-1/2}$ (respectively $Nu\sim Pr^{1/2}$) at later times. These different scalings can be explained via the boundary layer theory and are consistent with previous theoretical/numerical predictions.


Introduction
Transport of passive and active scalars in multiphase turbulence is very important in many industrial processes and natural phenomena, from vaporization of atomized fuel jets (Gorokhovski & Herrmann 2008;Ashgriz 2011;Gao et al. 2022;Boyd & Ling 2023), to rain formation and atmosphere-ocean heat/mass exchanges (Duguid & Stampfer Jr 1971;Deike 2022) or even to the uptake of nutrients and other biochemicals by cells in complex flows (Aksnes & Egge 1991;Magar & Pedley 2005).While the mixing of active or passive scalars in turbulent single-phase flows has been extensively analyzed (Antonia & Orlandi 2003;Kasagi et al. 1992;Kim & Moin 1989;Warhaft 2000;Pirozzoli et al. 2016;Zonta et al. 2012a,b;Zonta & Soldati 2014), it remains a challenging task in turbulent multiphase flows (Gauding et al. 2022;Ni 2024), where most of the available studies considered the case of heat/mass transfer from/to isolated drops and bubbles (Boussinesq 1905;Levich 1962;Bird et al. 2002;Bothe et al. 2004;Figueroa-Espinoza & Legendre 2010), with some remarkable 2 F. Mangani, A. Roccon, F. Zonta, and A. Soldati exceptions -focusing on swarms of drops/bubbles -appeared only recently (Herlina & Wissink 2016;Albernaz et al. 2017;Dodd et al. 2021;Farsoiya et al. 2021;Shao et al. 2022;Scapin et al. 2022;Hidman et al. 2023;Farsoiya et al. 2023).One of the crucial aspects in turbulent multiphase flows, which makes these flows so difficult to analyze, is the presence of interfaces that dynamically move and deform in time and space according to the flow conditions, and that mediate heat/species transport, mixing and phase change phenomena (Deckwer 1980;Gvozdić et al. 2018;Liu et al. 2022;Pelusi et al. 2023).
In this work, we focus on the numerical simulation of the heat transfer process in a droplet laden turbulent channel flow, looking in particular at the role of the Prandtl number , i.e. the ratio between momentum and thermal diffusivity, in the process.Compared to single-phase turbulence, where the range of scales that must be resolved to perform a direct numerical simulation (DNS) is purely dictated by the smallest scales of turbulence (Kolmogorov scale), when the mixing of scalars in multiphase turbulence is analyzed, two further additional scales come into the picture.The first one is the Batchelor scale (Batchelor 1959;Batchelor et al. 1959), which determines the smallest scale of the temperature/concentration field.The second important scale is the Kolmogorov-Hinze scale (Kolmogorov 1941;Hinze 1955), and is linked to the multiphase nature of the flow.This scale can be used, perhaps with some limitations (Qi et al. 2022), to determine the critical size of a drop/bubble that will not undergo breakage in turbulence.These two scales -and their corresponding ratio to the Kolmogorov scale, i.e. the smallest length scale of the turbulent flow field -control the system dynamics and define the minimal grid requirements that must be satisfied to perform a DNS of scalar mixing in multiphase turbulence (keeping always in mind that performing a simulation that resolves the interface dynamics down to the molecular scale is at present almost unfeasible).In this context, the major constraint is usually posed by the Batchelor scale, which becomes smaller than the Kolmogorov length scale when Prandtl numbers large than unity are considered.Overall, the wide range of scales involved in the process makes simulations of scalar mixing in multiphase turbulence a challenging task and limits the space parameters that can be explored by means of direct numerical simulations.Our simulations are initialized injecting a swarm of large and deformable drops (initially warmer) inside a turbulent channel flow (initially colder).The system is described by coupling the direct numerical simulations of turbulent heat transfer with a phase-field method, employed to describe the drops topology (Zheng et al. 2015;Mirjalili et al. 2022).We simulate realistic values of the Prandtl number up to  = 8, similar to those obtained in liquid-liquid systems.We remark here that simulations of mass transfer problems in wall-bounded flow configurations, where the typical Schmidt number  (i.e. the mass transfer counterpart of ) is O (10 2 ∼ 10 3 ), e.g. ≃ 600 for  2 in freshwater (Wanninkhof 1992), are currently out of reach even using the most advanced computing.Indeed, the resulting Batchelor scale would be at least one order of magnitude smaller, thus requiring grid resolutions comparable or larger than those employed for state-of-the-art single-phase DNS (Lee & Moser 2015;Pirozzoli et al. 2021) but with a much larger computational cost as the systems of equations to be solved is more complex and restrictive (also from the temporal discretization point of view).
The present study has three main objectives.First, we want to investigate the macroscopic dynamic of the drops and of the heat transfer process by analyzing the drop size distribution and the mean temperature behavior of the two phases over time.Second, we want to characterize the influence of the Prandtl number, i.e. of the microscopic flow properties, on the macroscopic flow properties (mean temperature, heat transfer coefficient) and, building on top of the numerical results, we want to develop a physically-based model to explain the observed results.Third, we want to study the influence of the Prandtl number and of drop size on the temperature distribution inside the drops, so to evaluate the corresponding flow mixing/ homogenization.
The paper is organized as follows.In section §2, the governing equations, the numerical method, and the simulation setup are presented.In section §3, the simulation results, in terms of drop size distribution and mean temperature of the two phases and heat transfer coefficient are carefully characterized and discussed.A simplified model is also developed to explain the observed results.The temperature distribution inside the drops is then evaluate at different Prandtl numbers and drop sizes.Finally, conclusions are presented in §4.

Methodology
We consider a swarm of large and deformable drops injected in a turbulent channel flow.The channel has dimensions   ×   ×   = 4ℎ × 2ℎ × 2ℎ along the streamwise (), spanwise () and wall-normal direction ().To describe the dynamics of the system, we couple direct numerical simulation (DNS) of the Navier-Stokes and energy equations, used to describe the turbulent flow, with a phase-field method (PFM), used to describe the interfacial phenomena.The employed numerical framework is described more in detail in the following.

Phase-field method
To describe the dynamics of drops and the corresponding topological changes (e.g.coalescence and breakage), we employ an energy-based phase field method (Jacqmin 1999;Badalassi et al. 2003), which is based on the introduction of a scalar quantity, the phase field , required to identify the two phases.The phase field  has a uniform value in the bulk of each phase ( = +1 inside the drops;  = −1 inside the carrier fluid) and undergoes a smooth change across the thin transition layer that separates the two phases.The transport of the phase field variable is described by a Cahn-Hilliard equation, which in dimensionless form reads as: where u = (, , ) is the velocity vector, Pe is the Péclet number,  is the phase field chemical potential while   is a penalty-flux term which will be further discussed later.The Péclet number is where  *  is the friction velocity ( *  = √︁  *  / * , with  *  the wall-shear stress and  * =  *  =  *  the density of the fluids), ℎ * is the channel half-height, M * is the mobility and  * is a positive constant (the superscript * is used to denote dimensional quantities hereinafter).The chemical potential  is defined as the variational derivative of a Ginzburg-Landau free-energy functional, the expression of which is chosen to represent an immiscible binary mixture of fluids (Soligo et al. 2019a,b,c).The functional is the sum of two contributions: the first contribution,  0 , accounts for the tendency of the system to separate into the two pure stable phases, while the second contribution,   , is a mixing term accounting for the energy stored at the interface (i.e.surface tension).The mathematical expression of the functional in dimensionless form is: where Ω is the considered domain and ℎ is the Cahn number, which represents the dimensionless thickness of the thin interfacial layer between the two fluids: where  * is clearly the dimensional thickness of the interfacial layer.From equation (2.3), the expression of the chemical potential can be derived as the functional derivative with respect to the order parameter: (2.5) At equilibrium, the chemical potential is constant throughout all the domain.The equilibrium profile for a flat interface can thus be obtained by solving ∇  = 0, hence: where  is the coordinate normal to the interface.As anticipated before, the last term in the right-hand side of the Cahn-Hilliard equation (equation 2.1) is a penalty-flux term employed in the profile-corrected formulation of the phase-field method, and is used to overcome some potential drawbacks of the standard formulation of the method, e.g.mass leakages among the phases and misrepresentation of the interfacial profile (Yue et al. 2007;Li et al. 2016).This penalty flux is defined as: where  = 0.0625/ℎ (Soligo et al. 2019c).

Hydrodynamics
To describe the hydrodynamics of the multiphase system, the Cahn-Hilliard equation is coupled with the Navier-Stokes equations.The presence of a deformable interface (and of the corresponding surface tension forces) is accounted for by introducing an interfacial term in the Navier-Stokes equations.Recalling that in the present case we consider two fluids with the same density ( * =  *  =  *  ) and viscosity ( * =  *  =  *  ), continuity and Navier-Stokes equations in dimensionless form read as: is the pressure field, while T c is the Korteweg tensor (Korteweg 1901) used to account for the surface tension forces and defined where I is the identity matrix and ⊗ represents the dyadic product.This approach is the continuum surface stress approach (Lafaurie et al. 1994;Gueyffier et al. 1999) applied in the context of PFM, and is analytically equivalent to the chemical potential forcing (Mirjalili et al. 2023).The dimensionless groups appearing in the Navier-Stokes equations are the shear Reynolds number   (ratio between inertial and viscous forces), and the Weber number Focus on Fluids articles must not exceed this page length   (ratio between inertial and surface tension forces), which are defined as: (2.11) where  * is the surface tension.Note that, consistently with the employed adimensionalization,   is defined using the half-channel height (and not the drop diameter).

Energy equation
The time evolution of the temperature field is obtained by solving the energy equation using a one-scalar model approach (Zheng et al. 2015).To avoid the introduction of further complexity in the system, we consider two fluids with the same thermophysical properties, i.e. same thermal conductivity  * , same specific heat capacity  *  and therefore same thermal diffusivity  * =  * / *  *  (since the density of the two phases is also the same).These properties have been evaluated at a reference temperature  *  = ( * ,0 +  * ,0 )/2, i.e. the average between the initial drops temperature and the carrier fluid temperature, and are assumed to be constant and uniform.Within these assumptions, the energy equation written in dimensionless form reads as: where  is the Prandtl number defined as with  * =  * / * the kinematic viscosity (i.e.momentum diffusivity).From a physical viewpoint,  represents the momentum-to-thermal diffusivity ratio.

Numerical discretization
The governing equations (2.1), (2.8), (2.9) and (2.12) are solved using a pseudo-spectral method, which uses Fourier series along the periodic directions (streamwise and spanwise) and Chebyshev polynomials along the wall-normal direction.The Navier-Stokes and continuity equations are solved using a wall-normal velocity-vorticity formulation: equation (2.9) is rewritten as a 4  ℎ order equation for the wall-normal component of the velocity  and a 2  order equation for the wall-normal component of the vorticity   (Kim et al. 1987;Speziale 1987).The Cahn-Hilliard equation (2.1), which in its original form is a 4  ℎ order equation is split into two 2  order equations using the splitting scheme proposed in Badalassi et al. (2003).Using this scheme, the governing equations are recasted as a coupled system of Helmholtz equations, which can be readily solved.
The governing equations are time advanced using an implicit-explicit scheme.For the Navier-Stokes, the linear part is integrated using a Crank-Nicolson implicit scheme, while the non-linear part is integrated using an Adams-Bashforth explicit scheme.Similarly, for the Cahn-Hilliard and energy equations, the linear terms are integrated using an implicit Euler scheme, while the non-linear terms are integrated in time using an Adams-Bashforth scheme.The adoption of the implicit Euler scheme for the Cahn-Hilliard equation helps to damp unphysical high-frequency oscillations that could arise from the steep gradients of the phase field.
As the characteristic length scales of the flow and temperature fields, represented by the Kolmogorov scale,  +  , and the Batchelor scale,  +  , are different when non-unitary Prandtl numbers are employed (being these two quantities linked by the following relation , a dual grid approach is employed to reduce the computational cost of the simulations and at the same time to fulfill the DNS requirements.In particular, when super-unitary Prandtl numbers are simulated, a finer grid is used to resolve the energy equation.Spectral interpolation is used to upscale/downscale the fields from the coarse to the refined grid and vice-versa when required (e.g.upscaling of the velocity field to compute the advection terms in the energy equation).
This numerical scheme has been implemented in a parallel Fortran 2003 MPI in-house proprietary code.The parallelization strategy is based on a 2D domain decomposition to divide the workload among all the MPI tasks.The solver execution is accelerated using openACC directives and CUDA Fortran instructions (solver execution) while the Nvidia cuFFT libraries are used to accelerate the execution of the Fourier/Chebyshev transforms.Overall, the computational method adopted allows for the accurate resolution of all the governing equations and the achievement of an excellent parallel efficiency thanks to the fine-grain parallelism offered by the numerical method used.The equivalent computational cost of the simulations is about 25 million CPU hours and the resulting dataset has a size of about 16 TB.

Boundary conditions
The system of governing equations is complemented by a set of suitable boundary conditions.For the Navier-Stokes equations, no-slip boundary conditions are enforced at the top and bottom walls (located at  = ±ℎ): (2.14) For the Cahn-Hilliard equation, no-flux boundary conditions are applied at the two walls, yielding to the following boundary conditions: (2.15) Likewise, for the energy equation, no-flux boundary conditions are applied at the two walls (i.e.adiabatic walls).  ( = ±ℎ) = 0 . (2.16) Along the streamwise and spanwise directions ( and ), periodic boundary conditions are imposed for all variables (Fourier discretization).The adoption of these boundary conditions leads to the conservation of the phase field and temperature fields over time: (2.17 where Ω is the computational domain.Regarding the phase-field, equation (2.17) enforces mass conservation of the entire system but does not guarantee the conservation of the mass of each phase (Yue et al. 2007;Soligo et al. 2019c), as some leakages between the phases may occur.This drawback is rooted in the phase-field method and more specifically in the curvature-driven flux produced by the chemical potential gradients (Kwakkel et al. 2020;Mirjalili & Mani 2021).This issue is here successfully mitigated with the adoption of the profile-corrected formulation that largely reduces this phenomenon.In the present cases, mass leakage between the phases occurs only in the initial transient when the phase field is initialized (see the section below for details on the initial condition) and is limited to 2% of the initial mass of the drops.After this initial transient, the mass of each phase keeps constant.
Case Table 1: Overview of the simulation parameters.For a fixed shear Reynolds number   = 300 and Weber number   = 3, we consider a single-phase flow case and four non-isothermal drop-laden flows characterised by different Prandtl numbers: from  = 1 to  = 8.The grid resolution is modified accordingly so as to satisfy DNS requirements.

Simulation set-up
The turbulent channel flow, driven by an imposed constant pressure gradient in the streamwise direction, has shear Reynolds number   = 300.The computational domain has dimensions × 600 wall units.The value of the Weber number is kept constant and is equal to   = 3.00, so to be representative of liquid/liquid mixtures (Than et al. 1988).To study the influence of the Prandtl number  on the heat transfer process, we consider four different values of :  = 1,  = 2,  = 4 and  = 8.These values cover a wide range of real-case scenarios: from low-Prandtl number fluids to water-toluene mixtures.
The grid resolution used to resolve the continuity, Navier-Stokes and Cahn-Hilliard equations is equal to   ×   ×   = 1024 × 512 × 513 for all the cases considered in this work.For the energy equation, the same grid used for the flow field and phase field is employed at the lower Prandtl numbers ( = 1 and  = 2), while a more refined grid, with   ×   ×   = 2048 × 1024 × 513 points, is used when the larger Prandtl numbers are considered ( = 4 and  = 8).The computational grid has uniform spacing in the homogeneous directions, while Chebyshev-Gauss-Lobatto points are used in the wallnormal direction.We refer the reader to table 1 for an overview of the main phzsical and computational parameters of the simulation.For the employed grid resolution, the Cahn number is set to ℎ = 0.01 while, to achieve convergence to the sharp interface limit, the corresponding phase field Péclet number is Pe = 1/ℎ = 50.
All simulations are initialized releasing a regular array of 256 spherical drops with diameter  = 0.4ℎ (corresponding to  + = 120 ..) inside a fully-developed turbulent flow field (obtained from a preliminary simulation).To ensure the independence of the results from the initial flow field condition, each case is initialized with a slightly different flow field realization.Naturally, the fields are equivalent in terms of statistics as they are all obtained from a statistically steady turbulent channel flow.The volume fraction of the drops is Φ =   /(  +   ) = 5.4%, with   and   the volume of the drops and carrier fluid, respectively.
The initial condition for the temperature field is such that all drops are initially warm (initial temperature  ,0 = 1), while the carrier fluid is initially cold (initial temperature  ,0 = 0).To avoid numerical instabilities that might arise from a discontinuous temperature field, the transition between drops and carrier fluid is initially smoothed using a hyperbolic tangent kernel.Figure 1 (which is an instantaneous snapshot captured at  + = 1000, for  = 1) shows a volume rendering of the temperature field (blue-cold, red-hot), inside which deformable drops (whose interface, iso-level  = 0, is shown in white) are transported.

Results
Results obtained from the numerical simulations will be first discussed from a qualitative viewpoint, by looking at instantaneous flow and drops visualizations, and then analyzed from a more quantitative viewpoint, by looking at the drop size distribution (DSD), and at the effect of the Prandtl number  on the average drops and fluid temperature.To explain the numerical results, and to offer a possible parametrization of the heat transfer process in drop-laden flows, we will also develop a simplified phenomenological model of the system.Finally, we will characterize the temperature distribution inside the drops, elucidating the effects of  and of the drop size on it.Note that, unless differently mentioned, results are presented using the wall-unit scaling system but for the temperature field, which is made dimensionless using the initial temperature difference as a reference scale (which is a natural choice in the present case).

Qualitative discussion
The complex dynamics of drops immersed in a non-isothermal turbulent flow is visualized in figure 1, where the drops (identified by isocontour of  = 0) are shown together with a volume-rendered distribution of temperature in the carrier fluid.Also shown in figure 1 is a close-up view of the temperature distribution inside the drop.
Once injected into the flow, each drop starts interacting with the flow and with the neighbouring drops.The result of the drop-turbulence and drop-drop interactions is the occurrence of breakage and coalescence events.A breakage event happens when the flow vigorously stretches the drop, leading to the formation of a thin ligament that breaks and generates two child drops.Upon separation, surface tension forces tend to retract the broken filaments and to restore the drop spherical shape.A coalescence event is observed when two drops come close each other.The small liquid film that separates the drops start to drain, and a coalescence bridge is formed.Later, surface tension forces enter the picture, reshaping the drop and completing the coalescence process.The dynamic competition between breakage and coalescence events, and their interaction with the turbulent flow, determines the number of drops, their size distribution, and their shape/morphology (i.e., curvature, interfacial area, etc.).
In the present case, drops not only exchange momentum with the flow and with the other drops, but also heat.Starting from an initial condition characterized by warm drops (with uniform temperature) and cold carrier fluid, and because of the imposed adiabatic boundary conditions, the system evolves towards an equilibrium isothermal state.During the transient to attain this thermal equilibrium state, heat is transported by diffusion and advection inside each of the two phases, and across the interface of the drops (see the temperature field inside and outside the drops, figure 1).The picture is then further complicated by the occurrence of breakage and coalescence events.This is represented in figure 2. When a breakage occurs (figure 2, top row), a thin filament is generated (figure 2a-c), which then leads to the formation of a smaller satellite drop (figure 2d-e).The filament and the satellite drop, given the large surface-to-volume ratio, exchange heat very efficiently, and become rapidly colder.In contrast, when a coalescence occurs (figure 2, bottom row), two drops having different temperature merge together.This induces an efficient mixing process, during which cold parcels of one drop become warmer and viceversa, warm parcles of the other drops become colder.Overall, breakup and coalescence events induce heat transfer modifications that are in general hard to predict a priori, since they do depend on the relative size of the involved parents/child drops.
Naturally, the problem of heat transfer in drop-laden turbulence is strongly influenced by the Prandtl number of the flow.This can be appreciated by looking at figure 3, where we show the instantaneous temperature field, together with the shape of the drops, at a certain instant in time ( + = 1500) and at the different Prandtl numbers:  = 1 (panel ),  = 2 (panel ),  = 4 (panel ) and  = 8 (panel ).In each panel, the temperature field is shown on a wall-parallel  + - + plane located at the channel center ( + = 0), and is visualized with a blue-red scale (blue-low, red-high).We observe that the temperature field changes significantly with .In particular, we notice an increase in the drop-to-fluid temperature difference for increasing , going from  = 1 (top panel) where this difference is small, to  = 8 (bottom panel) where this difference is large.The heat transfer from the drops to the carrier fluid becomes slower as  increases, consistently with a physical situation in which the  number is increased by reducing the thermal diffusivity of the fluid, while keeping the momentum diffusivity constant (i.e.constant kinematic viscosity, and hence shear Reynolds number).Also, the temperature structures, both inside and outside the drop, become thinner and more complicated at higher , since their characteristic lengthscale, the Bachelor scale  +  ∝  −1/2 , becomes smaller for increasing  (Batchelor 1959;Batchelor et al. 1959).In addition, smaller drops have, on average, a lower temperature compared to larger drops, regardless of the value of .All these aspects will be discussed in more detail in the next sections.

Drop Size Distribution
To characterize the collective dynamics of the drops, we compute the drop size distribution (DSD) at steady-state conditions, averaging over a time window Δ + = 3000, from  + = 3000 to 6000.It is worth mentioning that a quasi-equilibrium DSD, very close to the steady one, is already achieved around  + ≃ 750, and only minor changes occur to the DSD afterwards.
Figure 4 shows the DSD obtained for the different cases considered here:  = 1 (dark violet),  = 2 (violet),  = 4 (pink), and  = 8 (light pink).The distributions have been computed considering, for each drop, the diameter of the equivalent sphere computed as: where  + is the volume of the drop.Also reported in figure 4 is the Kolmogorov-Hinze scale,  +  , which can be computed as (Perlekar et al. 2012;Roccon et al. 2017;Soligo et al. 2019a): where   is the turbulent dissipation, here evaluated at the channel center where most of the drops collect because of their deformability (Lu & Tryggvason 2007;Soligo et al. 2020;Mangani et al. 2022).The Kolmogorov-Hinze scale identifies the critical diameter below which drop breakage is unlikely to occur (Kolmogorov 1941;Hinze 1955).Separated by the Kolmogorov-Hinze scale, we observe the emergence of two different regimes.For drops smaller than the Kolmogorov-Hinze scale, we find the coalescence-dominated regime (left, gray area), in which drops, which are smaller than the critical scale, are generally not prone to break.For drops larger than the Kolmogorov-Hinze scale, we find the breakup-dominated F. Mangani, A. Roccon, F. Zonta, and A. Soldati Kolmogorov-Hinze scale We note that, for equivalent diameters above the Hinze scale, our results follow quite well the theoretical scaling law and match the drops/bubbles size distributions obtained in literature considering similar flow instances (Deike et al. 2016;Di Giorgio et al. 2022;Soligo et al. 2021;Deike 2022;Crialesi-Esposito et al. 2023).Below the Hinze scale, for equivalent diameters in the range 25 <  +  <  +  our results match reasonably well the theoretical scaling law.For equivalent diameters  +  < 25 .., we observe an underestimation of the DSD compared to the proposed scaling.This is linked to the grid resolution, and in particular to the problem in describing very small drops (Soligo et al. 2021).

Mean temperature of drops and carrier fluid
We now focus on the average temperature of the drops and of the carrier fluid.We consider the ensemble of all drops as one phase, and the carrier fluid as the other phase (using the value of the phase field as a phase discriminator), and we compute the average temperature for each phase.The evolution in time of the drops and carrier fluid temperature,   and   respectively, is shown in figure 5, for the different values of .Together with the results obtained by current direct numerical simulations, filled symbols, in figure 5 we also show the predictions obtained by a simplified phenomenological model (solid lines), the details of which will be described and discussed later (see section 3.4).We start considering the DNS results only.As expected, we observe that the average temperature of the drops (violet to pink symbols) decreases in time, while the average temperature of the carrier fluid (blue to cyan symbols) increases in time, until the thermodynamic equilibrium, at which both phases have the same temperature, is asymptotically reached.For this reason, simulations have been run long enough for the average temperature of both phases to be sufficiently close to the equilibrium temperature.In particular, we stopped the simulations at  + ≃ 6000, when the condition with  ,0 the initial temperature of the drops, is satisfied by all simulations.The equilibrium temperature,   , can be easily estimated a-priori: since the two walls are adiabatic, and the homogeneous directions periodic, the energy of the system is conserved over time.The energy balance can be written as: where  *  and  *  is the mass of the carrier fluid and of the drops;  * ,0 and  * ,0 the physical value of the initial temperatures of the two phases and  *  the specific heat capacity.Considering that the two phases have equal density and specific heat capacity, we obtain: Recalling the definition of volume fraction, Φ =   /(  +   ), and making the equation dimensionless, we finally get: which is represented by the horizontal dashed line in Figure 5. Figure 5 provides also a clear indication that the higher the Prandtl number, the larger the time it takes for the system to reach the equilibrium temperature,   .The trend can be observed for both the drops and carrier fluid, as the two phases are mutually coupled (the heat released from the drops is adsorbed by the carrier fluid).This result confirms our previous qualitative observations, see figure 3 and discussion therein, that a large  (small thermal diffusivity) reduces the heat released by the drops.It is also interesting to observe that the behavior of the mean temperature of the two phases appears self-similar at the different .

A phenomenological model for heat transfer rates in droplet laden flows
In an effort to provide a possible interpretation of the previous results -and in particular to explain the average temperature behavior shown in figure 5 -, we develop a simple physicallysound model of the heat transfer process in droplet-laden turbulence.We start by considering the heat transfer mechanisms from a single drop of diameter  * to the surrounding fluid: t + P r = 1 P r = 2 P r = 4 P r = 8 P r = 1 P r = 2 P r = 4 P r = 8 Final equilibrium temperature Model (Eq.3.18) Figure 5: Time evolution of the mean temperature of drops (violet to pink colors, different symbols) and carrier fluid (blue to cyan colors, different symbols) for the different Prandtl numbers considered.DNS results are reported with full circles while the predictions obtained from the model are reported with continuous lines.The equilibrium temperature of the system,   , is reported with a horizontal dashed line.In general, it can be observed that by increasing the Prandtl number (corresponding to a decrease of the thermal diffusivity), the heat transfer process between the two phases slows down and more time is required to achieve the equilibrium condition.
Reportedly (Schlichting & Gersten 2016, 218), the thermal boundary layer thickness,  *  , can be expressed as  *  =  *  −  where  * is the momentum boundary layer thickness, and  is an exponent that depends on the flow condition in the proximity of the boundary where the boundary layer evolves.In particular, the exponent  ranges from  = 1/3 for no slip conditions, usually assumed for solid particles, to  = 1/2, usually assumed for clean gas bubbles.For an in-depth discussion on the topic, we refer the reader to appendix A. As a consequence, the heat transfer rate observed from drops/bubbles is expected to be larger than that observed from solid particles, since the no-slip boundary condition generally weakens the flow motion near the interface (Levich 1962;Herlina & Wissink 2016;Bird et al. 2002).We can now rewrite the equation of the model in dimensionless form, using the initial dropto-carrier fluid temperature difference Δ * =  * ,0 −  * ,0 as reference temperature, and  * / * 2  as reference time: where  + is the drop diameter in wall units, while   =  *   * / * is the Reynolds number based on the boundary layer thickness (which can be assumed constant among the different cases).Equation 3.10 can be rewritten as: where C is a constant whose value depends only on the flow structure, i.e. on Re  .Equation 3.11 describes the heat released by a single drop of dimensionless diameter  + .Assuming now that the turbulent flow is laden with drops of different diameters the general equation describing the heat released by the -th drop of diameter  +  becomes: where F  is the lumped-parameters representation of the right-hand side of the temperature evolution equation for the -th drop.As widely observed in literature (Deane & Stokes 2002;Soligo et al. 2019c), and also confirmed by the present study (figure 4), we can hypothesize an equilibrium drops-size-distribution (DSD) by which the number density of drops scales as  + −3/2 , in the sub-Hinze range of diameters (10 <  + < 110), and as  + −10/3 in the super-Hinze range of diameters (110 <  + < 240).With this assumption, and considering 7 classes of drops diameter for the sub-Hinze range, and 4 classes for the super-Hinze range, we can integrate equation 3.12 to obtain the time evolution of the temperature of each drop in time: (3.13) From a weighted averaged of the temperature (based on the number of drops in each class, as per the theoretical DSD), we obtain the average temperature of the drops,   .
To obtain the mean temperature of the carrier fluid, we consider that (adiabatic condition at the walls), the heat released by the drops is entirely absorbed by the carrier fluid.The heat released by the drops with a certain diameter  *  , can be computed as: where  *  () is the number of drops for that specific diameter (as per the DSD).The overall heat released by all drops can be calculated as the summation over all different classes of diameters: where   is the employed number of classes.Finally, the mean temperature of the carrier fluid is In dimensionless form (dividing by the initial drop-to-carrier fluid temperature Δ * ) equation 3.16 becomes: (3.17)The results of the model are shown in figure 5. Interstingly, under the simplified hypothesis of the model (chiefly, spherical shape of the drops, constant drop-size-distribution evaluated at the equilibrium), we observe that the behavior of the mean temperature is very well captured by the model (represented by the solid lines in figure 5) i.e. when  = 1/3 -typical of boundary layers around solid objects (i.e.solid particles).
Reasons for this might be the presence of wakes/sheltering effects between drops, but also the fact that drops are strongly advected by the mean flow, and the flow condition at the drop surface can be different from the slip one, and is in general not of simple evaluation.Given the relationship   / ∼  −2/3 postulated by the model (equation 3.18), which provides results in very good agreement with the numerical ones, it seems reasonable to rescale the F. Mangani, A. Roccon, F. Zonta, and A. Soldati t + = t + /P r 2/3 P r = 1 P r = 2 P r = 4 P r = 8 P r = 1 P r = 2 P r = 4 P r = 8 Final equilibrium temperature Figure 6: Time evolution of the mean temperature of drops (violet to pink colors) and carrier fluid (blue to cyan colors) for the different Prandtl numbers considered obtained from DNS and reported against the dimensionless time t+ =  + / 2/3 .The equilibrium temperature of the system,   , is reported with a horizontal dashed line.The DNS results reported over the new dimensionless time nicely collapse on top of each other, highlighting the self-similarity of the  , profiles.
time variable as: A representation of the DNS results in terms of the rescaled time, equation 3.19, is shown in figure 6.We observe a nice collapse of the two sets of curves -drops and carrier fluid (red and blue) -for the different values of , which clearly demonstrates the self-similar behavior of .For this reason, the rescaling of time t+ =  + / 2/3 , will be also used in the following.

Heat transfer from particles and drops/bubbles
It is now important to discuss the behavior of the heat transfer coefficient (and its dimensionless counterpart, the Nusselt number ), also in the context of available literature results.Naturally, similar considerations can be made to evaluate the mass transfer coefficient, in particular at liquid/gas interfaces (Levich 1962;Bird et al. 2002).
For solid particles, a balance between the convective time scale near the surface, and the diffusion time scale, gives an heat transfer coefficient (Krishnamurthy & Subramanian 2018): (3.20) and the corresponding Nusselt number: where  is an exponent that depends on the flow conditions and links the boundary layer thickness to the particle Reynolds number.Usually,  = 1/3 for small Reynolds numbers (Krishnamurthy & Subramanian 2018) while  = 1/2 for large Reynolds numbers (Ranz 1952;Whitaker 1972;Michaelides 2003).
Using similar arguments (balance between convective and diffusion time scales), but considering now that at the surface of a drop/bubble a slip velocity, and therefore a certain degree of advection, can be observed (Levich 1962;Bird et al. 2002;Herlina & Wissink 2016), the heat transfer coefficient is found to scale as: (3.22) and the corresponding Nusselt number as: where also in this case the exponent  does depend on the considered Reynolds number.Two regimes are usually defined (Theofanous et al. 1976): a low Reynolds number regime, for which  = 1/2, and a high Reynolds number regime, for which  = 3/4.An alternative approach, which gives similar predictions, is to use the penetration theory of Higbie (1935), in which turbulent fluctuations are invoked to estimate a flow exposure (or contact) time, to compute the heat/mass transfer coefficient.Such approach has been widely used in bubbleladen flows (Colombet et al. 2011;Herlina & Wissink 2014, 2016;Farsoiya et al. 2021).
We can now evaluate the heat transfer coefficient from our DNS at different , and compare it to the proposed scaling laws.Note that the heat transfer coefficient is obtained as: where the numerator represents the temperature difference of the drops between the time steps  and  + 1, while the denominator represents the temperature difference between the drop and the carrier fluid evaluated halfway in time between step  and  + 1 (i.e. at  + 1/2).The quantity  is the total interfacial area between drops and carrier fluid, while Δ is the time step used to evaluate the heat transfer.Here, we have evaluated the heat transfer coefficient taking the heat released by the drops as reference; an equivalent result, but with opposite sign, can be obtained using the heat absorbed by the carrier fluid as reference, and taking into account the different volume fraction of the two phases.
The dimensionless heat transfer coefficient, equation 3.24, is shown as a function of , and at different time instants, in figure 7.For a better comparison, the results are normalized by the value of the heat transfer coefficient for  = 1.The two reference scaling laws, H ∼  −2/3 obtained for  = 1/3 and H ∼  −1/2 obtained for  = 1/2 are also shown by a dotted and a dashed line.We note that at the beginning of the simulations (see for example  + = 250), the heat transfer coefficient is close to H ∼  −2/3 , while at later times it tends towards H ∼  −1/2 , hence approaching the scaling law proposed for heat/mass transfer in gas-liquid flows (Levich 1962;Magnaudet & Eames 2000;Bird et al. 2002;Herlina & Wissink 2014, 2016;Colombet et al. 2018;Farsoiya et al. 2021).
A possible explanation is that, as time advances, the shape of the drops becomes complex, and coalescence/breakups more frequent, thus inducing an overall surface decrease that is associated to an heat transfer increase.This reflects into an heat transfer process that is slower at the beginning, H * ∼  −2/3 , and faster at later times H * ∼  −1/2 .

Influence of the drop size on the average drop temperature
In the previous sections, we have studied the behavior of the mean temperature field of the drops and of the carrier fluid considered as single entities.However, while this description is perfectly reasonable for the carrier fluid -which can be considered a continuum -it can be questionable for the drops, which are not a continuum phase by nature.We now take the dispersed nature of the drops into account and we evaluate, for each drop, the equivalent diameter and the corresponding mean temperature.This is sketched in figure 8, where the average temperature of each drop (represented by a dot) is shown as a function of its equivalent diameter, at different time instants (between  + = 1050 and  + = 2400).Each panel refers to a different Prandtl number.Note that, at  + = 2400, the case  = 1 has almost reached the thermodynamic equilibrium (figure 5).It is clearly visible that, regardless of the considered time, small drops have an average temperature close to the equilibrium one.This is particularly visible at smaller Prandtl numbers, i.e. when heat transport is faster, but it can be observed also at larger .In contrast, the average temperature of larger drops is larger.Hence, the average temperature of the drops seems directly proportional to the drop size, as can be argued considering that the heat released by the drop, and hence its temperature reduction, is   / ∝  −1 (equation 3.16).It is therefore not surprising that the scatter plot at a given time instant is characterized by dots distributing in stripes-like fashion, with slope that decreases with time.This behavior is observed at all , although the range of drops temperature (y axis) at small  is definitely narrower (because of their larger heat loss) compared to that at large .It is also interesting to note -in particular at  = 4 and  = 8 (panels  and ) -the presence of drops with a temperature smaller than the equilibrium one (dots falling below the horizontal line that marks the equilibrium temperature).We can link this behavior to the small relaxation time of small drops that therefore adapt quickly to the local temperature of the carrier fluid, which can be smaller than the equilibrium one for two main reasons.First, at the early stages of the simulations, and at high Prandtl numbers, the temperature of the carrier fluid is lower than the equilibrium one.Second, temperature fluctuations (of both negative and positive sign) are present also in the carrier fluid.These fluctuations, in the form of hot/cold striations are more likely observed at large  (see the striation-like structures at  = 8 in figure 3𝑑).3.7.Temperature distribution inside the drops In many applications, in particular to evaluate mixing efficiency and flow homogeneity, not only the average temperature of drops is important, but also its space and time distribution inside the drops.To understand it, we now look at the PDF of the temperature fluctuations inside the drops, where   is the local temperature inside the drop, and   is the average temperature of all drops at a certain time (as per figure 5).Results are shown in figure 9.The first row of figure 9 shows the probability density function of  ′  at different , and at two different time instants:  + = 600 (panel , left) and  + = 1500 (panel , right).The second row of figure 9 shows the PDFs obtained at two different rescaled time instants, t+ =  + / 2/3 : t+ = 600 (panel , left) and t+ = 1500 (panel , right).
Considering first figure 9 ( + = 600), we notice that all PDFs have a rather regular shape, characterized by the presence of both positive and negative fluctuations (with respect to the average temperature), with a slight asymmetry towards the positive ones (positive fluctuations are more likely observed).A comparison between the curves obtained at different  shows that the range of temperature fluctuations is wider at larger .This is due to the small thermal diffusivity at large , which allows temperature fluctuations in the drop to survive much longer before they are damped and spread by diffusion.Naturally, at later times (figure 9𝑏,  + = 1500), the range of temperature fluctuations reduces.Indeed, as heat is transferred from the drops to the carrier fluid, the maximum temperature of drops reduces, and so does the range of temperature fluctuations inside the drop.This trend is more pronounced for negative fluctuations, as the minimum temperature inside the drops is somehow bounded by the temperature of the carrier fluid (which increases only a little, from  ,0 = 0 to   = 0.054, during the simulation).This latter observation is visible in the shape of the PDFs at  = 1, 2 and 4, since the system is closer to the thermal equilibrium at this time instant (the thermal equilibrium is identified in panel  by a vertical dashed line and marked with a label,    ): a sharp drop of the PDF, which does not significantly trespass the    limit, is observed.In contrast, positive temperature fluctuations are subject to relatively weaker constraints (they are only bounded by the maximum initial temperature of the drops).This results into a PDF that gets asymmetric, positively-skewed.It is also interesting to observe the development of a pronounced peak about the equilibrium temperature    , which corresponds to the presence of small drops (generated by breakages events) that -given their small thermal relaxation time and heat capacity -almost immediately adapt to the equilibrium temperature (see also figure 2𝑑,  ).However, a discussion on the temperature fluctuations, captured from flows at different  and after the same time  + from the initial condition, could be misleading, because it puts in contrast flows at different thermal states (i.e.different average temperature, and different temperature gradients, see figure 5).To filter out this effect, we compute the PDFs of the temperature fluctuations at the same rescaled time instants t+ =  + / 2/3 .By doing this, all cases can be considered at similar thermal conditions (see also figure 6).The resulting PDFs, at t+ = 600 and t+ = 1500, are shown in figure 9-.Note that, for the sake of clarity, the corresponding  + , which is different from case to case, is reported between brackets in the legend.In the rescaled time units, the collapse between the different curves is quite nice.The slight difference between the curves is due to the fact that, although the system is at the same thermal state (same t+ ), it is at different flow state (different  + ), i.e. the instantaneous drop size distributions are different.This gives the slightly larger negative fluctuations at larger  (which, being at a later stage, is characterized by the presence of smaller and colder drops), and slightly larger positive fluctuations at smaller  (which, being at an earlier flow state, is characterized by the presence of larger and warmer drops).
From a closer look at figure 9 (t + = 1500), we note very clearly the constraint set by the thermal equilibrium condition: the PDF cannot significantly trespass the limit represented by   (vertical dashed line), which is very similar for all , given the similar thermal state.Also visible is the peak, already discussed in figure 9𝑏, that emerges very close to the equilibrium temperature   , and which is due to the presence of small drops that adapt quickly to the local temperature of the carrier fluid.As previously noticed in figure 9𝑐, the higher probability of finding small drops at lower  is also responsible for the narrowing of the PDF (reduction of positive temperature fluctuations).

Conclusions
In this work, we studied heat transfer in a turbulent channel flow laden with large and deformable drops.The drops are initially warmer than the carrier fluid and as the simulations advance, heat is transferred from the drops to the carrier fluid.Simulations considered a fixed value of the Reynolds number, Re  = 300, and Weber number,   = 3, and analyze different Prandtl number values, from  = 1 to  = 8.The Prandtl number is changed by changing the thermal diffusivity.The investigation is based on the direct numerical simulation of turbulent heat transfer, coupled with a phase-field method, used to describe interfacial phenomena.First, we focused on the drops dynamics, observing that after an initial transient (up to  + = 1000), the drop size distribution (DSD) reaches a quasi-equilibrium condition where it follows the scaling  +  −3/2 in the coalescence-dominated regime and  +  −10/3 in the breakage-dominated regime.The threshold between the coalescence-dominated and the breakage-dominate regimes is represented by the Kolmogorov-Hinze scale.Then, we characterize the behavior of the average temperature of the drops and of the carrier fluid: as expected, the average temperature of drops decreases in time, while the average temperature of the carrier fluid increases in time, until reaching the equilibrium condition of uniform temperature in the whole system.We clearly observed that the higher the Prandtl number, the larger the time it takes for the system to reach the equilibrium temperature.Interestingly, the time behavior of the temperature profiles of both drops and carrier fluid is self-similar.Building on top of these numerical results, we developed a phenomenological model that can accurately reproduce the time evolution of the mean temperatures at all Prandtl numbers considered here.This model gave us the opportunity to introduce a new self-similarity variable (time, t+ ) that accounts for the Prandtl number effect, and by which all results collapse on a single curve.In addition, we also computed the heat transfer coefficient H (and its dimensionless counterpart, the Nusselt number ) and showed that it scales as H ∼  −2/3 (which corresponds to a Nusselt number scaling  ∼  1/3 ) at the beginning of the simulation, and tends to H ∼  −1/2 (or, alternatively,  ∼  1/2 ) at later times.These different scalings are consistent with previous literature predictions, and can be explained via the boundary layer theory (appendix A).The effects of the Prandtl number on the temperature distribution inside the drops has been investigated.We observe that by increasing the Prandtl number, the PDFs become wider and thus large temperature fluctuations are more likely to be observed.Interestingly, when the PDFs are compared at a the same rescaled time t+ (i.e.accounting for the Prandtl number effect), all curves collapse on top of each other, with only minor differences possibly due to the different instantaneous drop size distribution.The effect of the drop size was also discussed: small drops adapt faster to the equilibrium temperature, thanks to their small heat capacity, compared to larger drops.Finally, it must be pointed out that, since the different phases of a multiphase flow can have different thermophysical properties, Prandtl numbers can be also different from phase to phase.This aspect, which was not considered in the present work, will be the topic of a future study.In addition, in the present work we have assumed a constant and uniform surface tension.However, in many circumstances, surface tension does depend on temperature, therefore inducing thermocapillary effects.This will be also the subject of a future investigation.

Figure 1 :
Figure1: Rendering of the computational setup employed for the simulations.A swarm of large and deformable drops is released in a turbulent channel flow.The flow goes from left to right driven by a constant pressure gradient.The temperature field is volume-rendered (blue-low, red-high) and the drops interface is shown in white (iso-level  = 0).As can be appreciated from the close-up view (top left, which shows the temperature field in a drop section), drops have a temperature higher than the carrier fluid and as a result, there is a heat flux from the drops to carrier fluid.The snapshot refers to  = 1 and  + = 1000.
Figure2: Influence of topology changes on heat transfer: Time sequence of a breakage event (top row, panels -) and of a coalescence event (bottom row, panels  -).In the top row, we can appreciate how during a breakage event in the pinch-off region heat is efficiently transferred from the drops to carrier fluid thanks to the high surface/volume ratio of this region.Likewise, the thermal relaxation time of the small drops generated during the breakage (bottom) is smaller than that of the parent drop (top).In the bottom row, we can appreciate how during a coalescence event, parcels of fluid with different temperatures mix.The two time sequences refer to the case  = 1 and snapshots are separated by Δ + = 15.

Figure 3 :
Figure 3: Instantaneous visualization of the temperature field (red-hot; blue-cold) on a  + −  + plane located at the channel center for  + = 1500.Drops interface (iso-level  = 0) are reported using white lines.Each panel refers to a different Prandtl number:  = 1 (panel ),  = 2 (panel ),  = 4 (panel ) and  = 8 (panel ).By increasing the Prandtl number (from top to bottom) and thus decreasing the thermal diffusivity, the heat transfer process slows down.As a consequence, the drop temperature is higher when larger Prandtl numbers are considered.The effects of the Prandtl number on the characteristic length scales can also be appreciated: as it increases, temperature structures become smaller.

Figure 8 :
Figure 8: Scatter plot of the drop equivalent diameter  + against the drop average temperature   .Each dot represents a different drop while its color (black to gray colormap) identifies different times, from  + = 1050 (black) up to  + = 2400 (gray).Each panel refers to a different Prandtl number.A sketch showing drops of different equivalent diameters is reported in the upper part of panel .As can be appreciated, larger drops have larger relaxation times, and thus a longer transient is required to reach thermal equilibrium.This effect becomes more pronounced as the Prandtl number is increased and thus the overall heat transfer process slows down (smaller heat transfer coefficient).

Figure 9 :
Figure 9: Probability density function (PDF) of the temperature fluctuations,  ′  =   −   inside the drops.Each case is reported with a different color (violet to light pink) depending on the Prandtl number.The first row shows the PDFs obtained at two different time instants:  + = 600 (panel ) and  + = 1500 (panel ).The second row shows the PDFs obtained at two rescaled time instants t+ = 600 (panel ) and t+ = 1500 (panel ), where the rescaled time is computed as t+ =  + / 2/3 .For panels -, the corresponding  + is reported between brackets.
are the mass, external surface, and specific heat of the drop, H * is the heat transfer coefficient, while  *  and  *  are the drops and carrier fluid temperature.The heat transfer coefficient can be estimated as the ratio between the thermal conductivity of the external fluid,  * , and a reference length scale, here represented by the thermal boundary layer thickness  * : H * ∼  * / *  .(3.8)With this assumption, and recalling that  * =  *  =  *  , equation 3.7 becomes:  *   * = 6   *  *  * Figure7: Time behavior of the dimensionless heat transfer coefficient for the different Prandtl numbers considered.Heat transfer coefficients are reported normalized by the value of the heat transfer coefficient obtained for  = 1 (at the same time instant).In this way, results obtained at different time instants can be conveniently compared.The two scaling laws that refer to  = 2/3 and  = 1/2 are also reported as references.