3.1. For $\textit{Re}_0 \lt 212$ : axisymmetric melting regime

In the absence of melting, the wake behind a fixed sphere remains axisymmetric up to the first bifurcation at $\textit{Re}_{\textit{c1}}\approx 212$ (Johnson & Patel Reference Johnson and Patel1999; Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012). For a melting sphere with $\textit{Re}_0\lt 212$ , both the interface and the surrounding flow are therefore expected to preserve axisymmetry throughout the evolution. This is confirmed in figure 4(a) for $\textit{Re}_0=180$ , which displays the temperature field, mid-plane streamlines and the 3-D interface coloured by $v_{{\varGamma }}$ . At all times shown, i.e. $t/t_{\!f} = 0.010$ , $0.412$ , and $0.722$ with $t_{\!f}$ denoting the moment at which the solid vanishes, the wake exhibits a symmetric recirculation zone and the interfacial shape remains axisymmetric. On the upstream side, $v_{{\varGamma }}$ decreases monotonically with the polar angle $\varphi$ , as quantified in figure 4(b). This trend arises from boundary-layer development: starting from the stagnation point, the thermal boundary layer thickens with increasing arc length, thereby reducing the local heat flux. Since $v_{{\varGamma }} \propto \delta _\theta ^{-1}$ (equation (2.5)), the melting rate correspondingly decays with $\varphi$ . By contrast, on the downstream side, flow separation isolates the rear surface from direct contact with the warmer inflow. The separated recirculation region traps relatively cold fluid, suppressing convective transport into the wake. As a result, the temperature gradient across the rear interface remains weak and $v_{{\varGamma }}$ exhibits an abrupt decline near the separation point, followed by persistently small values further downstream, as evident in figure 4(b).

Figure 4.

Melting dynamics at $\textit{Re}_0=180$ . (a) Snapshots of the temperature field $\theta$ and streamlines on the mid-plane at $z=0$ , together with the 3-D melting interface coloured by the local melting rate $v_{{\varGamma }}$ , shown at $t/t_{\!f} = 0.010,0.412,0.722$ (top to bottom). The corresponding effective Reynolds numbers $\textit{Re}_e(t)$ are indicated in each panel, and the same notation is used in the subsequent figures illustrating the melting process. (b) Angular distribution of $v_{{\varGamma }}$ along the interface as a function of the polar angle $\varphi$ , from $t/t_{\!f}=0.052$ to $0.845$ in increments of $\Delta t/t_{\!f}=0.072$ . The inset at the lower left shows the polar angle $\varphi$ ; the origin is the body’s instantaneous mass centre $(x_c,0,0)$ , and $\varphi$ is measured from the front stagnation point. The inset at the upper right shows the interface-averaged melting rate $\bar {v}_{\varGamma }(t)$ as a function of the normalised remaining volume $V(t)/V_0$ , revealing a clear $-1/6$ scaling.

As shown in figure 4(a), the front of the melting sphere approximately retains a spherical arc throughout the process, whereas the rear interface progressively flattens. This behaviour is consistent with previous experimental and numerical observations (Huang et al. Reference Huang, Jinzi, M.Nicholas and Ristroph2015; Yang et al. Reference Yang, Howland, Liu, Verzicco and Lohse2024a ; Hao & Tao Reference Hao and Tao2001), underscoring the close coupling between wake structure and rear-interface evolution. The flattening originates from the recirculating wake, which establishes a secondary stagnation point at the rear centre ( $\varphi = \pi$ ). At this location, the recirculating wake produces a local stagnation point where the thermal boundary layer reaches its minimum thickness, yielding the peak melting rate $v_{\varGamma }$ . This enhanced recession at $\varphi =\pi$ diminishes the rear curvature and drives the interface towards a planar shape. According to classical boundary-layer theory for axisymmetric stagnation-point flow (Schlichting & Gersten Reference Schlichting and Gersten2000, (5.70)), a planar rear interface leads to an almost uniform viscous boundary-layer thickness. At $\textit{Pr}=7$ , the thermal layer lies within the viscous one and obeys the same scaling, so its thickness is likewise nearly uniform. This, in turn, implies an approximately uniform melting rate along the entire rear surface. Hence, once the interface becomes nearly flat, the boundary-layer structure naturally enforces the spatial uniformity of $v_{\varGamma }$ at the rear interface observed later in figure 4(b) at $t/t_{\!f}=0.845$ , stabilising the planar configuration. This tendency mirrors the self-similar retreating shapes documented for dissolving or eroding bodies in high-Reynolds-number flows ( $\textit{Re}\sim 10^4$ ), where nearly uniform wall shear or solute concentration profiles control the late-stage morphology (Ristroph et al. Reference Ristroph, Moore, Childress, Shelley and Zhang2012; Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013; Huang et al. Reference Huang, Jinzi, M.Nicholas and Ristroph2015). The present results thus suggest that similar boundary-layer mechanisms operate at the rear interface in the moderate- $\textit{Re}_0$ regime. Notably, as shown later in § 3.4, analogous uniformity also emerges at the front interface for high Reynolds number of $\textit{Re}_0=1000$ .

In addition, figure 4(b) shows that the overall melting rate increases over time. This trend is further illustrated by the upper-right inset, which presents the interface-averaged melting rate as a function of the normalised remaining volume. The observed behaviour can be rationalised using classical boundary-layer scaling arguments (Johnson & Patel Reference Johnson and Patel1999; Schlichting & Gersten Reference Schlichting and Gersten2000; Huang et al. Reference Huang, Jinzi, M.Nicholas and Ristroph2015). For a shrinking body, the viscous boundary-layer thickness satisfies the scaling $\delta _u / V^{1/3} \sim (V^{1/3})^{-1/2}$ . Consequently, the interface-averaged melting rate scales as

(3.1) \begin{equation} \bar {v}_{{\varGamma }} \sim \delta _\theta ^{-1} \sim \delta _u^{-1} \sim V^{-1/6}, \end{equation}

which agrees with the $-1/6$ scaling observed in the inset at the upper right of figure 4(b).

The influence of the initial Reynolds number $\textit{Re}_0$ on melting dynamics within the axisymmetric regime is illustrated in figure 5. For non-melting spheres, increasing $\textit{Re}_0$ reduces the separation angle $\varphi _s$ and enlarges the toroidal wake vortex (Johnson & Patel Reference Johnson and Patel1999). A comparable trend emerges in the melting problem: larger $\textit{Re}_0$ produces stronger wake recirculation, which progressively alters the rear-interface morphology. Figure 5(a–d) shows the mid-plane interface ( $z=0$ ) at successive times for different $\textit{Re}_0$ , with the local melting rate $v_{{\varGamma }}$ indicated by colour. To enable direct comparison, figure 5(e) superposes the interfaces at the moment when the solid volume has decreased to $V(t)/V_0=0.1$ . At moderate Reynolds numbers ( $\textit{Re}_0=25$ and $50$ ), the bodies preserve an approximately ellipsoidal form throughout melting, and the rear interface remains convex. In contrast, at higher $\textit{Re}_0$ ( $100 \lesssim \textit{Re}_0 \leqslant 200$ ), the rear interface undergoes marked flattening, induced by stronger and more persistent recirculating wakes. This contrast stems from the reduced extent of separation at low-to-moderate $\textit{Re}_0$ , where the separation angle remains large ( $\varphi _s \approx 140^\circ$ for $\textit{Re}_0 \leqslant 50$ ; Johnson & Patel Reference Johnson and Patel1999), and from the decrease in effective Reynolds number $\textit{Re}_e$ as the body shrinks. A more quantitative description of these shape changes is provided in § 4, where aspect-ratio-based measures are introduced. The implications of interface morphology for global melting rates and hydrodynamic forces are analysed in §§ 4 and 5.

Figure 5.

Effect of the initial Reynolds number $\textit{Re}_0$ on interface evolution in the axisymmetric melting regime. (a–d) Time evolution of the interface on the mid-plane at $z=0$ for $\textit{Re}_0=25,\ 100,\ 150$ and $200$ , coloured by the local melting rate $v_{{\varGamma }}$ . The sequences span nearly the entire melting process, with frames shown at uniform time intervals up to $t/t_{\!f} \approx 0.9$ . (e) Comparison of the interfaces on the mid-plane at $z=0$ when the remaining solid volume reaches $V(t)/V_0=0.1$ , highlighting the progressive flattening of the rear surface as $\textit{Re}_0$ increases.

3.2. For $212 \lt \textit{Re}_0 \lt 273$ : steady planar-symmetric melting regime

For $\textit{Re}_0 \gt 212$ , the wake behind a translating sphere loses axisymmetry but retains a reflectional symmetry plane, up to the second bifurcation threshold at $\textit{Re}_{\textit{c2}} \approx 273$ (Johnson & Patel Reference Johnson and Patel1999; Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012). Within this steady planar-symmetric regime, we examine the case $\textit{Re}_0 = 270$ . Since the flow spontaneously selects a symmetry plane, rather than being imposed, we analyse it using a rotated coordinate system $x$ – $\eta$ – $ \zeta$ , where $x$ – $\eta$ defines the symmetry plane and $\zeta$ is the transverse direction. Note that, in the present study, the symmetry plane is entirely determined by numerical noise. As a result, the rotated coordinate system $x$ – $\eta$ – $\zeta$ used for the symmetry-breaking regimes throughout the paper is $\textit{Re}_0$ –dependent. Figure 6(a) presents the 3-D streamlines in the vicinity of the sphere at $t/t_{\!f} = 0.502$ , viewed along the $\zeta$ direction. The upstream fluid is first entrained by the lower spiral, then transported azimuthally towards the upper spiral and finally guided around the lower spiral before moving downstream. Figure 6(b) shows the 3-D isosurfaces of the streamwise vorticity $\omega _x$ at the same instant, with red and blue surfaces representing $\omega _x = \pm 0.25$ . These illustrate the canonical features of the planar-symmetric wake: the recirculating vortex ring behind the sphere becomes tilted, with unequal upper and lower spirals, as illustrated in figure 6(c) where the snapshots of the melting body and wake are presented. The stronger spiral entrains fluid into the weaker one, thereby breaking the closed toroidal bubble of the axisymmetric melting regime. This imbalance generates a pair of counter-rotating streamwise vortices, which induce a transverse lift force in the $\eta$ direction. Still in figure 6(c), the front interface remains rounded, while the rear evolves asymmetrically: the lower side melts faster, producing an inclined planar surface that ultimately stabilises at an angle of $-9.2^\circ$ relative to the $\eta$ direction. This asymmetry arises because the upper spiral traps colder recirculated fluid, suppressing local melting, whereas the lower spiral is replenished by warmer inflow, enhancing melting. The rear stagnation point (red markers) shifts laterally during the evolution, initially located at negative $\eta$ ( $\eta \lt 0$ or $\varphi \gt \pi$ ) but crossing into positive $\eta$ ( $\eta \gt 0$ or $\varphi \lt \pi$ ) at late times of $t/t_{\!f} = 0.920$ , accompanied by a reversal of spiral dominance, as highlighted in the magnified view. This reversal can be attributed to the inclination of the melting body: the streamline path connecting the front stagnation point to the rear is shorter in the lower half, resulting in stronger vortex formation on that side. Such a scenario is also consistent with vortex structures observed behind inclined non-melting bodies such as spheroids (Wang et al. Reference Wang, Yang, Andersson, Zhu, Liu, Wang and Lu2021b ) and cylinders (Kharrouba, Pierson & Magnaudet Reference Kharrouba, Pierson and Magnaudet2021). The implications of this spiral reversal for the lift force are examined in § 5.2.

Figure 6.

Melting dynamics at $\textit{Re}_0=270$ in the steady planar-symmetric regime. (a) Three-dimensional streamlines viewed along the $\zeta$ direction. (b) Three-dimensional isosurfaces of the streamwise vorticity $\omega _x$ viewed along the $\eta$ direction, with red and blue surfaces indicating $\omega _x=\pm 0.25$ . The results shown in (a,b) correspond to $t/t_{\!f} = 0.502$ , at which the effective Reynolds number is $\textit{Re}_e(t)\approx 199.56$ . (c) Snapshots of the temperature field $\theta$ , streamlines on the symmetry plane ( $x$ – $\eta$ ) and the 3-D interface coloured by the local melting rate $v_{{\varGamma }}$ , shown at $t/t_{\!f}=0.008$ , $0.335$ , $0.669$ and $0.920$ . Red points mark the rear stagnation position, highlighting its lateral migration along the $\eta$ direction and eventual reversal of spiral dominance in the wake.

Figure 7 examines in detail the melting characteristics at $t/t_{\!f}=0.335$ for $\textit{Re}_0=270$ . Figure 7(a) displays the temperature field $\theta$ and the local melting rate $v_{{\varGamma }}$ in the reflectional symmetry plane ( $x$ – $\eta$ ). A pronounced thermal asymmetry is observed: the lower spiral contains warmer fluid than the upper, consistent with the enhanced melting on the lower rear surface. The strongest melting occurs at the rear stagnation point, where the temperature gradient is steepest. This ‘melting kernel’ is highlighted in all panels by red circles. The rear-view distribution of $v_{{\varGamma }}$ , shown in figure 7(b), further emphasises the asymmetry, with the lower hemisphere retreating more rapidly. Figure 7(c) shows the temporal evolution of $v_{{\varGamma }}$ as a function of the polar angle $\varphi$ (defined in the inset of figure 7 a), from $t/t_{\!f}=0.084$ to $0.837$ at intervals of $\Delta t/t_{\!f}=0.084$ . The peak melting rate consistently coincides with the stagnation point, whose angular location gradually migrates towards larger $\eta$ with time, reflecting the evolving wake geometry we described in figure 6(b). In addition, the spatial variation of $v_{{\varGamma }}$ across the rear interface becomes progressively smoother, reminiscent of the trend towards uniform boundary-layer-driven melting observed in the axisymmetric melting regime. This suggests that, even in the presence of reflectional asymmetry, the system dynamically reorganises to reduce rear-interface thermal gradients.

Figure 7.

Rear-interface melting characteristics at $\textit{Re}_0=270$ . (a) Temperature field $\theta$ and local melting rate $v_{{\varGamma }}$ on the symmetry plane ( $x$ – $\eta$ ). (b) Distribution of $v_{{\varGamma }}$ viewed from the rear (positive $x$ axis). The results shown in (a,b) correspond to $t/t_{\!f} = 0.335$ , at which the effective Reynolds number is $\textit{Re}_e(t)\approx 220.66$ . (c) Temporal evolution of $v_{{\varGamma }}$ as a function of the polar angle $\varphi$ , from $t/t_{\!f}=0.084$ to $0.837$ at constant time intervals. Red circles mark the rear stagnation point in all panels.

Figure 8 summarises the influence of the initial Reynolds number $\textit{Re}_0$ on the morphological evolution of the rear interface within the steady planar-symmetric melting regime, considering cases with $\textit{Re}_0 = 240$ , $250$ , $260$ and $270$ . As $\textit{Re}_0$ increases, the wake asymmetry intensifies, producing progressively larger inclination angles of the rear interface relative to the $x$ direction. This trend is illustrated in figure 8(a), which overlays the interfaces at a fixed volume fraction $V(t)/V_0=0.1$ . While the front hemispheres collapse onto nearly identical spherical arcs, the diverging rear interfaces clearly demonstrate the increasing anisotropy of wake-induced melting at higher Reynolds numbers. To quantify this inclination, we define the rear-interface angle as $\textrm{tan}^{-1}(\textrm{sgn}(\bar {n}_\eta )\sqrt {\bar {n}^2_\eta +\bar {n}^2_\zeta }/|\bar {n}_x|)$ , where $(\bar {n}_x,\bar {n}_\eta ,\bar {n}_\zeta )$ are the components of the area-averaged outward normal vector $\bar {\boldsymbol{n}}$ on the rear interface. The time evolution of the rear-interface angle, shown in figure 8(b), confirms this trend: the inclination increases from about $-6.2^\circ$ at $\textit{Re}_0=240$ to $-9.2^\circ$ at $\textit{Re}_0=270$ . At late stages ( $t/t_{\!f} \gt 0.8$ ), the inclination saturates, consistent with the nearly uniform rear melting rates observed earlier (see figure 7 c).

Figure 8.

Influence of $\textit{Re}_0$ on rear-interface evolution in the steady planar-symmetric melting regime. (a) Comparison of interfaces on the symmetry plane ( $x$ – $\eta$ ) when the remaining solid volume reaches $V(t)/V_0=0.1$ . (b) Time evolution of the inclination angle of the rear interface relative to the $x$ axis for different $\textit{Re}_0$ . Increasing $\textit{Re}_0$ enhances the inclination of the rear interface.

3.3. For $273 \lt \textit{Re}_0 \lt 355$ : periodic planar-symmetric melting regime

When the Reynolds number surpasses the second critical threshold, $\textit{Re}_{\textit{c2}} \approx 273$ , the wake undergoes a Hopf bifurcation, transitioning from steady to unsteady behaviour. The resulting flow exhibits time-periodic vortex shedding while retaining planar reflectional symmetry (Johnson & Patel Reference Johnson and Patel1999). In this regime, the wake generates oscillatory lift forces with a non-zero mean, yet the vortex separation line remains nearly fixed and the rear stagnation point fluctuates only weakly within $\pm 9^\circ$ (Johnson & Patel Reference Johnson and Patel1999). Thus, despite the periodic shedding, the rear stagnation point is effectively quasi-stationary.

Accordingly, we anticipate that once melting commences, the rear interface will continue to evolve into a planar and inclined surface, in close analogy with the behaviour observed in the steady planar-symmetric melting regime (§ 3.2). This expectation is confirmed in figure 9(a), which presents representative snapshots of the streamline topology and temperature field $\theta$ on the symmetry plane ( $x$ – $\eta$ ) for $\textit{Re}_0 = 300$ , together with the 3-D interface coloured by the local melting rate $v_{{\varGamma }}$ , over the interval $0 \lt t/t_{\!f} \lt 0.910$ . Once melting is activated, the rear interface begins to incline, and the wake gradually loses its temporal periodicity, reorganising into a new steady asymmetric configuration. By $t/t_{\!f} = 0.910$ , the lower spiral has become larger than the upper, a reversal analogous to that observed in the steady planar-symmetric melting regime (§ 3.2). The temporal development of this melting-induced transition is further illustrated in figure 9(b), which shows that a tilted and nearly planar rear interface persists in the later stage. The dynamical suppression of wake periodicity is quantified in figure 9(c), where the time histories of the polar angle of the rear stagnation point are plotted. These quantities initially exhibit small oscillations, consistent with the underlying unsteady wake. As melting progresses, the effective Reynolds number $\textit{Re}_e(t)$ decreases, leading to a decay of these oscillations and a monotonic shift towards smaller $\varphi$ values, corresponding to the positive $\eta$ direction. This behaviour demonstrates that interfacial evolution suppresses vortex shedding and re-establishes a quasi-steady wake morphology.

Figure 9.

Melting processes of $\textit{Re}_0=300$ in the periodic planar-symmetric melting regime. (a) Snapshots of the temperature field $\theta$ and streamline at the symmetry plane ( $x$ – $\eta$ ) and 3-D interface coloured by the melting rate $v_{{\varGamma }}$ at $t=1,40,80,115$ ( $t/t_{\!f}=0.008,0.317,0.633,0.910$ ), from top to bottom. The red points indicate the rear stagnation points at each moment. (b) Time evolution of the interface coloured by the melting rate $v_{{\varGamma }}$ at the reflectional symmetry plane ( $x$ – $\eta$ ) from $t=0$ to $t=110$ ( $t/t_{\!f}=0.871$ ) with a constant time interval of $\Delta t=10$ ( $\Delta t/t_{\!f}=0.079$ ). (c) Time evolution of polar angle of the rear stagnation point, where the dynamical suppression of wake periodicity is observed as the effective Reynolds number $\textit{Re}_e(t)$ decreases.

Figure 10 presents a phase-resolved view of the melting dynamics during a single oscillation cycle in the early stage of melting. Four characteristic phases are selected within one shedding period, $\psi = 0$ , $\pi /2$ , $\pi$ and $3\pi /2$ , corresponding to $t = 64.1$ , $65.5$ , $66.9$ and $68.3$ ( $t/t_{\!f} = 0.507, 0.519, 0.530, 0.540$ ). Figure 10(a) shows the temperature field and streamlines on the symmetry plane ( $x$ – $\eta$ ), together with the interface coloured by $v_{{\varGamma }}$ . While the global wake remains quasi-steady, small periodic modulations are evident, reflecting the influence of vortex shedding on local convective transport. Figure 10(b) illustrates the rear-face distribution of $v_{{\varGamma }}$ , which exhibits clear phase-dependent variations: the peak melting rate oscillates in synchronisation with the shedding cycle. This behaviour arises from periodic changes in the thickness of the thermal boundary layer at the rear stagnation zone, which governs local heat flux. Figure 10(c) confirms this interpretation by showing the polar distribution of $v_{{\varGamma }}$ on the symmetry plane ( $x$ – $\eta$ ). At $\psi = 0$ ( $t=64.1$ ), the lower spiral is in the process of detachment, which weakens convective transport near the rear surface and reduces the melting rate. In contrast, at $\psi = \pi$ ( $t=66.9$ ), the lower spiral is replenished by inflow and remains attached to the surface, intensifying local convection and enhancing the melting rate.

Figure 10.

Phase-resolved melting dynamics during one oscillation cycle at $\textit{Re}_0=300$ . Four phases are shown, $\psi = 0, \pi /2, \pi , 3\pi /2$ , corresponding to $t = 64.1, 65.5, 66.9, 68.3$ ( $t/t_{\!f}=0.507, 0.519, 0.530, 0.540$ ). During this interval, the effective Reynolds number $\textit{Re}_e(t)$ decreases from about $194.36$ to $189.02$ . (a) Temperature field $\theta$ and streamlines on the symmetry plane ( $x$ – $\eta$ ), together with the 3-D interface coloured by the local melting rate $v_{{\varGamma }}$ . (b) Distribution of $v_{{\varGamma }}$ on the rear face of the solid, viewed from the rear. (c) Polar distribution of $v_{{\varGamma }}$ at the symmetry plane ( $x$ – $\eta$ ) for the four phases, where the inset defines the polar angle $\varphi$ for this regime. Red points denote the rear stagnation point in all panels.

3.4. For $\textit{Re}_0 \gt 355$ : chaotic melting regime

When the initial Reynolds number $\textit{Re}_0$ exceeds the third critical threshold, $\textit{Re}_{\textit{c3}} \approx 355$ , the wake behind a non-melting sphere transitions into a fully 3-D chaotic state (Mittal Reference Mittal1999). For a melting sphere in this regime, the rear stagnation point, typically associated with the maximum local melting rate, hence no longer remains spatially or temporally fixed, but undergoes irregular fluctuations. Consequently, the rear interface cannot evolve into a stable planar surface, in contrast with the behaviours observed in the steady and periodic planar-symmetric melting regimes.

Figures 11 and 12 confirm these expectations. For both $\textit{Re}_0=500$ and $1000$ , figures 11(a) and 12(a) show mid-plane temperature and streamline fields superimposed on the 3-D interface coloured by $v_{{\varGamma }}$ , while figures 11(b) and 12(b) display the corresponding rear-face melting rate distributions. Here we adopt the Cartesian $x$ – $y$ – $z$ coordinates rather than $x$ – $\eta$ – $\zeta$ , since no clear planar symmetry persists. Unlike at $\textit{Re}_0=300$ (figure 9), the wakes now exhibit strongly disordered vortex shedding and intense small-scale turbulence, characteristic of chaotic flow. At $\textit{Re}_0=500$ , the rear surface remains rounded only up to $t/t_{\!f} \leqslant 0.387$ , whereas at $\textit{Re}_0=1000$ a rounded rear persists throughout. This difference reflects the role of small-scale vortices generated by Kelvin–Helmholtz-type instabilities at boundary-layer separation (Tomboulides & Orszag Reference Tomboulides and Orszag2000). These vortices disrupt large-scale coherence, suppress organised recirculation and render the rear stagnation point highly unsteady in both space and time. Consequently, the peak in $v_{{\varGamma }}$ is spatially diffuse, and the interface never reorganises into a planar rear – in sharp contrast to the lower- $\textit{Re}_0$ regimes (§§ 3.1–3.3), where a relatively stable stagnation point drives the formation of a well-defined planar interface.

Figure 11.

Melting dynamics for $\textit{Re}_0 = 500$ . (a) Snapshots of the temperature field $\theta$ and streamlines on the mid-plane at $z=0$ , together with the 3-D interface coloured by the local melting rate $v_{{\varGamma }}$ , at $t=1,62,123$ ( $t/t_{\!f}=0.006,0.387,0.767$ ). (b) Corresponding melting rate distributions on the rear surface.

Figure 12.

As in figure 11, but for $\textit{Re}_0=1000$ . Snapshots are shown at $t=1,90,180$ ( $t/t_{\!f}=0.004,0.388,0.776$ ).

An exception to this chaotic pattern is observed at $t/t_{\!f} = 0.767$ for the case $\textit{Re}_0 = 500$ (figure 11 a), where a weakly inclined rear interface emerges, oriented towards the negative $y$ axis. This phenomenon can be attributed to the progressive reduction of the effective Reynolds number $\textit{Re}_e(t)$ . As $\textit{Re}_e(t)$ falls below the chaotic threshold, the wake transitions back into a periodic, and eventually steady regime, thereby restoring partial symmetry and enabling planar interface development. Despite the irregularity of the wake in the chaotic regime, figures 11(b) and 12(b) reveal a gradual increase in the spatial uniformity of $v_{{\varGamma }}$ across the rear surface as melting progresses. This observation reinforces the general principle, established at lower Reynolds numbers, that the system dynamically reorganises to reduce gradients in the melting rate – a trend that persists even in high- $\textit{Re}_0$ turbulent environments.

Figure 13 provides a quantitative view of this tendency for $\textit{Re}_0=1000$ , with figures 13(a) and 13(b) showcasing the evolution of the melting rate along the front ( $0\lt \varphi \lt \pi /2$ ) and the back ( $\pi /2\lt \varphi \lt \pi$ ) interfaces, respectively. Here $\bar {\bar {v}}_{{\varGamma }}(t)$ is evaluated by averaging the local $v_{{\varGamma }}(t)$ over a narrow angular band ( $\varphi \pm 1.5^\circ$ ) and a time window ( $t \pm 10$ ). We also note that the results are nearly identical for the opposite side ( $-\pi \lt \varphi \lt 0$ ) due to the uniformity of the melting in this chaotic flow regime, and hence they are not shown here. The resulting distributions in both figures confirm that, even amid turbulence, the interface evolution systematically promotes homogenisation of $v_{{\varGamma }}$ over time. More specifically, the chaotic wake dynamics does not preclude self-organisation of melting: instead, the melting rate evolves towards greater spatial uniformity, with homogenisation extending from the rear to the front interface at $\textit{Re}_0=1000$ – a feature absent at lower Reynolds numbers ( $\textit{Re}_0 \leqslant 300$ ) but consistent with the dissolution process in the high- $Re_0$ regime reported by Huang et al. (Reference Huang, Jinzi, M.Nicholas and Ristroph2015). Besides confirming that a melting body tends to converge towards a shape with nearly uniform material removal along its surface, it is further demonstrated (Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013; Moore Reference Moore2017) that this principle inherently promotes the formation of a smooth, rounded front interface – a prediction that is fully consistent with our numerical observations.

Figure 13.

Local melting rate distribution for the case $\textit{Re}_0=1000$ from $t=10$ ( $t/t_{\!f}=0.043$ ) to $t=210$ ( $t/t_{\!f}=0.905$ ) with a constant time interval $\Delta t=20$ ( $\Delta t/t_{\!f}=0.086$ ). (a) Melting rate $v_{{\varGamma }}$ as a function of the polar angle $\varphi$ defined by the inset at the front on mid-plane at $z=0$ . (b) Tempo-spatially averaged melting rate $\bar {\bar {v}}_{{\varGamma }}$ on the rear interface as a function of the polar angle $\varphi$ . Details of the averaging procedure are given in the text.

At the end of this section, a natural question that arises is whether the wake bifurcations observed during melting can still be related to the classical critical Reynolds numbers established for non-melting rigid spheres, namely $\textit{Re}_{c1}$ , $\textit{Re}_{c2}$ and $\textit{Re}_{c3}$ in figure 2(a). Our results show that this correspondence no longer holds. For instance, the $\textit{Re}_0=300$ case loses vortex shedding only at $\textit{Re}_e \approx 138$ ( $\textit{Re}_{c2}\approx 273$ ), the $\textit{Re}_0=500$ case leaves the chaotic regime around $\textit{Re}_e \approx 206$ ( $\textit{Re}_{c3}\approx 355$ ) and the $\textit{Re}_0=1000$ case retains chaotic features throughout the melting process and only begins to develop a weakly time-periodic behaviour near $\textit{Re}_e \approx 268$ . This systematic reduction of critical Reynolds numbers can be traced to two coupled mechanisms. First, the wake exhibits a hysteresis-like response: in a small test where the melting process was artificially stopped in the $Re_0=300$ case, vortex shedding persisted for several additional periods even though the instantaneous $Re_e$ had fallen well below the corresponding rigid-sphere Hopf bifurcation threshold, demonstrating a delayed decay of unsteady modes. Second, as the solid melts, its shape becomes increasingly disk-like, and such an increase in aspect ratio is known to reduce the bifurcation Reynolds numbers (Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012). These combined effects explain why the classical critical Reynolds numbers cannot be directly applied to melting spheres.

5.1. Drag force

Figure 17(a) presents the transient drag response for the case $\textit{Re}_0 = 200$ . The dashed line indicates the well-established empirical correlation for rigid spheres (Clift, Grace & Weber Reference Clift, Grace and Weber2005),

(5.2) \begin{equation} C_D^H = \frac {24}{\textit{Re}} \big( 1 + 0.15 \textit{Re}^{0.687} \big) + \frac {0.42}{1 + 4.25 \times 10^4 \textit{Re}^{-1.16}}, \end{equation}

valid over $0 \lt \textit{Re} \lt 3\times 10^5$ . Before melting ( $t \leqslant 0$ ), the drag force is in close agreement with (5.2), as expected. Immediately after melting begins, however, $C_D(t)$ exhibits a short-lived but distinct drop for $0 \lt t/t_{\!f} \lt 0.002$ , highlighted in the inset. This transient decrease is analysed in detail very soon. Following this brief interval, $C_D(t)$ increases monotonically as the effective Reynolds number $\textit{Re}_e(t)$ decreases with the shrinking body size.

Figure 17.

(a) Time evolution of the drag coefficient $C_D$ for $\textit{Re}_0=200$ , with the inset providing a magnified view of the onset of melting. (b) Variation of $C_D$ with the effective Reynolds number $\textit{Re}_e(t)$ for different initial Reynolds numbers $\textit{Re}_0$ . In both panels, the dashed line denotes the empirical correlation for quasi-steady drag of rigid spheres given by (5.2).

The overall dependence of $C_D(t)$ on $\textit{Re}_e(t)$ is summarised in figure 17(b) for a range of $\textit{Re}_0$ . Because $\textit{Re}_e(t)$ decreases with time, the trajectories should be interpreted from right to left. At moderate Reynolds numbers ( $25 \leqslant \textit{Re}_0 \leqslant 200$ ), the simulated drag coefficients collapse closely onto the rigid-sphere correlation (5.2). In contrast, at higher $\textit{Re}_0$ systematic deviations arise: the computed drag exceeds the empirical prediction. This discrepancy originates from the evolving morphology of the melting body. While at low-to-moderate $\textit{Re}_0$ the body remains nearly spherical or slightly prolate, at higher $\textit{Re}_0$ the interface develops a bluff, cup-cap-like geometry (§ 3.1), which promotes earlier separation and enhanced pressure drag, thereby exceeding the rigid-sphere correlation.

Moreover, caution is required when interpreting the evolution of $C_D(t)$ as a function of $\textit{Re}_e(t)$ , since the body undergoes simultaneous variations in both size and shape. Such temporal changes can, in principle, induce unsteady hydrodynamic contributions through added-mass and history effects. To the best of our knowledge, no rigorous theory exists for the unsteady drag of a single melting or dissoluble body, even in the idealised case of a perfect sphere. Nevertheless, insight may be gained by analogy with isotropic vapour bubbles in superheated or subcooled liquids, as studied by Legendre, Borée & Magnaudet (Reference Legendre, Borée and Magnaudet1998). Importantly, the added-mass coefficient depends only on body shape and not on the interfacial condition (Mougin & Magnaudet Reference Mougin and Magnaudet2002). Thus, for a spherical body, whether rigid or slippery, the well-known added-mass coefficient is $C_M = 1/2$ . This equivalence strongly suggests that both bubble collapse and sphere melting experience the same added-mass correction. For a sphere with time-varying diameter, the added-mass correction to the drag coefficient may be written as $C_I(t) = 8C_M\mathcal{U}(t)$ , where $\mathcal{U}(t) = D_e'(t)$ is the ratio of the interfacial radial velocity to the free-stream velocity. In the present context, $\mathcal{U}(t)$ can be estimated from the Stefan condition (2.5) as $\mathcal{U}(t) \approx -\textit{St}(\textit{Re}_0\textit{Pr})^{-1}\delta _\theta ^{-1}(t)$ . With $\delta _\theta \approx 1.13\textit{Pr}^{-1/2}\textit{Re}_0^{-1/2}D_e^{1/2}(t)$ (Schlichting & Gersten Reference Schlichting and Gersten2000), this scaling reduces to $\mathcal{U}(t) \approx -0.885St \textit{Pr}^{-1/2}\textit{Re}_e^{-1/2}(t)$ . Substituting $\textit{St} = 0.25$ and $\textit{Pr}=7$ , we obtain $\mathcal{U}(t) \approx -0.0836\textit{Re}_e^{-1/2}(t)$ , and therefore $C_I(t) \approx -0.334\textit{Re}_e^{-1/2}(t)$ . This correction is modest: $C_I \approx -0.106$ for $\textit{Re}_e=10$ , $-0.0236$ for $\textit{Re}_e=200$ and $-0.0106$ for $\textit{Re}_e=1000$ .

In all cases, the added-mass correction remains below $4\,\%$ of the quasi-steady drag coefficient (5.2), and can therefore be neglected for the present simulations at $\textit{St}=0.25$ . At higher Stefan numbers, however, the increased interfacial recession velocity suggests that added-mass effects may become significant. This expectation is confirmed in figure 18. Figure 18(a) shows the evolution of $C_D(t)$ versus effective Reynolds number $\textit{Re}_e(t)$ at fixed $\textit{Re}_0=50$ for $\textit{St}$ between $0.1$ and $2$ . A clear trend emerges: larger $\textit{St}$ systematically reduces $C_D(t)$ at identical $\textit{Re}_e(t)$ . Figure 18(b) quantifies this effect, showing that at $\textit{Re}_e=40$ the instantaneous drag decreases from $1.71$ at $\textit{St}=0.25$ to $1.47$ at $\textit{St}=2$ , consistent with the scaling $\mathcal{U}(t)\sim -St$ derived above. Since figure 18(c) demonstrates that interface morphology is nearly unchanged across cases, the drag reduction cannot be attributed to shape effects but instead reflects the growing role of added-mass contributions at higher melting rates.

Figure 18.

Effect of Stefan number $\textit{St}$ on the temporal evolution of drag at $\textit{Re}_0=50$ . (a) Drag coefficient $C_D$ as a function of effective Reynolds number $\textit{Re}_e(t)$ for varying $\textit{St}$ . (b) Dependence of $C_D$ on $\textit{St}$ at fixed $\textit{Re}_e=40$ . (c) Mid-plane interface ( $z=0$ ) at $\textit{Re}_e=40$ for different $\textit{St}$ , showing negligible morphological differences. Panels (a,c) share the same legend.

As for the history force, Legendre et al. (Reference Legendre, Borée and Magnaudet1998) and Magnaudet & Legendre (Reference Magnaudet and Legendre1998) derived analytical expressions for the history coefficient $C_H(t)$ for bubbles with time-varying radii in the creeping-flow limit. In this regime, $C_H(t)$ originates from the diffusion of vorticity around the bubble and becomes relevant only when $\textit{Re}_e(t)\ll 1$ and $\mathcal{U}\textit{Re}_e(t)\ll 1$ . At higher Reynolds numbers ( $\textit{Re}_e(t)\gg 1$ ), the history effect reduces to a viscous correction consistent with viscous potential flow, scaling as $C_H(t) = 48/\textit{Re}_e(t) - C_D(t)$ . This correction is significant at very early times but decays rapidly as $\textit{Re}_e$ grows. In the present study, $\textit{Re}_e(t)\gg 1$ throughout most of the melting process, and although the viscous-potential assumption is not strictly valid for rigid solids, Legendre et al. (Reference Legendre, Borée and Magnaudet1998) showed that $C_H(t)$ is always smaller than the added-mass correction $C_I(t)$ for vaporising bubbles. Hence, it is reasonable to conclude that the history force is negligible here as well. Taken together, both $C_I(t)$ and $C_H(t)$ make only minor contributions to the total drag, which explains why for $25 \leqslant \textit{Re}_0 \leqslant 200$ the computed drag coefficient $C_D(t)$ follows the quasi-steady prediction (5.2) so closely.

We now return to the abnormal dip observed in figure 17(a) at the onset of melting, $t/t_{\!f} \to 0^+$ . A closer examination reveals that $C_D(t)$ undergoes a sudden reduction immediately after melting begins. Figure 19(a) shows this onset behaviour for $\textit{Re}_0=200$ at different Stefan numbers, $\textit{St}=0.1$ , $0.25$ and $1$ . The magnitude of the drop clearly intensifies with increasing $\textit{St}$ , from $\Delta C_D \approx 0.008$ at $\textit{St}=0.1$ to $\Delta C_D \approx 0.048$ at $\textit{St}=1$ . This trend indicates that faster melting rates induce stronger drag reduction. Additionally, to confirm that the sudden drag drop is solely caused by the melting process rather than any numerical artefact associated with activating the phase change, we conduct an auxiliary test at $\textit{St}=1$ . In this test, melting is switched on at $t=0$ and subsequently switched off at $t=0.5$ . As shown in figure 19(b), once melting is turned off, the reduced $C_D(t)$ rises back to its pre-melting value, as expected. Classical diffusion-controlled history forces are unlikely to explain this effect in the $\textit{Re}_e \gg 1$ regime (Magnaudet & Legendre Reference Magnaudet and Legendre1998). Instead, as Magnaudet & Legendre (Reference Magnaudet and Legendre1998) reported for bubbles, when $\textit{Re}_e\gg 1$ and $\mathcal{U}\textit{Re}_e\gg 1$ , surface recession alters the local radial velocity field, thereby modifying the viscous contribution to the drag. The effect strengthens with larger $\mathcal{U}$ (or $\textit{St}$ ) and is most pronounced at the very onset of melting. This interpretation is confirmed in figure 19(c), which decomposes the drag reduction (black line) for $\textit{St}=1$ into its viscous (blue line) and pressure (red line) components. The two contributions satisfy $\Delta C_{D,\mu } \approx 2\Delta C_{D,p}$ , consistent with the well-known balance that viscous pressure is approximately half of the viscous drag (Legendre & Magnaudet Reference Legendre and Magnaudet1998). This demonstrates that the observed drag reduction originates from viscous, rather than inertial effects. Further evidence is provided in figure 19(d), where the velocity profiles inside the boundary layer at $\varphi =63^\circ$ show that surface recession thickens the viscous boundary layer between $t=0$ and $t=0.2$ , thereby reducing $\partial |\boldsymbol{u}|/\partial \xi$ at the wall and hence the viscous shear stress. This mechanism is consistent with the experimental observations of Vakarelski, Chan & Thoroddsen (Reference Vakarelski, Chan and Thoroddsen2015), who reported drag reduction during the initial stages of sphere melting and attributed it to effective boundary-layer thickening due to interfacial recession.

Figure 19.

Quantitative results of the sudden drop in drag coefficient $C_D (t)$ immediately after melting begins, corresponding to $\textit{Re}_0=200$ . (a) Time evolution of $C_D (t)$ at different Stefan numbers, showing larger drag reductions at higher $\textit{St}$ . (b) Auxiliary test of the drag coefficient $C_D(t)$ at $\textit{St}=1$ , in which melting is switched on at $t=0$ and subsequently switched off at $t=0.5$ . Once melting ceases, $C_D(t)$ returns to its pre-melting value, confirming that the observed drag drop is entirely attributable to the melting process. (c) Decomposition of drag reduction ( $\Delta C_{D}$ , black solid line) into pressure ( $\Delta C_{D,p}$ , red line) and viscous ( $\Delta C_{D,\mu }$ , blue line) contributions for $\textit{St} = 1$ , with dashed line showcasing $\Delta C_{D,\mu }/2$ for reference. (d) Velocity magnitude profiles $|\boldsymbol{u}|$ at $\varphi = 63^\circ$ , comparing $t=0$ (red solid line) and $t=0.2$ (blue dashed line) for $\textit{St}=1$ . Surface recession suddenly thickens the boundary layer and thus reduces viscous shear.

5.2. Lift force

We now turn to the lift force acting on the melting sphere, focusing on the steady planar-symmetric and periodic planar-symmetric melting regimes ( $212 \lt \textit{Re}_0 \lt 355$ ). Figure 20 summarises the evolution of the lift coefficient $C_L$ , plotted against both dimensionless time $t/t_{\!f}$ and effective Reynolds number $\textit{Re}_e(t)$ . In the absence of melting ( $t/t_{\!f} \lt 0$ ), the lift response agrees with prior characterisations of the wake: $C_L$ is steady for $212 \lt \textit{Re}_0 \lesssim 270$ , but becomes oscillatory for $270 \lesssim \textit{Re}_0 \lesssim 300$ , consistent with the transition boundary reported by Johnson & Patel (Reference Johnson and Patel1999).

Figure 20.

Evolution of the lift coefficient $C_L$ in the steady planar-symmetric and periodic planar-symmetric melting regimes ( $212 \lt \textit{Re}_0 \lt 355$ ). (a) Coefficient $C_L$ as a function of dimensionless time $t/t_{\!f}$ , showing a sharp decline near $t/t_{\!f} \approx 0.7$ and eventual reversal at late times. (b) Coefficient $C_L$ versus effective Reynolds number $\textit{Re}_e$ , illustrating that larger $\textit{Re}_0$ accelerates the decline and reversal of lift due to stronger asymmetry in rear-interface melting.

Once melting commences, the early-to-intermediate stages ( $t/t_{\!f} \lt 0.7$ ) retain this pre-melting behaviour, indicating that modest interfacial recession does not immediately reorganise the lift dynamics. Then after that ( $t/t_{\!f} \gt 0.7$ ), $C_L(t)$ undergoes a sharp decline when the morphology change marks the transition from a spherical-cap-like to a cup-cap geometry, as previously displayed in figure 16(a). Physically, it corresponds to accelerated melting of the lower rear interface and the subsequent growth of the lower spiral, which restores approximate vorticity balance across the wake and suppresses the net lift force. For later times ( $t/t_{\!f} \gt 0.8$ ), $C_L$ not only vanishes but eventually reverses sign for all cases considered, implying that a freely translating body would migrate from the positive to the negative $\eta$ direction. The $C_L$ – $\textit{Re}_e$ representation in figure 20(b) highlights that the decay and reversal occur earlier for larger $\textit{Re}_0$ . This enhanced sensitivity arises because higher $\textit{Re}_0$ strengthens the asymmetry between upper and lower melting rates at the rear interface (cf. figure 8), accelerating the reorganisation of the wake.

To further clarify the mechanism underlying lift reversal in the late stage of melting, we analyse in detail the case $\textit{Re}_0=270$ , as presented in figure 21(a). The lift coefficient $C_L(t)$ is decomposed into its pressure and viscous components, $C_{L,p}(t)$ and $C_{L,\mu }(t)$ . The viscous contribution remains nearly constant throughout the process, indicating that shear stresses play only a secondary role in the lift reversal. By contrast, the pressure component dominates and is directly responsible for the reversal of $C_L$ . The inertial origin of this effect can be traced to the relative intensities of the lower and upper spirals, which reverse the sign of $\omega _\zeta ^{\eta \gt 0} - \omega _\zeta ^{\eta \lt 0}$ , and subsequently, such reversal also changes the orientation of the produced counter-rotating streamwise vortices in the wake, thus flipping the sign of $C_{L,p}$ . This process is visualised in figure 21(b), which shows isocontours of streamwise vorticity $\omega _x=\pm 0.25$ at representative time instants ( $t_1$ – $t_4$ shown in figure 21 a). The reversal of lift coincides precisely with a change in the dominant sign of $\omega _x$ , confirming that the transverse force arises directly from the counter-rotating vortex pair in the wake, as anticipated in the classical picture of bluff-body lift generation (Lighthill Reference Lighthill1956). Since $\omega _x$ itself is produced by stretching and tilting of azimuthal vorticity $\omega _\zeta$ , the exchange of dominance between the upper and lower spirals naturally reorganises the wake structure and drives the observed lift reversal.

Figure 21.

Lift-force decomposition and streamwise vorticities at $\textit{Re}_0=270$ . (a) Time evolution of the lift coefficient $C_L$ its pressure ( $C_{L,p}$ ) and viscous ( $C_{L,\mu }$ ) contributions, showing that the reversal of $C_L$ is controlled almost entirely by the pressure component. (b) Isocontours of streamwise vorticity $\omega _x$ on the $x$ – $\zeta$ plane (red: $\omega _x=0.25$ ; blue: $\omega _x=-0.25$ ) at four instants $t_1$ – $t_4$ marked in (a). The exchange in dominance between upper and lower spirals reverses the sign of $\omega _x$ and hence the direction of lift.

5.3. Torque

We now examine the torques acting on the melting sphere, since a freely descending body would rotate in response, thereby introducing an additional degree of freedom that could modify the melting dynamics. For clarity, we focus on the representative case $\textit{Re}_0=270$ , noting that other Reynolds numbers exhibit qualitatively similar behaviour.

Figure 22(a) shows that only the torque component $T_\zeta$ becomes dynamically relevant, growing noticeably as soon as the melt starts $t/t_{\!f} \gt 0$ , while $T_x$ and $T_\eta$ remain negligible. Decomposition of $T_\zeta$ into pressure and viscous contributions (figure 22 b) reveals that the $\zeta$ -component torque is almost entirely pressure-driven. Thus, the asymmetric pressure distribution generates a steadily increasing anticlockwise torque about the $\zeta$ axis, perpendicular to the $C_D$ – $C_L$ plane. Examination of instantaneous shapes (see figure 6) shows that this anticlockwise torque tends to reorient the inclined rear face of the body towards a vertical configuration. If free to rotate, the body would oscillate about this orientation, effectively homogenising the melting rate over its surface. This behaviour is consistent with the experiments of Machicoane et al. (Reference Machicoane, Bonaventure and Volk2013), who reported that freely moving ice spheres in turbulent flows preserved a nearly spherical shape for long durations owing to their ability to rotate.

Figure 22.

Torque evolution experienced by the melting sphere at $\textit{Re}_0=270$ . (a) Time evolution of torque coefficients $T_x$ , $T_\eta$ and $T_\zeta$ . Only the $\zeta$ component becomes significant over time, while $T_x$ and $T_\eta$ remain negligible. (b) Decomposition of $T_\zeta$ into pressure ( $T_{\zeta ,p}$ ) and viscous ( $T_{\zeta ,\mu }$ ) contributions, showing that the torque is almost entirely pressure-driven.

The origin of this anticlockwise torque is clarified in figure 23(a), where the body surface is partitioned into four regions at $t/t_{\!f}=0.502$ in an $x$ – $\eta$ projection. Regions 2 and 4 generate anticlockwise moments (red portions), while regions 1 and 3 produce clockwise moments (blue portions). Owing to the inclination of the body, regions 2 and 4 possess larger effective surface areas. Moreover, region 4 directly faces the oncoming flow, resulting in elevated pressure relative to region 3. These geometric and hydrodynamic factors combine to favour anticlockwise torque. The time evolution of centroid position in figure 23(b) confirms that melting drives a northeastward shift, further amplifying contributions from regions 2 and 4. Quantitative evidence is given in figure 23(c), which compares net contributions: anticlockwise torques dominate progressively, yielding the observed monotonic growth of $T_{\zeta ,p}$ and hence $T_\zeta$ .

Figure 23.

Clarification of the origin of the anticlockwise torque acting on the body for $\textit{Re}_0=270$ . (a) Illustration of the four regions of the interface at the symmetry plane ( $x$ – $\eta$ ) at $t/t_{\!f}=0.502$ , where region 1 and region 3 (blue portions) generate clockwise moments while region 2 and region 4 (red portions) generate anticlockwise moments. The black filled circle represents the mass centre. It is observed that region 4 directly faces the oncoming flow, resulting in elevated pressure relative to region 3. (b) Time evolution of the interface at the symmetry plane ( $x$ – $\eta$ ) and the mass centre represented by the black filled circle, which confirms that melting drives a northeastward shift of the mass centre. (c) Time evolution of the pressure component of the torque $T_{\zeta ,p}$ , with the net contribution from region 1 + region 3 (blue) and region 2 + region 4 (red). For ease of comparison, the red dotted curve represents the absolute values of the red solid curve.

6.1. Vertically translating sphere parallel to gravity

We first consider a melting sphere translating vertically in warm water. In this configuration, the colder melt released from the solid is subject to buoyancy acceleration in the streamwise direction. Depending on the orientation of the translation relative to gravity, two distinct flow responses arise. When the sphere rises against gravity ( $\textit{Ri} \gt 0$ ), the buoyant melt fluid is entrained downstream along the free-stream direction. This assisting buoyancy enhances momentum transport in the boundary layer, stabilises its evolution and delays separation. By contrast, when the sphere descends with gravity ( $\textit{Ri} \lt 0$ ), the cold melt sinks counter to the free stream. This opposing buoyancy destabilises the boundary layer, promotes earlier separation and strengthens the wake recirculation (Kotouč et al. Reference Kotouč, Bouchet and Dušek2009). These competing tendencies substantially alter the interfacial evolution compared with the purely forced-convection case.

To isolate buoyancy effects, we fix the initial Reynolds number at $\textit{Re}_0=200$ and vary the Richardson number in the range $\textit{Ri}\in [-0.5,0.5]$ . Figure 24(a) presents instantaneous cross-sections of the melting body on the vertical mid-plane at $y=0$ , coloured by the local melting rate $v_{{\varGamma }}$ , for four representative cases: $\textit{Ri}=-0.5,\ {-}0.1,\ 0.2$ and $0.5$ . The influence of buoyancy on the separation location is readily apparent. For example, at $\textit{Ri}=0.5$ the assisting buoyancy delays separation considerably, and the body retains an ellipsoidal morphology. In contrast, at $\textit{Ri}=-0.5$ the opposing buoyancy enhances the rear recirculation, producing an exaggeratedly flattened back interface. A more quantitative comparison is given in figure 24(b), which overlays the 2-D interfaces corresponding to different $\textit{Ri}$ when the remaining solid volume has decayed to $V(t)/V_0=0.1$ . The contrasting buoyancy effects are evident: negative $\textit{Ri}$ yields strongly flattened rears, while positive $\textit{Ri}$ produces more rounded geometries.

Figure 24.

Influence of the buoyancy on melting dynamics when the gravity is aligned with the streamwise direction. The initial Reynolds number is maintained at $\textit{Re}_0=200$ . (a) Time evolution of the interfacial shape, coloured by the local melting rate $v_{{\varGamma }}$ , on the mid-plane at $y=0$ for four representative Richardson numbers of $\textit{Ri}=0.5,\ 0.2,\ {-}0.1$ and $-0.5$ . The temporal sampling interval is fixed at $\Delta t =10$ for all cases. (b) Superposition of the mid-plane interface ( $y=0$ ) when the remaining solid volume reaches $V(t)/V_0=0.1$ , highlighting the contrasting influence of stabilising ( $\textit{Ri}\gt 0$ ) and destabilising ( $\textit{Ri}\lt 0$ ) buoyancy.

Furthermore, even in the absence of phase change, i.e. for an adiabatically heated sphere, Kotouč et al. (Reference Kotouč, Bouchet and Dušek2009) demonstrated that buoyancy strongly alters wake dynamics, giving rise to distinct vortical regimes as $\textit{Ri}$ varies. More recently, Li et al. (Reference Li, Jiang, Xu and Zhao2025) extended this analysis to freely translating spheres, identifying attached flow at $\textit{Ri}=0.5$ , axisymmetric wakes at $\textit{Ri}=0.1$ and $0.2$ , steady planar-symmetric flow at $\textit{Ri}=-0.1$ , quasi-periodic vortex shedding at $\textit{Ri}=-0.2$ and fully chaotic flow at $\textit{Ri}=-0.5$ . These canonical regimes are faithfully reproduced in the present simulations of melting spheres. Figure 25 illustrates the contrasting cases of $\textit{Ri}=0.5$ and $\textit{Ri}=-0.5$ : the assisting buoyancy stabilises the wake into an attached-flow configuration, whereas the opposing buoyancy promotes strong recirculation and chaotic vortex shedding. The same sequence of transitions is also evident in figure 24, which compares mid-plane interfacial morphologies across different $\textit{Ri}$ . These results confirm that the stabilising and destabilising effects of buoyancy on wake topology remain the dominant organising principle, even when coupled to phase change.

Figure 25.

Melting dynamics of a sphere at $\textit{Re}_0=200$ with gravity parallel to the streamwise direction. (a) Assisting buoyancy case, $\textit{Ri}=0.5$ . (b) Opposing buoyancy case, $\textit{Ri}=-0.5$ . Shown are snapshots of the temperature field $\theta$ and streamlines on the mid-plane at $z=0$ , together with the 3-D interface coloured by the local melting rate $v_{{\varGamma }}$ , at three representative instants: $t=1, 50, 80$ ( $t/t_{\!f}=0.009,0.469,0.844$ for $\textit{Ri}=0.5$ ; $t/t_{\!f}=0.010,0.497,0.795$ for $\textit{Ri}=-0.5$ ).

6.2. Horizontally translating sphere perpendicular to gravity

We now consider the case where gravity acts perpendicular to the sphere’s translational direction, corresponding to horizontal motion in warm water. In this configuration, the melt-induced cold fluid is always deflected downward ( $-y$ direction) by buoyancy, ensuring that the Richardson number remains negative ( $\textit{Ri}\lt 0$ ). Unlike the vertical case, buoyancy here does not simply assist or oppose the streamwise flow; instead, it perturbs the wake to select a preferred reflectional plane of bifurcation. As a result, flow asymmetry necessarily develops in the $x$ – $y$ plane. To isolate this effect, we fix the initial Reynolds number at $\textit{Re}_0=200$ and examine two representative values, $\textit{Ri}=-0.1$ and $-0.5$ , with $\textit{Ri}=0$ included as a reference.

The results, shown in figure 26, demonstrate that buoyancy acts asymmetrically on the boundary layer: the downward motion of the cold melt stabilises the upper hemisphere ( $y\gt 0$ ), delaying separation, while destabilising the lower hemisphere ( $y\lt 0$ ) and promoting earlier separation. Consequently, the wake and melting dynamics become increasingly asymmetric as $|\textit{Ri}|$ increases. At $\textit{Ri}=0$ , the wake remains axisymmetric throughout the melting process. At $\textit{Ri}=-0.1$ , however, a tilt develops in the wake, producing unequal spirals and a rear interface inclined backward. This configuration closely resembles the steady planar-symmetric regime observed at $\textit{Re}_0=270$ (see figure 6 c). Under stronger buoyancy forcing ( $\textit{Ri}=-0.5$ ), the downward convection nearly suppresses the lower spiral, resulting in a highly asymmetric wake. At later times (left panel of figure 26 c), the interface develops a rounded upper edge and a sharp lower edge, consistent with the experimental observations of Hao & Tao (Reference Hao and Tao2001) due to the asymmetric flow separations.

Figure 26.

Melting dynamics of a horizontally translating sphere at $\textit{Re}_0=200$ with gravity acting perpendicular to the streamwise direction. Cases shown: (a) $\textit{Ri}=0$ , (b) $\textit{Ri}=-0.1$ , (c) $\textit{Ri}=-0.5$ . Left panels: interfacial evolution on the mid-plane at $z=0$ , coloured by the local melting rate $v_{{\varGamma }}$ , with fixed temporal spacing $\Delta t=10$ . Right panels: snapshots of the temperature field $\theta$ and streamlines on the mid-plane at $z=0$ , together with the 3-D interface coloured by $v_{{\varGamma }}$ at $t=40$ ( $t/t_{\!f} = 0.397,0.400,0.408$ for $\textit{Ri}=0,-0.1,-0.5$ ). Increasingly negative $\textit{Ri}$ produces stronger wake asymmetry through stabilisation of the upper and destabilisation of the lower boundary layer.

Finally, we compare the survival time of the melting sphere under different buoyancy conditions, with the initial Reynolds number fixed at $\textit{Re}_0=200$ . Figure 27(a) shows the case of vertical translation. For moderate values $-0.2 \leqslant \textit{Ri} \leqslant 0.3$ , the complete melting time $t_{\!f}$ increases monotonically with more positive $\textit{Ri}$ , reflecting a reduced mean melting rate $\bar {v}_{{\varGamma }}$ . This trend is consistent with the shape effect reported by Yang et al. (Reference Yang, Howland, Liu, Verzicco and Lohse2024a ), who showed that cup-cap morphologies melt faster than slightly prolate ones at $\textit{Re}_0 \sim \mathcal{O}(100)$ . In our simulations, assisting buoyancy ( $\textit{Ri}\gt 0$ ) delays separation and promotes a prolate morphology with lower melt rates, whereas opposing buoyancy ( $\textit{Ri}\lt 0$ ) strengthens recirculation, producing a cup-cap shape with faster decay (see figure 24). At the extremes, deviations from this monotonic trend appear. For $\textit{Ri} \gt 0.3$ , suppression of the wake bubble and the progressively tighter attached flow enhance rear-side melting, thereby reducing $t_{\!f}$ . For $\textit{Ri} \leqslant -0.5$ , strong disordered wake dynamics drives the body towards a more rounded morphology with a smaller aspect ratio (see figure 24 b), which slows the overall melting. As $\textit{Ri}$ decreases further below $-0.5$ , the aspect ratio changes only slightly, and $t_{\!f}$ reaches a plateau. Figure 27(b) shows the case of horizontal translation. Here, $t_{\!f}$ decreases monotonically with more negative $\textit{Ri}$ . The downward buoyancy plume destabilises the lower boundary layer, narrowing the lower spiral, and drives warmer liquid into the wake, thereby enhancing melting of the lower rear surface and accelerating overall decay.

Figure 27.

Dependence of the complete melting time $t_{\!f}$ on the Richardson number $\textit{Ri}$ at fixed $\textit{Re}_0=200$ . (a) Vertical translation (gravity parallel to the streamwise direction): within $-0.2 \leqslant \textit{Ri} \leqslant 0.3$ , $t_{\!f}$ increases with $\textit{Ri}$ , corresponding to a transition from a cup-cap shape to a prolate shape. For $\textit{Ri} \gt 0.3$ , the wake is increasingly suppressed, and melting is enhanced by progressively tighter attached flow. For $\textit{Ri} \leqslant -0.5$ , buoyancy induces chaotic recirculation, and the global melting rate evolution reaches a plateau. (b) Horizontal translation (gravity perpendicular to the streamwise direction): $t_{\!f}$ decreases monotonically with increasingly negative $\textit{Ri}$ , as buoyancy destabilises the lower wake, suppresses the lower spiral and accelerates rear melting.

The analysis of buoyancy effects presented thus far has relied on the assumption of a linear density–temperature relation. However, as noted earlier, the ice–water system exhibits a pronounced density anomaly, for which a quadratic relation is required to capture the local maximum density at $4\,^\circ \textrm{C}$ . In Appendix A, we assess how adopting this quadratic model qualitatively modifies the melting dynamics. The key observations can be summarised as follows.

First, when the ambient water is at $20\,^\circ \textrm{C}$ , the cold fluid released from the melting interface must pass through $4\,^\circ \textrm{C}$ , where the density peaks within the thermal boundary layer. Nevertheless, the overall flow remains dominated by sinking cold water, consistent with the behaviour reported by Weady et al. (Reference Weady, Tong, Zidovska and Ristroph2022) for warmer ambient conditions such as $T_\infty = 8\,^\circ \textrm{C}$ . In this regime, the density anomaly does not fundamentally alter the flow morphology; instead, it primarily shifts the separation position and modifies the size of the recirculation region. Consequently, the qualitative conclusions drawn from the linear density–temperature model remain valid for $20\,^\circ \textrm{C}$ ambient water, as is typical for fluids without density anomalies.

Second, as the ambient temperature decreases and approaches $4\,^\circ \textrm{C}$ , the anomalous rising buoyancy associated with cold fluid near the melting interface gradually overtakes the sinking component. Revisiting the upward-translation configuration in figure 25(a), this increasing dominance of anomalous upward motion leads to earlier flow separation and the formation of a larger recirculation zone. A comparable transition of flow regimes can be reproduced by varying $\textit{Ri}$ from $0.5$ to $-0.5$ , which modifies the buoyancy forcing through a similar underlying mechanism.