Universal expression for the drag on a fluid sphere

An expression was developed for prediction of drag coefficients for any spherical particle, drop or bubble in an infinite, homogeneous liquid. The formula reproduces the limiting cases for gas bubbles and solid spheres, as well as the exact Hadamard-Rybczynski solution. The accuracy of the expression, which is valid for Reynolds numbers up to a few hundred, was confirmed by comparison with published numerical predictions of the drag coefficient for a range of physical circumstances.


Introduction
Bubbles, drops and particles are widespread in science and engineering phenomena. Knowledge of the behavior of single bubbles and drops is not only directly relevant to many applications; it also supports understanding of the corresponding swarms, e.g., [1].
Wegener et al. [2] recently provided a comprehensive review of theory, experimental data and pertinent approximations describing the dynamics of single drops in fluid systems. The steady rate of movement of spherical particle, drops and bubbles is quantified by the drag coefficient, C d . Of relevance here is Wegener et al.'s summary of drag formulas (their Tables 1 and  2), and the ranges of (drop) Reynolds numbers over which the different formulas apply.
We consider the drag on a single spherical solid particle, liquid drop or gas bubble (collectively referred to as a particle) moving in an otherwise quiescent, infinite, homogeneous liquid, without interphase mass transfer (i.e., Sherwood number of zero). Spherical particles will occur if surface tension is dominant, or if inertia is negligible. Clift et al. [3] and Bhaga and Weber [4] provide systematic characterizations of particle shape as it varies with Reynolds, Eötvös and Morton numbers. These authors show, for instance, that spherical shapes for Reynolds numbers of about 100 or higher are found for sufficiently small Eötvös and Morton numbers.
Below, we develop and test a new expression for C d applicable to spherical particles. Our approach is to build an interpolation based on known limiting cases (e.g., small Reynolds number), as well as a validated large Reynolds number drag expression.

Drag formula interpolations between gas bubbles and solid spheres
There are several approximations for C d in the literature for the limiting cases of bubbles or solid spheres [5][6][7][8]. As mentioned, Wegener et al. [2] provide a summary of these approximations, remarking that they "approximate only certain intervals of the standard drag curve." A well-known interpolation is that of Rivkind and Ryskin [9]: with where ρ is the density, μ the viscosity, U the (steady) velocity of the sphere, R 0 the Reynolds number, X the viscosity ratio (X ! 0 for a gas bubble and 1 for a solid sphere) and a the sphere radius. Subscript "0" refers to the material outside the sphere and "1" to that within. Although Eq (1) is straightforward, it does not approach known limits, as given below.

Analytical drag results for large Reynolds number
Unlike the situation for R 0 ( 1, exact results for R 0 ) 1 do not exist. For laminar flow (i.e., oscillations of the particle do not occur), Harper and Moore [13] as well as Parlange [14] obtained approximate expressions for the drag for this case, however. In both [13] and [14], it was observed that, to a first approximation, flow inside the particle is described by a Hill's vortex, and outside by a potential flow. Both approaches give drag predictions that are "numerically indistinguishable" [14]. This conclusion follows from the minor effect on the drag made by slightly different assumptions in [13] and [14]. Barry and Parlange [15] compared predictions of both theories to experimental results on recirculation within the particle [4], and found that the theory of Parlange [14] is more accurate. Thus, this theory is the starting point for the developments presented below. The drag formula of Parlange [14] is: where P = ρ 1 r À 1 0 is the density ratio. The two constants, α and β, are defined by integrals. Parlange [14] simplified Z by taking α = 2 and β = 1. Numerical evaluations (Appendix) give α and β as: In the following, we use α and β as given by Eq (7) For later convenience in manipulating Eq (6), we define A as: New drag formula Eq (6) holds for R 0 ) 1, so it is not surprising that Eq (3) (the Hadamard-Rybczynski limit) is not obtained for R 0 ( 1. However, we can force Eq (6) to do so by using standard Padé approximations [16]. First, we rewrite Eq (6) as: which is identical in order to Eq (6) but additionally approaches Eq (3) as R 0 ! 0. Next, we modify Eq (9) so that it will reduce to Eq (5) (with χ = 2), in which case it must be corrected for R 0 ( 1 without affecting predictions for R 0 ) 1. This occurs with a Padé approximant that maintains the first two orders of Eqs (4) and (6) for R 0 small and large, respectively. We proceed stepwise and satisfy the different limiting cases given above. First, to ensure that the predictions of Eq (9) are not affected for R 0 ) 1, we take corrections that are exponentially small in that limit, and so replace Eq (9) with: which reduces to Eq (9) for R 0 ) 1 (or for τ = 0) and Eq (3) for R 0 ! 0. Second, the parameters τ, λ and ω (all > 0) are chosen to ensure that Eq (10) reduces to the other limits given above. For this purpose, we observe that Eq (5) requires a Padé approximant in powers of R 0 , rather than R 1=2 0 , appearing in Eq (10) The initial appearance of R 1=2 0 is removed from the small-R 0 expansion of Eq (10) if: In addition, to satisfy Eq (5) (with χ = 2) for R 0 ( 1 requires that τ is given by: Once τ is obtained using Eq (12), λ and ω follow straightforwardly from Eq (11). The only arbitrary element in the deviation of Eq (10) is the form of the exponential corrections. Alternative functional forms were investigated, with the most promising being bð1 þ lR 1=2 0 Þ À 1 . However, based on comparisons with published numerical results (below), we found that only the exponential form reduced to Eq (9) rapidly enough for R 0 ) 1. Eq (10) constitutes a new, fully analytical, expression for the drag of a spherical particle. It reduces quickly to Eq (9) for R 0 ) 1. All coefficients are determined by the behavior of C d under different conditions, i.e., no empirical coefficient is determined by curve fitting of numerical results.
Because the bubble (X = 0) and the solid sphere (X ! 1) are oft-investigated special cases, we present Eq (10) for these limits.

Comparison with numerical results
Eq (10) is compared with numerical results from the literature. For convenience, numerous comparisons are collected in S1 Text. Specifically, tables of results from Eq (10) are compared with those from [3,11,[17][18][19][20][21][22][23][24][25][26][27]. Representative results are presented below. First, we consider in Fig 1 the gas bubble (X = 0), where results from the literature are collected. The figure compares Eq (10) with values given by Clift et al. [3], Oliver and Chung [11], Brabston and Keller [18] and Magnaudet et al. [25]. For the data of Clift et al. [3], the predictions agree well with the published values except for R 0 ! 200. This is not surprising as Clift et al. [3] obtained values at R 0 = 300 and 400 by interpolation with higher Reynolds numbers, when Eq (10) does not apply. Eqs (1) and (5) are accurate over their reported ranges of validity. The numerical simulations of Brabston and Keller [18] for R 0 up to 200 are in close agreement with Eq (10). The accuracy of the numerical results of Brabston and Keller is confirmed from their agreement with the numerical results of Oliver and Chung [11] over the narrower range of 1/2 < R 0 50. The results of Magnaudet et al. [25] agree well with Eq (10) for the entire range of R 0 considered. Fig 2 makes the same comparison for a solid sphere. Interestingly, the disagreement with the Clift et al. [3] values, and the predictions of the Rivkind and Ryskin [9] formula, Eq (1), is very small and is limited to R 0 between approximately 10 and 100. Again, the drag values reported by Oliver and Chung [11], Chang et al. [23] and Chang and Maxey [24] are all similar, and agree more closely with Eq (10) than the results reported by Clift et al., although the results from Eq (10) seem slightly high. This figure includes results calculated from the formula of Flemmer and Banks ( [28], Eq (7), which tends to be slightly lower than the numerical results. Besides the Flemmer and Banks [28] formula, there are several expressions for the solid sphere drag coefficient available. However, as shown by Mikhailov and Silva Freire [29], the largest variations between them occur at about R 0 = 100, with a maximum deviation of about 5%, so the Flemmer and Banks formula can be taken as representative.
Comparisons for various X are given in Table 1, which lists drag values computed by Oliver and Chung [11] for X = 0, 0.333, 1, 3, 1 and 1/2 R 0 50. The results for X = 0 and X ! 1 were already presented in the figures, of course. As expected, the agreement for intermediate X values is excellent and is similar to that shown in Figs 1 and 2. Table 1 includes the estimates calculated with Eq (1), which are less accurate than those from Eq (10). Table 2 lists numerical results from Juncu [30], who considered different density ratios, P, for R 0 = 100. We mention, in passing, that C d is largely insensitive to changes in P, and if P is not specified, typically P = X is assumed (as done here). The agreement between the numerical results and Eq (10) is excellent. Again, the estimates of Eq (10) are usually slightly above the numerical values. Table 2 includes predictions from Eq (1), although this expression does not account for variations of C d with P.

Conclusion
We developed a formula, Eq (10), to predict the drag for a spherical particle, for all viscosity ratios between gas bubbles and solid spheres. It was derived as a Padé approximant that interpolates between known analytical results at low and moderate Reynolds numbers assuming that the particle does not oscillate or wobble. The formula can be used to predict the drag coefficient for spherical bubbles, drops and particles for any viscosity and density ratios. Surface tension is assumed sufficient to maintain the spherical shape of the particle. Both Fig 2. 5 of Clift et al. [3] and Fig [8] of Bhaga and Weber [4] show how the particle shape changes (spherical, ellipsoidal, spherical cap, etc.) with the Eötvös, Morton and Reynolds numbers, and thus indicate the range of conditions for which Eq (10) applies. Eq (10) is more theoretically based than, for instance, the formula of Rivkind and Ryskin [9] and appears to be more accurate as well, especially at low Reynolds numbers.