The dissolution of a solid spherical particle is a canonical problem that finds many areas of application. In this work, we provide a generalized theory that takes into account the role of forced convection in the solvent (or, alternatively, the relative motion of the particle in the solvent), showing that the problem can be formulated in terms of four dimensionless parameters. Focusing on the case when one of these (the Reynolds number) is small, we consider asymptotic and numerical approaches to the problem, with a key outcome being a numerical method, implemented in the finiteelement software Comsol Multiphysics, that is able to solve the resulting axisymmetric movingboundary problem, even when over 90% of the particle has dissolved and its shape is far from spherical. We also demonstrate how this approach relates to and supersedes earlier efforts, providing a quantitative assessment of the often unquestioningly used Ranz–Marshall correlation for mass transfer from a sphere. In particular, we find that this correlation may overpredict the dissolution time by a factor of four, whereas a correlation by Clift et al. that is cited and used less often performs considerably better, even in the highly convectiondominated regime for which it was not originally intended.
I. INTRODUCTION
The dissolution of a solid spherical particle is a canonical problem that is found in many industrial and consumer applications, ranging from pharmaceutical and food products, to chemicals, detergents, and paints;^{9} for example, in the pharmaceutical industry, it is wellestablished that the dissolution of drug particles is among the most crucial steps underlying the performance of solid dosage forms.^{13,31} As a consequence, mathematical modeling constitutes an indispensable complement to experimental studies of dissolution.
The modeling of solid particle dissolution in a background solvent can arguably be traced back to the late 19th century, when Noyes and Whitney^{26} postulated that the dissolution rate is linearly proportional to the difference between the particle's solubility, $ C s ,$ and the concentration of dissolved solid in the bulk solution, C_{t}. Soon after, this evolved into the wellknown Nernst–Brunner equation by treating dissolution as a diffusioncontrolled process using Fick's first law of diffusion;^{7,25} in fact, the Nernst–Brunner formalism remains the dominant model describing drug particle dissolution behavior in the pharmaceutical industry to date.^{32}
A key assumption underlying the Nernst–Brunner formalism is the presence of a static layer immediately adjacent to the surface of the particle, in which the dissolved substance diffuses. While the thickness of the layer, h, cannot be measured in reality, it is a manifestation of the resistance to mass transfer of molecules in dissolution. Based on fluid dynamics models, for a system where convective mass transfer takes place, such as drug particles suspended in fluid, the diffusion layer thickness is calculated as the ratio of the characteristic length, L, to the Sherwood number, Sh,^{34} which itself can be thought of as the ratio of convective mass transfer to diffusive mass transfer. For spherical particles, L is considered to be equal to the diameter of the particle, and Sh can be determined by the Ranz–Marshall equation,^{28} for example. When there is no convective mass transfer, h is taken to be equal to the radius of the particle.
Although h is computed by considering the mass transfer in threedimensional (3D) space, as described by the Ranz–Marshall correlation,^{28} the Nernst–Brunner equation itself was derived under the assumption that diffusion takes place under onedimensional (1D) conditions. This is because only in 1D can the concentration gradient across the diffusion layer be a constant value, denoted as $ ( C s \u2212 C t ) / h$. However, realistic solid particle dissolution is not a 1D phenomenon, as actual dissolution occurs in 3D; this limitation was raised by Wang and Flanagan.^{45,46} By solving Fick's first law using a nonconstant, distancedependent concentration gradient, they showed that the widely used “cubic root law” in drug dissolution modeling, which was derived from the Nernst–Brunner formalism, is valid only when the particle size is substantially greater than the diffusion layer thickness.^{45} In other words, modeling of particle dissolution using the classical 1D approach, though widely adopted, may misrepresent the dissolution behavior in reality.
More recently, in light of these constraints, So et al.^{32} presented a model capable of describing drug particle dissolution under 3D conditions, but still following the same fundamental principle underpinning the Nernst–Brunner formalism, i.e., that dissolution is driven by diffusion within the diffusion layer. Specifically, they proposed that particle dissolution behavior can be modeled by solving Fick's second law of diffusion in 3D space. Solving Fick's second law gives rise to the solute concentration as a function of time and position within the diffusion layer. Once a solution is identified, the particle dissolution rate can be subsequently obtained by subjecting the solute concentration to Fick's first law of diffusion. This approach has an advantage over Wang and Flanagan's method,^{45} in which only Fick's first law was employed. The advantage arises from the fact that Fick's first law is constrained to describing the steadystate diffusion only, whereas Fick's second law covers both steady and nonsteady state processes.
Against this background, the purpose of this paper is to derive a model for particle dissolution without many of the assumptions mentioned above. In contrast to the principles underlying the Noyes–Whitney and Nernst–Brunner equations,^{7,25,26} we will not assume the presence of a diffusion layer adjacent to the surface of the particle; neither do we assume the validity or otherwise of the Ranz–Marshall correlation. Instead, the model development is driven by conservation laws for mass, momentum and dissolved species concentration as dissolution proceeds. Although the initial focus will be for a particle dissolving in solvent that streams past it, i.e., forced convection, the formulation will also encompass the case of a stationary fluid, as well as being able to delineate when natural convection may become dominant; this can occur because the dissolved species can lead to local variations in the density of the solvent. In this way, the derived framework gives a description that is arguably more general and transparent than those documented previously, which have frequently focused on diffusiondominated dissolution.^{8,10,12,14,18,29,37,47} Moreover, although we consider a particle that is initially spherical, our framework is able to account for vast deviations from this shape, which are a natural consequence of dissolution as the effect of convection increases; by extension, the framework can also account for particles that are not initially spherical. In addition, we note that, although the literature review given above is far from complete as regards all aspects of earlier experimental and theoretical work on the dissolution of solid particles, it is sufficient for our purposes; for a more complete recent review, see Cao et al.^{9}
In what follows, we formulate dissolution as a movingboundary problem that involves the interplay of the hydrodynamic flow of the solvent and the advection and diffusion of solute in the solvent. The resulting equations are then tackled by employing a combination of asymptotic analysis—both regular and singular perturbation theory  and numerical simulation; for the latter, the finiteelement software Comsol Multiphysics is used. The results from the asymptotic analysis help to confirm the accuracy of the numerical simulations in certain limits, while the latter is used to obtain solutions in parameter regimes when the former cannot be used. The advantage of this twopronged approach is that it enables us to obtain an overall view of this canonical problem, in contrast to the majority of earlier work, which considers either the diffusiondominated limit, in which case the particle can be considered as a shrinking sphere, or includes convection, but does not allow for the fact that the dissolving particle will not remain spherical. Indeed, with our numerical simulations, we are able compute the shape of the remaining particle during dissolution; an original finding is that the greater the convection, the greater the degree of asymmetry in dissolution with respect to the front and rear of the particle.
The structure of the paper is as follows. In Sec. II, we formulate the problem; in Sec. III, the model equations are nondimensionalized and analyzed asymptotically in two differing limits: one in which diffusion of the dissolved species dominates convection (Sec. III B) and the other in which convection dominates diffusion (Sec. III C). In Sec. IV, we describe the method that we use to solve the full equations numerically. In Sec. V, the results of the asymptotic analysis and the full numerical solutions are compared; in general, very good agreement is obtained. Conclusions are drawn in Sec. VI. In addition, three Appendixes provide supplementary details of the analysis: Appendix A gives the corresponding equations that prevail if natural convection dominates forced convection; Appendix B considers the diffusiondominated regime, providing a previously overlooked analytical solution in the limit of slow dissolution; and Appendix C discusses simplified models for forcedconvectiondominated dissolution.
II. MODEL EQUATIONS
Consider an initial configuration at time t = 0 consisting of a solid sphere of radius a that is fixed at the origin, past which a solvent of infinite expanse flows at a constant speed, U, in the xdirection; the situation can, therefore, be assumed to be axisymmetric. For $ t > 0 ,$ the sphere is able to dissolve; because of the configuration, the resulting body will also be axisymmetric, and it will later on prove to be convenient to work in spherical polar $ ( r , \theta , \varphi )$ coordinates. Denoting the resulting surface by $ r = R ( \theta , t ) ,$ the proposed evolution is as shown in Fig. 1. Let c_{0} be the initial solid concentration in the sphere, where we shall suppose that $ c 0 > c sat ,$ with c_{sat} denoting the saturation concentration of the dissolving solid in the liquid. In addition, let c denote the concentration of the solute in the solvent; in view of the above discussion regarding axisymmetry, we will have that $ c = c ( r , \theta , t ) .$ We will also assume that $ c = c \u221e ( \u2265 0 )$ for $ r > R ( \theta , t )$ at $ t = 0 ,$ where $ c \u221e < c sat ;$ even though we will ultimately take $ c \u221e = 0 ,$ we have allowed for the possibility that $ c \u221e > 0 ,$ as this is done elsewhere,^{8,18,29} although this will have no impact on the analysis.
A final detail that we wish to note is that, in the formulation above, $ \rho s = c 0 ,$ if one decides to interpret concentration in terms of mass per unit volume, as is indeed often done in the pharmaceutical literature.^{15,16,24} Of course, if one interprets the concentration in terms of the number of moles per unit volume, then $ \rho s = c 0 M ,$ where M is the known relative molecular mass of the particle material. However, although ρ_{s} and c_{0} are clearly related to each other in this particular problem, we will nevertheless adopt distinguishing notation, in view of potential future extensions of the model: for example, in a pharmaceutical context, the particle is likely to consist of both drug and excipient which can dissolve at different rates, meaning that the density of the particle cannot be directly related to a single concentration in a trivial way.
III. ANALYSIS
A. Nondimensionalization
In general, it would be necessary to solve Eqs. (42)–(44), subject to (28)–(36), numerically, and this will be done later. However, we can already note that there is a partial decoupling, in that c does not appear in Eqs. (42) and (43). Nevertheless, the system of equations does not decouple entirely, since the solution for u will depend on the location of the dissolution front; this depends on c, which clearly depends on u via (44).
In the ensuing Secs. III B and III C, we consider asymptotic solutions that can be obtained analytically or quasianalytically for $ P e \u226a 1$ and $ P e \u226b 1.$ On the other hand, when $ P e \u223c O ( 1 ) ,$ no such simplified analysis is possible; consequently, we do not consider it, other than solely numerically in Sec. V.
B. $ P e \u226a 1$
The analysis would now proceed differently, depending on whether $  \rho \u2212 1  / P e \u2009 \u2272 \u2009 1$ and $  \rho \u2212 1  / P e \u226b 1.$
1. $  \rho \u2212 1  / P e \u2272 1$
For the case when $  \rho \u2212 1  / P e \u226a 1 ,$ we need simply to set ν = 0 in (59).
2. $  \rho \u2212 1  / P e \u226b 1$
However, strictly speaking, Eq. (97) is erroneous, as it was obtained by assuming that $ R ( 0 )$ is constant; consequently, (97) will not satisfy (91), and the problem can, in general, only be solved numerically. However, a justification for (97) is given by Feng:^{18} namely that, typically, $ S \u226a 1 ;$ this will be reassessed later in Sec. IV when we compare this solution against two numerical solutions obtained via different methods.
C. $ P e \u226b 1$

the exponent for Re is $ 1 / 3 ,$ rather than 1/2;

the constant of proportionality is different;

the original correlation suggests that $ S h R M \u2192 2$ as $ R e \u2192 0$ if $ S c \u226a R e \u2212 3 / 2 ,$ but (147) gives that this should be $ S c \u226a R e \u2212 1$ instead;
 Eq. (147) is based on using the characteristic length scale as the radius of the sphere, rather than the diameter. This becomes evident when using Eq. (143) when $ P e \u226a 1 :$ in this case, the solution for c is given by$ c = 1 r ,$leading to $ S h = 1.$ In fact,$ S h R M ( P e diam ) = 2 S h ( P e ) ,$
IV. NUMERICAL IMPLEMENTATION
There remain two numerical tasks:
We describe each of these in turn.
A. Full equations
Equations (42)–(44), subject to (28)–(36), were solved using the 2D axial symmetry mode and transient solver in the finite elementbased software, Comsol Multiphysics. In combination with the software's deformed mesh mode, whereby an arbitrary Lagrangian–Eulerian (ALE) formulation is employed in order to solve free or movingboundary problems,^{38,39,43} Lagrangian P2–P1 elements were used for Eqs. (42) and (43), whereas Lagrangian secondorder elements were used for Eq. (44). With regard to the size of the computational domain, its outer edge was taken to be ten times that of the initial radius of the sphere; this was the extent used in several earlier works involving slow flow past a sphere,^{40–42} albeit without dissolution, and was found to be adequate here also. However, as dissolution proceeds, considerable mesh distortion was found to occur, leading to some numerical challenges. Comsol Multiphysics offers the options of Laplace smoothing or Winslow smoothing^{48} in these situations; however, it is wellestablished that the latter is more effective,^{20} although there were nonetheless some difficulties in being able to model complete dissolution. Even so, there were typically no difficulties until around 90% dissolution, and using further measures, such as remeshing, enabled us to compute up until 99% dissolution.
$ N el$ .  $ N dof$ .  $ N s$ . 

6571  56 134  80 
8397  74 655  120 
10 738  98 247  160 
12 572  117 538  200 
$ N el$ .  $ N dof$ .  $ N s$ . 

6571  56 134  80 
8397  74 655  120 
10 738  98 247  160 
12 572  117 538  200 
The accuracy of the numerical method was verified in a number of ways. First, we checked the method for large values of Pe, assuming that the interface does not move. Figure 2 shows $ \u2212 P e \u2212 1 / 3 \u2202 c / \u2202 r$ at r = 1 as a function of $ \pi \u2212 \theta $ for $ P e = 10 2 , 10 4 , 10 6 ,$ as determined numerically. We see that the code is able to adequately reproduce the analytical result given by Eq. (141) for high values of Pe.
As is evident, good agreement is obtained for all of the meshes, although nonetheless the mesh for which $ N s = 200$ was used for all subsequent computations.
B. Diffusiondominated problem
For this case, the onedimensional version of the arbitrary Lagrangian–Eulerian (ALE) formulation was used; this constitutes an approach that is independent of that explained in Sec. IV A to finding solutions for the $ P e \u226a 1$ regime, since the velocity field does not have to be solved for. Lagrangian secondorder elements were used to discretize Eq. (91) on a domain that extended to ten times the sphere radius and criterion (153) was again used; needless to say, far fewer elements and degrees of freedom were required, with meshindependent solutions being obtain for $ N dof$ on the order of several hundreds.
Figure 4 shows a comparison for $ f ( t )$ in the limit of low Péclet number. As seen in the figure, the agreement between the two approaches developed here is very good; for completeness, we have also included the corresponding result obtained using the method from Rice and Do,^{29} which clearly gives a consistent overestimation of the remaining fraction.
V. RESULTS
A. Effect of Re, Sc, and $C$
Figure 5 shows the extent of remaining particle for $ R e = 10 \u2212 3$ and $ P e = 10 \u2212 3 , 10 4 ,$ i.e., $ S c = 1 , 10 7$, after 0%, 30%, 60%, and 90% of the original volume has dissolved; for these results, $ c 0 / c sat = 290 ,$ corresponding to ibuprofen in 0.05 M of pH 6.6 phosphate buffer.^{24} For Fig. 5(b), it is evident that dissolution occurs increasingly nonspherically, with a cusplike form developing toward the rear of the particle. Moreover, considerably more dissolution occurs at the front end of the particle than at the rear: when 90% has dissolved, the particle has dissolved almost all of the way to its original center along $ \theta = \pi ,$ whereas very little dissolution has occurred along $ \theta = 0.$ This result is in line with Fig. 2.
B. Comparison with simplified models
First of all, note that the different approaches indicate different numbers of governing dimensionless parameters. In our approach, these are $ R e , S c , \rho $, and $ C ,$ with the list reducing to $ P e ( = ReSc ) , \rho $, and $ C ,$ when $ R e \u226a 1.$ On the other hand, the approach using correlations appears to oversimplify the parametric dependence: as evident from Eq. (C5), the nondimensional problem will not depend on ρ and $ c 0 / c sat$ at all, but only Re and Sc as conveyed through the Sherwood number. These comments notwithstanding, Fig. 9 shows the time taken for 90% of the particle to dissolve as a function of Pe for $ R e = 10 \u2212 3 ,$ ρ = 1, $ c 0 / c sat = 2.9 , 290 , 2900.$ Although the agreement between the different approaches appears reasonable on this log –log plot as a whole, this is not the case for low to intermediate values of Pe. This is seen even more clearly in Table II, which shows the actual time, in seconds, taken to dissolve 90% of the initial load for different values of Pe, as computed via the different methods and for the same values of $ R e , \rho $, and $C$ as in Fig. 9. Several features are notable here:

for high values of Pe, the Ranz–Marshall correlation overpredicts the dissolution time by as much as a factor of four, which may be deemed very inaccurate in pharmaceutical contexts, in particular;

the correlation due to Clift et al.^{11} works well for all values of Pe, even though it was originally suggested to be of use for Pe only as high as 70;

the correlation due to Acrivos and Taylor^{3} works best at the highest value of Pe, overestimating the dissolution time only slightly. Note also that we have not given any value when $ P e \u226a 1 ,$ as the correlation was derived based on the assumption that $ P e \u226b 1.$ On the other hand, one might ask why the correlation is still accurate even though the initial particle decreases in size during dissolution, meaning effectively lower values of Pe. In fact, this is because a 90% reduction in particle volume results in approximately a halving of the particle radius, meaning a still high value of Pe. Nevertheless, were we to compute the time taken for complete dissolution, we might expect this correlation to be inaccurate. This is further explored in Appendix C.
.  Pe .  

Method .  10^{−3} .  10^{−1} .  10^{2} .  10^{4} . 
Ranz and Marshall^{28}  4.40 $ \xd7 10 \u2212 3$  4.27 $ \xd7 10 \u2212 1$  2.93 $ \xd7 10 2$  1.30 $ \xd7 10 4$ 
Clift et al.^{11}  3.90 $ \xd7 10 \u2212 3$  3.78 $ \xd7 10 \u2212 1$  1.23 $ \xd7 10 2$  3.03 $ \xd7 10 3$ 
Acrivos and Taylor^{3}  ⋯  ⋯  1.47 $ \xd7 10 2$  3.20 $ \xd7 10 3$ 
Full numerical  3.70 $ \xd7 10 \u2212 3$  3.63 $ \xd7 10 \u2212 1$  1.23 $ \xd7 10 2$  3.05 $ \xd7 10 3$ 
.  Pe .  

Method .  10^{−3} .  10^{−1} .  10^{2} .  10^{4} . 
Ranz and Marshall^{28}  4.40 $ \xd7 10 \u2212 3$  4.27 $ \xd7 10 \u2212 1$  2.93 $ \xd7 10 2$  1.30 $ \xd7 10 4$ 
Clift et al.^{11}  3.90 $ \xd7 10 \u2212 3$  3.78 $ \xd7 10 \u2212 1$  1.23 $ \xd7 10 2$  3.03 $ \xd7 10 3$ 
Acrivos and Taylor^{3}  ⋯  ⋯  1.47 $ \xd7 10 2$  3.20 $ \xd7 10 3$ 
Full numerical  3.70 $ \xd7 10 \u2212 3$  3.63 $ \xd7 10 \u2212 1$  1.23 $ \xd7 10 2$  3.05 $ \xd7 10 3$ 
VI. CONCLUSION
Furthermore, although this work has focused primarily on dissolution, it is worth taking stock of the underlying flow physics. Essentially, we have considered a quasisteady slow flow past a bluff body that it is initially spherical, but whose shape changes with time. Although solute concentration does not enter the equations for the velocity field, the latter does enter the equation for the former; thence, dissolution alters the shape of the bluff body, which in turns affects the velocity field. Thus, although slow flow is by now a wellresearched topic in fluid mechanics, we are unaware of any other work that has explored this particular coupling between fluid mechanics, solute concentration and dissolution.
Numerous directions are now possible for extending this work. One would be to consider how the results differ when natural convection dominates; this could occur as a consequence of dissolution altering the local density of the surrounding solvent, as indicated in Appendix B. Another direction would be to consider what happens for higher values of Re, which will introduce an even greater coupling between the fluid mechanics and dissolution than when $ R e \u226a 1.$ A further extension would be to consider a particle that is allowed to move in the solvent, requiring an additional forcebalance condition to be imposed on the particle; although some attempts to do this already do exist,^{15,16,24} this has only been done in the context of the Ranz–Marshall correlation. Still further field would be to consider dissolution from porous granules; this would require the development of a diffusioncontrolled dissolution model^{22,23} that is coupled to the flow of a solvent through a porous body.^{19,27,33}
ACKNOWLEDGMENTS
This publication has emanated from research supported in part by a grant from Science Foundation Ireland cofunded under the European Regional Development Fund under Grant No. 12/RC/2275_P2. For the purpose of Open Access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission. Also, the authors would like to thank the anonymous reviewers for their important comments and suggestions.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Milton Assunção: Conceptualization (supporting); Investigation (lead); Software (lead); Validation (lead); Writing – original draft (supporting); Writing – review & editing (supporting). Michael Vynnycky: Conceptualization (equal); Formal analysis (lead); Project administration (supporting); Writing – original draft (lead); Writing – review & editing (lead). Kevin Moroney: Conceptualization (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX A: $ U \u226a \rho sat B g a 2 / \mu $
If (41) is not satisfied, e.g., in situations where the solvent is assumed to be stagnant, a different nondimensionalization is necessary.
APPENDIX B: SOLUTION TO EQUATION (98) WHEN $ S > 2 \pi $
APPENDIX C: SIMPLIFIED MODELS
A simplified approach is to assume that the particle remains spherical during dissolution, which is proportional to a rate, k, that is itself proportional to the Sherwood number that would prevail under steadystate conditions;^{15,16,24} this allows the use of wellknown correlations for the Sherwood number as a function of the Reynolds and Schmidt numbers, most notably the Ranz–Marshall relation,^{28} although others are available also.^{11} Mathematically, this approach leads to a firstorder ordinary differential equation for the particle radius, $ r d ,$ as a function of time, as follows.
Figure 10 compares the results of Eqs. (C10), (C18), and (C21); note, however, that we have only plotted results using Eq. (C21) for Pe as low as 10, as it is strictly speaking only valid for $ P e \u226b 1.$ Note also that this figure can also be compared directly with three of the curves in Fig. 9. The differences, if any, are very small, indicating, that computing the time to 90% dissolution gives a very good approximation to the total dissolution time.