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 Reference Verzicco and Orlandi1996; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015a). It has already been validated extensively, and applied to studies of turbulent flows (van der Poel et al. Reference van der Poel, Ostilla-Mónico, Verzicco, Grossmann and Lohse2015b; Yang, Verzicco & Lohse Reference Yang, Verzicco and Lohse2016; Yang et al. Reference Yang, Chong, Wang, Verzicco, Shishkina and Lohse2020; Wang, Lohse & Shishkina Reference Wang, Lohse and Shishkina2021a). To simulate the phase transition process, we use AFiD and the phase-field method presented by Favier et al. (Reference Favier, Purseed and Duchemin2019). 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. Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993; Favier et al. Reference Favier, Purseed and Duchemin2019). The implementation and validation of the phase-field model are shown in previous work (Liu et al. Reference Liu, Ng, Chong, Lohse and Verzicco2021; Yang et al. Reference Yang, Chong, Liu, Verzicco and Lohse2022).

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 Reference Couston2021), 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. Reference Haumann2020) as well as industrial applications such as the PCMs, where the undercooling effect is also widely observed (Zahir et al. Reference Zahir, Mohamed, Saidur and Al-Sulaiman2019; Shamseddine et al. Reference Shamseddine, Pennec, Biwole and Fardoun2022).

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. Reference Favier, Purseed and Duchemin2019).

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. Reference Yang, Chong, Wang, Verzicco, Shishkina and Lohse2020), and $0.1\leq \theta _c\leq 1.4$, which is a common temperature ratio for ice and PCM (Dietsche & Müller Reference Dietsche and Müller1985). 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. Reference Sharma, Tyagi, Chen and Buddhi2009), 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 Reference James1968). 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 Reference van der Poel, Stevens and Lohse2013). 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. (Reference Yang, Chong, Wang, Verzicco, Shishkina and Lohse2020) 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. Reference Yang, Chong, Wang, Verzicco, Shishkina and Lohse2020), 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. (Reference Yang, Chong, Wang, Verzicco, Shishkina and Lohse2020).

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.