1. Introduction
Sub-micron pollutants can evade the body’s respiratory defences and deposit in the lung’s alveoli, with exposure linked to stroke, heart disease and cancer (Donaldson et al. Reference Donaldson, Aitken, Tran, Stone, Duffin, Forrest and Alexander2006; Poland et al. Reference Poland, Duffin, Kinloch, Maynard, Wallace, Seaton, Stone, Brown, MacNee and Donaldson2008; Kwon, Ryu & Carlsten Reference Kwon, Ryu and Carlsten2020). But particulate size is not the only issue – mounting evidence suggests that particle shape also plays a critical role in health outcomes (Donaldson et al. Reference Donaldson, Aitken, Tran, Stone, Duffin, Forrest and Alexander2006; Poland et al. Reference Poland, Duffin, Kinloch, Maynard, Wallace, Seaton, Stone, Brown, MacNee and Donaldson2008; Thu et al. Reference Thu, Haider, Khan, Sohail and Hussain2023). This article 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 characterised by the Knudsen number,
${K\hspace {-.05cm}n}={\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$
nanometres. (The Reynolds and Mach numbers are negligible for sub-micron particulates.)
There is a wealth of experimental and theoretical data across the
$K\hspace {-.05cm}n$
scale for the drag on slowly translating spheres (Millikan Reference Millikan1923a
; Sone & Aoki Reference Sone and Aoki1977; Allen & Raabe Reference Allen and Raabe1982, Reference Allen and Raabe1985; Beresnev, Chernyak & Fomyagin Reference Beresnev, Chernyak and Fomyagin1990; Kalempa & Sharipov Reference Kalempa and Sharipov2020) (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 (Bird Reference Bird1994); however, DSMC is extremely computationally intensive, particularly for low-speed, low-
$K\hspace {-.05cm}n$
, external flows (Livi et al. Reference Livi, Di Staso, Clercx and Toschi2022; Clercx et al. Reference Clercx, Livi, Di Staso and Toschi2024; Tatsios et al. Reference Tatsios, Vasileiadis, Gibelli, Borg and Lockerby2025). 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 characterise motion.
Given these challenges, pragmatic means have been developed in the aerosol community to predict particle behaviour across the
$K\hspace {-.05cm}n$
scale. For example, there is a large body of work on predicting the mobility of aggregates of spheres at arbitrary
$K\hspace {-.05cm}n$
(Sorensen Reference Sorensen2011; Corson, Mulholland & Zachariah Reference Corson, Mulholland and Zachariah2017), as these can be a good model for many aerosol particles (e.g. soot). The leading heuristic approach for spanning the
$K\hspace {-.05cm}n$
scale for arbitrary geometries is the adjusted sphere method (ASM) proposed by Dahneke (Reference Dahneke1973b
). The method’s effectiveness has been demonstrated by both experiments (Shapiro et al. Reference Shapiro, Vainshtein, Dutcher, Emery, Stolzenburg, Kittelson and McMurry2012; Thajudeen, Jeon & Hogan Reference Thajudeen, Jeon and Hogan2015) and DSMC simulations (Zhang et al. Reference Zhang, Thajudeen, Larriba, Schwartzentruber and Hogan2012; Livi et al. Reference Livi, Di Staso, Clercx and Toschi2022). The ASM predicts the drag for an arbitrary particle by interpolating between the known ‘continuum limit’ (
${K\hspace {-.05cm}n}=0$
) and the behaviour in ‘free-molecular flow’ (
${K\hspace {-.05cm}n}\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 exist).
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 correctly reproduces the way in which a non-spherical particle’s drag approaches its limiting value as
${K\hspace {-.05cm}n}\!\rightarrow \!\infty$
, it fails to do the same as
${K\hspace {-.05cm}n}\!\rightarrow 0$
; it reproduces the correct limit (the continuum drag force), but not the correct approach to it (discussed more later). Finally, the method can be conceptually confusing; it requires defining a notional sphere and ‘adjusted Knudsen number’ for each unique 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 gap, then, for a simpler and more general heuristic approach to predicting the drag on slow-moving non-spherical particles across the
$K\hspace {-.05cm}n$
scale.
2. Revisiting the Cunningham correction factor
In 1910, Ebenezer Cunningham published a short theoretical article on the settling velocity of spherical particles in a gas (Cunningham Reference Cunningham1910). 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\hspace {-.05cm}n$
.
To this day, the ‘Cunningham correction factor’ is a fundamental tool in a vast array of applications, for example: understanding the transport of atmospheric pollutants, designing filtration systems, modelling particle deposition in the lungs, analysing powder technologies, controlling contaminants and monitoring combustion emissions.
For slow-moving spheres, the drag force (
$D$
) is linearly related to particle speed (
$V$
)
\begin{equation} D= \frac {6 \pi \mu R V }{C} , \end{equation}
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
\begin{equation} C_{\textit{Cunningham}}= 1+{K\hspace {-.05cm}n}\,(\mathcal{A} + B e^{-c/{K\hspace {-.05cm}n}}), \end{equation}
where
$\mathcal{A}$
,
$B$
and
$c$
are fitted coefficients and
${K\hspace {-.05cm}n}=\lambda /R$
.
Although Cunningham is widely credited with (2.2), the exponential term does not feature at all in Cunningham’s article. It was originally proposed by Knudsen & Weber (Reference Knudsen and Weber1911), and only given a theoretical basis/interpretation in the 1920s (Millikan Reference Millikan1923b ).
Cunningham’s correction was originally of a simpler form
\begin{equation} C_{\textit{original}}= 1+\mathcal{B}\,{K\hspace {-.05cm}n} , \end{equation}
where the coefficient
$\mathcal{B}$
can be obtained theoretically by considering the free-molecular limit. Cunningham’s intention was to propose a heuristic expression capable of predicting the particle settling velocity across the
$K\hspace {-.05cm}n$
scale, between the known limit in continuum conditions (i.e. Stokes drag;
$C=1$
) and the correct asymptotic behaviour in the limit of free-molecular flow
\begin{equation} C \sim \mathcal{B}\, {K\hspace {-.05cm}n}\quad \text{as}\quad {K n} \rightarrow \infty . \end{equation}
Over a decade later, Millikan (Reference Millikan1923b
) pointed out a problem with Cunningham’s expression (2.3): it does not reproduce the correct asymptotic behaviour at small Knudsen numbers (
$ C \nsim (1+ \mathcal{B}\, {K\hspace {-.05cm}n}) \hbox{ as } {K\hspace {-.05cm}n} \rightarrow 0$
). The form of the asymptotic approach to the continuum limit is correct, but the coefficient, which can be extracted from the analytical work of Basset (Reference Basset1888) on slip flow, is not equal to
$\mathcal{B}$
. Instead,
\begin{equation} C \sim (1+\mathcal{A}\, \textit{Kn})\quad \text{as}\quad \textit{Kn} \rightarrow 0 . \end{equation}
For this reason, Millikan (Reference Millikan1923b
) proposed a pragmatic modification that would allow both asymptotic results to be captured. Millikan modified the prefactor to
$K\hspace {-.05cm}n$
in Cunningham’s expression (2.3) so that it would smoothly transition from
$\mathcal{A}$
to
$\mathcal{B}$
in the appropriate limits
\begin{equation} C_{\textit{Cunningham}} =1+{K\hspace {-.05cm}n}\,(\mathcal{A}+(\mathcal{B}-\mathcal{A})e^{-c/{K\hspace {-.05cm}n}}), \end{equation}
where
$c$
is a parameter that Millikan described as the ‘rapidity of shift’, which must be fitted to experimental data. Note, (2.2) and (2.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 (Davies Reference Davies1945; Allen & Raabe Reference Allen and Raabe1982, Reference Allen and Raabe1985); and from its original form (2.3) to the one actually due to Knudsen & Weber (Reference Knudsen and Weber1911).
2.1. A new correction factor
Millikan’s sound argument in 1923 was that the Cunningham heuristic needed to be modified so that the asymptotic behaviour of (2.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 (2.4) and (2.5), without need for a fitting parameter
\begin{equation} C_{\textit{ne}w}=\mathrm{e}^{-(\mathcal{B}-\mathcal{A}){K\hspace {-.05cm}n}}+\mathcal{B}{K\hspace {-.05cm}n} . \end{equation}
2.2. Initial test case: spherical particles
The new form of correction factor (2.7), being free of empirical parameters, can serve as a powerful predictive heuristic. Although the focus of this article is non-spherical particles, a sphere is the logical first benchmark for evaluating its accuracy, as reliable data are readily available.
For a sphere, the coefficient
$\mathcal{A}$
can be obtained by solving the Stokes equations with a Maxwell (Reference Maxwell1879) slip boundary condition (as done by Basset Reference Basset1888), with the result taken to first order in
$K\hspace {-.05cm}n$
(as done by Millikan Reference Millikan1923b
). In this way it is easy to show, for the case of the sphere, that
$\mathcal{A}=\beta =(2-\sigma )/\sigma$
, where
$\beta$
is the ‘slip’ coefficient, which Maxwell’s condition gives in terms of
$\sigma$
, the accommodation coefficient (
$\sigma =1$
for diffuse molecular reflection at the particle surface, and
$\sigma =0$
for purely specular reflection). For all the results in this paper, the Maxwell model and diffuse molecular reflection will be adopted (
$\sigma =\beta =1$
). Note, Maxwell derived his slip boundary condition for planar surfaces, but it is applicable to any particle in the limit of
${K\hspace {-.05cm}n}\rightarrow 0$
(the limit at which
$\mathcal{A}$
is evaluated). The coefficient
$\mathcal{B}$
comes directly from the free-molecular result of Epstein (Reference Epstein1924). For diffuse reflection of gas molecules at the sphere surface,
$\mathcal{B}=18/(8+\pi )$
.
Figure 1 compares the prediction of (2.7) with the data from Millikan (Reference Millikan1923b ), relatively recent experiments (Allen & Raabe Reference Allen and Raabe1985) and modern kinetic theory (Beresnev et al. Reference Beresnev, Chernyak and Fomyagin1990). Given (2.7) is heuristic, with no fitting parameters, its prediction is extremely good: the maximum difference between the prediction of (2.7) and the experimental data of Allen & Raabe (Reference Allen and Raabe1985) is less than 4 % of the continuum value.
Drag on a slowly translating sphere against
$K\hspace {-.05cm}n$
. Comparison of Millikan’s data (Millikan Reference Millikan1923b
) (
$+$
), experiments of Allen & Raabe (Reference Allen and Raabe1985) (
$\boldsymbol{\cdots }$
), kinetic theory of Beresnev et al. (Reference Beresnev, Chernyak and Fomyagin1990) (
$\bigcirc$
) and the proposed heuristic (
$1/C_{\textit{ne}w}$
, —), (2.7).
3. Application to non-spherical particles
The forms of the asymptotic limits (2.4) and (2.5) are the same for particles of arbitrary shape. In the continuum limit (
${K\hspace {-.05cm}n}\rightarrow 0$
), a small departure from Stokes’ law can be expressed as a series expansion truncated to first order in
$K\hspace {-.05cm}n$
; a correction factor in the form of (
$1+ \mathcal{A}{K\hspace {-.05cm}n}$
) 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\hspace {-.05cm}n}\rightarrow \infty$
becomes proportional to the density of the gas, and the correction factor proportional to
$K\hspace {-.05cm}n$
; as in (2.4).
The core idea of this article is: if
$\mathcal{A}$
and
$\mathcal{B}$
can be found for a particle, (2.7) can be used to approximate the gas’s resistance to its motion across the entire
$K\hspace {-.05cm}n$
scale – whatever the particle shape.
3.1. 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 (Happel & Brenner Reference Happel and Brenner1983)
\begin{equation} \boldsymbol{f}=- \unicode{x1D646} \boldsymbol{\cdot } \boldsymbol{v} , \end{equation}
where
$\boldsymbol{f}$
is the force on the particle,
$\boldsymbol{v}$
is the particle velocity (relative to the gas) and
$\unicode{x1D646}$
is a 3
$\times$
3 tensor, whose components depend (potentially independently) on
$K\hspace {-.05cm}n$
. It is relatively easy to calculate the resistance tensor in certain limits
\begin{equation} \unicode{x1D646}^{\,0} = \lim _{{K\hspace {-.05cm}n} \to 0}\hspace {-.1cm}\left ( \unicode{x1D646}\,\right ), \quad \unicode{x1D646}^{1} = \lim _{{K\hspace {-.05cm}n} \to 0}\left (\frac {\mathrm{d} \unicode{x1D646}}{\mathrm{d}{K\hspace {-.05cm}n}}\right )\!, \quad \unicode{x1D646}^{\,\infty } = \lim _{{K\hspace {-.05cm}n} \to \infty }\hspace {-.07cm}({K\hspace {-.05cm}n} \unicode{x1D646}\,). \end{equation}
From these, following the spirit of Cunningham, but using (2.7), it is possible to propose a simple element-wise correction to the continuum resistance tensor
\begin{equation} K_{\textit{ij}}=K_{\textit{ij}}^{\mathrm{0}}/C_{\textit{ij}} , \end{equation}
where
$C_{\textit{ij}}=e^{-(\mathcal{B}_{\textit{ij}}-\mathcal{A}_{\textit{ij}}){K\hspace {-.05cm}n}}+\mathcal{B}_{\textit{ij}}{K\hspace {-.05cm}n}$
, and where the correct asymptotic behaviour at both limits is ensured by setting
${\mathcal{A}_{\textit{ij}}} =-K^{\mathrm{1}}_{\textit{ij}}/K^{\mathrm{0}}_{\textit{ij}}$
and
${\mathcal{B}_{\textit{ij}}} = K^{\mathrm{0}}_{\textit{ij}}/K^{\mathrm{\infty }}_{\textit{ij}}$
.
3.2. Evaluating the correction tensor
To evaluate
$ \unicode{x1D63E}$
for a given particle, all that is required is
$ \unicode{x1D646}^{\,\mathrm{0}}$
,
$ \unicode{x1D646}^{\mathrm{1}}$
and
$ \unicode{x1D646}^{\,\mathrm{\infty }}$
. The continuum resistance tensor
$ \unicode{x1D646}^{\,\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 (Pozrikidis Reference Pozrikidis1992) or the method of fundamental solutions (MFS) (Jordan & Lockerby Reference Jordan and Lockerby2025). For some geometries, analytical results exist (for example, the spheroid (Oberbeck Reference Oberbeck1876; Happel & Brenner Reference Happel and Brenner1983); see Appendix A).
The free-molecular resistance tensor,
$ \unicode{x1D646}^{\,\mathrm{\infty }}$
, for arbitrary convex particles, can be obtained by evaluating the following integral over the particle surface,
$S$
(Dahneke Reference Dahneke1973a
; Halbritter Reference Halbritter1974):
\begin{equation} \unicode{x1D646}^{\,\infty } = \frac {\mu }{L} \int _{S} \left ( \frac {\sigma }{2} \unicode{x1D644} + \gamma \, \boldsymbol{n}\boldsymbol{n}\right )\, {\rm d}S , \end{equation}
where
$ \unicode{x1D644}$
is the identity tensor,
$\boldsymbol{n}$
is an outward-facing surface normal,
$\gamma = (8+ \pi \sigma -6\sigma )/4$
and
$\sigma$
is the accommodation coefficient. The analytical result for the spheroid (Dahneke Reference Dahneke1973a
) is included in Appendix A. For more complex convex particles, the surface integral (3.4) can be performed numerically, for example using the MFS (Lockerby Reference Lockerby2022). For non-convex particles, DSMC can be adopted (Gallis, Torczynski & Rader Reference Gallis, Torczynski and Rader2001; Chinnappan et al. Reference Chinnappan, Kumar, Arghode and Myong2019), which becomes very efficient approaching free-molecular flow.
The first-order resistance tensor,
$ \unicode{x1D646}^{\mathrm{1}}$
, is not so straightforward. For a sphere, a simple analytical solution to the Stokes equations with slip boundary conditions exists (Basset Reference Basset1888). This can be truncated to first order in
$K\hspace {-.05cm}n$
, as done by Millikan (Reference Millikan1923a
). 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’ (Keh & Chang Reference Keh and Chang2008).
Fortunately, a very convenient technique has been developed, relatively recently, for evaluating
$ \unicode{x1D646}^{\mathrm{1}}$
directly (Ramachandran & Khair Reference Ramachandran and Khair2009; Stone Reference Stone2010). It exploits the reciprocal theorem (Happel & Brenner Reference Happel and Brenner1983; Masoud & Stone Reference Masoud and Stone2019), and has been employed to derive expressions for the first-order slip/
$K\hspace {-.05cm}n$
effect on drag around spheroids (Masoud & Stone Reference Masoud and Stone2019), Janus spheres (Ramachandran & Khair Reference Ramachandran and Khair2009) and a range of other configurations (Lockerby Reference Lockerby2025). Adapting the form in Lockerby (Reference Lockerby2025) we can write
\begin{equation} \unicode{x1D646}^{\mathrm{1}}=K_{\textit{ij}}^{\mathrm{1}}=-\frac {\beta \, L}{\mu } \int _{S} \boldsymbol{\tau ^{i}}\boldsymbol{\cdot }\boldsymbol{\tau ^{j}}\, \, {\rm d}S , \end{equation}
where
$\boldsymbol{\tau ^i}$
is the surface shear-stress vector generated by a unit particle velocity in the
$i^{{}}$
th direction, from a no-slip solution to the Stokes equations. The closed-form expression for the spheroid (Masoud & Stone Reference Masoud and Stone2019) is provided in Appendix A.
3.3. Results
In the absence of comprehensive experimental data for drag on slow-moving non-spherical particles across the
$K\hspace {-.05cm}n$
scale, DSMC simulations (stochastic solutions to the Boltzmann equation) represent the accepted and most reliable benchmark. Besides DSMC simulations for aggregations of spheres (Zhang et al. Reference Zhang, Thajudeen, Larriba, Schwartzentruber and Hogan2012), only spheroids have been properly studied (Livi et al. Reference Livi, Di Staso, Clercx and Toschi2022; Clercx et al. Reference Clercx, Livi, Di Staso and Toschi2024; Tatsios et al. Reference Tatsios, Vasileiadis, Gibelli, Borg and Lockerby2025; Zhang et al. Reference Zhang, Chang, Wang and Xia2025); the most careful and detailed studies being by Livi et al. (Reference Livi, Di Staso, Clercx and Toschi2022) and Clercx et al. (Reference Clercx, Livi, Di Staso and Toschi2024).
Figure 2 compares data from Clercx et al. (Reference Clercx, Livi, Di Staso and Toschi2024) with (3.3), for drag on various aspect ratio, prolate and oblate spheroids, as a function of
$K\hspace {-.05cm}n$
. 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.
Resistance tensor components for prolate (a,b) and oblate (c,d) spheroids of aspect ratios 4 (a,c) and 10 (b,d). Motion parallel (
$\triangle$
,
$K_{x\hspace {-.01cm}x}$
) and perpendicular (
$\bigcirc$
,
$K_{y\hspace {-.01cm}y}$
) to the polar axis. Comparison of DSMC (Clercx et al. Reference Clercx, Livi, Di Staso and Toschi2024) with (3.3).
Figure 3 compares (3.3) with the data of Livi et al. (Reference Livi, Di Staso, Clercx and Toschi2022) for aspect ratio 2 spheroids, at various
$K\hspace {-.05cm}n$
, as a function of orientation. The benefit of the resistance tensor description (3.1), is that the drag for any particle orientation can be determined by a simple transformation:
$ \unicode{x1D646}'= \unicode{x1D64D}\boldsymbol{\cdot }\unicode{x1D646} \boldsymbol{\cdot }\unicode{x1D64D}^\top$
, where
$ \unicode{x1D64D}$
is a rotation matrix. For figure 3, the flow direction (
$x'$
) is rotated from the spheroid’s axis of revolution (
$x$
) by a single-axis rotation
$\theta$
about
$z$
:
$ \unicode{x1D64D}= \unicode{x1D64D}_z(\theta )$
. Again, the general level of agreement is very good, with the greatest discrepancy at
${K\hspace {-.05cm}n}=1$
, where the heuristic is farthest from either limit and DSMC is hardest to perform accurately.
Resistance component in the direction of motion (
$x'$
) for prolate (a) and oblate (b) spheroids of aspect ratio 2, against
$\theta$
, the angle (in degrees) between
$x'$
and
$x$
(where
$x$
is the spheroid’s axis of revolution). Shown are DSMC data (Livi et al. Reference Livi, Di Staso, Clercx and Toschi2022) at
${K\hspace {-.05cm}n}=1(\square ),\,5(\triangle ),\,7(\Diamond ),\,9($
$\bigcirc$
$)\,$
and
$10(\triangleleft )$
; (3.3) (—).
The final test case involves predicting the drag on an infinitely thin circular disc of radius
$R$
, moving perpendicular to its surface along the
$x$
-axis; recent kinetic-theory results, based on the Bhatnagar-Gross-Krook (BGK) model, provide a benchmark for comparison across the entire
$K\hspace {-.05cm}n$
-scale (Tomita, Taguchi & Tsuji Reference Tomita, Taguchi and Tsuji2025). The disc is the high-aspect-ratio limit of an oblate spheroid, so, again, we can use the results in Appendix A (with
$e=1$
and
$b=R$
). This case is particularly challenging for the heuristic, because the Maxwell slip model assumes a planar surface on the scale of the mean free path, and the disc’s edge has infinite curvature. Given this, the agreement shown in figure 4 is remarkably good; the maximum difference is 4 % of the continuum resistance value.
4. Discussion and future work
The results in figures 1–4 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 Appendix A, 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\lt {K\hspace {-.05cm}n}\lt 10$
). It is only an approximation, and no replacement for future, much-needed, experiments and Boltzmann-equation solutions for non-spherical particles.
In this work, the Maxwell slip boundary condition has been adopted throughout. Modern kinetic theory predicts a slightly higher value for
$\beta$
than Maxwell’s model (for
$\sigma =1$
), typically in the range
$\beta =1.11$
–
$1.15$
depending on the specific molecular-collision model (Sharipov Reference Sharipov2004). However, the results presented here are not particularly sensitive to the precise value of
$\beta$
chosen. For example, if
$\beta =1.13$
(a working average proposed by Young (Reference Young2011)) the maximum difference in the results presented in figures 1, 2 and 4 would be around 3 % of the corresponding continuum value; occurring at
${K\hspace {-.05cm}n}\approx 1$
.
The predictions for spheroids have not been tested in the low transition regime, due to the compounded difficulties in performing DSMC there (Tatsios et al. Reference Tatsios, Vasileiadis, Gibelli, Borg and Lockerby2025). It is fair to point out that the 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 cannot ensure
$\lim _{{K\hspace {-.05cm}n}\rightarrow 0}(\mathrm{d} \unicode{x1D646}/\mathrm{d}{K\hspace {-.05cm}n})= \unicode{x1D646}^{\mathrm{1}}$
, which (3.3) does by design. For an oblate spheroid of aspect ratio 3, the ASM underpredicts the steepness of the
$K_{xx}({K\hspace {-.05cm}n})$
curve at
${K\hspace {-.05cm}n}=0$
by approximately 20 %; this increases to nearly 30 % for aspect ratio 10. It is also fair to point out here that Livi et al. (Reference Livi, Di Staso, Clercx and Toschi2022) and Clercx et al. (Reference Clercx, Livi, Di Staso and Toschi2024) 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 (2.7) assumes that
$\mathcal{B}-\mathcal{A}\gt 0$
. This is the case for all the results presented in this article, but not true in general. For spheroids, it is true for all possible aspect ratios and orientations; except for very high-aspect-ratio oblate spheroids when travelling perpendicular to their axes of revolution (
$b/a \gtrsim 45$
). The quantity
$\mathcal{B}-\mathcal{A}$
can be thought of as reflecting the compactness of the transition regime: when it becomes negative, the separation between the near-continuum (2.5) and free-molecular (2.4) limits is very large, and both asymptotes cannot be satisfied using a correction factor of the form in (2.7). In the event that
$\mathcal{B}\lt \mathcal{A}$
, a practical fix is to set the exponent (
$\mathcal{B}-\mathcal{A}$
) to zero. This will sacrifice accurate asymptotic behaviour as
${K\hspace {-.05cm}n}\rightarrow 0$
for a sensible correction across the
$K\hspace {-.05cm}n$
scale.
Brenner proved that the continuum resistance tensor (
$ \unicode{x1D646}^{\,\mathrm{0}}$
) is symmetric (Happel & Brenner Reference Happel and Brenner1983). Quick inspection of (3.4) and (3.5) reveals that
$ \unicode{x1D646}^{\mathrm{1}}$
and
$ \unicode{x1D646}^{\,\mathrm{\infty }}$
are also symmetric; it follows, then, that
$\boldsymbol{C}$
, and any predicted resistance tensor, will also be symmetric. As discussed in Happel & Brenner (Reference Happel and Brenner1983), an alternative proof exists for the symmetry of
$ \unicode{x1D646}^{\,\mathrm{0}}$
due to Landau & Lifshitz (Reference Landau and Lifshitz1958) that does not 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\hspace {-.05cm}n$
scale, as predicted by (3.3). One of the consequences of this is that the ‘sine-squared drag law’ should hold at all
$K\hspace {-.05cm}n$
, which is supported by DSMC simulations (Bird Reference Bird1994; Livi et al. Reference Livi, Di Staso, Clercx and Toschi2022; Clercx et al. Reference Clercx, Livi, Di Staso and Toschi2024; Zhang et al. Reference Zhang, Chang, Wang and Xia2025).
Resistance tensor component in the direction of motion (
$x$
) for an infinitely thin circular disc of radius
$R$
(
$=L$
) moving perpendicular to its surface. Comparison of the BGK-Boltzmann solution of Tomita et al. (Reference Tomita, Taguchi and Tsuji2025) (
$\bigcirc$
) and (3.3) (—).
The form of the correction tensor tacitly assumes that the direction of the principal axes remains constant across the entire
$K\hspace {-.05cm}n$
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 (Brenner Reference Brenner1965). Such particles will almost certainly require numerical integration and full Stokes-flow simulations to evaluate the required resistance tensors in each limit. Preliminary investigation suggests that the rotational resistance matrix may behave similarly to the translational one in the three limits of (3.2), so the same heuristic could be used. The coupling tensor is more complicated, however, and its components may even switch sign between the continuum and free-molecular limits.
Data availability
A short Matlab subroutine for evaluating the spheroid resistance tensor as a function of
$K\hspace {-.05cm}n$
, used in figures 2–4, is downloadable at https://doi.org/10.6084/m9.figshare.30264322.
Declaration of interests
The author reports no conflict of interest.
Appendix A.
$ \unicode{x1D646}^{\,\boldsymbol{0}}$
,
$ \unicode{x1D646}^{\kern1pt\boldsymbol{1}}$
and
$ \unicode{x1D646}^{\,\boldsymbol{\infty}}$
for spheroids
The continuum (Happel & Brenner Reference Happel and Brenner1983), first-order (Masoud & Stone Reference Masoud and Stone2019) and free-molecular (Dahneke Reference Dahneke1973a
) resistance tensors for a spheroid defined by
$(x/a)^2 + (y/b)^2 + (z/b)^2 = 1$
are
\begin{align*} \unicode{x1D646}^{\,\mathrm{0}} &=16\,\pi \mu \,b \,e^3 \left ( \frac {\boldsymbol{i}_x \boldsymbol{i}_x}{2 e\, r+\left (4 e^2-2\right ) \sin ^{-1}(e)} + \frac {\boldsymbol{i}_y\boldsymbol{i}_y+\boldsymbol{i}_z \boldsymbol{i}_z}{(2\,e^2 + 1)\,\mathrm{sin}^{-1}(e) -e\,r} \right ) \! , \nonumber\\ \unicode{x1D646}^{\mathrm{1}} &=-32\,\pi \mu L \,e^3 \beta \left (\boldsymbol{i}_x \boldsymbol{i}_x \frac {e-r^2 \,\mathrm{tanh}^{-1}(e)}{2\,{{\left (e\,r+\mathrm{sin}^{-1}(e)\,{\left (2\,e^2 -1\right )}\right )}}^2 } \right . \nonumber\\ &\quad \left .+( \boldsymbol{i}_y\boldsymbol{i}_y+\boldsymbol{i}_z \boldsymbol{i}_z) \frac {\mathrm{tanh}^{-1}(e)\,{\left (e^2 +1\right )-e}}{{{\left (e\,r-\mathrm{sin}^{-1}(e)\,{\left (2\,e^2 +1\right )}\right )}}^2 }\right ) \hbox{ and} \\ \unicode{x1D646}^{\,\mathrm{\infty }} &=\frac {\mu }{L} \frac {\pi b^2}{e^3} \biggl (\boldsymbol{i}_x \boldsymbol{i}_x \left [ \sigma e^3 + 2e \gamma + r^2 (\sigma e^2 - 2\gamma ) \,\mathrm{tanh}^{-1}(e)\right ]\nonumber \\ & \quad +(\boldsymbol{i}_y\boldsymbol{i}_y+\boldsymbol{i}_z \boldsymbol{i}_z)\Bigl [-\gamma e r^2 + e^3 \sigma + r^2 (\gamma (1 + e^2) + e^2 \sigma ) \mathrm{tanh}^{-1}(e)\Bigr ] \biggr ) , \end{align*}
where
$\boldsymbol{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).
Open access
