4.1. Heat-transfer mechanisms and temperature field in the decoupled regime

Figure 3.

Effect of the viscosity ratio $m$ on (a) $D_{2}$ (in the relative reference frame moving at the mean speed of the annulus, $V=U_{2}$ ) and (b) the source term related to viscous dissipation $W_{2}^{\star }/\beta$ , as functions of the volume fraction of the outer phase $\beta$ . The thin film region (TFR) is highlighted by the grey area. For both the inner ( $j=1$ ) and the outer ( $j=2$ ) phase: (c) dimensionless velocity profiles $u_{j}(y)$ for different values of the viscosity ratio $m$ and a fixed volume fraction $\beta =0.3$ ; (d) average speed $U_{j}$ as a function of $\beta$ and for different values of $m$ .

As described in § 3.2, the decoupled model (3.10) applies to cases where the thermal conductivity of the outer phase is much bigger compared with that of the inner phase, $\mathcal{K}\leqslant \mathcal{O}(\sqrt {\varepsilon })$ . This implies that the inner and the outer phases are thermally decoupled and the model reduces to a ‘one-side’ approach. Such a scenario is typical of the majority of applications, including gas–liquid and liquid–liquid systems (see table 1).

By inspection of (3.10), we see that the temperature evolves in time and space due to three different mechanisms: advection equal to the average speed of the outer phase $U_2$ , effective thermal diffusion through the coefficient $D_2^\star$ , and a constant source term $S_2^\star$ due to the imposed heat flux at the channel wall and viscous dissipation. Specifically, $D_2^\star$ has the typical structure of the Aris–Taylor dispersion coefficient (it scales with the square of the Péclet number), being the sum of heat diffusivity, equal to the unity, and axial shear-induced diffusion, $D_{2}(V,\,\beta, \,m)$ whose expression is given in (C40b ). The latter is an apparent diffusion mechanism for the averaged temperature and has its physical origin in the non-uniformity of the velocity in the outer phase: from the physical point of view, local advection enhances/reduces the axial diffusion in terms of the leading-order variables.

Figure 3(a) shows the evolution of $D_2$ in a coordinate system moving at the mean speed of the outer phase, $V=U_{2}$ , according to Aris (Aris Reference Aris1959). In fact, in the decoupled model, thermal dispersion is exclusively driven by the annulus. The presence of the other phase, which acts as a thermal insulator, reduces the effective heat-diffusion compared with the single-phase scenario where $D_2=2/105$ ; such a limit is reached only when the core vanishes, i.e. for $\beta \to 1$ . The maximal shear-induced diffusion is obtained in the rigid core limit ( $m\to \infty$ ), where the velocity profile of the annulus is almost linear (see figure 3 a), maximising the spread of temperature in the axial direction. In the free surface limit ( $m\to 0$ ), instead, the flow in the outer phase is so slow that the shear-induced diffusion becomes negligible.

Due to the many applications that involve confined thin-films (e.g. Craster & Matar (Reference Craster and Matar2009); Lavalle et al. (Reference Lavalle, Mergui, Grenier and Dietze2021)), we highlight the thin-film region (TFR) in figure 3(a), defined arbitrarily as $\beta \leqslant \sqrt {\varepsilon }\ll 1$ . In the TFR, the outer phase is ‘thin’ enough that $D_{2}$ can be expanded as

(4.1) \begin{equation} D_{2}(V=U_{2},\,\beta \to 0,\,m)= \begin{cases} \dfrac {3\,m^2}{40}\,\beta ^{4}+\mathcal{O}(\beta ^{5}) & \text{if}\,\, m\,\,\text{finite}, \\[12pt] \dfrac {1}{120} \beta ^{2}+\mathcal{O} (\beta ^{3}) & \text{if}\,\, m\to \infty, \end{cases} \end{equation}

showing that the shear-induced diffusion becomes negligibly small in thin films. This can be explained by looking at figure 3(d): when $\beta \to 0$ , the speed of the liquid becomes so small that the system is dominated by diffusion.

In general, advection and diffusion compete in the transient dynamics described by (3.10). To quantify this, similarly to Dejam et al. (Reference Dejam, Hassanzadeh and Chen2014)and Ling et al. (Reference Ling, Tartakovsky and Battiato2016), we introduce the ratio between the effective diffusivity (3.11a ) and its single phase limit (which corresponds to the analogous of the classical Aris–Taylor dispersion coefficient for the thermal problem), $\mathcal{D}_{{AT}}=1+\frac {2}{105}\,{\textit{Pe}}^{2}$ (Horne & Rodriguez Reference Horne and Rodriguez1983; Berkowitz & Zhou Reference Berkowitz and Zhou1996; Wang et al. Reference Wang, Cardenas, Deng and Bennett2012):

(4.2) \begin{equation} \Pi _{2}=\dfrac {\varepsilon ^{-1}\,{\textit{Pe}}\,D_{2}^{\star }}{\mathcal{D}_{{AT}}}= \dfrac {1+{\textit{Pe}}^{2}\,D_{2}}{1+\dfrac {2}{105}\,{\textit{Pe}}^{2}}. \end{equation}

This normalisation facilitates the identification of different regimes, based on the competition between diffusion and advection, as shown in figure 4, where the evolution of $\Pi _{2}$ with respect to the Péclet number and the viscosity ratio at two fixed volume fractions, $\beta =\{0.2,\,0.8\}$ , is presented. Interestingly, for small Péclet numbers ( ${\textit{Pe}}\lt 1$ ), $\Pi _{2}\to 1$ and it is independent of $m$ and $\beta$ , meaning that advective mixing of thermal energy in the outer phase is negligible compared with heat diffusion. When the Péclet number is sufficiently large, the advection dominates and $\Pi _{2}$ assumes always values lower than unity. This means that the core-annular flow pattern does not enhance heat diffusion compared with single-phase flow, in particular, when the volume fraction is small. Also, the mechanism of shear-induced diffusion becomes negligible in the free-surface limit ( $m\to 0$ ) since the velocity in the outer layer is so small that, regardless of the film thickness $\beta$ , $\Pi _{2}$ approaches zero without reaching a plateau as ${\textit{Pe}}\to \infty$ . The regime map in figure 4 can be used in transient multiphase heat-transfer applications to identify a specific set of operational conditions based on fluid and/or flow properties to enhance or reduce heat diffusion.

Figure 4.

Normalised coefficient of shear-induced thermal diffusivity $\Pi _{2}$ for a decoupled system ( $\mathcal{K}\to 0$ ) – see definition in (4.2) – against the flow Péclet number ${\textit{Pe}}$ , for different values of volume fraction $\beta$ and viscosity ratio $m$ .

The source term $S_{2}^{\star }$  (3.11b ) is affected by the dimensionless flux $q_{w}$ and the Brinkman number ${Br}$ . In particular, it scales with the inverse of the thickness of the outer phase $S_2^\star \sim q_w/\beta$ , meaning that heating a thin film is very effective since thermal energy is delivered to a small volume of fluid. The viscous dissipation $W_2^\star =W_2^\star (\beta, \,m )$ defined in (C28b, c ), see (C46c ), contributes to the source term, accounting not only for the viscous heating produced in the annulus, but also for the production into the core that is exchanged across the interface. As shown in figure 3(b), in the free-surface limit ( $m\to 0$ ), the velocity in the annulus is so slow that viscous dissipation is negligible, while for iso-viscous fluids ( $m=1$ ), $W_2^\star /\beta =3$ as for single-phase flow. However, in the rigid-core limit ( $m\to \infty$ ), only the outer phase is dissipative and the heat generated due to friction becomes unbounded as the film gets thinner. Due to this singularity around $\beta =0$ , the scaling in the thin-film region can be obtained using Laurent series for $\beta \to 0$ as

(4.3) \begin{equation} \dfrac {W_{2}^{\star }}{\beta }\left (\beta \to 0,\,m\right )= \begin{cases} 3\,m+\mathcal{O}(\beta ) & \text{if}\,\, m\,\,\text{finite}, \\[6pt] \beta ^{-1}+1+\mathcal{O}(\beta ) & \text{if}\,\, m\to \infty . \end{cases} \end{equation}

The ADHT equation (3.10) for the decoupled regime admits an analytical solution, setting the initial temperature of the outer phase to a uniform value over the entire channel half-length $x\gt 0$ and imposing the same temperature at the inlet $x=0$ for $t\gt 0$ , i.e. $\left \langle \vartheta _{2}\right \rangle (x,\,0 )=\left \langle \vartheta _{2}\right \rangle (0,\,t )\equiv 0$ $\forall \,x,\,t\gt 0,$ we get – see Carslaw & Jaeger (Reference Carslaw and Jaeger1959) and Van Genuchten & Alves (Reference Van Genuchten and Alves1982) –

(4.4) \begin{eqnarray} \left \langle \vartheta _{2}\right \rangle \left (x,\,t\right ) &=& \, S_{2}^{\star }\, \left \{{\vphantom{\left [\frac {x-U_{2}\,t}{2\left (D_{2}^{\star }\,t\right )^{\frac {1}{2}}}\right ]}} t+\frac {(x-U_{2}\,t)}{2\,U_{2}} \operatorname {erfc}\!\left [\frac {x-U_{2}\,t}{2\left (D_{2}^{\star }\,t\right )^{\frac {1}{2}}}\right ] + \right. \nonumber\\[6pt] && \left. -\frac {(x+U_{2}\,t)}{2\,U_{2}} \exp \!\left (\frac {U_{2}\,x}{D_{2}^{\star }}\right ) \operatorname {erfc}\!\left [\frac {x+U_{2}\,t}{2\left (D_{2}^{\star }\,t\right )^{\frac {1}{2}}}\right ]\right \} \,\!. \end{eqnarray}

At long times, the diffusion smears out the temperature profile and the steady-state temperature $\left \langle \vartheta _{2}^{\infty }\right \rangle$ becomes a linear function of $x$ (like in classical internal forced-convection problems, see Incropera (Reference Incropera2007)) with slope equal to $S_{2}^{\star }/ U_{2}$ , i.e.

(4.5) \begin{equation} \left \langle \vartheta _{2}^{\infty }\right \rangle \left (x\right ) = \lim _{t\to +\infty }\left \langle \vartheta _{2}\right \rangle \left (x,\,t\right ) = \frac {S_{2}^{\star }}{U_{2}}\,x. \end{equation}

This solution will be used in § 4.2 to derived a closed-form model for the heat-transfer coefficient. Note that, given the solution for the mean temperature (i.e. its leading order), higher order corrections can be easily calculated plugging (4.4) into (C15), (C23).

4.2. Nusselt number in the decoupled regime

In the decoupled regime, since $\mathcal{K}\to 0$ and $\mathcal{K}\,\mathcal{A}\lesssim \varepsilon$ , the bulk temperature (3.15) simplifies to

(4.6) \begin{align} \vartheta _{b}\left (x,\,t\right )= \left \langle \vartheta _{2}\right \rangle +\sqrt {\varepsilon }\;U_{2}^{-1}\,{\Big \langle u_{2}(y)\,\vartheta _{2}^{(1)}(y)\Big \rangle }+\mathcal{O}(\varepsilon ). \end{align}

Upon substitution of (4.6) in (3.14), the local Nusselt number (3.14) yields

(4.7) \begin{equation} {\textit{Nu}}\left (x,\,t\right )=\frac {4\,q_{w}/ \sqrt {\varepsilon }}{ \left .\vartheta _{2}^{(1)}\right |_{y=1}-U_{2}^{-1}\,{\left \langle u_{2}(y)\,\vartheta _{2}^{(1)}(y)\right \rangle } }. \end{equation}

Equation (4.7) shows how the Nusselt number evolves in space and time and is calculated by plugging the analytical solution (4.4) into the first-order correction (C15) and using the velocity profile for $u_2(y)$ , given in (A5a ). Here, the Nusselt number is studied with respect to the fixed reference frame, in accordance with the literature, where the bulk temperature is always defined in systems with $V=0$ .

The Nusselt number in laminar decoupled core-annular flows with uniform wall heat flux is a function of space, time, the volume fraction, the viscosity ratio, the dimensionless heat-flux, the Brinkman and the Péclet numbers only, i.e. ${\textit{Nu}}={\textit{Nu}} (x,\,t,\,\beta, \,m,\,q_{w},\,{Br}^{\prime },\,{\textit{Pe}} )$ . The Nusselt number and, as a consequence, the heat-transfer coefficient $h$ depend on the Péclet number only in the transient regime, while at the steady-state, ${\textit{Nu}}^{\infty }={\textit{Nu}}^{\infty } (\beta, \,m,\,{Br}^{\prime } )$ . Note that, as typical in the heat transfer community (Shah & London Reference Shah and London1978), the Péclet number can be seen as the product of the Prandtl and Reynolds numbers, i.e. ${\textit{Pe}}={Pr}\,{Re}$ , where ${Pr}=\mu _{2}/\! (\rho _{2}\,\alpha _{2} )$ and ${Re}=\rho _{2}\,\hat {U}\,\hat {H}/\mu _{2}$ , respectively.

When the thermal response reaches the steady state, we get the following expression for the asymptotic Nusselt number ${\textit{Nu}}^{\infty }$ , obtained by combining (4.5), (C15) and (4.7),

(4.8) \begin{align} {\textit{Nu}}^{\infty }(\beta, \,m,\,{Br}^{\prime })= \frac {280(3-\beta )^{2}}{\beta \left (\chi \,{Br}^{\prime }+45\,\beta ^2-245\,\beta +336\right )}, \end{align}

where $\chi$ accounts for the contribution of viscous dissipation as

(4.9) \begin{equation} \chi =\frac {18\,m \left [3\,m\,\beta \left (7 - 3\,\beta \right )^{2} - \left (5\,\beta ^{2} - 35\,\beta + 56\right )\left (1 - \beta \right )^{3}\right ]}{{[1 + \beta (\beta ^{2} - 3\,\beta + 3) (m - 1)]}^{2}}. \end{equation}

Equation (4.8) shows that, in the absence of viscous dissipation ( ${Br}^{\prime }=0$ ) or – equivalently – in the inviscid limit ( $m\to 0$ ), the asymptotic Nusselt number is uniquely a function of the volume fraction of the outer phase $\beta$ . Such an expression embeds the effects of viscous flow in both phases and heat-transfer (advection and diffusion) in the outer phase and it has been obtained starting from first principles. In figure 5(a), the Nusselt number in the non-dissipative regime is represented by the solid black line, which converges to the single phase limit of $140/17$ as $\beta \to 1$ . Instead, in the dissipative regime, the single-phase limit of ${\textit{Nu}}^{\infty }$ approaches $140/\! (108\,{Br}^{\prime }+17 )$ according to (278) of Shah & London (Reference Shah and London1978).

Figure 5.

(a) Dependence of the asymptotic Nusselt number ${\textit{Nu}}^{\infty }$ (absolute value) on the volume fraction $\beta$ , for fixed viscosity ratios $m$ and ${Br}^{\prime }=1$ . Evolution of the normalised Nusselt number ${\textit{Nu}}/{\textit{Nu}}^{\infty }$ : over space at fixed times and Péclet numbers – (b) ${\textit{Pe}}=1$ , (c) ${\textit{Pe}}=0.1$ , (d) ${\textit{Pe}}=0.01$ – over time and across the transport regimes at fixed axial locations – (e) $x=20$ , (f) $x=40$ , with $\varepsilon =0.01$ , ${Br}^{\prime }=0$ , $q_{w}=\beta =0.1$ , $m=1$ .

The viscosity ratio impacts $ {\textit{Nu}}^{\infty }$ only in the dissipative regime, as shown in figure 5(a) for ${Br}^{\prime }=1$ . Interestingly, when the viscosity ratio $m$ is sufficiently large (but finite), there exists conditions where viscous heating ( ${Br}^{\prime }\gt 0$ ) prevails over the wall heat flux ( $q_{w}\gt 0$ ), and the bulk temperature becomes greater than the wall temperature, $\left .\vartheta _{2}\right |_{y=1}-\vartheta _{b}\lt 0$ , reversing the direction of the heat exchange. When this happens, the heat-transfer coefficient $h$ becomes negative, see (3.14), and the singularity of the Nusselt number represents the limiting situation where $\left .\vartheta _{2}\right |_{y=1}=\vartheta _{b}$ . The critical value of volume fraction $\overline {\beta } = \beta (m,\,{Br}^{\prime } )$ which determines the position of the vertical asymptote can be calculated by finding the zeros of the denominator in (4.8) numerically. To ensure that the Nusselt number remains positive in the entire film thickness range ( $0\lt \beta \lt 1$ ), the viscosity ratio and the Brinkman number have to satisfy the inequality $\chi \,{Br}^{\prime }\gt -91/36$ .

Finally, it is interesting to observe that the asymptotic Nusselt number in the TFR scales as follows:

(4.10) \begin{equation} {\textit{Nu}}^{\infty }(\beta \to 0,\,m,\,{Br}^{\prime }\neq 0)=\! \begin{cases} \dfrac {15}{2 \left (1 - 3\,m\,{Br}^{\prime }\right )}\,\beta ^{-1} + \mathcal{O}(1)\,+ \ldots \!& \text{if}\,\, m\, \text{finite}, \\[12pt] \dfrac {60}{7\,{Br}^{\prime }} + \mathcal{O}(\beta ) & \text{if}\,\, m\to \infty, \end{cases} \end{equation}

while for non-dissipative flows ( ${Br}^{\prime }=0$ ), we get that ${\textit{Nu}}^{\infty }\sim 15/\! (2\,\beta )$ . In the finite-viscosity regimes, the asymptotic Nusselt number ${\textit{Nu}}^{\infty }$ shows a singularity. Specifically, increasing $m$ shifts the position of the vertical asymptote towards lower values of $\beta$ , and a larger number of terms in the expansion (4.10) is required to ensure the convergence of the series to (4.8).

The time and spatial evolution of the local Nusselt number, normalised with respect to its asymptotic value for $t\to \infty$ , i.e. ${\textit{Nu}}/{\textit{Nu}}^{\infty }$ , is shown in figure 5(b–d) for different Péclet numbers. The resulting plots have the form of S-shaped breakthrough curves and they are bounded between two horizontal asymptotes: the upper one denotes the asymptotic Nusselt number, while the lower one corresponds to the Nusselt number predicted by (4.7) at the initial time instant, $t=0$ . At sufficiently high Péclet numbers, the S-shaped curves are almost rigidly transported by the flow over time, indicating an advective-dominated behaviour; a reduction of the Péclet number leads to a diffusion-dominated regime (see figure 4).

To determine how the advective and diffusive mechanisms compete in the transient dynamics, we plot the normalised Nusselt number for the different transport regimes considered in figure 5(b–d) over an extended time interval at two fixed positions along the channel: $x=20$ (half-length of the channel) and $x=40$ (end of the channel). At high Péclet numbers ( $ {\textit{Pe}}=1$ ), the heat is primarily advected downstream and the system quickly reaches the fully developed condition. This means that the heat transfer is quite efficient due to the strong influence of advection. Conversely, decreasing the Péclet number ( $ {\textit{Pe}}=\{0.1,\,0.01\}$ ) results in a slower increase of the Nusselt number due to the effect of diffusion. An interesting feature in both figures 5(e) and 5(f) is the existence of a common intersection point between all the curves at $t=t_{e}$ that is independent of the value of the Péclet number. For $t\gt t_{e}$ , the Nusselt number increases with the Péclet number, whereas for $t\lt t_{e}$ , the behaviour of the system is reversed. At a fixed axial position, the curves appear to pivot around this point.

4.3. Hierarchy of time scales in thermal dispersion

The trends illustrated in figure 5 also reflect the development of the forced-convective thermal boundary layer, similarly to the single-phase scenario (Siegel & Sparrow Reference Siegel and Sparrow1959; Siegel & Perlmutter Reference Siegel and Perlmutter1963; Fakoor-Pakdaman et al. Reference Fakoor-Pakdaman, Andisheh-Tadbir and Bahrami2014). As time elapses, the thermal boundary layer progressively penetrates and fills up the thickness $\beta$ occupied by the annulus. As a result, increasingly larger distances after $x=0$ become thermally fully developed. In fact, since a finite amount of time, $t^{\ast }=x^{\ast }/ U_{2}$ , is required for the entrance fluid to be advected downstream and reach the axial position $x^{\ast }$ , beyond this coordinate, $x\geqslant x^{\ast }$ , there has not been any penetration of the fluid which was originally outside the channel before the start of the transient (at $t=0$ ). For $x\lt x^{\ast }$ , the thermal problem reaches its long-time thermal behaviour (linear temperature profile and constant Nusselt number).

Similarly to mass transfer (Allaire, Mikelić & Piatnitski Reference Allaire, Mikelić and Piatnitski2010; Feder, Flekkøy & Hansen Reference Feder, Flekkøy and Hansen2022), the thermal mixing is enhanced in regimes where diffusion has had time to even out the transverse variations of temperature (across the layer) while there are still longitudinal variations at the large scale (along it). Referring to the time scales introduced in § 2.2 and choosing the advective characteristic time scale $\hat {\tau }_{a}$ as the reference time, we obtain the following time scale hierarchy valid for the decoupled model (with $j=2$ and $\mathcal{K}\to 0$ ):

(4.11) \begin{equation} \tau _{\hat {H},\,2} \ll \tau _{a,\,2} \ll \tau _{\hat {L},\,2}\qquad \Rightarrow \qquad \varepsilon \,{\textit{Pe}}\,\beta ^2 \ll U_{2} \ll \dfrac {{\textit{Pe}}}{\varepsilon }. \end{equation}

By inspection of (4.11), we can conclude that if the channel is sufficiently shallow ( $\varepsilon \to 0$ ), these three times scales are sharply separated over a wide range of Péclet numbers, as shown in figure 6. In particular, since the mean velocity of the annulus is bounded between $0$ and $1$ , i.e. $0\leqslant {U_{2}}\leqslant 1$ (see figure 3 d), thermal dispersion is observed for times of the order of $\hat {\tau }_a$ whenever ${\textit{Pe}}\gtrsim \varepsilon$ . For smaller Péclet numbers, the problem approaches its purely diffusive limit, where mixing is not influenced by shear-induced mechanisms. Interestingly, in the inviscid limit ( $m\to 0$ ), the outer phase is almost at rest ( $U_{2}\to 0$ ) and does not experience any shear-induced mixing, while in the TFR limit, the outer phase is so thin that transverse diffusion is almost instantaneous ( $\tau _{\hat {H},2}\to 0$ as $\beta \to 0$ ) at any $ {\textit{Pe}}$ .

Finally, it is worth recalling that our analysis cannot describe the early-time dynamics of the system – which may be referred to as pre-asymptotic dispersion regime (Young & Jones Reference Young and Jones1991; Taghizadeh et al. Reference Taghizadeh, Valdés-Parada and Wood2020) – and the upscaled model (3.10) holds only at times that are much greater than $\tau _{\hat {H},\,2}$ , i.e. $t\gg \varepsilon \,{\textit{Pe}}\,\beta ^2$ , holding no memory from the initial conditions (Ananthakrishnan, Gill & Barduhn Reference Ananthakrishnan, Gill and Barduhn1965; Sankarasubramanian, Gill & Benjamin Reference Sankarasubramanian, Gill and Benjamin1973; Fallon & Chauhan Reference Fallon and Chauhan2005).

Figure 6.

Time scale hierarchy and corresponding dominant heat-transfer mechanisms for the decoupled model ( $\mathcal{K}\to 0$ ) as a function of the volume fraction $\beta$ at different viscosity ratio $m$ and Péclet number ${\textit{Pe}}$ . $\tau _{2}=\{\tau _{\hat {H},\,2},\,U_{2},\,\tau _{\hat {L},\,2}\}$ are the characteristic times of transverse diffusion, advection (average speed of the annulus) and longitudinal diffusion, see (4.11), respectively. The dispersion regime is characterised by the dynamical competition between longitudinal advection and diffusion. $\varepsilon =0.01$ .

4.4. Ramifications for the modelling of thin films and Taylor bubbles

In some circumstances, the modelling of flow patterns, such as intermittent slug flow, relies on an idealisation of the liquid film region, borrowing closure relations from core-annular flows (see for example, Balestra, Zhu & Gallaire (Reference Balestra, Zhu and Gallaire2018); Picchi et al. (Reference Picchi, Ullmann and Brauner2018)). Taylor bubble flow, in fact, consists of a sequence of liquid slugs and elongated bubbles, and – if the bubble is sufficiently long – a region of uniform film thickness forms, which can be considered, as a first approximation, a (local) region of fully developed core-annular flow. This approach neglects the evolution of the thin-film typical of Bretherton’s problem (Bretherton Reference Bretherton1961), while preserving model simplicity. In the context of heat transfer, the widely used three-zone model developed by Thome et al. (Reference Thome, Dupont and Jacobi2004) and Dupont, Thome & Jacobi (Reference Dupont, Thome and Jacobi2004) assumes that the film region of an evaporating/condensing Taylor bubble can be treated in this way.

Figure 7.

(a) Transient evolution of the first-order correction of the temperature profile of the annulus $\vartheta _{2}^{(1)}$ in the decoupled regime at $x=40$ , for $\beta =0.1$ , $m=0.01$ , $\varepsilon =0.01$ , ${\textit{Pe}}=1$ , $q_{w}=0.1$ and ${Br}={Br}^{\prime }=0$ . (b) Right ordinate: time evolution of the wall temperature ( $\vartheta _{2}^{(1)}|_{y=1}$ , black solid lines with circles) and the bulk temperature ( $U_{2}^{-1}\,\langle u_{2}\,\vartheta _{2}^{(1)}\rangle$ , black solid line). Left ordinate: Nusselt number (4.7). The lower asymptote (dotted line) identifies the starting Nusselt number, ${Nu}|_{t\to 0}=16(3-\beta )\beta ^{-1}(8-3\,\beta )^{-1}$ ; the upper asymptote (dashed line) denotes the fully developed Nusselt number ${Nu}^{\infty }$  (4.8).

Specifically, the heat transfer coefficient between the thin film surrounding a Taylor bubble and the channel wall is modelled as $h={\kappa _{2}}/{\delta }$ (where $\delta$ is the film thickness) considering only the heat conduction in the transverse direction (across the film) at the steady-state and neglecting the impact of advection (see Dupont et al. (Reference Dupont, Thome and Jacobi2004); Thome et al. (Reference Thome, Dupont and Jacobi2004); Dai et al. (Reference Dai, Guo, Fletcher and Haynes2015); Magnini & Thome (Reference Magnini and Thome2017); Zhang & Nikolayev (Reference Zhang and Nikolayev2023)). In this way, the heat-transfer coefficient is obtained by assuming a linear temperature profile in the film. The results presented in § 4.2 can be used to check the validity of this hypothesis. Combining (4.10) with (3.14) and (3.19), we obtain the following scaling law for the steady-state heat-transfer coefficient of core-annular flows in a planar geometry:

(4.12) \begin{equation} h=\dfrac {15}{8}\dfrac {\kappa _2}{\beta \,\hat {H}}. \end{equation}

At the steady-state, $h$ is a function of the fluid conductivity and the film thickness only, i.e. $h=h(\kappa _{2},\,\beta \hat {H})$ . Our analysis confirms that the heat-transfer is primarily driven by conduction in the film, and the temperature profile remains almost linear along the transverse direction: when the film is thin enough, the dynamic effects due to the flow are negligible, see figure 3(d), $U_2\to 0$ as $\beta \to 0$ . This can be easily shown combining (C15) with the steady-state axial derivative of the averaged temperature (4.4) for ${Br}=0$ ,

(4.13) \begin{equation} \left .\vartheta _{2}^{(1)}\right |_{t\to \infty } = \frac {q_{w}}{\sqrt {\varepsilon }}\, \left [\frac {5\left (3 + y\right )\left (1 - y\right )^3 - \beta ^3\left (5-\beta \right )} {20\,\beta ^2\left (3 - \beta \right )} + y + \frac {\beta }{2} - 1\right ], \end{equation}

as plotted in figure 7(a). The temperature profile is almost linear in proximity of the wall, while it flattens out at the interface where the heat-flux is zero; this behaviour is embedded in the numerical factor $15/8$ in (4.12).

However, in regimes where the film thickness is not small compared with the channel size (Aussillous & Quéré Reference Aussillous and Quéré2000), the flow in the film is not negligible, and the Nusselt number should be estimated using (4.8) instead of (4.12). In addition, accounting for transient effects would introduce the dependence of the Péclet number into the problem – see figure 7(b).

4.5. Heat-transfer mechanisms and temperature field in the coupled regime

Figure 8.

Effective coefficients of advection as functions of the volume fraction $\beta$ , for different viscosity ratios $m$ , when $\mathcal{K}=1$ in the absolute reference frame ( $V=0$ ). (a) Advective coefficients normalised with the mean flow speed, i.e. $a_{jj}^{\star }/{U_{j}}$ : $j=1$ (left) and $j=2$ (right). (b) Coefficients of coupled advection $a_{12}^{\star }/{U_{1}}$ (left) and $a_{21}^{\star }/{U_{2}}$ (right). Those coefficients are independent of the speed of the reference frame.

When the thermal conductivities of the two phases are of the same order, $\mathcal{K}=\mathcal{O}(1)$ , the heat transfer problem is described by two coupled ADHT equations (3.8), (3.9). Differently from the decoupled regime, the phases can exchange energy through the fluid–fluid interface. By inspection of (3.8), (3.9), we see that the temperatures evolve in time and space due to three different mechanisms: advection, diffusion and storage/source terms. Specifically, each ADHT equation shows (i) a canonical and a coupled effective advection term; (ii) a canonical and a cross-coupling effective diffusion term; (iii) storage and source terms. In the following, we elucidate the physical interpretation of those effective coefficients to describe the main heat-transfer mechanisms of the coupled regime.

The evolution of the advection coefficient $a_{jj}^{\star }=U_{j}+\omega _{jj}^{\star }$ is shown in figure 8(a) in the fixed reference frame ( $V=0$ ). For each phase, the advection can be written as the sum of the mean phase velocity and an extra contribution due to phase coupling. When $\beta$ is small, the lubrication effect of the thin annulus enhances the effective advection of the core that scales as $\beta ^{-1}$ and it is not affected by the viscosity ratio. The coupling advective coefficients $a_{ij}^{\star }$ appear in each ADHT due to the continuity of temperature at the fluid–fluid interface (see § C.5.2) and link one phase with the advection of the other one. Those coupling coefficients are plotted in figure 8(b): when $\beta$ is small (in the TFR), the coupled advection in the core becomes negligible, $a_{12}^{\star }\to 0$ , and, therefore, axial temperature gradients of the annulus do not directly affect the core. Instead, regardless of the value of $\beta$ , $a_{21}^{\star }\to 0$ in the rigid-core limit (for $m\to \infty$ ), where the velocity profile in the core is almost flat, see figure 3(c), and does not trigger any coupled advection mechanism.

Figure 9.

(a) Conductance ratios $G_{j}/ G_{{eq}}$ as functions of the volume fraction $\beta$ , for $\mathcal{K}=1$ . (b) Effective coefficients of diffusion as functions of $\beta$ , for different values of the viscosity ratio $m$ , when $\mathcal{K}=1$ , in the reference frame moving at the mean flow speed ( $V=1$ ).

To investigate the diffusive mechanisms, we can think of heat transfer in the flowing phases via an electrical circuit analogy, as shown in figure 9(a). Specifically, the core-annular flow is represented by a parallel connection of resistors with a specific thermal conductance equal to $G_{1}$ and $G_2$ ; the equivalent conductance is then $G_{\text{eq}}=G_{1}+G_{2}$ . In this framework, we can recast the diffusive coefficients in analogy with the Aris–Taylor formalism to highlight the contribution of heat- and shear-induced diffusion as

(4.14a–b) \begin{equation} d_{11}^{\star }=\underbrace {\frac {G_{1}}{G_{\text{eq}}}}_{\textrm {(i)}}+{\textit{Pe}}^{2}\underbrace {\left (D_{1}+\delta _{jj}^{\star }\right )}_{\textrm {(ii)}}, \quad \quad d_{22}^{\star }=\underbrace {\frac {G_{2}}{G_{\text{eq}}}}_{\textrm {(i)}}+\mathcal{A}^{2}\,{\textit{Pe}}^{2}\underbrace {\left (D_{2}+\delta _{22}^{\star }\right )}_{\textrm {(ii)}}, \end{equation}

The first contribution, (i), is a volume-weighted average of the thermal conductivities of the individual phases and recalls the typical effective description proposed by Maxwell for the heat diffusion coefficient of composites (Maxwell Reference Maxwell1873); its evolution with the volume fraction is shown figure 9(a), where we set $\mathcal{K}=1$ . Note that the ratio $G_{j}/ G_{\text{eq}}$ also appears in front of each time term in (3.8), (3.9), representing the thermal inertia of the phases, i.e. their capacity to store heat and delay its transmission. In the limit of $\mathcal{K}\to 0$ , the ratio $G_{2}/ G_{\text{eq}}$ equals unity and the coefficient reduces to one discussed for the decoupled model. The second contribution in the effective diffusive coefficients, (ii), accounts for the shear-induced diffusion, $D_1$ and $D_2$ , and the coupling between the phases via the terms $\delta _{jj}^{\star }$ . Both the effective diffusion coefficients scale with the square of the Péclet number, while the extra contribution due to the interaction is specific of this class of coupled-layers problems since it enter the ADHT equations by imposing the continuity of temperatures at the fluid–fluid interface (see the derivation in § C.5.2). Those terms are quadratic expressions in the speed of the moving reference frame, $V$ , see (C41), consistent with previous works in the context of mass diffusion, e.g. see Aris (Reference Aris1959).

The evolution of $D_{2}$ with $\beta$ and $m$ has been discussed in § 4.1, figure 3(a), while the shear-induced diffusion for the core, $D_{1}$ , is plotted in figure 9(b) setting $V=1$ . Specifically, in the rigid-core limit, $m\to \infty$ , the shear-induced diffusion becomes negligible since the inner phase has a flat velocity profile, see figure 3(c), preventing any extra spreading across the axial direction rather than heat diffusion. Conversely, in the free-surface limit, $m\to 0$ , the outer phase is almost at rest and behaves like a static coating: $D_{1}$ increases linearly with the thickness $\beta$ as

(4.15) \begin{equation} D_{1}(V=1,\,m\to 0)=\dfrac {2}{105} + \dfrac {\beta }{15}, \end{equation}

showing that a thinner core will result in higher transverse variations in its velocity profile, promoting the shear-induced axial diffusion mechanism. As expected, the single-phase limit of $2/105$ is recovered when $\beta \to 0$ .

Figure 10.

Normalised coefficient of shear-induced thermal conductivity $\Pi _{2}$ for the outer phase of core-annular flows ( $\mathcal{K}=1$ ) against the flow Péclet number ${\textit{Pe}}$ , for different values of volume fraction $\beta$ and viscosity ratio $m$ .

To explore the influence of the flow parameters ( $m$ , $\beta$ , ${\textit{Pe}}$ , $\mathcal{K}$ , $\mathcal{A}$ ) on the heat-transfer mechanisms in the coupled regime, we set $V=1$ and study the ratio between the coefficient of shear-induced thermal diffusivity in (4.14) and its single-phase limit, $\Pi _{j}={d_{jj}^{\star }}/{\mathcal{D}_{\text{AT},\,j}}$ , similarly to what done in § 4.1:

(4.16a–b) \begin{equation} \Pi _{1} = \dfrac {1+\dfrac {1-\beta }{K\,\beta }+\mathcal{A}^{2}{\textit{Pe}}^{2} \left (D_{1}+\delta _{11}^{\star }\right )} {1+\dfrac {2}{105}\,\mathcal{A}^{2}\,{\textit{Pe}}^{2}}, \quad \quad \Pi _{2} = \dfrac {1+\dfrac {K\,\beta }{1-\beta }+{\textit{Pe}}^{2} \left (D_{2}+\delta _{22}^{\star }\right )} {1+\dfrac {2}{105}\,{\textit{Pe}}^{2}}. \end{equation}

The evolution of $\Pi _{2}$ against the Péclet number is shown in figure 10, varying the viscosity ratio $m$ for three different volume fractions $\beta =\{0.05,\,0.2,\,0.8\}$ ; the conductivity ratio $\mathcal{K}$ has been set to unity. The normalised diffusion coefficient of the outer phase has the form of an S-shaped curve and allows the identification of a diffusion-dominated and an advection-dominated regime. Differently from the decoupled model, the purely diffusive regime is characterised by a limiting value that exceeds unity, i.e. $\Pi _{2}\to 1+K\,\beta /{ (1-\beta )}$ for ${\textit{Pe}}\to 0$ , and equal to the conductance ratio $G_{2}/{G_{\text{eq}}}$ , see the expression of $t_{2}^{\star }$ given in (C38). The advective limit is strongly affected by the viscosity ratio and the volume fractions. Specifically, when the thickness $\beta$ is small, the combined effect of advection and the coupling between phases reduces the diffusion coefficients compared with the case of single-phase flow. This effect is more pronounced in the free-surface than in the rigid-core limit. When the annulus is thick enough, instead, $\Pi _{2}$ can overcome the purely diffusive limit for sufficiently small values of the viscosity ratio $m$ . The shear-induced diffusion ratio for the inner phase $\Pi _{1}$ shows a similar behaviour, as presented in figure 11: at low Péclet numbers, $\Pi _{1}\to G_{1}/{G_{\text{eq}}}=1+ (1-\beta )/ (K\,\beta )$ , see the expression of $t_{1}^{\star }$ given in (C38). However, in the advective regime, the shear-induced diffusion in the core is minimal in the rigid-core limit, where transverse variations of velocity are negligible.

Figure 11.

Normalised coefficient of shear-induced thermal diffusivity $\Pi _{1}$ for the inner phase of core-annular flows ( $\mathcal{K}=1$ ) against the flow Péclet number $\mathcal{A}\,{\textit{Pe}}$ , for different values of viscosity ratio $m$ , setting the volume fraction to (a) $\beta =0.05$ and (b) $\beta =0.4$ .

A similar analysis can be done for the cross-coupling diffusive terms ${\textit{Pe}}^{2}\,d_{12}^{\star }$ and $\mathcal{A}^{2}\,{\textit{Pe}}^{2}\,d_{21}^{\star }$ in (3.8), (3.9), introducing the following ratios:

(4.17a-b) \begin{equation} \Pi _{12}=\dfrac {{\textit{Pe}}^{2}\,d_{12}^{\star }}{1+\dfrac {2}{105}\,{\textit{Pe}}^{2}},\qquad \quad \Pi _{21}=\dfrac {\mathcal{A}^{2}\,{\textit{Pe}}^{2}\,d_{21}^{\star }}{1+\dfrac {2}{105}\,\mathcal{A}^{2}\,{\textit{Pe}}^{2}}. \end{equation}

The physical origin of these coupling terms is due to advection only, as can be seen in figure 12. Both coefficients, in fact, tend to zero in the purely diffusive regime ( ${\textit{Pe}}\to 0$ ) and play a role only when advection is important, i.e. at high Péclet numbers. Note that the coupling terms enter the model when the temperature continuity is enforced, see Appendix C.5.2, suggesting that, at high Péclet numbers, heat diffusion in one phase is influenced by diffusion in the other one and vice versa. This mechanism can either play as negative diffusion (with respect to a reference frame that moves with the mean flow $V=1$ , and, therefore, the core is always faster than the annulus), as in the case of the annulus, see figure 12(a), or either enhance or reduce diffusion depending on $m$ and $\beta$ in the case of the core, see figure 12(b).

Figure 12.

Normalised coefficients of cross-coupling diffusivity, (a) $\Pi _{12}$ and (b) $\Pi _{21}$ , for a core-annular system ( $\mathcal{K}=1$ ) against the corresponding flow Péclet number, for different values of viscosity ratio $m$ and volume fraction $\beta$ .

Finally, we study the source terms in (3.8), (3.9) to interpret how the external source of energy $q_w$ affects the model. To do so, observing that the source terms are identical up to a factor equal to $K^2$ , i.e. $K^{2}\,e_{1}^{\star }=e_{2}^{\star }$ , see (C45), and are opposite-signed, we consider the linear combination, $K^{2}$  (3.8) + (3.9), so that the exchange terms cancel out to obtain

(4.18) \begin{equation} \frac {q_{w}}{\varepsilon }\left (K^{2} g_{1}^{\star }+g_{2}^{\star }\right )= \frac {q_{w}\,\hat {H}}{\varepsilon \,\kappa _{2}}\left (G_{1}+G_{2}\right )= \frac {q_{w}\,\hat {H}}{\varepsilon \,\kappa _{2}}\,G_{{eq}}. \end{equation}

The resulting source term is the imposed heat flux at the channel wall multiplied by the sum of the two thermal conductances, see figure 9(a), confirming the consistency of our model. In other words, in the equation for the annulus, the source term is positive and greater with respect to the decoupled model by a factor ${3\,K}/{2 (1-\beta )}\gt 0$ – see (C44) – while in the core, the source term is negative, $-K/{2 (1-\beta )}\lt 0$ , so as to preserve the continuity of the temperature and the thermal flux across the interface. Overall, the energy entering both layers is equal to (4.18) and consistent with the imposed boundary conditions (2.16).

4.6. Transient and steady-state temperature field in the coupled regime

In the coupled regime, the system of coupled ADHT equations (3.8), (3.9) can be solved numerically and admits an analytical solution only at the steady state.

At the steady state, in fact, diffusion smears out any sharp differences of temperature, and the solution $\langle \vartheta _{j}^{\infty }\rangle (x)$ , $j=\{1,\,2\}$ , can be derived by setting to zero the time derivatives and diffusive terms in (3.8), (3.9), yielding an inhomogeneous system of linear first-order ordinary differential equations (ODEs), which can be recast via a linear combination as

(4.19a,b) \begin{align} \left \{ \begin{array}{ll}\dfrac{\textrm {d}{\left \langle \vartheta _{1}^{\infty }\right \rangle }}{\textrm{d}x} = \eta _{1} \left(\left \langle \vartheta _{1}^{\infty }\right \rangle -\left \langle \vartheta _{2}^{\infty }\right \rangle \right ) + \gamma _{1},\\[12pt] \dfrac{\textrm {d}{\left \langle \vartheta _{2}^{\infty }\right \rangle }}{\textrm{d} x} = \eta _{2} \left (\left \langle \vartheta _{1}^{\infty }\right \rangle -\left \langle \vartheta _{2}^{\infty }\right \rangle \right ) + \gamma _{2},\end{array} \right . \end{align}

where

(4.20a) \begin{gather} \lambda =\varepsilon \left (a_{11}^{\star }\,a_{22}^{\star }-a_{12}^{\star }\,a_{21}^{\star }\right ),\quad \eta _{1}=-\frac {a_{22}^{\star }\,e_{1}^{\star }+a_{12}^{\star }\,e_{2}^{\star }}{\mathcal{A}\,{\textit{Pe}}\,\lambda },\quad \eta _{2}=\frac {a_{21}^{\star }\,e_{1}^{\star }+a_{11}^{\star }\,e_{2}^{\star }}{{\textit{Pe}}\,\lambda }, \end{gather}

(4.20b) \begin{gather} \gamma _{1}=\frac {\left (a_{22}^{\star }\,w_{1}^{\star }-a_{12}^{\star }\,w_{2}^{\star }\right ){Br}+ \left (a_{22}^{\star }\,g_{1}^{\star }-a_{12}^{\star }\,g_{2}^{\star }\right )q_{w}}{\mathcal{A}\,{\textit{Pe}}\,\lambda }, \end{gather}

(4.20c) \begin{gather} \gamma _{2}=\frac {\left (-a_{21}^{\star }\,w_{1}^{\star }+a_{11}^{\star }\,w_{2}^{\star }\right ){Br}+ \left (-a_{21}^{\star }\,g_{1}^{\star }+a_{11}^{\star }\,g_{2}^{\star }\right )q_{w}}{{\textit{Pe}}\,\lambda }. \end{gather}

A particular solution to system (4.19) (see Kamke (Reference Kamke1977); Polyanin & Zaitsev (Reference Polyanin and Zaitsev2017)) is represented by two parallel lines of slope $M$ and intercepts $Q_{j}$ :

(4.21) \begin{equation} \big\langle \vartheta _{j}^{\infty }\big\rangle \left (x\right )=Mx+Q_{j},\quad \text{with} \quad M=\frac {\eta _{1}\,\gamma _{2}-\eta _{2}\,\gamma _{1}}{\eta _{1}-\eta _{2}}, \quad Q_{j}=-\frac {\eta _{j}\left (\gamma _{1}-\gamma _{2}\right )}{\left (\eta _{1}-\eta _{2}\right )^{2}}, \end{equation}

meaning that, at the steady state, the temperature of both phases increases linearly in the axial direction with slope equal to $M(q_{w}/{\textit{Pe}},\,\beta, \,m,\,\mathcal{K}\,\mathcal{A},\,{Br}^{\prime })$ given by

(4.22) \begin{equation} M = \frac {2\,q_{w}}{\varepsilon \,{\textit{Pe}}}\left \{\frac {\left (1-\beta \right )^3+m\left [\beta \left (\beta ^2 - 3\,\beta +3\right )+12\,{Br}^{\prime }\right ]}{m\,\beta ^2(3-\beta )+\mathcal{K}\,\mathcal{A}\left [2\left (1-\beta \right )^3+3\,m\,\beta \left (1-\beta \right )(2-\beta)\right ]}\right \}, \end{equation}

where the dimensionless group $\mathcal{K}\,\mathcal{A}$ represents the volumetric heat capacity ratio. Remarkably, in the free-surface limit, $m\to 0$ , $M$ is independent of $\beta$ and equals $q_{w}/(\varepsilon \,{\textit{Pe}}\,\mathcal{K}\,\mathcal{A})$ , whereas the single-phase limits are

(4.23) \begin{equation} \frac {\varepsilon \,{\textit{Pe}}\,M}{q_{w}} =\! \begin{cases} \dfrac {1+12\,m\,{Br}^{\prime }}{\mathcal{K}\,\mathcal{A}} \!& \text{if}\,\, \beta \to 0,\\[10pt] 1+12\,{Br}^{\prime } & \text{if}\,\, \beta \to 1, \end{cases} \end{equation}

recovering, in the absence of viscous dissipation ( ${Br}^{\prime }=0$ ), the expected steady-state thermal behaviour, i.e. $M\propto \{1/(\mathcal{K}\,\mathcal{A});\,1\}$ , see figure 13(a). The difference between the intercepts $Q_{2}-Q_{1}$ corresponds to the steady-state thermal lag between the outer and the inner phase. Interestingly, in the free-surface limit, $m\to 0$ , the normalised thermal lag increases linearly with $\beta$ as $ [4 (1-\beta )+5\,K\,\beta ]/ (10\,K )$ .

Figure 13.

Normalised (a) slope and (b) difference in intercepts for the steady-state solutions (4.21) to the coupled model, as a function of the volume fraction of the outer phase $\beta$ and for fixed values of the viscosity ratio $m$ , corresponding to a liquid–liquid scenario with $\mathcal{A} = 0.51$ and $\mathcal{K} = 5.18$ , see table 1.

Figure 14.

Transient evolution of the averaged temperature $\langle \vartheta _{j}\rangle (x,\,t)$ in the coupled regime, (a) $j=1$ , (b) $j=2$ . (c) Absolute value of the difference between the temperatures of the two phases at the interface $\vartheta _{j}|_{y=1-\beta }$ , including corrections up to the second order. Dashed horizontal lines represent the $\mathcal{O}(\varepsilon ^{1/2},\,\varepsilon, \varepsilon ^{3/2})$ tolerances choosing $\varepsilon =0.01$ . (d) Time evolution of the Nusselt number $\boldsymbol{Nu}(x,\,t)$ against the axial coordinate $x$ . For this liquid–liquid scenario (see table 1), the simulation parameters are: $\boldsymbol{Pe}=1$ , $\mathcal{A}=0.51$ , $\mathcal{K}=5.18$ , $m=0.625$ , $q_{w}=0.1$ , $\beta =0.5$ , $\textit{Br}^{\prime }=0$ . Mesh resolution: $\Delta {x}=0.02$ .

The transient formulation of (3.8), (3.9) can be solved numerically, as shown in figure 14(a, b), choosing a liquid–liquid system from table 1; the numerical solution has been obtained using the pdepe solver of MATLAB, setting an initial temperature equal to zero in both phases and imposing the steady-state heat flux $M$ , given in (4.22), at the right boundary of the domain ( $x=20$ ). A second-order accurate spatial discretisation scheme is employed to convert the original problem into a set of ODEs, which are then integrated to obtain approximate solutions at the specified times (Skeel & Berzins Reference Skeel and Berzins1990). Time discretisation is performed using a multistep variable-step variable-order (VSVO) solver based on numerical differentiation formulae (NDFs) of orders 1–5 (Shampine & Reichelt Reference Shampine and Reichelt1997; Shampine, Reichelt & Kierzenka Reference Shampine, Reichelt and Kierzenka1999). We use equally spaced meshes in the intervals $0\leqslant {x}\leqslant 20$ and $0\leqslant {t}\leqslant 18$ , with 1001 and 501 points, respectively. To ensure highly accurate results, the relative and absolute errors at each integration step have been set to $10^{-12}$ , with norm control enabled to manage the overall error in the solution vector. As time elapses, the solution evolves and an increasing region of the channel reaches the steady state.

Interestingly, we can use the numerical solution to check the model consistency and verify that the continuity of temperature at the fluid–fluid interface has been properly enforced. This issue represents a key aspect of the model, since such a boundary condition has been imposed in asymptotic terms, see § C.5.2. From figure 14(c), we can see that the temperature difference, $\vartheta _{1}-\vartheta _{2}= (\left \langle \vartheta _{1}\right \rangle +\sqrt {\varepsilon }\,\vartheta _{1}^{(1)}+\varepsilon \,\vartheta _{1}^{(2)} )- (\left \langle \vartheta _{2}\right \rangle +\sqrt {\varepsilon }\,\vartheta _{2}^{(1)}+\varepsilon \,\vartheta _{2}^{(2)} )$ evaluated at the interface and keeping the first three terms of the expansions is within the $\mathcal{O}(\varepsilon \sqrt {\varepsilon })$ error, consistent with the two-scale asymptotic expansion in (3.3). Our model is in fact accurate up to the order $\mathcal{O}(\varepsilon )$ .

4.7. Nusselt number in the coupled regime

In the coupled regime, the bulk temperature $\vartheta _{b}$ , see (3.15), simplifies to

(4.24) \begin{equation} \vartheta _{b}\left (x,\,t\right )=\frac {\mathcal{K}\,\mathcal{A}\left (1-\beta \right )U_{1}\left (\left \langle \vartheta _{1}\right \rangle +\sqrt {\varepsilon }\,\vartheta _{1}^{(1)}\right )+\beta \,U_{2}\left (\left \langle \vartheta _{2}\right \rangle +\sqrt {\varepsilon }\,\vartheta _{2}^{(1)}\right )}{\mathcal{K}\,\mathcal{A}\left (1-\beta \right )U_{1}+\beta \,U_{2}} +\mathcal{O}(\varepsilon ), \end{equation}

and the Nusselt number can be computed though the definition (3.19), using the numerical solution of (3.8), (3.9) for the averaged temperatures $\left \langle \vartheta _{j}\right \rangle$ to estimate the first-order terms in (4.24) via (C12) and (C15). The Nusselt number in laminar coupled core-annular flows with uniform wall heat-flux is a function only of ${\textit{Nu}}={\textit{Nu}} (x,\,t,\,\beta, \,m,\,\mathcal{A},\,\mathcal{K},\,q_{w},\,{Br}^{\prime },\,{\textit{Pe}} )$ and, based on the definition of the heat-transfer coefficient, see (3.12), it allows for obtaining the two-phase heat-transfer coefficient.

An example of the time and spatial evolution of ${\textit{Nu}} (x,\,t )$ is given in figure 14(d): the Nusselt number evolves according to an S-shaped trend and is bounded between a lower and an upper horizontal asymptote, corresponding respectively to its limiting values for $t\to \infty$ and $t\to 0$ . The impact of the Péclet number (i.e. the competition between advection and diffusion) is qualitatively similar to what is described in figure 5(b–f) for the decoupled regime and it is not shown here only for the sake of brevity.

Figure 15.

Limiting two-phase Nusselt number for a core-annular flow of unitary conductivity ratio, $\mathcal{K}=1$ , as a function of the volume fraction $\beta$ , for fixed values of the viscosity ratio $m$ . Panels (a) to (c) refer to a non-dissipative core-annular flow ( ${Br}^{\prime }=0$ ) with increasing diffusivity ratios as increasing powers of the small-scale parameter, $\varepsilon =0.01$ : (a) $\mathcal{A}=\varepsilon$ , (b) $\mathcal{A}=\sqrt {\varepsilon }$ , (c) $\mathcal{A}=1$ . In (d), $\mathcal{A}={Br}^{\prime }=1$ .

At the steady state, the axial gradients of temperature $\partial \langle \vartheta _{j}\rangle /\partial x$ can be replaced by the slope $M$ given in (4.21) and the Nusselt number $ {\textit{Nu}}^{\infty }= {\textit{Nu}}^{\infty } (\beta, m, \mathcal{A}, \mathcal{K}, {Br}^{\prime } )$ can be written in a closed form (its full expression looks quite cumbersome and it is not reported here only for the sake of brevity). Figure 15 shows the evolution of $ {\textit{Nu}}^{\infty }$ with respect to the volume fraction and the viscosity ratio for a liquid–liquid case, see table 1, keeping the conductivity ratio as $\mathcal{K}=1$ ; the panels ( $a$ )–( $c$ ) show the effect of the diffusivity ratio $\mathcal{A}$ . First, the single-phase limits of the steady-state Nusselt number for a vanishing core (i.e. the channel section is entirely occupied by the outer phase) yields to

(4.25) \begin{equation} {\textit{Nu}}^{\infty }\left (\beta \to 1,\,m,\,\mathcal{A},\,\mathcal{K},\,{Br}^{\prime }\right )=\! \begin{cases} 4 + \mathcal{O}(1-\beta )& \text{if}\,\, m \to 0,\\[6pt] \dfrac {140}{108\,{Br}^{\prime }+17} + \mathcal{O}(1-\beta ) & \text{otherwise}, \end{cases} \end{equation}

confirming that the result obtained as $m\not \to 0$ is consistent with the single-phase benchmark expression given by (278) of Shah & London (Reference Shah and London1978) ( ${\textit{Nu}}^{\infty }=140/17$ ) for a uniform heat flux boundary condition imposed in a planar channel. Note that in the free-surface limit, $m\to 0$ , the outer phase moves so slowly (see figure 3 d) that the heat flux cannot be convected downstream and the Nusselt number converges to the value of $4$ , in agreement with the case of two different temperatures specified at each boundary, see (268) of Shah & London (Reference Shah and London1978).

Interestingly, when $m\to 0$ , the convective heat-transfer is less efficient compared with the single-phase flow and ${\textit{Nu}}^{\infty }\lt 140/17$ in the whole range of volume fractions, see figure 15(a–c). In the free-surface limit, in fact, the annulus moves so slowly (see figure 3 c) that the asymptotic Nusselt number is unaffected by the thermal diffusivity ratio $\mathcal{A}$ :

(4.26) \begin{equation} {\textit{Nu}}^{\infty }\left (\beta, \,m\to 0,\,\mathcal{K}\right )=\frac {140\,K}{35\,K\,\beta +17\left (1-\beta \right )}. \end{equation}

Notice that, in the coupled regime ( $k=0$ ), the use of $K$ or $\mathcal{K}=K\,\varepsilon ^{k}$ is equivalent in (4.26). In the rigid core limit, instead, convective heat-transfer is favoured and is always enhanced compared with the single-phase flow, namely ${\textit{Nu}}^{\infty }\gt 140/17$ , in the whole range of $\beta$ . This can be explained by the shape of the velocity profile (linear in the annulus and flat in the core, see figure 3 c), that maximises heat transfer ensuring the most efficient replacement of fluid over the heated surface. This effect is amplified by lowering the diffusivity ratio $\mathcal{A}$ . Specifically, the thermal diffusivity represents the promptness of a material in dissipating a temperature inhomogeneity relative to its tendency to store thermal energy, and, therefore, when $\mathcal{A}$ is less than one, the diffusivity of the core is greater than that of the annulus. In our case, a lower value of $\mathcal{A}$ results in a greater tendancy of the outer phase to accumulate the heat received from the surroundings rather than diffusing it, enhancing the heat-transfer coefficient. This behaviour is favoured at low volume fractions, i.e. when the fluid in contact with the channel wall is sufficiently thin. In contrast, increasing $\beta$ leads to a larger thermal resistance in the annulus and gives a greater difference between the wall and the bulk temperature defined in (4.24), reducing the heat-transfer coefficient $h$ and the Nusselt number.

In other words, two-phase flows enhance convective heat transfer only under certain conditions. Specifically, only if the viscosity ratio is finite and the diffusivity ratio is small enough, we observe an enhanced heat-transfer coefficient, in particular at low volume fractions.

Finally, the effect of viscous dissipation on the Nusselt number is shown in figure 15(d), setting ${Br}^{\prime }=1$ while keeping $\mathcal{A}=1$ . When the two fluids have the same viscosity, i.e. $m=1$ , the Nusselt number does not depend on the volume fraction and viscous dissipation lowers the value of $140/(108\,{Br}^{\prime }+17)$ . The curves obtained for values of $m$ larger than $1$ are located below this threshold, while those where $0\lt m\lt 1$ lie above it. In any case, viscous heating ( ${Br}^{\prime }\gt 0$ ) reduces the efficiency of convective heat-transfer and the dependence on the viscosity ratio is flipped compared with the non-dissipative scenario. This can be attributed to the combination of two factors (Shah & London Reference Shah and London1978): (i) a reduction in the wall temperature gradient near the wall region due to viscous heating; and (ii) a slower rise of the bulk temperature along the channel axis, due to a reduced amount of heat transferred through the wall.