A correction tensor for approximating drag on slow-moving particles of arbitrary shape

TL;DR

This paper introduces a novel correction tensor for accurately estimating drag on slow-moving, arbitrarily shaped particles without experimental fitting, validated against theoretical and experimental data.

Contribution

It presents a new correction tensor that generalizes the Cunningham correction factor for non-spherical particles, enabling better predictions of their transport properties.

Findings

Accurately predicts drag on non-spherical particles.

Validated against experiments and kinetic theory.

Provides a simple, general tool for micro/nano particle transport.

Abstract

A new form of the Cunningham correction factor is presented that requires no experimental fitting. It is expanded to provide a predictive heuristic for non-spherical particles, via definition of a "correction tensor''. Its accuracy is tested against experiments and kinetic theory for the sphere, and stochastic solutions to the Boltzmann equation for a range of spheroids. It represents a simple, general tool for approximating transport properties of non-spherical micro/nano particles in a gas.

Equations

$$D=\frac{6\pi \mu RV}{C}$$

$$C_{\mathrm{Cunningham}}=1+Kn(\mathcal{A}+B{e}^{-c/Kn})$$

$$C_{\mathrm{original}}=1+\mathcal{B}Kn,$$

$$C\sim \mathcal{B}KnasKn\to \infty .$$

$$C\sim (1+\mathcal{A}Kn)asKn\to 0.$$

$$C_{\mathrm{Cunningham}}=1+Kn(\mathcal{A}+(\mathcal{B}-\mathcal{A})e^{-c/Kn}),$$

$$C_{\mathrm{new}}=e^{-(\mathcal{B}-\mathcal{A})Kn}+\mathcal{B}Kn.$$

$$\mathbf{f}=-\mathbf{K}\cdot \mathbf{v}$$

$${\mathbf{K}}^0,{\mathbf{K}}^1=\underset{Kn\to 0}{\lim }\left(\mathbf{K},\frac{\mathrm{d}\mathbf{K}}{\mathrm{d}Kn}\right),{\mathbf{K}}^{\infty }=\underset{Kn\to \infty }{\lim }(Kn\mathbf{K}).$$

$$K_{ij}=K_{ij}^0/C_{ij}$$

$${\mathbf{K}}^{\infty }=\frac{\mu }{L}{\int }_S\left(\frac{\sigma }{2}\mathbf{I}+\gamma \mathbf{n}\mathbf{n}\right)dS$$

$${\mathbf{K}}^1=K_{ij}^1=-\frac{L}{\mu }{\int }_S{\mathbf{\tau }}^{\mathbf{i}}\cdot {\mathbf{\tau }}^{\mathbf{j}}dS,$$

$${\mathbf{K}}^0=16\pi \mu b{e}^3\left(\frac{{\mathbf{i}}_x{\mathbf{i}}_x}{2er+\left(4e^2-2\right){\sin }^{-1}(e)}+\frac{{\mathbf{i}}_y{\mathbf{i}}_y+{\mathbf{i}}_z{\mathbf{i}}_z}{(2e^2+1){\sin }^{-1}(e)-er}\right),$$

$${\mathbf{K}}^1=-32\pi \mu L{e}^3\beta \left({\mathbf{i}}_x{\mathbf{i}}_x\frac{e-r^2{\tanh }^{-1}(e)}{2{\left(er+{\sin }^{-1}(e)\left(2e^2-1\right)\right)}^2}+({\mathbf{i}}_y{\mathbf{i}}_y+{\mathbf{i}}_z{\mathbf{i}}_z)\frac{{\tanh }^{-1}(e)\left(e^2+1\right)-e}{{\left(er-{\sin }^{-1}(e)\left(2e^2+1\right)\right)}^2}\right)and$$

$$\displaystyle\bm{K}^{\mathrm{\infty}}=\frac{\mu}{L}\frac{\pi b^{2}}{e^{3}}\biggl{(}\bm{i}_{x}\bm{i}_{x}\left[\sigma e^{3}+2e\gamma+r^{2}(\sigma e^{2}-2\gamma)\,\mathrm{tanh}^{-1}(e)\right]$$

$$\displaystyle\hskip 56.9055pt+(\bm{i}_{y}\bm{i}_{y}+\bm{i}_{z}\bm{i}_{z})\Bigl{[}-\gamma er^{2}+e^{3}\sigma+r^{2}(\gamma(1+e^{2})+e^{2}\sigma)\mathrm{tanh}^{-1}(e)\Bigr{]}\biggr{)}\,,$$

Full text

††preprint: APS/123-QED

Introduction – Sub-micron pollutants can evade the body’s respiratory defenses and deposit in the lung’s alveoli, with exposure linked to stroke, heart disease, and cancer [1, 2, 3]. But particulate size is not the only issue — mounting evidence suggests that particle shape also plays a critical role in health outcomes [4, 2, 5, 3]. This Letter proposes a simple and general method for predicting the transport properties of such slow-moving non-spherical sub-micron particles.

The scale of a particulate is characterized by the Knudsen number, ${K\hskip-1.42271ptn}={\lambda}/{L}$, a dimensionless quantity representing the ratio of the gas’s mean free path ($\lambda$) to the particle’s characteristic length scale ($L$). For air at sea level, $\lambda\approx 65$ nanometers. (The Reynolds and Mach numbers are negligible for sub-micron particulates).

There is a wealth of experimental and theoretical data across the ${K\hskip-1.42271ptn}$ scale for the drag on slowly translating spheres [6, 7, 8, 9, 10, 11] (which is related to their diffusivity via the Stokes-Einstein relation), but none, or very little, for the drag on non-spherical particles. Experiments are extremely challenging, and numerical simulations must involve solving the Boltzmann equation, or some approximation to it. Some progress in this has been made using the direct-simulation Monte Carlo (DSMC) method [12]; however, DSMC is extremely computationally intensive, particularly for low-speed, low-${K\hskip-1.42271ptn}$, external flows [12, 13, 14, 15]. Furthermore, studying non-spherical particles not only expands the parameter space but also complicates the analysis of resistance, as a single drag coefficient is insufficient. Instead, a “resistance tensor” must be determined to characterize motion.

Given these challenges, pragmatic means have been developed in the aerosol community to predict particle behavior across the ${K\hskip-1.42271ptn}$ scale. For example, there is a large body of work on predicting the mobility of aggregates of spheres at arbitrary ${K\hskip-1.42271ptn}$ [16, 17], as these can be a good model for many aerosol particles (e.g. soot).

A popular heuristic approach for spanning the ${K\hskip-1.42271ptn}$ scale, applied to arbitrary geometries (not just agglomerates), is the Adjusted Sphere Method (ASM) proposed by Dahneke [18] in 1973. The method’s effectiveness has been demonstrated by both experiments [19, 20, 21, 22] and DSMC [13, 23]. The ASM predicts drag on an arbitrary particle by interpolating between the known “continuum limit” (${K\hskip-1.42271ptn}=0$) and the behavior in “free-molecular flow” (${K\hskip-1.42271ptn}\rightarrow\infty$). The main assumption being that, when scaled appropriately, the drag’s transition between these limits, for any given particle, mirrors that of a simple sphere (for which reliable data exists).

Despite its success, ASM has significant drawbacks. The model is semi-empirical and based on results for a translating sphere, which makes its application to very non-spherical shapes and other resistance coefficients problematic. Furthermore, while the method accurately reproduces the asymptotic behavior for a non-spherical particle as ${K\hskip-1.42271ptn}\!\rightarrow\!\infty$, it fails to do so as ${K\hskip-1.42271ptn}\!\rightarrow 0$; it reproduces the correct limit value (the continuum one), but not the overall asymptotic behavior (discussed later). Finally, the method can be conceptually confusing; if used correctly, it requires defining a notional sphere and “adjusted Knudsen number” for each resistance coefficient. For a completely arbitrary geometry, with a full resistance tensor, this would result in a matrix of spheres and Knudsen numbers representing the same particle.

There is a need, then, for a simpler and more general heuristic approach for non-spherical particles.

The History of the Cunningham Correction Factor – In 1910, Ebenezer Cunningham published a short theoretical article on the settling velocity of spherical particles in a gas [24]. The article’s lasting legacy is an empirically-fitted expression that describes how drag on a spherical particle departs from Stokes law for non-zero ${K\hskip-1.42271ptn}$.

To this day, the “Cunningham correction factor” is a fundamental tool in aerosol science and used in a vast array of associated technologies. For example, to understand the transport of atmospheric pollutants, for the design of filtration systems, in powder technologies, in lung deposition modeling, in contaminant control, in monitoring combustion emissions, and more.

For slow-moving spherical particles, the drag force ($D$) is linearly related to particle speed ($V$):

[TABLE]

where $\mu$ is the dynamic viscosity, $R$ is the radius of the sphere, and $C$ is a correction factor, the form of which (normally attributed to Cunningham) is

[TABLE]

where $\mathcal{A}$, $B$ and $c$ are fitted coefficients and ${K\hskip-1.42271ptn}=\lambda/R$.

Although Cunningham is widely credited with (2), the exponential term doesn’t feature at all in Cunningham’s article. It was originally proposed by Knudsen and Weber [25], and only given a theoretical basis/interpretation in the 1920s by Millikan [26].

Cunningham’s correction was originally of a simpler form:

[TABLE]

where the coefficient $\mathcal{B}$ can be obtained theoretically. Cunningham’s intention was to propose a heuristic expression capable of predicting particle settling velocity across the ${K\hskip-1.42271ptn}$ scale, between the known limit in continuum conditions (i.e. Stokes drag; $C=1$) and the correct asymptotic behavior in the limit of free-molecular flow:

[TABLE]

Millikan [26], over 10 years later, pointed out a failing in Cunningham’s expression (3), in that it did not reproduce the correct asymptotic behavior at small Knudsen numbers ($C\nsim(1+\mathcal{B}\,{K\hskip-1.42271ptn})\hbox{ as }{K\hskip-1.42271ptn}\rightarrow 0$). It is correct in form, but the coefficient, which can be extracted from the analytical work of Basset [27] on slip flow, is not equal to $\mathcal{B}$. Instead,

[TABLE]

For this reason, Millikan [26] proposed a pragmatic modification that would allow both asymptotic behaviors to be captured. Millikan modified the prefactor to ${K\hskip-1.42271ptn}$ in Cunningham’s expression (3) so that it would smoothly transition from $\mathcal{A}$ to $\mathcal{B}$ in the appropriate limits:

[TABLE]

where $c$ is a parameter that Millikan described as the “rapidity of shift”, which must be fitted to experimental data. Note, (2) and (6) are equivalent to each other; $B=\mathcal{B}-\mathcal{A}$.

It was Millikan’s intervention then, that changed the course of history for the Cunningham correction factor: from its origin as a predictive heuristic to its modern use as a fitting function [28, 7, 8]; and from its original form (3) to the one actually due to Knudsen and Weber [25]. Millikan’s analysis was partly an attempt to claim some intellectual ownership of Knudsen and Weber’s equation by providing a physical basis for it. It is ironic, then, that neither Millikan, Knudsen nor Weber typically gets the credit.

A New Correction Factor – Millikan’s sound argument in 1923 was that the Cunningham heuristic needed to be modified so that the analytical limit of (5) could also be captured. However, for this purpose, his introduction of two new parameters into Cunningham’s expression (to satisfy one more condition), is needlessly complicated. It is also less general, of course, because the spare coefficient $c$ needs to be fitted to experimental data.

Millikan did not spot (or at least did not mention) that a simpler expression exists that satisfies both (4) and (5), without need for a fitting parameter:

[TABLE]

Test Case: Spherical Particles — The new form of correction factor (7), being free of empirical parameters, can serve as a powerful predictive heuristic. Although the focus of this Letter is non-spherical particles, a sphere is the logical first benchmark for evaluating its accuracy, as reliable data is readily available.

For a sphere, the coefficient $\mathcal{A}$ can be obtained by solving the Stokes equations with a Maxwell slip boundary condition [29] (as done by Basset [27]), with the result taken to first order in ${K\hskip-1.42271ptn}$ (as done by Millikan [26]). In this way it is easy to show, for the case of the sphere, that $\mathcal{A}=\beta$ where $\beta$ is the boundary condition’s ‘slip’ coefficient. Note, Maxwell derived his slip boundary condition for planar surfaces, but it’s applicable to any particle in the limit of ${K\hskip-1.42271ptn}\rightarrow 0$ (the limit at which $\mathcal{A}$ is evaluated). For diffuse molecular reflection, Maxwell found $\beta=1$, but modern kinetic theory predicts a slightly higher value, typically in the range $\beta=1.11$–$1.15$, depending on the specific molecular-collision model [30]. Young proposed a working average of $\beta=1.13$ [31], which will be adopted here and for the rest of this Letter. The coefficient $\mathcal{B}$ comes directly from Epstein’s free-molecular result [32]. For diffuse reflection of gas molecules at the sphere surface, $\mathcal{B}=18/(8+\pi)$.

Figure 1 compares the prediction of equation (7) with the data from Millikan [26], relatively recent experiments [8], and modern kinetic theory [10]. Given equation (7) is heuristic, with no fitting parameters, its prediction is exceptionally good.

Application to Non-spherical Particles – The forms of the asymptotic limits (4) and (5) are the same for particles of arbitrary shape. In the continuum limit (${K\hskip-1.42271ptn}\rightarrow 0$), a small departure from Stokes law can be expressed as a series expansion truncated to first order in ${K\hskip-1.42271ptn}$; a correction factor in the form of ($1+\mathcal{A}{K\hskip-1.42271ptn}$) follows. In the free-molecular limit, there are no gas-molecule interactions (i.e. no collisions): the effect of a gas molecule on the resistance to particle motion is independent of any other. Given the low particle speed (relative to the average gas molecule), the resistance generated by the gas (for a given particle motion) as ${K\hskip-1.42271ptn}\rightarrow\infty$ becomes proportional to the density of the gas, and the correction factor proportional to ${K\hskip-1.42271ptn}$; as in (4).

The core idea of this Letter is: if $\mathcal{A}$ and $\mathcal{B}$ can be found for a particle, equation (7) can be used to approximate the gas’s resistance to its motion across the entire ${K\hskip-1.42271ptn}$ scale — whatever the particle shape.

A Correction Tensor – When a particle moves slowly enough that the gas flow responds linearly to its motion, the resistance to the particle’s translation is most generally described by a resistance tensor [33]:

[TABLE]

where $\bm{f}$ is the force on the particle, $\bm{v}$ is the particle velocity (relative to the gas), and $\bm{K}$ is a 3$\times$3 tensor, whose components depend (potentially independently) on ${K\hskip-1.42271ptn}$. It is relatively easy to calculate the resistance tensor in certain limits:

[TABLE]

From these, following the spirit of Cunningham, but using equation (7), it is possible to propose a simple element-wise correction to the continuum resistance tensor:

[TABLE]

where $C_{ij}=e^{-(\mathcal{B}_{ij}-\mathcal{A}_{ij}){K\hskip-1.42271ptn}}+\mathcal{B}_{ij}{K\hskip-1.42271ptn}$, and where the correct asymptotic behavior at both limits is ensured by setting ${\mathcal{A}_{ij}}=-K^{\mathrm{1}}_{ij}/K^{\mathrm{0}}_{ij}$ and ${\mathcal{B}_{ij}}=K^{\mathrm{0}}_{ij}/K^{\mathrm{\infty}}_{ij}$.

Evaluating the Correction Tensor — To evaluate $\bm{C}$ for a given particle, all that is required is $\bm{K}^{\mathrm{0}}$, $\bm{K}^{\mathrm{1}}$ and $\bm{K}^{\mathrm{\infty}}$. The continuum resistance tensor $\bm{K}^{\mathrm{0}}$ can be determined by solving the Stokes equations with a no-slip boundary condition; for complex geometries this can be done, for example, using the boundary element method [34] or the Method of Fundamental Solutions (MFS) [35]. For some geometries, analytical results exist (for example, the spheroid [36, 33]; see the appendix).

The free-molecular resistance tensor, $\bm{K}^{\mathrm{\infty}}$, for arbitrary convex particles, can be obtained by evaluating the following integral over the particle surface, $S$ [37, 38]:

[TABLE]

where $\bm{I}$ is the identity tensor, $\bm{n}$ is an outward-facing surface normal, $\gamma=\left(8+\pi\sigma-6\sigma\right)/4$ and $\sigma$ is the accommodation coefficient (typically $\sigma\approx 1$). The analytical result for spheroids [37] is included in the appendix. For more complex convex particles, the surface integral (11) can be performed numerically (e.g., using the MFS [39]). For non-convex particles, DSMC can be adopted [40, 41] which is very efficient in free-molecular flow.

The first-order resistance tensor, $\bm{K}^{\mathrm{1}}$, is not so straightforward. For a sphere, a simple analytical solution to the Stokes equations with slip boundary conditions exists [27]. This can be truncated to first order in ${K\hskip-1.42271ptn}$, as Millikan did [6]. However, for non-spherical particles, such solutions generally do not exist. Even the spheroid slip solution requires an “infinite-series form of semi-separation of variables” [42].

Fortunately, a very convenient technique has been developed, relatively recently, for evaluating $\bm{K}^{\mathrm{1}}$ directly [43, 44]. It exploits the reciprocal theorem [33], and has been employed to derive expressions for the first-order slip/${K\hskip-1.42271ptn}$ effect on drag around spheroids [45], Janus spheres [43]; and a range of other configurations [46]. Adapting the form in [46] we can write:

[TABLE]

where $\bm{\tau^{i}}$ is the surface shear-stress vector generated by a unit particle velocity in the $i^{\mathrm{th}}$ direction, from a no-slip solution to the Stokes equations. Masoud and Stone’s closed-form expression for the spheroid [45] is provided in the appendix.

Results – In the absence of comprehensive experimental data for drag on slow-moving non-spherical particles across the ${K\hskip-1.42271ptn}$ scale, DSMC simulations (stochastic solutions to the Boltzmann equation) represent the accepted and most reliable benchmark. Besides DSMC simulations for aggregations of spheres [23], only spheroids have been properly studied [13, 14, 15, 47]; the most careful and detailed being by Livi et al. [13] and Clercx et al. [14].

Figure 2 compares data from Clercx et al. [14] to equation (10), for drag on various aspect ratio, prolate and oblate spheroids, as a function of ${K\hskip-1.42271ptn}$. Here the spheroid is defined by $(x/a)^{2}+(y/b)^{2}+(z/b)^{2}=1$, where $x$ is along its axis of revolution; for all cases $L=\sqrt[3]{ab^{2}}$. The general level of agreement in Figure 2 is excellent.

Figure 3 compares (10) with the data of Livi et al. [13] for aspect ratio 2 spheroids, at various ${K\hskip-1.42271ptn}$, as a function of orientation. The benefit of the resistance tensor description (8), is that the drag for any particle orientation can be determined by a simple transformation: $\bm{K}^{\prime}=\bm{R}\cdot\bm{K}\cdot\bm{R}^{\top}$, where $\bm{R}$ is a rotation matrix. For Figure 3, the flow direction ($x^{\prime}$) is rotated from the spheroid’s axis of revolution ($x$) by a single-axis rotation $\theta$ about $z$: $\bm{R}=\bm{R}_{z}(\theta)$. Again, the general level of agreement is very good, with the greatest discrepancy at ${K\hskip-1.42271ptn}=1$, where the heuristic is farthest from either limit and DSMC is hardest to perform accurately.

Discussion and Future Work – The results in Figures 1-3 demonstrate that the heuristic proposed has considerable predictive power. Importantly, there are no fitting parameters: the predictions are all based on known analytical results (the ones for a spheroid are provided in the appendix, for convenience). However, it is also important to stress that the prediction is based on an interpolation across the highly complex “transition regime” ($0.1<{K\hskip-1.42271ptn}<10$). It is only an approximation, and no replacement for future, much-needed, experiments and Boltzmann-equation solutions for non-spherical particles.

The predictions have not been tested in the low transition regime, due to the compounded difficulties in performing DSMC there [15]. It is fair to point out that the Adjusted Sphere Method (ASM) also performs well in the high transition regime. However, a disadvantage of the ASM, other than its conceptual issues and need for empirical input, is that it can’t ensure $\lim_{{K\hskip-1.42271ptn}\rightarrow 0}(\mathrm{d}\bm{K}/\mathrm{d}{K\hskip-1.42271ptn})=\bm{K}^{\mathrm{1}}$, which equation (10) does by design. For an oblate spheroid of aspect ratio 3, the ASM underpredicts the steepness of the $K_{xx}({K\hskip-1.42271ptn})$ curve at ${K\hskip-1.42271ptn}=0$ by about 20%; this increases to nearly 30% for aspect ratio 10. It is also fair to point out here that Livi et al. [13] and Clercx et al. [14] also proposed useful interpolation schemes through their data, but these used spheroid DSMC data for fitting, and are thus only applicable to the spheroid.

The correction factor (7) assumes that $\mathcal{B}-\mathcal{A}>0$. This is the case for all the results presented here, but not true in general. In the event that $\mathcal{B}<\mathcal{A}$, a practical fix is to set the exponent ($\mathcal{B}-\mathcal{A}$) to zero. This will sacrifice accurate asymptotic behavior as ${K\hskip-1.42271ptn}\rightarrow 0$ for a sensible correction across the ${K\hskip-1.42271ptn}$ scale.

Brenner proved that the continuum resistance tensor ($\bm{K}^{\mathrm{0}}$) is symmetric [33]. Quick inspection of (11) and (12) reveals that $\bm{K}^{\mathrm{1}}$ and $\bm{K}^{\mathrm{\infty}}$ are also symmetric; it follows, then, that $\bm{C}$, and any predicted resistance tensor, will also be symmetric. As discussed in [33], an alternative proof exists for the symmetry of $\bm{K}^{\mathrm{0}}$ due to Landau and Lifshitz [48, 49] that doesn’t require invocation of the fluid equations at all, and is purely based on thermodynamic considerations. It thus seems likely that the resistance tensor retains symmetry across the ${K\hskip-1.42271ptn}$ scale as predicted by (10). One of the consequences of this is that the “sine-squared drag law” should hold at all ${K\hskip-1.42271ptn}$, which is also supported by DSMC [13, 14, 47, 12].

The form of the correction tensor tacitly assumes that the direction of the principal axes remains constant across the entire ${K\hskip-1.42271ptn}$ scale. This assumption may not hold true, and the extent of its validity requires further investigation. Future work should also include applying the heuristic correction to the full 6$\times$6 resistance tensor for the study of particles with rotational-translational coupling [50].

Data availability – A short script for evaluating the spheroid resistance tensor as a function of ${K\hskip-1.42271ptn}$, used in Figures 2 and 3, will be made available on publication.

Appendix A $\bm{K}^{\mathrm{0}}$, $\bm{K}^{\mathrm{1}}$ and $\bm{K}^{\mathrm{\infty}}$ for spheroids

The continuum [33], first-order [45], and free-molecular [37] resistance tensors for a spheroid defined by $(x/a)^{2}+(y/b)^{2}+(z/b)^{2}=1$ are

[TABLE]

where $\bm{i}_{x,y,z}$ are unit vectors, $r=a/b$, and $e=\sqrt{1-r^{2}}$ (which is real for oblate spheroids and imaginary for prolate).