3.1. The Laurent hierarchy

The mapping function can be expanded in a Laurent series (Galin Reference Galin1945; Polubarinova-Kochina Reference Polubarinova-Kochina1945; Howison Reference Howison1992; Bazant et al. Reference Bazant, Choi and Davidovitch2003; Rycroft & Bazant Reference Rycroft and Bazant2016)

(3.11)\begin{equation} f(\omega, t) = a_1(t) \omega + a_0(t) + a_{-1}(t) \omega^{-1} + \cdots. \end{equation}

In addition to providing a relatively simple solution strategy for the PG equation (3.10), the coefficients $a_m$ offer physical insight into the size and shape of the evolving solid body. Reflection symmetry about the real ($x$) axis implies that the mapping function from the unit circle must also be symmetric, $\textrm {Im} [\,f(\bar \omega )]=-\textrm {Im} [\,f(\omega )]$, which requires that the Laurent coefficients $a_m$ are all real. The coefficient $a_1$ relates the overall size of the solid domain to the radius of the unit circle, and is usually referred to as the conformal radius. The conformal centre $a_0$ approximates the shift in the centre of mass when mapping to the real domain. Points within the solid are bounded by the condition $|z-a_0| \le 2a_1$ (Pommerenke Reference Pommerenke1975), while far from the solid object ($|\omega | \rightarrow \infty$) $z \sim a_1 \omega$. The flux from (3.7) corresponds to a renormalized flow field in the physical domain, with a far field velocity $a_1 U_\infty$. In the mathematical domain the Péclet number must therefore be replaced by the rescaled value $Pe\, a_1$, reflecting the decreasing size of the dissolving cylinder; the flux is therefore explicitly time dependent (Bazant et al. Reference Bazant, Choi and Davidovitch2003).

Substituting the Laurent expansion for $f(\omega )$ (3.11), the evolution of the boundary can be found from (3.10):

(3.12)\begin{equation} \textrm{Re} \left[Pe \sum_{m^\prime=-\infty}^1 \sum_{m^{\prime\prime}=-\infty}^1 {m^{\prime} a_{m^\prime} \dot a_{m^{\prime\prime}} \,\textrm{e}^{-\textrm{i}({m^{\prime}}-{m^{\prime\prime}})\theta}} \right] = -\sigma(\theta, t). \end{equation}

The mapping function can be approximated by a finite number of terms, with coefficients found by taking moments $\int _{-{\rm \pi} }^{{\rm \pi} } \textrm {e}^{\textrm {i}m\theta } \cdots \textrm {d}\theta$ of the evolution equation. For a Laurent series including terms up to $\omega ^{-M}$, there are $M+2$ unknowns, matched by $M+2$ cosine moments

(3.13)\begin{equation} \sum_{m^\prime=-M}^1 \sum_{m^{\prime\prime}=-M}^1\,m^{\prime} a_{m^\prime} \dot a_{m^{\prime\prime}} \left[\delta(m-m^\prime+m^{\prime\prime}) + \delta(m+m^\prime-m^{\prime\prime}) \right] = - \frac{\sigma_m(t)}{{\rm \pi} Pe}, \end{equation}

where $\sigma _m$ are moments of the flux

(3.14)\begin{equation} \sigma_m(t) = \int_{-{\rm \pi}}^{\rm \pi} \cos(m\theta) \sigma(\theta, t)\,\textrm{d}\theta. \end{equation}

The coefficients satisfy the differential equations

(3.15)\begin{equation} \sum_{m^{\prime}=-M}^{1-m} m^{\prime}a_{m^\prime} {\dot a}_{m^\prime+m} + (m^{\prime}+m)a_{m^\prime+m} {\dot a}_{m^\prime} = s_m, \end{equation}

with $s_m(t) = -\sigma _m(t)/({\rm \pi} Pe)$ and $0 \le m \le (M+1)$. Similar equations have been derived in the low Péclet limit (Rycroft & Bazant Reference Rycroft and Bazant2016), but without assuming a reflection symmetry. However, for symmetric objects we recover the same differential equations as Rycroft & Bazant (Reference Rycroft and Bazant2016), but with right-hand sides that are valid at any Péclet number, not just the asymptotic limits. In this work we examine high Péclet (§ 4) and low Péclet (§ 5.2) asymptotics as well as the general case (§ 5). For small Péclet numbers ($Pe \ll 1$) the Laurent series (3.11) can be truncated after the first two terms (§ 5.2) and closed form solutions obtained (Rycroft & Bazant Reference Rycroft and Bazant2016). At $Pe=0$, (3.15) reduces to the hierarchy describing Laplacian growth ($s_m=\delta _{1,m}$), analysed (for example) by Gustafsson & Vasil'ev (Reference Gustafsson and Vasil'ev2006).

Given an expression for the flux around the unit circle – (4.2), (5.5) or (5.6) – the evolution of the boundary can be found by solving a system of differential algebraic equations (DAEs) of the form ${\boldsymbol {F}}({\boldsymbol {y}},\dot {\boldsymbol {y}}) = 0$, with initial conditions $a_1=1$, $a_{m\ne 1} = 0$. Equation (3.15) was implemented as a residual to the IDA solver (https://computing.llnl.gov/projects/sundials), using Python bindings (https://pypi.org/project/Assimulo/) provided by the Assimulo suite (Andersson, Führer & Åkesson Reference Andersson, Führer and Åkesson2015). The IDA implements backwards differencing of variable orders and variable step sizes to maintain the solution within a specified tolerance ($10^{-8}$ in this work). Consistent initial conditions, meaning that ${\boldsymbol {y}}(0)$ and $\dot {\boldsymbol {y}}(0)$ satisfy the DAEs, were implemented by noting that at $t=0$ all the coefficients except $a_1$ are zero, corresponding to a circular shape. The initial derivatives of the Laurent coefficients are then

(3.16)\begin{equation} \dot a_{1-m}(0) = \frac{s_m}{1+\delta_{m,0}}. \end{equation}

Event handling in Assimulo was used to terminate the simulation once the conformal radius was less than twice the specified tolerance.

The area and centre-of-mass position of the dissolving solid can be related to the Laurent coefficients by applying Green's theorem to the complex function $F(z)$ (Rycroft & Bazant Reference Rycroft and Bazant2016)

(3.17)\begin{equation} \iint_A F'(z)\,\textrm{d} x\,\textrm{d} y = -\frac{1}{2\textrm{i}}\oint_C F(z_b)\,\textrm{d}\bar z_b, \end{equation}

where $A$ is the (two-dimensional) domain bounded by the closed curve $C$ and the contour is taken in the counter-clockwise direction. Mapping to the unit circle $\omega _b = \textrm {e}^{\textrm {i}\theta }$, the contour integral can be converted to a definite integral,

(3.18)\begin{equation} \iint_A F^\prime(z)\,\textrm{d} x \,\textrm{d} y = \tfrac{1}{2}\int_{-{\rm \pi}}^{\rm \pi} F(\,f(\omega_b)) \overline{\omega_b f^\prime(\omega_b)}\,\textrm{d}\theta. \end{equation}

We can use this result to calculate the undissolved area $A(t)$ (with $F(z) = z$), and the location of the centre of mass ${\rm \Delta} x(t)$ (with $F(z) = z^2/2$), using the Laurent expansion for $f(\omega )$ (3.11)

(3.19)\begin{gather} A(t) = {\rm \pi}\sum_{{m}=-M}^1\,m a_{m}^2(t), \end{gather}

(3.20)\begin{gather}{\rm \Delta} x(t) = \frac{\rm \pi}{2A(t)} \sum_{{m^\prime}=-M}^1 \sum_{{m^{\prime\prime}}=-M}^1 (m^\prime+m^{\prime\prime}) a_{m^\prime+m^{\prime\prime}}(t)a_{m^\prime}(t)a_{m^{\prime\prime}}(t). \end{gather}

4.1. Comparison with experiment

Numerical solutions of the hierarchy (3.15), with the high Péclet flux moments (4.2), lead to a surprisingly close match with experimentally observed shapes (figure 1), despite the differences in fluid flow (viscous versus potential) and boundary conditions. However, solutions from the PG equation were matched to experimental images with the same undissolved area, rather than the same time. Figure 4 includes a comparison of the area of the undissolved solid (3.19) with experimental measurements from Dutka et al. (Reference Dutka, Starchenko, Osselin, Magni, Szymczak and Ladd2020). The time axis is scaled by the collapse time $t_0$, which is the time taken for complete dissolution of the sample. In the scaled time the areas match to within the uncertainties of the experimental measurements, so that the different underlying physics is reflected only in a change of time scale. The experimental conditions result in a time scale $\tau = 57.3 \ \textrm {h}$ (§ 2), while the sample took $1000 \ \textrm {h}$ to dissolve completely. The dimensionless collapse time from experiment is therefore $t_0 = 17.5$, whereas the collapse time from the PG equation is $t_0 = 14.5$. The size and shape of the undissolved solid can therefore be mapped to the experimental observations by a change in time scale of less than 20 %.

Figure 4.

The area and shift in centre of mass of a dissolving disk as a function of time. The solid lines are a fully convergent solution ($M=128$) of the PG equation (3.15), while experimental measurements (Dutka et al. Reference Dutka, Starchenko, Osselin, Magni, Szymczak and Ladd2020) are shown as solid circles. The dashed lines are single moment approximations to (3.19) and (3.20), $A/A_0 \approxeq a_1^2$ and ${\rm \Delta} x \approxeq a_0$ (§ 4.3). The time scale $t_0$ is the collapse time of the disk, when the dissolving solid vanishes entirely: $t_0 = 14.5$ from the PG equation versus $17.5$ from experiment.

The centre of mass of the disk shifts downstream as it dissolves, but the theoretical prediction is approximately twice as large as the experimentally observed shift (figure 4). Discrepancies in dissolution time scale and drift of the centre of mass have been investigated by numerical simulations and will be discussed in § 4.2. The assumption of potential flow and transport-limited kinetics leads to faster dissolution at the leading edge of the disk than in the experiment, but nevertheless the shapes at matching areas (figure 1) are very similar. We discuss the independence of shape and shift in terms of the Laurent coefficients in §§ 4.3 and 5.3.

We note that the experiment used for figure 4 is not the same as the one shown in figure 1. For figure 1 we selected the experiment with the most regular shapes, but halfway through that experiment the flow rate was quadrupled, changing the time scale. We therefore used another experiment, with a constant flow rate but slightly less regular shapes, for figure 4. A complete description of the experimental results can be found in Dutka et al. (Reference Dutka, Starchenko, Osselin, Magni, Szymczak and Ladd2020); it includes videos of all the experiments, and examples of the images utilized for figures 1 and 4 (figures 11 and 10 respectively from that work).

4.2. Limitations of potential flow and transport-limited kinetics

We have shown in § 4.1 that the differences between the theoretical model described by (3.1)–(3.5a–c) and microfluidic experiments can largely be accounted for by a change in time scale and a shift in the position of the centre of mass; the solid shapes are remarkably similar and the outlines of the boundary overlap almost exactly. Finite-volume simulations of equations enable us to track the effects of the different physics that separates the solutions of (3.1)–(3.5a–c) from the laboratory experiments: Stokes versus potential flow, finite versus infinite Damköhler number and confined versus unbounded systems. The finite-volume simulations are second-order approximations to (2.9a,b)–(2.11a–c), laid out in § 2; the methods have been described in detail elsewhere (Starchenko et al. Reference Starchenko, Marra and Ladd2016; Starchenko & Ladd Reference Starchenko and Ladd2018; Dutka et al. Reference Dutka, Starchenko, Osselin, Magni, Szymczak and Ladd2020).

Figure 5 shows the shape of a dissolving disk from finite-volume simulations (red lines) at different dimensionless times ($t=4, 8, 12$). Images in each row correspond to a particular physical model indicated by the title – potential flow versus Stokes flow and transport-limited dissolution versus a finite reaction rate; the dimensionless time in each column is the same. The black dashed lines are from a conformal mapping of potential flow and transport-limited kinetics in an infinite domain, at times chosen so that the area of the undissolved solid matches that in the finite-volume simulation; those areas are indicated in each panel.

Figure 5.

Finite-volume simulations (red lines) of a dissolving disk (Dutka et al. Reference Dutka, Starchenko, Osselin, Magni, Szymczak and Ladd2020) at $Pe=250$; the simulations are matched with solutions of the Polubarinova–Galin equation (3.15) at the same area. Dimensionless times for the three columns are $t=4$, $t=8$, and $t=12$ respectively, with the time scale from (2.8c). (a) Finite-volume simulations with potential flow and transport-limited kinetics, as used in conformal mapping. A finite-volume simulation with Stokes flow instead of potential flow is shown in (b). (c) Results with the reaction rate characterizing the experiments ($k=4.5 \times 10^{-3} \ \textrm {mm} \ \textrm {s}^{-1}$); the Damköhler number $Da_\mathrm {II} = k/U_\infty = 66.7$. Solutions of the PG equation (black dashed lines) have been scaled in time and shifted to coincide with the finite-volume simulations. The time scales and shifts are shown in figure 6.

Results from a potential flow simulation, with a slip boundary condition on the surface of the dissolving disk, are shown in figure 5(a). Transport-limited kinetics, $c({\boldsymbol {r}} \in S)=0$, were used to match the solution from the PG equation. In this case solutions of the PG equation (3.15) need no rescaling of the time or shift in centre of mass to map onto the finite-volume results. The time scale ($t_C$) and shift (${\rm \Delta} x$) exactly match the finite-volume simulation, as indicated by the red line (with a slope of 1) in figure 6(a) and the red line in figure 6(b).

Figure 6.

Time scale and centre of mass for different flow and transport equations. (a) The rescaled time ($t_C$) of the conformal mapping (potential flow, $Da_\mathrm {II} \rightarrow \infty$, infinite domain) is compared with the dissolution time scale for different physical models ($t$): potential flow, $Da_\mathrm {II} \rightarrow \infty$, unbounded domain (red); Stokes flow, $Da_\mathrm {II} \rightarrow \infty$, unbounded domain (blue); Stokes flow, $Da_\mathrm {II} = 66.7$, unbounded domain (magenta); Stokes flow, $Da_\mathrm {II} = 66.7$, finite domain (black). For potential flow (with $Da_\mathrm {II} = \infty$ and an unbounded domain) the time scales $t_C$ and $t$ are the same, so the slope is unity. The decreased slopes for Stokes-flow simulations reflect the slower dissolution in these cases. (b) The shift in centre of mass for different physical models. The black circles are the experimental results; the corresponding simulations are shown by the black lines.

In Stokes flow the fluid velocity on the disk vanishes, which sets up a viscous boundary layer around the disk with a thickness of order $d$, the gap between the plates of the Hele-Shaw cell. The disk dissolves more slowly than in the potential flow case, as indicated in figure 6(a). In order to match shapes, rescaled times $t_C$ were determined by matching areas between solutions of the PG equation and finite-volume simulations

(4.4)\begin{equation} A(t_C) = A_{sim}(t). \end{equation}

The linear dependence of $t_C$ on $t$ (figure 6a) shows that there is a single rescaling of the time that enables the shapes from conformal mapping to match those from Stokes-flow simulations. The shift in centre of mass is reduced by switching from potential to Stokes flow (figure 6b) and is more in line with the experimental observations, although the shift (blue line) is still too large.

A closer match with the experimental time scale can be achieved by including a finite surface dissolution rate, which in dimensionless form is given by (2.11a–c). Now the shift in centre of mass closely aligns with the experimental measurements shown in figure 4, although the dissolution is slightly slower (magenta line in figure 6a). The remaining discrepancy can be resolved by using a similar sized Hele-Shaw cell, with dimensions $3.8 R \times 3.3 R$ as opposed to an unbounded system. This produces a small change in the centre of mass position (black line in figure 6b) but reduces the time for complete dissolution to $t_0 = 17.5$ (figure 6a) or $1000 \ \textrm {h}$, which precisely matches the experiment.

4.3. Analytic solutions for the low-order Laurent coefficients

Analytical approximations to the Laurent coefficients can be obtained by first truncating the expansion of the mapping function (3.11), beyond the third term. The truncated moment hierarchy from (3.15) becomes

(4.5)\begin{gather} -2a_{-1}^0 \dot a_{-1}^0 + 2a_1^0 \dot a_1^0 = -\frac{8}{\sqrt{{\rm \pi}^3 Pe}}\sqrt{a_1^0}, \end{gather}

(4.6)\begin{gather}-a_{-1}^0 \dot a_0^0 + a_1^0 \dot a_0^0 = \frac{8}{3\sqrt{{\rm \pi}^3 Pe}}\sqrt{a_1^0}, \end{gather}

(4.7)\begin{gather}-a_{-1}^0 \dot a_1^0 + a_1^0 \dot a_{-1}^0 = \frac{8}{15\sqrt{{\rm \pi}^3 Pe}}\sqrt{a_1^0} . \end{gather}

The zero superscript denotes approximate solutions for $a_m$ from the three-term expansion of $f$. Noticing that $a_1$ (solid black line in figure 7) is significantly larger that $a_{-1}$ (solid blue line in figure 7), we can solve for $a_1^0$ by neglecting the small term $2a_{-1}^0\dot a_{-1}^0$ in (4.5),

(4.8)\begin{equation} a_1^0(t) = \left(1 - \frac{t}{t_0^0}\right)^{2/3}. \end{equation}

The collapse time, when $a_1^0$ vanishes,

(4.9)\begin{equation} t_0^0 = \frac{\sqrt{{\rm \pi}^3 Pe}}{6} = 14.674 \end{equation}

is within $1\,\%$ the exact ($M = 128$) value, $t_0 = 14.506$. There are only small differences between the exact numerical solution for $a_1$ (solid black line) and the approximate result $a_1^0$ (dashed black line), mostly near the collapse point; the exact coefficient $a_1$ vanishes at $t_0$, while the analytical approximation $a_1^0$ vanishes at $t_0^0$, which is approximately 1 % larger.

Figure 7.

Time-dependent Laurent coefficients at $Pe=250$. The solid lines are numerical results with $M = 128$. The time has been scaled by the analytical collapse time $t_0^0 = \sqrt {{\rm \pi} ^3 Pe}/6$; the numerical collapse time $t_0 = 0.98857t_0^0$. The dashed lines are the analytical approximations from (4.8) and (4.12a,b).

The remaining coefficients can be found by replacing the time variable in (4.6) and (4.7) by $a_1^0$, noting that $\dot a_1^0 = -2/(3t_0^0 \sqrt {a_1^0})$

(4.10)\begin{gather} \left(-a_{-1}^0+a_1^0\right) \frac{\textrm{d} a_0^0}{\textrm{d} a_1^0} = -\frac{2}{3}a_1^0, \end{gather}

(4.11)\begin{gather}-a_{-1}^0+a_1^0 \frac{\textrm{d} a_{-1}^0}{\textrm{d} a_1^0} = -\frac{2}{15}a_1^0. \end{gather}

Solving first for $a_{-1}^0$ and then for $a_0^0$, the solutions in terms of $a_1^0$ are

(4.12a,b)\begin{equation} a_0^0 = 5 \exp(-15/2)\left[\text{Ei}(15/2)-\text{Ei}(15/2+\ln a_1^0) \right],\quad a_{-1}^0 = -\tfrac{2}{15}a_1^0 \ln a_1^0, \end{equation}

where ${\text{Ei}(x)} = -\int _{-x}^{\infty } t^{-1}\,\textrm {e}^{-t}\,\textrm {d} t$ is the exponential integral function. The analytical solutions from (4.12a,b) are compared with numerical results ($M = 128$) in figure 7. There are only small differences in centre of mass (red lines) between the numerical (solid) and analytical (dashed) results. At the collapse point ($t=t_0^0$) $a_1^0 = a_{-1}^0 = 0$, while the shift in centre of mass reaches a constant value close to the numerical one, $a_0^0 = 5 \exp(-15/2)\text{Ei}(15/2) = 0.800$. The discrepancy in $a_{-1}^0$ (blue lines) in figure 7 comes from the neglect of higher-order Laurent coefficients ($m < -1$), rather than the approximation for $a_1^0$.

4.4. Time-dependent area and centre of mass

In § 4.2, we noted that the impact of potential flow and transport-limited kinetics was largely restricted to the time dependence of the area of the solid domain and the shift in its centre of mass (figure 6); the shapes at the same area were similar for the different physical models (figure 5). Here we examine the time dependence of the conformal area ($a_1$) and centre ($a_0$), which are closely connected to the area and centre of mass of the disk. Figure 4 shows that approximations to the area and centre of mass derived from the first two coefficients in the Laurent expansion (dashed lines) are almost indistinguishable from the exact results (solid lines). Taking the terms involving just $m = 0$ and $m=1$ from (3.19) and (3.20), we have:

(4.13a,b)\begin{equation} \frac{A}{A_0} \approxeq a_1^2,\quad {\rm \Delta} x \approxeq a_0. \end{equation}

The area corresponds approximately to the square of the conformal radius, whereas the centre-of-mass position corresponds to the conformal centre.

Equations (4.13a,b) and (4.8) imply that the area of the object will decrease as $(1-t/t_0^0)^{4/3}$. The scaling is not exact, but a comparison of the solid and dashed black lines in figures 4 and 7 shows that it is very close. Interestingly a similar scaling is predicted for mechanical erosion at large Reynolds numbers (Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013), despite the different mechanisms for material removal. Erosion of friable materials such as clay (Ristroph et al. Reference Ristroph, Moore, Childress, Shelley and Zhang2012) is due to the large shear stresses on the solid surface caused by the flowing fluid. On the other hand, dissolution is a result chemical reactions at the solid surface. However, the interplay of convection and diffusion, momentum in one case (erosion) and mass in the other (dissolution), may offer an explanation for the similarity in the scaling of the area.

In contrast to the shapes in figure 1, with the cusp pointing downstream, erosion produces objects with the cusp pointing upstream (Ristroph et al. Reference Ristroph, Moore, Childress, Shelley and Zhang2012). Nevertheless, in both cases the flux scales as $\sqrt {U_\infty /a_1}$, which is also the case for dissolution at large $Re$ (Huang et al. Reference Huang, Moore and Ristroph2015). This gives a general scaling relation for the evolution of the conformal radius $\dot a_1 \sim a_1^{-1/2}$ which leads to (4.8). The area of two-dimensional bodies then scales as $(1-t/t_0)^{4/3}$ (Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013), while the volume of three-dimensional bodies scales as $(1-t/t_0)^2$ (Huang et al. Reference Huang, Moore and Ristroph2015). However, at low Péclet numbers the area decays approximately linearly in time (§ 5.2) for both dissolution (Rycroft & Bazant Reference Rycroft and Bazant2016) and erosion (Quaife & Moore Reference Quaife and Moore2018).

4.5. Motion of the leading and trailing edges

The position and curvature of the boundary near the leading ($\theta = {\rm \pi}$) and trailing ($\theta = 0$) edges can be written in terms of cumulants of the Laurent coefficients,

(4.14)\begin{equation} S_n^\pm(t) = \sum_{m=-M}^1 (\pm 1)^m m^n a_m(t), \end{equation}

where the minus sign refers to the leading edge, and the plus sign to the trailing edge. The positions and curvatures are given by

(4.15a,b)\begin{equation} z_b^\pm = S_0^\pm,\quad C_b^\pm = \frac{\textrm{d}^2x/\textrm{d}\theta^2}{(\textrm{d} y/\textrm{d}\theta)^2} = - \frac{S_2^\pm}{(S_1^\pm)^2}. \end{equation}

The velocity of the leading edge increases slowly at first and then more rapidly, as shown in figure 8. Similarly to a parabolic tip (Cummings et al. Reference Cummings, Hohlov, Howison and Kornev1999), the velocity of the leading edge is proportional to the square-root of the curvature. The multiplicative factor, ${2/\sqrt {Pe\,{\rm \pi} }}$, can be found by summing the time derivatives of the coefficients

(4.16)\begin{equation} \dot {S}_0^-(0) = \sum_{m=-M}^1 (-1)^m \dot{a}_m(0) = \frac{4}{3 t_0^0}\left(\tfrac{1}{2}+\sum_{m=1}^\infty (-1)^{m-1}\frac{1}{4m^2-1}\right) =\frac{2}{\sqrt{Pe\,{\rm \pi}}}, \end{equation}

where $\dot {a}_m(0)$ given by (3.16). It is larger than the velocity of a parabolic tip (Cummings et al. Reference Cummings, Hohlov, Howison and Kornev1999) by a factor of $\sqrt {2}$, but consistent with the flux at the leading edge of the circle (4.2). The connection between velocity and curvature is approximately maintained at later times, as shown in figure 8.

Figure 8.

Velocity of the leading and trailing edges, $\dot {z}_b^\pm$, at $Pe=250$ ($M = 128$); the plus sign refers to the trailing edge and the minus sign to the leading edge. The curvature of the leading edge (4.15b) is shown as the dashed black line.

By contrast, the trailing edge is essentially stationary ($\dot {S}_0^+(0)=0$) until late in the dissolution, when the effective Péclet number, $a_1(t) Pe$, is small. In reality the trailing edge has a small velocity even at $t=0$ (§ 5); however, a limitation of the asymptotic approximation is that the flux at the trailing edge is predicted to vanish (4.1), whereas in fact it takes a constant value of $1/{\rm \pi}$ (Choi et al. Reference Choi, Margetis, Squires and Bazant2005). However, these differences are not significant when making comparisons with the experiment, given the uncertainties in the data.

5.1. Shapes of a dissolving disk at different Péclet numbers

The shapes of a dissolving disk have been calculated for different initial Péclet numbers ($Pe$); the results are summarized in figure 9. At sufficiently large Péclet numbers ($Pe > 100$), the shapes derived from the asymptotic flux (3.8) are almost indistinguishable from the exact shapes from (5.5) (figure 9a). The sequence of shapes is independent of $Pe$, with a dimensionless time scale proportional to $\sqrt {Pe}$. Since time was made non-dimensional by the fluid velocity (2.8c), the physical time scale is proportional to $Pe^{-1/2}$.

Figure 9.

Dissolution in potential flow by an exact conformal mapping (red solid lines) and the high Péclet asymptotic approximation (black dashed lines). Results are shown at different dimensionless times and at Péclet numbers: $Pe = 100$ (a), $Pe = 10$ (b) and $Pe = 1$ (c). The time for complete dissolution $t_0$ derived from the exact conformal mapping is shown in parentheses.

The disk develops a cusp at the trailing edge, which becomes more pronounced as it dissolves. This shape is characteristic of dissolution experiments and simulations at low $Re$, but high $Pe$ (Soulaine et al. Reference Soulaine, Roman, Kovscek and Tchelepi2017; Dutka et al. Reference Dutka, Starchenko, Osselin, Magni, Szymczak and Ladd2020). The fluid near the trailing edge is nearly saturated with mineral ions and the disk dissolves much more slowly here than on the leading edge. The slow dissolution of the trailing edge also accounts for the positive shift in the centre of mass (figure 4). There does not seem to be an asymptotic shape where only the size ($a_1$) and centre ($a_0$) are time varying. Instead, the tip is constantly evolving in a history dependent manner. Interestingly, high Péclet precipitation does have an invariant shape (Davidovitch et al. Reference Davidovitch, Choi and Bazant2005), although it is not known what initial conditions lead to that shape. We were unable to find a similar factorization for the mapping function of a dissolving object, and numerical results show that an invariant shape does not evolve from an initially circular disk.

By contrast, dissolution at high Reynolds numbers results in a flat trailing edge, due to enhanced mixing in the wake behind the object (Huang et al. Reference Huang, Moore and Ristroph2015). The leading edge of a dissolving cylinder (or sphere) has a finite radius of curvature (Huang et al. Reference Huang, Moore and Ristroph2015), as it does at low $Re$. On the other hand, erosion at high Reynolds number leads to a triangular-shaped tip with straight edges meeting at approximately $90^\circ$ (Ristroph et al. Reference Ristroph, Moore, Childress, Shelley and Zhang2012; Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013).

At lower $Pe$, (3.8) incorrectly predicts a cusp at the trailing edge, while the exact solution is more egg shaped (figure 9b). Furthermore, the asymptotic flux does not capture the persistent circular shape when $Pe < 1$ (figure 9c). The circular shape at low Péclet (${Pe < 1}$) has been found previously (Rycroft & Bazant Reference Rycroft and Bazant2016), using a low Péclet asymptotic expansion. The numerical approach outlined in this section is uniformly valid over the whole range of Péclet numbers, and smoothly connects the low and high-$Pe$ asymptotes. In § 5.2 we improve on the low Péclet expansion in Rycroft & Bazant (Reference Rycroft and Bazant2016) by including logarithmic (in $Pe$) contributions to the flux.

5.2. Low Péclet limit

At low Péclet numbers there is an asymptotic expansion of (5.4) for which the leading terms are (Choi et al. Reference Choi, Margetis, Squires and Bazant2005)

(5.6)\begin{equation} \sigma(\theta, t) = \frac{1 + (1-C)Pe\, a_1(t)\cos \theta + {\textit{O}}(Pe^2)}{C}; \end{equation}

here, $C(t) = -\gamma - \ln (Pe\, a_1(t)/4)$ and $\gamma$ is Euler's constant.

The flux (5.6) has only two moments up to linear order in $Pe$: from (3.14)

(5.7a,b)\begin{equation} \sigma_0 = 2{\rm \pi} C^{-1},\quad \sigma_1 = {\rm \pi}Pe\, a_1 (C^{-1}-1). \end{equation}

With the low-$Pe$ flux (5.6), (4.5)–(4.7) are

(5.8)\begin{gather} -2a_{-1} \dot a_{-1} + 2a_1 \dot a_1 = -\frac{2C^{-1}(t)}{Pe}, \end{gather}

(5.9)\begin{gather}-a_{-1} \dot a_0 + a_1 \dot a_0 = -a_1(C^{-1}(t)-1), \end{gather}

(5.10)\begin{gather}-a_{-1} \dot a_1 + a_1 \dot a_{-1} = 0 . \end{gather}

If the disk is initially circular, $a_{-1} = 0$ for all times (since $\dot a_{-1} \propto a_{-1}$); in fact, if we ignore contributions to the flux from order $Pe^2$ and above, all the Laurent coefficients other than $a_0$ and $a_1$ are zero. Taking $m = -M$ in (3.15), we find from an equation similar to (5.10) that $a_{-M} = 0$. Working up through the hierarchy (3.15), with successively decreasing values of $m$, it can be seen that all the coefficients up to $a_0$ vanish.

In order to analyse the dissolution of complex shapes, Rycroft & Bazant (Reference Rycroft and Bazant2016) argued that the logarithmic term in $C$ could be assumed constant, $C(t) = C_0 = -\gamma - \ln (Pe/4)$. In that case, the area and shift depend linearly on time (Rycroft & Bazant Reference Rycroft and Bazant2016), with a collapse time and final shift in centre of mass

(5.11a,b)\begin{equation} t_0 = \frac{Pe}{2}C_0,\quad {\rm \Delta} x(t_0) = \frac{Pe}{2} (C_0-1). \end{equation}

Taking into account the time scale defined in (2.8c), the dissolution time (in units of $R^2/D$) diverges as $\log Pe$ in the Smoluchowski (diffusion) limit, where the solution corresponds to Laplacian growth.

Taking the time dependence of $C$ into account, and solving (5.8) with the initial condition $a_1(0) = 1$ ($a_{-1}=0$), we obtain an implicit equation for $a_1$,

(5.12)\begin{equation} t = \frac{Pe}{4}\left[\left(1-2\gamma-2\ln\frac{Pe}{4}\right)\left(1-a_1^2\right)+ a_1^2 \ln a_1^2\right]. \end{equation}

Next, we solve for $a_0$ with $a_1$ as the independent variable, as in (4.10) and (4.11). Equation (5.9) can then be rewritten as

(5.13)\begin{equation} a_1 \frac{\textrm{d} a_0}{\textrm{d} a_1}\frac{\textrm{d} a_1}{\textrm{d} t} = -\frac{C^{-1}(a_1)}{Pe}\frac{\textrm{d} a_0}{\textrm{d} a_1} = -a_1(C^{-1}(a_1)-1). \end{equation}

Integrating with respect to $a_1$, and applying the boundary condition $a_0(1) = 0$ ($t=0$), we have

(5.14)\begin{equation} a_0 = \frac{Pe}{4}\left[\left(-1-2\gamma-2\ln\frac{Pe}{4}\right)\left(1-a_1^2\right)+ a_1^2 \ln a_1^2\right]. \end{equation}

We then obtain exact asymptotic expressions for the collapse time and final shift at low $Pe$,

(5.15a,b)\begin{equation} t_0 = \frac{Pe}{2}\left(C_0 + \frac{1}{2}\right),\quad {\rm \Delta} x(t_0) = \frac{Pe}{2}\left(C_0 - \frac{1}{2}\right). \end{equation}

The approximate collapse time (5.11a) is shorter than the exact result (5.15a) by $Pe/4$ and the final shift by $Pe/4$ also ((5.11b) and (5.15b)).

5.3. Dependence of area, shift and shape on time and Péclet number

A log–log plot, shown in figure 10(a), indicates that the undissolved area is approximately a power of the time remaining before collapse $A/A_0 = (1-t/t_0)^\alpha$. At high Péclet numbers $\alpha \approx 4/3$ as indicated in (4.8); however, contributions from Laurent coefficients with ${m < 0}$ modify the power, making it slightly less than $4/3$. Results at $Pe=100$ (solid black line) are asymptotic and show no discernible change at larger Péclet numbers. At lower Péclet numbers the exponents decrease towards a linear decay $A/A_0 = (1-t/t_0)$ (dotted line), which was predicted (Rycroft & Bazant Reference Rycroft and Bazant2016) from a low Péclet asymptotic expansion of the flux (§ 5.2). However, if the Péclet number is reduced below $Pe = 0.01$, there is only a minimal change from the blue line, suggesting that the predicted asymptotic limit, shown by the dotted line in figure 10(a), is never reached. Rycroft & Bazant (Reference Rycroft and Bazant2016) argued that the logarithmic term in $C$ could be assumed to be constant, but this approximation fails to capture the correct evolution of the area near the collapse point. In § 5.2 we developed an exact asymptotic (low $Pe$) calculation for the dissolution of an initially circular disk, which shows that there are logarithmic corrections to the linear decay that occurs for constant $C$.

Figure 10.

Exact solutions for the time-dependent area (a) and collapse time (b). The $x$ axis in panel (a) indicates the time remaining before collapse; note that time is decreasing in the positive direction. The areas are plotted logarithmically to show the approximate power-law scaling; the dashed line is the linear decay at low $Pe$ from Rycroft & Bazant (Reference Rycroft and Bazant2016), while the dotted line is the high $Pe$ asymptote $A(t)/A_0 = (1-t/t_0)^{4/3}$. The numerically determined collapse times are shown in panel (b) as a function of Péclet number. Asymptotic expressions for high and low Péclet numbers, from (4.9) and (5.15a,b), are shown as solid and dashed lines; the approximate low $Pe$ asymptote (5.11a,b) is shown as the dotted line.

The numerically determined collapse times $t_0(Pe)$ are shown in figure 10(b). The collapse time is scaled by the high-$Pe$ asymptotic value $t_0^0 = \sqrt {Pe\,{\rm \pi} ^3/36}$. Since $t_0^0$ is an analytical approximation (4.9), it does not exactly capture the collapse time, even at large $Pe$, but is slightly too large; nevertheless, the collapse time at $Pe=100$ is within 1 % of $t_0^0$. At lower Péclet numbers the disk dissolves faster than predicted by the asymptotic flux (4.2), as indicated in figures 9(b) and 9(c). When $Pe < 1$ the disk remains circular (figure 9c), and all the Laurent coefficients apart from $a_1$ (radius) and $a_0$ (shift) are vanishingly small; an initially circular disk simply shrinks and shifts downstream (Rycroft & Bazant Reference Rycroft and Bazant2016). The exact low $Pe$ asymptote (§ 5.2) is indicated by the dashed line; the approximate asymptote assuming a constant $C = C(0)$ is shown as the dotted line.

The centre of mass of the dissolving disk ${\rm \Delta} x(t)$ is shown in figure 11 over a range of Péclet numbers. At low Péclet numbers the shift is nearly linear in time, but there are small logarithmic corrections (5.14), as discussed in § 5.2. At higher Péclet numbers the shift increases more rapidly with time and undergoes a sharp increase in slope near the collapse point. Nevertheless, there is a well-defined value for the shift at the collapse point ${\rm \Delta} x(t_0)$, which we have determined numerically as ${\rm \Delta} x(t_0, Pe\rightarrow \infty ) = 0.769$. The approximate value from (4.12a,b) is $0.800$, which is $4\,\%$ larger. The asymptotic approximations at high and low Péclet number are shown by the dashed lines in figure 11(b) and the approximate low $Pe$ shift at the collapse point (Rycroft & Bazant Reference Rycroft and Bazant2016) by the dotted line.

Figure 11.

Exact solutions for the shift in centre of mass (${\rm \Delta} x$). (a) The time-dependent shift ${\rm \Delta} x(t)$ for different Péclet numbers. (b) The final shift ${\rm \Delta} x(t_0)$ versus Péclet number. Asymptotic expressions for high and low Péclet numbers, from (4.12a,b) and (5.14), are shown as solid and dashed lines; the approximate low $Pe$ asymptote from (5.11a,b) is shown as the dotted line.

The shape of the dissolving object can be characterized by the aspect ratio

(5.16)\begin{equation} A_R(t) = \frac{4A(t)}{{\rm \pi} L(t)^2} \end{equation}

where the area of the object $A(t)$ is given in (3.19) and the length $L$ is the distance between the leading and trailing edges (4.14),

(5.17)\begin{equation} L(t) = S_0^{+}(t) - S_0^{-}(t) = 2\sideset{}{^\prime}\sum_{m=-M}^1 a_m(t), \end{equation}

here the sum is restricted to the odd coefficients $a_1, a_{-1}, a_{-3}, \ldots$. For an elliptical object, $A_R$ is equal to the ratio of minor to major axes. The aspect ratio varies from one at $t=0$ (circular disk), to a smaller value at $t_0$ (elliptical disk) resulting from the asymmetry introduced by the flow (figure 12a). For small Péclet numbers ($Pe < 0.1$) the aspect ratio remains constant, meaning that the disk is always circular. However, at $Pe = 1$ there is a growing asymmetry, not easily visible in figure 9(c), with $A_R \approx 0.9$ near the collapse point. At higher Péclet numbers the asymmetry becomes much larger, with a terminal aspect ratio $A_R(t_0) = 0.027$ at $Pe = 1000$ (figure 12b).

Figure 12.

Exact solutions for the aspect ratio of the dissolving disk. (a) The time-dependent aspect ratio $A_R(t)$ for different Péclet numbers. (b) The final aspect ratio $A_R(t_0)$ versus Péclet number.

5.4. The interplay of fluid velocity and conformal radius

The effective Péclet number $Pe\,a_1(t)$ decreases with time as the solid dissolves, and eventually the dynamics is controlled by low $Pe$ asymptotics (Rycroft & Bazant Reference Rycroft and Bazant2016). Figure 13 shows the shapes of a dissolving cylinder at high initial Péclet number ($Pe=1000$). The surface coordinates have been scaled by the conformal radius $a_1$ so that the shapes near the collapse point remain visible; the coordinates are shifted along the $x$-axis to keep the centre of mass at the origin. The time $\delta t = t_0 - t$ indicates the time remaining before the solid vanishes.

Figure 13.

Dissolution at different times to collapse, $\delta t = t_0 - t$, for a Péclet number $Pe=1000$. Outlines of the shape of the undissolved solid (red lines) are scaled by the conformal radius ($a_1$), which is indicated in the legend. The blue lines show the region within 0.05 of the tip on a scale expanded sixteen fold. The inset plots (blue) have been shifted a small distance downstream to separate them from the main figures.

At high Péclet number ($Pe=1000$) the disk develops a distinct cusp at the trailing edge, with slightly concave regions near the tip. By contrast, at lower Péclet numbers (figure 9) the trailing tip is always convex. However, although the tip looks sharp on the scale of the plot, an expanded scale shows that it remains slightly rounded (blue inset figures).

Starting with a high Péclet number means that the undissolved solid must be very small before the low-$Pe$ limit is approached. The collapse to the finite-time singularity ($\delta t \rightarrow 0$) is then very abrupt, as illustrated by the logarithmic time scale in the last row of figure 13. In high velocity erosion, an invariant triangular shape develops on the leading edge, which shrinks and shifts downstream with time (Moore et al. Reference Moore, Ristroph, Childress, Zhang and Shelley2013). Here, the shape of the solid body is evolving all the way to the end; there is no invariant shape.

Once $Pe\,a_1$ drops below unity, all Laurent coefficients other than $a_1$, $a_0$, and $a_{-1}$ become small. The undissolved solid then takes up a roughly elliptical shape, which persists until collapse; the final shape can be seen in the last panel of figure 13. It may seem surprising in view of § 5.2 and figure 9(c) that the final shape is elliptical not circular. However, when the initial Péclet number is large, the object enters into the low $Pe$ asymptotic regime with a non-circular shape, so that the higher-order Laurent coefficients are decaying in time, but remain non-zero (Rycroft & Bazant Reference Rycroft and Bazant2016).