Close-contact melting on hydrophobic textured surfaces: confinement and meniscus effects

Abstract

We investigate the dynamics of close-contact melting (CCM) on ‘gas-trapped’ hydrophobic surfaces, with specific focus on the effects of geometrical confinement and the liquid–air meniscus below the liquid film. By employing dual-series and perturbation methods under the assumption of small meniscus deflections, we obtain numerical solutions for the effective slip lengths associated with velocity $\lambda$ and temperature $\lambda _t$ fields, across various values of aspect ratio $\Lambda$ (defined as the ratio of the film thickness $h$ to the structure’s periodic length $l$) and gas–liquid fraction $\phi$. Asymptotic solutions of $\lambda$ and $\lambda _t$ for $\Lambda \ll 1$ and $\Lambda \gg 1$ are derived and summarised for different surface structures, interface shapes and $\Lambda$, which reveal a different trend of $\lambda$ for $\Lambda \ll 1$ and depending on the presence of a meniscus. In the context of constant-pressure CCM, our results indicate that longitudinal grooves can enhance heat transfer under the effects of confinement and a meniscus when $\Lambda \lesssim 0.1$ and $\phi \lt 1 - 0.5^{2/3} \approx 0.37$. For gravity-driven CCM, the parameters of $l$ and $\phi$ determine whether the melting rate is enhanced, reduced or nearly unaffected. We construct a phase diagram based on the parameter matrix $(\log _{10} l, \phi )$ to delineate these three regimes. Lastly, we derive two asymptotic solutions for predicting the variation in time of the unmelted solid height.

1. Introduction

In the domain of melting, a phenomenon known as close-contact melting (CCM) stands out. Unlike conventional convection-driven melting with expanding molten liquid spaces, CCM involves a distinct process where the surface of an unmelted solid and a heating element are pressed together under external force. This results in the formation of a thin molten film flow between the solid and heating surfaces, facilitating solid–liquid phase change heat transfer at a low-Reynolds-number flow of the liquid. Close-contact melting can be observed in two primary modes: heat source-driven (constant pressure) and unmelted solid-driven (gravity-driven) (Hu et al. 2022), depending on the external force applied to the object. The former mode refers to the scenario where the heated surface penetrates into an unmelted solid, which includes glacier drilling by heating (Schuller, Kowalski & Raback 2016), subtractive machining (Mayer & Moaveni 2008) and the ‘melt down’ of nuclear fuel (Emerman & Turcotte 1983). The latter mode can be observed readily in daily life, such as when ice, butter or chocolate melts on a heated substrate. Also, it is found in latent heat thermal energy storage (Kozak, Rozenfeld & Ziskind 2014) and thermal management systems (Fu et al. 2022), where a considerable heat flux can be achieved easily. In the area of thermal energy storage, the melting rate is equivalent to the charging rate.

Given the widespread applications of CCM, it has been studied extensively for decades (Hu et al. 2022). These studies explore the effects of geometry (Moore & Bayazitoglu 1982; Sparrow & Geiger 1986; Dong et al. 1991; Fomin, Wei & Chugunov 1995) and boundary conditions for heating (Moallemi & Viskanta 1985; Saito et al. 1985a , b ; Bejan 1992; Fomin, Saitoh & Chugunov 1997; Groulx & Lacroix 2007; Schueller & Kowalski 2017), with a primary focus on heat transfer characteristics. Recently, attention has shifted towards understanding the mechanisms of flow and heat transfer within the liquid film. Hu et al. measured the thickness variation of the thin film for the first time, confirming a magnitude of $O(10^2)$ $\unicode{x03BC}$ m (Hu et al. 2019). Non-Newtonian features, including Bingham (Kozak et al. 2019), power-law (Kozak 2023) and Carreau properties for the viscosity variations (Hu & Fan 2023), as well as temperature-dependent properties (Cregan, Williams & Myers 2020), have been considered to study their impact on film thickness variation and flow behaviour. Strategies such as the use of permeable surfaces (Turkyilmazoglu 2019), slip surfaces (Aljaghtham, Premnath & Alsulami 2021; Kozak 2022) and pressure-enhanced conditions (Chen et al. 2022; Fu et al. 2022) have been proposed to accelerate melting.

In addition, hydrophobic or superhydrophobic surfaces with microscale-trapped air can provide an effective slip length, of the order of $O(10^2)$ $\unicode{x03BC}$ m, for liquid transport (Lee, Choi & Kim 2008; Quéré 2008). This characteristic is considered a promising strategy for reducing liquid film thickness, thereby potentially enhancing heat transfer of CCM. Recently, Kozak theoretically studied postpatterned hydrophobic surfaces and showed that a temperature slip length $\lambda ^{*}_{t}$ (the definition is given in § 2.2.1) greater than velocity slip length $\lambda ^{*}$ (Kozak 2022) leads to a reduced heat transfer; the calculation utilised the formulae $ \lambda ^{*} / l^{*}=3 \sqrt {\pi /(1-\phi )}/16-3\ln (1+\sqrt {2})/(2 \pi )$ (Davis & Lauga 2010) and $\lambda ^{*}=3\lambda ^{*}_t/4$ (Enright et al. 2014), where $^{*}$ represents dimensional quantities, $l^{*}$ is the period of a post pattern and $\phi$ is the liquid–gas area fraction. Although this work established a framework for analysing velocity and temperature slip effects on CCM, the expressions for $\lambda ^*$ rely on three assumptions: (i) shear flow in the liquid above the surface; (ii) a flat gas–liquid interface; (iii) perfect slip on the gas–liquid interface. Assumption (iii) is satisfied when $ h_s^{*} {\mu _l^{*}}/{\mu _g^{*}} \gg l^{*}$ (Belyaev & Vinogradova 2010; Asmolov & Vinogradova 2012), where $h_s^{*}$ , $ {\mu _l^{*}}$ and $\mu _g^{*}$ denote microstructure height, liquid viscosity and gas viscosity, respectively. However, assumptions (i) and (ii) may not apply to CCM due to the features of Poiseuille flow and deformation of the meniscus, leading to a complicated effective slip. The assumption of a constant ratio $\lambda _t^*/ \lambda ^*$ is valid for large aspect ratios $\Lambda \equiv h^*/l^*$ (Enright et al. 2014), where $h^*$ is the liquid film thickness. Furthermore, the slip lengths are governed by the geometrical patterns of the microstructure, which should be evaluated and compared.

In the literature on slip, since the early works by Philip (1972a , b ), it has been shown that the velocity slip length $\lambda ^{*}$ is influenced by confinement effects (Lauga & Stone 2003; Sbragaglia & Prosperetti 2007; Teo & Khoo 2008; Feuillebois, Bazant & Vinogradova 2009; Schnitzer & Yariv 2017; Landel et al. 2020). The asymptotic limit of slip length for flow past transverse slip regions inside a tube (radius $R^*$ ) has been determined (Lauga & Stone 2003), specifically $\lambda ^{*}_{\perp } / l^{*} \sim \ln (\sec (\phi \pi /2 ) )/(2\pi )$ when $l^{*}/R^{*} \ll 1$ , while $\lambda ^{*}_{ \perp }/R^{*} \sim \phi /(4-4\phi )$ when $l^{*}/R^{*} \gg 1$ ; this result indicates that the velocity slip length can depend on the tube radius $R^{*}$ instead of the pattern period $l^{*}$ . Sbragaglia & Prosperetti (2007) numerically obtained for longitudinal grooves a variation of velocity slip length $\lambda _{\parallel } = \lambda _{\parallel }^{*}/l^{*}$ for various aspect ratios $\Lambda = h^{*} / l^{*}$ , where longitudinal means grooves oriented parallel to the flow, demonstrating that $\lambda _{\parallel }$ remains constant when $\Lambda \geqslant 10$ , but shows a scaling law $\lambda _{\parallel } \sim \Lambda$ when $\Lambda \lt 10$ . Similar confinement effects on the effective slip length can be found for either longitudinal or transverse grooves with symmetrical or asymmetrical boundary conditions (Teo & Khoo 2008; Schnitzer & Yariv 2017) or on interfacial velocity at the liquid–gas interface (Landel et al. 2020). An interesting fact related to confinement effects is that longitudinal grooves provide the largest slip lengths while transverse grooves yield the smallest ones in the limit of $\Lambda \ll 1$ (Feuillebois et al. 2009). This contrasts with predictions for $\Lambda \gg 1$ , where arrays of posts have a larger slip length than grooves. Given the special feature of CCM that film thickness $h^{*}$ (also the channel height) is variable and sensitive to conditions, we can consider the examples from Kozak (2022), where $l^{*} \sim 50 \ \unicode{x03BC}$ m and $h^{*} \sim 50 \ \unicode{x03BC}$ m ( $\phi =0.99$ ), as well as $l^{*} \sim 16 \ \unicode{x03BC}$ m and $h^{*} \sim 80 \ \unicode{x03BC}$ m ( $\phi =0.84$ ). In this case, confinement effects may be non-negligible due to $1 \leqslant \Lambda \leqslant 5$ .

At the same time, another issue arises concerning the curved gas–liquid interface due to the wide variation of pressure exerted on the liquid film in CCM. It has been demonstrated that meniscus protrusion into cavities will either increase or decrease the slip length $\lambda ^{*}$ depending on the aspect ratio $\Lambda$ (Sbragaglia & Prosperetti 2007; Teo & Khoo 2010; Crowdy 2016; Game et al. 2017; Kirk, Hodes & Papageorgiou 2017). When $\Lambda \gg 1$ , both slip lengths $\lambda ^{*}_{\perp }$ and $\lambda ^{*}_{\parallel }$ decrease along with a larger meniscus angle $\theta$ that protrudes into grooves (Teo & Khoo 2010). However, for $ \Lambda \ll 1$ , the slip length $\lambda ^{*}$ can be increased by a meniscus for symmetrical boundaries (Sbragaglia & Prosperetti 2007; Kirk et al. 2017). Besides the velocity slip length $\lambda ^{*}$ , confinement and meniscus effects also play an important role in the temperature slip length $\lambda ^{*}_{t}$ (Enright et al. 2014; Kirk et al. 2017; Hodes et al. 2023; Tomlinson et al. 2024). The primary conclusions drawn are that the confinement effect results in a larger temperature slip length, while meniscus influences tend to reduce it. However, for CCM studies, the latter has not been considered (Kozak 2022).

1.1. Scope and structure of this work

Since constant-pressure and gravity-driven CCM are characterised by constant and growing liquid film thicknesses, respectively, as well as also corresponding to two different application scenarios, melt drilling and thermal energy storage, respectively, we focus on parallel-groove textured surfaces as a promising candidate for both modes of CCM, owing to the theoretical limit of a significant slip length for a given gas–liquid contact ratio (Feuillebois et al. 2009). Furthermore, by incorporating double re-entrant structures, these surfaces can readily maintain the Cassie state for a wide range of liquids (Liu & Kim 2014; Wilke et al. 2022). In § 2, we introduce the theoretical framework for CCM, considering the influences of the velocity slip length $\lambda$ and the temperature slip length $\lambda _t$ (noting that quantities without a $^*$ are dimensionless), where the subscript $\parallel$ represents longitudinal grooves; also, $f$ and $c$ denote flat interfaces and curved interfaces (menisci), respectively. Specifically, dimensional analysis of the problem is presented in § 2.1, followed by the derivation of the effective thermal slip length $\lambda _{t,f}$ and $\lambda _{t,c}$ , corresponding to a flat gas–liquid interface (§ 2.2.1) and a curved surface (§ 2.2.2), respectively. The effective velocity slip lengths $\lambda _{\parallel, f}$ and $\lambda _{\parallel, c}$ are derived for flat (§ 2.3.1) and the curved (§ 2.3.2) liquid–gas interface, respectively. Slip lengths for all cases are solved by the standard dual-series method, whose details are introduced in the appendices. Based on the obtained effective slip lengths, a generalised one-dimensional description is developed for both constant-pressure and gravity-driven CCM in § 2.4.

In § 3.1, we discuss the influences of the dimensionless film height $\Lambda$ on the slip lengths for the flat (§ 3.1.1) and curved (§ 3.1.2) gas–liquid interfaces. Asymptotic solutions of $\lambda$ and $\lambda _t$ for $\Lambda \gg 1$ and $\Lambda \ll 1$ are summarised in table 1. In addition, a modified asymptotic solution for $\lambda _{\parallel, c}$ is proposed for a wide range of $\Lambda \gtrsim 0.2$ . We then discuss the slip effects on constant-pressure CCM in § 3.2, revealing the critical conditions of $\Lambda$ and $\phi$ to achieve a faster melting rate. On the other hand, we analyse gravity-driven CCM in § 3.3, under limiting conditions derived from the dimensionless periodic length $l$ , and construct a $l$ - $\phi$ phase diagram to illustrate the influences of slip.

Table 1. Dimensional asymptotic formula of velocity ( $\lambda _{\parallel }^{*}$ or $\lambda _{\perp }^{*}$ ) and thermal ( $\lambda _{t}^{*}$ ) slip lengths for various confinement and meniscus effects.

2. Physical model and theoretical framework

A typical CCM process on a microtextured hydrophobic surface is sketched in figure 1(a). Here and in the following, the notation $^{*}$ represents a dimensional quantity. A cuboid-shaped solid phase change material (PCM) with length $L^{*}$ , width $W^{*}$ and initial height $H_{0}^{*}$ is heated from below by a textured (microgrooved) surface maintained at wall temperature $T_{w}^{*}$ and characterised by a periodic length $l^{*}$ . The solid PCM, initially at its melting point $T_{m}^{*}\lt T_{w}^{*}$ , will continuously melt from the bottom under a certain pressure, which could be either a specific constant value $\mathcal{P}^{*}$ or the self-weight pressure of the PCM, $\rho _s^{*} g^{*} (H^{*}-h^{*})$ , depending on the mode of CCM. Given the microgrooved structure (periodic length $l^*$ and gas fraction $\phi$ ), the mixed interface, composed of solid–liquid and gas–liquid contacts, when approximated with a one-dimensional description, introduces effective coefficients of temperature slip $\lambda _{t}^{*}$ and velocity slip $\lambda ^{*}$ , which are described below and significantly affect the fluid flow and heat transfer in the CCM process.

Figure 1. (a) Close-contact melting of a cuboid-shaped unmelted solid with initial height $H_{0}^{*}$ , length $L^{*}$ and width $W^{*}$ pressed downwards by constant pressure $\mathcal{P}^{*}$ or self-weight $\rho _s^*g^*(H^*-h^*)$ on a heated microgrooved hydrophobic surface with characteristic periodic length $l^{*}$ , where $h^*$ is the liquid film thickness; quantities with $^{*}$ are dimensional. Schematic diagrams of boundary conditions for (b) two-dimensional descriptions, or (c) a one-dimensional description, for the temperature distribution and temperature slip length $\lambda _{t}$ ; similarly, (d)–(e) characterise the velocity boundary conditions and velocity slip length $\lambda$ on the longitudinal grooves. It is noted that length is scaled by $l^{*}$ for convenience. (f) The dimensionless schematic diagram ( $L^* \ll W^*$ ) by scaling as $X = x^*/L^*$ , $Y = y^*/h_0^*$ and $T = (T^*-T_m^*)/(T_w^*-T_m^*)$ , showing the remaining solid height $H$ and film thickness $h$ is influenced by the effective velocity slip length $\lambda$ and temperature slip length $\lambda _{t}$ .

2.1. Scaling analysis and non-dimensionalization

Typically, a thin melted film of liquid separates the substrate from the solid PCM. Hence, a lubrication-style analysis is suggested. The relevant physical parameters are the gravitational acceleration $g^{*}$ , viscosity $\mu ^{*}$ , solid density $\rho _s^{*}$ , liquid density $\rho _l^{*}$ , liquid specific heat capacity $c_{p,l}^{*}$ , thermal conductivity of liquid $k_l^{*}$ and latent heat of fusion $\mathcal{H}^*$ . By scaling pressure and stresses with $p_c^* = \rho _s^*g^*(H_0^*-h_0^*)$ for the gravity-driven mode or $p_c^* = \mathcal{P}^*$ for the constant-pressure mode, velocity with $u_c^* = {h_0^*}^2p_c^*/(\mu ^*L^*)$ , time with $t_c^* = L^*/u_c^*$ and temperature with $T_w^{*}-T_m^{*}$ , we can define dimensionless variables as

(2.1) \begin{equation} \begin{gathered} t =\frac {t^*}{t_c^*}, \quad X=\frac {x^*}{L^*}, \quad Y=\frac {y^*}{h_0^*}, \quad Z=\frac {z^*}{l^*}, \quad u = \frac {u^*}{u_c^*}, \quad v = \frac {v^*L^*}{u_c^*h_0^*}, \quad w = \frac {w^*L^*}{u_c^*l^*},\\ \quad P = \frac {P^*}{p_c^*}, \quad T=\frac {T^*-T_m^*}{T_w^*-T_m^*},\quad h = \frac {h^*}{h_0^*}, \quad H =\frac {H^*}{H_0^*}, \quad \varDelta = \frac {H^*-h^*}{H_0^*-h_0^*}, \\ \mathscr{h} = \frac {h_0^*}{L^*}, \quad A = \frac {H_0^*}{L^*}, \quad \mathscr{l}=\frac {l^*}{L^*}, \quad l =\frac {l^*}{h_0^*},\quad \Lambda =\frac {h^*}{l^*}=\frac {h}{l}, \end{gathered} \end{equation}

where $h_0^*$ is the initial film thickness identified below as (2.51) in § 2.4. The aspect ratio $\Lambda$ is a key parameter related to the confinement effect in the following discussion. Hence, the dimensionless variables are denoted as the velocity vector $\boldsymbol{V}\equiv u\boldsymbol{e}_X+v\boldsymbol{e}_Y+w\boldsymbol{e}_Z$ , temperature $T$ , pressure $P$ , spatial gradient operator $\boldsymbol{\nabla }\equiv \boldsymbol{e}_X\partial _X+\boldsymbol{e}_Y\partial _Y +\boldsymbol{e}_Z\partial _Z$ and $\partial _t$ is the time derivative. The Reynolds number $\textrm{Re}$ , Péclet number $\textrm{Pe}$ , Prandtl number $\textrm{Pr}$ , Eckert number $\textrm{Ec}$ , Stefan number $\textrm{St}$ , density ratio $\rho$ and hydrostatic pressure ratio $p_h$ are, respectively, defined as

(2.2a) \begin{align}&\qquad\qquad\textrm{Re} \equiv \frac {\rho _l^*L^*u_c^*}{\mu ^*}, \quad \textrm{Pe} \equiv \frac {u_c^*L^*\rho _l^*c_{p,l}^*}{k_l^*}, \quad \textrm{Pr} \equiv \frac {\mu ^{*} c_{p,l}^{*}}{k_l^{*}}, \end{align}

(2.2b) \begin{align}&\textrm{Ec} \equiv \frac {{u_c^*}^2}{ c_{p,l}^{*} \left (T_w^{*}-T_m^{*}\right )}, \quad \textrm{St} \equiv \frac {c_{p,l}^{*} \left (T_w^{*}-T_m^{*}\right )}{\mathcal{H}^*}, \quad \rho \equiv \frac {\rho _s^{*}}{\rho _l^{*}}, \quad p_h = \frac {\rho _l^*g^*h_0^*}{p_c^*}. \end{align}

In this case, the equations of continuity, momentum (in the $X$ , $Y$ -and $Z$ directions), and energy, including the dissipation of mechanical to thermal energy in the bulk liquid and an energy balance at the melting front, are, respectively,

(2.3a) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad \boldsymbol{\nabla } \cdot \boldsymbol{V}=0, \end{align}

(2.3b) \begin{align}&\qquad\qquad \textrm{Re}\mathscr{h}^2\partial _t u + \textrm{Re}\mathscr{h}^2 ({\boldsymbol{V}} \cdot \boldsymbol{\nabla }) {u} = -\partial _X P + \mathscr{h}^2 \partial _X^2 u + \partial _Y^2 u + l^{-2} \partial _Z^2 u, \end{align}

(2.3c) \begin{align}& {}\qquad\textrm{Re}\mathscr{h}^4\partial _t v + \textrm{Re}\mathscr{h}^4 ({\boldsymbol{V}} \cdot \boldsymbol{\nabla }) v = -\partial _Y P + \mathscr{h}^4 \partial _X^2 v + \mathscr{h}^2 \partial _Y^2 v + l^{-2} \mathscr{h}^2 \partial _Z^2 v + p_h, \end{align}

(2.3d) \begin{align}& {}\qquad\,\,\,\textrm{Re}\mathscr{h}^2\mathscr{l}^2 \partial _t w + \textrm{Re}\mathscr{h}^2 \mathscr{l}^2 ({\boldsymbol{V}} \cdot \boldsymbol{\nabla }) w = -\partial _Z P + \mathscr{h}^2 \mathscr{l}^2 \partial _X^2 w + \mathscr{l}^2 \partial _Y^2 w + \mathscr{h}^2 \partial _Z^2 w, \end{align}

(2.3e) \begin{align}& \textrm{Pe}\mathscr{h}^2\partial _t T + \textrm{Pe}\mathscr{h}^2 ({\boldsymbol{V}} \cdot \boldsymbol{\nabla }) {T} = \mathscr{h}^2 \partial _X^2 T + \partial _Y^2 T + l^{-2} \partial _Z^2 T + \textrm{PrEc} \left [ \left ( \mathscr{h} \partial _X u \right )^2 + \left ( \partial _Y u \right )^2 \right . \nonumber\\ &\qquad\qquad\qquad\qquad\qquad \left . + \left ( l^{-1} \partial _Z u \right )^2 + \left ( \mathscr{h}^2 \partial _X v \right )^2 + \left ( \mathscr{h} \partial _Y v \right )^2 + \left ( l^{-1} \mathscr{h} \partial _Z v \right )^2 \right . \nonumber\\ &\qquad\qquad\qquad\qquad\qquad \left . + \left ( \mathscr{h} \mathscr{l} \partial _X w \right )^2 + \left ( \mathscr{l} \partial _Y w \right )^2 + \left ( \mathscr{h} \partial _Z w \right )^2 \right ], \end{align}

(2.3f) \begin{align}&\qquad\qquad\qquad\qquad \frac {\textrm{St}}{\mathscr{h}(A-\mathscr{h}) \textrm{Pe} \rho } \partial _Y T(X,h) = \partial _t \varDelta, \end{align}

where $\varDelta$ is the remaining height of the solid to be melted and $\partial _t\varDelta$ represents the melting rate.

In accordance with previous experimental and theoretical results (Moallemi, Webb & Viskanta 1986; Hu & Fan 2023), the following well-validated assumptions are adopted: (i) the lubrication approximation is valid due to $ \mathscr{h}^2 \ll 1$ ; (ii) the flow is quasisteady, $\partial /\partial {t}=0$ , as a consequence of a low effective Reynolds number, $\textrm{Re}\mathscr{h}^2 \ll 1$ ; (iii) thermophysical properties are constant; (iv) the solid–liquid interface is flat, i.e. $ h=h(t)$ ; (v) the hydrostatic pressure in the liquid film is neglected ( $p_h \ll 1$ ); (vi) heat convection is negligible $\textrm{Pe}\mathscr{h}^2 \ll 1$ when $\textrm{St} \leqslant 0.1$ for CCM (Hu & Fan 2023; Ezra & Kozak 2024); (vii) viscous dissipation is negligible due to $ \textrm{PrEc} \ll 1$ for most PCMs and thermal conditions; (viii) the initial temperature of the PCM solid remains $T_m$ ; (ix) the heat flux across the gas–liquid interface is negligible due to the much low thermal conductivity of the gas compared with the liquid. Additionally, we consider (x) the initial geometrical conditions $ \mathscr{l}^2 \ll 1$ and $\mathscr{h} \ll A$ . Consequently, with $l =l^*/h_0^*$ providing the corresponding transverse dimensions of the liquid flow channel, the governing equations can be simplified considerably to

(2.4a) \begin{align}&\qquad \partial _Xu+\partial _Yv=0, \\[-9pt] \nonumber \end{align}

(2.4b) \begin{align}& 0= -\partial _XP+\partial _Y^2u+ l^{-2}\partial _Z^2u, \\[-9pt] \nonumber \end{align}

(2.4c) \begin{align}&\qquad\quad\,\, 0= \partial _YP, \\[-9pt] \nonumber \end{align}

(2.4d) \begin{align}&\qquad\quad\,\, 0 = \partial _Z P, \\[-9pt] \nonumber \end{align}

(2.4e) \begin{align}&\qquad 0 = \partial _Y^2 T+l^{-2}\partial _Z^2T, \\[-9pt] \nonumber \end{align}

(2.4f) \begin{align}&\qquad\, \partial _YT(X,h)= \frac {d\Delta }{d \tau }, \end{align}

where $\tau \equiv t\textrm{St}/ [\rho \textrm{Pe}\mathscr{h}(A-\mathscr{h}) ]$ is a rescaled time, which recognises that a balance of heat conduction and the heat of fusion control the dynamics. Additionally, the force balance on the remaining solid is satisfied as

(2.5) \begin{equation} \int _{-1/2}^{1/2} P {\textrm d}X = \varDelta ^{\mathscr{c}}, \end{equation}

where $\mathscr{c} = 0$ represents the constant-pressure CCM mode and $\mathscr{c} = 1$ represents a gravity-driven configuration. Then, the boundary conditions for the velocity and temperature fields within a single periodic unit cell are given as

(2.6a) \begin{align}& u(Y=h)=0,\ \partial _Yu(Z=\pm 1/2)=0, \ \boldsymbol{n}\cdot \nabla u(Y=0,|Z|\lt \phi /2)=0,\ u(Y=0,\notag\\&\quad |Z|\gt \phi /2)=0, \end{align}

(2.6b) \begin{align}&\qquad\qquad\qquad\quad v(Y=0)=0,\ v(Y=h) = \rho A\mathscr{h}^{-1}\partial _{t} H - \partial _{ t} h, \end{align}

(2.6c) \begin{align}& T(Y=h)=0,\ \partial _YT(Z=\pm 1/2)=0, \ \boldsymbol{n}\cdot \nabla T(Y=0,|Z|\lt \phi /2)=0,\ T(Y=0,\notag\\&\quad |Z|\gt \phi /2)=1. \end{align}

The goal of the subsequent parts of § 2 is to determine the flow and temperature fields so as to find how the PCM height changes with time, i.e. $\varDelta (\tau )$ or $H(\tau )$ . A significant computational difficulty is to do this calculation for different groove configurations, i.e. gas–liquid fractions $\phi$ and dimensionless periodic lengths $l$ and liquid film thicknesses $\Lambda$ . The specific results are given in § 3.

2.2. Average heat flux and effective thermal slip length

In this section, a dual-series approach similar to previous work (Lauga & Stone 2003; Sbragaglia & Prosperetti 2007; Kirk et al. 2017) was adopted to give exact solutions for an effective thermal slip length $\lambda _{t}^{*}$ in the two-dimensional configuration. For convenience and consistency with previous work, the characteristic length $l^{*}$ of the periodic topography of the substrate is used for the length scale throughout this section, leading to dimensionless coordinates $(x,y,z)\equiv (x^*,y^*,z^*)/l^*$ and channel height $\Lambda \equiv h^*/l^*$ . Consequently, the schematic diagram of the configuration is depicted in figure 1(b). The dimensionless variables $X$ and $Y$ introduced in the previous subsection will return when we obtain a one-dimensional model of the problem in § 2.4.

In addition to the temperature distribution, our primary interest lies in the average heat flux $\left \langle q^{\prime \prime }\right \rangle \equiv \int _{-1/2}^{1/2} -\partial _yT(\Lambda, z) dz$ through the upper boundary of the liquid film, as this parameter determines the melting rate as indicated by (2.4f ). Consequently, the introduction of an effective thermal slip length $\lambda _t \equiv \lambda _t^*/l^*$ can reduce the two-dimensional problem $T(y,z)$ in the $y-z$ plane to a one-dimensional problem $\bar T(y)$ in the $y$ direction, as illustrated in the schematic of the thermal transport in figure 1(b). In the following, $\lambda _{t,f}$ and $\lambda _{t,c}$ represent, respectively, the thermal slip lengths of a flat and curved gas–liquid interface.

2.2.1. On the flat interface, $\lambda _{t,f}$

Based on (2.4e ), the dimensionless form of the energy equation is rewritten as

(2.7) \begin{equation} \partial _{yy} T + \partial _{zz} T = 0. \end{equation}

As depicted in figure 1(b), the boundary conditions are

(2.8a) \begin{align}&\qquad\qquad {T}(\Lambda, z)=0, \end{align}

(2.8b) \begin{align}& \begin{cases} {T}(0,z)=1, & \phi /2 \lt |z| \leqslant 1/2 \\ \partial _y {T}(0,z)= 0, & |z| \leqslant \phi /2 \end{cases} \end{align}

(2.8c) \begin{align}&\,\, \partial _z T(y,-1/2)= \partial _zT(y,1/2) = 0, \end{align}

we assume that the heat flux vanishes at the gas–liquid boundary owing to the smallness of the ratio of the gas and liquid thermal conductivities. We can observe that the boundary conditions for this two-dimensional problem introduce the fraction of the boundary that involves the gas–liquid fraction $\phi$ . Due to the linearity of the boundary-value problem, the temperature $T$ can be divided into two components $T_b$ and $T_s$ , namely $T(y,z)=T_b(y)+T_s(y,z)$ , where $T_b(y) = 1 - y/\Lambda$ satisfies the boundary conditions on temperature in the $y$ direction, and $T_s(y,z)=\tilde {T}(y,z)/\Lambda$ accounts for the influence of the hydrophobic surface. Hence, $\tilde {T}$ also satisfies the Laplace equation as

(2.9) \begin{equation} \partial _{yy} \tilde {T} + \partial _{zz} \tilde {T} = 0. \end{equation}

The boundary conditions for $\tilde {T}$ are

(2.10) \begin{align}&\qquad\qquad\qquad \tilde {T}(\Lambda, z)=0, \end{align}

(2.11) \begin{align}& \begin{cases} \tilde {T}(0,z)=0, & \phi /2 \lt |z| \leqslant 1/2; \\ -1 + \partial _y\tilde {T}(0,z) = 0, & |z| \leqslant \phi /2. \end{cases} \end{align}

The general solution $\tilde {T}$ , accounting for condition (2.8c ), is

(2.12) \begin{equation} \begin{aligned} \tilde {T} (y, z) = c_0 \frac {y}{\Lambda } + d_0 + \sum _{n=1}^\infty \left [c_n \cosh \left (2\pi n y\right ) + d_n \sinh \left (2\pi n y\right )\right ] \cos \left (2\pi n z\right ), \end{aligned} \end{equation}

where $c_n$ and $d_n$ for $n \in [0, \infty )$ are constants to be determined for given $\phi$ and $\Lambda$ . Following the standard procedure of a dual-series method to numerically solve (2.12) as detailed in Appendix A, $c_n$ and $d_n$ for $n \in [0, N-1]$ can be obtained with a numerical truncation for large enough $N$ . The choice of $N$ should be chosen to ensure convergence of the numerical results (Teo & Khoo 2008). Here, we adopt $N=500$ as there is no significant deviation with $N=1000$ , which is also consistent with Landel et al. (2020). Therefore, the average heat flux $\left \langle q^{\prime \prime }\right \rangle$ can be calculated by substituting $c_0$ as

(2.13) \begin{equation} \left \langle q^{\prime \prime }\right \rangle = -\int _{-1/2}^{1/2} \partial _yT(\Lambda, z) dz = \frac {1}{\Lambda } - \frac {c_0}{\Lambda ^2}. \end{equation}

The definition of an effective thermal slip length $ \lambda _{t}$ representative of the composite flat liquid–gas surface and solid–liquid surface (Enright et al. 2014) is

(2.14) \begin{equation} \lambda _{t} \equiv \frac {T_b(y=0)- \left \langle T(y=0)\right \rangle }{\left \langle \partial _y T(y=0)\right \rangle } = \frac {1-\int _{-1/2}^{1/2}\left (1+\tilde {T}(y=0)/\Lambda \right ){\textrm d} z}{\int _{-1/2}^{1/2} \partial _yT(0, z) {\textrm d}z}. \end{equation}

With the results of $1-\int _{-1/2}^{1/2}(1+\tilde {T}(y=0)/\Lambda ) {\textrm d}z= c_0/\Lambda$ due to $c_0 = -d_0$ from Appendix A and $\int _{-1/2}^{1/2} \partial _yT(0, z) {\textrm d}z= \left \langle q^{\prime \prime }\right \rangle$ , we find the effective thermal slip length $ \lambda _{t,f}$ for a flat, composite gas–liquid interface:

(2.15) \begin{equation} \begin{aligned} \lambda _{t,f} = \frac {{\Lambda }c_0}{\Lambda -c_0}. \end{aligned} \end{equation}

2.2.2. On the curved interface, $\lambda _{t,c}$

We next consider how the deformation of the gas–liquid interface affects the thermal slip length. A curved liquid–air interface is formed by the pressure difference between liquid $ P_l^{*}$ and gas $ P_g^{*}$ , and described approximately by the radius of curvature $ R^{*} = \sigma ^{*}/(P_l^{*} - P_g^{*})$ , where $\sigma ^{*}$ is the surface tension. Hence, the shape of the meniscus can be described by the geometric relationship

(2.16) \begin{equation} y = \sqrt {R^2 - \phi ^2 / 4} - \sqrt {R^2 - z^2}, \end{equation}

where $R = R^{*} / l^{*}$ . For $R \gg \phi / 2$ , this expression simplifies to

(2.17) \begin{equation} y = \frac {1}{8R} \left (4z^2 - \phi ^2\right ) + {O}\left (\frac {1}{R^3}\right ) = - \epsilon \eta (z) \quad \text{for} \quad |z| \leqslant \phi / 2, \end{equation}

after defining $\epsilon = 1 / (8R)$ , which will be assumed $\ll 1$ , and $\eta (z) = \phi ^2 - 4z^2$ . It should be noted that the protrusion angle $\theta$ satisfies $R \sin \theta = \phi / 2$ , which leads to a maximum $\epsilon = \sin \theta /4\phi = 0.217$ for $\phi = 0.2$ and $\theta = 10^{\circ }$ . Again, neglecting the influence of the conductivity of the gas phase, then we have the condition of no heat flux across the meniscus, and the temperature field should satisfy

(2.18) \begin{equation} \boldsymbol{n} \cdot \boldsymbol{\nabla }T = 0 \quad \text{at} \quad y =-\epsilon \eta (z), \end{equation}

which is equivalent to

(2.19) \begin{equation} \partial _y T(-\epsilon \eta, z)+\epsilon \eta ^{\prime } \partial _z T(-\epsilon \eta, z)=0, \end{equation}

where $\eta ^{\prime } \equiv \textrm{d} \eta /\textrm{d}z =-8z$ . Since $\epsilon \ll 1$ , $T(-\epsilon \eta, z)$ can be obtained from $T(0, z)$ by Taylor expansion as $T(-\epsilon \eta, z)=T(0, z)-\epsilon \eta \partial _y T(0, z)+O (\epsilon ^2 )$ , leading to

(2.20) \begin{equation} \begin{aligned} & \partial _y T(-\epsilon \eta, z)=\partial _y T(0, z)-\epsilon \eta \partial _{yy} T(0, z)+O(\epsilon ^2), \\ & \partial _z T(-\epsilon \eta, z)=\partial _z T(0, z)+O(\epsilon ). \end{aligned} \end{equation}

Then, substitution of (2.20) into (2.19) yields

(2.21) \begin{equation} \partial _y T(0, z)-\epsilon \eta \partial _{yy} T(0, z)+\epsilon \eta ^{\prime } \partial _z T(0, z)+O(\epsilon ^2)=0. \end{equation}

Hence, the boundary conditions are

(2.22) \begin{equation} \begin{aligned} & T(\Lambda, z)=0 \\ & \left \{\begin{array}{ll} T(0, z)=0, & \phi /2\lt |z| \leqslant 1 / 2; \\ \partial _y T(0, z)-\epsilon \eta \partial _{y y} T(0, z)+\epsilon \eta ^{\prime } \partial _z T(0, z)=0, & |z| \leqslant \phi /2. \end{array} \right . \end{aligned} \end{equation}

After a perturbation analysis by substituting $T=T^{(0)}+\epsilon T^{(1)}+O (\epsilon ^2 )$ (details provided in Appendix B), we find that the presence of a curved gas–liquid interface has no influence on the average heat flux and thermal slip length, which means

(2.23) \begin{equation} \lambda _{t,f}=\lambda _{t,c}\equiv \lambda _t. \end{equation}

We have now completed the necessary aspects of the thermal problem inside the thin film to estimate the time variations of the PCM, which will be described in § 2.4.

2.3. Flow rate – pressure gradient relationship and velocity slip length

We now consider the flow in the thin film since we need the pressure distribution to complete the force balance on the PCM. Based on (2.4b )–(2.4d ), the main liquid flow is along the $x$ direction as indicated in figure 1(d). Similar to § 2.2, in addition to obtaining the velocity field $u$ , we also focus on the variation of the flow rate-pressure gradient relationship $Q-\partial _x P$ due to the presence of trapped gas, which will influence the pressure distribution within the liquid film and the liquid film thickness during the CCM process. By introducing the velocity slip length $\lambda \equiv \lambda ^*/l^*$ to enable an equivalent flow rate $\bar Q = Q$ as demonstrated in figure 1(e), these steps allow a dimensional reduction from $u(y,z)$ to $\bar u(y)$ . In the following, $\lambda _{\parallel, f}$ and $\lambda _{\parallel, c}$ represent the velocity slip lengths by flat and curved gas–liquid interfaces, respectively, while the main liquid flow is parallel to the orientation of the grooves.

2.3.1. Velocity slip length on the flat interface, $\lambda _{\parallel, f}$

With the definition of variables introduced in the previous section, the dimensionless form of the momentum equation is

(2.24) \begin{equation} \partial _{yy}u+\partial _{zz}u=\partial _xP, \end{equation}

where $u={u}^{*} \mu ^{*} L^{*}/({{l}^{*}}^2p_c^{*})$ . Due to linearity, the velocity $u(y,z)$ can be divided into two components $u_b(y)$ and $u_{s}(y,z)$ , where $u_b(y) = -\partial _xP (\Lambda y-y^2 )/2$ accounts for the flow without slip and $u_s(y,z)=-\partial _xP\Lambda \tilde {u}/2$ accounts for slip along the hydrophobic surface. Hence, substituting $u(y,z)= -\partial _xP (\Lambda y-y^2+\Lambda \tilde {u}(y,z) )/2$ leads to a boundary-value problem

(2.25) \begin{equation} \partial _{yy} \tilde {u}+\partial _{zz} \tilde {u}=0, \end{equation}

with corresponding homogeneous boundary conditions

(2.26) \begin{align}&\qquad\qquad\quad \tilde {u}(\Lambda, z)=0, \end{align}

(2.27) \begin{align}& \begin{cases}\tilde {u}(0, z)=0, & \phi /2 \lt |z| \leqslant 1/2; \\ 1+\partial _{y} \tilde {u}(0, z)=0, & |z|\leqslant \phi /2.\end{cases} \end{align}

This problem is effectively identical to the thermal problem in § 2.2.1. The general solution of $\tilde {u}(y,z)$ with periodic boundary conditions at $z = \pm 1/2$ is

(2.28) \begin{equation} \begin{aligned} \tilde {u} (y,z) = r_0 +q_{0}\frac {y}{\Lambda }+\sum _{n=1}^\infty \left [r_n\cosh \left (2\pi n y\right )+q_n \sinh \left (2\pi n y\right )\right ] \cos \left (2\pi n z\right ), \end{aligned} \end{equation}

where constants $r_n$ and $q_n$ ( $n \in [0, N-1]$ ), dependent on $\phi$ and $\Lambda$ , also can be numerically determined similar to the dual-series approach in Appendix C. Consequently, the increased flow rate $Q_d$ due to the presence of slip can be obtained via

(2.29) \begin{equation} Q_d = \int _{0}^{\Lambda }\int _{-1/2}^{1/2} u_s \textrm{d}z\textrm{d}y =\int _{0}^{\Lambda }\int _{-1/2}^{1/2} -\frac {{\textrm d}P}{{\textrm d}x}\frac {\Lambda \tilde {u}}{2} \textrm{d}z\textrm{d}y= -\frac {{\textrm d}P}{{\textrm d}x}\frac {{\Lambda }^2r_0}{4}. \end{equation}

On the other hand, for a model pressure-driven channel flow with one boundary of slip length $\lambda$ , the increased flow rate $Q_d$ can be written as

(2.30) \begin{equation} Q_d = -\frac {{\textrm d}P}{{\textrm d}x}\frac {{\Lambda ^3\lambda }}{4\Lambda +4\lambda }. \end{equation}

Comparing the last two expressions, the effective velocity slip length $\lambda _{\parallel, f}$ on the longitudinal grooves for a flat liquid–gas interface is

(2.31) \begin{equation} \lambda _{\parallel, f} = \frac {\Lambda r_0}{\Lambda -r_0}. \end{equation}

A reader should notice that the detailed boundary value problems for the thermal and velocity fields, and corresponding dual integral equations solutions, show that the two problems are identical, $r_n=-c_n$ and $r_0=c_0$ . Thus, $\lambda _{\parallel, f}=\lambda _{t,f}$ as evident from (2.31) and (2.15).

2.3.2. Velocity slip length on the curved interface, $\lambda _{\parallel, c}$

To account for a curved gas–liquid interface, we consider the configuration similar to § 2.2.2. In order to satisfy the condition of zero shear stress at the steady-state gas–liquid interface, the velocity should satisfy

(2.32) \begin{equation} \boldsymbol{n} \cdot \boldsymbol{\nabla }u = 0 \quad \text{at} \quad y =-\epsilon \eta, \end{equation}

which is equivalent to

(2.33) \begin{equation} \epsilon \eta ^{\prime } \partial _z u(-\epsilon \eta, z)+\partial _y u(-\epsilon \eta, z)=0. \end{equation}

Since $\epsilon \ll 1$ , $u(-\epsilon \eta, z)$ can be obtained from $u(y, z)$ by Taylor expansion as $u(-\epsilon \eta, z)=u(0, z)-\epsilon \eta \partial _y u(0, z)+O (\epsilon ^2 )$ as

(2.34a) \begin{align}& \partial _y u(-\epsilon \eta, z)=\partial _y u(0, z)-\epsilon \eta \partial _{yy} u(0, z)+O(\epsilon ^2), \end{align}

(2.34b) \begin{align}&\qquad\quad \partial _z u(-\epsilon \eta, z)=\partial _z u(0, z)+O(\epsilon ). \end{align}

Then, substitution of (2.34a ) and (2.34b ) into (2.33) yields

(2.35) \begin{equation} \partial _y u(0, z)-\epsilon \eta \partial _{yy} u(0, z)+\epsilon \eta ^{\prime } \partial _z u(0, z)+O(\epsilon ^2)=0. \end{equation}

Hence, the boundary conditions are

(2.36) \begin{equation} \begin{aligned} & u=0 \quad \text{ at } \quad y=\Lambda \\ & \left \{\begin{array}{ll} u=0, & \phi /2\lt |z| \leqslant 1 / 2 \\ \partial _y u-\epsilon \left ( \eta \partial _{y y} u+ \eta ^{\prime } \partial _z u\right )=0, & |z| \leqslant \phi /2, \end{array} \quad\text{ at } y=0\right . . \end{aligned} \end{equation}

After substituting the perturbation expansion for the velocity field, $u=u^{(0)}+\epsilon u^{(1)}+O (\epsilon ^2 )$ , we solve for the first-order velocity $\tilde {u}^{(1)}$ (see Appendix D) as

(2.37) \begin{equation} \tilde {u}^{(1)}=r_0^{(1)} +q_0^{(1)}\frac {y}{\Lambda }+\sum _{n=1}^{\infty }\left [r_n^{(1)} \cosh \left (2\pi n y\right )+q_n^{(1)} \sinh \left (2\pi n y\right )\right ] \cos \left (2\pi n z\right )\!; \end{equation}

the coefficients $r_0^{(0)}$ , $r_n^{(0)}$ , $r_0^{(1)}$ , and $r_n^{(1)}$ can be obtained numerically. Consequently, the total flow rate $Q$ can be calculated as

(2.38) \begin{equation} \begin{aligned} Q & =\int _{-1/2}^{ 1/2} \int _{0}^{\Lambda } u^{(0)} \textrm{d} y \textrm{\ d} z+\epsilon \left (\int _{-1/2}^{ 1/2} \int _{0}^{\Lambda } u^{(1)} \textrm{d} y \textrm{\ d} z+\int _{-\phi /2}^{\phi /2} u^{(0)}(0, z) \eta \textrm{\ d} z\right )+O(\epsilon ^2) \\ &=Q^{(0)}+\epsilon Q^{(1)}+O(\epsilon ^2), \end{aligned} \end{equation}

where

(2.39a) \begin{align}&\qquad\qquad\qquad\qquad\quad Q^{(0)} = -\frac {dP}{dx}\left [\frac {{\Lambda }^3}{12}+\frac {{\Lambda }^2r_0^{(0)}}{4}\right ], \end{align}

(2.39b) \begin{align}& Q^{(1)} = -\frac {dP}{dx}\left [\frac {{\Lambda }^2r_0^{(1)}}{4}+ \frac {\phi ^3\Lambda }{3}r_0^{(0)}+\Lambda \sum _{n=1}^{\infty } r_n^{(0)} \frac {\sin (n\pi \phi )-n\pi \phi \cos (n\pi \phi )}{n^3\pi ^3} \right ]. \end{align}

Based on $Q=-\partial _xP\Lambda ^3/12+Q_d = -\partial _xP\Lambda ^3/12 -\partial _xP\Lambda ^3\lambda /(4\Lambda +4\lambda )$ , the slip length $\lambda _{\parallel, c}$ along the curved interface is

(2.40) \begin{equation} \begin{aligned} \lambda _{\parallel, c} &= \lambda _{\parallel }^{(0)} + \epsilon \lambda _{\parallel }^{(1)} = \frac {r_0^{(0)}\Lambda }{\Lambda - r_0^{(0)}} + \epsilon \frac {4\left (\Lambda + \lambda _{\parallel }^{(0)}\right )^2}{{\Lambda }^4}\frac {Q^{(1)}}{-\partial _X P} = \frac {r_0^{(0)}\Lambda }{\Lambda - r_0^{(0)}} \\ & \quad + \epsilon \frac {4\left (\Lambda + \lambda _{\parallel }^{(0)}\right )^2}{{\Lambda }^4}\!\left [ \! \frac {{\Lambda }^2 r_0^{(1)}}{4} + \frac {\phi ^3 \Lambda }{3} r_0^{(0)} + \Lambda \! \sum _{n=1}^{\infty } r_n^{(0)} \frac {\sin (n\pi \phi )-n\pi \phi \cos (n\pi \phi )}{n^3\pi ^3} \! \right ]\!. \end{aligned} \end{equation}

2.4. A one-dimensional description of CCM with slip

The velocity slip length $\lambda \equiv \lambda ^*/l^*$ and temperature slip length $\lambda _{t}\equiv \lambda _t^*/l^*$ , describing average effects from the mixed boundary conditions of the two-dimensional problem, enable us to simplify the CCM problem into an effective one-dimensional problem describing flow and heat exchange along the direction of liquid flow, shown in the schematic diagram in figure 1(f), where the notation $ \bar {()}$ implies that the variables are averaged in the $z$ direction. The corresponding governing equations (2.4a ), (2.4b ), (2.4e ) and (2.4f ) can be replaced, respectively, by

(2.41a) \begin{align}&\,\, \partial _X \bar u +\partial _Y \bar v = 0, \end{align}

(2.41b) \begin{align}& 0=-\partial _X \bar P+\partial _Y^2 \bar u, \end{align}

(2.41c) \begin{align}&\qquad 0=\partial _Y^2 \bar T, \end{align}

(2.41d) \begin{align}& \partial _Y \bar T(X,h)= \frac {{\textrm d}\Delta }{{\textrm d} \tau } . \end{align}

Hence, the corresponding boundary conditions in figure 1(f) are

(2.42a) \begin{align}&\qquad\qquad\qquad \bar u(X,0) = l\lambda \partial _Y \bar u(X,0), \quad \bar u(X,h)=0, \end{align}

(2.42b) \begin{align}&\qquad\qquad 1-\bar T(X,0)=-l\lambda _{t}\partial _Y \bar T(X,0), \quad \bar T(X,h)=0, \end{align}

(2.42c) \begin{align}& \bar v(X,0)=0,\quad \bar v(X,h)=\left [ \left (A-{\mathscr{h}}\right )\frac {\rho }{\mathscr{h}}\frac {{\textrm d} \Delta }{{\textrm d} \tau }+\left (\rho -1\right )\frac {{\textrm d} h}{{\textrm d} \tau }\right ]\frac {\textrm{St}}{\rho \textrm{Pe}\mathscr{h}\left (A-\mathscr{h}\right )}, \end{align}

(2.42d) \begin{align}&\qquad\qquad\qquad\qquad \partial _X \bar P(0,Y)=0, \quad \bar P(1/2,Y)=0. \end{align}

It is worth noting that the upper velocity boundary $\bar v(X,h)$ is derived from $\bar v^*(x^*,h^*) = \rho \partial _{t^*} H^* - \partial _{ t^*} h^*$ based on mass conservation at the melting front.

The velocity profile $\bar u(X,Y)$ can be obtained by integrating (2.41b ) twice with respect to $Y$ and applying boundary conditions (2.42a ), yielding

(2.43) \begin{equation} \bar u = \frac {h^2}{2}\frac {\partial \bar P}{\partial X}\left (\frac {Y^2}{h^2}-\frac {Y/l+\lambda }{\Lambda +\lambda }\right ), \end{equation}

where $ h/l = \Lambda$ . Then, substituting (2.43) into the continuity equation (2.41a ) and integrating along $y$ from 0 to $h$ leads to

(2.44) \begin{equation} \frac {h^2}{2}\frac {\partial ^2 \bar P}{\partial X^2}\int _{0}^{h}\left (\frac {Y^2}{h^2}-\frac {Y/l+\lambda }{\Lambda +\lambda }\right ) {\textrm d}Y+ \underbrace {\int _{0}^{h}\frac {\partial v}{\partial Y} {\textrm d}Y}_{\bar v|_{Y=h}-\bar v|_{Y=0}}=0. \end{equation}

Noting that $\bar v(X,h) \approx \textrm{St} \partial _\tau \varDelta /(\textrm{Pe}\mathscr{h}^2)$ due to the typical property $\rho -1 \lt {O}(10^{-1})$ , then substituting (2.42c ) into (2.44) leads to

(2.45) \begin{equation} \frac {\partial ^2 \bar P}{\partial X^2} = \frac {d \Delta }{d\tau } \frac {\textrm{St}}{\textrm{Pe}\mathscr{h}^2}\frac {6}{h^3}\frac {\Lambda +\lambda }{\Lambda +4\lambda } . \end{equation}

The pressure distribution can be obtained by integrating (2.45) twice with respect to $x$ , and using conditions (2.42d ), to obtain

(2.46) \begin{equation} \bar P = \frac {d \Delta }{d\tau } \frac {\textrm{St}}{\textrm{Pe}\mathscr{h}^2}\frac {12}{h^3}\frac {\Lambda +\lambda }{\Lambda +4\lambda } \left (X^2-\frac {1}{4}\right ). \end{equation}

Then, we can substitute (2.46) into the force balance (2.5) to obtain a relationship between $\varDelta$ and $h$ :

(2.47) \begin{equation} \frac {d \Delta }{d\tau } = -\frac {\textrm{Pe}\mathscr{h}^2}{\textrm{St}}h^3\frac {\Lambda +4\lambda }{\Lambda +\lambda } \varDelta ^{\mathscr{c}}. \end{equation}

The temperature profile can be obtained by substituting (2.42b ) into (2.41c ), yielding

(2.48) \begin{equation} \bar T = \frac {\Lambda -Y/l}{\Lambda +\lambda _{t}}. \end{equation}

Substitution of (2.48) into (2.41d ) gives

(2.49) \begin{equation} \frac {d\Delta }{d\tau }=-\frac {1}{l(\Lambda +\lambda _{t})}. \end{equation}

The combination of (2.47) and (2.49) leads to

(2.50) \begin{equation} h^4\frac {\Lambda +4\lambda }{\Lambda +\lambda }\left (1+\frac {\lambda _{t}}{\Lambda }\right )=\frac {\textrm{St}}{\textrm{Pe} \mathscr{h}^2}\frac {1}{\varDelta ^{\mathscr{c}}}. \end{equation}

Equations (2.49) and (2.50) are coupled expressions between the solid height $\varDelta (\tau )$ that remains and the film thickness $h$ . Note that because of the definition $\Lambda =h/l$ where $l$ is fixed, $h$ varies with $\tau$ .

Recalling $\mathscr{h}=h_0^*/L^*$ and $h=h^*/h_0^*$ in (2.1) indicates that $h_0^*$ has not been specified. Usually, the initial height satisfies $ H_0^* \gg h^*$ (i.e. $\varDelta \approx H$ ) and we can further find a characteristic initial film thickness $h_0^*$ that satisfies the $\textrm{St}/\textrm{Pe}\mathscr{h}^2=1$ , corresponding to the initial solid height $H=1$ and no-slip surface $\lambda =\lambda _t=0$ . This enables us to define $h_0^*$ as

(2.51) \begin{equation} h_0^* = \left [\frac {(T_w^*-T_m^*)k_l^*{L^*}^2\mu ^*}{\mathcal{H}^*\rho _l^*p_c^*}\right ]^{\frac {1}{4}}, \end{equation}

from which (2.50) and (2.49) are transformed into the equations ( $\Lambda =h/l$ )

(2.52a) \begin{align}& h^4\frac {\Lambda +4\lambda }{\Lambda +\lambda }\left (1+\frac {\lambda _{t}}{\Lambda }\right ) = \frac {1}{H^{\mathscr{c}}}, \end{align}

(2.52b) \begin{align}&\quad \frac {{\textrm d}H}{{\textrm d}\tau } = -\frac {1}{l(\Lambda +\lambda _{t})}. \end{align}

Based on (2.52a ), we can conclude that a constant film thickness $h$ is maintained for the constant-pressure mode ( $\mathscr{c}=0$ ) though $\lambda$ and $\lambda _t$ are dependent on $h$ . On the other hand, for the gravity-driven mode ( $\mathscr{c}=1$ ) there is a coupling between $h$ and $H$ , making the solutions of $H(\tau )$ and $h(\tau )$ complicated.

Furthermore, we introduce the Nusselt number, $\textrm{Nu}$ , based on the instantaneous height of the liquid film $h^*(t^*)$ , to demonstrate the heat transfer capability under different conditions. $\textrm{Nu}$ is defined as

(2.53) \begin{equation} {\textrm{Nu}} \equiv -\frac {k_l^{*}\partial _{y^{*}} \bar T^{*}(0,h^{*})}{\left (T_w^{*}-T_m^{*}\right )}\frac {h^{*}}{k_{l}^{*}} = \frac {1}{l(\Lambda + \lambda _{t})}. \end{equation}

At this stage, we have completed all of the analysis of the thin-film flow and heat transfer necessary for modelling CCM processes, which may be enhanced using a corrugated substrate. The main variables that dictate performance are the slip lengths $\lambda$ and $\lambda _t$ , geometric ratios $l$ and $\Lambda$ , the latter which changes in time, and the resulting $\textrm{Nu}$ , which determines the melting rate, i.e. the power density of the CCM.

3. Results and discussions

The goal of this section is to calculate the melting rate (i.e. $\textrm{Nu}$ ) or the instantaneous solid height (i.e. $H(\tau )$ ), which is a complicated function of geometric parameters ( $l$ ) and slip lengths ( $\lambda$ , $\lambda _t$ ). It also depends on the operation controlled by hydrostatic pressure or a constant applied pressure. The subsections that follow summarise the main results for different cases.

3.1. Confinement and meniscus effects on the slip length

3.1.1. Velocity slip length $\lambda$ and temperature slip length $\lambda _t$ at a flat interface

Figure 2. Slip lengths as a function of $\Lambda$ . (a) Temperature slip length $\lambda _t$ , (b) velocity slip length $\lambda _{\parallel, f}$ on longitudinal grooves and (c) $\lambda _{\perp, f}$ on transverse grooves versus the aspect ratio $\Lambda =h^*/l^*$ when the liquid–gas interface is flat. (d) Comparison of the ratio of the velocity to temperature slip lengths, $\lambda _{,f}/\lambda _t$ , between longitudinal and transverse grooves. In figures (a)–(c) the solid lines are numerical results. Dashed lines and dot–dash lines are asymptotic solutions.

The temperature $\lambda _t$ (§ 2.2.1) and velocity $\lambda _{\parallel, f}$ slip lengths (§ 2.3.1) as a function of the aspect ratio $\Lambda =h^*(t^*)/l^*$ (see figure 1 a), for a given gas–liquid fraction $\phi$ , are presented as solid lines in figures 2(a) and 2(b), respectively. For comparison, we also calculate the slip length $\lambda _{\perp, f}$ variation on a surface of transverse grooves in Appendix F, and present it in figure 2(c). It can be observed that across a wide range of gas–liquid fractions $\phi$ , each slip length of $\lambda _t$ , $\lambda _{\parallel, f}$ and $\lambda _{\perp, f}$ maintains a specific fixed value for $\Lambda \gg 1$ , and exhibits a scaling law of $\sim \Lambda$ when $\Lambda \ll 1$ .

Asymptotic solutions for these slip lengths are readily available. These solutions are also provided in figure 2 where they are drawn as dashed lines for $\Lambda \gg 1$ and dot–dash lines for $\Lambda \ll 1$ . In the limit of $\Lambda \gg 1$ , the asymptotic solutions of $\lambda _{\parallel, f}$ and $\lambda _{\perp, f}$ are well known (Lauga & Stone 2003; Teo & Khoo 2008) as

(3.1a) \begin{align}&\,\, \lambda _{\parallel, f}= \frac {1}{\pi } \ln \left (\sec \left (\frac {\phi \pi }{2}\right )\right ), \quad \Lambda \gg 1, \end{align}

(3.1b) \begin{align}& \lambda _{\perp, f}= \frac {1}{2\pi } \ln \left (\sec \left (\frac {\phi \pi }{2}\right )\right ), \quad \Lambda \gg 1. \end{align}

In the limit of $\Lambda \ll 1$ , we found the asymptotic solutions proposed by Teo & Khoo (2008) are incorrect due to mistakes in calculating $Q_d$ . The velocity slip length $\lambda _f$ should be

(3.2a) \begin{align}&\,\,\,\,\, \lambda _{\parallel, f}= \frac {\phi }{1-\phi }\Lambda, \quad \Lambda \ll 1, \end{align}

(3.2b) \begin{align}& \lambda _{\perp, f}= \frac {\phi }{4\left (1-\phi \right )}\Lambda, \quad \Lambda \ll 1, \end{align}

which is consistent with the other related results (Lauga & Stone 2003; Landel et al. 2020). Furthermore, asymptotic solutions for the temperature slip length are obtained as

(3.3a) \begin{align}& \lambda _{t} = \frac {1}{\pi } \ln \left (\sec \left (\frac {\phi \pi }{2}\right )\right ), \quad \Lambda \gg 1, \end{align}

(3.3b) \begin{align}&\qquad \lambda _{t} = \frac {\phi }{1-\phi }\Lambda, \quad \Lambda \ll 1. \end{align}

As already highlighted at the end of § 2.3.1, based on the results reported in figures 2(a) and 2(b), it is observed that $\lambda _{\parallel, f}=\lambda _t$ , while $\lambda _{\perp, f}$ remains a constant value proportional to them in both asymptotic ranges. As shown in figure 2(d), the ratio $\lambda _{\perp, f}/\lambda _t$ varies between 0.25 and 0.5 depending on $\phi$ and $\Lambda$ , consistent with $\lambda _{\parallel, f} = \lambda _t = 2\lambda _{\perp, f} = \ln (\sec (\phi \pi /2) )/\pi$ for $\Lambda \gg 1$ , and $\lambda _{\parallel, f} = \lambda _t = 4\lambda _{\perp, f} = \phi \Lambda /(1-\phi )$ for $\Lambda \ll 1$ .

3.1.2. Velocity slip length $\lambda _{\parallel, c}$ at a curved gas–liquid interface

Figure 3. Slip lengths at a curved interface. (a–b) Variation of the first-order velocity slip length $\lambda ^{(1)}$ for different $\Lambda$ , where a dotted line denotes (3.4), a dashed line (3.9) and a dot–dash line (3.10). (c) Total velocity slip length $\lambda _{\parallel, c}=\lambda ^{(0)}+\epsilon \lambda ^{(1)}$ versus $\Lambda$ , with a dashed line denoting (3.12) and a dot–dash line (3.11). (d) The ratio of the total velocity and temperature slip lengths, $\lambda _{\parallel, c}/\lambda _t$ , versus $\Lambda$ .

We next consider the velocity slip length $\lambda _{\parallel, c}$ at a curved gas–liquid interface as a function of $\Lambda$ and $\phi$ . The semilog and log–log results of the first-order slip length $\lambda ^{(1)}$ , computed as described in § 2.3.2, are shown in figures 3(a) and 3(b), respectively. For different $\phi$ , $\lambda ^{(1)}$ exhibits two plateau intervals as $\Lambda \ll 1$ and $\Lambda \gg 1$ , with a decreasing trend as $\Lambda$ increases. Notably, when $\Lambda \gtrsim 2$ , $\lambda ^{(1)}$ becomes negative. This occurs because the increased area for flow caused by the meniscus is negligible, while there is a reduced velocity gradient at the boundary $y=0$ due to the meniscus.

The asymptotic solution of the first-order slip length $\lambda ^{(1)}$ , with ${O}(\Lambda ^{-1})$ error, has already been proposed (Sbragaglia & Prosperetti 2007; Kirk et al. 2017),

(3.4) \begin{equation} \lambda _{\|}^{(1)}= -\phi ^3\int _{0}^{1} \frac {\left [1-\cos \left (\phi \pi s\right )\right ]\left (1-s^2\right )}{\cos \left (\phi \pi s\right )-\cos \left (\phi \pi \right )}\textrm{d}s + {O}(\Lambda ^{-1}), \quad \Lambda \gg 1, \end{equation}

which only accounts for the changed velocity profile and neglects effects of increased flow rate associated with an expanded area, corresponding to the slip length

(3.5) \begin{equation} \lambda _{\|, c}=\ln (\sec (\phi \pi / 2)) / \pi -\epsilon \phi ^3 \mathcal{F}(\phi ),\quad \Lambda \gg 1; \end{equation}

the function $\mathcal{F}(\phi )$ is defined shortly. Based on (2.39b ) in the limit of $\Lambda \gg 1$ , we can give the flow rate $Q^{(1)}$ as

(3.6) \begin{equation} Q^{(1)}= -\frac {\phi ^3\Lambda ^2l}{4} \mathcal{F}(\phi ) +\pi \phi ^4\Lambda l^2 \mathcal{G}(\phi ),\quad \Lambda \gg 1, \end{equation}

where $\mathcal{F(\phi )}$ and $\mathcal{G(\phi )}$ are

(3.7) \begin{equation} \mathcal{F}(\phi )=\int _{0}^{1} \frac {\left [1-\cos \left (\phi \pi s\right )\right ]\left (1-s^2\right )}{\cos \left (\phi \pi s\right )-\cos \left (\phi \pi \right )}\textrm{d}s \end{equation}

and

(3.8) \begin{equation} \mathcal{G}(\phi )=\int _{0}^{1} \frac {s\left (1-s^2/3\right )\sin \left (\phi \pi s/2\right )}{\sqrt {\cos \left (\phi \pi s\right )-\cos \left (\phi \pi \right )}}\textrm{d}s. \end{equation}

These two integral functions can be referred to in previous work (Sneddon 1966; Sbragaglia & Prosperetti 2007). Hence, a modified asymptotic solution of $\lambda _{\parallel }^{(1)}$ for $\Lambda \gtrsim 0.2$ is proposed as

(3.9) \begin{equation} \lambda _{\|}^{(1)}=\left (-\phi ^3 \mathcal{F}(\phi ) +\frac {4}{\Lambda }\phi ^4 \mathcal{G}(\phi )\right )\left (1+ \frac {1}{4\Lambda \pi } \ln \left (\sec \left (\frac {\phi \pi }{2}\right )\right )\right )^2,\quad \Lambda \gtrsim 0.2. \end{equation}

In the limit of $\Lambda \ll 1$ , an asymptotic solution can be derived by substituting $\lambda _{\parallel }^{(0)}=l\phi \Lambda /(1-\phi )$ and $Q^{(1)} = 2\phi ^3\Lambda ^2l/3$ into the second term of (2.40), leading to

(3.10) \begin{equation} \lambda _{\parallel }^{(1)}=\frac {8\phi ^3}{3\left (1-\phi \right )^2},\quad \Lambda \ll 1. \end{equation}

As depicted in figures 3(a) and 3(b), the asymptotic solutions (3.4) and (3.10) for the plateaus can accurately predict the first-order velocity slip length $\lambda ^{(1)}$ for different $\phi$ . However, their applicability is limited ( $\Lambda \lesssim 0.01$ or $\Lambda \gtrsim 10^2$ ), while the newly proposed modified asymptotic solution (3.9) can accurately predict the first-order velocity slip for $\Lambda \gtrsim 0.2$ .

Based on the velocity slip length of zeroth-order $\lambda ^{(0)}$ and first-order $\lambda ^{(1)}$ , the slip length for a curved liquid surface, $\lambda _{\parallel, c}$ , can be finally obtained according to the equation $\lambda _{\parallel, c}=\lambda ^{(0)}+\epsilon \lambda ^{(1)}$ , which is plotted for different $\phi$ as solid lines in figure 3(c). The corresponding asymptotic solutions of $\lambda _{\parallel, c}$ are also drawn as dot–dash and dashed lines, corresponding, respectively, to the equations

(3.11) \begin{equation} \lambda _{\parallel, c}=\frac {\phi }{1-\phi }\Lambda + \epsilon \frac {8\phi ^3}{3\left (1-\phi \right )^2},\quad \Lambda \lesssim 0.01, \end{equation}

(3.12) \begin{align} \lambda _{\parallel, c}&=\frac {1}{\pi } \ln \left (\sec \left (\frac {\phi \pi }{2}\right )\right )\nonumber\\&\quad + \epsilon \left (-\phi ^3 \mathcal{F}(\phi ) +\frac {4}{\Lambda }\phi ^4 \mathcal{G}(\phi )\right )\left (1+ \frac {1}{4\Lambda \pi } \ln \left (\sec \left (\frac {\phi \pi }{2}\right )\right )\right )^2,\quad \Lambda \gtrsim 0.2. \end{align}

The presence of a meniscus leads to a constant value of $\lambda _{\parallel, c}$ when $\Lambda \ll 1$ . Hence, the results for $\lambda _{\parallel, c}/\lambda _t$ indicate potentially large values up to ${O}(10^2)$ compared with $\lambda _{\parallel, f}/\lambda _t=1$ , when $\Lambda \ll 1$ for arbitrary $\phi$ , as plotted in figure 3(d). The results imply that the existence of a meniscus may enhance CCM due to a significant effect of velocity slip compared with thermal slip.

The asymptotic solutions obtained in this section (§ 3.1) are summarised in table 1 for convenient comparison. The effects of gas–liquid fraction $\phi$ , aspect ratio $\Lambda$ and the presence of a meniscus on the velocity slip length $\lambda$ and temperature slip length $\lambda _t$ can be summarized as follows. (i) Increasing $\phi$ always leads to an increase in both $\lambda$ and $\lambda _t$ , as the proportion of the gas–liquid interfaces $\phi$ that allows slip increases. (ii) In the absence of deformation of the meniscus, when $\Lambda \gg 1$ , the slip lengths become independent of $\Lambda$ because for such aspect ratios there is no change on the local velocity or temperature distribution near the gas–liquid interface (figure 8 g–i); under such conditions, a no-slip top boundary is effectively equivalent to a free-shear or thermally insulating one. In contrast, when $\Lambda \ll 1$ , the no-slip condition at the top boundary significantly affects the flow or thermal field near the gas–liquid interface (figure 8 a–f), causing both $\lambda$ and $\lambda _t$ to decrease and approach zero as $\Lambda \to 0$ . (iii) When deformation of the meniscus is included, the temperature distribution is modified, but since the meniscus is assumed to be adiabatic, it does not contribute to the net heat flux, and the effective temperature slip length $\lambda _t$ remains unchanged. (iv) However, for velocity slip, the meniscus provides an additional flow area, and this contribution becomes more pronounced as $\Lambda$ decreases. Consequently, $\lambda$ does not vanish as $\Lambda \to 0$ but instead approaches a finite value that depends solely on $\phi$ .

3.2. Characterizing constant-pressure CCM

We now use the main results of § 3.1 to calculate the melting rate. Recalling (2.52a ) with $\mathscr{c}=0$ , we observe that a constant dimensionless liquid film thickness $h=h^*/h^*_0$ can be determined due to a constant exerted pressure, yielding

(3.13) \begin{equation} h = \left [\frac {1+ \lambda /\Lambda }{\left (1+4 \lambda /\Lambda \right )\left (1+ \lambda _{t}/\Lambda \right )}\right ]^{\frac {1}{4}}. \end{equation}

Substitution of (3.13) into the definition of $\textrm{Nu}$ (2.53) leads to a prediction for the melting rate

(3.14) \begin{equation} {\textrm{Nu}} = \left [\frac {1+4\lambda /\Lambda }{\left (1+ \lambda /\Lambda \right )\left (1+ \lambda _{t}/\Lambda \right )^3}\right ]^{\frac {1}{4}}. \end{equation}

Note that ${\textrm{Nu}} = 1$ represents the heat transfer performance of a smooth, unstructured surface ( $\lambda =\lambda _t=0$ ). Therefore, ${\textrm{Nu}} \gt 1$ indicates that the effects of velocity slip and temperature slip enhance heat transfer. Conversely, ${\textrm{Nu}} \lt 1$ implies a reduction in heat transfer.

3.2.1. Flat liquid–gas interface

Figure 4. Variation of $\textrm{Nu}_f$ along with (a) aspect ratio $\Lambda$ by numerical results and (b) liquid–gas fraction $\phi$ by asymptotic formula for flat gas–liquid interface. Here $\parallel$ and $\perp$ represent longitudinal and transverse grooves, respectively.

Substituting the numerical results of $\lambda _{\parallel, f}$ and $\lambda _{\perp, f}$ , respectively, obtained from § 2.2.1 and Appendix F into (3.14) yields the results of $\textrm{Nu}_f$ versus $\Lambda$ for a flat interface shown in figure 4(a). For both longitudinal and transverse groove structures, the existence of slip always reduces heat transfer. When $\Lambda \gg 1$ , $\textrm{Nu}_f$ asymptotically equals 1 due to negligible slip effect. For the case of $\Lambda \ll 1$ , $\textrm{Nu}_f$ decreases and then asymptotically equals a constant value for either the longitudinal or transverse grooves. It is obvious that $\textrm{Nu}_{\parallel, f} \gt \textrm{Nu}_{\perp, f}$ is valid for any $\phi$ . Therefore, we can determine that $\textrm{Nu}_f$ varies only as a function of $\phi$ within the asymptotic region ( $\Lambda \ll 1$ as indicated in figure 2). By substituting the asymptotic solutions from § 3.1.1, we find

(3.15) \begin{align}& \textrm{Nu}_{\parallel, f}=\left [\left (1+3\phi \right )\left (1-\phi \right )^3\right ]^{\frac {1}{4}}, \quad \Lambda \ll 1; \end{align}

(3.16) \begin{align}&\quad\,\, \textrm{Nu}_{\perp, f}=\left [\frac {4\left (1-\phi \right )^3}{4-3\phi }\right ]^{\frac {1}{4}},\quad \Lambda \ll 1, \end{align}

which are plotted in figure 4(b). All curves generally show a decreasing trend of $\textrm{Nu}$ as $\phi$ increases. For $\Lambda \ll 1$ , the parallel and perpendicular values of $\textrm{Nu}$ start at 1 and decrease with increasing liquid–gas fraction $\phi$ . For larger values of $\Lambda =$ 0.2 and 1, the curves for both $\textrm{Nu}_{\parallel, f}$ and $\textrm{Nu}_{\perp, f}$ exhibit more pronounced decreases and diverge more significantly as $\phi$ approaches 1. The curves suggest that the parameter $\Lambda$ significantly impacts $\textrm{Nu}_f$ , with higher values of $\Lambda$ leading to more rapid decreases in $\textrm{Nu}_f$ .

3.2.2. A curved liquid–gas interface on longitudinal grooves

Figure 5. Variation of (a) $\lambda /\Lambda$ and (b) $\lambda _t/\Lambda$ versus $\Lambda$ , where solid lines represent numerical results and dot–dash lines represent asymptotic solutions listed in table 1. (c) Variation of $\textrm{Nu}$ versus $\Lambda$ for curved gas–liquid interface, dot–dash lines are asymptotic solutions (3.17). (d) Map of $\textrm{Nu}$ at the various combinations of $\Lambda$ , $\phi$ and transverse/longitudinal surface structures for constant-pressure mode. Symbols denote , $\phi =0.2$ and $\Lambda \gg 1$ ; , $\phi =0.9$ and $\Lambda \gg 1$ ; , $\phi =0.2$ and $\Lambda \ll 1$ ; , $\phi =0.9$ and $\Lambda \ll 1$ .

For a curved liquid–gas interface on longitudinal grooves, first we observe that $\lambda _{\parallel, c}/\Lambda$ and $\lambda _t/\Lambda$ can be considered as the functions of $\Lambda$ as shown in figures 5(a) and 5(b), respectively. Similarly, by substituting the numerical solutions for $\lambda _{\parallel, c}$ and $\lambda _{t}$ into (3.14), we can obtain the numerical results of $\textrm{Nu}_c$ for considering a meniscus on longitudinal grooves, as plotted in figure 5(c). Compared with the flat interface, $\textrm{Nu}_c$ is also close to 1 for $\Lambda \gg 1$ . However, as $\Lambda \to 0$ , there exists a critical $\phi$ at which $\textrm{Nu}_c$ can be larger than 1 and heat transfer is enhanced. By substituting the asymptotic slip length solution (3.11) at $\Lambda \ll 1$ into (3.14), we find

(3.17) \begin{equation} \textrm{Nu}_{c}=\left (\frac {(1-\phi )^3\left (3 \Lambda \left (1+2 \phi -3 \phi ^2\right )+8 \phi ^2 \sin {\theta }\right )}{3 \Lambda (1-\phi )+2 \phi ^2 \sin {\theta }}\right )^{\frac {1}{4}} \overset {\Lambda \rightarrow 0}{=}\left [4(1-\phi )^{3}\right ]^{\frac {1}{4}}, \end{equation}

which is plotted as dot–dash lines for different $\phi$ in figure 5(c). Then, by setting ${\textrm{Nu}}_c =1$ , we can obtain the critical fraction $\phi _{cr} = 1-2^{-2/3} \approx 0.37$ . Therefore, we infer that when $\phi \gt 0.37$ , heat transfer will be reduced regardless of the value of $\Lambda$ . However, when $\phi \lt 0.37$ , heat transfer can be enhanced to accelerate melting, with $\textrm{Nu}_c$ approaching a constant value of $4^{1/4}(1-\phi )^{3/4}$ as $\Lambda$ decreases.

In order to visually compare the effects of surface structure, a meniscus, aspect ratio ( $\Lambda$ ) and gas–liquid ratio ( $\phi$ ) on $\textrm{Nu}$ , a two-dimensional phase diagram is presented in figure 5(d). In this diagram, the contours correspond to varying $\textrm{Nu}$ , with white dashed lines representing $\textrm{Nu}$ isoclines and other lines indicating results derived from asymptotic solutions. It can be observed that for a flat interface, ${\textrm{Nu}}\lt 1$ always remains for surfaces with transverse groove structures, indicating a reduction in heat transfer effectiveness. In contrast, surfaces with longitudinal grooves maintain ${\textrm{Nu}} \rightarrow 1$ for small $\phi$ values, but $\textrm{Nu}$ decreases monotonically as $\phi$ increases. When a meniscus (curved interface) is introduced, ${\textrm{Nu}}\gt 1$ can be achieved for longitudinal grooves initially at small $\phi$ and $\Lambda \lt 1$ , while $\textrm{Nu}$ shows a non-monotonic trend versus $\phi$ . The enhancement of $\textrm{Nu}$ at $\phi \lt 0.37$ is more pronounced as $\Lambda$ decreases.

Overall, under the condition of constant pressure, longitudinal grooves ( $\parallel$ ) exhibit the best heat transfer performance among all surfaces. Furthermore, only with the presence of a meniscus and significant confinement effects ( $\Lambda \ll 1$ ), can one achieve an enhancement in heat transfer at an appropriate gas–liquid fraction ( $\phi \lt 0.37$ ).

3.3. Gravity-driven CCM

In contrast with the constant pressure mode just discussed, when the PCM is supported only by hydrostatic pressure, i.e. the gravity-driven mode (i.e. $\mathscr{c}=1$ ), the film thickness $h^*(t^*)$ changes in time due to the continuously decreasing size of the PCM solid.

3.3.1. Numerical results of solid height $H(\tau )$ and film thickness $h(\tau )$

Figure 6. Variation of $H$ and $h$ versus $\tau$ for different conditions of (a) $l=10^3$ , (b) $l=10^2$ and (c) $l=1$ . In all figures, the black line represents the remaining height $H$ , the blue line represents the film thickness $h$ , the red dotted line () represents the magnitude of $l$ , red dashed line () and green dot–dash line() represent the no-slip solution of $H(\tau )$ (3.19a ) and $h(\tau )$ (3.19b ), respectively. Four values of $\phi = \{0.1, 0.2, 0.3, 0.5\}$ are chosen for calculating each case. Experimental data (black squares) are replotted (Moallemi et al. 1986) for comparison with the case $\phi = 0$ .

Recalling (2.52a ) and (2.52b ) with $\mathscr{c}=1$ , we can rewrite them by using $h(\tau )=l\Lambda (\tau )$ from (2.1), yielding

(3.18a) \begin{align}& \Lambda ^4\frac {1+4\lambda /\Lambda }{1+\lambda /\Lambda }\left (1+\frac {\lambda _{t}}{\Lambda }\right ) =\frac {1}{Hl^4}, \end{align}

(3.18b) \begin{align}&\qquad \frac {{\textrm d}H}{{\textrm d}\tau }=-\frac {1}{l(\Lambda +\lambda _{t})} . \end{align}

This set of equations consists of a first-order ordinary differential equation and an algebraic equation, which can be numerically solved as detailed in Appendix G. For comparison, the analytical solution for a no-slip condition can be obtained by substituting $\lambda =\lambda _t=0$ into the above equations, yielding

(3.19a) \begin{align}& H_{no-slip} = \left (1-\frac {3}{4}\tau \right )^{\frac {4}{3}}, \end{align}

(3.19b) \begin{align}& h_{no-slip} = \left (1-\frac {3}{4}\tau \right )^{-\frac {1}{3}}. \end{align}

We estimated that the magnitude of $l$ ranges from $10^{-1}$ to $10^{3}$ , based on the representative values of the following dimensionless parameters: $\textrm{St} = 0.01-0.1$ and ${A}^2 = 0.1-10$ , while $H_0^{*}$ varies from 0.01 m to 1 m and $l^{*}$ varies from 10 $\unicode{x03BC}$ m to 1 mm, with representative PCMs chosen as water and tetradecanol (Moallemi et al. 1986; Feuillebois et al. 2009); by examining their thermophysical properties, liquid metals such as gallium or caesium are inadequate due to $\textrm{Re}\mathscr{h}^2 \gtrsim 1$ . It is worth noting that $l=l^* / h_0^*=l^* [ (T_w^*-T_m^* ) k_l^* L^{* 2} \mu ^* ]^{-1/4}(\mathcal{H}^* \rho _l^* p_c^*)^{1/4}$ includes the information of the heating condition, thermophysical properties of materials and sizes of grooves and the PCM solid, which are convenient indicators. Hence, varying $\phi$ and $l$ allows us to investigate the influences of various materials and surface structures.

The representative results for $l \in \{10^{3}, 10^{2}, 10^{0}\}$ are drawn in figure 6 from figures 6(a) to 6(c), respectively. Since we have shown in § 3.2 that a flat interface ( $\parallel, f$ ) always reduces the melting rate and similar results occur for the gravity-driven CCM here, only results considering a meniscus ( $\parallel, c$ ) are discussed in the following. For large $l=10^3$ and $10^2$ , the results show that the melting rate ${\textrm d}H/{\textrm d}\tau$ decreases with increased $\phi$ (black line) as illustrated in figures 6(a) and 6(b), though the film thickness $h$ also decreases along with larger $\phi$ (blue line). This indicates that the temperature slip is more pronounced than the velocity slip. Compared with the case of no-slip (red dashed line), $\phi \leqslant 0.3$ exhibits an accelerated melting rate due to slip, similar to the pressure-driven CCM in figure 5(c). When $l=1$ in figure 6(c), it can be demonstrated that all slip cases gradually approach those of no-slip (figure 6 d and 6 e). Although no experimental data are currently available to validate the model predictions of $H(\tau )$ or $h(\tau )$ under varying gas–liquid fraction $\phi$ and structure parameter $l$ , experimental data for the limiting case of $\phi = 0$ were reported in two-dimensional rectangular CCMs by Moallemi et al. (1986), allowing for direct comparison. As shown in figure 6(c), the model prediction agrees well with the experimental results, confirming the validity of the limiting solution (i.e. the no-slip case) presented in this study.

In general, an increase in the value of $\phi$ will invariably result in a reduction in the thickness of the liquid film $h$ . However, the melting rate is decreased due to the more pronounced effect of temperature slip. The potential for slip to accelerate the melting process depends on the specific combination of $l$ and $\phi$ . When $l$ is sufficiently large (e.g. $10^2$ or $10^3$ ), even a modest increase in $\phi$ can enhance heat transfer, thereby accelerating the melting process. Conversely, as $l$ diminishes, the influence of slip decreases correspondingly, ultimately becoming negligible.

Figure 7. Phase diagram of melting time ratio $\tau _r$ of gravity-driven CCM with various combinations of gas–liquid fraction $\phi$ and periodic length $l$ on longitudinal grooves, where white dotted lines are contour lines for various $\tau _r$ : black solid line, $\phi = 1 - 2^{-2/3}$ ; black dashed line, the asymptotic limit of $l\lesssim 0.1$ for no-slip solutions (3.20); black dot–dash line, the asymptotic limits of (3.22a ) and (3.22b ) for solution (3.21).

3.3.2. Phase diagram and asymptotic solution

To illustrate the effects of $\phi$ and $l$ , we selected 288 and 110 values from the ranges $\phi$ from 0.1–0.9 and $l$ from $10^{-2}$ – $10^{3}$ , respectively. These values were used to construct the parameter matrix (log $_{10}l$ , $\phi$ ), which was then substituted into the numerical solution calculations to obtain the total melting time $\tau _{end}$ defined as $H(\tau _{end})=0$ for each case of (log $_{10}l$ , $\phi$ ). The resulting phase diagram is shown in figure 7, where the colour code represents $\tau _r = \tau _{end}/\tau _{no-slip}=3\tau _{end} /4$ , meaning the melting time ratio of slip surface to no-slip surface. $\tau _{no-slip}=4/3$ is obtained from (3.19a ), which corresponds to $\tau _r=1$ . Therefore, $\tau _r \lt 1$ means heat transfer enhancement, while $\tau _r \gt 1$ means heat transfer reduction.

It can be observed that the region where enhancement effects can be achieved exists only in the lower right-hand corner of the phase diagram. The upper limit boundary is at $\phi = 1 - 2^{-2/3}$ , which is consistent with the conclusions drawn from the constant pressure CCM in § 3.2.2. The left-hand critical point of the region of heat transfer enhancement is approximately located at $(1.17, 0.125)$ . Unlike the region of heat transfer reduction, where $\tau _r$ changes monotonically with $\phi$ when $l$ is fixed, in the enhancement region, $\tau _r$ first decreases and then increases as $\phi$ decreases. This behaviour indicates that for a given $l$ , there is always a specific value of $\phi$ that maximises heat transfer enhancement, which is plotted as yellow scatters in figure 7.

Based on this phase diagram, we can further find two regions corresponding to two asymptotes (see details in Appendix H). One approximation is related to the no-slip solution for $\Lambda \geqslant 10$ , as

(3.20) \begin{equation} H(\tau ) = \left (1-\frac {3}{4}\tau \right )^{\frac {4}{3}}, \quad l \lesssim 0.1. \end{equation}

Another approximate solution satisfies $\Lambda \lesssim 0.01$ , as

(3.21) \begin{equation} H(\tau ;\phi )= \left (1-\frac {3\sqrt {2}}{4}\tau (1-\phi )^{\frac {3}{4}} \right )^{\frac {4}{3}}, \end{equation}

whose approximate limit conditions should be satisfied as

(3.22a) \begin{align}&\qquad\quad 100\left (\frac {1-\phi }{4}\right )^{\frac {1}{5}} \lesssim l, \end{align}

(3.22b) \begin{align}& l\left [{\frac {4}{3}\sin \theta \frac {\phi ^2}{(1-\phi )(1-2\phi )}}\right ] \gtrsim 10. \end{align}

Therefore, the upper right-hand region enclosed by the two dash–dotted lines in figure 7 represents the parameter range where asymptotic solution (3.21) is valid.

Comparison of numerical and asymptotic results from $l = 1-10^3$ are conducted in Appendix H. It shows a good agreement with predictions in the phase diagram by selecting $\phi =\{0.1,0.2,0.3,0.4,0.5,0.6\}$ for cases.

4. Conclusions

In this work, we established a theoretical framework to investigate the dynamics of CCM on gas-trapped hydrophobic surfaces, with a particular focus on the effects of liquid film confinement and the meniscus at the gas–liquid interface. Utilizing dual-series and perturbation methods under the assumption of small meniscus deflections, we obtained numerical results for the velocity slip length $\lambda$ and temperature slip length $\lambda _t$ across various aspect ratios $\Lambda$ and gas fraction $\phi$ , representative of the liquid film thickness $h$ relative to the period $l$ of the microstructure. For $\Lambda \gg 1$ , the scaling laws are $\lambda _{\parallel, f} = \lambda _t = 2\lambda _{\perp, f} = \ln (\sec (\phi \pi /2) )/\pi$ , while for $\Lambda \ll 1$ , they are $\lambda _{\parallel, f} = \lambda _t = 4\lambda _{\perp, f} = \phi \Lambda /(1-\phi )$ , where $f$ , $\parallel$ and $\perp$ denote a flat interface, longitudinal grooves and transverse grooves, respectively. Accounting for the meniscus at the gas–liquid interface, the velocity slip length changes to (3.5) and (3.11), showing a slight decrease for $\Lambda \gg 1$ and a significant increase for $\Lambda \ll 1$ compared with the flat interface. This behaviour results from the larger area fraction and pronounced slip effects induced by the meniscus at $\Lambda \ll 1$ , whereas at $\Lambda \gg 1$ , the meniscus contributes negligible additional area, and the velocity gradient is reduced at the bottom. Additionally, a modified asymptotic solution for $\lambda _{\parallel, c}$ (3.12) is proposed, which is broadly applicable in the range $\Lambda \gtrsim 0.2$ .

For constant-pressure CCM, the dimensionless film thickness $h$ remains constant for a specific pressure $\mathcal{P}$ and both slip lengths, $\lambda$ and $\lambda _t$ , lead to a slip-effective $\textrm{Nu}$ that determines whether heat transfer is enhanced ( ${\textrm{Nu}} \gt 1$ ) or decreased ( ${\textrm{Nu}} \lt 1$ ). When the gas–liquid interface is flat, we found that transverse grooves surfaces always deteriorate heat transfer regardless of $\Lambda$ and $\phi$ , while longitudinal grooves exhibit behaviour identical to a smooth surface at small $\phi$ but show deteriorated performance at large $\phi$ . In the presence of a meniscus at the gas–liquid interface, we observed an enhanced melting rate when $\Lambda \lesssim 0.1$ and $\phi \lt 1-\sqrt [3]{2}/2 \approx 0.37$ . It is important to note that $\textrm{Nu}$ does not vary monotonically with $\phi$ ; instead, maximum values of $\textrm{Nu}$ occur under specific conditions within the range $\phi \lt 0.37$ seen in figure 5(c).

In gravity-driven CCM, the film thickness $h$ increases monotonically in time. Slip effects have a minor influence on $h$ and become negligible when $\Lambda \gtrsim 10$ or when there is a combination of small $\phi$ and $\Lambda \sim 10^{-1} $ – $10$ . Conversely, significant slip effects are observed for smaller $\Lambda$ , particularly for $\Lambda \lesssim 10^{-2}$ . We developed a two-dimensional phase diagram based on $(\log _{10} l, \phi )$ to identify regions of enhancement, reduction, or negligible impact of the gas–liquid boundaries (the generator of slip) relative to the no-slip case. Additionally, we derived asymptotic solutions and their limiting conditions.

The results reveal enhanced heat transfer and accelerated melting power in CCM are achievable only when a meniscus is present and there are significant confinement effects. The critical conditions depend on the gas–liquid interface fraction being less than 0.37. It should be noted that the quantitative evaluation of heat transfer enhancement or reduction presented in this study is valid primarily for small meniscus deflections. However, we believe that our conclusions regarding the critical conditions for heat transfer enhancement or reduction remain applicable even for large meniscus deformations. This is because, under the assumed boundary conditions, a highly curved interface at small $\Lambda$ increases the meniscus area ratio, which is expected to further amplify the observed trends. Nonetheless, a dedicated quantitative analysis is required in future work to confirm this effect.

From a practical perspective, this work has indicated that the conditions for utilizing gas-trapped hydrophobic surfaces to enhance CCM melting rates are quite stringent. Future work may explore the potential for heat transfer enhancement using liquid-infused slippery surfaces (Hardt & McHale 2022) or polymer brush surfaces (Chen et al. 2023) to significantly reduce the impact of effective thermal slip while maintaining velocity slip.