Melting of olive oil in immiscible surroundings: experiments and theory

Abstract

We report on the melting dynamics of frozen olive oil in quiescent water for Rayleigh numbers up to $10^9$. The density difference results in an upward buoyancy-driven flow of liquid oil forming a thin film around the frozen oil. We experimentally investigate flat, cylindrical and spherical shapes and we derive theoretical expressions for the local film thickness, velocity and the local melt rate for these three canonical geometries. Our theoretical models predict the correct order of magnitude and the correct scaling as compared with our experimental findings.

1. Introduction

Understanding the complicated dynamics of phase change is relevant to predicting and controlling many natural and industrial processes. Melting and dissolution are examples of the classical Stefan problem, where the boundary is defined by the phase of the material and the evolution of the boundary follows from the material undergoing phase change. Common examples include the freezing of water to make ice cubes, novel developments in phase change materials, where the latent heat of fusion is used as a temporary energy storage (Dhaidan & Khodadadi 2015), and the melting of ice around the Earth's North and South Poles (Holland, Bitz & Tremblay 2006; Feltham 2008; Cenedese & Straneo 2023).

During melting, the cold melt generally flows along the body, giving rise to non-uniform melting, and therefore changing the shape of the object, which then feeds back on the flow. This shape change or self-sculpting process of objects subject to melting, erosion or dissolution has been a topic of recent interest. The evolution of eroding clay balls and cylinders in a flow have been studied by Ristroph et al. (2012). More recently, more studies have been done with quiescent surroundings, as in Cohen et al. (2016), Davies Wykes et al. (2018), Pegler & Davies Wykes (2020) and Cohen et al. (2020), where they studied the pattern formation due to natural convection and dissolution of hard candy and salt, submerged in water. Pattern formation was also studied by Guérin et al. (2020) who, using experiments, reveal the dynamics of karst geomorphology and rillenkarren formations. Further insights into the emergence of rock formations due to dissolution are provided by Davies Wykes et al. (2018) and Huang et al. (2020), who emphasise the importance of the directionality of the shaping process. In Huang & Moore (2022) a class of exact solutions is given for the shape of the pinnacles, which show that the tip curvature is large, but finite. Recently there have been direct numerical simulations by Yang et al. (2023a,b), who use the phase-field method to study the morphology of melting ice in a Rayleigh–Bénard geometry and stratification of salt concentration around a melting cylinder.

Hitherto, most studies have focused on the melting of miscible fluids, i.e. a frozen object submerged in the same substance in liquid phase, or a similar miscible liquid or solution (e.g. melting of ice in salty water). The case of immiscible melting has received little attention. Immiscible melting can be achieved in multiple ways: either we have organic compounds like oils and waxes and combine those with water, or we have metals (e.g. gallium) inside water or oils. The most experimentally accessible option is to use an oil with a freezing point around $0\,^\circ$C and water. Motivated in the context of nuclear meltdown accidents in nuclear reactors involving molten core material, Taghavi-Tafreshi, Dhir & Catton (1979) looked at the case of melting of an immiscible liquid where the frozen oil is below warmer water. In the melting of a horizontal wall the process is governed by the Rayleigh–Taylor instability (de Gennes et al. 2004), as the melting material (oil) is positively buoyant in the water. The interface in their case is therefore heavily undulated by the pinching-off droplets. Nevertheless, they have found a scaling of the Nusselt number ($Nu$) with the Rayleigh number ($Ra$) of $Nu \sim Ra^{1/4}$, despite all the intricacies associated with the pinching-off droplets and the non-uniform liquid oil layer. These findings were followed up by Farhadieh & Epstein (1982) where they looked at a variety of oily or waxy substances melting in a variety of watery solutions. They observe from their experimental measurements that their data are bounded by two different $Nu \propto Ra^{1/4}$ scaling laws, but that their data follow more closely a ${Nu} \propto {Ra}^{1/5}$ scaling set forth by Gerstmann & Griffith (1967) in terms of absolute agreement, although they did not make any definitive conclusion on the scaling exponent due to the scatter in their experimental findings. In this work we study the melting process of frozen olive oil in water around room temperature, as this is most accessible experimentally. We will, however, focus on configurations where the Rayleigh–Taylor instability only plays a minor role.

For thermal convection problems with phase change the three dimensionless control parameters are the Rayleigh (${Ra}_o = {g(\rho _w-\rho _o)L^3}/{(\nu _o\alpha _o\rho _o)}$), Stefan (${Ste} = {c_p\Delta T}/{\mathcal {L}}$) and Prandtl (${Pr}=\nu _o/\alpha _o$) numbers, where $g$ is gravity, $\rho _w$ is the density of the water, $\rho _o$ is the density of the olive oil, $L$ is a relevant length scale, $\nu$ is the kinematic viscosity, $\alpha$ is the thermal diffusivity, $c_p$ is the specific heat, $\Delta T = T_\infty - T_o$ is the temperature difference between the ambient and the melting temperature and $\mathcal {L}$ is the latent heat of fusion. We have used $o$ and $w$ subscripts to denote quantities related to oil and water, respectively. Here, the Rayleigh number is the time scale associated with thermal transport due to diffusion, as compared with the time scale of thermal transport due to convection. Whereas in classical Rayleigh–Bénard convection buoyancy is created by (generally) small density changes due to temperature changes ($\beta \Delta T$), in our case a large density difference is immediately created due to the different substances $((\rho _w-\rho _o)/\rho _o)$. A high Rayleigh number means intense thermal driving of the system. The Stefan number describes the ratio of specific heat vs the latent heat of fusion, where the latent heat is the heat needed or released by a phase change. A higher Stefan number means that phase change happens faster. The Prandtl number is a material property describing the ratio of the momentum diffusivity to the thermal diffusivity, and determines whether the thermal boundary layer is embedded in the momentum boundary layer or vice versa.

The present work has the following structure: in § 2 we describe the experimental set-up. In § 3, results for the melt rate of frozen olive oil are shown for three different geometries: a vertical wall, a cylinder and a ball. For the cylinder we show two different initial Rayleigh numbers (initial sizes). The obtained local melt rates are compared with theoretical models that are derived in detail in § 4. The effect of the assumption of constant viscosity is discussed, and a correction for the variation of the viscosity with temperature is derived. Lastly, we discuss our findings in detail in § 5 and finish with our conclusions in § 6.

2. Experimental set-up

The schematic in figure 1 shows the experimental set-up we use to study the melting of frozen olive oil. We choose olive oil as our working fluid since the freezing/melting temperature is reasonable to achieve in the laboratory ($T_o = -8\,^\circ$C), the resulting contact temperature is above the freezing temperature for the ambient water and it is readily available. We use a rectangular glass tank of $400\ {\rm mm} \times 500\ {\rm mm} \times 800\ {\rm mm}$, filled with water. During the melting process, the water is quiescent and kept at a constant room temperature $T_\infty = 20\,^\circ$C. We only investigate here for ambient temperatures equal to this room temperature, so as to avoid any spurious flows created by natural convection caused by heat transfer from the surroundings to the bath. A small PVC holder (poor thermal conductivity) is included in the frozen olive oil during the freezing process and is attached to a support that is submerged in water. The olive oil is cooled down to a temperature of $T_i = -14\,^\circ$C. Further material properties of both substances can be found in table 1, including the surface tension. Since $\rho _{oil} < \rho _{water}$ the melted olive oil will rise and collect at the apex of the objects where it periodically pinches off. There, the surface tension is relevant for the pinch-off dynamics. We do not study this effect here and will ignore the surface tension, which we elaborate on in the discussion section. For studies of droplet deformation and breakup in viscous fluids we refer to e.g. Stone (1994), Cohen et al. (1999) and Cohen & Nagel (2001). We study three different canonical geometries: a vertical wall, a horizontal cylinder and a ball; see figure 1. The melting process is recorded through interval imaging. For this, a DSLR camera (Nikon D850) with a 100 mm macro objective (Zeiss Makro Planar $T^*$ 2/100) is used, resulting in a resolution of $30\ \mathrm {\mu }{\rm m}\ {\rm px}^{-1}$. An LED light source and light diffuser are used to create a uniformly lit background. The original movies can be found in the supplementary material available at https://doi.org/10.1017/jfm.2024.843. The images are binarised after which we find the contour, area and, for the cylinders and balls, the centroid of the object, see figures 2(b) and 2(c). From the evolution of the contour we find the local melt rate.

Figure 1. Schematic of the experimental set-up. The frozen olive oil object is submerged in quiescent water. The dimensions of the glass tank are $400\ {\rm mm} \times 500\ {\rm mm} \times 800\ {\rm mm}$. A PVC holder is incorporated in the frozen oil during freezing, and connects the frozen oil to a support. A white LED light source is used with a light diffuser to create a uniform background illumination. Three canonical geometries are shown: vertical wall (photographed), horizontal cylinder and ball. A camera periodically photographs the melting objects from the side (vertical wall) or the front (cylinder and ball).

Table 1. Material properties of water and olive oil. Properties for water are taken at $T = 10^\circ$C, following Bejan (1993). Thermal properties for olive oil are taken from Carbajal Valdez et al. (2006) and Turgut, Tavman & Tavman (2009). Values are taken for a constant temperature of $3\,^\circ$C; the mean temperature in the olive oil melt layer. The value for the latent heat of fusion is taken from EasyCalculation.com (2005). Comparing with values for other vegetable oils from Gudheim (1944) we find that they are of comparable magnitude. For quantities given without error estimates in their source we have assumed the last digit to indicate the precision. Note that the thermal properties of olive oil are not entirely self-consistent: $\alpha \rho c_p/\lambda = 0.82 \pm 0.04$. This is due to the fact that not many different sources are available, and measurements were done on different olive oils, possibly with different compositions. Wherever one of these quantities is needed, we opt for the ‘most direct’ quantity. $^*$The surface tension $\sigma$ was measured using the pendant drop method for olive oil submerged in water at a temperature of $20\,^\circ$C.

Figure 2. (a) Melting of a horizontal cylinder in ambient water showing several stages of the droplet pinch-off at the top where the melted olive oil collects. Note that we consider the melting cylinder only in two-dimensional cross-section in this work. The distance between the pinching-off droplets is $22.5\ {\rm mm} \pm 2\ {\rm mm}$. (b) Example image from a melting ball experiment with the contour that is obtained from image analysis overlaid on the original image. Note that the region where the liquid olive oil melt collects and periodically pinches off from the ball is ignored. (c) Contours of a melting ball for $t=200$ to 3200 s, with intervals of 200 s. Note that the tracking on the apex of the ball is hindered by the collecting and periodic pinching-off of liquid olive oil melt. Therefore, data are ignored in the grey shaded region.

The image processing is applied to all images that are taken during an experiment, typically with an interval time of $\Delta t = 20$ s. In figure 2(c) contours from a single experiment are shown at different times. Such an image shows a qualitative description of the melting process of a ball. The contours at the top are more closely spaced, whereas contours at the bottom are more distant, revealing that the melt rate at the bottom of the ball is higher than the melt rate at the top.

3. Results

Here we will look at the melt rates obtained from the experiments. We will then compare these with analytical expressions – derived in the next section – and discuss the applicability of the theory for the vertical wall, for the two horizontal cylinders and lastly for a ball.

3.1. Vertical wall

We first look at the case of the melting of a block of olive oil of 30 cm height; see figure 6 for the definitions of the axes and other quantities. We calculate the horizontal local melting rate from the evolving contours; see figure 3. For the profiles in the early stages of the melting process, the shape, despite slightly changing over time, can be regarded as vertical. At later times, the profile of the initially rectangular block has sculpted itself away from its rectangular shape. While the total process of melting takes approximately 40 min, here, we just show melt rates obtained during the first 5 min. In the lower regions of the vertical wall ($z < 2$ cm), finite-size effects of the object are observed (the bottom corner is rounded over time). We do not show the upper edge region of the melting wall, since the results are heavily influenced by the accumulation and detachment of oil droplets. Away from the top and bottom corners it can be seen that the melt rates are remarkably constant over time, in both scaling ($U\propto z^{-1/4}$) and magnitude. The grey dotted line shows the theoretical model with constant viscosity in the melt layer, whereas the black dashed line shows the theoretical model where the viscosity $\mu (T)$ varies in the melt layer (see § 4.2). We find that our analytical expression predicts the correct order of magnitude and the correct scaling in this region. Note that our model does not contain any fitting parameters, and even though several approximations have been made, there is order of magnitude and scaling agreement between the model and the experiments. The inclusion of the variable viscosity lowers the prediction of the melting rate and therefore improves the prediction. Henceforth, we only show models with temperature-dependent variable viscosity in the melting layer.

Figure 3. Horizontal melt rate as a function of the height for the case of a vertical wall for ${Ra}=O(10^9)$. We show three theoretical approximations that are derived in a later section. Here, the theoretical approximation with constant viscosity $\mu _o$ is shown as a grey dotted line, the black dashed line assumes a linear viscosity profile in the melt layer and the grey dash-dotted line assumes variable viscosity resulting from a linear temperature profiles in the melt layer. The intermediate part of the profile shows a scaling of $-1/4$ with the height. Profiles for times from 100 until 380 s from the start of the melting process are shown. At later times the vertical wall is no longer an appropriate approximation. See the supplementary material for a movie of this experiment.

3.2. Cylinder

We perform experiments for horizontal cylinders with initial radii 25 mm and 60 mm corresponding to Rayleigh numbers ${Ra} \approx 10^7$ and ${Ra} \approx 10^8$, see figure 10 for the definition of the axes and other measurements. Figure 4 shows the radial melt rates $U(\theta ) = {\rm d} r(\theta )/ {\rm d} t$ for the small and large cylinders, as a function of the polar angle $\theta$, where $\theta = 0^\circ$ is at the bottom of the cylinder. For the analysis we restrict ourselves to the early stages of melting, when the assumption of a circular cross-section is still valid, since this is one of the key geometrical assumptions in the theoretical model. The theoretical model, without any fitting parameter, for the cylindrical geometry is included as a black dashed line. Here, the temperature-dependent viscosity is included in our model. A reasonable agreement is found between our experiments and our model for both cylinders in terms of order of magnitude and shape. There are some notable differences between the model and the experimental observations. The melting process of the cylindrical shape has a maximum melt rate along the surface at an angle of $\theta \approx 60^\circ$ from the bottom, whereas the theory predicts a monotonic decrease with increasing $\theta$, such that the predicted maximum melt rate is at the bottom. For $\theta \geq 160^\circ$ the theory and experiments do not match due to the collection of melt at the top before pinching off (Shi, Brenner & Nagel 1994) at the top of the cylinder and subsequent rise to the water's surface. The small cylinder has a higher melt rate than the large cylinder, which we can precisely predict from our theory since $U \propto R^{-1/4}$ ((4.58)) gives us a ratio of 1.24 in the melt rates which we also observe in our experiments within the experimental error (the value on the ordinate is just below $10\ \mathrm {\mu }{\rm m}\ {\rm s}^{-1}$ for the $R_0 = 25\ {\rm mm}$ cylinder and $\approx 8\ \mathrm {\mu }{\rm m}\ {\rm s}^{-1}$ for the $R_0 = 60\ {\rm mm}$ cylinder, giving a value of $\approx 1.2$, very close to our predicted value).

Figure 4. (a) Cross-sectional melt rate as a function of the angle $\theta$ (polar angle starting from the bottom, see figure 10) for a small cylinder ($R_0 = 25\ {\rm mm}$) ${Ra}=O(10^7)$. The black dashed line shows the theoretical model for the cylindrical geometry. (b) Large cylinder $(R_0 = 60\ {\rm mm} )\ {Ra}=O(10^8)$. See the supplementary material for movies of these experiments.

3.3. Ball

Finally, we look at the melting of a ball with initial radius $R=60\ {\rm mm}$; see figure 2(b) for the experiment and figure 10 for the definition of the axes and other measurements. Figure 5 shows the melt rate for this ball and compares the experiments with the theory. Again, we restrict ourselves to the early stages of melting, when the assumption of a circular cross-section is still valid. Note that we also included a temperature-dependent viscosity in our model. As for the cylinders, the model shows reasonable agreement in terms of scaling and shape. Deviations for $\theta \leq 75^\circ$ may be caused by the ambient water, which we will discuss in § 5. For high angles ($\theta \geq 160^\circ$) the oil layer is much thicker and the flow is influenced by periodically detaching droplets. In the dissolution-based problem discussed in Davies Wykes et al. (2018), the interface evolution differs greatly from what we show here in figure 2. Their sugar balls dissolve into water, creating a Rayleigh–Taylor-type instability similar to our work. In our present work the instability is delayed due to the thin film of very viscous olive oil (as compared with water), the surface tension and the curved surface, all modifying the instability. It is still unstable but with different wavelengths and growth rates, see e.g. de Gennes et al. (2004), Trinh et al. (2014) and Balestra, Brun & Gallaire (2016). From the dispersion relation for semi-infinite bodies from e.g. Mikaelian (1996) we get a positive growth rate only for length scales from 5 cm and larger. While our objects are bigger than $\approx$5 cm, the instability does not have sufficient time to grow to be apparent. We strictly see the pinch-off at the apex of our objects. For very large objects that are submerged in baths with smaller temperature differences (slow flow rates) this effect might be important.

Figure 5. Cross-sectional inward melt rate of a ball as a function of the angle $\theta$ for ${Ra}= O(10^8)$. Note that $\theta = 0^\circ$ is at the bottom of the ball. The solid green lines are results from the experiments. The black dashed line is the theoretical model for the spherical geometry. Note that, now, we only show the model with temperature-dependent viscosity. The reproducibility of this type of experiment is shown in Appendix B. See the supplementary material for a movie of this experiment.

4. Theory

In this section we will derive the analytical models for the melt rate of the frozen olive oil objects that were shown in figures 3–5. We start with the theory for the vertical wall. After deriving this model we see, from comparison with the experimental results, that we need to include a temperature-dependent viscosity, and we incorporate this in the analysis. Using this as a starting point we then derive the theory for the cylindrical geometry and the spherical geometry.

4.1. Vertical wall

We start with the melting of the vertical wall, as shown in figure 6. With horizontal and vertical axes $x$ and $z$, respectively, the frozen wall is along the $z$-axis at $x=0$, the liquid melt layer is in the range $0 \leq x \leq h(z)$ while the ambient water, causing the oil to melt, stretches from $x=h(z)$ to infinity in the $x$-direction. The surrounding water is cooling down as the oil is melting and therefore flows downward under the influence of gravity, reaching velocities of the order of cm s$^{-1}$. These velocities are relatively low and, as such, we will assume that the ambient water is stationary. Inside the melt layer we have the horizontal velocity $u$ and the vertical velocity $w$ in the $x$ and $z$ directions, respectively. Gravitational acceleration is $g$ in the negative $z$ direction. The most important are the densities $\rho _o$ and $\rho _w$ and the dynamic viscosities $\mu _o$ and $\mu _w$. Under these circumstances, where the oil is very viscous, and assuming we are in a steady state, the $w$-component of the Navier–Stokes equations simplifies such that we have a balance between the pressure gradient due to buoyancy and the viscous forces:

(4.1)\begin{equation} \mu_o \left(\frac{\partial^2w}{\partial x^2} + \frac{\partial^2w}{\partial z^2}\right) - g(\rho_o - \rho_w) = 0.\end{equation}

Here, we ignore surface tension effects, which we show in the discussion section to be a valid assumption. At the interface between solid and liquid oil we have $w(x=0) = 0$. The film has a thickness of $O$(mm), which is very small with respect to the height (30 cm), such that we are in the thin film limit (${\partial }/{\partial x} \gg {\partial }/{\partial z}$) and we can therefore neglect the derivatives in the $z$-direction in (4.1). At the oil–water interface the stresses have to be balanced across the interface:

(4.2)\begin{equation} \mu_o \left. \frac{\partial w}{\partial x} \right|_{x=h^-} = \mu_w \left. \frac{\partial w}{\partial x} \right|_{x=h^+}, \end{equation}

where $h^-$ ($h^+$) indicates the gradient at the interface on the oil (water) side of the oil–water interface. Since our oil is much more viscous as compared with our water (table 1: $\mu _o/\mu _w = O(100)$) and there is no reason to expect any sharp gradients in the water, we have for the velocity gradient at the oil–water interface

(4.3)\begin{equation} \left. \frac{\partial w}{\partial x} \right|_{x=h^-} = \frac{\mu_w}{\mu_o } \left. \frac{\partial w}{\partial x} \right|_{x=h^+} \approx 0. \end{equation}

The solution of a simplified (4.1) obeying these boundary conditions is

(4.4)\begin{equation} w(x,z) = \beta x(2h(z)-x) \quad \text{with } \beta = \frac{g (\rho_w - \rho_o)}{2\mu_o} = \frac{g\Delta \rho}{2\mu_o}.\end{equation}

From the continuity equation follows the horizontal velocity:

(4.5)\begin{equation} \frac{\partial u}{\partial x} ={-}\frac{\partial w}{\partial z} ={-}2\beta x\frac{{\rm d} h}{{\rm d} z}. \end{equation}

We can now find the melt rate $U(z)$ at which the wall melts by integrating the previous equation to obtain

(4.6)\begin{equation} u(x,z) = U(z) - \beta x^2\frac{{\rm d} h}{{\rm d} z}. \end{equation}

From mass conservation, see figure 7, we can relate the melt rate and the film thickness:

(4.7)\begin{equation} U(z) = \frac{{\rm d} }{{\rm d} z}\left(\int_0^h w\, {\rm d} x\right) = \frac{{\rm d} }{{\rm d} z}\left(\frac{2}{3}\beta h(z)^3\right).\end{equation}

To find $h(z)$ we need to consider the thermal transport. The ambient water, with a temperature of $T_\infty = 20\,^\circ$C far away, transfers heat to the melt layer. This, in turn, transfers heat to the solid oil, causing this to melt further. We first consider the advection–conduction equation inside the thin melt layer. Assuming stationarity and using the thin film approximation we arrive at

(4.8)\begin{equation} u\frac{\partial T}{\partial x} + w \frac{\partial T}{\partial z}= \alpha\frac{\partial^2T}{\partial x^2}, \end{equation}

where $\alpha = (\lambda /(\rho c_p))_o$ is the thermal diffusivity, $\lambda$ is the thermal conductivity and $c_p$ is the specific heat capacity. In order to make the following analysis tractable (the Mises transformation (Schlichting & Gersten 2016, p. 183), where instead of $x$ and $z$ the coordinates $z$ and the streamfunction $\psi$ are used, leads to ${\partial T}/{\partial z} = \alpha ({\partial }/{\partial \psi }) ( w ({\partial T}/{\partial \psi }))$, which is no more tractable), we focus on the first term and neglect the second term of the left-hand side of (4.8). For the first term we approximate $u$ by only taking the first term of (4.6):

(4.9)\begin{equation} U(z)\frac{\partial T}{\partial x} = \alpha\frac{\partial^2T}{\partial x^2}.\end{equation}

After the result for $T(x,z)$ is obtained we will estimate the involved error in these approximations. One boundary condition is that $T$ equals the melting temperature $T_o$ at $x = 0$. At the interface of the melt layer and the ambient water, $x=h$, temperature and heat flux must be continuous. Since the velocity is small in both the oil and water, heat conduction is prominent. At the boundary, the water flowing down is only in contact with the wall for a short time (see the discussions in § 5). Therefore, we approximate here the temperature at the interface between the water and oil with the so-called contact temperature, occurring when two semi-infinite media with different temperatures are brought in contact. With material properties $f = \sqrt {\lambda \rho c_p}$, this contact temperature is then (see e.g. § 5.7, (5.63) from Incropera & De Witt 1990)

(4.10)\begin{equation} T^* = \frac{f_o T_o + f_w T_\infty}{f_o + f_w}, \end{equation}

where $T_\infty$ is the temperature of the water far away. Filling in the values for water and oil we get

(4.11)\begin{equation} T^* = T_o + 0.75(T_\infty - T_o). \end{equation}

The solution of (4.9) with $T_o$ at $x = 0$ and $T= T^*$ at $x=h$ results in

(4.12)\begin{equation} T(x,z) = T_o + 0.75 (T_\infty - T_o) \frac{\displaystyle \int_0^x {\rm e}^{{U(z)x}/\alpha}\,{\rm d} x}{\displaystyle \int_0^h {\rm e}^{{U(z)x}/\alpha}\,{\rm d} x}.\end{equation}

The heat flux at the wall, $\lambda ({\rm d} T/{{\rm d} x})_{x=0}$, results in the melting rate. With latent heat $\mathcal {L}$ this means that

(4.13)\begin{equation} \lambda({\rm d} T/{{\rm d} x})_{x=0} = U(z)\rho_o\mathcal{L}.\end{equation}

Combining (4.12) and (4.13) gives

(4.14)\begin{equation} \frac{0.75\lambda(T_\infty - T_o)}{\displaystyle \int_0^h {\rm e}^{{U(z)x}/\alpha}\,{{\rm d} x}} = U(z) \rho_o \mathcal{L}. \end{equation}

We now rewrite the expression and introduce the non-dimensionalisation $\xi = x/h$:

(4.15)\begin{equation} 0.75\lambda(T_\infty - T_o) = U(z) \rho_o \mathcal{L} h \int_0^1 {\rm e}^{{U(z)h \xi}/\alpha}\, {\rm d} \xi. \end{equation}

Working out the integral and rearranging gives

(4.16)\begin{gather} 0.75\lambda(T_\infty - T_o) = \rho_o \mathcal{L} \alpha [ {\rm e}^{{U(z)h}/\alpha} -1 ], \end{gather}

(4.17)\begin{gather}\alpha \ln\left( 1+ \frac{0.75\lambda(T_\infty - T_o)}{\rho_o \mathcal{L} \alpha} \right) = U(z) h =\frac{\beta}{2}\frac{{\rm d}}{{\rm d} z}(h^4), \end{gather}

where in the last step we used (4.7). We can now solve (4.17) for $h$ without further approximation:

(4.18)\begin{gather} h(z) = z^{1/4} \left[ \frac{2\alpha}{\beta} \ln(1+\varLambda) \right]^{1/4}, \end{gather}

(4.19)\begin{gather}\varLambda = 0.75(T_\infty - T_o)\frac{c_p}{\mathcal{L}} = 0.75 Ste. \end{gather}

We note that the dependence of the bath temperature on the layer thickness $h$ is very weak, not only is there a logarithmic dependence, but also the fourth root would make experimentally verifying the dependence of $T_\infty$ on $h$ very challenging. Now that we have found an expression for $h$, $T$ and $U$ we can consider the approximations we have made before. We first focus on neglecting the term $w({\partial T}/{\partial z})$ in (4.8). Since $Uh/\alpha$ is in all our experiments small with respect to unity we may, from (4.12), for an estimate we take $T-T_o$ as $C x/h$ with $C$ a constant. The main term, from (4.6), $U(z)=2\beta h^2 ({{\rm d} h}/{{\rm d} z})$. Hence, the term $u({\partial T}/{\partial x})$ connected herewith is $2C\beta h ({{\rm d} h}/{{\rm d} z})$. The term $w({\partial T}/{\partial z})$ is smaller, using (4.4), by an amount $\xi ^2(2-\xi )/2$. This is significantly, although not negligibly, smaller; the error is $18.75\,\%$ at $\xi =1/2$. To obtain an estimate of the temperature distribution this is acceptable. We next consider the approximation of $u$ in (4.6). We now use in (4.12), instead of $U$, the full expression in (4.6) and carry out the integration of the thus obtained version of (4.14) numerically, inserting the solution ((4.18)) for $h$. We find that the value of the integral then differs only $0.6\,\%$, owing to the fact that the argument of the exponent is $\ll 1$ for both cases (with and without the second term in (4.6)), such that the exponent is roughly constant $\approx 1$.

Figure 6. Schematic overview of melted olive oil (yellow) at a flat plate of frozen olive oil (orange), submerged in water. Here, $U(z)$ is the melt rate at the frozen/liquid oil interface.

Figure 7. Illustration of the mass flows in a small control volume (dashed grey) of melting olive oil at a vertical wall. Oil is ‘injected’ into the volume from the wall since the left of the control volume is following the interface, and oil enters the volume from below due to buoyancy forces. The total ingress from these two contributions equals the egress at the top of the control volume.

Finally, we consider the neglect of convection in the transfer of heat between the oil film and water. In this connection we can consider the melt oil film as a solid since the velocity, as figures 3 and 5 show, is of order of $10\ \mathrm {\mu }{\rm m}\ {\rm s}^{-1}$, which is negligibly small. Water is cooled by the cold wall and flows downward. We can make an estimate of the velocity and the associated heat transfer by making use of the data on the opposite case, a hot wall and a cold fluid, since that is extensively dealt with in the literature. Following Bejan (1993) (pp. 345–346) for the case of the vertical wall, we define the Rayleigh number by ${Ra}_w = \epsilon _w g ( T_\infty - T^* ) (L/2)^3 {Pr}_w/\nu _w^2$, where, in addition to the quantities in table 1, $\epsilon =1.7\times 10^{-4}\ {\rm K}^{-1}$ is the thermal expansion coefficient of water at $16.5\,^\circ {\rm C}, \nu _w=\mu _w/\rho _w$ the kinematic viscosity of water and ${Pr}=\nu /\alpha$ the Prandtl number of water. Here, ${Ra}_w$ depends on the location along the wall, for which we have taken the half-way height $L/2=0.15$ m. Using these quantities and table 1, we find ${Ra}_w=2.5\times 10^8$. With (7.43) from Bejan (1993) (p. 345) and $G$ from figure 7.5 (Bejan 1993, p. 346) the downward velocity is 1.04 cm s$^{-1}$. In this calculation the estimate of a water velocity of the order of cm s$^{-1}$, made earlier in this paper, was used. The corresponding Reynolds number (${Re}=vL/\nu = 2800$) is far below $O(10^5)$ where the flow becomes turbulent. Bejan (1993) ((7.51)) gives, for the above mentioned ${Ra}_w$, a Nusselt number ${Nu}=KL/(2\lambda _w)=58$, where $K$ is a heat transfer coefficient, with a value of approximately $225\ {\rm W}\ ({\rm m}^{2}{\rm K})^{-1}$. The associated heat flux from the oil film to the water is then $1.6\ {\rm kW}\ {\rm m}^{-2}$. This is small with respect to the heat flux by conduction which is $\lambda _o (T^*-T_o)/h$. With $h\approx 0.3\ {\rm mm}$ from ((4.18)) this is ${\approx }12\ {\rm kW}\ {\rm m}^{-2}$. An additional and strong argument for the neglect of free convection is the following: consider (4.17). The quantity $0.75(T_\infty - T_o)$ is $T^* - T_o$. In our theory this is a constant. If free convection where important $T^*$ would depend on $z$ and, as (4.17) shows, $h$ would not vary as $z^{1/4}$, but then deviate from that.

For comparison with experiments, $U(z)$ is the best quantity. From (4.7) (or, alternatively, (4.17)) and (4.18) we have

(4.20)\begin{equation} U(z) = \frac{\alpha}{h}\ln{(1+\varLambda)} \propto z^{{-}1/4}, \end{equation}

which scales with the height $z$ to a power of $-1/4$. This exponent has been shown before in similar configurations by Wagner (1949), Ostrach (1953), Merk & Prins (1954) and Wells & Worster (2011). We note that the melt rate only has a weak dependence on the bath temperature $T_\infty$. Combining (4.18) and (4.20), we see that $U$ scales as the logarithm of $\varLambda$ taken to the power $3/4$. To achieve a doubling of our $U$ therefore means increasing our bath temperature to ${\approx }71\,^\circ {\rm C}$, which would create all kinds of problems with spurious flows created by heat leak with the surroundings causing natural convection. The predicted melt rate in (4.20) is drawn in figure 3 as the grey dotted line. This shows that the slope agrees satisfactorily with the experimental data, but the values are considerably higher. We have now estimated the effect of three approximations viz. the reduction of $u$ to $U$, the neglect of $w\partial T/\partial z$ in (4.8) and the neglect of free convection at the water side of the melt film. None of these seems dominant as the source of the discrepancy between the theory and the measurements. Initially, we suspected the neglect of the vertical heat transport to be the main cause. We made another estimate of the involved error by first writing (4.8) as

(4.21)\begin{equation} \frac{\partial}{\partial x}(u T ) + \frac{\partial}{\partial z}(w T ) = \alpha \frac{\partial^2 T}{\partial^2 x}, \end{equation}

and subsequently integrating this term by term over the melt layer, using the solution of (4.12) for the temperature and $U$ and (4.4) for $u$ and $w$. Whereas the first term on the left and the conduction term both give $U(T^*-T_o)$, the second term on the left, the neglected term, results in $(U/8)(T^* + 3T_o )$. With $T^*=13\,^\circ$C and $T_o=-8\,^\circ$C, this term has a value of $-11/8\ {\rm Km}\ {\rm s}^{-1}$, in absolute magnitude only 6 % of $U(T^*-T_o)$. This again does not mark the neglect of vertical heat transport as causing the discrepancy.

So far we have approximated the viscosity of the oil as constant inside the melt layer, however, the strong temperature gradient inside the melt layer does not allow us to model the viscosity as constant. The temperature dependence of the viscosity of the oil is shown in figure 8 for the relevant temperature range. The temperature at the wall is $-8\,^\circ {\rm C}$ and the temperature at the oil–water interface is $T^* = 13\,^\circ$C. The viscosity varies from 380 to 100 cP over this interval. In the following section we calculate again the vertical velocity $w$ in the melt, however, now taking a variable viscosity into account.

Figure 8. Olive oil viscosity as a function of the temperature. Measurements are done using an Anton Paar MCR502 rheometer with Peltier cooling. Measurements are done starting at a temperature of $20\,^\circ$C and decreasing step-wise to a minimum of $-8\,^\circ$C. We estimate an error of ${\pm }5$ cP.

4.2. Variation in viscosity

In table 1 we state a value for the viscosity of olive oil of 170 cP. This value is taken at a mean olive oil temperature of $3\,^\circ {\rm C}$. An Anton Paar MCR502 rheometer was used to measure temperature dependence of the viscosity of the olive oil, see figure 8. As the oil is cooling down and approaching its freezing temperature its viscosity is increasing substantially. To account for the variation in the viscosity, we recalculate the vertical velocity $w$. In order to do this, we have to consider the Navier–Stokes equations where constant viscosity is not assumed and from that we obtain (different from (4.1))

(4.22)\begin{equation} \frac{\partial}{\partial x}\left(\mu_o(x)\frac{\partial w}{\partial x} \right) ={-}g\Delta \rho. \end{equation}

Crucially, here, we see that the viscosity is now inside the first derivative. After integrating once with respect to $x$ and rearranging we obtain

(4.23)\begin{equation} \frac{\partial w}{\partial x} = \frac{-g\Delta \rho x}{\mu_o(x)} + C.\end{equation}

The boundary conditions remain unchanged:

(4.24)\begin{gather} w|_{x=0} = 0, \end{gather}

(4.25)\begin{gather}\left.\frac{\partial w}{\partial x}\right|_{x=h} = 0. \end{gather}

We find the integration constant $C$ from using the boundary equation (4.25)

(4.26)\begin{equation} \frac{\partial w}{\partial x} = \frac{g\Delta \rho (h-x)}{\mu_o(x)}. \end{equation}

This cannot be further solved analytically without assuming a particular shape for $\mu _o(x)$. We will follow two approaches. For the first, we will assume that the viscosity varies linearly inside the thin film to find the velocity profile $w$, which can be done fully analytically. For the second approach we will assume a linear temperature profile inside the thin film, and by numerically integrating (including the viscosity curve of figure 8), we find the velocity profile.

For the first approach we will assume that the viscosity is varying linearly inside the thin film:

(4.27)\begin{equation} \mu_o(x) = \mu_o(T_o)\left(1-\frac{x}{h}\right) + \mu_o(T^*)\frac{x}{h}. \end{equation}

We introduce the dimensionless quantity $\tilde {\mu }$:

(4.28)\begin{equation} \tilde{\mu}= \frac{\mu_o(T^*) - \mu_o(T_o)}{\mu_o(T_o)}. \end{equation}

Using this definition, (4.26) can be rewritten and we obtain

(4.29)\begin{equation} \frac{\partial w}{\partial x} = \frac{2 \beta^* h (h-x)}{h+\tilde{\mu} x}, \end{equation}

where $\beta ^*$ is similar to $\beta$ ((4.4)) but with $\mu _o = \mu _o(T_o)$. After integrating and applying the boundary condition (4.24) we get

(4.30)\begin{equation} w(x,z) = 2 \beta^* h^2 \left(\frac{\tilde{\mu}+1}{\tilde{\mu}^2} \ln \left(\frac{\tilde{\mu} x}{h}+1\right)-\frac{x}{h \tilde{\mu}}\right). \end{equation}

We see by expansion in $\tilde \mu$ that, for $\tilde \mu = 0$, $w$ is the same as in (4.4). Note that $w$ is still a function of both $x$ and $z$, since $h$ is a function of the height. As before, to get the melt rate $U(z)$ we integrate over a control volume in the liquid melt layer (see (4.7)):

(4.31)\begin{equation} U(z) = \frac{{\rm d}}{{\rm d} z}\left(\beta^* \frac{2}{3}h(z)^3\right)\left(\frac{3 (1+\tilde{\mu})^2}{\tilde{\mu}^3} \ln (1+\tilde{\mu}) -\frac{6 + 9 \tilde{\mu}}{2 \tilde{\mu}^2}\right) .\end{equation}

Comparing this result with the previously found expression for melt rate $U(z)$ (right-hand side of (4.7)), it is seen that the varying viscosity introduces a correction to the melt rate that is dependent on the values for the viscosity at the wall and at the oil–water interface. Using (4.17) we see that the ratio between the melt rates for the constant viscosity vs the linearly varying viscosity is $\beta ^*/\beta$ times the fourth root of the last bracketed term in (4.31).

For the second approach we will assume that the temperature is varying linearly inside the thin film:

(4.32)\begin{equation} T(\xi) = \xi T^* + (1-\xi) T_o. \end{equation}

We can now fill this into (4.26) and integrate to find $w$:

(4.33)\begin{gather} \frac{\partial w}{\partial x} = \frac{g\Delta \rho (h-x)}{\mu_o(T(\xi))}, \end{gather}

(4.34)\begin{gather}\frac{w(\xi)}{h^2} = g\Delta \rho \int \frac{1-\xi}{\mu_o(T(\xi))}\, {\rm d} \xi +C. \end{gather}

The curve $\mu _o(T)$ from figure 8 is used to numerically integrate the profile and, using the boundary condition ((4.24)), the integration constant can be determined. The comparison of all these velocity profiles can be found in figure 9. Going from a constant viscosity to a viscosity that varies linearly dramatically decreases the velocity $w$ and stays relatively close to a parabolic profile. If we now vary the viscosity according to the curve of figure 8 (on average this viscosity is lower than the linear approximation) the curve is slightly higher but the shape of the profile is still roughly parabolic for the viscosities under consideration.

Figure 9. Vertical velocity profiles assuming constant viscosity (grey dotted, (4.4)), linear viscosity (black dashed, (4.30)) and linear temperature profile and viscosity following figure 8 (grey dash-dotted, (4.34)).

We can now find $U$ for the general case of $\mu (T(\xi ))$ by rewriting (4.7):

(4.35)\begin{equation} U(z) = \frac{\partial}{\partial z}\left(h^3\int_0^1 \frac{w}{h^2}\, {\rm d} \xi\right) = \frac{\partial}{\partial z}(h^3 ) \int_0^1 \frac{w}{h^2}\, {\rm d} \xi, \end{equation}

where we note that the integral is only a function of $\xi$ (see (4.34)) and evaluates to a constant, and can therefore be taken outside the derivative. We evaluate now $hU(z)$,

(4.36)\begin{equation} h U(z) = \frac 34 \frac{\partial}{\partial z}(h^4 ) \int_0^1 \frac{w}{h^2}\, {\rm d} \xi, \end{equation}

which we can now insert into (4.15) to find an equation similar to (4.17):

(4.37)\begin{equation} \alpha \ln\left( 1+ \frac{0.75\lambda(T_\infty - T_o)}{\rho_o \mathcal{L} \alpha} \right) = U(z) h = \frac 34 \frac{\partial}{\partial z}(h^4 ) \int_0^1 \frac{w}{h^2}\, {\rm d} \xi. \end{equation}

We can solve for $h$ (again realising that the last integral evaluates to a number and does not depend on $z$ nor $h$):

(4.38)\begin{equation} h(z) = z^{1/4} \left[ \frac{\alpha\ln(1+\varLambda)}{\displaystyle \dfrac 34 \int_0^1 \dfrac{w}{h^2}\, {\rm d} \xi} \right]^{1/4}. \end{equation}

We now plug this in (4.35) to obtain

(4.39)\begin{gather} U(z)= \frac{\partial}{\partial z} z^{3/4} \left[ \frac{\alpha\ln(1+\varLambda)}{\displaystyle \dfrac 34 \int_0^1 \dfrac{w}{h^2}\, {\rm d} \xi} \right]^{3/4} \int_0^1 \frac{w}{h^2}\, {\rm d} \xi, \end{gather}

(4.40)\begin{gather}U(z)= \frac{1}{z^{1/4}} [\alpha\ln(1+\varLambda) ]^{3/4} \left[\frac 34 \int_0^1 \frac{w}{h^2}\, {\rm d} \xi \right]^{1/4} . \end{gather}

One can see that (4.40) simplifies to (4.20) once we take $\mu _o$ constant, which allows us to integrate equation (4.34) analytically (to obtain a half-parabolic profile), and we integrate that once more with respect to $\xi$ to obtain the vertical flux.

4.3. Cylinder

For the horizontal cylinder (see figure 2a), we need to perform a coordinate transform to polar coordinates, see figure 10. The buoyancy term $\beta$ now depends on the angle $\theta$ with the vertical direction, and the coordinate transform is applied to the governing equations. The melt rate $U$ and film thickness $h$ now depend on the angle $\theta$ instead of height $z$. Since $h \ll R$, with $R$ the radius of the cylinder, the tangential velocity (with $x = r-R$) has the same profile as $w$ for the vertical wall. Equation (4.41) shows the continuity equation in polar coordinates where the axial dependence has been assumed absent, where $u_r$ is the radial velocity and $u_\theta$ is the tangential velocity, with $u_\theta = \bar {\beta }x(2h(\theta )-x)$ analogous to the vertical velocity $w$ for the vertical wall, where $\mu _o$ is taken as constant, conforming to (4.4), and $\bar {\beta } = \beta \sin {\theta }$ is the buoyancy term adapted to the geometry:

(4.41)\begin{equation} \frac{\partial}{\partial x}((R+x)u_r ) = \frac{\partial}{\partial\theta}(u_\theta).\end{equation}

Substituting $u_\theta$ in (4.41) and using an analogous boundary condition to the vertical wall case, $u_r(x=0) = U(\theta )$, results in an expression for $u_r$:

(4.42)\begin{equation} u_r(\theta,x) = U(\theta) + \frac{1}{R}\left(-\frac{1}{3}\frac{{\rm d} }{{\rm d} \theta}(\bar\beta)x^3 +\frac{{\rm d} }{{\rm d} \theta}(\bar\beta h)x^2 \right).\end{equation}

An expression for the melt rate can be obtained by considering mass conservation:

(4.43)\begin{equation} U(\theta) = \frac{1}{R}\frac{\partial}{\partial\theta}\int_0^{h(\theta)}u_\theta(\theta,x)\,{{\rm d} x} = \frac{1}{R}\frac{{\rm d} }{{\rm d} \theta}\left(\frac{2}{3}\bar\beta h^3\right). \end{equation}

The advection–conduction equation, analogous to (4.9), with $\xi = {x}/{h}$, follows

(4.44)\begin{equation} u_r\frac{\partial T}{\partial \xi} = \frac{\alpha}{h}\frac{\partial^2 T}{\partial \xi^2}. \end{equation}

For the temperature we make an approximation analogous to the case of the vertical wall, i.e. $u_r \approx U(\theta )$. The boundary conditions are analogous to the case of a vertical wall:

(4.45)\begin{gather} T|_{\xi=0} = T_o, \end{gather}

(4.46)\begin{gather}T|_{\xi=1} = T_o + 0.75\Delta T, \end{gather}

where $\Delta T = T_\infty - T_o$. The solution for the temperature profile is

(4.47)\begin{equation} T(\xi) = T_o + 0.75\Delta T \frac{\displaystyle \int_0^\xi \exp{ \left[\dfrac{2h}{3R\alpha} \dfrac{{\rm d}}{{\rm d} \theta} ( \bar\beta h^3 ) \xi \right] {\rm d}\xi } }{\displaystyle \int_0^1 \exp{\left[ \dfrac{2h}{3R\alpha} \dfrac{{\rm d}}{{\rm d} \theta} (\bar\beta h^3)\xi \right]{\rm d} \xi}}.\end{equation}

Analogous to (4.13) the heat flux balance is

(4.48)\begin{equation} \rho_oU(\theta)\mathcal{L} = \frac{0.75\lambda \Delta T}{h} \frac{1}{\displaystyle \int_0^1\exp\left[\dfrac{2 \bar\beta h}{3R\alpha}\dfrac{{\rm d}}{{\rm d} \theta}\left(\sin(\theta) h^3\right)\xi \right] {\rm d}\xi}.\end{equation}

We define the quantities

(4.49a,b)\begin{equation} q^4 = \frac{R\alpha}{\beta} \quad \text{and}\quad H = \frac{h^4}{q^4}. \end{equation}

Then the integral in (4.48) becomes

(4.50)\begin{equation} \int_0^1 \exp \left[ \left(\frac{1}{2}\frac{{\rm d} H}{{\rm d} \theta}\sin(\theta) + \frac{2}{3}H \cos(\theta)\right)\xi \right] {\rm d} \xi . \end{equation}

We define a function $G(\theta )$:

(4.51)\begin{equation} G(\theta) = \frac{2}{3}H\cos(\theta) + \frac{1}{2}\frac{{\rm d} H}{{\rm d} \theta}\sin(\theta), \end{equation}

such that we can write the integral as

(4.52)\begin{equation} \int_0^1\exp{[ G(\theta)\xi ]}\, {\rm d} \xi. \end{equation}

From (4.43), (4.49a,b) and (4.51) we deduce that

(4.53)\begin{equation} hU(\theta) = \alpha G(\theta). \end{equation}

Inserting this into (4.48) and using (4.52) results in

(4.54)\begin{equation} 0.75 Ste = \exp(G(\theta)) - 1. \end{equation}

With (4.19) and taking logarithms we arrive at

(4.55)\begin{equation} G(\theta) = \frac{2}{3}H\cos(\theta) + \frac{1}{2}\frac{{\rm d} H}{{\rm d} \theta}\sin(\theta) = \ln{(1+\varLambda)}. \end{equation}

Solving and requesting regularity at $\theta = 0$ gives

(4.56)\begin{equation} H = \left( \frac{h}{q}\right)^4 = 4\ln{(1+\varLambda)}\frac{\displaystyle \int_0^\theta \sin(\theta^*)^{1/3}\,{\rm d} \theta^*}{\sin(\theta)^{4/3}}.\end{equation}

The final solution for the melt film thickness and melt rate can now be found by rewriting equation (4.56), using (4.49a,b), and substituting the result for the melt film thickness in (4.53):

(4.57)\begin{gather} h_{cyl}(\theta) = R^{1/4}\left[ \frac{2\alpha}{\beta}\ln{(1+\varLambda)} \right]^{1/4}f_{cyl}(\theta) , \end{gather}

(4.58)\begin{gather}U_{cyl}(\theta) = \frac{\alpha}{h_{cyl}}\ln{(1+\varLambda)}, \end{gather}

where $f_{cyl}(\theta )$ is a shape function

(4.59)\begin{align} f_{cyl} &= \frac{\left(\displaystyle \int_0^\theta \sin(\theta^*)^{1/3}\, {\rm d} \theta^*\right)^{1/4}}{\sin(\theta)^{1/3}}, \end{align}

(4.60) \begin{align} &= \begin{cases} \frac{\sqrt[4]{3{}_2F_1\left(\dfrac{1}{2},\dfrac{2}{3}; \dfrac{5}{3};\sin(\theta)^2\right)}}{\displaystyle\sqrt{2}} & 0 \leq \theta \leq \dfrac{\rm \pi}{2}\\ \frac{\sqrt[4]{\dfrac{3 \sqrt{3} \varGamma \left(\dfrac{2}{3}\right)^3}{2^{2/3} {\rm \pi}} -\dfrac{3}{4}\sin(\theta)^{{4}/{3}} {}_2F_1\left(\dfrac{1}{2},\dfrac{2}{3};\dfrac{5}{3}; \sin(\theta)^2\right)}}{\displaystyle\sqrt[3]{\sin(\theta)}} & \dfrac{\rm \pi}{2} \leq \theta \leq {\rm \pi}\end{cases}, \end{align}

where ${}_2F_1$ is Gauss's hypergeometric function and $\varGamma$ the complete gamma function. Compared with the expressions that were found for the vertical wall, the difference is in this shape function, which compensates for the geometry varying when following the boundary of the wall, and the dependence on the radius $h_{cyl} \propto R^{1/4}$. Similar expressions occur in Acrivos (1960a,b) and also resemble the solutions for a dissolving vertical cylinder found by Pegler & Davies Wykes (2020).

Figure 10. Schematic overview of molten olive oil (yellow) over a circular shape of frozen olive oil (orange), submerged in water (blue). Here, $U(\theta )$ is the melt rate at the frozen/liquid oil interface. Definitions are used for the horizontal cylinder (polar coordinates) and for the ball (spherical coordinates).

4.4. Ball

The problem of a melting ball is very similar to the cylinder described above. Here, $\theta = 0$ is again defined on the bottom side of the object, the azimuthal angle $\phi$ is defined positive in clockwise direction and $x = r-R$, with $x=0$ at the surface, is defined in the same manner as the horizontal cylinder; see figure 10. An important difference is a flow focusing due to the varying circumference of the ball with changing polar angle $\theta$. The continuity equation in spherical coordinates, where $u_\theta = \bar \beta x(2h-x)$ (again assuming $\mu _o$ is constant), is unchanged and azimuthal symmetry is assumed:

(4.61)\begin{equation} \frac{1}{r^2}\frac{\partial}{\partial r}(r^2u_r) +\frac{1}{r\sin(\theta)}\frac{\partial}{\partial\theta}(u_\theta\sin(\theta)) = 0. \end{equation}

From this, $u_r$ can be found as before:

(4.62)\begin{equation} u_r = U(\theta) - \bar \beta \left(2 h\cos(\theta) + \frac{{\rm d} h}{{\rm d} \theta}\sin(\theta)\right)\frac{x^2}{R} + \frac{2}{3}\bar\beta \cos(\theta)\frac{x^3}{R}. \end{equation}

The control volume over the film thickness is now taken in three dimensions:

(4.63)\begin{equation} 2{\rm \pi} R^2\sin(\theta) U(\theta) = \frac{{\rm d} }{{\rm d} \theta}\int_0^h2{\rm \pi} R\sin(\theta)u_\theta \,{{\rm d} x}, \end{equation}

which after substitution of $u_\theta$ gives

(4.64)\begin{equation} U(\theta) = \frac{1}{R\sin(\theta)}\frac{\partial}{\partial \theta}\left(\frac{2}{3}\bar\beta h^3\sin(\theta)^2\right). \end{equation}

Following a similar procedure as before, we find a new function $G(\theta )$:

(4.65)\begin{equation} G(\theta) = \frac{4}{3}\cos(\theta) H + \frac{1}{2}\sin(\theta)\frac{{\rm d} H}{{\rm d} \theta} = \ln{(1+\varLambda)}. \end{equation}

This can be solved to obtain the final solutions for the film thickness and the melt rate:

(4.66)\begin{gather} h_{ball}(\theta) = R^{1/4}\left[ \frac{2\alpha}{\bar\beta}\ln{(1 + \varLambda)} \right]^{1/4}f_{ball}, \end{gather}

(4.67)\begin{gather}U_{ball}(\theta) = \ln{(1 + \varLambda)}\frac{\alpha}{h_{ball}}, \end{gather}

where $f_{ball}$ is a shape function for the spherical geometry, different from the one for the cylinder,

(4.68)\begin{align} f_{ball} &= \frac{\displaystyle \left(\int_0^\theta \sin(\theta^*)^{5/3}\, {\rm d} \theta^*\right)^{1/4}}{\displaystyle\sin(\theta)^{2/3}}, \end{align}

(4.69) \begin{align} &= \begin{cases} \frac{\sqrt[4]{3 {}_2F_1\left(\dfrac{1}{2},\dfrac{4}{3}; \dfrac{7}{3};\sin(\theta)^2\right)}}{\displaystyle 2^{3/4}} & 0 \leq \theta \leq \dfrac{\rm \pi}{2}\\ \frac{\sqrt[4]{\dfrac{\sqrt{3} \varGamma \left(\dfrac{1}{3}\right)^3}{5 \sqrt[3]{2} {\rm \pi}} -\dfrac{3}{8}\sin(\theta)^{{8}/{3}} {}_2F_1\left(\dfrac{1}{2},\dfrac{4}{3};\dfrac{7}{3};\sin(\theta)^2\right) }}{\displaystyle\sin(\theta)^{{2}/{3}}} & \dfrac{\rm \pi}{2} \leq \theta \leq {\rm \pi}\end{cases}. \end{align}

We would like to highlight that the solutions for the thickness $h$ ((4.18), (4.57) and (4.66)) and for the melt rate $U$ ((4.20), (4.58) and (4.67)) for all three geometries have a very similar form.

5. Discussion

In § 4 we already discussed some approximations that we made. We have shown that our models match relatively well with our experimental findings, especially considering that our model does not contain any free (fitting) parameters. During the derivation we made several approximations, and the model does not include all effects. We will now go through the remaining approximations and assess their validity.

First, we have made use of the thin film approximation, which seems like a reasonable approximation since our layer thickness is of $O$(mm) while our objects are of $O(100\ {\rm mm})$.

Second, we had made the assumption that ${\partial w}/{\partial x} |_{x=h} \approx 0$, which also seems justified since the ratio of the viscosities, even in the worst case, is $\mu _o/\mu _w \approx 75 \gg 1$.

Third, the assumption of a constant contact temperature $T^*$ along the wall turns out to be realistic for the vertical wall and for the cylinder, as can be seen from figures 3 and 4. In the case of the ball, the agreement with the model is good between $60^\circ$ and $130^\circ$ but there is a significant difference in the bottom region; see figure 5. The reason for that becomes clear from the following analysis. When two semi-infinite media of different temperatures are brought into contact the interface assumes a temperature, the contact temperature $T^*$, given in (4.10), which remains constant thereafter. In our case one of the media, the melt film, is of finite extent $h$. If we take, for convenience, equal material properties at both sides, the contact temperature changes in time according to (see Appendix A)

(5.1)\begin{equation} T^* = T_o + \frac{T_\infty - T_o}2 \operatorname{erf}\left(\frac h{\sqrt{\alpha t}}\right),\end{equation}

where $\operatorname {erf}(\tau ) = (2/{\sqrt {{\rm \pi} }}) \int _0^\tau {\rm e}^{-t^2}\, {\rm d} t$ is the error function. Given a typical time of 20 s of contact of a water element from top to bottom of the ball, and a $h = 2\ {\rm mm}$, with $\alpha = 8 \times 10^{-8}\ {\rm m}^2\ {\rm s}^{-1}$, this means that, after 20 s, the error function in (5.1) has still 96 % of its initial value. However, near the bottom, where $\theta =0$ and $h$ is very small, of the order of tenths of a millimetre, and where probably the contact time is longer, the error function decreases from its initial value. This means a drop in the contact temperature and thereby of the melting rate. For the vertical wall the film thickness is of the order of a $1\ {\rm mm}$ along most of the wall. Unfortunately, we are unable to locally measure the temperature profiles since the scales are too small (and probes too big).

Four, throughout the analysis we had assumed that the problem is time-independent. Since our freezing temperature is $T_i = -14\,^\circ$C and the melting point is $T_o = -8\,^\circ$C all the matter has to warm up 6 K before it melts. At $t=0$ a skin layer of the temperature grows inside the material. The typical dimensionless similarity variable $\eta = x/\sqrt {\alpha t}$ can then be used to find the temperature profile inside the material which goes like $\operatorname {erf}(\eta /2)$. The typical penetration depth of the temperature is thus given by $\sqrt {\alpha t}$, such that the speed at which this front moves is $U_{{skin\ layer}} = \sqrt {{\alpha }/{(4t)}}$. If we equate this to our melt speed of $U \approx 10\ \mathrm {\mu }{\rm m}\ {\rm s}^{-1}$ we find a typical time scale of $t \approx 200$ s at which the speed of the penetrating skin layer reaches the same speed as the melting boundary. In other words, for times below a few minutes, there is energy spent on heating up material that is not melted in this time. For larger times, the speeds of the skin layer and the melting boundary are the same, and energy is only spent on heating material that is also melted. This thus means that for small times we overpredict the melting rate, and the actual melting rate is slightly lower since we heat more material up than we melt. Note that the energy spent on heating is relatively small as compared with the energy spent on the phase transition $c_p(T_i - T_o)/\mathcal {L} = 4.4\,\%$, such that the effect is comparatively small. This is, however, a clear case in which our prediction causes an over-prediction of the melt rate. A small correction could be made by introducing an effective latent heat. The idea is to include not only the energy for the phase change (the latent heat), but also the energy needed to heat up the material from its initial temperature of $-14\,^\circ$C to the melting point $-8\,^\circ$C. Heating up the oil by $\Delta T = -8\,^\circ {\rm C} - (-14\,^\circ {\rm C}) = 6\ {\rm K}$ prior to melting requires an energy of $c_p\Delta T = 11.82\ {\rm kJ}\ {\rm kg}^{-1}$. Comparing this with the latent heat of olive oil $267\ {\rm kJ} {\rm kg}^{-1}$, this amounts to roughly 4.4 %. An effective latent heat could then be defined that takes into consideration: $\mathcal {L}_{eff} = \mathcal {L}+c_p(T_i - T_o)$.

Another experimental issue not yet discussed in detail might be that, as we remove the frozen oil from the metallic mould and then place it in our water tank the oil has slightly heated up. We are not sure whether all objects were $-14\,^\circ$C throughout. Whereas the melting profiles $U$ are more or less constant for the vertical wall and the large horizontal cylinder, for the ball and the small cylinder the melting profiles $U$ change a bit over time. The reason between those could be the varying time between releasing the olive oil from the mould and placing it in our aquarium; see figure 1. We hypothesise that, for the vertical wall and the large cylinder, the object was left (relatively) long in air and would already have started forming the temperature skin layer. For the experiments with the small cylinder and the ball we quickly used the object after releasing it from the mould, not allowing the temperature skin layer to develop. This could then explain the differences in steadiness between the melting profiles’ $U$ for the various geometries.

Fifth, up to now we have ignored surface tension effects. To justify this assumption, we compare the hydrostatic pressure and the Laplace pressure. The typical difference in hydrostatic pressures (for the vertical wall at mid-height) is given by $P_{hydrostatic} = g\Delta \rho L/2 \approx 200\ {\rm Pa}$. Now the Laplace pressure can be approximated as $P_{Laplace} = \sigma |{\rm d}^2h/{\rm d} z^2 | = 0.2\ {\rm mPa}$. Since the Laplace pressure is smaller by 6 orders of magnitude, we are confident in ignoring the effects of surface tension. This ratio is substantially less for the cylindrical and spherical experiments; ${\approx }200$ and ${\approx }500$, respectively, but still high enough for us to be confident in ignoring the effects of surface tension. However, surface tension might be important for the spherical and cylindrical geometries for $\theta > 90^\circ$ to stabilise the oil film.

Sixth, we observe no turbulence in the ambient fluid nor in the thin oil film, such that we expect our results to be valid at least up to ${Ra}=10^8$. Since our Stefan number is relatively low (0.21) the heating up of the oil is small as compared with the heat needed to melt. For high Stefan numbers we expect an additional conduction problem inside the objects that, for small $t$, delays the melting. So we can therefore conclude that our results should also be valid for Stefan numbers at least up to ${Ste}=0.2$.

Lastly, throughout our derivation we have considered the ambient water to be stationary, which we have discussed before. The velocity in the oil layer along the wall is of order 1 mm s$^{-1}$ and the velocity in water 1 cm s$^{-1}$. This is so low that heat transfer is still mainly by conduction and will influence $T^*$ by a small amount. However, for applications with much larger Rayleigh and Stefan numbers this situation might change, and our analysis and assumptions will no longer hold. For example increased velocities in the surrounding water may affect the pinch-off at the apex due to increased shear forces, or the attached oil film may show instabilities which might result in detachment of the oil film from the object prior to reaching the apex.