Drag Coefficient Correction for Spherical and Nonspherical Particles

Article pubs.acs.org/IECR

Drag Coefficient Correction for Spherical and Nonspherical Particles Suspended in Square Microducts Zachary D. Hensley and Dimitrios V. Papavassiliou* School of Chemical, Biological, and Materials Engineering, University of Oklahoma, Norman, Oklahoma 73019-1004, United States ABSTRACT: Even though drag in creeping flow has been extensively studied in the past, recent developments in microfluidic technologies have led to a reevaluation of the drag equations that were theoretically derived for ideal cases. Stokes’ Law is the bulwark of equations describing drag in creeping flow, but it applies in cases of uniform flow around a spherical object. Empirical and theoretical corrections to the Stokes drag for other particle shapes have been proposed, as well as corrections that account for the presence of walls that affect the shape of the velocity profile around the particle. The present work utilizes computational fluid dynamics to investigate the validity and the differences between results obtained with these equations at low Reynolds number flows. In addition, simpler equations using the ratio of the diameter of the particle to the height of the channel and the aspect ratio of nonspherical particles are used to predict drag for particles that are spherical, ellipsoidal, and cylindrical in square microducts. simulations.11 Simulations with DPD have been used to provide flow behavior around individual nanorods and to calculate quantities such as drag and torque in order to compare the calculated results to results obtained from equations available in the literature.12 The theoretically expected results were used mostly to validate the DPD results. The contribution of the present work is centered on the analysis of the effects of the size of the microchannel relative to the size of a particle on the drag coefficient and on the lift forces. Effects of the particle location relative to the microchannel wall are also investigated. In addition to spherical particles, prolate spheroids (ellipsoidal particles) and cylindrical particles of different aspect ratios are examined, and corrections to the Stokes equation are provided.

1. INTRODUCTION Research in microfluidics has become increasingly important to the advancement of nanomaterials and micromachining technologies. Industrial processing, biological applications, and colloids science can benefit from a deeper understanding of microflows. It has also been suggested that microfluidics holds promise for large-scale industrial processes to be miniaturized and run in parallel.1 In bioinspired applications, microchannels have been used to study low Reynolds numbers blood flow. Experiments have established how particles order and separate by lift forces caused by the flow lamina.2 Other bio applications include drug delivery with carbon nanotubes (CNTs) and treatment of cancer with hyperthermia.3,4 Recently, gold-infused nanoparticles have been used to detect cancer cells and may serve as a backbone for molecularly targeted photodiagnostics and therapy.5 For biology, the concept of a “lab on a chip” is a goal that is gradually approached with the improving manipulation and understanding of microscale properties.6 Furthermore, multiphase flow experiments in microducts have been designed to determine the regions of slug, slug/annular, bubbly, and pulse flow.7 Some researchers have tackled microfluidic flows using simulations based on atomistic, Lattice Boltzmann, and continuum approaches.8 Simulations for particle settling and interactions for colloidal mixtures, however, can prove to be computationally intensive due to the solid−fluid interactions and varying shapes of particles. Models were often based on simplifications, such as the assumption that all particles were spheres in low Reynolds number flow, and generalized equations were developed to describe the particle behavior.9 Some current methods involve fluid particle dynamics (FPD) to simplify the problem and provide an explanation for hydrodynamic interactions and particle aggregation.10 Another simulation technique that may be appropriate for examining flows over nanoscale particles is dissipative particle dynamics (DPD). DPD is a coarse-graining computational technique that allows a set of governing equations to be solved for time scales that are much larger than the time scales used in molecular dynamics © 2014 American Chemical Society

2. BACKGROUND Stokes derived a fundamental equation for drag in creeping flow around a spherical particle, where the drag coefficient was given as C D,Stokes =

24 Rep

(1)

The Stokes solution applies for spherical particles in uniform flow in an infinite domain at low Rep. Correction factors, however, must be applied on the basis of velocity profile exposure and on particle shape. Fidleris and Whitmore designed an elegant experiment involving several fluids and particles of different densities to study the effect of the presence of walls. Both the laminar and turbulent regime were studied, and it was conclusively shown that a correction to account for the distance from the wall is needed.13 These experiments confirmed that the ratio of the diameter of the particle to the diameter of the circular Received: Revised: Accepted: Published: 10465

February 21, 2014 May 19, 2014 May 22, 2014 May 22, 2014 dx.doi.org/10.1021/ie5007646 | Ind. Eng. Chem. Res. 2014, 53, 10465−10474

Industrial & Engineering Chemistry Research

Article

⎡ CD 24 ⎢ = 1 + 0.1118(RepK1K 2)0.6567 K2 RepK1K 2 ⎢⎢ ⎣

pipe is important for drag, as proposed by the Landenburg equation:14 u∞ = umax

1 dp

(1 + 2.4 ) D

(2)

+

The experiments of Fidleris and Whitmore12 also helped to produce equations for the turbulent flow regime. The accuracy of Ladenburg’s correction factor for a sphere in a pipe was first questioned by Faxen and later by Wakiya.15 Faxen proved that Ladenburg’s equation was the first approximation of a series and it could be expanded to give a more accurate value for the drag coefficient. In the early 1950s, Wakiya, using the method of reflections, stated that the coefficient appearing in the denominator of eq 2 should be 2.1 instead of 2.4 for a sphere placed in the center of a pipe.16 Until the 1960s, particles of shapes other than spherical had not been considered in-depth. Happel and Brenner provided fundamental relationships for ellipsoidal and cylindrical particles in an infinite domain at low Reynolds numbers using asymptotic theory. (In the present work, the word ellipsoidal is used to indicate a prolate spheroid particle that occurs with the rotation of an ellipse around its long axis, following Happel and Brenner.) The concept of an equivalent radius of a particle was used to provide a correction for such nonspherical particle shapes, given as follows:17 r=

8c 3 ⎡ 2φ ⎢ − φ2 − 1 + ⎣

2φ − 1 2

3/2

(φ − 1)

⎛φ+ ln⎜ ⎝φ−

1

⎛ G ⎞2 um = 0.343ru∞⎜ ⎟ ⎝v⎠

⎞⎤ ⎟⎥ − 1 ⎠⎦

φ

(3)

1

FL = 6πμu∞r(0.343Rev 2 )

The radius calculated by this equation was then used as the radius of a sphere and was substituted into the well-known Stokes equation for the force acting on a sphere

FD = 6πμu∞r

CL =

(4)

24 C (1 + ARep B) + E Rep 1 + Re

p

(8)

2FL 2 Ap ρu∞

(9)

The premise of this work is to analyze how the correction to the drag coefficient because of the presence of walls varies in the Stokes regime for a microduct and to compare this correction with previous corrections for a pipe. While the velocity profile for fully developed laminar flow in a pipe is given as

For finite Rep, Levenspiel, in the late 1980s, introduced a fourparameter empirical equation of the form18 CD =

(7)

Later works analyzed the cases of multiple spheres and the effect that the walls have on lift. The importance of such experiments was the finding that particles always moved away from the wall. In the case of flow between two walls, a particle was found to move to a position of equilibrium.22 McLaughlin proposed a more complex equation to better predict lift for a larger range of Reynolds numbers, resulting in a formula that includes the tensor component of the lift.23 Miyazaki continued to expand on the work of Saffman, analyzing the lift force as a tensor for simple shear flow.24 This leads to the overall force equation for lift on a sphere and the coefficient of lift, both of which are defined below:

φ2 − 1 2

(6)

Although drag is the main focus of this article, we discuss here the importance of the lift force. While particles were known to experience a sort of lift in low Reynolds flows, it was not until Saffman produced correlations in the 1960s that this phenomenon could be explained. Saffman’s relation depends strictly on the assumption that Rep < Rev < 1. This eventually led to the following equation for velocity migration:21

1 2

⎤ 0.4305 ⎥ 3305 ⎥ 1 + Re K K ⎥ p 1 2 ⎦

⎡ ⎛ r ⎞2 ⎤ u(r ) = umax ⎢1 − ⎜ ⎟ ⎥ ⎝R⎠ ⎦ ⎣

(5)

This model can be applied to both spherical and nonspherical particles. Levenspiel and Haider defined each of the coefficients in eq 5 as a function of sphericity. (Sphericity is a parameter that can be used to describe how close a particle shape is to that of a sphere and is defined as the ratio of the surface area of a sphere with the same volume as the particle divided by the actual surface area of the particle.) An issue with this equation is that it fails to take into account wall effects. A shape factor scruple was introduced recently in an attempt to describe the drag of irregularly shaped particles. This scruple is based on sphericity. Combined with Levenspiel’s four-parameter model, the number of types of particles for which the equation could be applied was further expanded.19 More recent works describe how to quantify drag for a large range of Reynolds numbers using two parameters, K1 and K2, that are functions of sphericity, dp/D, and φ. Ganser developed the following empirical equation to describe drag for all particles with Rep K1K2 < 105:20

(10)

the velocity profile for flow in a duct is given as an infinite series,25 u(y , z) =

1 ΔP 2 (z − b 2 ) + 2 μX cosh

∞

∑ A n cos (2n + 1)πz n=0

(2n + 1)πy 2b

2b

(11)

ΔP

An =

16b2 μX ( −1)n [(2n + 1)π ]3 cosh

(2n + 1)πa 2b

(12)

We considered the case of a square microduct for this work; however, it should be noted that varying one side of the channel can have a profound impact on the velocity profile. Lastly, the coefficient of lift was also calculated for each simulation run and the results are presented. 10466

dx.doi.org/10.1021/ie5007646 | Ind. Eng. Chem. Res. 2014, 53, 10465−10474


Article

the microducts was 40 μm, there were no wake effects for all cases simulated. Velocity measurements were taken along the center streamline. This ensured that the velocity behind the particle reached u∞ (no wake) again before the periodic condition was applied. The simulations would terminate when the residuals, as calculated by the software, were less than 10−5. In total, more than 150 separate simulation cases were run to obtain the data presented herein.

3. METHODOLOGY AND SIMULATION The simulation runs were performed using ANSYS Fluent 14.0, a finite volume solver. Geometries and meshing were created in ICEM, with the exception of elliptical geometries that were generated in Design Modeler.26 Each mesh was composed of strictly hexahedral elements. Curved surfaces were o-gridded to reduce skewness. To validate the numerical approach and to explore the numerical mesh size that is adequate for the simulations, we considered the case of a microduct without suspended particles. These initial simulations were conducted to confirm that the velocity profile in the microduct agreed with eq 11. The number of computational cells was changed between 40 000 and 200 000 and the root-mean-square (rms) of the error analyzed. The rms error for all of these simulations was on the order of 10−3; however, subsequent runs were carried out with roughly 200 000 cells in order to obtain more descriptive velocity contour maps. The solver is a solution of the integral form of the Navier−Stokes equation with finite volume methods. This is equivalent to a direct numerical simulation for laminar flow. The forces on a particle were directly calculated by the software as surface integrals on the surface of the particles. Equation 9 was used for the computation of the lift coefficient, and the drag coefficient was computed as CD = 2 FD/(ρu∞2Ap). The case of a spherical particle in an infinite domain was simulated next. The spherical particle case was run with the computational box having zero shear at the outer boundaries, and the drag force on the particle was compared to the results of Stokes Law. Ellipsoidal and cylindrical particle cases were run with a computational box having no shear, and the results were compared to Stokes Law using Happel and Brenner’s equivalent radius.17 To verify that the numerical grid does not affect the simulation answer, a comparison of drag force between runs was used for spherical particles and Rep = 0.1. If the difference in drag was less than 3%, then the current simulation answer was used. When the drag force converged, the ratio of the drag coefficient obtained by the simulation divided by the drag force calculated by Stokes’ law was 1.03. The observed 3% deviation from Stokes law is within expectations for the case of Rep = 0.1, which is at the upper limit of the applicability of the Stokes equation. An example of a particle in flow that could have a Rep at about 0.1 is a red blood cell in a saline solution with fluid velocity on the order of 0.01 cm/s. A Rep that is 2 orders of magnitude smaller could be obtained for particles like multiwalled carbon nanotubes in aqueous flows with velocity on the order of 0.01 cm/s. The height of the channel was varied to obtain different dp/H aspect ratios. The computational mesh was nonuniform; it was finer on the surface of the particles, where the grid spacing was 0.1 μm, and increased with a ratio of 1.2 from the particle surface. The mass flow rate was varied to ensure a constant bulk velocity (uavg = 0.1 m/s) through the microducts between simulation runs. This kept the Rep = 0.1 when particles were located in the center of the microduct. For the case of generating data to construct Moody type plots, the mass flow rate changed to obtain different bulk velocities, as needed. Periodic boundary conditions were specified in the streamwise direction. For the periodic boundary conditions to make physical sense, the node arrangement at the entrance and exit of the computational box had to be the same. The channel length was initially varied to ensure that the channel was sufficiently long to prevent the velocity of the fluid around the particle to be affected by the wake of the particle that is located upstream in the periodic computational box. It was determined that when the length for

4. RESULTS 4.1. Spherical Particles. As Ladenburg14 proved roughly 100 years ago, the ratio of the diameter of the particle divided by

Figure 1. Drag coefficient for spheres of three different diameters centered in a microchannel and with varying the size of the microduct. (Rep = 0.1)

Figure 2. Drag coefficient for a spherical particle with a diameter of 0.5 μm placed off-center at three different locations in the microduct. (Rep = 0.1 for the particle placed in the center of the channel.)

the diameter of the pipe leads to a simple drag correction. It was therefore assumed that drag on a sphere in a micro channel should follow a similar trend. For the case of a square duct, the hydraulic diameter is equal to the height of the duct. Runs were performed for spheres of various sizes dp, while the height of the duct was modified to obtain sets of data with the same dp/H but different dp (Figure 1). Note that as dp/H → 0, the value of the drag coefficient approaches the value of the drag coefficient for Stokes flow because the particle diameter becomes very small relative to the height of the microduct. In that case, the flow is very similar to flow around a particle in an infinite domain. It was important to also verify that the placement of the sphere at different heights in the channel did not play a major role in the drag force. The sphere was placed at four different distances from 10467



Article

Figure 3. All data points for spherical particles plotted against the parameter dp/H in order to estimate the deviation from Stokes’ Law. (Rep = 0.1)

Figure 5. Drag coefficients for an ellipsoidal particle with φ = 8 that was placed at different locations within the microduct to determine the difference in drag. (Rep = 0.1 for the particle located at the center of the duct.)

Figure 4. Moody type plot for different values of dp/H. The mass flow rate was adjusted to obtain different Rep values. All particles were centered in the channel. The lines are drawn on the basis of the results obtained using eq 13.

Figure 6. Several ellipsoidal particles with varying φ were placed in the center of the microduct to determine the importance of the particle aspect ratio, φ. (Rep = 0.1)

the duct wall: center of the duct (ΔL = 0), shifted closer to the wall by one micrometer (ΔL = 1 μm), shifted toward the wall by two micrometers (ΔL = 2 μm), and shifted toward the wall by three micrometers (ΔL = 3 μm). In each case, the flow rate was changed through the duct so that the average fluid velocity entering the duct was 0.1 m/s. This ensured that if a particle was placed at the center of the channel its Reynolds number would be Rep = 0.1. Drag increased slightly for the same dp/H as the particle approached the wall. In Figure 2, we present the slight variations in drag that were obtained from the simulations for various particle locations. A linear regression was performed for all spherical particle data, and the results are presented in Figure 3.

The value of the drag coefficient using a linear correlation is slightly higher than that established by Ladenburg14 and later corrected by Wakiya.16 The reason for this is the difference in the shape of the velocity profile in a square duct versus a pipethe velocity profile in a square duct is dependent on both the y and z positions, while fully developed velocity in a pipe is a function of only the radial position. Furthermore, the velocity fields in pipes and square ducts have different symmetry, in particular, full rotational symmetry in circular pipes does not apply in square ducts. It is, thus, clearly demonstrated that a significant enough difference occurs when comparing drag for particles suspended in ducts to the drag for particles in pipes. This leads to the

Table 1. Error Analysis for Drag Coefficient for Currently Established Equations in the Literature and Spherical Particlesa

dp/H 0.0625 0.1 0.125 a

Fidleris equation ⎡ ⎤ ⎢ ⎥ 24 ⎢ 0.407 ⎥ 0.681 CD = (1 + 0.15Rep ) + ⎛ ⎞ ⎢ Rep 8710 ⎥ ⎜1 + ⎟ ⎢⎣ Re p ⎠ ⎥ ⎝ ⎦ 14.3% 18.8% 7.64%

Ganser equation (eq 6)

present correlation (eq 13)

10.0% 12.3% 6.9%

1.1% 0.1% 2.1%

The percent differences are provided for cases of different parameter dp/H. 10468



Article

Figure 7. All data for ellipsoidal particles plotted against the proposed parameter (φ)2/5(dp/H). (Rep = 0.1)

Figure 8. A Moody type plot generated with changes in the parameter (φ)2/5(dp/H). The lines are drawn on the basis of the results obtained using eq 14.

suggestion of a simple linear correction for spheres in a square duct, as follows: CD C D,Stokes

⎛ dp ⎞ = ⎜1 + 2.5 ⎟ H⎠ ⎝

(13)

Additional runs were conducted at different Reynolds numbers with dp/H held constant to obtain data for a Moody type diagram, represented by Figure 4. The case of Stokes’ Law is represented when the ratio of dp/H approaches infinity. Table 1 is a summary of the deviation of established drag coefficient relations and the data obtained through the simulations. The proposed eq 13 exhibits dramatically smaller deviations than other equations. 4.2. Ellipsoidal Particles. With sufficient theory provided by Happel and Brenner 17 for ellipsoidal particles (prolate spheroids) in an infinite medium, it seemed reasonable to establish the correction factor in a duct. Again, we had to determine if the particle placement away from the duct center contributed greatly to drag, so the particle was placed at different locations away from the center of the duct while maintaining a constant bulk velocity (0.1 m/s). With the exception of two points, a linear relation between drag and dp/H is exhibited (Figure 5). These two points correspond to the case where dp was large relative to the size of the channel and the particle almost touched the wall. While it appears that using an equivalent radius would alleviate the need to supply a correction for φ, this hypothesis was tested. Different values of φ (φ = 8, 16, 32) were tested for ellipsoidal particles centered in a channel. The aspect ratio in this case is defined as the end-to-end distance of the particle divided by its diameter. Figure 6 is an indication that the functionality of the aspect ratio is not precisely captured by the Happel and Brenner17 relationship. The correction factor was

Figure 9. Drag coefficients for a cylindrical particle with φ = 8 that was placed at different locations within the microduct. (Rep = 0.1 for the particle paced at the center of the duct.)

then assumed to be a function of both φ and dp/H. A simple linear fitting does not work for this case, so an exponential relationship was pursued, and the appropriate exponent for φ was found to be 2/5. Figure 7 is an evaluation of the CD ratio versus the new parameter (φ)2/5(dp/H). The correction for the drag coefficient for elliptical particles was determined as follows: CD C D,H&B

⎡ 2 ⎛ d p ⎞⎤ = ⎢1 + 2.33(φ)5 ⎜ ⎟⎥ ⎢⎣ ⎝ H ⎠⎥⎦

(14)

Table 2 is a presentation of the deviation of results obtained using different drag coefficient relations and the data obtained through the present simulations. It provides a quantitative measure of the agreement between eq 14 and the data. Similar to the case of

Table 2. Error Analysis for the Drag Coefficient Obtained for Ellipsoidal Particles When Using Different Equationsa

(φ)2/5(dp/H) 0.144 0.230 0.287 a



Happel and Brenner equations (eqs 3 and 4)

present correction (eq 14)

39.9% 28.3% 9.93%

22.6% 34.3% 40.0%

5.21% 3.09% 1.09%

The table provides the percent differences for different values of (φ)2/5(dp/H). 10469



Article

Figure 10. Results for runs with a different aspect ratio, φ, evaluated to determine the correction that must be applied for drag on cylindrical particles. (Rep = 0.1)

Figure 12. Moody type plot generated for cylindrical particles and different values of (φ)2/5(dp/H).

Figure 11. All data for cylinders plotted versus the same parameter that was generated from the data for ellipsoids. (Rep = 0.1)

Figure 13. Coefficient of lift for a sphere as a function of ΔL/H compared to the analytically obtained value.

spherical particles, a Moody type diagram was generated with different (φ)2/5(dp/H), as seen in Figure 8. 4.3. Cylindrical Particles. Due to the Stokes paradox, drag around cylindrical particles is a more difficult case to analyze theoretically. No exact analytical solution has been generated for drag around a finite cylinder because of the presence of the flat faces of the cylinders. For consistency purposes, Happel and Brenner’s17 expression for ellipsoidal particles in an infinite medium is again used to calculate the drag for cylinders. A ratio of φ = 8 was chosen, in analogy to the aspect ratio of 8 used for ellipsoidal particles, and the particle was placed at different locations in the duct with its long axis parallel to the direction of the flow. The particle aspect ratio in this case is defined as the ratio of the length of the cylinder divided by its diameter. The results are presented in Figure 9. Because it has already been

found (on the basis of our results for the ellipsoidal particles) that the φ aspect ratio is important, runs testing φ = 16 and 32 were also performed. Simulations with different placement relative to the duct wall were also completed for the multiple φ to demonstrate that regardless of the magnitude of the aspect ratio, the height is not an important factor when the bulk velocity is constant. Figure 10 is the presentation of the results, depicting this finding in graphical form. We note again that the effect of the aspect ratio φ on the drag coefficient is not linear, so we employed an exponential fit. The drag coefficient varies with the aspect ratio φ raised to the power of two-fifths for both cylindrical and ellipsoidal particles. Figure 11 is an evaluation of the drag on cylindrical particles versus φ2/5(dp/H). The correction for cylindrical particles is as follows:

Table 3. Overview of the Results for Different Cases of (φ)2/5(dp/H) with the Percent Deviation from the Simulation Findingsa

2/5

(φ) (dp/H) 0.115 0.23 0.287 a



Happel and Brenner equations (eqs 3 and 4)

Morgan and Green, 200327 F = 8πμLu∞/(2 ln(φ) − 1)

Xu et al., 200928 F = 4πμLu∞/(ln(φ))

present correction (eq 15)

66.4% 67.1% 65.1%

50.2% 14.0% 2.98%

2.43% 25.9% 33.1%

48.2% 12.5% 1.65%

0.85% 1.77% 1.23%

Several equations for the calculation of the drag coefficient have been derived in literature for the behavior of cylindrical particles. 10470



Article

Figure 14. Contours of the velocity magnitude in the xz-plane for flow around a sphere that is originally centered in the channel and then shifted down by 1 μm for each image. The legend is given in units of m/s.

CD CD,H&B

⎡ ⎛ d p ⎞⎤ = ⎢1 + 4.20(φ)2/5 ⎜ ⎟⎥ ⎢⎣ ⎝ H ⎠⎥⎦

lift force calculated through our simulations. For all spherical particle cases, a linear trend is observed; however, the Miyazaki equation is found to overestimate the lift observed in each run. In the case of lift, the parameter dp/H does not appear to be important in determining the lift. In fact, it is the displacement of the particle relative to the center of the channel, given as ΔL/H, that seemed to be of particular importance in describing lift behavior. The reason is that the difference in the velocity above and below the particle, which results in varying pressures on the top and bottom of the sphere, is what contributes to lift, and this difference would depend on the y-location of the particle. Contour plots are provided for spherical, ellipsoidal, and cylindrical particles (Figures 14−16, respectively) for each height placement within the duct. The distortion of the velocity profile due to the presence of the particle is evident from these figures. The velocity above the particle becomes slightly higher and the velocity below the particle (between the particle and the duct wall) becomes smaller, leading to lift phenomena. The

(15)

The deviation between the obtained data points from the simulations and different equations for the drag coefficient available in the literature for cylindrical particles are presented in Table 3. While the relation obtained by Xu et al.28 works well for a particular value of φ2/5(dp/H), eq 15 leads to consistently smaller deviations from the data. Lastly, Figure 12 is an evaluation of the effects of the new (φ)2/5(dp/H) parameter for various low Reynolds number flow cases. 4.4. Lift Forces on the Microparticles. The lift force on the microparticles was calculated for each simulation run. Prior work on lift forces has mainly focused on the case of spherical particles. Miyazaki24 and other researchers investigated the migration velocity, while Saffman’s initial equation provided a direct calculation of the lift force generated. Figure 13 is a comparison of the empirical correlation of Miyazaki (eq 8) compared to the 10471



Article

Figure 16. Contours of the velocity magnitude in the xz-plane for flow around a cylinder with aspect ratio 8 that is originally centered in the channel and then shifted down by 1 μm for each image. The legend determining the velocity magnitude is the same as in Figure 14.

Figure 15. Contours of the velocity magnitude in the xz-plane for flow around an ellipsoid with an aspect ratio of 8 that is originally centered in the channel and then shifted down by 1 μm for each image. The legend determining the velocity magnitude is the same as in Figure 14.

coefficient CL is found to vary linearly with ΔL/H, leading to the following simple relation:

Figure 17. Coefficient of lift for an ellipsoidal particle is plotted vs (ΔL/ H)(dp/L).

ΔL (16) H where the coefficient is found to be A = 5.52. Ellipsoidal particles behave different than spheres. There should be a strong dependence on ΔL/H for the lift on the particle; however, this is not sufficient enough to describe the lift fully. Because lift is in the direction perpendicular to the flow direction, the cross-sectional area of the particle in the flow direction relative to the cross-sectional area of the particle in vertical direction also plays a role. The ratio of the cross-sectional area perpendicular to the flow over the cross-sectional area parallel to the flow reduces to dp/L for the ellipsoidal particle case, while it is equal to 1 for a sphere. In Figure 17, we plot the lift coefficient as a function of (ΔL/H)(dp/L). The following equation was obtained on the basis of the plot: CL = A

Figure 18. Coefficient of lift for a cylindrical particle is plotted vs (ΔL/ H)(πdp/4L).

ΔL d p CL = B (17) H L 2 where the coefficient is B = 10.17, with an R value of 0.5.

The lift for cylindrical particles is expected to be a function of similar parameters as for ellipsoidal particles. The cross-sectional 10472



■

area parallel to the flow is different in this case, and the ratio of the cross-sectional areas reduces to πdp/4L. Simulations for cylinders with larger aspect ratios than that of the ellipsoidal particles were generated, and the lift coefficient was calculated. In Figure 18, we plot CL vs (ΔL/H)(πdp/4L), leading to the following equation: ΔL d pπ CL = C H 4L

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: (405) 325-5811. Fax: (405) 325-5813. Notes

The authors declare no competing financial interest.

■

(18)

where the coefficient is C = 15.07, with an R2 value of 0.52. The results in Figures 17 and 18, as well as the R2 values offered above, indicate that there could be room for developing a different set of dimensionless parameters in the equations for the calculation of the lift coefficient. Other parameters were indeed considered; however, this was the best option among many tried. The physical justification, as already mentioned, is to use the cross-sectional area of the particle that is perpendicular to the direction of the lift. To make this quantity dimensionless, the cross-sectional area in the flow direction was used. Other parameters resulted in much lower R2 values, and the data did not fall on a line, even for the spherical particles (as observed in Figure 13).

ACKNOWLEDGMENTS

Part of this work was done while D.V.P. was serving at the National Science Foundation (NSF). Any opinion, findings, conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF. This work was supported by the Advanced Energy Consortium (http://www.beg.utexas.edu/aec/). Member companies include BP America Inc., BG Group, Petrobras, Schlumberger, Statoil, Shell, Repsol, and Total.

■

5. CONCLUDING REMARKS As discussed earlier, different fluid velocity profiles result in differences to the drag experienced by a microparticle placed in a microfluidic device. The velocity field varies in pipes in the radial direction, while the velocity changes in a microduct with respect to two directions. Examining contours of velocity in ducts, one can see the significant change that occurs just by doubling the length of one of the sides of the duct cross section while keeping the bulk fluid velocity constant. It should seem reasonable that as one side continues to grow so that a ≫ b, slit flow is achieved. Further work is needed to provide an accurate representation of drag coefficient correction for the case of a ≠ b. Particle shape corrections for drag in low Reynolds number flows have been developed in the past for various particles; most of these tend to be tedious and do not account for the presence of walls. The present work has confirmed that drag for spherical particles is affected by the ratio dp/H, while for ellipsoidal and cylindrical particles, it is affected by both φ and dp/H. At first glance, one may ask why the correction coefficients for ellipsoidal particles and spherical particles are not the same. This is because Happel and Brenner’s17 equivalent radius is indeterminate for the case of φ = 1, even though as φ approaches zero the equivalent radius is equal to the sphere radius. In the denominator of eq 3, the logarithmic term reduces to 0 for φ = 1, while it is multiplied by a term that is undefined in this case. Because our results are for particles oriented with their long axes parallel to the wall, further work is needed to examine how particle orientation can affect the drag on ellipsoidal and cylindrical particles. The effect of the presence of multiple particles on lift and drag is another example of future work. The present results provide corrections to the Stokes equation that could be applied in the design of microducts where particles flow or to update simulation models that use an equation of motion to predict the trajectories of microparticles with nonspherical shapes (e.g., colloidal particles, microfibers, carbon nanotubes, etc.) in networks of microducts.

NOMENCLATURE a = Semiaxis of the microduct in the y direction Ap = Cross-sectional area perpendicular to the flow field A, B, C = Constants used in Levenspiel equation (eq 5) b = Semiaxis of the microduct in the z direction c = Semiaxis of the perpendicular to the flow field CD = Drag coefficient CD,Fluent = Drag coefficient from the simulations CD,H&B = Drag coefficient based on the analytical result of Happel and Brenner17 CD,Stokes = Drag coefficient from the Stokes equation CL = Coefficient of lift D = Pipe diameter dp = Diameter of the particle E = Constant used in Levenspiel18 equation (eq 5) FD = Drag force FL = Lift force G = Velocity gradient H = Height of the microduct [=] 2b K1, K2 = Constants used in Ganser20 equation (eq 6) L = End-to-end length of a particle (cylindrical or ellipsoidal) r = Radius of the particle R = Radius of the pipe Rep = Reynolds number based on the particle diameter [=] (ρu∞2dp/μ) Rev = Reynolds number based on shear rate [=] (βr2/v) u∞ = Streamline velocity at long distance upstream from the particle center um = Migration velocity of the particle umax = Maximum fluid velocity in a duct or pipe (located at the center line) X = Length of the duct y = Location in the y direction relative to the origin (located at the center of the duct) z = Location in the z direction relative to the origin (located at the center of the duct)

Greek Characters

β = Shear rate ΔL = Distance of particle center of gravity from the center of the duct ΔP = Pressure drop over the length of the duct μ = Dynamic viscosity of water v = Kinematic viscosity of water φ = Aspect ratio (length of particle/diameter of particle) 10473



■

Article

Customer. Retrieved from ANSYS Simulation Driven Product Development. (27) Morgan, H., Green, N.. AC Electrokinetics: Colloids and Nanoparticles; Research Studies Press: Philadelphia, PA, 2003. (28) Xu, Y.-Q.; Barnard, A.; McEuen, P. Bending and Twisting of Suspended Single-Walled Carbon Nanotubes in Solution. Nano Lett. 2009, 1−6.

REFERENCES

(1) Squires, T.; Quake, S. Microfluidics: Fluid physics at the nanoliter scale. Rev. Mod. Phys. 2005, 77, 977−1026. (2) Carlo, D.; Irimia, D.; Tompkins, R.; Toner, M. Continuous inertial focusing, ordering, and separation of particles in microchannels. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 18892−18897. (3) Kam, N.; OConnell, M.; Wisdom, J.; Dai, H. Carbon nanotubes as multifunctional biological transporters and near-infrared agents for selective cancer cell destruction. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 11600−11605. (4) Duong, H. M.; Papavassiliou, D. V.; Mullen, K. J.; Wardle, B. L.; Maruyama, S. A numerical study of effective thermal conductivity of biological fluids containing single-walled carbon nanotubes. Int. J. Heat Mass Transfer 2009, 52, 5591−5597. (5) Huang, X.; El-Sayed, I.; Qian, W.; El-Sayed, M. Cancer cell imaging and photothermal therapy in the near-infrared region by using gold nanorods. J. Am. Chem. Soc. 2005, 128, 2115−2120. (6) Beebe, D.; Mensing, G.; Walker, G. Physics and applications of microfluidics in biology. Annu. Rev. Biomed. Eng. 2002, 4, 261−286. (7) Tiggelaar, R.; Van der Schaaf, J.; De Loos, S. Gas-liquid dynamics at low Reynolds numbers in pillared rectangular micro channels. Microfluid. Nanofluid. 2010, 9, 131−144. (8) Worner, M. . Numerical modeling of multiphase flows in microfluidics and micro process engineering: A review of methods and applications. Microfluid. Nanofluid. 2012, 12, 841−886. (9) Tam, C. The drag on a cloud of spherical particles in low Reynolds number flow. J. Fluid Mech. 1969, 38, 537−546. (10) Hajime, T.; Takeaki, A. Simulation method of colloidal suspensions with hydrodynamic interactions: Fluid particle dynamics. Phys. Rev. Lett. 2010, 85, 1338−1341. (11) Groot, R. D.; Warren, P. B. Dissipative particle dynamics: Bridging the gap between atomistic and mesocopic simulation. J. Chem. Phys. 1997, 107, 4423−4435. (12) Chen, S.; Phan-Thien, N.; Khoo, B.; Fan, X. J. Flow around spheres by dissipative particle dynamics. Phys. Fluids 2006, 18, 103605− 103619. (13) Fidleris, V.; Whitmore, R. L. The physical interaction of spherical particles in suspensions. Rheologica Acta 1961, 1, 573−580. (14) Ladenburg, A. Uber Partielle Racemie. Ber. Dtsch. Chem. Ges. 1908, 966−969. (15) Faxen, H. The resistance against the movement of a rigour sphere in viscous fluids, which is embedded between two parallel layered barriers. Ann. Phys. (Berlin, Ger.) 1922, 68, 89−119. (16) Wakiya, S. A spherical obstacle in the flow of a viscous fluid through a tube. J. Phys. Soc. Jpn. 1953, 8, 254−257. (17) Happel, J., Brenner, H. Low Reynolds Number Hydrodynamics; Prentice-Hall: Englewood Cliffs, NJ, 1965. (18) Haider, A.; Levenspiel, O. Drag coefficient and terminal velocity of spherical and nonspherical particles. Powder Technol. 1989, 58, 63− 70. (19) Thompson, T. L.; Clark, N. N. A holistic approach to particle drag prediction. Powder Technol. 1991, 67, 57−66. (20) Ganser, G. H. A rational approach to drag prediction of spherical and nonspherical particles. Powder Technol. 1993, 77, 143−152. (21) Saffman, P. G. The lift on a small sphere in a slow shear flow. J. Fluid Mech. 1964, 22, 385−400. (22) Vasseur, P.; Cox, R. G. The lateral migration of spherical particles sedimenting in a stagnant bounded fluid. J. Fluid Mech. 1977, 80, 561− 591. (23) McLaughlin, J. B. Inertial migration of a small sphere in linear shear flows. J. Fluid Mech. 1991, 224, 261−274. (24) Miyazaki, K.; Bedeaux, D.; Avalos, J. B. Drag on a sphere in slow shear flow. J. Fluid Mech. 1995, 296, 373−390. (25) Sarvepalli, D.; Schmidtke, D.; Nollert, M. Design consideration for a microfluidic device to quantify the platelet adhesion to collagen at physiological shear rates. Ann. Biomed. Eng. 2009, 37, 1331−1341. (26) ANSYS INC. http://www1.ansys.com/customer/content/ documentation/121/wb2_help.pdf (accessed Nov 2009). ANSYS 10474


Drag Coefficient Correction for Spherical and Nonspherical Particles

Recommend Documents