2.1. Problem description

We illustrate the theoretical framework using a minimal yet representative enclosure of PCM (Dhaidan & Khodadadi Reference Dhaidan and Khodadadi2015) – a rectangular cavity of width $W$ and height $H$ (see figure 1 a). It is filled with solid PCM initialised at the ambient temperature $\tilde T=T_c$ , which is below the fusion (melting) temperature $T_f$ . Keeping its top and bottom walls insulated, the cavity’s left wall absorbs heat from a hot working device at temperature $\tilde T=T_h\gt T_f$ , with its right boundary subject to a cold ambient environment at $\tilde T=T_c\lt T_f$ . This configuration finds typical application in enhancing thermal comfort in buildings (Lachheb et al. Reference Lachheb, Younsi, Youssef and Bouadila2024), where PCM is encapsulated within bricks, with $T_f$ lying between the indoor ( $T_c$ ) and outdoor ( $T_h$ ) temperatures. The ideal configuration – where $T_f$ always equals $T_c$ , preventing any heat transfer in the solid phase – has been discussed in previous studies (Esfahani et al. Reference Esfahani, Hirata, Berti and Calzavarini2018; Favier, Purseed & Duchemin Reference Favier, Purseed and Duchemin2019; Li et al. Reference Li, Jiao and Jia2022; Yang et al. Reference Yang, Chong, Liu, Verzicco and Lohse2022).

Figure 1.

(a) A rectangular cavity filled with subcooled PCM, maintaining temperatures at $T_h$ (hot) and $T_c$ (cold) at the left and right boundaries, respectively. (b) Time evolution of the domain-averaged velocity $\langle |{\boldsymbol{ u}}| \rangle =(g\alpha \Delta T_l H)^{-1/2}\int |\boldsymbol{\tilde u}| \mathrm{d} {\mathcal{V} } /{\mathcal{V}}_0$ , where ${\mathcal{V}}_0=\textit{WH}$ . Here, ${{Ra}}=10^7$ , ${Pr}=0.1$ , ${{St}}=0.1$ , ${{S}}_T=1.0$ and $\gamma =1.0$ . The inset shows the instantaneous temperature field and streamlines at ${\textit {Fo}}{{St}}=0.0046$ .

Using tildes to indicate dimensional variables, the thermo-hydrodynamic equations for the liquid phase read

(2.1a) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot} \boldsymbol{\tilde u} & =0,\\[-12pt]\nonumber \end{align}

(2.1b) \begin{align} \frac {\partial {\boldsymbol{ \tilde {u}} }}{\partial \tilde t}+\boldsymbol{{ \tilde {u}} }\boldsymbol{\cdot} \boldsymbol{\nabla }\boldsymbol{ \tilde {u}} & =\boldsymbol{g}\alpha ({ \tilde {T}} -T_0)+ {\nu \nabla ^2 \boldsymbol{{ \tilde {u}} }}-\frac {\boldsymbol{\nabla }{{ \tilde {p}} }}{\rho _0},\\[-12pt]\nonumber \end{align}

(2.1c) \begin{align} \frac {\partial {\tilde T}}{\partial \tilde t}+\boldsymbol{\tilde u}\boldsymbol{\cdot} \boldsymbol{\nabla } \tilde {T} & = {\kappa \nabla ^2 \tilde {T}}, \end{align}

where $\boldsymbol{\tilde u}$ , $\tilde T$ and $\tilde p$ denote the velocity, temperature, and pressure of the liquid, respectively; its kinematic viscosity is $\nu$ , while the thermal expansion coefficient and diffusivity are $\alpha$ and $\kappa,$ respectively. Here, $\rho _0$ denotes the reference density measured at the reference temperature $T_0=(T_h+T_f)/2$ and $\boldsymbol{g}$ indicates the gravitational acceleration. At the melting interface between the solid ( $s$ ) and liquid ( $l$ ) phase, the Stefan condition $\mathcal{L} \boldsymbol{\tilde V} \boldsymbol{\cdot} \boldsymbol{n}/C_p=[(\kappa \boldsymbol{\nabla } \tilde T)_s-(\kappa \boldsymbol{\nabla } \tilde T)_l]\boldsymbol{\cdot} \boldsymbol{n}$ holds, where $\mathcal{L}$ and $C_p$ are the specific latent heat and specific heat, respectively; $\boldsymbol{\tilde V}$ represents the velocity of the interface, and $\boldsymbol{n}$ denotes its liquid-facing unit normal.

The above governing equations denote that the material properties affecting the system performance include $\alpha$ , $T_f$ , $\nu$ , $\kappa$ , $C_p$ and $\mathcal{L}$ . Their effects, along with those of the geometry ( $W$ and $H$ ), and the working conditions ( $T_h$ and $T_c$ ) are characterised by five dimensionless numbers: the Rayleigh number ${Ra}$ , Prandtl number $Pr$ , Stefan number ${St}$ , subcooling strength ${{S}}_T$ and the cavity’s aspect ratio $\gamma$ , defined as

(2.2) \begin{equation} {{Ra}}=\frac {g\alpha \Delta T_l H^3}{\nu \kappa },\ {Pr}=\frac {\nu }{\kappa },\ {St}=\frac { C_p \Delta T_l}{\mathcal{L}}, {{S}}_T=\frac {\Delta T_s}{\Delta T_l},\ \gamma =\frac {W}{H}. \end{equation}

Here, $\Delta T_s=T_f-T_c$ and $\Delta T_l=T_h-T_f$ are the temperature differences across the solid and liquid phases, respectively. Additionally, the melting behaviour of PCMs is characterised by the Nusselt number $\textit {Nu}$ , Fourier time $\textit {Fo}$ and liquid fraction $f$

(2.3) \begin{equation} {\textit {Nu}}=-\frac {H}{\Delta T_l} \left \langle \frac {\partial \tilde T}{\partial \tilde y} \right \rangle _{\tilde y=0},\ {\textit {Fo}}=\frac {\kappa \tilde t}{W^2},\ f=\frac {\tilde w}{W}, \end{equation}

where $\tilde w$ is the average thickness of the liquid layer (figure 1 a), and $\langle \ \rangle _{\tilde y=0}$ denotes spatial averaging along the hot wall at $\tilde y = 0$ .

2.2. Linear stability analysis

We obtain the transition criterion through linear stability analysis, choosing the endpoint of the conduction-dominated regime as the base state. Since the horizontal temperature gradient imposes a heterogeneous buoyancy on the fluid (figure 1 a), the corresponding base flow $\boldsymbol{\tilde {u}}_b$ is not stationary, as evidenced by the temporal evolution of the domain-averaged velocity $\langle |\boldsymbol{u}| \rangle$ in figure 1(b). This non-stationary base flow $\boldsymbol{\tilde {u}}_b$ differs from the quiescent counterpart of the canonical Rayleigh–Bénard convection, and cannot be derived analytically. However, the base liquid temperature $\tilde {T}_b$ is readily calculable, because the convective heat transfer is negligible in the conduction-dominated regime. By solving the heat conduction equation, we obtain $\tilde {T}_b = T_h-{\Delta T_l} \tilde y /{\tilde w_b}$ , where $\tilde w_b = W/(1 + {{S}}_T^{\ast })$ results from the energy balance between the liquid and solid phases, $\Delta T_l/\tilde w_b= \Delta T_s/ ( W-\tilde w_b )$ . Here, the superscript $^{\ast }$ indicates that ${{S}}_T^\ast$ is the critical subcooling strength required for the base state to reach the onset of instability.

To linearise (2.1), we decompose the flow and temperature fields into the base state $ ( \tilde {\boldsymbol{u}}_b, \tilde {p}_b, \tilde T_b )$ and the disturbance state $ ( \tilde {\boldsymbol{u}}^{\prime }, \tilde {p}^{\prime }, \tilde T^{\prime } )$ , such that $ ( \tilde {\boldsymbol{u}}, \tilde {p}, \tilde T ) = ( \tilde {\boldsymbol{u}}_b, \tilde {p}_b, \tilde T_b )+ ( \tilde {\boldsymbol{u}}^{\prime }, \tilde {p}^{\prime }, \tilde T^{\prime } )$ . Selecting $\tilde w_b$ as the length scale, $U_0$ and $U$ the base and perturbation velocity scales, respectively, we derive the dimensionless linearised equations for $ ( \tilde {\boldsymbol{u}}^{\prime }, \tilde {p}^{\prime }, \tilde T^{\prime } )$

(2.4a) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot} \boldsymbol u' &=0 ,\\[-10pt]\nonumber \end{align}

(2.4b) \begin{align} \frac {\partial {\boldsymbol{ u}' }}{\partial t} + \frac {U_0 \tilde w_b}{\kappa } (\boldsymbol{ u }_b \boldsymbol{\cdot} \boldsymbol{\nabla }{\boldsymbol{u}'}+\boldsymbol{ u }' \boldsymbol{\cdot} \boldsymbol{\nabla }{ \boldsymbol{u}_b} ) &= \frac { {{Ra}}{Pr}\gamma ^3}{\big(1+{{S}}_T^{\ast }\big)^3} T' \boldsymbol{e_x}+ {Pr} {\nabla ^2} \boldsymbol{ u }' -\boldsymbol{\nabla }{{ p'} },\\[-10pt]\nonumber \end{align}

(2.4c) \begin{align} \frac {\partial {T'}}{\partial t} - v' + \frac {U_0 \tilde w_b}{\kappa }\boldsymbol{u}_b\boldsymbol{\cdot} \boldsymbol{\nabla } T' &= {\nabla ^2} T', \end{align}

where we have also chosen $\tilde w_b^2/\kappa$ as the characteristic time, $\rho _0 \kappa U/\tilde w_b$ as the characteristic pressure and $\Delta T_l U \tilde w_b /\kappa$ as the characteristic temperature. Here, $v^{\prime }$ denotes the $y$ -component of the dimensionless perturbation velocity $\boldsymbol{u}^{\prime }$ .

We note that $U_0\tilde w_b/\kappa$ in (2.4b ) and (2.4c ) represents the Péclet number, which measures the relative strength of convection to conduction. Since this number is inherently small in the conduction-dominated base state, we simplify the theoretical analysis by neglecting all terms involving $U_0\tilde w_b/\kappa$ in (2.4). This assumption proves reasonable, as will be verified a posterior by numerical simulations (see Appendix B). With this simplification, we obtain

(2.5a) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot} \boldsymbol u' &=0 ,\\[-10pt]\nonumber \end{align}

(2.5b) \begin{align} \frac {\partial {\boldsymbol{ u}' }}{\partial t} &= \frac { {{Ra}}{Pr}\gamma ^3}{\big(1+{{S}}_T^{\ast }\big)^3} T' \boldsymbol{e_x}+ {Pr}\nabla ^2 \boldsymbol{ u }' -\boldsymbol{\nabla }{{ p'} },\\[-10pt]\nonumber \end{align}

(2.5c) \begin{align} \frac {\partial {T'}}{\partial t} - v' &= \nabla ^2 T'. \end{align}

Using the normal mode approach, a dimensionless perturbation field $\phi$ of either the velocity, temperature or pressure is expressed as $\phi ' (x,y,t)=\hat \phi (y)\exp (\mathrm{i}\mathcal{K}x+\sigma t)$ , where $\hat \phi$ presents the complex magnitude of perturbation, $\mathcal{K}$ and $\sigma$ stand for the wavenumber and growth rate, respectively; here, $\mathrm{i}=\sqrt {-1}$ . Substituting the normal modes into (2.5), we obtain the governing equations for the perturbation magnitudes

(2.6a) \begin{align} 0&=\mathrm{i}\mathcal{K}\hat u +\frac {\partial \hat v}{\partial y},\\[-10pt]\nonumber \end{align}

(2.6b) \begin{align} \sigma \hat u &=\frac {{{Ra}}{Pr}\gamma ^3}{\big(1+{{{S}}_T^{\ast }}\big)^3} \hat T +{Pr}\left (\frac {\partial ^2 \hat u}{\partial y^2} -\mathcal{K}^2 \hat u\right )-\mathrm{i}\mathcal{K}\hat p,\\[-10pt]\nonumber \end{align}

(2.6c) \begin{align} \sigma \hat v &= {Pr} \left (\frac {\partial ^2 \hat v}{\partial y^2}-\mathcal{K}^2\hat v \right )-\frac {\partial \hat p}{\partial y},\\[-10pt]\nonumber \end{align}

(2.6d) \begin{align} \sigma \hat T &= \left ( \frac {\partial ^2 \hat T}{\partial y^2}-\mathcal{K}^2\hat T \right )+\hat v. \end{align}

The corresponding boundary conditions at the hot wall $y=0$ are

(2.7) \begin{align} \hat u = \hat v =0,\ \hat T =0. \end{align}

Those at the melting interface $y=1$ read

(2.8) \begin{align} \hat u = \hat v =0,\ \mathcal{K} \hat T \mathrm{coth}\big(\mathcal{K} {{{S}}_T^{\ast }}\big) + \frac {\partial \hat T}{\partial y}=0, \end{align}

where the condition for $\hat {T}$ is derived from the Stefan condition (Davis, Müller & Dietsche Reference Davis, Müller and Dietsche1984; Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2019; Lu et al. Reference Lu, Zhang, Luo, Wu and Yi2022).

Equations (2.6), (2.7) and (2.8) indicate that the transition depends on $Pr$ , ${Ra}$ , $\gamma$ and ${{S}}_T^{\ast }$ . To investigate the influence of these parameters, we simplify (2.5) into a single equation representing the marginal state ( $\sigma =0$ ). Taking the divergence of (2.5b ) and using (2.5a ), we derive

(2.9) \begin{equation} 0 =\frac {{{Ra}}{Pr}\gamma ^3}{\big(1+{{S}}_T^{\ast }\big)^3} \frac {\partial T'}{\partial x}-{\nabla ^2}{{ p'} }. \end{equation}

Next, differentiating (2.9) with respect to $y$ results in

(2.10) \begin{equation} 0 =\frac {{{Ra}}{Pr}\gamma ^3}{\big(1+{{S}}_T^{\ast }\big)^3} \frac {\partial ^2 T'}{\partial x \partial y}-{\nabla ^2}{\frac {\partial p'}{\partial y} }. \end{equation}

We then take the Laplacian of the $y$ -component of (2.5b ) and obtain

(2.11) \begin{equation} \frac {\partial }{\partial t} {\nabla ^2} {v' }={Pr}{\nabla ^4} v'-{\nabla ^2}\frac {\partial p'}{\partial y}. \end{equation}

Further subtracting (2.10) from (2.11) yields

(2.12) \begin{equation} \frac {\partial }{\partial t} {\nabla ^2} {v' }={Pr}{\nabla ^4} v'-\frac {{{Ra}}{Pr}\gamma ^3}{\big(1+{{S}}_T^{\ast }\big)^3} \frac {\partial ^2 T'}{\partial x \partial y}. \end{equation}

Incorporating the normal modes into (2.5c ) and (2.12), the perturbation magnitudes at the marginal state ( $\sigma =0$ ) are governed by

(2.13a) \begin{align} \mathcal{K}^2\hat T-\hat v-\hat T^{(2)} & = 0, \end{align}

(2.13b) \begin{align} -\mathrm{i}\mathcal{K}\gamma ^3{{Ra}}\hat T^{(1)} +\mathcal{K}^2\big(1+{{S}}_T^\ast \big)^3\big(\mathcal{K}^2\hat v-2\hat v^{(2)}\big)+\big(1+{{S}}_T^\ast \big)^3\hat v^{(4)} & =0, \end{align}

where $\hat T^{(n)} = {\partial ^n \hat T}/{\partial y^n}$ and $\hat v^{(n)} = {\partial ^n \hat v}/{\partial y^n}$ . By eliminating $\hat v$ , we derive

(2.14) \begin{equation} \mathcal{K}^6\hat T - \frac {\mathrm{i}{{Ra}}\gamma ^3 \mathcal{K}\hat T^{(1)}}{\big(1+{{{S}}_T^\ast }\big)^3} =3\mathcal{K}^4\hat T^{(2)}-3\mathcal{K}^2\hat T^{(4)}+\hat T^{(6)}, \end{equation}

which suggests that the transition is characterised solely by a critical effective Rayleigh number ${{Ra}}_e^{\ast }={{Ra}} \gamma ^3/(1+{{{S}}_T^\ast })^3$ .

2.3. Results and discussions

Realising the critical liquid fraction $f^{\ast }=1/ ( 1+{{S}}_T^{\ast } )$ upon substituting $\tilde w_b=W/(1+{{S}}_T^{\ast })$ into (2.3), we further obtain

(2.15) \begin{equation} f^{\ast }=\left (\frac {{{Ra}}_e^{\ast }}{{{Ra}} \gamma ^3}\right )^{1/3}. \end{equation}

Applying a spectral method (Weideman & Reddy Reference Weideman and Reddy2000) to (2.6)–(2.8), we calculate $f^{\ast }=1/ ( 1+{{S}}_T^{\ast } )$ for various configurations within the ranges of $10^6 \leqslant {Ra} \leqslant 10^{10}$ , $10^{-5} \leqslant Pr \leqslant 10^3$ and $0.5 \leqslant \gamma \leqslant 2$ , further yielding ${{Ra}}_e^{\ast }$ via (2.15). The results consistently exhibit a similar trend across these ranges. This trend is revealed by a representative set of data presented in figure 2.

Figure 2(a) shows that ${{Ra}}_e^{\ast }$ remains constant with respect to ${Ra}$ . Additionally, ${{Ra}}_e^{\ast }$ is independent of $Pr$ when ${Pr}\lt 10^{-3}$ and ${Pr}\gt 10$ , plateauing at $6248$ and $5200$ , respectively (figure 2 b). Notably, $f^{\ast }$ values corresponding to both the highest and lowest values of ${{Ra}}_e^{\ast }$ are closely matched, as indicated by $ ( 6248/5200 )^{1/3}\approx 1.06$ . Consequently, in the subsequent theoretical modelling, $f^{\ast }$ can be approximated using the median value of ${{Ra}}_e^{\ast }=5724$ . However, such an approximation fails for narrow (small $\gamma$ ) cavities, where ${{Ra}}_e^{\ast }$ increases rapidly when $\gamma$ decreases from $0.75$ , as evidenced in figure 2(c). This failure stems from the stronger confinement effect, as commonly observed in low-aspect-ratio canonical Rayleigh–Bénard convection (Wang et al. Reference Wang, Ma, Chen and Sun2012; Huang et al. Reference Huang2013b ; Huang & Xia Reference Huang and Xia2016; Yu et al. Reference Yu, Goldfaden, Flagstad and Scheel2017; Shishkina Reference Shishkina2021; Ahlers et al. Reference Ahlers2022).

Figure 2.

Typical dependence of the critical liquid fraction $f^{\ast }$ and critical effective Rayleigh number ${{Ra}}_e^{\ast }$ on (a) $Ra$ with ${Pr}=10$ and $\gamma =1$ , (b) $Pr$ with ${Ra}=10^6$ and $\gamma =1$ and (c) $\gamma$ with ${Ra}=10^6$ and ${Pr}=10$ .

Having derived the transition criterion $f^{\ast }$ , we employ it to develop the mathematical model for the melting curve. Following Esfahani et al. (Reference Esfahani, Hirata, Berti and Calzavarini2018), Favier et al. (Reference Favier, Purseed and Duchemin2019) and Yang et al. (Reference Yang, Chong, Liu, Verzicco and Lohse2022), we start from the dimensional energy equation for both phases (Voller, Cross & Markatos Reference Voller, Cross and Markatos1987; Huang et al. Reference Huang, Wu and Cheng2013a )

(2.16) \begin{equation} \frac {\partial {\tilde T}}{\partial \tilde t}+\boldsymbol{\tilde u} \boldsymbol{\cdot} \boldsymbol{\nabla } \tilde {T} ={\kappa \nabla ^2 \tilde {T}} -\frac {\mathcal{L}}{C_p}\frac {\partial {f_{loc}}}{\partial \tilde t}, \end{equation}

where $f_{loc}$ represents the local liquid fraction (see (A6)) and its volume average is equivalent to the macroscopic counterpart $f$ in (2.3). The first and last terms represent the change rates of the sensible and latent heats, respectively. The ratio of the two terms, i.e. the Stefan number, is typically far below unity; namely, ${St}=C_p\Delta T_l/\mathcal L\ll 1$ , because of the high $\mathcal L$ of PCMs. This condition as commonly satisfied in prior studies (Beckermann & Viskanta Reference Beckermann and Viskanta1989; Regin, Solanki & Saini Reference Regin, Solanki and Saini2009; Arasu & Mujumdar Reference Arasu and Mujumdar2012; Ho & Gao Reference Ho and Gao2013; Kean, Sidik & Kaur Reference Kean, Sidik and Kaur2019; Khan & Khan Reference Khan and Khan2019), allowing for neglect of the first term.

Further applying Gauss’s theorem, we obtain the integral form of (2.16)

(2.17) \begin{equation} \frac {\mathcal{L}}{C_p}\frac {\partial f}{\partial \tilde t} W =\kappa \left [\left\langle \frac {\partial \tilde T}{\partial \tilde y}\right\rangle _{\tilde y=W} - \left\langle \frac {\partial \tilde T}{\partial \tilde y}\right\rangle _{\tilde y=0} \right ], \end{equation}

where $\langle \partial \tilde T/ \partial \tilde y \rangle _{\tilde y=W}$ , the average temperature gradient at the cold wall, varies over time generally due to the advancing melting interface. However, in the limit of ${{St}}\ll 1$ as revealed above, the interfacial velocity is negligible compared with the conduction rate, as implied by the normalised Stefan condition $\boldsymbol{\tilde V} \boldsymbol{\cdot} \boldsymbol{n}/({\kappa }/H)={{St}} (\boldsymbol{\nabla } T_{s}+ \boldsymbol{\nabla } T_{l})\boldsymbol{\cdot} \boldsymbol{n}$ . Consequently, the conduction process can be considered quasi-steady on the time scale of the interfacial evolution, leading the temperature gradient at the cold wall to be

(2.18) \begin{equation} \left\langle \frac {\partial \tilde T}{\partial \tilde y}\right\rangle _{\tilde y=W}=-\frac { \Delta T_s}{W-\tilde {w}}, \end{equation}

where $W{-}\tilde w$ quantifies the average thickness of the solid phase (see figure 1 a).

At the hot wall, the temperature gradient $\langle \partial \tilde T/ \partial \tilde y \rangle _{\tilde y=0}$ depends on the conduction and convection states of the liquid phase. When conduction dominates, this gradient becomes $-\langle \partial \tilde T/\partial \tilde y\rangle _{\tilde y=0}=\Delta T_l/\tilde w$ , which combining (2.17) and (2.18) yields the conduction melting curve in the implicit form

(2.19) \begin{equation} \frac {f(1+{{S}}_T)(f-2{{S}}_T+f{{S}}_T)-2{{S}}_T\ln [1-f(1+{{S}}_T)]}{2(1+{{S}}_T)^3} = \textit {Fo}{St}. \end{equation}

Due to the subcooling effects, this solution differs from the classical Stefan solution $f \sim \sqrt {\textit {Fo}}$ (Esfahani et al. Reference Esfahani, Hirata, Berti and Calzavarini2018; Fu et al. Reference Fu, Yan, Gurumukhi, Garimella, King and Miljkovic2022), but recovers to it when ${{S}}_T=0$ .

Figure 3.

The present model is compared with (a–c) experimental measurements (squares) of n-octadecane ( $Pr=56$ ) in a cubic cavity (Ho & Gao Reference Ho and Gao2013), and (d–f) our two-dimensional numerical data (circles) for a rectangular cavity. The numerical data in (d–f) are represented on a log–log scale in (g–i) to highlight the three melting regimes.

As convection dominates, the temperature gradient can be modelled as $-\langle \partial \tilde T/\partial \tilde y \rangle _{\tilde y=0}={\textit {Nu}} \Delta T_l/H$ following (2.3), resulting in the convection melting curve

(2.20) \begin{equation} \frac {f}{{\textit {Nu}}\gamma }-\frac {{{S}}_T}{({\textit {Nu}}\gamma )^2 }\ln \left [ {\textit {Nu}}\gamma (f-1)+{{S}}_T \right ]= \textit {Fo} {St} +C, \end{equation}

where $C$ is a constant to be determined. Because both melting curves concatenate at $ ( {\textit {Fo}}^{\ast } {St},f^{\ast } )$ , subtracting (2.19) from (2.20) yields the constant

(2.21) \begin{align} C &=\frac {f^{\ast }}{{\textit {Nu}}\gamma }-\frac {{{S}}_T}{({\textit {Nu}}\gamma )^2 }\ln \left [ {\textit {Nu}}\gamma (f^{\ast }-1)+{{S}}_T \right ] \nonumber \\[4pt] & \quad -\frac {f^{\ast }(1+{{S}}_T)(f^{\ast }-2{{S}}_T+f^{\ast }{{S}}_T)-2{{S}}_T\ln [1-f^{\ast }(1+{{S}}_T)]}{2(1+{{S}}_T)^3}, \end{align}

where the Nusselt number is expressed using established models, specifically $\textit {Nu} \approx 0.5{{Ra}}^{1/4}{Pr}^{1/4}$ for ${Pr}\leqslant 0.1$ and $\textit {Nu} \approx 0.24{{Ra}}^{0.26}{Pr}^{0}$ for $ Pr\gt 0.1$ (Beckermann & Viskanta Reference Beckermann and Viskanta1989; Shishkina Reference Shishkina2016; Wang et al. Reference Wang, Liu, Verzicco, Shishkina and Lohse2021).

By substituting (2.21) into (2.20), the domain of the logarithmic function necessitates

(2.22) \begin{equation} f \lt f_s=1-\frac {{{S}}_T}{{\textit {Nu}} \gamma }. \end{equation}

This indicates that melting will cease at a saturated liquid fraction $f_s$ , leaving a portion of solid PCM unmelted.

Thus far, we have derived the critical liquid fraction $f^{\ast }$ and established a mathematical model for the melting curve. As shown in figure 3(a–c), the theoretical melting curves agree well with the prior experimental data (Ho & Gao Reference Ho and Gao2013). Besides, these curves also capture the trend revealed in our two-dimensional simulations, see figure 3(d–f). In the long-time limit, these curves asymptotically flatten, converging to $f_s$ as expected. To better characterise the transition, we depict these curves on a log–log scale in figure 3(g–i). The temporal evolution of $f$ exhibits three distinct phases: the first two follow a power-law relationship $f \sim (\textit {Fo} {St})^\beta$ with exponents $\beta \approx 0.5$ and $\beta \approx 1.0$ , respectively, followed by a saturated phase. Notably, this power-law scaling is exact in no-subcooling cases ( ${{S}}_T=0$ ), where $\beta =0.5$ and $1.0$ correspond to the conduction-dominated and convection-dominated regimes, respectively (Esfahani et al. Reference Esfahani, Hirata, Berti and Calzavarini2018; Favier et al. Reference Favier, Purseed and Duchemin2019; Li et al. Reference Li, Jiao and Jia2022). For present subcooling scenarios ( ${{S}}_T\gt 0$ ), the first two regimes with different $\beta$ are precisely demarcated by the predicted transition criterion $f=f^{\ast }$ (grey dashed line). Our theory’s predictive ability is further validated through benchmark comparisons with both our simulations and published data (Beckermann & Viskanta Reference Beckermann and Viskanta1989) as detailed in Appendix A.

Despite the generally satisfactory agreement between theoretical predictions and the experimental/numerical data, slight discrepancies are observed. They arise from two approximations: (i) neglecting the time variation of sensible heat in (2.16) and (ii) assuming a steady temperature gradient at the cold wall to obtain (2.18). Both approximations hold when ${St} \ll 1$ . At ${St} = 0.1$ , figure 3 indicates that the maximum relative error between the theoretical and numerical predictions of $f$ is ${\approx}15\,\%$ . This error increases with $St$ , reaching ${\approx}20\,\%$ when ${St} \in ( 0.2, 0.4 )$ , see Appendix C.

We now exploit the theoretical melting curve $f(\textit {Fo})$ to predict the key performance indicators. We start with the energy density $\textit {ED}$ as the cumulative energy stored per unit volume of PCM (Woods et al. Reference Woods, Mahvi, Goyal, Kozubal, Odukomaiya and Jackson2021; Fu et al. Reference Fu, Yan, Gurumukhi, Garimella, King and Miljkovic2022), which depends on the liquid fraction $f$ through

(2.23) \begin{equation} \textit {ED}(\tilde {t})=\rho [C_p(T_f -T_c)+\mathcal{L}] f . \end{equation}

Here, we have neglected the liquid’s sensible heat that plays a negligible role, as demonstrated in Fu et al. (Reference Fu, Yan, Gurumukhi, Garimella, King and Miljkovic2022) and confirmed by our experimentally verified calculations (see Appendix E). Notably, the maximum value of $\textit {ED}$ , reached at $f=f_s$ , represents the energy storage capacity, whereas the temporal derivative of $\textit {ED}$ serves as another indicator, the power density (Woods et al. Reference Woods, Mahvi, Goyal, Kozubal, Odukomaiya and Jackson2021; Fu et al. Reference Fu, Yan, Gurumukhi, Garimella, King and Miljkovic2022).

Practically, (2.23) characterises how design parameters of a PCM system affect its performance via the dimensionless numbers encoded in the melting curve. This characterisation, requiring only simple algebraic calculations, enables the rapid determination of parameters customised to meet specified performance criteria, thereby streamlining the system’s design.

Besides precisely characterising performance, the theoretical solutions also provide general guidance for optimising design parameters. For example, the critical liquid fraction $f^{\ast }$ signifies the duration of the conduction-dominated regime, which features inefficient heat transfer. Hence, reducing $f^{\ast }$ can enhance the melting rate (Sun et al. Reference Sun, Zhang, Medina and Lee2016; Azad et al. Reference Azad, Groulx and Donaldson2021a ), i.e. the power density. In the current configuration, this can be achieved by increasing the aspect ratio $\gamma$ , see (2.15).