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 finite-element software Comsol Multiphysics, that is able to solve the resulting axisymmetric moving-boundary 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 convection-dominated regime for which it was not originally intended.
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 well-established 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 Whitney26 postulated that the dissolution rate is linearly proportional to the difference between the particle's solubility, and the concentration of dissolved solid in the bulk solution, Ct. Soon after, this evolved into the well-known Nernst–Brunner equation by treating dissolution as a diffusion-controlled 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 three-dimensional (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 one-dimensional (1D) conditions. This is because only in 1D can the concentration gradient across the diffusion layer be a constant value, denoted as . 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 non-constant, distance-dependent 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 steady-state diffusion only, whereas Fick's second law covers both steady and non-steady 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 diffusion-dominated 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 moving-boundary 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 finite-element 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 two-pronged 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 diffusion-dominated 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 diffusion-dominated regime, providing a previously overlooked analytical solution in the limit of slow dissolution; and Appendix C discusses simplified models for forced-convection-dominated 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 x-direction; the situation can, therefore, be assumed to be axisymmetric. For 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 coordinates. Denoting the resulting surface by the proposed evolution is as shown in Fig. 1. Let c0 be the initial solid concentration in the sphere, where we shall suppose that with csat 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 We will also assume that for at where even though we will ultimately take we have allowed for the possibility that 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, 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 where M is the known relative molecular mass of the particle material. However, although ρs and c0 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.
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).
The analysis would now proceed differently, depending on whether and
For the case when we need simply to set ν = 0 in (59).
However, strictly speaking, Eq. (97) is erroneous, as it was obtained by assuming that 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, this will be reassessed later in Sec. IV when we compare this solution against two numerical solutions obtained via different methods.
the exponent for Re is rather than 1/2;
the constant of proportionality is different;
the original correlation suggests that as if but (147) gives that this should be 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 in this case, the solution for c is given byleading to In fact,where Moreover, (149) gives
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 element-based 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 moving-boundary problems,38,39,43 Lagrangian P2–P1 elements were used for Eqs. (42) and (43), whereas Lagrangian second-order 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 smoothing48 in these situations; however, it is well-established 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 re-meshing, enabled us to compute up until 99% dissolution.
|10 738||98 247||160|
|12 572||117 538||200|
|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 at r = 1 as a function of for 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 was used for all subsequent computations.
B. Diffusion-dominated problem
For this case, the one-dimensional 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 regime, since the velocity field does not have to be solved for. Lagrangian second-order 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 mesh-independent solutions being obtain for on the order of several hundreds.
Figure 4 shows a comparison for 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.
A. Effect of Re, Sc, and
Figure 5 shows the extent of remaining particle for and i.e., , after 0%, 30%, 60%, and 90% of the original volume has dissolved; for these results, 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 non-spherically, with a cusp-like 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 whereas very little dissolution has occurred along 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 , and with the list reducing to , and when 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 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 ρ = 1, 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 , and 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 Taylor3 works best at the highest value of Pe, overestimating the dissolution time only slightly. Note also that we have not given any value when as the correlation was derived based on the assumption that 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.
|Method .||10−3 .||10−1 .||102 .||104 .|
|Ranz and Marshall28||4.40||4.27||2.93||1.30|
|Clift et al.11||3.90||3.78||1.23||3.03|
|Acrivos and Taylor3||⋯||⋯||1.47||3.20|
|Method .||10−3 .||10−1 .||102 .||104 .|
|Ranz and Marshall28||4.40||4.27||2.93||1.30|
|Clift et al.11||3.90||3.78||1.23||3.03|
|Acrivos and Taylor3||⋯||⋯||1.47||3.20|
Furthermore, although this work has focused primarily on dissolution, it is worth taking stock of the underlying flow physics. Essentially, we have considered a quasi-steady 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 well-researched 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 A further extension would be to consider a particle that is allowed to move in the solvent, requiring an additional force-balance 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 diffusion-controlled dissolution model22,23 that is coupled to the flow of a solvent through a porous body.19,27,33
This publication has emanated from research supported in part by a grant from Science Foundation Ireland co-funded 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.
Conflict of Interest
The authors have no conflicts to disclose.
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).
The data that support the findings of this study are available from the corresponding author upon reasonable request.
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
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 steady-state conditions;15,16,24 this allows the use of well-known 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 first-order ordinary differential equation for the particle radius, 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 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.