Melting and solidification in periodically modulated thermal convection

Abstract

Melting and solidification in periodically time-modulated thermal convection are relevant for numerous natural and engineering systems, for example, glacial melting under periodic sun radiation and latent thermal energy storage under periodically pulsating heating. It is highly relevant for the estimation of melt rate and melt efficiency management. However, even the dynamics of a solid–liquid interface shape subjected to a simple sinusoidal heating has not yet been investigated in detail. In this paper, we offer a better understanding of the modulation frequency dependence of the melting and solidification front. We numerically investigate periodic melting and solidification in turbulent convective flow with the solid above and the melted liquid below, and sinusoidal heating at the bottom plate with the mean temperature equal to the melting temperature. We investigate how the periodic heating can prevent the full solidification, and the resulting flow structures and the quasi-equilibrium interface height. We further study the dependence on the heating modulation frequency. As the frequency decreases, we found two distinct regimes, which are ‘partially solid’ and ‘fully solid’. In the fully solid regime, the liquid freezes completely, and the effect of the modulation is limited. In the partially solid regime, the solid partially melts, and a steady or unsteady solid–liquid interface forms depending on the frequency. The interface height can be derived based on the energy balance through the interface. In the partially solid regime, the interface height oscillates periodically, following the frequency of modulation. Here, we propose a perturbation approach that can predict the dependency of the oscillation amplitude on the modulation frequency.

1. Introduction

Periodic thermal modulations are common in our daily life (Berger & Wille 1972; Davis 1976), from low frequencies (in geophysical applications, e.g. daily and seasonal sunlight cycles) to high frequencies (in industrial applications, e.g. electrical pulses). Periodic thermal modulation with phase transitions also plays an important role in this broad range of problems, in both nature and industrial applications, such as sea ice melting with tidal warm current and melt ponds with sun radiation (Perovich & Polashenski 2012; Kim et al. 2018; Popović et al. 2018; Yang et al. 2023a), tidal heating in Enceladus (Meyer & Wisdom 2007), seasonal heating and cooling, and phase change materials (PCMs) with cycles of storing and releasing energy (Sharma et al. 2009). Understanding this frequency-dependent nonlinear thermal response to time-periodic boundary conditions is necessary for evaluating the heat transport in these systems and dynamics of the melting front, which are critical to e.g. estimate the glacier melt rate in the geophysical context and the optimization of PCM-based thermal management systems in industrial applications.

The effect of periodic thermal boundary conditions on single-phase flow has been studied in depth (Jin & Xia 2008; Yang et al. 2020; Urban et al. 2022). Contrary to the general thought that the time-averaged global quantities are unchanged by modulation because the net force averaged over a cycle vanishes, a significant enhancement (up to 25 %) of heat transfer was found with a moderate period of thermal oscillation (${\sim }100$ free-fall time units) from the bottom wall, due to the perturbation of the thermal boundary layers. Due to the significant effect of modulation, different ways of modulation are also studied, such as temporally modulated temperature at the boundary (Jin & Xia 2008; Yang et al. 2020; Urban et al. 2022), modulated rotation (Geurts & Kunnen 2014; Sterl, Li & Zhong 2016) and modulated gravity (Gresho & Sani 1970; Rogers et al. 2000), which all have important effects on the heat transport.

With a solid–liquid phase transition, the problem of periodic heating and cooling is substantially more complicated due to the presence of a freely moving phase boundary (Stefan problem). It becomes a two-way coupled problem, i.e. the interface shape changes due to the flow beneath while such change of shape has feedback to the flow structures. Previously, phase change problems with steady forcing have been studied both numerically (Favier, Purseed & Duchemin 2019; Purseed et al. 2020; Couston et al. 2021; Ravichandran & Wettlaufer 2021; Yang et al. 2023b,c) and experimentally (Davis, Müller & Dietsche 1984; Dietsche & Müller 1985), where the focus was on the interface height and the roughness evolution. Similar to the modulation effect on single-phase flow, thermal modulations also have a significant and relevant influence on phase evolution, but they have hardly been explored so far.

Some previous work investigated melting and solidification under modulated temperature boundary conditions, for situations with pure conduction (Shamberger et al. 2020) and with weak convection (Mazzeo et al. 2015). These previous studies focus mainly on the phase change without turbulent flow. However, the coupling dynamic between the solid–liquid interface and the turbulent convection has been shown to have a significant effect on the melting and solidification processes (Favier et al. 2019; Couston et al. 2021), as well effects on fluid properties, such as density anomalies (Wang et al. 2021c,d; Wang, Calzavarini & Sun 2021b; Yang et al. 2022) and salinity (Du et al. 2023; Yang et al. 2023b). A comprehensive review of phase change with flows can be found in a recent study (Du, Calzavarini & Sun 2024). Moreover, the evolution of the solid–liquid interface over different thermal modulation frequencies is even less explored.

In order to better understand the solid–liquid interface dynamics for modulated heating and cooling, and give a full picture of the parameter space, we select turbulent Rayleigh–Bénard convection (RBC) as a model system, where a fluid is heated from below and cooled from above. The RBC is a paradigmatic example in the study of global heat transport in thermally driven turbulent flow (Ahlers, Grossmann & Lohse 2009; Lohse & Xia 2010; Chillà & Schumacher 2012; Xia 2013; Shishkina 2021; Lohse & Shishkina 2023, 2024), as it shares characteristics common to many systems of interest for natural and industrial applications, and is also widely applied to investigate the dynamics of multiphase flow (Zhong, Funfschilling & Ahlers 2009; Lakkaraju et al. 2013; Wang, Mathai & Sun 2019; Liu et al. 2022a,b) and phase changes driven by convective heat transfer (Davis et al. 1984; Dietsche & Müller 1985; Favier et al. 2019; Purseed et al. 2020). We add a simple harmonic heating temperature boundary condition at the bottom plate and a constant cooling temperature boundary condition at the top plate. We model the melting and solidification process with the phase-field method, which is used widely for phase-change problems (Favier et al. 2019; Purseed et al. 2020; Couston et al. 2021; Ravichandran & Wettlaufer 2021). The objective of this paper is to show how a pre-existing fluid layer resists full solidification while periodically receiving or losing heat through the bottom boundary, and to quantify the solid–liquid interface height.

The paper is organized as follows. The set-up and numerical methods are described in § 2. The main results are presented in §§ 3–5. The flow and solid–liquid interface structure, and their temporal evolution under different modulation frequencies, are discussed in § 3. The dependence of the average solid–liquid interface height and the heat transfer on the modulation frequency is discussed in § 4. The oscillation amplitude of the solid–liquid interface is discussed in § 5. The paper ends with the conclusions and an outlook in § 6.

2. Governing equations and control parameters

The flow in RBC is confined between two parallel plates separated by a distance $H$, with gravitational acceleration $g$ acting vertically to these plates. We numerically solve the velocity field $\boldsymbol {u}$ and the temperature field $\theta$ in the liquid phase from the Navier–Stokes equations within the Oberbeck–Boussinesq approximation with the direct numerical simulations solver AFiD, which is a second-order staggered finite difference, open-source code from our research group (Verzicco & Orlandi 1996; van der Poel et al. 2015a). It has already been validated extensively, and applied to studies of turbulent flows (van der Poel et al. 2015b; Yang, Verzicco & Lohse 2016; Yang et al. 2020; Wang, Lohse & Shishkina 2021a). To simulate the phase transition process, we use AFiD and the phase-field method presented by Favier et al. (2019). In this method, the phase-field variable $\phi$ is continuous in time and space, and transitions smoothly from value $1$ in the solid to value $0$ in the liquid. The applied phase-field model was initially derived based on entropy conservation, which guarantees the thermodynamic consistency, and also satisfies the Gibbs–Thompson relation (Wang et al. 1993; Favier et al. 2019). The implementation and validation of the phase-field model are shown in previous work (Liu et al. 2021; Yang et al. 2022).

The complete governing equations of the flow include the continuity equation, the momentum equation and the heat transfer equation:

(2.1)\begin{gather} \frac{\partial u_i}{\partial x_i}=0, \end{gather}

(2.2)\begin{gather} \frac{\partial u_i}{\partial t}+u_j\,\frac{\partial u_i}{\partial x_j} ={-}\frac{\partial p}{\partial x_i} + \theta\delta_{iz}+\sqrt{\frac{Pr}{Ra}}\,\frac{\partial^2 u_i}{\partial x_j^2}-\frac{(1-\phi)^2u_i}{\eta}, \end{gather}

(2.3)\begin{gather} \frac{\partial \theta}{\partial t} ={-}u_i\,\frac{\partial \theta}{\partial x_i}+(1+(\lambda-1)(1-\phi))\sqrt{\frac{1}{Ra\,Pr}}\,\frac{\partial^2 \theta}{\partial x^2_j}-\frac{1}{St}\,\frac{{\rm d}Q(\phi)}{{\rm d}\phi}\,\frac{{\rm d}\phi}{{\rm d}t}, \end{gather}

where $\delta _{iz}$ is the Kronecker delta function.

The phase change process is modelled by the phase-field equation

(2.4)\begin{equation} \frac{\partial\phi}{\partial t} = M\,\nabla^2\phi +St\,\frac{\alpha M}{\epsilon}\,(\theta-\theta_m)\,\frac{{\rm d}p}{{\rm d}\phi}-\frac{M}{4\epsilon^2}\,\frac{{\rm d}G(\phi)}{{\rm d}\phi}. \end{equation}

The independent control parameters in these equations are the Rayleigh number $Ra$ (measuring the strength of the thermal driving), the Prandtl number $Pr$ (intrinsic material property of the liquid), the Stefan number $St$ (ratio of sensible heat and latent heat), and the dimensionless cooling temperature at the top plate $-\theta _c$:

(2.5a––d)\begin{equation} Ra=\frac{\beta g (T_h-T_m) H^3}{\nu\kappa_l} ,\quad Pr=\frac{\nu}{\kappa_l} ,\quad St=\frac{c_p(T_h-T_m)}{\mathcal{L}} ,\quad \theta_c=\frac{T_m-T_c}{T_h-T_m}, \end{equation}

where $\beta$ is the thermal expansion coefficient, $\nu$ is the kinematic viscosity of the liquid, $\kappa _l$ is the thermal diffusivity in the liquid phase, $c_p$ is the specific heat capacity, and $\mathcal {L}$ is the latent heat. Here, $T_h$, $T_m$ and $T_c$ are the magnitudes of the heating temperature at the bottom plate, the melting temperature (which may be variable due to e.g. pressure effects; Couston 2021), and the cooling temperature at the top plate, respectively.

In the governing equations, the length is rescaled by the domain height $H$, the temperature $\theta$ is a dimensionless representation of temperature $T$, relative to the melting point at atmospheric pressure and scaled by the temperature difference $T_h - T_m$ across the liquid phase, and the velocity is rescaled by the associated free-fall velocity $U_f = \sqrt {g\beta H (T_h-T_m) }$ and corresponding free-fall time $t_f=H/U_f$. Note that the free-fall time scale is much shorter than the diffusive time scale $t_d$, with the relation $t_f=t_d/\sqrt {Ra\,Pr}$, so the modulation period in our parameter space is always shorter than the diffusive time scale. The dimensionless melting temperature of the solid is set as $\theta _m=0$. For the periodically modulated thermal boundary condition, we take a sinusoidal modulation signal to the dimensionless bottom temperature as

(2.6)\begin{equation} \theta_{h} = \sin(2{\rm \pi} ft), \end{equation}

where we introduce the modulation frequency $f$ non-dimensionalized by the free-fall time unit, as an additional control parameter. Equation (2.6) implies that the actual temperature difference between the bottom and the solid–liquid interface varies sinusoidally between $-(T_h-T_m)$ and $(T_h-T_m)$. Note that the bottom boundary can go below the freezing point without solidification because solidification cannot occur directly at the bottom plate where there is only the liquid phase ($\phi =0$) in the method that we use. Physically, this can be regarded as the undercooling effect, where no nuclei formation and phase transition occurs, in spite of being below the freezing point. This also has relevance to the natural environment (e.g. the temperature of cold polar waters sometimes drops below the freezing point; Haumann et al. 2020) as well as industrial applications such as the PCMs, where the undercooling effect is also widely observed (Zahir et al. 2019; Shamseddine et al. 2022).

In (2.4), the phase-field method includes the dimensionless mobility $M$, the dimensionless measurement of the interface thickness $\epsilon$, the coupling parameter $\alpha$, and the penalty parameter $\eta$. The function $G(\phi )=\phi ^2(1-\phi )^2$ is a double-well potential function, and $Q(\phi )=\phi ^3(10-15\phi +6\phi ^2)$ is a smoothing function to ensure a smooth transition between the solid and liquid phases. The choices of these parameters are $M=1$ and $\alpha =St/\epsilon$, $\epsilon$ is proportional to the grid size, and $\eta$ is equal to the time step, which is the same as in the previous study (Favier et al. 2019).

Our main focus in this paper is the effect of the modulation frequency $f$ and the dimensionless undercooling temperature $\theta _c$ on the melting and solidification dynamics. We will vary $10^{-5}< f\leq 1$, following the range in a previous study of thermal modulated convection (Yang et al. 2020), and $0.1\leq \theta _c\leq 1.4$, which is a common temperature ratio for ice and PCM (Dietsche & Müller 1985). Frequency $f$ is typically low in geophysical contexts due to factors like ocean currents and seasonal sunlight, while the heating frequency of PCM can span across various frequency regimes (Sharma et al. 2009), making the range of $f$ explored in this study applicable to different scenarios. If $\theta _c=0$, then the solid will finally melt completely, and $\theta _c\to \infty$ means that the liquid will finally freeze completely. We fix $Ra=10^8$, $Pr=10$ and $St=10$, which represents water at approximately $8\,^\circ {\rm C}$. The thermal diffusivity ratio is $\lambda =\kappa _{s}/\kappa _{l}=7$ for ice (James 1968). When the modulation is absent, one would expect that for any $\theta _c>0$, the system will finally freeze completely.

We conduct two-dimensional simulations in a domain with aspect ratio $\varGamma =W/H=2$, where $W$ and $H$ are the horizontal and vertical lengths of the domain, respectively, as shown in figure 1. Note that in the RBC for $Pr\gtrsim 1$, there are close similarities between two- and three-dimensional RBCs (van der Poel, Stevens & Lohse 2013). We impose no-slip boundary conditions for the top and bottom plates, and periodic boundary conditions in horizontal directions. Initially, the velocity field is set to $\boldsymbol {u}=0$, the temperature field $\theta$ in the liquid is set as a linear profile with some random fluctuations to trigger the convective flow, and we have a flat solid–liquid interface at $H/2$, as shown in figure 1. We also conducted a resolution dependence check of our simulation at $Ra=10^8$, $f=0.01$, $\theta _c=0.1$, as shown in figure 2.

Figure 1. Illustration of the two-dimensional set-up and the boundary conditions. The black curve shows a qualitative temperature profile.

Figure 2. (a) The averaged melt front $\bar {h}$ evolution as a function of time for different resolutions $N$ in the vertical direction at $Ra=10^8$, $f=0.01$, $\theta _c=0.1$. (b) The averaged melt front height as a function of $N$. One can see the convergence of $\bar {h}$, and our final choice of resolution is $N=192$.

One of the key response parameters in the system is the interface height ${h}$. The dimensionless local interface height $h(x,t)/H$ is defined by the location $\phi =1/2$. In this paper, we focus mainly on the temporal evolution of the horizontally averaged interface height $\bar {h}(t)=(1/W)\int ^W_0h(x,t)\,{{\rm d}\kern 0.7pt x}$. Another key response is the dimensionless heat flux $Nu$, defined as

(2.7)\begin{equation} Nu=\frac{{Q}\bar{h}}{\theta_h-\theta_m}, \end{equation}

where ${Q}$ is the dimensionless vertical temperature gradient at the bottom plate.

3. Flow structures with temperature oscillation

In figure 3, we show the representative series of temperature field variations during one period of sinusoidal heating temperature oscillations for different frequencies and $\theta _c=0.1$. Two distinct flow regimes can be identified, as follows.

1. (i) The ‘fully solid’ regime for high frequencies (e.g. $f=1$ as shown in figure 3a i–iv). Here, the influence of the modulation is negligible because it is too fast to be sensed by the system. Given that the temperature at the top is lower than the melting temperature, the system freezes completely in the final state.

2. (ii) The ‘partially solid’ regime for low frequency (e.g. $f=10^{-2}$ as shown in figure 3b i–iv). In this regime, the heating time in one period is long enough to initiate convective flow, which breaks the symmetry between heating and cooling. There are two stages in one period: one is the heating phase (figure 3b ii) when $\theta _h>0$, and the other is the cooling phase (figure 3b iv) when $\theta _h<0$. In the heating phase, frequent plume emissions are observed near the bottom plate, while in the cooling phase, there is no plume emission from the bottom plate due to the stable stratification. Although there is still convective flow in the cooling phase, it is much weaker than that in the heating phase. Therefore the solid partially melts. The ‘partially solid’ regime can be further distinguished in two sub-regimes: one with convection always active for intermediate frequency (see figure 3(b); another with intermittent convection, i.e. switching on and off with the heating and cooling phases, for low frequency (see figure 3c).

Figure 3. (a) Snapshots of the temperature field in the liquid phase for $f=1$ (note that $t_f=\sqrt {Ra\,Pr}\,t_d$) as time evolves (i–iv). The system finally freezes completely. (b,c) Snapshots of the temperature field in the liquid phase corresponding to four different phases of the sinusoidal period, for (b) $f=10^{-2}$, (c) $f=10^{-4}$, once a statistical steady state is achieved. The colour map ranges from $\theta =-1$ to $\theta =1$. The solid phase is represented in white, and the black solid line represents the interface.

For the heating phase, when the convective heat flux is strong and the solid undergoes significant melting, the conductive heat flux within the solid also intensifies. This is described by the conduction relation $\theta _c/(1-\bar {h})$, where $\bar {h}$ increases as the solid melts. Consequently, a balance is achieved between the increased heat flux in both the liquid and the solid. For the cooling phase, when the heat flux is weak and the solid undergoes more freezing, the conductive heat flux within the solid weakens. This results in another balance of heat flux between the liquid and the solid. Therefore, in the final equilibrium stage, the total amounts of melting and freezing during the heating and cooling phases become equal, leading the solid–liquid interface to reach an equilibrium height.

As the frequency decreases further (e.g. $f=10^{-4}$ as shown in figure 3c i–iv), with very long heating and cooling times in one period, we found the solid–liquid interface to become more unsteady, i.e. the interface height increases in the heating phase (figure 3c ii) and then decreases in the cooling phase (figure 3c iv). During the heating phase, larger temperature differences result in stronger convective flow and thus a higher solid–liquid interface. Moreover, multiple plumes merge into a single strong plume near the bottom plate, which significantly deforms the solid–liquid interface. During the cooling phase, smaller temperature differences suppress the convective flow, and the height of the solid–liquid interface becomes lower. Thus flow is completely damped out, and only purely diffusive heat flux remains.

In figure 4, we show the temporal evolutions of the mean temperature profile for different frequencies, where we can see the layer movement in the ‘partially solid’ regime. In figures 4(a,b), the interface is independent of time. When the modulation frequency $f$ decreases (figures 4c,d), the oscillation magnitude of the interface gradually grows, and the oscillation period of the interface stays the same as that of the thermal modulation at the bottom plate. We can also see the asymmetric solid–liquid interface evolution: as the temperature of the bottom plate increases in the heating phase, the solid–liquid interface rises quickly due to strong convective heat flux, while as the temperature of the bottom plate decreases in the cooling phase, the solid–liquid interface moves downwards slowly due to the weak convective heat flux. In figure 5, we show the temporal evolutions of the mean temperature profile for different $\theta _c$, where we can see that $\theta _c$ affects mainly the equilibrium layer height.

Figure 4. Temporal evolutions of mean temperature profiles for (a) $f=10^{-2}$, (b) $f=10^{-3}$, (c) $f=4\times 10^{-4}$, (d) $f=10^{-4}$ at $\theta _c=0.1$. The horizontal axis represents time, with total length $6T_0$, where $T_0$ is the time of one period in each case. The colour map ranges from $\theta =-1$ to $\theta =1$. The solid phase is represented in white. One can see that the period of the solid–liquid interface matches the period of thermal modulation at the bottom plate, and the oscillation amplitude of the solid–liquid interface increases as $f$ decreases.

Figure 5. Temporal evolutions of mean temperature profile for (a) $f=10^{-2}$, $\theta _c=0.1$; (b) $f=10^{-2}$, $\theta _c=0.4$; (c) $f=4\times 10^{-4}$, $\theta _c=0.1$; (d) $f=4\times 10^{-4}$, $\theta _c=0.4$. The horizontal axis represents time, with total length $6T_0$, where $T_0$ is the time of one period in each case. The colour map ranges from $\theta =-1$ to $\theta =1$. The solid phase is represented in white.

In figure 6, we present the explored parameter space (in the parameter space spanned by $f$ and $\theta _c$) and classify two different regimes, which again we identify qualitatively. The ‘fully solid’ regime exists for high $f$ and high $\theta _c$ because the convective heat transfer is limited under fast modulation, and high undercooling temperature at the top plate also tends to freeze the liquid. The ‘partially solid’ regime exists for low $f$ and low $\theta _c$, where convective heat transfer is strong enough to break the symmetry between the heating phase and the cooling phase so that the net heat melts the solid instead of it being completely frozen. For very low $f$, when the period is long enough, the solid–liquid interface oscillates within one period. When $f$ is low enough, the solid can freeze completely in one cooling phase, which requires a much lower $f$ than that explored in our parameter range.

Figure 6. Explored parameter space in the $f$ versus $\theta _c$ parameter plane, displaying the different flow regimes. These are indicated by different colours: the ‘fully solid’ regime is black, and the ‘partially solid’ regime colours represent the mean height $\bar {h}$.

4. Quasi-equilibrium interface height

To distinguish the two regimes quantitatively, we plot the equilibrium solid–liquid interface height $\bar {h}$ as a function of $f$ at different $\theta _c$ in figure 7(a). The dependence of $\bar {h}$ on $f$ shows a similar trend for different $\theta _c$. When $f$ is large, $\bar {h}=0$ because the liquid freezes completely (‘fully solid’ regime). As $f$ decreases to a certain value, $\bar {h}$ starts to increase after reaching equilibrium (‘partially solid’ regime). With $f$ decreasing further, the mean heating $\bar {h}$ reaches a local maximum and then starts to decrease again. For frequencies slightly smaller than the optimal one, the solid–liquid interface starts to oscillate obviously within the period. The oscillation amplitude of the solid–liquid interface is shown as the shaded region, which is the region between the maximum and minimum $h$ in one period.

Figure 7. (a) Mean equilibrium solid–liquid interface height $\bar {h}$ as a function of the modulation frequency $f$ at different $\theta _c$. The shaded regions represent the region between the maximum height and minimum height of the solid–liquid interface during one period. Two regimes are classified: ‘fully solid’ and ‘partially solid’ as $f$ decreases. (b) Plots of $Nu$ at the bottom plate as a function of the modulation frequency $f$ at different $\theta _c$. Here, $Nu$ is calculated based on (2.7). In all cases, $Nu(f)$ shows a pronounced maximum at medium frequencies.

Based on the Stefan boundary condition, the solid–liquid interface is related directly to the heat flux through the interface. Note that in this section, we consider only the time-averaged interface height and heat flux; the temporal oscillation of the interface and the heat flux will be considered in § 5. In figure 7(b), we plot the time-averaged heat flux $Nu$ (averaged over the whole simulation time) at the bottom plate as a function of the frequency $f$ for different top plate temperatures $\theta _c$. For all these $\theta _c$, $Nu$ shows a trend similar to that of $\bar {h}$: $Nu$ keeps constant at high $f$; as $f$ decreases, $Nu$ increases and reaches an optimal point, after which $Nu$ decreases again. The trend of $Nu$ versus $f$ obtained here behaves similarly to that of single-phase modulated thermal convection in Yang et al. (2020) because of the same mechanism of perturbing the boundary layers by the temporal modulation. The effect of modulation can be explained by introducing the Stokes thermal boundary layer (Yang et al. 2020), which affects $Nu$ by disturbing the thermal boundary layer and velocity boundary layer. Depending on the thickness of these boundary layers, different regimes of modulation frequency can be classified; see again Yang et al. (2020).

Another explanation for the observed different regimes can be attributed to the separation of time scales. The equations governing the motion are normalized by the free-fall time scale $t_f$, where the characteristic velocity is determined by the buoyancy difference between the melt temperature and the bottom temperature, and the total depth of the domain. It can be inferred that when modulation is faster than the free-fall time scale, the frequency of boundary layer oscillations is on the same time scale as the free-fall time scale. When this occurs, the boundary condition is evolving too fast for the convective plumes to fully develop, thus the system enters the ‘solid’ regime, where the bottom boundary condition is effectively zero by a separation of time scales. When the modulation is slower than the diffusive time scale ($\sqrt {Ra\,Pr}\,t_f$), there can be an obvious movement of the front over a long enough time in one period. For modulation periods falling between the free-fall time scale and the diffusive time scale, modulation influences heat transfer and results in melting, yet there is insufficient time for the melt front to adequately respond within a single period.

The relation between heat flux and the equilibrium height can be derived based on the non-dimensionalized Stefan boundary condition:

(4.1)\begin{equation} St\,\frac{{\rm d}h}{{\rm d}t}=\frac{1}{\sqrt{Ra\,Pr}}\,(Q_l-Q_s)=\frac{1}{\sqrt{Ra\,Pr}}\left(Nu\,\frac{\theta_h-\theta_m}{h}-\lambda\,\frac{\theta_m-(-\theta_c)}{1-h}\right), \end{equation}

where $Q_l$ and $Q_s$ represent the heat flux through the liquid phase and the solid phase, respectively. We have set the melting temperature as $\theta _m=0$. At equilibrium state, i.e. ${{\rm d}h}/{{\rm d}t}=0$, (4.1) is reduced to

(4.2)\begin{equation} Nu\,\frac{1}{\bar{h}}=\lambda\,\frac{\theta_c}{1-\bar{h}}. \end{equation}

From (4.2), we can obtain the expression for the equilibrium height $\bar {h}$:

(4.3)\begin{equation} \bar{h}=\frac{{Nu}/{\theta_c}}{\lambda+{Nu}/{\theta_c}}. \end{equation}

We plot $\bar {h}$ versus $Nu/\theta _c$ in figure 8. The numerical data show good agreement with (4.3). There are some deviations at large $\bar {h}$, corresponding to very low modulation frequency. The reason is that in this regime, the time-dependent term ${\rm d}h/{\rm d}t$ cannot be neglected. However, most of the simulation data, especially for high frequency, follow the steady solution quite well.

Figure 8. Mean equilibrium solid–liquid interface height $\bar {h}$ as a function of $Nu/\theta _c$ at different $\theta _c$. The dashed line represents the derived model from the steady energy balance equation (4.2), which agrees well with the simulation data in the ‘partially solid’ regime. The hollow circles represent the cases where an obvious oscillation of interface in one period is observed, e.g. the case in figure 3(c).

5. Dependence of interface oscillation amplitude on control parameters

To obtain a relation between the oscillation amplitude $A$ of the solid–liquid interface and the control parameters $f$ and $\theta _c$ in the ‘partially solid’ regime, we look into (4.1) with the fluctuation of $h$ taken into account. The solid–liquid interface height can be expanded as

(5.1)\begin{equation} h(t)=\bar{h}+h_1(t), \end{equation}

where $h_1(t)\ll \bar {h}$. We assume that the modulation of the heat flux follows the thermal modulation at the bottom plate with a certain phase delay $\psi$, so that we can represent the heat flux modulation by a sinusoidal function $Q_l(t)=\bar {Q}(1+\varepsilon \sin (2{\rm \pi} f t+\psi ))$, where $\bar {Q}=\langle \partial _n\theta \rangle _{t,S}$ is the time-averaged normal heat flux over the melt front and time. Then by applying (5.1) and the full equation of heat flux (4.1), we obtain

(5.2)\begin{equation} St\left(\frac{{\rm d}\bar{h}}{{\rm d}t}+\frac{{\rm d}h_1(t)}{{\rm d}t}\right) =\frac{1}{\sqrt{Ra\,Pr}}\left(Q_l(t)-\theta_c\lambda\left(\frac{1}{1-\bar{h}} +\frac{h_1(t)}{(1-\bar{h})^2}+O(h^2_1(t))\right)\right). \end{equation}

If we cancel out the time-independent terms using (4.2), then the equation for $h_1$ can be obtained as

(5.3)\begin{equation} \sqrt{Ra\,Pr}\,St\,\frac{{\rm d}h_1}{{\rm d}t}=\varepsilon\bar{Q}\sin(2{\rm \pi} f t+\psi)-\frac{\bar{Q}^2}{\theta_c\lambda}\,h_1, \end{equation}

which can be rewritten as

(5.4)\begin{equation} \frac{{\rm d}h_1}{{\rm d}t}=B_1\sin(2{\rm \pi} f t+\psi)-B_2h_1, \end{equation}

where $B_1=\varepsilon \bar {Q}/(\sqrt {Ra\,Pr}\,St)$ and $B_2=\bar {Q}^2/(\theta _c\lambda \,\sqrt {Ra\,Pr}\,St)$ are both independent of $t$. By solving (5.4), and substituting the analytical solution for $h_1(t)$ into (5.1), we can obtain the full analytical solution for $h(t)$:

(5.5)\begin{align} h(t)=\bar{h}+\left(\frac{B_1}{B_2^2+(2{\rm \pi} f)^2}\,({B_2\sin(2{\rm \pi} f t + \psi)} - {(2{\rm \pi} f)\cos(2{\rm \pi} f t + \psi)}) + c\,{\rm e}^{{-}B_2t}\right)+O(\varepsilon^2), \end{align}

where $c$ is a constant and depends on the initial condition. Since we focus only on the equilibrium state, $t\rightarrow \infty$ implies ${\rm e}^{-B_2t}\rightarrow 0$. Based on (5.5), the oscillation amplitude of $h(t)$ is

(5.6)\begin{equation} A(\theta_c,f)=\frac{B_1}{\sqrt{B_2^2+(2{\rm \pi} f)^2}}\sim\frac{1}{f}. \end{equation}

Thus for $f\gg B_2/2{\rm \pi}$, we obtain the scaling relation $A\sim f^{-1}$, independent of $\theta _c$. A rough estimate from our simulation results gives $B_2\sim 10^{-4}$; therefore the assumption $f\gg B_2/2{\rm \pi}$ is valid for the large $f$ in our simulations. We plot $A$ as a function of $f$ at different $\theta _c$ from our simulations, and compare them to the scaling relation (5.6), as shown in figure 9. The amplitude from simulations shows good agreement with the prediction. The difference between simulation results and model prediction at low $f<10^{-4}$ is because $f$ is close to $B_2$. As $f$ decreases further, we expect that $A$ will reach the asymptotic value $A={B_1}/{B_2}$, which is independent of $f$, but depends on $\theta _c$. Note that the assumption that $Q_l$ is sinusoidal with a phase lag relative to the bottom plate temperature as well as the dependence of $\varepsilon$ on $f$ remains to be validated for even lower frequencies where asymmetry is observed, e.g. figure 4(d).

Figure 9. The oscillation amplitude $A=(h_{max}-h_{min})/2$ of the solid–liquid interface as a function of the modulation frequency $f$ at different $\theta _c$. The dashed line represents the model prediction from the perturbation solution (5.6) as $A=B_1f^{-1}$ with fitting $B_1=7.9\times 10^{-6}$, since $\varepsilon$ cannot be determined theoretically, which agrees well with the simulation data.

6. Conclusion and outlook

In conclusion, this study presents a two-dimensional numerical investigation of the dynamics and thermal responses of an oscillating melt–solidification front under simple harmonic heating at the bottom plate. In general, we classified two regimes (‘fully solid’ and ‘partially solid’) for the thermal response of the melting front at different modulation frequencies $f$ and different $\theta _c$ at the upper plate. We also quantify the equilibrium interface height $\bar {h}$ based on the energy balance in the system, where $Nu$ follows a trend similar to that of $Nu$ in single-phase modulated Rayleigh–Bénard convection (RBC) (Yang et al. 2020). We further derived a solution for the oscillation amplitude of the solid–liquid interface from a perturbation solution, which shows good agreement with the simulation data. Finally, we identified the two regimes in the $f\unicode{x2013}\theta _c$ plane: ‘fully solid’ and ‘partially solid’.

We presented a comprehensive investigation of phase transition dynamics under a harmonic thermal modulation coupled with turbulent thermal convection. Our work expands upon the recent numerical study of melting/freezing with steady forcings (Couston et al. 2021; Ravichandran & Wettlaufer 2021) to unsteady forcing. The frequency regime of the thermal modulation is crucial and depends on specific applications. For high-frequency applications such as pulsed electronic devices (Yang, Khandekar & Groll 2009), the regime is expected to be fully solid or partially solid with a steady layer due to the rapid pulsing frequency ($\,f$). In geophysical contexts such as melt ponds and glacier melting influenced by factors such as ocean currents (Ding, He & Xia 2022) and seasonal sunlight (Perovich & Polashenski 2012), the modulation frequency is much lower, placing basal melting in the partially solid regime with an unsteady layer. In this regime, the quasi-equilibrium solid–liquid interface undergoes deformation and oscillation within each period. Phase change material finds applications across various frequency regimes (Sharma et al. 2009), making the framework developed in this study applicable to different scenarios. The framework can also be extended to investigate other free-boundary problems, including dissolution (Davies Wykes et al. 2018; Mac Huang et al. 2020) and erosion (Ristroph et al. 2012; Amin et al. 2019), where the free-boundary condition depends on the concentration gradient and the tangential velocity.

Numerous open questions remain. A theoretical relation describing the dependence of the regime transitions on $f$ and $\theta _c$ is currently lacking, considering the complexity of turbulent convection with changing heating temperature and deformable solid–liquid interface. Exploring the flow dynamics in three dimensions and comparing it to two-dimensional RBC (van der Poel et al. 2013) would be of interest. However, conducting three-dimensional direct numerical simulations remains computationally demanding. The impact of $Ra$, $Pr$ and domain aspect ratio on the melting topography and the equilibrium statistics is still unclear and warrants further investigation; see the Appendix for a preliminary test. Note that Purseed et al. (2020) found multi-stability in RBC with a melting boundary, which has implications for the choice and impact of initial conditions (here $H/2$). We ran a series of simulations with different initial heights from $0.1H$ to $0.9H$ at $Ra=10^8$, $f=0.1$, $\theta _c=0.1$. In our results, we did not see multiple flow states, which could be because the choice of $Ra$ in our case is too large for the existence of the state of pure conduction. Multiple equilibrium states could appear for larger $\theta _c$, which decreases the interface height and the effective $Ra$. Moreover, with different starting phases of the heating force (start with heating or cooling), and initial temperature profiles (start with a dynamical state, with convection active, or with a static stable bottom boundary layer), the results can be different. Answering these questions will require a thorough exploration with simulations of different initial conditions with different control parameters. Future studies should also consider multi-component liquids such as seawater, where both temperature and salinity play significant roles in the flow structure and phase change process. The effect of fluid properties and salinity effect is also explored by a series of studies (Wang et al. 2021c,d; Wang, Calzavarini & Sun 2021b; Du et al. 2023), considering the density anomaly and mushy-layer-induced convection. Additionally, variations in salt concentration alter the temperature corresponding to the density maximum, thereby influencing flow structures substantially.