2.1. Initial and boundary conditions

The mathematical description of the problem is concluded by prescribing initial and boundary conditions. Initially, the fluid region is filled with quiescent solvent, mathematically

(2.4) \begin{equation} \boldsymbol{u}(\boldsymbol{x}, 0) = \boldsymbol{0}, \quad c(\boldsymbol{x}, 0) = 0 \qquad \mbox{in} \;\;\; \varSigma _l(0). \end{equation}

Fluid velocity at the boundaries satisfies

(2.5) \begin{equation} \boldsymbol{u} = V \hat {n}^\varGamma = (V_{\!x}, V_{\!y} ) \quad \mbox{ on } \quad \varGamma (t), \quad \quad \boldsymbol{u} = (\varOmega _0 y, -\varOmega _0 x ) \quad \mbox{ on } \quad \varGamma _w, \end{equation}

where $V_{\!x}$ and $V_{\!y}$ are the horizontal and vertical components of the interface velocity, respectively, and $\hat {n}^\varGamma = (n_x^\varGamma , n_y^\varGamma )$ is the unit normal to the interface into the fluid region. Although the fluid velocity at the interface is assumed to be the same as the interface velocity, one can also use $V = 0$ (Huang et al. Reference Huang, Shelley and Stein2021). This assumption is relevant because the interface velocity is typically very small ( ${\sim} \boldsymbol{O}(10^{-5})$ ) in dissolution problems. On the other hand, the fluid velocity at the wall $\varGamma _w$ is specified based on a fixed clockwise rotation speed $\varOmega _0 (\geqslant 0)$ . This indicates that the wall is treated as rotating with a constant angular velocity, without any acceleration or oscillation.

At the interface $\varGamma (t)$ , concentration equals the saturation concentration ( $c_{\textit{sat}}$ ) of the solution:

(2.6) \begin{equation} c_{\varGamma } := c(\boldsymbol{x},t) = c_{\textit{sat}} \quad \mbox{ on } \quad \varGamma (t). \end{equation}

As the solute dissolves, the interface recedes as a necessary consequence of mass conservation, described by the flux balance (Wells & Worster Reference Wells and Worster2011; Pegler & Davies Wykes Reference Pegler, Wykes and Megan2021, and references therein):

(2.7) \begin{equation} \rho _s (c_s - c_{\varGamma } ) V = - \rho _0 D \boldsymbol{\nabla }c \boldsymbol{\cdot }\boldsymbol{n}^\varGamma \quad \mbox{ on } \quad \varGamma (t), \end{equation}

where $\rho _s$ is the density of the solid, $c_s$ is the concentration of the undissolved solute and $V$ is the normal velocity of the interface. The condition (2.7) indicates that the solute dissolution rate/speed depends on the concentration flux at the interface as well as the mass diffusivity of the solute in the solvent. Interestingly, in nature, the mass diffusivity of most solute materials is relatively low. As a result, we can expect the dissolution process to be slower than other phase-change processes, such as melting or solidification. However, the current investigation does not focus on any specific solute or solvent. At the outer boundary $\varGamma _w$ , we assume a Neumann condition:

(2.8) \begin{equation} -\!D \boldsymbol{\nabla }c \boldsymbol{\cdot }\boldsymbol{n}^w = 0 \quad \mbox{ on } \quad \varGamma _w \end{equation}

representing no diffusive flux, wherein $\boldsymbol{n}^w$ represents a unit normal out of the fluid region.

2.2. Non-dimensionalisation

Introducing the annular gap width $L$ as the characteristic length scale, we define the characteristic velocity as

(2.9) \begin{equation} u_c = \left ( g \beta _c L \Delta c \right )^{1/2}\!, \end{equation}

which is analogous to the free-fall velocity commonly used in natural convection, obtained by balancing the buoyancy and convective terms appearing in the momentum equation. Using this velocity scale, we render the non-dimensional variables as

(2.10) \begin{equation} \boldsymbol{x}^\prime = \frac {\boldsymbol{x}}{L},\quad t^\prime = \frac {t u_c}{L},\quad \boldsymbol{u}^\prime = \frac {\boldsymbol{u}}{u_c}, \quad p^\prime = \frac {p}{\rho u_{c}^2}, \quad c^\prime = \frac {c}{\Delta c}, \quad \mbox{ where } \quad \Delta c = \max _{t = 0}{(\lvert c - c_{\varGamma }\rvert )}. \end{equation}

The corresponding dimensionless equations read (after dropping the prime symbols)

(2.11) \begin{align} & \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} = 0, \quad \mbox{ in } \quad \varSigma _{l}(t) \times (0, \infty ), \end{align}

(2.12) \begin{align} &\frac {\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u} = - \boldsymbol{\nabla \!} p + \frac {1}{\textit{Re}} {\nabla} ^2 \boldsymbol{u} - c \hat {j}, \quad \mbox{ in } \quad \varSigma _{l}(t) \times (0, \infty ), \end{align}

(2.13) \begin{align} \frac {\partial c}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla } c &= \frac {1}{\textit{Pe}} {\nabla} ^2 c, \quad \mbox{ in } \quad \varSigma _{l}(t) \times (0, \infty ). \\[12pt] \nonumber \end{align}

Equations (2.11)–(2.13) are associated with the initial conditions

(2.14) \begin{equation} \boldsymbol{u}(\boldsymbol{x}, 0) = \boldsymbol{0}, \quad c(\boldsymbol{x}, 0) = 0 \quad \mbox{ in } \quad \varSigma _l(0), \end{equation}

boundary conditions

(2.15) \begin{align} & c = 1 \quad \mbox{ on } \quad \varGamma (t), \qquad \qquad \frac {\partial c}{\partial n^w} = 0 \quad \mbox{ on } \quad \varGamma _w, \end{align}

(2.16) \begin{align} & \boldsymbol{u} = V\! \hat {n^w} \quad \mbox{ on } \quad \varGamma (t), \qquad \qquad \boldsymbol{u} = (\varOmega y, -\varOmega x ) \quad \mbox{ on } \quad \varGamma _w, \end{align}

and the Stefan condition

(2.17) \begin{align} & V = \frac {St}{\textit{Pe}} \frac {\partial c}{\partial n^\varGamma } \quad \mbox{ on } \quad \varGamma (t). \end{align}

A close look into (2.11)–(2.17) yields that the convective dissolution problem can be studied in terms of four dimensionless parameters: Péclet number ( $\textit{Pe}$ ), Reynolds number ( $\textit{Re}$ ) (or, two derived parameters – solutal Rayleigh number ( $Ra = \textit{Re} \boldsymbol{\cdot }\textit{Pe}$ ) and Schmidt number ( $Sc = \textit{Pe}/\textit{Re}$ )); Stefan number ( $St$ ); and $\varOmega$ . These parameters are defined as

(2.18) \begin{align} & \textit{Re} = \frac {u_c L}{\nu }, \;\; \textit{Pe} = \frac {u_c L}{D}, \;\; Ra = \frac {g \beta _c \Delta c L^3}{\nu D}, \;\; Sc = \frac {\nu }{D}, \;\; St = \frac {\rho _0 \Delta c}{\rho _s \Delta c_s}, \;\; \varOmega = \frac {\varOmega _0 L}{u_c}, \end{align}

where $\Delta c_s = c_s - c_\varGamma$ .

2.3. Stream function–vorticity formulation

The stream function–vorticity formulation of the Navier–Stokes equations in a two-dimensional domain facilitates a more comprehensive and intuitive understanding of flow mechanisms than the conventional representation using primitive variables ( $\boldsymbol{u}, p$ ). In a two-dimensional Cartesian framework, the relationships between the stream function $(\psi )$ and vorticity $(\omega )$ are defined as follows:

(2.19) \begin{equation} (u,v) = \left ( \frac {\partial \psi }{\partial y}, - \frac {\partial \psi }{\partial x} \right )\!, \quad \omega = \frac {\partial v}{\partial x} - \frac {\partial u}{\partial y}, \end{equation}

where $ u,v$ denote the velocity components of the fluid in $x,y$ directions, respectively. Using these relations, the governing flow system becomes

(2.20) \begin{align} &\dfrac {\partial ^2 \psi }{\partial x^2} + \dfrac {\partial ^2 \psi }{\partial y^2} = - \omega , \quad \mbox{ in } \quad \varSigma _{l}(t) \times (0, \infty ), \end{align}

(2.21) \begin{align} & \dfrac {\partial \omega }{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }) \omega = \frac {1}{\textit{Re}} {\nabla} ^2 \omega - \frac {\partial c}{\partial x}, \quad \mbox{ in } \quad \varSigma _{l}(t) \times (0, \infty ). \end{align}

The corresponding boundary conditions for $\psi$ are

(2.22) \begin{align} &\psi = 0, \quad \dfrac {\partial \psi }{\partial n^\varGamma } = V_{\!x} n_x^\varGamma + V_{\!y} n_y^\varGamma , \quad \mbox{ on } \quad \varGamma (t), \end{align}

(2.23) \begin{align} & \psi = 0, \quad \dfrac {\partial \psi }{\partial n^w} = \varOmega (x n_x^w + y n_y^w), \quad \mbox{ on } \quad \varGamma _w. \end{align}

To complete the system, boundary conditions for the vorticity $\omega$ are required to be computed from (2.20) utilising conditions (2.22) and (2.23). These are discussed/specified later.

2.4. Gibbs–Thomson effect

The Stefan condition (2.17) indicates that the interface velocity can be determined solely from the parameter ratio $St/\textit{Pe}$ and the concentration flux into the fluid region at the interface. Notably, this movement/velocity may also be influenced by the shape of the interface surface. In convex regions, where the local curvature exceeds the mean curvature, the interface moves more rapidly, enhancing the dissolution process. Conversely, in concave regions, the movement is steadier, resulting in a slower dissolution process. Although the initial circular interface of this model is smooth, it may develop irregularities in areas of both low and high curvature during dissolution. In such instances, it is necessary to equilibrate the interface motion between the concave and convex regions. The Gibbs–Thomson effect effectively addresses this requirement by establishing equilibrium in mass exchange between the solid and fluid regions due to interface effects. In this context, the interface velocity is adjusted as (Perez Reference Perez2005; Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013; Huang et al. Reference Huang, Shelley and Stein2021)

(2.24) \begin{equation} V = \dfrac {St}{\textit{Pe}} \dfrac {\partial c}{\partial n} - \epsilon\! \left ( \kappa - \frac {2 \pi }{S} \right ) \quad \mbox{ on } \quad \varGamma , \end{equation}

where $\epsilon$ is a smoothing term in the interface velocity to enhance numerical stability without affecting the physical dissolution process (Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013; Huang et al. Reference Huang, Shelley and Stein2021). Here, $\kappa$ and $S$ represent the local curvature and total arc length of the interface, respectively. When the interface is circular, the local curvature $\kappa$ is equivalent to the mean curvature given by $2\pi /S$ . As a result, no adjustment to the interface velocity is required for a steadily advancing circular interface.

3.1. Transformed model in computational domain

The initial annular domain of the problem evolves due to the dynamic characteristics of the interface, which may lead to complex shapes. Therefore, in order to solve the problem on a uniform Cartesian grid numerically, the physical domain is transformed at every time step into a fixed unit square $([0,1]\times [0,1])$ using a boundary-fitted transformation (see figure 2). This transformation is implemented so that the physical boundaries of the problem align with the boundary curves of the $(\xi ,\eta )$ coordinates. In the present scenario, we map the inner boundary $\varGamma (t)$ (the interface) to the line $\xi = 0$ and the outer rotating wall $\varGamma _w$ to the line $\xi = 1$ . The other two lines, $\eta = 0$ and $\eta = 1$ , act as artificial boundaries, allowing the computational domain to be defined without resorting to the use of physical boundaries.

Figure 2.

Layout of a typical $10\times 24$ grid in the (a) physical $(x,y)$ domain and (b) computational $(\xi ,\eta )$ domain.

The corresponding Cartesian grid in the $(\xi ,\eta )$ domain is generated by solving the following equations (see Thompson, Thames & Mastin (Reference Thompson, Thames and Mastin1974) and Nandi & Sanyasiraju (Reference Nandi and Sanyasiraju2021) for further details):

(3.1) \begin{align} & \alpha x_{\xi \xi } - 2\beta x_{\xi \eta } + \gamma x_{\eta \eta } + J^{2} \big (Px_{\xi } + Qx_{\eta } \big ) = 0, \end{align}

(3.2) \begin{align} & \alpha y_{\xi \xi } - 2\beta y_{\xi \eta } + \gamma y_{\eta \eta } + J^{2} \big (Py_{\xi } + Qy_{\eta } \big ) = 0, \end{align}

with

(3.3) \begin{align} &\alpha = x_{\eta }^{2} + y_{\eta }^{2}, \quad \beta = x_{\xi } x_{\eta } + y_{\xi } y_{\eta } , \quad \gamma = x_{\xi }^{2} + y_{\xi }^{2}, \quad J = x_{\xi } y_{\eta } - y_{\xi } x_{\eta }. \end{align}

Here, the functions $P$ and $Q$ are used to control the spacing between grid lines. To avoid the interpolation at every time step of calculations, time derivative is calculated in the $(\xi , \eta )$ domain using the following relation:

(3.4) \begin{align} \left (\dfrac {\partial g}{\partial t}\right )_{x, y} & = \dfrac {\partial (x, y, g)}{\partial (\xi , \eta , t)}\big/ \dfrac {\partial (x, y, t)}{\partial (\xi , \eta , t)} \nonumber \\ & = \left (\dfrac {\partial g}{\partial t}\right )_{\xi , \eta } - \dfrac {(g_{\xi }y_{\eta } - g_{\eta }y_{\xi })}{J} \left (\dfrac {{\rm d}x}{{\rm d}t}\right )_{\xi , \eta } + \dfrac {(g_{\xi }x_{\eta } - g_{\eta }x_{\xi })}{J} \left (\dfrac {{\rm d}y}{{\rm d}t}\right )_{\xi , \eta }\!, \end{align}

where $g \in \{ c, \omega \}$ .

Using (3.4), the governing equation for the concentration field $c$ in the computational plane reads as

(3.5) \begin{equation} \left (\dfrac {\partial c}{\partial t}\right )_{\xi , \eta } = A c_{\xi \xi } + B c_{\xi \eta } + C c_{\eta \eta } + D c_{\xi } + E c_{\eta } - \big (u^{\xi }c_{\xi } + u^{\eta } c_{\eta } \big ), \end{equation}

where

(3.6) \begin{align} & A = \dfrac {\alpha }{J^{2}\textit{Pe}},\quad B = - \dfrac {2 \beta }{J^{2}\textit{Pe}},\quad C = \dfrac {\gamma }{J^{2}\textit{Pe}}, \nonumber\\ & D = \dfrac {P}{\textit{Pe}} - \dfrac {1}{J}\! \left (x_{\eta }\frac {{\rm d}y}{{\rm d}t} - y_{\eta } \frac {{\rm d}x}{{\rm d}t}\right ), \quad E = \frac {Q}{\textit{Pe}} - \frac {1}{J}\! \left ( y_{\xi } \frac {{\rm d}x}{{\rm d}t} - x_{\xi } \frac {{\rm d}y}{{\rm d}t}\right )\!. \end{align}

Similarly, the equation for the vorticity field can be expressed as

(3.7) \begin{equation} \left (\dfrac {\partial \omega }{\partial t}\right )_{\xi , \eta } = A_1 \omega _{\xi \xi } + B_1 \omega _{\xi \eta } + C_1 \omega _{\eta \eta } + D_1 \omega _{\xi } + E_1 \omega _{\eta } - \big(u^{\xi } \omega _{\xi } + u^{\eta } \omega _{\eta } \big) + R, \end{equation}

where

(3.8) \begin{align} & A_1 = \dfrac {\alpha }{J^{2}\textit{Re}},\quad B_1 = - \dfrac {2 \beta }{J^{2}\textit{Re}},\quad C_1 = \dfrac {\gamma }{J^{2}Re}, \quad R = - \dfrac {1}{J} ( c_\xi y_\eta - c_\eta y_\xi ), \nonumber\\ & D_1 = \dfrac {P}{\textit{Re}} - \frac {1}{J}\! \left ( x_{\eta } \frac {{\rm d}y}{{\rm d}t} - y_{\eta } \frac {{\rm d}x}{{\rm d}t}\right ), \quad E_1 = \dfrac {Q}{\textit{Re}} - \frac {1}{J}\! \left ( y_{\xi } \frac {{\rm d}x}{{\rm d}t} - x_{\xi } \frac {{\rm d}y}{{\rm d}t} \right )\!. \end{align}

Here, the contravariant fluid velocity components $u^{\xi }$ and $u^{\eta }$ in the $\xi$ and $\eta$ directions are expressed as follows:

(3.9) \begin{equation} u^{\xi } = \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla } \xi = \frac {\psi _\eta }{J}, \quad u^{\eta } = \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla } \eta = -\frac {\psi _\xi }{J}, \end{equation}

which satisfy the continuity equation (2.11). Therefore, in computations, techniques such as upwinding or downwinding can be applied using these contravariant components instead of Cartesian ones. This approach facilitates an accurate numerical representation of fluid flow across different coordinate systems. The stream function and the interface normal velocity in the $\xi \eta$ plane are described by

(3.10) \begin{align} & \alpha \psi _{\xi \xi } - 2\beta \psi _{\xi \eta } + \gamma \psi _{\eta \eta } + J^{2} P \psi _{\xi } + J^{2}Q \psi _{\eta } = - J^{2}\omega , \end{align}

(3.11) \begin{align} & V = - \dfrac {St}{J \textit{Pe} \sqrt {\alpha }} (\alpha c_{\xi } - \beta c_{\eta } ) - \epsilon\! \left ( \kappa - \frac {2 \pi }{S} \right )\!, \quad \mbox{ on } \quad \xi =0 . \end{align}

At the artificial boundaries $\eta = 0$ and $\eta = 1$ , the flow variables are periodic. On the other hand, at the other boundaries ( $\xi = 0, 1$ ), concentration $c$ and stream function $\psi$ are computed using (2.15), (2.22) and (2.23). Finally, the boundary conditions for $\omega$ at these boundaries are derived from the conditions for $\psi$ and (3.10). For example, at the boundary $\xi = 1$ , (2.23) yields

(3.12) \begin{align} & \psi _{\eta } = 0, \quad \psi _{\eta \eta } = 0, \end{align}

(3.13) \begin{align} & \psi _\xi = \dfrac {J \varOmega }{\alpha } ( xy_\eta - yx_\eta ) = p_1 \quad \text{(say)}. \end{align}

This leads to the equation for $\omega$ as follows:

(3.14) \begin{equation} \omega = -\frac {1}{J^2} ( \alpha \psi _{\xi \xi } - 2\beta \psi _{\xi \eta } + J^2 P p_1 ). \end{equation}

A similar approach can be used to obtain the condition for $\omega$ at $\xi = 0$ . It is important to note that the derivatives involved in the $\omega$ boundary condition are computed using a five-point ghost cell technique to incorporate the specified lower-order non-zero derivative condition efficiently.

3.2. Numerical scheme implementation for flow and concentration

While implementing the numerical scheme, the spatial derivatives are discretised using standard second-order central difference approximations, while the convection terms are approximated using a third-order QUICK scheme (Johnson & MacKinnon Reference Johnson and MacKinnon1992; Leonard Reference Leonard1995), as described below:

(3.15) \begin{equation} \big(u^{\xi } c_{\xi }\big)_{i} = \dfrac {u^{\xi }_{i} + \big\lvert u^{\xi }_{i}\big\rvert }{2} \phi ^{+} + \dfrac {u^{\xi }_{i} - \big\lvert u^{\xi }_{i}\big\rvert }{2} \phi ^{-}, \end{equation}

where

(3.16) \begin{align} & \phi ^{+} = \frac {1}{6h} \left ( \phi _{i-2} - 6\phi _{i-1} + 3\phi _{i} + 2\phi _{i+1} \right ), \end{align}

(3.17) \begin{align} & \phi ^{-} = \frac {1}{6h} \left ( -2\phi _{i-1} - 3\phi _{i} + 6\phi _{i+1} - \phi _{i+2} \right ). \end{align}

To compute the stream function $\psi$ , Poisson’s equation (3.10) is discretised with the boundary condition $\psi = 0$ at $\xi = 0, 1$ . The resulting system of algebraic equations is solved using the bi-conjugate gradient-stabilised method with incomplete LU factorisation as a preconditioner, up to a convergence tolerance of $10^{-10}$ of the residual vector. On the other hand, for the concentration $(c)$ and vorticity $(\omega )$ fields, the advection–diffusion equations (3.5) and (3.7) are split in the time direction using the ADI scheme (In’t Hout & Welfert Reference In’t Hout and Welfert2009):

(3.18a) \begin{align} &Y^{(0)} = c^{(n)} + \Delta tf(c^{(n)}), \end{align}

(3.18b) \begin{align} &Y^{(1)} = Y^{(0)} + \lambda \Delta t [ f_{1}(Y^{(1)}) - f_{1}(c^{(n)}) ], \end{align}

(3.18c) \begin{align} &Y^{(2)} = Y^{(1)} + \lambda \Delta t [ f_{2}(Y^{(2)}) - f_{2}(c^{(n)}) ], \end{align}

(3.18d) \begin{align} &Y^{(3)} = Y^{(0)} + \zeta \Delta t [f (Y^{(2)}) - f(c^{(n)}) ], \end{align}

(3.18e) \begin{align} &Y^{(4)} = Y^{(3)} + \lambda \Delta t [ f_{1}(Y^{(4)}) - f_{1}(Y^{(2)}) ], \end{align}

(3.18f) \begin{align} &c^{(n+1)} = Y^{(4)} + \lambda \Delta t [ f_{2}(Y^{(5)}) - f_{2}(Y^{(2)}) ], \end{align}

where $f$ is defined as the sum of three functions: $f_1$ , which consists of derivatives in the $\xi$ direction; $f_2$ , which involves derivatives in the $\eta$ direction; and $f_{12}$ , which includes mixed derivatives and additional terms. The constant parameters $\lambda$ and $\zeta$ determining the stability of the ADI scheme are chosen as $\lambda \geqslant {1}/{2} + {\sqrt {3}}/{2}$ and $\zeta = {1}/{2}$ , which are in accordance with the stability analysis (In’t Hout & Welfert Reference In’t Hout and Welfert2009; Nandi & Sanyasiraju Reference Nandi and Sanyasiraju2021).

The ADI splitting (3.18) involves six intermediate steps for computing the required variables for the next time step, comprising two explicit and four implicit steps. The internal variables $Y^{(0)}$ and $Y^{(3)}$ can be obtained directly from the first and fourth explicit steps, respectively. The other variables $Y^{(1)},\ Y^{(2)},\ Y^{(4)}$ and $Y^{(5)}$ are computed by solving the following four tri-diagonal systems:

(3.19a) \begin{align} &(I - \lambda \Delta t S_{\xi })[Y^{(1)}] = [Y^{(0)}] - \lambda \Delta t S_{\xi }[c^{(n)}] + \lambda \Delta t [B_{1}], \end{align}

(3.19b) \begin{align} &(I - \lambda \Delta t S_{\eta })[Y^{(2)}] = [Y^{(1)}] - \lambda \Delta t S_{\eta }[c^{(n)}] + \lambda \Delta t [B_{2}], \end{align}

(3.19c) \begin{align} & (I - \lambda \Delta t S_{\xi })[Y^{(4)}] = [Y^{(3)}] - \lambda \Delta t S_{\xi }[Y^{(2)}] + \lambda \Delta t [B_{3}], \end{align}

(3.19d) \begin{align} & (I - \lambda \Delta t S_{\eta })[Y^{(5)}] = [Y^{(4)}] - \lambda \Delta t S_{\eta }[Y^{(2)}] + \lambda \Delta t [B_{4}], \end{align}

which are deduced from the remaining four implicit steps. Here, $S_\xi$ and $S_\eta$ are two tri-diagonal matrices derived from discretised equations involving functions $f_{1}(c)$ and $f_{2}(c)$ . The column vectors $B_{l}$ (where $l = 1, \ldots , 4$ ) appearing in (3.19a – d ) represent the boundary conditions. The boundary conditions for the variables at the intermediate time steps are treated the same as those at the working time steps.

It is worth noting that although there are four tri-diagonal systems, two distinct tri-diagonal matrices are sufficient to handle the computations for the required variables at a particular time step. The ADI scheme is preferred over other implicit or explicit methods, such as the Crank–Nicolson and total variation diminishing Runge–Kutta schemes, for solving the convection–diffusion equation for vorticity and concentration owing to the former’s computational economy and unconditional stability.

3.3. Interface tracking

Recall from § 3.1 that $\xi = 0$ corresponds to the interface. Therefore, the unit normal at the interface can be expressed as $\displaystyle \boldsymbol{n}^\varGamma = ( {y_\eta }/{\sqrt {\alpha }}, - {x_\eta }/{\sqrt {\alpha }} )$ . Thus, given the normal velocity  $V$ from (3.11), the interface in the $xy$ plane evolves according to

(3.20) \begin{equation} \dfrac {{\rm d}x}{{\rm d}t} = \dfrac {V\! y_\eta }{\sqrt {\alpha }}, \qquad \dfrac {{\rm d}y}{{\rm d}t} = -\dfrac {V\! x_\eta }{\sqrt {\alpha }}, \end{equation}

which can be expressed in a compact form:

(3.21) \begin{equation} \dfrac {{\rm d} \boldsymbol{X}}{{\rm d}t} = \boldsymbol{g}(\boldsymbol{X}, t), \end{equation}

where $\boldsymbol{X} = (x, y)^\top$ and $\displaystyle \boldsymbol{g} = ( V\! y_\eta /\sqrt {\alpha }, -V\! x_\eta /\sqrt {\alpha } )^\top$ . Equations (3.21) are solved using the third-order total variation diminishing Runge–Kutta scheme:

(3.22a) \begin{align} & \boldsymbol{X}^{(1)} = \boldsymbol{X}^{(n)} + \Delta t \boldsymbol{g}(\boldsymbol{X}^{(n)}, t), \end{align}

(3.22b) \begin{align} & \boldsymbol{X}^{(2)} = \frac {3}{4} \boldsymbol{X}^{(n)} + \frac {1}{4} X^{(1)} + \frac {\Delta t}{4} \boldsymbol{g}(\boldsymbol{X}^{(1)}, t), \end{align}

(3.22c) \begin{align} & \boldsymbol{X}^{(n+1)} = \frac {1}{3} \boldsymbol{X}^{(n)} + \frac {2}{3} \boldsymbol{X}^{(2)} + \frac {2\Delta t}{3} \boldsymbol{g}(\boldsymbol{X}^{(2)}, t). \end{align}

Here, we have used the total variation diminishing Runge–Kutta method instead of the previously used ADI scheme as the ordinary differential equation system (3.20) involves only unidirectional derivative terms along the $\eta$ direction. However, it is worth noting that a first-order explicit scheme can also be used for solving these boundary update equations rather than higher-order methods, as its effects on the accuracy of the numerical solutions are negligible (Udaykumar, Mittal & Shyy Reference Udaykumar, Mittal and Shyy1999).

3.4. Overview of computational framework

Numerical simulations at a given time step consist of several key steps starting with grid generation using (3.1)–(3.2). Next, the stream function $\psi$ is calculated from (3.10) based on the known values of vorticity $\omega$ . Following this, the fluid velocity components $u^{\xi }$ and $u^{\eta }$ are derived from (3.9) using the stream function. The boundary conditions for $\omega$ are then set according to (3.14). Afterward, the concentration $c$ and vorticity $\omega$ are solved separately using the ADI scheme outlined in (3.18). Once these values are obtained, the moving interface is updated using (3.20) with the computed concentration field $c$ , and the grid is regenerated to reflect the new interface location. This cycle of computation is repeated until convergence is achieved. The numerical solution is assumed to converge when the infinity norm of the absolute errors of the variables $\omega , c, x^\varGamma\!, y^\varGamma$ satisfies the specified tolerance of $10^{-8}$ . Once the process converges, the flow field, concentration field and interface are all stored at the current time step.

The code validation and grid independence studies for our computations are presented in Appendices A and B, respectively. The CPU time taken to complete one cycle of computation within the convergence loop is approximately 0.6 s for $50 \times 300$ spatial grid points using MATLAB 2022b in an Intel Xeon(R) Gold 6326 processor with a clock speed of 2.90 GHz and 64 GB of RAM.

4.1. Qualitative impacts of rotation and buoyancy on flow and solute dissolution

In this section, we discuss the spatio-temporal dynamics of the dissolved solute concentration and the flow field for different values of $Ra$ and $\varOmega$ . Figure 3 depicts concentration distribution overlaid by streamlines at four instances corresponding to $A_d(t) = 10\,\%,\ 30\,\%,\ 60\,\% \text{ and } 90\,\%$ (where $A_d$ corresponds to the dissolved area, defined in (B3)) for $Ra = 10^5$ and $\varOmega = 0$ to 2 with an increment of 0.5. At the early stages ( $A_d(t)$ approximately up to 10 %), irrespective of the rotation speed considered here, dissolution and hence the phase-change process are driven by diffusion only (see the first column of figure 3). Thus, the dissolved solute is localised near the interface, and its concentration reduces as we move away from the interface and eventually decays to zero near the wall. However, rotation breaks the left–right symmetry about the vertical axis $x = 0$ , as evident from the first column of figure 3.

Figure 3.

Streamlines (white contours) overlaying the concentration contours of the dissolved solute at $A_d(t) = 10\,\%$ , $30\,\%$ , $60\,\%$ and $90\,\%$ (left to right) for $Ra = 10^5$ and $\varOmega =$ 0–2 with an increment of 0.5 (top to bottom). The corresponding colour map illustrates the concentration variations in the domain at these instances.

In the absence of rotation, the primary vortices, which are observed initially on either side of the solute along the horizontal central line, gradually move towards the bottom of the solute over time and remain there until the dissolution process is complete. On the other hand, with a rotation speed of $0.5$ , a primary vortex forms on the right-hand side of the solute’s surface and eventually settles there even when dissolution reaches the $90\,\%$ stage. However, at rotation speeds of $1,\;1.5$ and $2$ , this vortex ceases to exist after certain stages of dissolution.

Figure 3 also shows that, during dissolution with no rotation, the interplay between diffusion and buoyancy-driven convection aids in developing a laminar velocity boundary layer near the interface except at the topmost part. Rotation induces an additional laminar boundary layer near the cylinder wall. The boundary layer near the interface forms due to the combined effects of buoyancy and centripetal forces. At low rotational speeds, this layer is accompanied by a primary vortex to the right of the interface. As the rotation speed increases, the vortex is pushed away from the interface and eventually disappears. As the cylinder rotates, the fluid in immediate contact with the wall is dragged along, causing the fluid at the surface to move faster than the fluid away from the surface, which leads to the development of an additional boundary layer. An increasing rotation speed induces more shear and larger velocity gradients near the wall, affecting the nature of the boundary-layer thickness.

Figure 4.

Streamlines (white contours) overlaying the concentration contours of the dissolved solute at $A_d(t) = 10\,\%$ , $30\,\%$ , $60\,\%$ and $90\,\%$ (left to right) for $\varOmega =$ 1 and $Ra = 10^5$ , $10^6$ and $10^7$ (top to bottom). The corresponding colour map illustrates the concentration variations in the domain at these instances.

Another important observation is that the flow patterns at rotation speeds greater than  $1$ begin to resemble concentric circular layers around both the solute and the cylindrical wall after a certain amount of dissolution. As the rotation speed increases, the gap between these layers becomes smaller. The flow also loses its closed trajectories and tends to become more uniform and structured.

A similar analysis of the flow regimes is performed by varying the Rayleigh number ( $Ra$ ) while keeping the rotation speed fixed at $\varOmega = 1$ , which corresponds to the case when the rotation-induced laminar velocity of the cylinder is the same as that of the buoyancy-induced characteristic velocity. Figure 4 depicts that the flow varies significantly as the Rayleigh number increases. For $Ra \geqslant 10^6$ , flow tends to become more irregular, and dissolved solute is localised in the right half of the cylinder.

As discussed earlier, in the absence of rotation, buoyancy causes the dissolved solute-laden heavier fluid to sink near the solute boundary and lighter fluid to rise near the cylinder, respecting mass conservation. As dissolution increases, this convective flow progressively becomes confined within the lower half of the cylinder (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.11009). When rotation is added to the cylinder, the upward flow near the cylinder is enhanced on the left. At the early stages of the dissolution, the qualitative feature of the flow remains similar on the left, barring fluid parcels moving faster near the cylinder, supported by the clockwise rotation of the cylinder. On the other hand, on the right, the upward flow of the fluid parcels near the cylinder is opposed by the clockwise-rotating cylinder. As a result, a region of upward-moving fluid parcels is formed, confined between two regions with downward-moving fluid parcels. The former carries fluid parcels atop the undissolved solute, leading to a forced convective flow and more dissolution on the left. With increasing dissolution of the solute, the clockwise flow (on the right) near the solute becomes progressively weaker, and the solute appears to be tilted on the right. We further noticed that the clockwise flow near the solute becomes weaker as the rotation speed increases, leading to only an anticlockwise flow near the solute and a clockwise flow near the cylinder (see supplementary movies 2–4).

Overall, we notice that rotation alters the flow topology that leads to the spreading of the dissolved solute around the interface instead of being localised only in the bottom half of the cylinder ( $\varOmega = 0$ ). This causes the concentration gradient near the upper part of the interface to be shallow, and, hence, the dissolution slows down.

To further explore the behaviour of the flow regimes, we trace trajectories of passive tracers released initially at the solute–fluid interface. These trajectories represent the respective pathlines of individual fluid particles coinciding with the tracer particle and reveal how the flow dynamics changes accordingly. The equations governing the trajectory of a particle $P$ with position $(X_{\!P}(t), Y_{\!P}(t))$ are

(4.1) \begin{equation} \dfrac {{\rm d}X_{\!P}}{{\rm d}t} = u(X_{\!P},Y_{\!P},t), \quad \dfrac {{\rm d}Y_{\!P}}{{\rm d}t} = v(X_{\!P},Y_{\!P},t), \end{equation}

where the Eulerian velocity at its location gives the Lagrangian velocity of a particle. These equations are associated with initial conditions $X_{\!P}(t = 0) = X_{\!P}^0, \; Y_{\!P}(t = 0) = Y_{\!P}^0$ .

Figure 5.

Pathlines of two passive tracers, released initially at $(0.0745, -1.0115)$ and $(-0.0204, -1.9441)$ , for rotation speeds ranging from $0$ to $1.5$ (shown left to right with an increment of $0.5$ ). The respective trajectories of these particles are shown in green and blue. The initial positions of the tracers are marked with open markers, while their final positions are marked with the corresponding filled markers. Arrow markers along the pathlines denote the direction of particle motion. The red contour denotes the final shape of the undissolved solute, while $+$ represents the centre of the cylinder.

Figure 5 shows the trajectories of two tracer particles initially placed at $(X_{\!P}^0, Y_{\!P}^0) = (0.0745$ , $-1.0115)$ and $(X_{\!P}^0, Y_{\!P}^0) = (-0.0204, -1.9441)$ for $Ra = 10^5$ and various rotation speeds ranging from 0 to 1.5. The initial positions of the tracers are chosen judiciously to capture the key features of the flow. The former (subsequently called ‘P1’) starts near the solute, where the flow variation is primarily dominated by buoyancy. On the other hand, the latter (subsequently called ‘P2’) is released near the cylinder, where the flow is dominated by rotation. The trajectories are tracked until 90 % dissolution of the solute is achieved. In the absence of rotation, P1 is driven by natural convection and follows a trajectory vertically downward before it moves upward along the cylinder wall, whereas P2 experiences upward motion only (see supplementary movie 5).

Cylinder rotation significantly alters the particle trajectories. Particles P1 and P2 primarily traverse in the anticlockwise and clockwise directions, respectively. However, for each $\varOmega$ , we observe a unique and distinct characteristic of the particle trajectories. When buoyancy dominates rotation ( $\varOmega = 0.5$ ), the trajectory of P1 is confined to the right half of the cylinder for the majority of the dissolution, except at a later stage when the solute shrinks significantly and the particle has to traverse a small distance to end on the left half of the cylinder. Nonetheless, it eventually ends up on the right half (filled green marker) of the cylinder after spending a small duration on the left. For the same flow parameters, P2 initially rotates in the clockwise direction in accordance with the clockwise motion of the cylinder. However, before completing a full revolution around the solute, the particle is caught in the anticlockwise-rotating fluid parcels on the right, attracting P2 towards the solute. Subsequently, the particle traverses around the solute in the anticlockwise direction and re-enters the clockwise-rotating boundary layer adjacent to the cylinder (see supplementary movie 6).

For $\varOmega \gt 1$ , when rotation dominates buoyancy, P1 remains confined within a neighbourhood of the initial interface position. In particular, for $\varOmega = 1.5$ , it is evident from figure 5 that the particle exhibits a cyclic motion around the solute. The broken symmetry of the trajectory is attributed to the buoyancy that causes a departure of the particle from a perfect cyclic motion. However, as $\varOmega$ increases, the asymmetry becomes less prominent, and the number of revolutions of the trajectory around the solute increases (not shown for brevity). See supplementary movie 7.

Finally, for $\varOmega = 1$ , we observe that, initially, P1 follows a trajectory similar to that for $\varOmega = 0.5$ . However, instead of being confined within the right half of the cylinder, it revolves around the solute as time progresses. After a couple of revolutions around the solute, the particle is slowly pushed away from the interface by a strong buoyant force below the solute. This process continues for a while before the particle is eventually trapped within the boundary layer near the cylinder wall and traverses in the clockwise direction. As anticipated, for both $\varOmega = 1$ and $\varOmega = 1.5$ , P2 rotates in the clockwise direction confined within the boundary layer near the cylinder (see supplementary movie 8).

The above observations of particle paths indicate that the fluid adjacent to the cylinder wall flows in the clockwise direction, while the fluid adjacent to the interface experiences an overall anticlockwise motion. For all the cases under consideration in figure 5, the pathlines are smooth, regular and do not exhibit any disorderly behaviour, a signature of laminar flow. This behaviour is particularly evident at higher rotation speeds, for which the trajectories exhibit smooth circular motions around the solute.

4.2. Quantitative impacts of rotation and buoyancy on solute transport

In this section, we quantify dissolution, mixing and interface shape due to nonlinear interactions between fluid flow and solute dissolution over a wide range of rotation speed ( $\varOmega$ ) and Rayleigh number ( $Ra$ ).

4.2.1. Dissolution time of the solute

Our numerical experiments are carried out until the solute experiences $95\,\%$ dissolution. Using these simulated data, we predict the time required for complete solute dissolution using spline extrapolation. Table 1 lists the time for the complete dissolution for Rayleigh numbers varying by two orders of magnitude and $\varOmega$ ranging from 0 to 1.5. Interestingly, rotation-induced forced convection does not enhance dissolution. Rather, for a fixed Rayleigh number, the dissolution process slows down as the rotation speed increases due to a shallow concentration gradient at the interface. We further note, for $Ra = 10^5$ and $10^6$ , a rotation speed $\varOmega = 1.5$ can significantly delay the complete dissolution of the solute with an almost 85 %–94 % increase in complete dissolution as compared with the case at $\varOmega = 0$ . Whereas, for $Ra = 10^7$ , the time taken for complete dissolution with $\varOmega = 1.5$ increases only up to 8 %–9 % relative to the case of no rotation. In summary, the relative strength of rotation to buoyancy plays a critical role in determining the dissolution of the solute. It appears that there exists a critical rotation speed $\varOmega _{\textit{cric}}(Ra)$ such that for $\varOmega \lt \varOmega _{\textit{cric}}(Ra)$ , dissolution has a weak dependence on rotation. Our prediction further indicates that $\varOmega _{\textit{cric}}$ increases with $Ra$ .

Table 1.

Time required for solute dissolution at various Rayleigh numbers and rotation speeds.

Figure 6.

Concentration of dissolved and undissolved solute in the absence of rotation and buoyancy: (a) initially ( $t = 0$ ) and (b) at a later time ( $t \gt 0$ ). Dissolution recedes the interface (red) radially. Due to this rotational symmetry, the Stefan problem associated with the dissolution can be studied through a one-dimensional model. (c) The concentration profile in the radial direction is shown from top to bottom at $t = 0$ , $t = 7.5$ (early time), $t = 30$ (intermediate time) and $t = 60$ (later time).

We now aim to theoretically estimate a relationship for the solute dissolution time. To achieve this, we first consider a scenario where the effects of buoyancy and rotational forces are neglected. In such cases, the chosen problem is driven solely by diffusion forces, which reduces it to a classical one-dimensional phase-change problem in the radial direction due to the symmetric nature of the system. Let $r_s(t)$ and $r_d(t)$ denote the radius of the undissolved solute and the distance traversed by the solid–liquid interface from its initial position, respectively, at an instantaneous time $t$ (see figure 6). Thus, we have

(4.2) \begin{equation} r_d(t) + r_s(t) = 1 \quad \text{with} \quad r_d(0) = 0 \quad \text{and} \quad r_s(0) = 1. \end{equation}

Overall, the problem can be viewed as the dissolution of a solute block in a fluid, with an expanding fluid region. Resorting to $c(r = r_d(t), t) = 1$ and following the analysis of Stefan (Reference Stefan1891), we obtain $r_d(t)\propto \sqrt {t}$ . Thus, from (4.2) we derive

(4.3) \begin{equation} 1-r_s(t) = k \sqrt {t} \quad \mbox{which yields} \quad A_d(t) = 2k\sqrt {t} - k^2t, \end{equation}

where $k$ is a constant of proportionality. The above-mentioned scaling relation is valid until the dissolved solute reaches the cylinder wall. Therefore, in the present study, the finite-size effect of the fluid region leads to a breakdown of the said power-law behaviour. We performed numerical simulations by neglecting buoyancy and rotation in our model formulation and determined the interface evolution, which agrees excellently with $k \sqrt {t}$ for $k = 0.08$ (correct up to two decimal places) up to $t \lesssim 30$ (see figure 7 a).

Figure 7.

Time required for the gradual dissolution of the solute: (a) without buoyancy or rotation, (b) for various rotation speeds at a fixed Rayleigh number $Ra \approx 10^5$ and (c) for various Rayleigh numbers at fixed rotation speed $\varOmega = 1$ . Solid curves represent numerical results up to 90 % dissolution. The dashed line, proportional to $\sqrt {t}$ , is shown as a guide to the eye. Square markers denote the predicted time for complete dissolution given in table 1.

Effects of $Ra$ and $\varOmega$ on the evolution of $r_d(t)$ and hence $A_d(t)$ are investigated next. We numerically compute $A_d(t)$ , defined in (B3), for different values of $Ra$ and $\varOmega$ . For a fixed $Ra$ , the evolution of $A_d(t)$ is independent of $\varOmega$ at the early stages of the dissolution. For example, with $Ra = 10^5$ , the effects of rotation become evident at $A_d(t) \gtrsim 0.25$ (see figure 7 b). Provided the evolution of the dissolved solute area satisfies the relation (4.3), one infers that $k$ depends on $\varOmega$ in such a way that with increasing $\varOmega$ , $k^2 t$ is dominated by $2 k \sqrt {t}$ in (4.3). This is attributed to the fact that at higher rotation speeds, the rotation-induced effects dominate, and the buoyancy effect becomes negligible, causing the system to behave more like a purely diffusion-driven problem. On the other hand, for a fixed $\varOmega$ , the effects of $Ra$ are evident throughout the dissolution process. However, the overall evolution of $A_d(t)$ can be described by a power-law relation $A_d(t) \sim t^a$ , for $a \in (0.5278, 0.6354)$ .

In an attempt to explain this nonlinear dependence on $Ra$ and $\varOmega$ , we revisit the non-dimensional parameters $Ra$ and $\varOmega$ . In a forced convection, the velocity field is proportional to $\varOmega$ . Thus, we write

(4.4) \begin{equation} \textit{Re} \sim \varOmega , \quad \textit{Pe} \sim \varOmega \quad \mbox{resulting in} \quad Ra \sim \varOmega ^2. \end{equation}

Subsequently, we observe that $A_d(t)$ follows a $\sqrt {t}$ relation provided

(4.5) \begin{equation} \sqrt {Ra_{\varOmega }} \lesssim 250, \quad \mbox{where } Ra_\varOmega = \frac {Ra}{\varOmega ^2}. \end{equation}

Within this specific range for $Ra_\varOmega$ , the flow is dominated by rotation, and the interface is consistently observed to maintain its nearly initial circular shape throughout the dissolution process. On the other hand, when $\sqrt {Ra_{\varOmega }} \gt 250$ , the interface loses its circular shape due to a dominant buoyancy force (see second and third rows of figure 3 and second and third panels of figure 10).

4.2.2. Mixing of the dissolved solute

Here, we quantify mixing of dissolved solute in the solvent and explore the roles of buoyancy and rotation, varying $Ra$ and $\varOmega$ . The degree of mixing is defined as (Jha, Cueto-Felgueroso & Juanes Reference Jha, Cueto-Felgueroso and Juanes2011)

(4.6) \begin{equation} \chi (t) = 1 - \dfrac {\sigma ^2(t)}{\sigma ^2_{\textit{max}}}, \end{equation}

where the variance

(4.7) \begin{equation} \sigma ^2(t) = \langle c^2 \rangle - {\langle c \rangle }^2 \end{equation}

of the concentration field can be calculated in terms of the spatial average defined as

(4.8) \begin{equation} \langle \ast \rangle (t) = \int _{\varSigma _l(t)} \ast (x,y,t) \, {\rm d}x \, {\rm d}y. \end{equation}

Here, $\sigma ^2_{\textit{max}}$ represents the maximum concentration variance observed during the dissolution or mixing process. Equation (4.6) signifies that at the perfectly segregated state ( $\chi = 0$ ) variance attains its maximum value; whereas, for the perfectly mixed state, at which concentration is uniform throughout the fluid region, $\chi = 1$ corresponds to the zero variance.

Figure 8.

Effects of buoyancy and rotation on the variation of the degree of mixing with the amount of dissolved solute, $A_d$ : (a) $Ra = 10^5$ , $10^6$ and $10^7$ (light to dark), without rotation ( $\varOmega = 0$ , solid line) and with rotation ( $\varOmega = 1$ , dashed line); (b) $\varOmega = 0.5$ , $1$ and $2$ (light to dark) for $Ra = 10^5$ (solid line) and $Ra = 10^6$ (dashed line). The black line corresponds to pure diffusion, i.e. without buoyancy or rotation.

Due to the absence of dissolved solute in the fluid region initially ( $c_0(\boldsymbol{x}, 0) \equiv 0$ ), spatial variance is zero and hence $\chi (0) = 1$ . As solute dissolves and mixes with the solvent, $\sigma ^2(t)$ increases, resulting in $\chi (t)$ decreasing until it reaches a local minimum, $\chi _{\textit{min}}^{\textit{local}, 1}$ . We denote the corresponding dissolved solute area using $A_{d, \textit{min}}^{\textit{local}, 1}$ . For $A_d(t) \leqslant A_{d, \textit{min}}^{\textit{local}, 1}$ , $\chi (t)$ is dominated by diffusion with a decay rate that depends on $Ra$ and $\varOmega$ . We note that, depending on $Ra$ and $\varOmega$ , $\chi (t)$ may attain the global minimum value $\chi _{\textit{min}}^{\textit{global}} (= 0)$ . The smallest $A_d(t)$ for which $\chi (t) = 0$ is denoted as $A_{d, \textit{min}}^{\textit{global}}$ . For the reference case of no rotation and buoyancy (black line in figure 8), $A_{d, \textit{min}}^{\textit{local}, 1} = A_{d, \textit{min}}^{\textit{global}}$ , which may also occur for certain combinations of $Ra$ and $\varOmega$ . In the absence of rotation ( $\varOmega = 0$ ), $A_{d, \textit{min}}^{\textit{global}}$ increases as $Ra$ increases (solid lines in figure 8 a). On the other hand, for $\varOmega = 1$ , $A_{d, \textit{min}}^{\textit{global}}$ decreases as $Ra$ increases (dashed lines in figure 8 a). However, no such correlation is obtained on varying $\varOmega\ (\gt 0)$ for a fixed $Ra$ (see figure 8 b). Evolution of $\chi (t)$ beyond $A_{d, \textit{min}}^{\textit{local}, 1}$ is governed by nonlinear interactions of buoyancy and rotation forces. Nonetheless, for a given $Ra$ , at 90 $\,\%$ dissolution of the solute, a better mixed state is obtained through rotation as compared with that in the absence of rotation (figure 8 a). Likewise, for a given rotation speed $\varOmega$ , at 90 $\,\%$ dissolution of the solute, a better mixed state is obtained for smaller $Ra$ (figure 8 b). Thus, we conclude that there is an optimal pair of $Ra$ and $\varOmega$ that causes the dissolved solute to mix better in the solvent.

4.2.3. Symmetry breaking due to rotation

As discussed in § 4.1, in the absence of rotation, the dissolved solute shrinks due to buoyancy force, spreads symmetrically about the $y$ axis and primarily occupies the bottom half of the cylinder – the centre of mass of the dissolved solute is (0, $-y^\ast (t; Ra)$ ). The introduction of cylinder rotation breaks this symmetry. To investigate this symmetry breaking, we compute the centre of mass of the dissolved solute as follows:

(4.9) \begin{equation} \mu _X(t) = \int _{\varSigma _l(t)} x c_X(x,t) \, {\rm d}x, \quad \mu _Y(t) = \int _{\varSigma _l(t)} y c_Y(y,t) \, {\rm d}y, \end{equation}

where

(4.10) \begin{align} & \displaystyle c_X(x, t) = \int _{\varSigma _l(t)} c_{X,Y}(x, y,t) \, {\rm d}y, \end{align}

(4.11) \begin{align} & \displaystyle c_Y(y, t) = \int _{\varSigma _l(t)} c_{X,Y}(x, y,t) \, {\rm d}x \end{align}

and

(4.12) \begin{equation} c_{X,Y}(x, y,t) = \dfrac {c(x, y,t)}{\langle c \rangle (t)}. \end{equation}

The angular departure of the line joining ( $\mu _X(t), \mu _Y(t)$ ) to the centre of the cylinder from the negative $y$ axis, measured in the counterclockwise direction, can be computed as

(4.13) \begin{equation} \theta (t) = \tan ^{-1}\left (\dfrac {\mu _Y(t)}{\mu _X(t)}\right ) - \frac {3\pi }{2}. \end{equation}

As anticipated, in the absence of rotation, the inclination angle is always zero (see dashed line in figure 9 a). We discuss our observations for $Ra \in \{ 10^5, 10^6, 10^7 \}$ and $\varOmega \in \{ 0.5, 1, 2 \}$ . For all the parameter values considered here, $\theta$ remains non-negative, signifying that the centre of mass of dissolved solute always shifts to the opposite direction of rotation of the cylinder. Though the overall variation of $\theta (t)$ with $A_d(t)$ is non-monotonic, at the early stages of dissolution, growing angular departure is dominated by diffusion and rotation. As the amount of dissolved solute increases in the solvent, buoyancy forces become stronger and counter the rotation. Nonlinear interactions of these two counterbalancing forces on the spreading of the dissolved solute result in non-monotonic variations of $\theta (t)$ with $A_d(t)$ . Each local optimum value $\theta _{opt}^{local, i} \; (i = 1, 2,\ldots )$ corresponding to a specific dissolution volume $A_{d, \textit{opt}}^{local, i}$ signifies a reversal of dominance between rotation and buoyancy. To be more specific, for $A_{d, \textit{opt}}^{\textit{local}, 1} \lt A_d(t) \lt A_{d, \textit{opt}}^{\textit{local}, 2}$ buoyancy dominates rotation and $\theta (t)$ decreases, whereas, for $A_{d, \textit{opt}}^{\textit{local}, 2} \lt A_d(t) \lt A_{d, \textit{opt}}^{\textit{local}, 3}$ , rotation dominates buoyancy and $\theta {(t)}$ increases and so on.

Figure 9.

Variation of the angle of inclination, $\theta (t)$ , with the dissolved volume, $A_d(t)$ : (a) fixed Rayleigh number $Ra = 10^5$ and varying rotation speed $\varOmega$ ; (b) fixed $\varOmega = 1$ with varying $Ra$ .

For $Ra = 10^5$ , $A_{d, \textit{opt}}^{\textit{local}, 1}$ ( $\theta _{opt}^{\textit{local}, 1}$ ) decreases (increases) with $\varOmega$ (see figure 9 a). On the other hand, for $\varOmega = 1$ , both $A_{d, \textit{opt}}^{\textit{local}, 1}$ and $\theta _{opt}^{\textit{local}, 1}$ increase with $Ra$ (see figure 9 b). We further note that for $Ra = 10^5$ and $\varOmega = 0.5$ , $\theta (t)$ does not vary significantly as solute dissolves. As rotation increases $(\varOmega \geqslant 1)$ , competing buoyancy and rotation forces result in a rapid oscillation in $\theta$ as solute dissolves, and the frequency of oscillations increases with $\varOmega$ (see figure 9 a). At a higher rotation speed, $\theta (t)$ decreases as $A_d(t)$ increases, signifying a homogeneous distribution of the dissolved solute in the fluid region, which is in agreement with our observations in § 4.1.

4.2.4. Interface evolution

This section examines how the shape of the interface evolves during the dissolution process at various rotation speeds. Figure 10 shows the interface at 15 % dissolution intervals up to 90 % for each rotation speed ranging from $0$ to $2$ . Irrespective of the rotation speed, the interface at the initial stages (up to 15 % dissolution) maintains its circular shape and attains different shapes at later stages depending on the resultant of the forces.

Figure 10.

The shape dynamics of the interface for rotation speeds ranging from $0$ to $2$ (shown from left to right with an increment of $0.5$ ) and $Ra = 10^5$ . For each case, the positions of the interface evolutions are shown (inward with colour changing from light to dark) at every 15 % interval of dissolution with the outermost curve representing the initial interface and the innermost at 90 % dissolution. The filled black circle indicates the centre of the cylinder.

In the absence of rotation, the interface evolves to an egg-like shape as solute dissolves. For the same parameter values, a qualitatively similar evolution has been observed by Huang et al. (Reference Huang, Shelley and Stein2021), in which a circular solute placed at the centre of a rectangular domain developed an egg-like shape during dissolution. At these stages, it appears that the solute moves downward due to buoyancy. However, in reality, the solute does not move from its initial position; rather, the convection-induced preferential dissolution of the upper part of the solute results in an apparent movement of the solute. A similar phenomenon is observed at a lower rotation speed of the cylinder ( $\varOmega \leqslant 0.5$ ). In such cases, the solute is seen to be tilted noticeably in the direction opposite to that of rotation. However, as rotation speed increases ( $\varOmega \geqslant 1$ ), the interface regains its circular shape. When the cylinder rotates sufficiently fast (e.g. $Ra = 10^5$ , $\varOmega \geqslant 1.5$ ), the interface maintains almost a perfectly circular shape exhibiting the axisymmetry of the flow and dissolution. A closer look reveals a slight asymmetry of the interface, which is in agreement with a non-zero $\theta$ for such cases.

Figure 11.

The shape dynamics of the interface for $Ra = 10^5$ , $10^6$ and $10^7$ (left to right) at a fixed rotation speed $\varOmega = 1.0$ . For each case, the positions of the interface evolutions are shown (inward with colour changing from light to dark) at every 15 % interval of dissolution, with the outermost curve representing the initial interface and the innermost at 90 % dissolution. The filled black circle indicates the centre of the cylinder.

A similar analysis is of interface evolution for $\varOmega = 1$ and varying $Ra \in \{ 10^5, 10^6, 10^7 \}$ indicates that the axisymmetry of the interface breaks down as the strength of the buoyancy forces ( $Ra$ ) increases (see figure 11). As discussed in § 4.2.1, dissolution of the solute exhibits identical power-law behaviour for $R_\varOmega \lesssim 250$ . Motivated by this observation, we investigated the evolution of the interface in terms of $Ra_\varOmega$ . As discussed earlier, we have seen through several experiments that the interface always exhibits a circular shape when $\sqrt {Ra_{\varOmega }} \lesssim 250$ and tends to lose this shape for higher values of $Ra_{\varOmega }$ . Here, we summarise the results of nine different numerical experiments in terms of three different values of the non-dimensional group $Ra_\varOmega \gt 250$ .

Figure 12.

Shape dynamics of the interface at 90 % dissolution for three different modified Rayleigh numbers: (a) $\sqrt {Ra_{\varOmega }} = 1000$ , (b) $\sqrt {Ra_{\varOmega }} = 2000$ and (c) $\sqrt {Ra_{\varOmega }} = 4000$ . The filled black circle indicates the centre of the cylinder.

Figure 12 shows an egg-like interface for all nine cases considered here. However, the orientation of the nose (the pointed end of the egg-like shape) of the interface is primarily determined by $\varOmega$ . As expected, the same rotation speed results in a similar interface and a similar orientation of the nose of the interface. Additionally, as $Ra_{\varOmega }$ increases, the nose of the interface becomes sharper, exhibiting behaviour similar to the trend observed when $Ra$ increases, as discussed in figure 11. It is worth mentioning that the above-mentioned behaviour of the interface is limited to $Ra \lt 10^8$ , for which the overall dynamics of the flow field and the solute dissolution exhibit significantly different features, and they are briefly discussed in the next section.

4.3. Dissolution dynamics at $Ra=10^8$

Throughout this paper, we have discussed results for the parameter ranges $10^5 \leqslant Ra \leqslant 6.4 \times 10^7$ and $0 \leqslant \varOmega \leqslant 2.5$ , for which the flow remains regular and exhibits a laminar solutal boundary layer. We close this section by summarising numerical simulations for $Ra = 10^8$ and $\varOmega = 0.5$ in figure 13. We observed that the interface adopts an uneven shape such that the solute region is concave. The solutal boundary layer exhibits characteristics inconsistent with laminar boundary layers (see figures 13 a and 13 b). To better understand this phenomenon, we followed the trajectories of four tracer particles – placed initially on the interface at angular locations 0, $\pi /2$ , $\pi$ and $3 \pi /2$ – as the solute dissolves. The trajectories of these tracers are captured up to 90 % of the solute dissolution and shown in figure 13(c). The loopy and relatively non-smooth particle trajectories indicate an irregular flow. Further quantification of these irregular flow dynamics is beyond the scope of the current study and will be considered in our future research.

Figure 13.

(a) Concentration (undissolved and dissolved solute) distribution at $0\,\%$ (initial), $15\,\%$ , $30\,\%$ , $60\,\%$ and $80\,\%$ dissolution (from left to right) for $Ra=10^8$ and $\varOmega =0.5$ . In each panel, the current position of four particles released at four different positions (leftmost panel) is shown. (b) Shape of the interface at different dissolution volumes from $0\,\%$ to $90\,\%$ with an increment of $15\,\%$ (inward with colour changing from light to dark). (c) Pathlines of passive tracers released at four different positions. In each case, the brown marker represents the starting position, while the green marker indicates the position at $90\,\%$ dissolution.