On the dissolution of a solid spherical particle

Assunção, M.; Vynnycky, M.; Moroney, K. M.

doi:10.1063/5.0144883

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.

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;⁹ 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 Whitney²⁶ 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 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.³²

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,³⁴ 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,²⁸ 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,²⁸ 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 $(C_{s} - 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 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.⁴⁵ 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.³² 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,⁴⁵ 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.⁹

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 $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, θ, ϕ)$ coordinates. Denoting the resulting surface by $r = R (θ, t),$ the proposed evolution is as shown in Fig. 1. Let c₀ 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, θ, t) .$ We will also assume that $c = c_{\infty} (\geq 0)$ for $r > R (θ, t)$ at $t = 0,$ where $c_{\infty} < c_{sat};$ even though we will ultimately take $c_{\infty} = 0,$ we have allowed for the possibility that $c_{\infty} > 0,$ as this is done elsewhere,^8,18,29 although this will have no impact on the analysis.

FIG. 1.

View large Download slide

Schematic of axisymmetric solid particle dissolution: (a) at t = 0, the particle is a sphere of radius a; (b) at a later time, it has partially dissolved and occupies $r \leq R (θ, t)$ ⁠, where $R (θ, t) < a$ ⁠.

The governing equations for conservation of mass, momentum, and solute concentration for the solvent are given by

\frac{\partial ρ_{l}}{\partial t} + \nabla \cdot (ρ_{l} u) = 0,

(1)

ρ_{l} (\frac{\partial u}{\partial t} + (u \cdot \nabla) u) = - \nabla p + μ \nabla^{2} u + ρ_{l} g,

(2)

\frac{\partial c}{\partial t} + \nabla \cdot (c u) = D \nabla^{2} c,

(3)

respectively, where u is the velocity of the solvent, p is the pressure, ρ_l is the density of the solvent, μ is its dynamic viscosity, D is the diffusion coefficient for the solute in the solvent, which we will take to be constant for simplicity, and g is the gravity vector; for later use, we will write

g = - g k,

where k is the unit vector in the upward vertical direction and g is the acceleration due to gravity. Note that although we have marked the direction of gravity in Fig. 1 as being perpendicular to the flow direction of the solvent at infinity, this will be immaterial in this paper, as its effect will ultimately be neglected; nevertheless, it is necessary to include it for the purposes of discussing the formulation of the model equations.

We must now set boundary conditions. As

r \to \infty

⁠, we take

u_{r} \to U cos θ, u_{θ} \to - U sin θ, c \to c_{\infty},

(4)

where u_r and

u_{θ}

are the r- and θ-components of

u .

At

r = R (θ, t)

⁠, we will suppose that

c = c_{sat},

(5)

which assumes that the dissolution (reaction) kinetics at the interface are fast compared to the mass transfer from the interface; this is a commonly used assumption.^{8,18,21,29,32} The other conditions will be obtained by considering mass and solute conservation for the system as a whole; although this is often done by referring to early work by Scriven³⁰ and the first edition of the textbook by Bird et al.,⁵ here we adopt a more fundamental approach, as it obtains the same result, while highlighting the significance of often-neglected simplifications. Thus, for mass conservation, we have

\frac{4}{3} π a^{3} ρ_{s} + 2 π \int_{0}^{π} \int_{a}^{\infty} ρ_{l} r^{2} sin θ d r d θ = 2 π ρ_{s} \int_{0}^{π} \int_{0}^{R (θ, t)} r^{2} sin θ d r d θ + 2 π \int_{0}^{π} \int_{R (θ, t)}^{\infty} ρ_{l} r^{2} sin θ d r d θ,

(6)

where ρ_s is the density of the particle, which we will treat as constant. Differentiating with respect to t and applying the Leibniz rule, we have

ρ_{s} \int_{0}^{π} \frac{\partial R}{\partial t} R^{2} (θ, t) sin θ d θ = ρ_{l} \int_{0}^{π} \frac{\partial R}{\partial t} R^{2} (θ, t) sin θ d θ - \int_{0}^{π} \int_{R (θ, t)}^{\infty} \frac{\partial ρ_{l}}{\partial t} r^{2} sin θ d r d θ,

(7)

using Eq. and applying the divergence theorem to the remaining volume integral, we have

ρ_{s} \int_{0}^{π} \frac{\partial R}{\partial t} R^{2} (θ, t) sin θ d θ = ρ_{l} \int_{0}^{π} \frac{\partial R}{\partial t} R^{2} (θ, t) sin θ d θ + \int_{0}^{π} {ρ_{l} u \cdot n |}_{r = R (θ, t)} R^{2} (θ, t) sin θ d θ,

(8)

where n is the unit normal pointing out of the solvent domain. Applying on a pointwise basis at

r = R (θ, t)

⁠, we have

u \cdot n = (\frac{ρ_{s}}{ρ_{l}} - 1) \frac{\partial R}{\partial t} at r = R (θ, t) .

(9)

Note that an issue arises if we assume, as is often done,^18,21,30 that ρ_l is constant. In this case, gives that

ρ_{s} = ρ_{l};

however, it is often the case that

ρ_{s} \neq ρ_{l} .

Hence, a self-consistent approach is to assume that ρ_l is not constant, which then allows to be derived correctly, noting that ρ_l in is not a constant value, but rather just the value of ρ_l at

r = R (θ, t) .

Having established that ρ_l cannot be constant throughout the solvent region, it is then necessary to invoke an appropriate constitutive relation for it. In the context of drug dissolution and many others, a plausible one is a linear dependence of ρ_l on the solute concentration, i.e.,

ρ_{l} = ρ_{ref} (1 + β (c - c_{ref})),

(10)

where ρ_ref is reference solvent density, c_ref is a reference solute concentration, and β is the solutal expansion coefficient.

As regards solute conservation for the system, this is given by

\frac{4}{3} π a^{3} c_{0} + 2 π \int_{0}^{π} \int_{a}^{\infty} c_{\infty} r^{2} sin θ d r d θ = 2 π c_{0} \int_{0}^{π} \int_{0}^{R (θ, t)} r^{2} sin θ d r d θ + 2 π \int_{0}^{π} \int_{R (θ, t)}^{\infty} c r^{2} sin θ d r d θ .

(11)

Differentiating with respect to t and applying the Leibniz rule again, we have

- 2 π c_{0} \int_{0}^{π} \frac{\partial R}{\partial t} R^{2} sin θ d θ = 2 π (\int_{0}^{π} (\int_{R (θ, t)}^{\infty} \frac{\partial c}{\partial t} r^{2} sin θ d r - c_{r = R (θ, t)} \frac{\partial R}{\partial t} R^{2} sin θ) d θ),

(12)

and, in view of , we have

\int_{0}^{π} \int_{R (θ, t)}^{\infty} \frac{\partial c}{\partial t} r^{2} sin θ d r d θ = (c_{sat} - c_{0}) \int_{0}^{π} \frac{\partial R}{\partial t} R^{2} sin θ d θ .

(13)

Using and applying the divergence theorem on , we will end up with

\int_{0}^{π} {(D \nabla c - c_{sat} u) \cdot n} R^{2} (θ, t) sin θ d θ = (c_{sat} - c_{0}) \int_{0}^{π} \frac{\partial R}{\partial t} R^{2} (θ, t) sin θ d θ,

(14)

which, when also applied on a pointwise basis at

r = R (θ, t)

⁠, reduces to

(D \nabla c - c_{sat} u) \cdot n = (c_{sat} - c_{0}) \frac{\partial R}{\partial t} at r = R (θ, t),

(15)

since the integral of

D \nabla c - c u

at infinity will be negligibly small. Thence, gives, on using ,

D \nabla c \cdot n = {c_{sat} (\frac{ρ_{s}}{ρ_{l}} - 1) + (c_{sat} - c_{0})} \frac{\partial R}{\partial t} at r = R (θ, t),

(16)

i.e.,

D \nabla c \cdot n = {c_{sat} (\frac{ρ_{s}}{ρ_{l}}) - c_{0}} \frac{\partial R}{\partial t} at r = R (θ, t) .

(17)

Note that is the more general form of the condition used elsewhere⁴ and accounts for the possibility that the particle does not remain spherical. In addition, we will have no slip at the surface of the solid particle; hence,

u \cdot t = 0,

(18)

where t is the unit tangential vector at the surface.

The problem also requires initial conditions, and for these, we simply take, for

r > a,

c = c_{\infty},

(19)

u = u_{0},

(20)

where

u_{0}

is a prescribed initial velocity field, and

R (θ, 0) = a .

(21)

A sensible choice for

u_{0}

would be the steady solution to the Navier–Stokes equations for flow past a sphere, which satisfies the first two conditions in , as well as and ; indeed, this is what we shall take later.

A final detail that we wish to note is that, in the formulation above, $ρ_{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 $ρ_{s} = c_{0} M,$ where M is the known relative molecular mass of the particle material. However, although ρ_s and c₀ 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

Taking c_ref = c_sat and

ρ_{ref} = ρ_{sat},

we nondimensionalize through

\begin{matrix} r^{'} = \frac{r}{a}, R^{'} = \frac{R}{a}, t^{'} = \frac{t}{a / U}, c^{'} = \frac{c - c_{\infty}}{c_{sat} - c_{\infty}}, \\ u^{'} = \frac{u}{U}, u_{0}^{'} = \frac{u_{0}}{U}, p^{'} = \frac{p + ρ_{sat} g (1 - β (c_{sat} - c_{\infty})) z}{ρ_{sat} (1 - β (c_{sat} - c_{\infty})) U^{2}}, \end{matrix}

(22)

where z is the spatial coordinate in the opposite direction to g in Fig. 1. On dropping the primes, we will have

\frac{\partial}{\partial t} ((1 + B (c - 1))) + \nabla \cdot ((1 + B (c - 1)) u) = 0,

(23)

(1 + B (c - 1)) (\frac{\partial u}{\partial t} + (u \cdot \nabla) u) = - \nabla p + \frac{1}{R e} \nabla^{2} u - \frac{B}{(1 - B) F r^{2}} c k,

(24)

\frac{\partial c}{\partial t} + \nabla \cdot (c u) = \frac{1}{P e} \nabla^{2} c,

(25)

where Re, Pe, and Fr are the Reynolds, Péclet, and Froude numbers, which are given by

R e = \frac{ρ_{sat} (1 - B) U a}{μ}, P e = \frac{U a}{D}, F r = \frac{U}{\sqrt{g a}},

(26)

respectively. Also,

P e = ReSc,

where Sc denotes the Schmidt number and is given by

S c = \frac{μ}{ρ_{sat} (1 - B) D},

(27)

and

B = β (c_{sat} - c_{\infty});

this does not have a commonly used name, although it quantifies the maximum density variation due to concentration differences relative to the reference density. Equations will be subject to the following boundary conditions: at

r = R (θ, t),

u \cdot t = 0,

(28)

u \cdot n = (ρ - 1) \frac{\partial R}{\partial t},

(29)

\frac{1}{P e} \nabla c \cdot n = - C \frac{\partial R}{\partial t},

(30)

c = 1,

(31)

where

ρ = \frac{ρ_{s}}{ρ_{sat}}, C = (\frac{c_{0} - ρ c_{sat}}{c_{sat} - c_{\infty}}),

(32)

as

r \to \infty

⁠,

u_{r} \to cos θ, u_{θ} \to - sin θ, c \to 0.

(33)

The initial conditions now become, at

t = 0,

c = 0, r > 1, 0 \leq θ \leq π,

(34)

u = u_{0},

(35)

R (θ, 0) = 1, 0 \leq θ \leq π .

(36)

For pharmaceutical applications, the usual regime of interest is for

B ≪ 1, R e ≪ 1.

The second of these inequalities describes the well-known slow flow regime,¹ wherein viscous forces dominate inertia. In this case, it is appropriate to rescale p via

p = P / R e,

leading to

\nabla \cdot u = 0,

(37)

R e (\frac{\partial u}{\partial t} + (u \cdot \nabla) u) = - \nabla P + \nabla^{2} u - (\frac{BRe}{F r^{2}}) c k,

(38)

\frac{\partial c}{\partial t} + \nabla \cdot (c u) = \frac{1}{P e} \nabla^{2} c,

(39)

subject to –, where now

R e = \frac{ρ_{sat} U a}{μ}, S c = \frac{μ}{ρ_{sat} D};

(40)

in fact, the simplification associated with

B ≪ 1

is the well-known Boussinesq approximation and effectively means that the density of the solvent does not change appreciably from its initial reference value, even though material is dissolving into it. One might also expect

BRe / F r^{2} ≪ 1,

i.e.,

U ≫ ρ_{sat} B g a^{2} / μ,

(41)

giving just

\nabla \cdot u = 0,

(42)

R e (\frac{\partial u}{\partial t} + (u \cdot \nabla) u) = - \nabla P + \nabla^{2} u,

(43)

\frac{\partial c}{\partial t} + \nabla \cdot (c u) = \frac{1}{P e} \nabla^{2} c .

(44)

Note, however, that precludes the possibility of using the above formulation directly for the case of a stagnant fluid at infinity, i.e.,

U = 0,

apart from for the particular case when

ρ_{s} = ρ_{l};

the resulting analysis for when

U ≪ ρ_{sat} B g a^{2} / μ

and

ρ_{s} \neq ρ_{l}

is given in Appendix A.

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 quasi-analytically for $P e ≪ 1$ and $P e ≫ 1.$ On the other hand, when $P e \sim O (1),$ no such simplified analysis is possible; consequently, we do not consider it, other than solely numerically in Sec. V.

B. $P e ≪ 1$

In this case, which corresponds to the dominance of diffusion over convection, Eq. would reduce

\nabla^{2} c \approx 0,

meaning that c would not be able to satisfy the initial condition; thus, we need to rescale t according to

t = PeT .

Doing so, we have

\nabla \cdot u = 0,

(45)

R e (\frac{1}{P e} \frac{\partial u}{\partial T} + (u \cdot \nabla) u) = - \nabla P + \nabla^{2} u,

(46)

\frac{\partial c}{\partial T} + P e \nabla \cdot (c u) = \nabla^{2} c,

(47)

subject to the following: at

r = R (θ, T),

u \cdot t = 0,

(48)

u \cdot n = (\frac{ρ - 1}{P e}) \frac{\partial R}{\partial T},

(49)

\nabla c \cdot n = - C \frac{\partial R}{\partial T},

(50)

c = 1;

(51)

as

r \to \infty,

u_{r} \to cos θ, u_{θ} \to - sin θ, c \to 0;

(52)

at

T = 0,

c = 0, r > 1, 0 \leq θ \leq π,

(53)

R = 1, 0 \leq θ \leq π .

(54)

The analysis would now proceed differently, depending on whether $| ρ - 1 | / P e ≲ 1$ and $| ρ - 1 | / P e ≫ 1.$

1. $| ρ - 1 | / P e ≲ 1$

Considering the distinguished limit

| ρ - 1 | / P e \sim 1

⁠, and setting

ν = | ρ - 1 | / P e,

an

O (1)

constant, we have, at leading order,

\nabla \cdot u = 0,

(55)

\frac{1}{S c} \frac{\partial u}{\partial T} = - \nabla P + \nabla^{2} u,

(56)

\frac{\partial c}{\partial T} = \nabla^{2} c,

(57)

subject to the following: at

r = R (θ, T),

u \cdot t = 0,

(58)

u \cdot n = ν \frac{\partial R}{\partial T},

(59)

\nabla c \cdot n = - C \frac{\partial R}{\partial T},

(60)

c = 1;

(61)

as

r \to \infty,

u_{r} \to cos θ, u_{θ} \to - sin θ, c \to 0;

(62)

at

T = 0,

c = 0, r > 1, 0 \leq θ \leq π,

(63)

R = 1, 0 \leq θ \leq π .

(64)

For the case when $| ρ - 1 | / P e ≪ 1,$ we need simply to set ν = 0 in (59).

2. $| ρ - 1 | / P e ≫ 1$

In this case, we would need to rescale u and P. Setting

U = \frac{u}{(ρ - 1) / P e}, \bar{P} = \frac{P}{(ρ - 1) / P e},

(65)

we have, at leading order,

\nabla \cdot U = 0,

(66)

\frac{1}{S c} (\frac{\partial U}{\partial T} + (ρ - 1) (U \cdot \nabla) U) = - \nabla \bar{P} + \nabla^{2} U,

(67)

\frac{\partial c}{\partial T} + (ρ - 1) \nabla \cdot (c U) = \nabla^{2} c,

(68)

subject to the following: at

r = R (θ, T),

U \cdot t = 0,

(69)

U \cdot n = \frac{\partial R}{\partial T},

(70)

\nabla c \cdot n = - C \frac{\partial R}{\partial T},

(71)

c = 1;

(72)

as

r \to \infty,

U_{r} \to 0, U_{θ} \to 0, c \to 0;

(73)

at

T = 0,

c = 0, r > 1, 0 \leq θ \leq π,

(74)

R = 1, 0 \leq θ \leq π .

(75)

Now, if

| ρ - 1 | ≪ 1,

there is once again decoupling between the problems for c and U. When

| ρ - 1 | \sim O (1),

there appears to be a coupled spherically symmetric problem in which

U, \bar{P}

⁠, and c have to be solved for simultaneously. Anticipating a solution of the form

U_{r} = U_{r} (r, T), U_{θ} = 0, R = R (T) .

(76)

Equations give, respectively,

\frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} U_{r}) = 0,

(77)

\frac{1}{S c} (\frac{\partial U_{r}}{\partial T} + (ρ - 1) U_{r} \frac{\partial U_{r}}{\partial r}) = - \frac{\partial \bar{P}}{\partial r} + \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial U_{r}}{\partial r}) - \frac{2 U_{r}}{r^{2}},

(78)

\frac{\partial c}{\partial T} + (ρ - 1) U_{r} \frac{\partial c}{\partial r} = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial c}{\partial r}),

(79)

subject to the following: at

r = R (T),

U_{r} = \dot{R},

(80)

\frac{\partial c}{\partial r} = - C \dot{R},

(81)

c = 1,

(82)

where dots denote differentiation with respect to

T;

as

r \to \infty,

U_{r} \to 0, c \to 0;

(83)

at

T = 0,

c = 0, r > 1,

(84)

R = 1.

(85)

Note that Eqs. are obtained from –, whereas is satisfied automatically by virtue of . Equations are a more complete form of the equations considered in earlier work,^8,18,21 wherein the r-component of the momentum equation was not mentioned at all. Then, from , , and , we have

U_{r} = R^{2} \dot{R} / r^{2},

(86)

leading to

\frac{1}{S c} (\frac{R^{2} \ddot{R}}{r^{2}} - 2 (ρ - 1) \frac{R^{4} {\dot{R}}^{2}}{r^{5}}) = - \frac{\partial \bar{P}}{\partial r},

(87)

\frac{\partial c}{\partial T} + (ρ - 1) \frac{R^{2} \dot{R}}{r^{2}} \frac{\partial c}{\partial r} = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial c}{\partial r});

(88)

thus,

\bar{P} = \frac{1}{S c} (\frac{R^{2} \ddot{R}}{r} - \frac{1}{2} (ρ - 1) \frac{R^{4} {\dot{R}}^{2}}{r^{4}}) + P_{\infty},

(89)

where

P_{\infty}

is the value of

\bar{P},

possibly time-dependent, at infinity.

The leading-order problem for c will now be spherically symmetric, noting that we have the stagnant fluid case. Setting

\begin{matrix} c = c^{(0)} + P e c^{(1)} + O (P e^{2}), \\ R = R^{(0)} + P e R^{(1)} + O (P e^{2}), \end{matrix}

(90)

we have, at

O (1),

\frac{\partial c^{(0)}}{\partial T} + (ρ - 1) \frac{R^{{(0)}^{2}} {\dot{R}}^{(0)}}{r^{2}} \frac{\partial c^{(0)}}{\partial r} = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial c^{(0)}}{\partial r}),

(91)

subject to the following: at

r = R^{(0)} (T),

\frac{\partial c^{(0)}}{\partial r} = - C {\dot{R}}^{(0)},

(92)

c^{(0)} = 1;

(93)

as

r \to \infty,

c^{(0)} \to 0;

(94)

at

T = 0,

c^{(0)} = 0, r > 1,

(95)

R^{(0)} = 1.

(96)

In fact, the problem given by Eqs. is equivalent to the one considered by Rice and Do,²⁹ who claimed to find an analytical solution for

c^{(0)}

when

ρ = 1 :

in our notation, their solution would take the form

c^{(0)} = \frac{R^{(0)}}{r} erfc (\frac{r - R^{(0)}}{2 T^{1 / 2}}),

(97)

where

R^{(0)}

satisfies

{\dot{R}}^{(0)} = - S (\frac{1}{π^{1 / 2} T^{1 / 2}} + \frac{1}{R^{(0)}})

(98)

subject to , where

S = (c_{sat} - c_{\infty}) / (c_{0} - c_{sat}) .

Setting

R^{(0)} = T^{1 / 2} V,

Eq. becomes

T \dot{V} = - S (\frac{1}{π^{1 / 2}} + \frac{1}{V}) - \frac{1}{2} V,

(99)

leading to

\int^{T} \frac{d T}{T} = - 2 \int^{V} \frac{V d V}{{(V + \frac{S}{π^{1 / 2}})}^{2} + 2 S - \frac{S^{2}}{π}} .

(100)

In fact, Rice and Do²⁹ only gave the solution for when

2 S - S^{2} / π > 0,

i.e.,

0 < S < 2 π

⁠; however, there is another form when

2 S - S^{2} / π < 0,

i.e.,

S > 2 π .

First, for

0 < S < 2 π,

setting

λ^{2} = 2 S - S^{2} / π > 0,

so that

\int^{T} \frac{d T}{T} = - 2 \int^{V} \frac{V d V}{{(V + \frac{S}{π^{1 / 2}})}^{2} + λ^{2}},

(101)

we put

λ W = V + S / π^{1 / 2},

which leads, on integration, to

T (1 + \frac{1}{λ^{2}} {(\frac{R^{(0)}}{T^{1 / 2}} + \frac{S}{π^{1 / 2}})}^{2}) = C exp (\frac{2 S}{λ π^{1 / 2}} {tan}^{- 1} (\frac{1}{λ} (\frac{R^{(0)}}{T^{1 / 2}} + \frac{S}{π^{1 / 2}}))),

(102)

where C is a constant to be determined. Since

R^{(0)} (0) = 1,

we must have

C = \frac{exp (- S π^{1 / 2} / λ)}{λ^{2}},

(103)

whence

λ^{2} T + {(R^{(0)} + \frac{S T^{1 / 2}}{π^{1 / 2}})}^{2} = exp (\frac{2 S}{λ π^{1 / 2}} [{tan}^{- 1} (\frac{1}{λ} (\frac{R^{(0)}}{T^{1 / 2}} + \frac{S}{π^{1 / 2}})) - \frac{π}{2}]),

(104)

which corresponds to Rice and Do's equation ;²⁹ the corresponding result when

S > 2 π

is given in Appendix B.

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:¹⁸ namely that, typically, $S ≪ 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 ≫ 1$

In this limit, which corresponds to the dominance of convection over diffusion, we may suspect that reduces to

\frac{\partial R}{\partial t} \approx 0,

(105)

and hence

R \equiv 1

at leading order for all time. However, this assumes that

\nabla c \cdot n ≪ P e,

which might not necessarily be the case. To check this, first note that the leading-order problem for u would be the same as the steady-state problem. However, the leading-order problem for c is not the steady-state problem; at first sight, it appears to be

\frac{\partial c}{\partial t} + \nabla \cdot (c u) = 0

(106)

subject to . However, since the partial differential equation (PDE) in is of lower order than in , it will now not be possible to satisfy all of the initial and boundary conditions; thus, there would have to be a boundary layer near

R = 1.

In fact, it is well-known for the case of a sphere,^2,3 and we shall recap it shortly, that the steady-state problem, consisting of but without the time derivative and subject to and , would have a boundary layer of width

P e^{- 1 / 3},

and thence

\nabla c \cdot n \sim P e^{1 / 3},

but it seems at first sight that this cannot be the case here. If it were, the convection and diffusion terms in Eq. would be of

O (P e^{- 1 / 3})

in the boundary layer, meaning that would there reduce to

\frac{\partial c}{\partial t} \approx 0.

(107)

Anticipating this difficulty, we should first rescale t via

t = P e^{1 / 3} T

⁠. Now, Eq. becomes

\frac{1}{P e^{1 / 3}} \frac{\partial c}{\partial T} + \nabla \cdot (c u) = \frac{1}{P e} \nabla^{2} c,

(108)

and, at

r = R (θ, T),

Eq. becomes

\frac{1}{P e^{2 / 3}} \nabla c \cdot n = - C \frac{\partial R}{\partial T};

(109)

the other governing equations, boundary and initial conditions remain as before. Observe also that, since

\nabla c \cdot n \sim P e^{1 / 3}

at

r = R (θ, T),

will imply that

R \approx 1

at leading order and that the leading-order correction to R will be of

O (P e^{- 1 / 3});

in fact, this can be determined once the leading-order problem for c has been solved.

In the bulk, we now have

\nabla \cdot (c u) \approx 0;

(110)

then, in view of , we will have

u \cdot \nabla c \approx 0.

(111)

Hence, c will be constant on streamlines and, since all of these emanate from upstream where

c \to 0,

we will have that

c \approx 0

in the bulk.

Prior to considering the boundary layer, it is worth writing out the governing equations in full in

(r, θ, ϕ)

-coordinates under axisymmetric conditions. First, it is convenient to introduce the Stokes streamfunction,

ψ,

which satisfies

u_{r} = - \frac{1}{r^{2} sin θ} \frac{\partial ψ}{\partial θ}, u_{θ} = \frac{1}{r sin θ} \frac{\partial ψ}{\partial r} .

(112)

Thus, will be automatically satisfied, whereas taking the

curl

of Eq. will lead to

{(\frac{\partial^{2}}{\partial r^{2}} + \frac{sin θ}{r^{2}} \frac{\partial}{\partial θ} (\frac{1}{sin θ} \frac{\partial}{\partial θ}))}^{2} ψ = 0

(113)

with Eq. becoming

\frac{1}{P e^{1 / 3}} \frac{\partial c}{\partial T} + \frac{1}{r^{2} sin θ} (\frac{\partial ψ}{\partial r} \frac{\partial c}{\partial θ} - \frac{\partial ψ}{\partial θ} \frac{\partial c}{\partial r}) = \frac{1}{P e} {\frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial c}{\partial r}) + \frac{1}{r^{2} sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial c}{\partial θ})} .

(114)

Since the boundary conditions for ψ are

ψ = 0, \frac{\partial ψ}{\partial r} = 0 at r = 1,

(115)

ψ \sim - \frac{1}{2} r^{2} sin θ as r \to \infty,

(116)

we arrive at

ψ = \frac{1}{4} (3 r - \frac{1}{r} - 2 r^{2}) {sin}^{2} θ .

(117)

For the boundary layer, setting

r = 1 + P e^{- 1 / 3} \bar{r}, ψ = P e^{- 2 / 3} \bar{ψ},

(118)

Equation reduces to

\frac{\partial c}{\partial T} + \frac{\partial \bar{ψ}}{\partial \bar{r}} \frac{\partial c}{\partial s} - \frac{\partial \bar{ψ}}{\partial s} \frac{\partial c}{\partial \bar{r}} = \frac{\partial^{2} c}{\partial {\bar{r}}^{2}},

(119)

i.e.,

\frac{\partial c}{\partial T} + \frac{3}{2} \bar{r} s (2 - s) \frac{\partial c}{\partial s} - \frac{3}{2} {\bar{r}}^{2} (1 - s) \frac{\partial c}{\partial \bar{r}} = \frac{\partial^{2} c}{\partial {\bar{r}}^{2}},

(120)

where

s = 1 - cos θ,

with

0 \leq s \leq 2,

subject to

c = 1 at \bar{r} = 0,

(121)

c \to 0 as \bar{r} \to \infty,

(122)

c = 0 at s = 0,

(123)

c = 0 at T = 0.

(124)

Moreover, we can write

\bar{ψ} = {\bar{r}}^{2} u_{\infty} (s),

where

u_{\infty} (s) = \frac{3}{4} s (2 - s),

so that Eq. becomes

\frac{1}{4 {\bar{ψ}}^{1 / 2} u_{\infty}} \frac{\partial c}{\partial T} + \frac{\partial c}{\partial ϕ} = \frac{\partial}{\partial \bar{ψ}} ({\bar{ψ}}^{1 / 2} \frac{\partial c}{\partial \bar{ψ}})

(125)

subject to

c = 1 at \bar{ψ} = 0,

(126)

c \to 0 as \bar{ψ} \to \infty,

(127)

c = 0 at s = 0,

(128)

c = 0 at T = 0,

(129)

where

ϕ = 2 \int_{0}^{s} \sqrt{u_{\infty} (s^{'})} d s^{'};

for later use, note that

ϕ = \frac{\sqrt{3}}{2} (θ - \frac{1}{2} sin 2 θ),

(130)

with the interval

0 \leq θ \leq π

transforming to

0 \leq ϕ \leq π \sqrt{3} / 2.

Although the form of is more compact than that of , Eq. has singularities at

Ψ = 0

and

s = 0, 2;

it would therefore be more convenient to consider the solution of , subject to –. Introducing the transformations

η = \frac{\bar{r}}{T^{1 / 2}}, c = G (s, η, T) .

(131)

Equation gives

T \frac{\partial G}{\partial T} - \frac{1}{2} η \frac{\partial G}{\partial η} + \frac{3}{2} η T^{3 / 2} (s (2 - s) \frac{\partial G}{\partial s} - η (1 - s) \frac{\partial G}{\partial η}) = \frac{\partial^{2} G}{\partial η^{2}},

(132)

subject to

G = 1 at η = 0,

(133)

G \to 0 as η \to \infty,

(134)

G = 0 at s = 0,

(135)

G = 0 at T = 0.

(136)

As

T \to 0,

we obtain

- \frac{1}{2} η \frac{d G}{d η} = \frac{d^{2} G}{d η^{2}}

(137)

subject to and , giving

G = erfc (\frac{1}{2} η) .

(138)

Notice now that will also satisfy both and . For general values of T, one would now need to solve , subject to the boundary conditions and the initial condition . Although we will not do so here, this can in principle be done numerically using the Box scheme for unsteady boundary-layer equations.⁶

In the case of steady state, we set

\partial c / \partial T

in Eq. to zero and drop . Identifying

η = \bar{ψ} / ϕ^{2 / 3}

as the similarity variable, and then setting

η = Y^{2},

we arrive at

\frac{d^{2} c}{d Y^{2}} + \frac{4 Y^{2}}{3} \frac{d c}{d Y} = 0.

(139)

Hence, on applying the boundary conditions, we obtain

c = \frac{\int_{Y}^{\infty} exp (- \frac{4}{9} {\bar{Y}}^{3}) d \bar{Y}}{\int_{0}^{\infty} exp (- \frac{4}{9} {\bar{Y}}^{3}) d \bar{Y}},

(140)

whence

{(\frac{\partial c}{\partial \bar{r}})}_{r = 1} = - \frac{{(\sqrt{3} / 2)}^{2 / 3} sin θ}{{(θ - \frac{1}{2} sin 2 θ)}^{1 / 3} \int_{0}^{\infty} exp (- \frac{4}{9} {\bar{Y}}^{3}) d \bar{Y}} .

(141)

An important reference quantity is the average Sherwood number, Sh, which is given for a sphere by

S h = \frac{\int_{0}^{2 π} \int_{0}^{π} {(- \partial c / \partial r)}_{r = 1} sin θ d θ d ϕ}{\int_{0}^{2 π} \int_{0}^{π} sin θ d θ d ϕ};

(142)

in the case of axisymmetric flow, this reduces to

S h = \frac{1}{2} \int_{0}^{π} {(- \frac{\partial c}{\partial r})}_{r = 1} sin θ d θ .

(143)

Thence, after some manipulation of , we obtain

S h = {(\frac{81 π^{2}}{2048})}^{1 / 3} {(\int_{0}^{\infty} exp (- \frac{4}{9} {\bar{Y}}^{3}) d \bar{Y})}^{- 1} P e^{1 / 3};

(144)

now, with

\int_{0}^{\infty} exp (- \frac{4}{9} {\bar{Y}}^{3}) d \bar{Y} = {(\frac{1}{12})}^{1 / 3} Γ (\frac{1}{3}),

(145)

where Γ is the gamma function given by

Γ (z) = \int_{0}^{\infty} x^{z - 1} e^{- x} d x, Re (z) > 0,

(146)

we arrive at

S h = C^{*} P e^{1 / 3},

(147)

where

C^{*} = \frac{{(\frac{243 π^{2}}{512})}^{1 / 3}}{Γ (\frac{1}{3})} \approx 0.6245,

(148)

as was first derived by Acrivos and Taylor,³ with higher-order corrections given by Acrivos and Goddard,² albeit in the context of heat transfer.

As is evident, Eq. differs from the often-cited form of the Ranz–Marshall correlation,²⁸

S h_{R M} = 2 + 0.6 R e^{1 / 2} S c^{1 / 3},

(149)

in the following ways:

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} \to 2$ as $R e \to 0$ if $S c ≪ R e^{- 3 / 2},$ but (147) gives that this should be $S c ≪ R e^{- 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 ≪ 1 :$ in this case, the solution for c is given by

$c = \frac{1}{r},$

(150)

leading to $S h = 1.$ In fact,

$S h_{R M} (P e_{diam}) = 2 S h (P e),$

(151)

where $P e_{diam} = 2 P e .$ Moreover, (149) gives

$S h = 1 + 0.3 \sqrt{2} R e^{1 / 2} S c^{1 / 3} .$

(152)

IV. NUMERICAL IMPLEMENTATION

There remain two numerical tasks:

the solution of the full transient axisymmetric Eqs. (23)–(25), subject to (28)–(36);
the solution of the diffusion-dominated spherically symmetric problem given by Eq. (91) subject to (92)–(96).

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 smoothing⁴⁸ in these situations; however, it is well-established that the latter is more effective,²⁰ 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.

For all cases, the same convergence criterion at each time step was used, namely,

{(\frac{1}{N_{dof}} \sum_{i = 1}^{N_{dof}} {(\frac{| E_{i} |}{A + R | U_{i} |})}^{2})}^{\frac{1}{2}} < 1,

(153)

where E_i is the solver's estimate of the absolute error in the latest approximation to the ith component of the scaled solution vector, U_i, at each time step, A is the absolute tolerance and

R

is the relative tolerance; for simplicity, we take

R = A

and denote this common value by

ε .

Moreover,

N_{dof}

is the number of degrees of freedom (DOF), corresponding to a mesh having

N_{el}

elements, with

N_{s}

elements at the dissolution surface; Table I shows how

N_{el}, N_{dof}

⁠, and

N_{s}

are related for four different meshes that were experimented with.

TABLE I.

Number of elements (⁠ $N_{el}$ ⁠), degrees of freedom (⁠ $N_{dof}$ ⁠), and surface elements (⁠ $N_{s}$ ⁠) for four different meshes.

$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 $- P e^{- 1 / 3} \partial c / \partial r$ at r = 1 as a function of $π - θ$ 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.

FIG. 2.

Comparison of the analytical solution obtained from Eq. (141) and numerical solutions for − P e − 1 / 3 ∂ c / ∂ r at r = 1, with P e = 10 2 , 10 4 , 10 6 and R e = 10 − 3 .

View large Download slide

Comparison of the analytical solution obtained from Eq. (141) and numerical solutions for $- P e^{- 1 / 3} \partial c / \partial r$ at r = 1, with $P e = 10^{2}, 10^{4}, 10^{6}$ and $R e = 10^{- 3} .$

Next, in Fig. 3, which is basically a grid independence test, we consider the fraction of the remaining solid as a function of dimensionless time, t, for

P e = 10^{3}

for four different meshes; here, the remaining volume is given by

V (t) = 2 π \int_{0}^{R (θ, t)} r^{2} sin θ d θ,

whence the fraction remaining,

f (t),

is given by

f (t) = \frac{3 V (t)}{4 π} .

FIG. 3.

$Remaining fraction, f, as computed for four different meshes as a function of dimensionless time, t, for R e = 10 − 3 , P e = 10 3 , c 0 / c sat = 2 , ρ = 1 , c ∞ = 0.$

View large Download slide

Remaining fraction, f, as computed for four different meshes as a function of dimensionless time, t, for $R e = 10^{- 3}, P e = 10^{3}, c_{0} / c_{sat} = 2, ρ = 1, c_{\infty} = 0.$

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. 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 $P e ≪ 1$ 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 $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,²⁹ which clearly gives a consistent overestimation of the remaining fraction.

FIG. 4.

$Remaining fraction, f, as a function of dimensionless time, t, for R e = 10 − 3 , P e = 10 − 2 , c 0 / c sat = 2 , ρ = 1 , c ∞ = 0.$

View large Download slide

Remaining fraction, f, as a function of dimensionless time, t, for $R e = 10^{- 3}, P e = 10^{- 2}, c_{0} / c_{sat} = 2, ρ = 1, c_{\infty} = 0.$

V. RESULTS

At this stage, the full problem depends on four dimensionless parameters:

R e, S c, ρ, C .

In what follows, we consider results for which

ρ = 1, c_{\infty} = 0,

meaning that we focus on the influence of Re, Sc and

C = c_{0} / c_{sat} - 1 (> 0) .

In Sec. V A, we focus on the solutions obtained for different values of Sc and

C

at low values of Re; while in Sec. V B, we compare the results from this full-model approach against those obtained from simplified treatments of the problem.

A. Effect of Re, Sc, and $C$

Figure 5 shows the extent of remaining particle for $R e = 10^{- 3}$ and $P e = 10^{- 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.²⁴ 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 $θ = 0.$ This result is in line with Fig. 2.

FIG. 5.

View large Download slide

Particle shape after 0% (solid line), 30% (dashed line), 60% (dotted line), and 90% (dotted-dashed line) of the original volume has dissolved. $R e = 10^{- 3}, c_{0} / c_{sat} = 290$ ⁠: (a) $P e = 10^{- 3}$ and (b) $P e = 10^{4}$ ⁠.

Although it is now possible to generate many figures like Fig. 5 for different combinations of dimensionless parameters, we focus instead summarizing such results in terms of the instantaneous sphericity of the dissolving particle. As originally defined by Wadell,⁴⁴ the sphericity,

Ψ

⁠, of a particle is the ratio of the surface area of a sphere with the same volume as the given particle to the surface area of the particle, and is given by

Ψ = \frac{π^{1 / 3} {(6 V_{p})}^{2 / 3}}{A_{p}},

(154)

where V_p is volume of the particle and A_p is the surface area of the particle. The sphericity of a sphere is unity by definition and any particle which is not a sphere will have sphericity less than 1. Thence, Fig. 6 shows

Ψ

as a function of normalized dissolution time for different values of Pe, with

R e = 10^{- 3}

and

c_{0} / c_{sat} = 290

⁠; here, the normalization in time is carried out with the time taken for 90% of the particle to dissolve, which we denote by

t_{d, 90}

_%. From this figure, it is evident that sphericity decreases as Pe increases, although even when Pe is as high as 10⁴ the value of

Ψ

is still above 0.98, in spite of the fact that the corresponding profile in Fig. 5(b) appeared far from spherical. Figure 7 shows a corresponding plot for different values of

c_{0} / c_{sat}

⁠, with

P e = 10^{4}

and

R e = 10^{- 3} .

Thus, even though

c_{0} / c_{sat}

varies by several orders of magnitude, it is clear that the profiles are very similar, and are affected rather by the high value of Pe. We have not given a corresponding plot for a low value of Pe, as

Ψ = 1

to a good approximation in this case; similarly, the results for intermediate values of Pe can be expected to lie between the line

Ψ = 1

and the curves shown in Fig. 7.

FIG. 6.

View large Download slide

Sphericity of the granule geometry up to 90% of dissolution for different values of Pe, with $R e = 10^{- 3}$ and $c_{0} / c_{sat} = 290$ ⁠. Note that the curves for $P e = 10^{- 3}$ and 10⁻¹ are practically on top of each other.

FIG. 7.

View large Download slide

Sphericity as a function of normalized dissolution time for different values of $c_{0} / c_{sat}$ ⁠, with $P e = 10^{4}, R e = 10^{- 3}$ ⁠. Note that the curves for $c_{0} / c_{sat} = 290$ and 2900 are practically on top of each other.

B. Comparison with simplified models

As seen from the above results, the dissolving particle does not necessarily maintain its initially spherical shape. However, an often-used simplified approach to model dissolution is to assume that the particle behaves as a shrinking sphere during dissolution and that it shrinks at a rate that is proportional to the Sherwood number that would prevail under steady-state conditions;^15,16,24 this allows the use of already available correlations for the Sherwood number for flow past a sphere as a function of the Reynolds and Schmidt numbers.^11,28 Mathematically, this approach leads to a first-order ordinary differential equation for the particle radius,

r_{d},

as a function of time, as outlined in Appendix C. Moreover, whereas the Ranz–Marshall correlation is often used,²⁸ this may not necessarily be the most suitable. For example, Clift et al.¹¹ gave the correlation

S h = \frac{1}{2} (1 + {(1 + 2 P e)}^{1 / 3})

(155)

as being appropriate when

R e \leq 1;

note that when

P e ≫ 1,

it is evident that

S h \sim 2^{- 2 / 3} P e^{1 / 3},

which is in reasonable quantitative agreement with , noting that

2^{- 2 / 3} \approx 0.63

⁠. A comparison of the correlations due to Acrivos and Taylor,³ Ranz and Marshall²⁸ and Clift et al.¹¹—Eqs. , , and , respectively—is given in Fig. 8 for

R e = 10^{- 3};

these are also compared against the full numerical solution to the steady-state version of Eqs. , subject to , , , and , with ρ = 1 and

R \equiv 1

⁠. From this figure, it is clear that the Ranz–Marshall correlation severely underpredicts the Sherwood number for large values of Pe, as compared to the other three methods, which give more or less the same result. Thus, it is now of interest to see how well these approaches fare when applied to the actual dissolution problem, and how well their results agree with those given by our full numerical method; we turn to this next.

FIG. 8.

View large Download slide

Sherwood number (Sh) as a function of Péclet number (Pe) for steady flow past a sphere (⁠ $R e = 10^{- 3}$ ⁠).

First of all, note that the different approaches indicate different numbers of governing dimensionless parameters. In our approach, these are $R e, S c, ρ$ ⁠, and $C,$ with the list reducing to $P e (= ReSc), ρ$ ⁠, and $C,$ when $R e ≪ 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^{- 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, ρ$ ⁠, 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.¹¹ 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³ 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 ≪ 1,$ as the correlation was derived based on the assumption that $P e ≫ 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.

FIG. 9.

View large Download slide

Time taken to dissolve 90% of the particle, $t_{d, 90}$ _%, when $R e = 10^{- 3} .$

TABLE II.

Time, in seconds, taken to dissolve 90% of the initial drug load, according to different correlations and full time-dependent numerical simulation. $R e = 10^{- 3}, ρ = 1, c_{\infty} = 0, c_{0} / c_{sat} = 290.$

	Pe
Method	10⁻³	10⁻¹	10²	10⁴
Ranz and Marshall²⁸	4.40 $\times 10^{- 3}$	4.27 $\times 10^{- 1}$	2.93 $\times 10^{2}$	1.30 $\times 10^{4}$
Clift et al.¹¹	3.90 $\times 10^{- 3}$	3.78 $\times 10^{- 1}$	1.23 $\times 10^{2}$	3.03 $\times 10^{3}$
Acrivos and Taylor³	⋯	⋯	1.47 $\times 10^{2}$	3.20 $\times 10^{3}$
Full numerical	3.70 $\times 10^{- 3}$	3.63 $\times 10^{- 1}$	1.23 $\times 10^{2}$	3.05 $\times 10^{3}$

	Pe
Method	10⁻³	10⁻¹	10²	10⁴
Ranz and Marshall²⁸	4.40 $\times 10^{- 3}$	4.27 $\times 10^{- 1}$	2.93 $\times 10^{2}$	1.30 $\times 10^{4}$
Clift et al.¹¹	3.90 $\times 10^{- 3}$	3.78 $\times 10^{- 1}$	1.23 $\times 10^{2}$	3.03 $\times 10^{3}$
Acrivos and Taylor³	⋯	⋯	1.47 $\times 10^{2}$	3.20 $\times 10^{3}$
Full numerical	3.70 $\times 10^{- 3}$	3.63 $\times 10^{- 1}$	1.23 $\times 10^{2}$	3.05 $\times 10^{3}$

VI. CONCLUSION

In this paper, we have considered the canonical problem of the dissolution of a solid spherical particle under forced convection conditions. Although some aspects of this problem have been treated on previous occasions,^{8,12,17,18,36,47} we have attempted to provide an overarching generalized description that encompasses all earlier work. Starting with conservation laws for mass, momentum and dissolved species, we have formulated an axisymmetric moving-boundary problem for the evolving location of the surface of the dissolving particle. Nondimensionalization leads to a formulation that depends on four dimensionless parameters:

R e (= \frac{ρ_{sat} U a}{μ}), S c (= \frac{μ}{ρ_{sat} D}), ρ (= \frac{ρ_{s}}{ρ_{sat}}), C (= \frac{c_{0} - ρ c_{sat}}{c_{sat} - c_{\infty}}) .

(156)

However, when

R e ≪ 1,

the problem depends only on three:

P e (= ReSc), ρ

⁠, and

C

⁠. The subsequent development was for the case when

R e ≪ 1

⁠; this is of particular relevance to drug particle dissolution. The governing equations were then solved numerically using finite-element methods. The subsequent results enabled us to determine the accuracy and appropriateness for dissolution modeling of the often-used Ranz–Marshall correlation;²⁸ this was found to severely overestimate the dissolution in convection-dominated cases, i.e., when

P e ≫ 1

⁠. On the other hand, a less-often used correlation due to Clift et al.¹¹ appeared to be more accurate in general, even well beyond the range in Pe anticipated by the original authors. Nevertheless, the use of both correlations requires the assumption that the particle remains spherical as it dissolves; however, from the full numerical computations, this was found to be far from the case for the convection-dominated regime. A further correlation, based on a boundary-layer analysis due to Acrivos and Taylor,³ was found to work well for the convection-dominated regime; however, as the particle dissolves, it is inevitable that a diffusion-dominated regime is encountered, rendering the correlation as it stands inappropriate for determining the time to complete dissolution.

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 $R e ≪ 1.$ 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 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 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.

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 ≪ ρ_{sat} B g a^{2} / μ$

If (41) is not satisfied, e.g., in situations where the solvent is assumed to be stagnant, a different nondimensionalization is necessary.

Setting

\begin{matrix} r^{'} = \frac{r}{a}, R^{'} = \frac{R}{a}, t^{'} = \frac{t}{a^{2} / D}, c^{'} = \frac{c - c_{\infty}}{c_{sat} - c_{\infty}}, \\ u^{'} = \frac{u}{D / a}, u_{0}^{'} = \frac{u_{0}}{D / a}, p^{'} = \frac{p + ρ_{sat} g (1 - B) z}{μ D / a^{2}}, \end{matrix}

(A1)

and with

B ≪ 1

⁠, we would now have

\nabla \cdot u = 0,

(A2)

\frac{1}{S c} (\frac{\partial u}{\partial t} + (u \cdot \nabla) u) = - \nabla p + \nabla^{2} u - Rac k,

(A3)

\frac{\partial c}{\partial t} + \nabla \cdot (c u) = \nabla^{2} c,

(A4)

subject to

u \cdot t = 0,

(A5)

u \cdot n = (ρ - 1) \frac{\partial R}{\partial t},

(A6)

\nabla c \cdot n = C \frac{\partial R}{\partial t},

(A7)

c = 1,

(A8)

u \to P e (cos θ, - sin θ, 0), c \to 0,

(A9)

where Ra denotes the Rayleigh number and is given by

R a = \frac{ρ_{sat} β (c_{sat} - c_{\infty}) g a^{3}}{μ D} .

(A10)

We can now note that

U ≪ ρ_{sat} B g a^{2} / μ

implies that

P e ≪ R a .

Alternatively, this can be rewritten in terms of the Richardson number, Ri, which is more commonly used as a measure of whether forced or natural convection dominates:³⁵ with

R i = Bga / U^{2},

natural convection is dominant if

RiRe ≫ 1,

where Re is as defined in .

Thence, on setting

t = T / R a, u = R a U, p = RaP,

(A11)

we obtain

\nabla \cdot U = 0,

(A12)

\frac{R a}{S c} (\frac{\partial U}{\partial T} + (U \cdot \nabla) U) = - \nabla P + \nabla^{2} U - c k,

(A13)

\frac{\partial c}{\partial T} + \nabla \cdot (c U) = \frac{1}{R a} \nabla^{2} c,

(A14)

subject to the following: at the surface of the particle,

U \cdot t = 0,

(A15)

U \cdot n = (ρ - 1) \frac{\partial R}{\partial T},

(A16)

\nabla c \cdot n = R a C \frac{\partial R}{\partial T},

(A17)

c = 1;

(A18)

as

r \to \infty,

U \to 0, c \to 0.

(A19)

Observe that this problem is also axisymmetric, albeit about an axis that is parallel to the gravity vector in Fig. 1.

APPENDIX B: SOLUTION TO EQUATION (98) WHEN $S > 2 π$

For

S > 2 π,

we set

λ^{2} = S^{2} / π - 2 S > 0,

so that Eq. gives

\int^{T} \frac{d T}{T} = - 2 \int^{V} \frac{V d V}{{(V + \frac{S}{π^{1 / 2}})}^{2} - λ^{2}} .

(B1)

Integrating both sides gives

\frac{1}{λ^{2}} {(R^{(0)} + \frac{S T^{1 / 2}}{π^{1 / 2}})}^{2} - T = C {(\frac{\frac{R^{(0)}}{T^{1 / 2}} + \frac{S}{π^{1 / 2}} + λ}{\frac{R^{(0)}}{T^{1 / 2}} + \frac{S}{π^{1 / 2}} - λ})}^{S / λ π^{1 / 2}},

(B2)

where C is a constant to be determined. Since

R^{(0)} (0) = 1,

we must have

C = 1 / λ^{2};

so,

{(R^{(0)} + \frac{S T^{1 / 2}}{π^{1 / 2}})}^{2} - T = {(\frac{\frac{R^{(0)}}{T^{1 / 2}} + \frac{S}{π^{1 / 2}} + λ}{\frac{R^{(0)}}{T^{1 / 2}} + \frac{S}{π^{1 / 2}} - λ})}^{S / λ π^{1 / 2}} .

(B3)

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,²⁸ although others are available also.¹¹ Mathematically, this approach leads to a first-order ordinary differential equation for the particle radius, $r_{d},$ as a function of time, as follows.

Mass balance considerations give

\frac{d}{d t} (\frac{4}{3} ρ_{s} π r_{d}^{3}) = - 4 π r_{d}^{2} k (c_{sat} - c_{\infty}),

(C1)

where

\frac{k r_{d}}{D} = S h;

(C2)

combining, we have

ρ_{s} r_{d} \frac{d r_{d}}{d t} = - DSh (c_{sat} - c_{\infty}) .

(C3)

At this stage, it is convenient to nondimensionalize through

r_{d} = a R_{d}, t = (\frac{a^{2} ρ_{s}}{D (c_{sat} - c_{\infty})}) τ,

(C4)

giving

R_{d} \frac{d R_{d}}{d τ} = - S h,

(C5)

subject to

R_{d} (0) = 1.

Ranz–Marshall:²⁸ In this case, with

S h = 1 + 0.3 \sqrt{2} R e_{d}^{1 / 2} S c^{1 / 3} .

(C6)

Equation becomes

R_{d} \frac{d R_{d}}{d τ} = - (1 + λ R_{d}^{1 / 2}),

(C7)

where

λ = 0.3 \sqrt{2} R e^{1 / 2} S c^{1 / 3};

note here that

R e_{d} = ρ_{sat} U r_{d} / μ,

which we rewrite as

R e_{d} = R e R_{d},

on using . Setting

Y = 1 + λ R_{d}^{1 / 2},

we have

\frac{2}{λ^{4}} \int_{1 + λ R_{d}^{1 / 2}}^{1 + λ} \frac{{(Y - 1)}^{3} d Y}{Y} = \int_{0}^{τ} d τ,

(C8)

whence

τ = \frac{2}{λ^{4}} {\frac{1}{6} ((1 + λ) (2 λ^{2} - 5 λ + 11) - (1 + λ R^{1 / 2}) (2 λ^{2} R - 5 λ R^{1 / 2} + 11)) + ln (\frac{1 + λ R^{1 / 2}}{1 + λ})} .

(C9)

This gives the nondimensional dissolution time,

τ_{d},

on setting

R_{d} = 0,

as

τ_{d} = \frac{1}{3 λ^{4}} {(1 + λ) (2 λ^{2} - 5 λ + 11) - 11 - 6 ln (1 + λ)};

(C10)

the dimensional dissolution time,

t_{d, 100}

_%, is then given by

t_{d, 100 %} = (\frac{a^{2} ρ_{s}}{D (c_{sat} - c_{\infty})}) τ_{d} .

(C11)

Clift et al.,¹¹ Eq. (3-49), for

R e \leq 1

⁠: Now, instead of , we have

S h = \frac{1}{2} (1 + {(1 + 2 P e)}^{1 / 3}),

(C12)

giving

R_{d} \frac{d R_{d}}{d τ} = - \frac{1}{2} (1 + {(1 + μ R_{d})}^{1 / 3}),

(C13)

where

μ = 2 R eSc;

thence,

2 \int_{R_{d}}^{1} \frac{R_{d} d R_{d}}{1 + {(1 + μ R_{d})}^{1 / 3}} = \int_{0}^{τ} d τ .

(C14)

Putting

Z = 1 + {(1 + μ R_{d})}^{1 / 3},

we have

τ = \frac{6}{μ^{2}} \int_{1 + {(1 + μ R_{d})}^{1 / 3}}^{1 + {(1 + μ)}^{1 / 3}} \frac{[{(Z - 1)}^{5} - {(Z - 1)}^{2}] d Z}{Z},

(C15)

whence

τ = \frac{6}{μ^{2}} {F (μ) - F (μ R_{d}) - 2 ln (\frac{1 + {(1 + μ)}^{1 / 3}}{1 + {(1 + μ R_{d})}^{1 / 3}})},

(C16)

where

F (μ) = \frac{1}{5} {(1 + μ)}^{\frac{2}{3}} (μ - 4) + \frac{1}{4} {(1 + μ)}^{\frac{1}{3}} (7 - μ) + \frac{μ}{3} + \frac{247}{60} .

(C17)

From , the dissolution time is given by

τ_{d} = \frac{6}{μ^{2}} {F (μ) - \frac{76}{15} - 2 ln (\frac{1 + {(1 + μ)}^{1 / 3}}{2})} .

(C18)

Acrivos and Taylor,³ for

R e ≪ 1, P e ≫ 1

: Using , we have

R_{d} \frac{d R_{d}}{d τ} = - ν R_{d}^{1 / 3},

(C19)

where

ν = C^{*} {(ReSc)}^{1 / 3},

whence it is straightforward to obtain

R_{d} = {(1 - \frac{5}{3} ν τ)}^{3 / 5},

(C20)

leading to

τ_{d} = \frac{3}{5 ν} .

(C21)

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 ≫ 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.

FIG. 10.

View large Download slide

Time taken to dissolve 100% of the particle, $t_{d, 100}$ _%, when $R e = 10^{- 3} .$

REFERENCES

1.

D. J.

Acheson

,

Elementary Fluid Dynamics

(

Oxford University Press

,

Oxford

,

1990

).

Google Scholar

2.

A.

Acrivos

and

J. D.

Goddard

, “

Asymptotic expansions for laminar forced convection heat and mass transfer. Part 1. Low speed flows

,”

J. Fluid Mech.

23

,

273

–

291

(

1965

).

https://doi.org/10.1017/S0022112065001350

Google Scholar

Crossref

3.

A.

Acrivos

and

T. D.

Taylor

, “

Heat and mass transfer from single spheres in Stokes flow

,”

Phys. Fluids

5

,

387

–

394

(

1962

).

https://doi.org/10.1063/1.1706630

Google Scholar

Crossref

4.

R.

Asthana

and

S. K.

Pabi

, “

An approximate solution for the finite-extent moving-boundary diffusion-controlled dissolution of spheres

,”

Mater. Sci. Eng., A

128

,

253

–

258

(

1990

).

https://doi.org/10.1016/0921-5093(90)90233-S

Google Scholar

Crossref

5.

R.

Bird

,

W.

Stewart

, and

E.

Lightfoot

,

Transport Phenomena

,

1st ed.

(

John Wiley

,

New York

,

1960

).

Google Scholar

6.

P.

Bradshaw

,

T.

Cebeci

, and

J. H.

Whitelaw

,

Engineering Calculation Methods for Turbulent Flow

(

Academic Press

,

London

,

1981

).

Google Scholar

7.

E.

Brunner

, “

Reaktionsgeschwindigkeit in heterogenen Systemen

,”

Z. Phys. Chem.

47U

,

56

–

102

(

1904

).

https://doi.org/10.1515/zpch-1904-4705

Google Scholar

Crossref

8.

M.

Cable

and

J. R.

Frade

, “

The diffusion-controlled dissolution of spheres

,”

J. Mater. Sci.

22

,

1894

–

1900

(

1987

).

https://doi.org/10.1007/BF01132424

Google Scholar

Crossref

9.

H.

Cao

,

X.

Jia

,

Y.

Li

,

C.

Amador

, and

Y.

Ding

, “

CFD-DNS simulation of irregular-shaped particle dissolution

,”

Particuology

50

,

144

–

155

(

2020

).

https://doi.org/10.1016/j.partic.2019.08.003

Google Scholar

Crossref

10.

Y. W.

Chen

and

P. J.

Wang

, “

Dissolution of spherical solid particles in a stagnant fluid—An analytical solution

,”

Can. J. Chem. Eng.

67

,

870

–

872

(

1989

).

https://doi.org/10.1002/cjce.5450670522

Google Scholar

Crossref

11.

R.

Clift

,

J.

Grace

, and

M. E.

Weber

,

Bubbles, Drops, and Particles

(

Academic Press

,

1978

).

Google Scholar

12.

M.

Doichinova

,

O.

Lavrenteva

, and

C.

Boyadjiev

, “

Non-linear mass transfer from a solid spherical particle dissolving in a viscous fluid

,”

Int. J. Heat Mass Transfer

54

,

2998

–

3003

(

2011

).

https://doi.org/10.1016/j.ijheatmasstransfer.2011.02.051

Google Scholar

Crossref

13.

A.

Dokoumetzidis

and

P.

Macheras

, “

A century of dissolution research: From Noyes and Whitney to the biopharmaceutics classification system

,”

Int. J. Pharm.

321

,

1

–

11

(

2006

).

https://doi.org/10.1016/j.ijpharm.2006.07.011

Google Scholar

Crossref

PubMed

14.

J. L.

Duda

and

J. S.

Vrentas

, “

Heat or mass transfer-controlled dissolution of an isolated sphere

,”

Int. J. Heat Mass Transfer

14

,

395

–

407

(

1971

).

https://doi.org/10.1016/0017-9310(71)90159-1

Google Scholar

Crossref

15.

D. M.

D'Arcy

and

T.

Persoons

, “

Mechanistic modelling and mechanistic monitoring: Simulation and shadowgraph imaging of particulate dissolution in the flow-through apparatus

,”

J. Pharm. Sci.

100

,

1102

–

1115

(

2011

).

https://doi.org/10.1002/jps.22337

Google Scholar

Crossref

PubMed

16.

D. M.

D'Arcy

and

T.

Persoons

, “

Understanding the potential for dissolution simulation to explore the effects of medium viscosity on particulate dissolution

,”

AAPS Pharm. Sci. Tech

20

,

47

(

2019

).

https://doi.org/10.1208/s12249-018-1260-4

Google Scholar

Crossref

17.

T.

Elperin

and

A.

Fominykh

, “

Effect of solute concentration level on the rate of coupled mass and heat transfer during solid sphere dissolution in a uniform fluid flow

,”

Chem. Eng. Sci.

56

,

3065

–

3074

(

2001

).

https://doi.org/10.1016/S0009-2509(01)00011-2

Google Scholar

Crossref

18.

J. Q.

Feng

, “

Diffusion-controlled quasi-stationary mass transfer for an isolated spherical particle in an unbounded medium

,”

Chem. Eng. Commun.

200

,

65

–

76

(

2013

).

https://doi.org/10.1080/00986445.2012.690355

Google Scholar

Crossref

19.

J. J. L.

Higdon

and

M.

Kojima

, “

On the calculation of Stokes-flow past porous particles

,”

Int. J. Multiphase Flow

7

,

719

–

727

(

1981

).

https://doi.org/10.1016/0301-9322(81)90041-0

Google Scholar

Crossref

20.

P.

Knupp

, “

Winslow smoothing on two-dimensional unstructured meshes

,”

Eng. Comput.

15

,

263

–

268

(

1999

).

https://doi.org/10.1007/s003660050021

Google Scholar

Crossref

21.

A.

Kovács

,

C. J. W.

Breward

,

K. E.

Einarsrud

,

S. A.

Halvorsen

,

E.

Nordgård-Hansen

,

E.

Manger

,

A.

Münch

, and

J. M.

Oliver

, “

A heat and mass transfer problem for the dissolution of an alumina particle in a cryolite bath

,”

Int. J. Heat Mass Transfer

162

,

120232

(

2020

).

https://doi.org/10.1016/j.ijheatmasstransfer.2020.120232

Google Scholar

Crossref

22.

K. M.

Moroney

,

L.

Kotamarthy

,

I.

Muthancheri

,

R.

Ramachandran

, and

M.

Vynnycky

, “

A moving boundary model of dissolution from binary drug-excipient granules incorporating granule microstructure

,”

Int. J. Pharm.

599

,

120219

(

2021

).

https://doi.org/10.1016/j.ijpharm.2021.120219

Google Scholar

Crossref

PubMed

23.

K. M.

Moroney

and

M.

Vynnycky

, “

Mathematical modelling of drug release from a porous granule

,”

Appl. Math. Model.

100

,

432

–

452

(

2021

).

https://doi.org/10.1016/j.apm.2021.07.023

Google Scholar

Crossref

24.

M.

Navas-Bachiller

,

T.

Persoons

, and

D. M.

D'Arcy

, “

Exploring bulk volume, particle size and particle motion definitions to increase the predictive ability of in vitro dissolution simulations

,”

Eur. J. Pharm. Sci.

174

,

106185

(

2022

).

https://doi.org/10.1016/j.ejps.2022.106185

Google Scholar

Crossref

PubMed

25.

W.

Nernst

, “

Theorie der Reaktionsgeschwindigkeit in heterogenen Systemen

,”

Z. Phys. Chem.

47U

,

52

–

55

(

1904

).

https://doi.org/10.1515/zpch-1904-4704

Google Scholar

Crossref

26.

A. A.

Noyes

and

W. R.

Whitney

, “

The rate of solution of solid substances in their own solutions

,”

J. Am. Chem. Soc.

19

,

930

–

934

(

1897

).

https://doi.org/10.1021/ja02086a003

Google Scholar

Crossref

27.

B. S.

Padmavathi

,

T.

Amaranath

, and

S. D.

Nigam

, “

Stokes-flow past a porous sphere using Brinkman model

,”

Z. Angew. Math. Phys.

44

,

929

–

939

(

1993

).

https://doi.org/10.1007/BF00942818

Google Scholar

Crossref

28.

W. E.

Ranz

and

W. R.

Marshall

, “

Evaporation from drops. 1

,”

Chem. Eng. Prog.

48

,

141

–

146

(

1952

).

Google Scholar

29.

R. G.

Rice

and

D. D.

Do

, “

Dissolution of a solid sphere in an unbounded, stagnant liquid

,”

Chem. Eng. Sci.

61

,

775

–

778

(

2006

).

https://doi.org/10.1016/j.ces.2005.08.003

Google Scholar

Crossref

30.

L. E.

Scriven

, “

On the dynamics of phase growth

,”

Chem. Eng. Sci.

10

,

1

–

13

(

1959

).

https://doi.org/10.1016/0009-2509(59)80019-1

Google Scholar

Crossref

31.

J.

Siepmann

and

F.

Siepmann

, “

Mathematical modeling of drug dissolution

,”

Int. J. Pharm.

453

,

12

–

24

(

2013

).

https://doi.org/10.1016/j.ijpharm.2013.04.044

Google Scholar

Crossref

PubMed

32.

C.

So

,

P.-C.

Chiang

, and

C.

Mao

, “

Modeling drug dissolution in 3-dimensional space

,”

Pharm. Res.

39

,

907

–

917

(

2022

).

https://doi.org/10.1007/s11095-022-03270-6

Google Scholar

Crossref

PubMed

33.

A. C.

Srivastava

and

N.

Srivastava

, “

Flow past a porous sphere at small Reynolds number

,”

Z. Angew. Math. Phys.

56

,

821

–

835

(

2005

).

https://doi.org/10.1007/s00033-005-2006-1

Google Scholar

Crossref

34.

K.

Sugano

, “

Theoretical comparison of hydrodynamic diffusion layer models used for dissolution simulation in drug discovery and development

,”

Int. J. Pharm.

363

,

73

–

77

(

2008

).

https://doi.org/10.1016/j.ijpharm.2008.07.002

Google Scholar

Crossref

PubMed

35.

D. J.

Tritton

,

Physical Fluid Dynamics

(

Van Nostrand Rheinhold

,

UK

,

1977

).

Google Scholar

Crossref

36.

S.

Tseng

,

T.-H.

Lin

, and

J.-P.

Hsu

, “

Unsteady dissolution of particle of various shapes in a stagnant liquid

,”

Chem. Eng. Sci.

123

,

573

–

578

(

2015

).

https://doi.org/10.1016/j.ces.2014.11.023

Google Scholar

Crossref

37.

J. S.

Vrentas

,

C. M.

Vrentas

, and

H. C.

Ling

, “

Equations for predicting growth or dissolution rates of spherical-particles

,”

Chem. Eng. Sci.

38

,

1927

–

1934

(

1983

).

https://doi.org/10.1016/0009-2509(83)85049-0

Google Scholar

Crossref

38.

M.

Vynnycky

and

S.

Kimura

, “

An analytical and numerical study of coupled transient natural convection and solidification in a rectangular enclosure

,”

Int. J. Heat Mass Transfer

50

,

5204

–

5214

(

2007

).

https://doi.org/10.1016/j.ijheatmasstransfer.2007.06.036

Google Scholar

Crossref

39.

M.

Vynnycky

and

S.

Kimura

, “

Can natural convection alone explain the Mpemba effect?

,”

Int. J. Heat Mass Transfer

80

,

243

–

255

(

2015

).

https://doi.org/10.1016/j.ijheatmasstransfer.2014.09.015

Google Scholar

Crossref

40.

M.

Vynnycky

and

M. A.

O'Brien

, “

An asymptotic and numerical study of slow, steady ascent in a Newtonian fluid with temperature-dependent viscosity

,”

Appl. Math. Comput.

219

,

3154

–

3177

(

2012

).

https://doi.org/10.1016/j.amc.2012.09.049

Google Scholar

Crossref

41.

M.

Vynnycky

and

M. A.

O'Brien

, “

Slow, steady ascent in a power-law fluid with temperature-dependent viscosity

,”

J. Non-Newtonian Fluid Mech.

195

,

9

–

18

(

2013

).

https://doi.org/10.1016/j.jnnfm.2012.12.001

Google Scholar

Crossref

42.

M.

Vynnycky

and

M. A.

O'Brien

, “

The slow, steady ascent of a hot solid sphere in a Newtonian fluid with strongly temperature-dependent viscosity

,”

Appl. Math. Comput.

231

,

231

–

253

(

2014

).

https://doi.org/10.1016/j.amc.2013.12.078

Google Scholar

Crossref

43.

M.

Vynnycky

and

M.

Zambrano

, “

Towards a ‘moving-point’ formulation for the modelling of oscillation-mark formation in the continuous casting of steel

,”

Appl. Math. Model.

63

,

243

–

265

(

2018

).

https://doi.org/10.1016/j.apm.2018.06.029

Google Scholar

Crossref

44.

H.

Wadell

, “

Volume, shape, and roundness of quartz particles

,”

J. Geol.

43

,

250

–

280

(

1935

).

https://doi.org/10.1086/624298

Google Scholar

Crossref

45.

J. Z.

Wang

and

D. R.

Flanagan

, “

General solution for diffusion controlled dissolution of spherical particles. 1. Theory

,”

J. Pharm. Sci.

88

,

731

–

738

(

1999

).

https://doi.org/10.1021/js980236p

Google Scholar

Crossref

PubMed

46.

J. Z.

Wang

and

D. R.

Flanagan

, “

General solution for diffusion-controlled dissolution of spherical particles. 2. Evaluation of experimental data

,”

J. Pharm. Sci.

91

,

534

–

542

(

2002

).

https://doi.org/10.1002/jps.10039

Google Scholar

Crossref

PubMed

47.

Y.

Wang

,

B.

Abrahamsson

,

L.

Lindfors

, and

J. G.

Brasseur

, “

Comparison and analysis of theoretical models for diffusion-controlled dissolution

,”

Mol. Pharm.

9

,

1052

–

1066

(

2012

).

https://doi.org/10.1021/mp2002818

Google Scholar

Crossref

PubMed

48.

A.

Winslow

, “

Numerical solution of the quasilinear Poisson equations in a nonuniform triangle mesh

,”

J. Comput. Phys.

1

,

149

–

172

(

1966

).

https://doi.org/10.1016/0021-9991(66)90001-5

Google Scholar

Crossref

On the dissolution of a solid spherical particle

I. INTRODUCTION

II. MODEL EQUATIONS

III. ANALYSIS

A. Nondimensionalization

B. $P e ≪ 1$

1. $| ρ - 1 | / P e ≲ 1$

2. $| ρ - 1 | / P e ≫ 1$

C. $P e ≫ 1$

IV. NUMERICAL IMPLEMENTATION

A. Full equations

B. Diffusion-dominated problem

V. RESULTS

A. Effect of Re, Sc, and $C$

B. Comparison with simplified models

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: $U ≪ ρ_{sat} B g a^{2} / μ$

APPENDIX B: SOLUTION TO EQUATION (98) WHEN $S > 2 π$

APPENDIX C: SIMPLIFIED MODELS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

On the dissolution of a solid spherical particle

I. INTRODUCTION

II. MODEL EQUATIONS

III. ANALYSIS

A. Nondimensionalization

B. P e ≪ 1

1. | ρ − 1 | / P e ≲ 1

2. | ρ − 1 | / P e ≫ 1

C. P e ≫ 1

IV. NUMERICAL IMPLEMENTATION

A. Full equations

B. Diffusion-dominated problem

V. RESULTS

A. Effect of Re, Sc, and C

B. Comparison with simplified models

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: U ≪ ρ sat B g a 2 / μ

APPENDIX B: SOLUTION TO EQUATION (98) WHEN S > 2 π

APPENDIX C: SIMPLIFIED MODELS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

B. $P e ≪ 1$

1. $| ρ - 1 | / P e ≲ 1$

2. $| ρ - 1 | / P e ≫ 1$

C. $P e ≫ 1$

A. Effect of Re, Sc, and $C$

APPENDIX A: $U ≪ ρ_{sat} B g a^{2} / μ$

APPENDIX B: SOLUTION TO EQUATION (98) WHEN $S > 2 π$