Computation of drag and diffusion coefficient for coronavirus: I.

Monte Carlo simulations and integral equation techniques allow for the flexible and efficient computation of drag and diffusion coefficients for virus mimetic particles. We highlight a Monte Carlo method that is useful for computing the drag on biomimetic particles in the free-molecular regime and a numerical technique to compute a (boundary integral solution) of the Stokes equation in the hydrodynamic limit. The free-molecular and the continuum results allow the construction of an approximation for the drag applicable over the full range of Knudsen numbers. Finally, we outline how this work and will be useful in modeling viral transport in air and fluids and in viral morphology measurements and viral separations via electrospray-differential mobility analyzers (ES-DMA).

Introduction

The coronavirus is believed to be mainly transported in air as embedded in drops and droplets (see an earlier study on pathogen transport in air by Vejerano and Marr (2018) and a recent study by Lieber (2021). The evaporation and the drying of these drops and droplets can lead to viral (or clumps of them) transport in the later stages or through filters, masks, and human respiratory tract as partially wet or dry.

Although convective processes are primarily responsible for transporting viruses and other bioaerosols in bulk, diffusion plays a crucial role in their deposition and capture on surfaces. The morphology of the coronavirus and other bioaerosols (such as pollen) are quite complex; thus, analytical solutions or approximations for the drag coefficient or diffusion coefficients of such particles are of interest. Additionally, greater understanding of the impact complex shape/morphology has on diffusion also applies to related areas of study involving general (bio) particle transport.

An area of considerable current interest, is the use of electrospray differential mobility analyzers, referred to as ES-DMA (Allmaier, 2011; Fernández-García, 2019; Hogan, 2006; Pease, 2012; Pease, 2011; Perez-Lorenzo and Fernandez de la Mora, 2021; Weiss, 2019). In this technique, a fine spray from a liquid-virus solution is introduced into an ionizing chamber and then into a flowing gaseous system (air or other gasses at atmospheric pressure). As the droplets change size, the bioparticles accumulate charge and begin to approximate a Boltzmann distribution; and sensors can then measure the bioparticles mobility (DMA) and concentration e.g., with a Condensation Nuclei Counter. The technique is useful for virus sizes between 20 and 800 nm, and the mobility measurements can be quite precise. Since the mobility depends on the viral structure, the technique allows for both structural determination and separation of the viruses.

The inverse problem of computing viral structure from mobility measurements requires knowledge of mobility dependence on the structure itself. The diffusion coefficient, is related to the particle mobility through the well-known Einstein relation. The mobility itself can be obtained from the calculation or measurements of drag force on a particle, and this has been an area of considerable interest in aerosol science and engineering (Cercignani et al., 1968; Friedlander, 2000; Fuchs, 1964; Ivchenko et al., 2007; Lea & Loyalka, 1982; Loyalka, 1992; Takata et al., 1993; Williams, 1973; Williams & Loyalka, 1991). Although early work focused on spherical particles and simple non-spherical geometries, there has also been work on non-spherical particles and particle fractal aggregates through the use of Monte Carlo techniques (Chan & Dahneke, 1981; Mackowski, 2006; Zhang et al., 2012).

Realistic viral morphology is non-spherical; even approximately spherical capsids have myriad protein protrusions. Indeed, in coronaviruses the surface spike glycoproteins give the viruses their distinctive, eponymous shape. To our knowledge there has not been a thorough investigation of the drag force these spike proteins induce on viral motion in fluids.

Since viral particle diameters span a range of sizes compared to the mean free path of air molecules (~0.06 μm), it is useful to have available results that can span the entire size range (or alternatively, Knudsen number). We have explored a range of computations in this regard in both in the free molecular and the continuum limits and explored approximations over the entire Knudsen number range. For computational efficiency, we have initially approximated spikes on a virus surface by small spheres, however, we will be exploring spheroidal, conical, and other shapes in the future.

The structure of the paper is as follows. In section II, we describe the free molecular results where we have used a Monte Carlo approach. In section III, we describe the continuum results, where we have used a numerical technique to solve a boundary integral equation associated with the the 3D Stokes equation. In section IV, we discuss the Sherman and the ASA approximations of Dahneke that one may use to obtain drag for all Knudsen numbers from the free molecular and the continuum results. In section V, we discuss the gas-surface interaction models and then briefly discuss experimental techniques in section VI. Section VII is a summary of the paper.

Drag in the free molecular limit

Following earlier works, we consider a particle situated in an infinite expanse of a gas flowing with a macroscopic velocity (we will choose its direction as the x-direction). Then the molecular distribution far away from the particle is described by a Maxwellian (Ivchenko et al., 2007):

with nondimensional velocity scaling factors:

In which is the number density, is the mass of a gas molecule, is the Boltzmann constant and is the temperature. The vectors and are the molecular- and particle-velocities, respectively, and the starred variants represent nondimensionalizations of those quantities. Also,

is an absolute Maxwellian.

We will also assume that any molecules reflected from a surface follow Maxwell's diffuse-specular reflection model:Where is a Maxwellian (with a density normalized by conservation of mass), and is the reflection operator. Here is the normal unit vector to the surface pointing into the gas, and is the coefficient of diffuse specular reflection ( is one for purely diffuse reflection and zero for purely specular reflection).

We assume that each gas molecule infinitesimally contributes to the net drag, via the relation:

In the free molecular regime intermolecular collisions are exceedingly rare, and collisions with convex surfaces are non-reentrant. In this case, one may compute the incipient distribution directly, through eq. (1). Particularly, for a spherical particle of radius one obtains:

Note that in eq. (6), the scalar pressure is:

The mobility then is,

and the diffusion coefficient is,

Note that these results depend on both the molecular density and the molecular mass, among other parameters. We will express the nondimensionalized drag as:

Using eq. (6), the drag (for a sphere of radius) in eq. (10) becomes:

And has the values 13.16 for purely diffuse scattering and 9.45 for scattering that is purely specular. We note that Williams gave an alternate derivation of this result for the case of diffuse reflection, via a calculation of the collisional cross-sections. Williams also considered simultaneous heat transfer and additional gas-surface interaction models (Williams, 1973).

For particles much smaller than the molecular mean free path, intermolecular collisions are effectively disregarded, but for non-reentrant surfaces (meaning that molecules reflected from surface point to surface point; either directly, through concave geometry, or indirectly via sky-shine physics), especially ones with non-convex geometries (as viral surfaces can be), it is not possible to obtain helpful analytical expressions.

The reentrant case involves difficult integral equations that account for reflections, and most authors have either resorted to numerical techniques in the face of intractability or used numerical methods such as Monte Carlo algorithms for the flexibility they offer in exploring a wide range of shapes and surface distributions such as found on viral surfaces.

Although both the Test Particle Monte Carlo (TPMC) and the Direct Simulation Monte Carlo (DSMC) techniques can be used for calculating the drag in the Free Molecular regime, in this work we utilize the TPMC technique as it gives us greater flexibility (Bird, 1994, Chan and Dahneke, 1981, Mackowski, 2006;).

Geometric cells as generally used in DSMC, are not needed as intermolecular collisions are neglected Chan & Dahneke, 1981. We also note that both Chan and Dahneke and Mackowski have previously used the TPMC technique to calculate drag on straight chain of uniform particles and agglomerates, respectively. We have found the procedure of Mackowski particularly convenient.

We have constructed a Mathematica program using this procedure to explore the dependence of drag on a spherically shaped virus, with varying distributions of smaller spheres (to approximate spike-glycoproteins) on its surface. In constructing this program, we were able to repurpose Wolfram Language software we constructed previously to model condensational growth on aerosol chains and agglomerates (White et al., 2020; White & Loyalka, 2015).

Since Mackowski described the procedure in detail, we only note essential points. Following Mackowski, we enclose the subject particle(s) in an artificial rectangular (or square) box, and distinguish between surface areas that are parallel, , or perpendicular, , to the direction of motion, and sample the starting molecules from these surfaces in proportion to their particle fluxes.

We sample the molecular speeds (note that Mackowski instead used a mean speed approximation) to support generality and flexibility in computing different gas/surface molecular scattering kernels. We sampled the speed, , by simply using a root finding algorithm on an equation of the form:Where is a uniform random variate on the unit interval, and . This is (theoretically) equivalent to sampling from a particular scaling of the Probability Distribution Function (PDF) of a gamma distribution, or more specifically, the PDF of Nakagami's one-sided distribution, which is itself a special case of the gamma distribution. For an arbitrary mean and standard deviation this distribution is:

We experimented with writing our own sampling routines, as well as using the provided Wolfram Language random variate routines. Our experience was that the provided routines offered superior speed compared to our own, and for our test cases the difference between results obtained from using eqs. (12) or (13)was not significant (see also the Appendix).

We verified the program against known results for a single sphere first and have obtained results in agreement with those from eq. (10), for five different values of . Next, we considered a number, , of equal sized spheres, in a straight chain. We report these results in Table 1 .

Table 1

Nondimensionalized Drag for spheres in a straight chain

We report the nondimensionalized drag, defined via eq. (10), for a the case of a linear chain of spheres (i.e., a group of 2-spheres embedded in a 3 dimensional Euclidean space), similar to a straight line of croquet or billiard balls touching one another, but not resting on any surface. The results are as expected; note that the drag increases slightly with the accommodation coefficient for a given geometry.

NN	α
	0	¼	½	¾	1
1	9.48 ± 0.02	10.37 ± 0.02	11.32 ± 0.02	12.25 ± 0.02	13.15 ± 0.02
2	14.71 ± 0.04	16.32 ± 0.05	17.82 ± 0.05	19.51 ± 0.05	21.08 ± 0.05
3	19.83 ± 0.07	22.23 ± 0.08	24.61 ± 0.08	26.81 ± 0.08	29.23 ± 0.09
4	25.07 ± 0.12	27.98 ± 0.12	31.09 ± 0.13	33.90 ± 0.13	36.75 ± 0.14

Subtracting the single sphere values from the two sphere values we deduced results for the Basic Chain Unit (BCU). We report these results in Table 2 along with a benchmark comparison from Chan and Dahneke.

Table 2

Nondimensionalized drag on a BCU

We tabulate the value of a BCU for five values of the accommodation coefficient and compare our values to a known result. For comparison in the normalizations used in this paper, we multiply the results of Chan and Dahneke by a constant factor. The differences are within about 4%.

α	Test Particle MC	Chan and Dahneke
0	5.23 ± 0.04	5.18 ± 0.04
¼	5.95 ± 0.05	–
½	6.50 ± 0.05	–
¾	7.26 ± 0.05	–
1	7.93 ± 0.05	7.64 ± 0.05

In terms of the BCU, the drag for spheres in a straight chain is approximately:

And the results of Table 1 are in accord with this.

In this work we approximated coronavirus bioparticles by modeling spike glycoproteins as spherical protrusions from a primary body meant to approximate the viral capsid. We first took each of these protrusions to be 1/16th the radius of the primary sphere and considered sixteen different distributions of protrusions as shown in Fig. 1 .

Fig. 1

Coronavirus-mimetic geometries.

The ratio between the radius of the main body and the radius of the protrusions was selected based on some preliminary findings of a coronavirus structure that relate the length of peplomers to capsid body (Schoeman & Fielding, 2019). The choice of these configurations was incremental and meant to extend complexity from a bare central single sphere (in geometry one) to the fully covered, “bumpy” sphere (in geometry sixteen).

The geometries that explored here are the result of trial and error in creating an artistic rendition, of the novel coronavirus. The result of that process was a stochastic geometry that randomizes itself during each rendering; we show one rendition in Fig. 3. However, even if the randomized geometry resembles illustrations of coronaviruses that have appeared online and in popular media, this induces too much complexity for a first study.

Fig. 3

Randomized Coronavirus Geometry.

Shown is a coronavirus mimetic geometry with usage of stochastically placed cones tipped with spheres as computational phantoms for the spike glycoproteins, and small, hemispherical protrusions randomly placed to model the membrane glycoproteins embedded in the capsid surface. Unlike the geometries represented in this study, the location, unlike the geometries represented in this study, we placed the surface features non-deterministically via sampling random matrices, in a procedure we intend to further explore in a later paper.

Rather than performing hundreds or thousands of simulations and reporting statistical distributions that are themselves statistical results on Monte Carlo simulations, we opted to simulate a bare sphere, then iteratively add protrusions to incrementally examine the effect of adding these protrusions on the drag experienced by the particle system.

We report our corresponding TPMC results for drag in Table 3 . We additionally computed the nondimensionalized drag obtained by the Volume Equivalent Sphere (VES) approximation. The radius of the volume equivalent sphere is

Table 3

Monte Carlo normalized drag simulation results.

Results of the drag Monte Carlo simulation across five different coefficients of diffuse/specular reflection, for sixteen different geometries meant to mimic the novel coronavirus SARS-CoV-2. Please note that the simulation index, , is arbitrary, and the sphere number, , is the number of protrusions from the main body. Within each geometry set, the top lines are the results of our Monte Carlo Simulations, the second lines are the drag experienced by a volume equivalent sphere, and the third lines are the ratios of the first two, respectively. The drag for the volume equivalent sphere cases was computed with the aid of eqs. (5) - (11). The results show that for each simulation geometry, an increase in the quantity of the capsid surface corresponds to an increase in the drag experienced by the particle. The uncertainty in the drag, which as the result of the Monte Carlo simulation, is the less precise than values obtained via evaluation of analytic formulae, drives the uncertainty of the calculations in general.

N(a)	Coefficient of Diffuse/Specular Reflection (α)
(S)(b)	0	¼	½	¾	1
1 (0)	9.483 ± 0.018	10.374 ± 0.019	11.318 ± 0.019	12.249 ± 0.020	13.150 ± 0.020
	9.453 ± 0.001	10.381 ± 0.001	11.309 ± 0.001	12.238 ± 0.001	13.165 ± 0.001
	1.003 ± 0.002	0.999 ± 0.002	1.001 ± 0.002	1.001 ± 0.002	0.999 ± 0.002
2 (1)	9.432 ± 0.019	10.379 ± 0.019	11.318 ± 0.020	12.258 ± 0.020	13.160 ± 0.021
	9.455 ± 0.001	10.383 ± 0.001	11.311 ± 0.001	12.239 ± 0.001	13.167 ± 0.001
	0.998 ± 0.002	1.000 ± 0.002	1.001 ± 0.002	1.002 ± 0.002	0.999 ± 0.002
3 (1)	9.432 ± 0.018	10.353 ± 0.019	11.299 ± 0.020	12.228 ± 0.020	13.137 ± 0.021
	9.455 ± 0.001	10.383 ± 0.001	11.311 ± 0.001	12.239 ± 0.001	13.167 ± 0.001
	0.998 ± 0.002	0.997 ± 0.002	0.999 ± 0.002	0.991 ± 0.002	0.998 ± 0.002
4 (4)	9.472 ± 0.018	10.449 ± 0.019	11.312 ± 0.020	12.285 ± 0.021	13.156 ± 0.021
	9.459 ± 0.001	10.388 ± 0.001	11.317 ± 0.001	12.245 ± 0.001	13.174 ± 0.001
	1.001 ± 0.002	1.006 ± 0.002	1.000 ± 0.002	1.003 ± 0.002	0.999 ± 0.002
5 (8)	9.602 ± 0.019	10.526 ± 0.019	11.448 ± 0.002	12.371 ± 0.020	13.275 ± 0.021
	9.465 ± 0.001	10.395 ± 0.001	11.324 ± 0.001	12.253 ± 0.001	13.182 ± 0.001
	1.014 ± 0.002	1.013 ± 0.002	1.011 ± 0.000	1.010 ± 0.002	1.007 ± 0.002
6 (16)	9.720 ± 0.019	11.597 ± 0.020	11.531 ± 0.020	12.454 ± 0.021	13.364 ± 0.021
	9.478 ± 0.001	10.408 ± 0.001	11.339 ± 0.001	12.269 ± 0.001	13.200 ± 0.001
	1.026 ± 0.002	1.018 ± 0.002	1.017 ± 0.002	1.015 ± 0.002	1.0124 ± 0.002
7 (32)	9.771 ± 0.019	10.697 ± 0.020	11.609 ± 0.020	12.531 ± 0.020	13.433 ± 0.021
	9.502 ± 0.001	10.435 ± 0.001	11.368 ± 0.001	12.301 ± 0.001	13.234 ± 0.001
	1.028 ± 0.002	1.025 ± 0.002	1.021 ± 0.002	1.019 ± 0.002	1.015 ± 0.002
8 (6)	9.553 ± 0.019	10.472 ± 0.019	11.427 ± 0.019	12.378 ± 0.020	13.281 ± 0.021
	9.462 ± 0.001	10.391 ± 0.001	11.320 ± 0.001	12.249 ± 0.001	13.178 ± 0.001
	1.010 ± 0.002	1.008 ± 0.002	1.009 ± 0.002	1.011 ± 0.002	1.008 ± 0.002
9 (14)	9.743 ± 0.019	10.625 ± 0.020	11.562 ± 0.020	12.495 ± 0.020	13.422 ± 0.020
	9.628 ± 0.001	10.573 ± 0.001	11.518 ± 0.001	12.463 ± 0.001	13.409 ± 0.001
	1.012 ± 0.002	1.005 ± 0.002	1.004 ± 0.002	1.004 ± 0.002	1.001 ± 0.002
10 (30)	9.956 ± 0.019	10.885 ± 0.020	11.784 ± 0.020	12.713 ± 0.020	13.597 ± 0.021
	9.499 ± 0.001	10.432 ± 0.001	11.364 ± 0.001	12.297 ± 0.001	13.230 ± 0.001
	1.048 ± 0.002	1.044 ± 0.002	1.037 ± 0.002	1.034 ± 0.002	1.028 ± 0.002
11 (62)	10.129 ± 0.019	11.030 ± 0.020	11.951 ± 0.020	12.805 ± 0.021	13.709 ± 0.021
	9.548 ± 0.001	10.486 ± 0.001	11.423 ± 0.001	12.360 ± 0.001	13.298 ± 0.001
	1.061 ± 0.002	1.052 ± 0.002	1.046 ± 0.002	1.036 ± 0.002	1.031 ± 0.002
12 (8)	9.607 ± 0.019	10.546 ± 0.019	11.487 ± 0.020	12.397 ± 0.020	13.297 ± 0.021
	9.465 ± 0.001	10.395 ± 0.001	11.324 ± 0.001	12.253 ± 0.001	13.182 ± 0.001
	1.015 ± 0.002	1.015 ± 0.02	1.015 ± 0.002	1.012 ± 0.002	1.009 ± 0.002
13 (20)	9.874 ± 0.019	10.797 ± 0.020	11.699 ± 0.020	12.651 ± 0.021	13.533 ± 0.021
	9.483 ± 0.001	10.415 ± 0.001	11.346 ± 0.001	12.277 ± 0.001	13.208 ± 0.001
	1.041 ± 0.02	1.037 ± 0.002	1.031 ± 0.002	1.031 ± 0.002	1.025 ± 0.002
14 (44)	10.233 ± 0.019	11.186 ± 0.020	12.075 ± 0.020	12.944 ± 0.020	13.843 ± 0.021
	9.521 ± 0.001	10.455 ± 0.001	11.390 ± 0.001	12.325 ± 0.001	13.259 ± 0.001
	1.075 ± 0.002	1.070 ± 0.002	1.060 ± 0.002	1.060 ± 0.002	1.044 ± 0.002
15 (92)	10.465 ± 0.020	11.368 ± 0.020	12.273 ± 0.020	13.137 ± 0.020	13.999 ± 0.021
	9.594 ± 0.001	10.536 ± 0.001	11.478 ± 0.001	12.420 ± 0.001	13.362 ± 0.001
	1.091 ± 0.002	1.079 ± 0.002	1.069 ± 0.002	1.058 ± 0.002	1.048 ± 0.002
16 (114)	11.569 ± 0.020	12.441 ± 0.020	13.287 ± 0.021	14.075 ± 0.021	14.887 ± 0.021
	9.628 ± 0.001	10.573 ± 0.001	11.518 ± 0.001	12.463 ± 0.001	13.409 ± 0.001
	1.202 ± 0.002	1.177 ± 0.002	1.154 ± 0.002	1.129 ± 0.002	1.110 ± 0.002

and the corresponding drag,Where is the nondimensionalized drag as given by eq. (10), with a radius . Note that the volume equivalent drag is then:

and these results are straightforward to compute from the data in Table 3. We see that the VES approximation is quite good (we note however that for the straight chain of uniform spheres, the VES approximation is better when the number of spheres is small (say, four or less) and worsens when n is larger).

To explore the usefulness of the volume equivalent sphere approximation further, we have obtained added results by scaling the radii of the geometry by a factor of ¾ and performing the same computations for the geometric layouts defined by cases two through sixteen in Fig. 1. We report these results in Table 4 .

Table 4

Monte Carlo normalized drag simulation results for Shrunken geometry results

of the drag Monte Carlo simulation across five different coefficients of diffuse/specular reflection, for sixteen different geometries meant to mimic the novel coronavirus SARS-CoV-2. The geometries used here are identical to the geometries used to create Table 3, the primary difference being that we scaled down the spheres used to model the particles by a factor of ¾, to ease some comparison on the relationship between corona size and drag. Please note that we did not duplicate experiment one replicated, as it consisted of a single bare sphere. Otherwise, the descriptions from Table 3 are still applicable, as well as the pictures in panel graphic Fig. 1.

N	Coefficient of Diffuse/Specular Reflection (α)
(S)	0	¼	½	¾	1
2 (1)	9.444 ± 0.009	10.368 ± 0.010	11.294 ± 0.010	12.225 ± 0.010	13.171 ± 0.010
	9.454 ± 0.001	10.382 ± 0.001	11.310 ± 0.001	12.238 ± 0.001	13.166 ± 0.001
	0.999 ± 0.001	0.999 ± 0.001	0.999 ± 0.001	0.999 ± 0.001	1.000 ± 0.001
3 (1)	9.435 ± 0.009	10.377 ± 0.009	11.289 ± 0.010	12.228 ± 0.010	13.164 ± 0.010
	9.456 ± 0.001	10.384 ± 0.001	11.312 ± 0.001	12.241 ± 0.001	13.169 ± 0.001
	0.999 ± 0.001	0.999 ± 0.001	0.998 ± 0.001	0.999 ± 0.001	1.000 ± 0.001
4 (4)	9.481 ± 0.009	10.404 ± 0.010	11.327 ± 0.010	12.249 ± 0.010	13.205 ± 0.010
	9.458 ± 0.001	10.384 ± 0.001	11.312 ± 0.001	12.241 ± 0.001	13.169 ± 0.001
	1.003 ± 0.001	1.002 ± 0.001	1.001 ± 0.001	1.001 ± 0.001	1.003 ± 0.001
5 (8)	9.524 ± 0.009	10.434 ± 0.010	11.381 ± 0.010	12.301 ± 0.010	13.220 ± 0.010
	9.458 ± 0.001	10.387 ± 0.001	11.315 ± 0.001	12.244 ± 0.001	13.173 ± 0.001
	1.007 ± 0.001	1.005 ± 0.001	1.006 ± 0.001	1.005 ± 0.001	1.004 ± 0.001
6 (16)	9.597 ± 0.009	10.513 ± 0.010	11.440 ± 0.010	12.355 ± 0.010	13.272 ± 0.010
	9.463 ± 0.001	10.393 ± 0.001	11.322 ± 0.001	12.251 ± 0.001	13.180 ± 0.001
	1.013 ± 0.001	1.012 ± 0.001	1.011 ± 0.001	1.009 ± 0.001	1.007 ± 0.001
7 (32)	9.669 ± 0.009	10.591 ± 0.010	11.494 ± 0.010	12.410 ± 0.010	13.338 ± 0.010
	9.474 ± 0.001	10.404 ± 0.001	11.334 ± 0.001	12.264 ± 0.001	13.194 ± 0.001
	1.021 ± 0.001	1.018 ± 0.001	1.014 ± 0.001	1.012 ± 0.001	1.011 ± 0.001
8 (6)	9.509 ± 0.009	10.443 ± 0.010	11.360 ± 0.010	12.285 ± 0.010	13.233 ± 0.010
	9.457 ± 0.001	10.385 ± 0.001	11.314 ± 0.001	12.242 ± 0.001	13.171 ± 0.001
	1.006 ± 0.001	1.006 ± 0.001	1.004 ± 0.001	1.003 ± 0.001	1.005 ± 0.001
9 (14)	9.606 ± 0.009	10.547 ± 0.010	11.461 ± 0.010	12.375 ± 0.010	13.282 ± 0.010
	9.462 ± 0.001	10.391 ± 0.001	11.320 ± 0.001	12.249 ± 0.001	13.178 ± 0.001
	1.015 ± 0.001	1.015 ± 0.001	1.012 ± 0.001	1.010 ± 0.001	1.008 ± 0.001
10 (30)	9.759 ± 0.009	10.662 ± 0.010	11.584 ± 0.010	12.487 ± 0.010	13.408 ± 0.010
	9.473 ± 0.001	10.403 ± 0.001	11.333 ± 0.001	12.262 ± 0.001	13.192 ± 0.001
	1.030 ± 0.001	1.025 ± 0.001	1.022 ± 0.001	1.018 ± 0.001	1.016 ± 0.001
11 (62)	9.909 ± 0.009	10.822 ± 0.010	11.719 ± 0.010	12.599 ± 0.010	13.488 ± 0.010
	9.493 ± 0.001	10.425 ± 0.001	11.357 ± 0.001	12.289 ± 0.001	13.221 ± 0.001
	1.044 ± 0.001	1.038 ± 0.001	1.032 ± 0.001	1.025 ± 0.001	1.020 ± 0.001
12 (8)	9.544 ± 0.009	10.468 ± 0.010	11.391 ± 0.010	12.327 ± 0.010	13.233 ± 0.010
	9.458 ± 0.001	10.387 ± 0.001	11.315 ± 0.001	12.244 ± 0.001	13.173 ± 0.001
	1.009 ± 0.001	1.008 ± 0.001	1.007 ± 0.001	1.007 ± 0.001	1.005 ± 0.001
13 (20)	9.696 ± 0.009	10.613 ± 0.010	11.525 ± 0.010	12.458 ± 0.010	13.366 ± 0.010
	9.466 ± 0.001	10.395 ± 0.001	11.325 ± 0.001	12.254 ± 0.001	13.183 ± 0.001
	1.024 ± 0.001	1.021 ± 0.001	1.018 ± 0.001	1.017 ± 0.001	1.014 ± 0.001
14 (44)	9.908 ± 0.010	10.823 ± 0.010	11.720 ± 0.010	12.620 ± 0.010	13.518 ± 0.010
	9.482 ± 0.001	10.413 ± 0.001	11.343 ± 0.001	12.274 ± 0.001	13.205 ± 0.001
	1.045 ± 0.001	1.040 ± 0.001	1.033 ± 0.001	1.028 ± 0.001	1.024 ± 0.001
15 (92)	10.143 ± 0.010	11.026 ± 0.010	11.906 ± 0.010	12.803 ± 0.010	13.659 ± 0.010
	9.513 ± 0.001	10.447 ± 0.001	11.381 ± 0.001	12.314 ± 0.001	13.248 ± 0.001
	1.066 ± 0.001	1.055 ± 0.001	1.046 ± 0.001	1.040 ± 0.001	1.031 ± 0.001
16 (114)	10.662 ± 0.010	11.549 ± 0.010	12.419 ± 0.010	13.267 ± 0.010	14.107 ± 0.011
	9.527 ± 0.001	10.462 ± 0.001	11.398 ± 0.001	12.333 ± 0.001	13.268 ± 0.001
	1.119 ± 0.001	1.104 ± 0.001	1.090 ± 0.001	1.076 ± 0.001	1.063 ± 0.001

We present an image of the geometries used for the Monte Carlo simulations and various calculations throughout the paper. Note that case number one is not shown as it consisted of a single, isolated sphere. Geometries two and three are identical save for the orientation of their single protrusions. In each of the geometries pictured, under the nondimensionalizations used in this paper, the larger, central body-particle has a radius of one. The exception to this would be the geometries used in the construction of Table 3. Table 3 geometries have main spheres with a radius of ¾. Primarily, and exceptions are noted where appropriate, the small spheres surrounding the main bodies all have radii of 1/16.

We discuss here first the results relating to Fig. 1 and Table 3. We reference each computation by our so-called simulation number, , which is a completely arbitrary index, and the number of small spheres on the larger sphere as the index . However, not all geometries are alike and for visual purposes readers may reference Fig. 1. Exact coordinate chart data showing the specifications for each geometry are available online via the publisher.

We note from Table 3 that the addition of spheres to the surface increases the drag. The drag also increases as diffuse reflection begins to dominate over specular reflection. Comparing simulation five to simulation twelve (both have eight protrusions, with a different spatial distribution), we see trivial difference between the two drag results.

However, extending the number of spheres and the complexity we see differences increase such that disorder and spread in sphere location increases drag, as seen specifically in comparing eleven with sixty-two small spheres and fourteen with forty-four small spheres. Despite having less spheres, simulation fourteen had higher drag and deviates further from the VES approximation.

Finally, simulation sixteen with 144 spheres has significant deviation from VES approximation greater than 10%. This simulation acts as the extreme example that a spread-out layout of smaller spheres and highlights the advantage of performing the computations.

Comparing the geometries with the larger geometries (small-sphere radius of 1/16) and the smaller geometries (small-sphere radius of 3/64) is insightful. We see that in the smaller cases there is much more uniformity with respect to the number and distribution of the small spheres, compared to the larger case where there is a much more noticeable increase in drag force as the number of protrusions increases. We noticed the same increase in drag with the accommodation coefficient, , for the scaled-down simulation.

We have computed drag on a range of virus configurations in the free-molecular limit by constructing and using a versatile Test Particle Monte Carlo (TPMC) program. Our results show that increasing protrusions result in higher drag. Our results also show that diffuse molecular-reflection increases drag, but also has a higher sensitivity to the number (and placement) of protrusions. Finally, we find that the volume equivalent sphere approximation supplies a good estimation of the drag for the types of geometries we studied, under the assumptions we made. The volume equivalent sphere approximation however gives an increasingly worse estimate of the drag as the number of protrusions increases.

We should note that in our program, we used about four million particles and parallel processing (4 cores). The program can be slow for complicated geometry computations (partly due to incorporation of visualization and use of Mathematica). We have not explored its optimization as it is in the development stage. Its use on a hundred or thousand core processor (for example) does not require any change.

Further, we have reported above drag results only along the x-axis (in line with our assumption of gas flow on this axis only), and the geometries that we considered here are (almost) spherically symmetric. However, we also computed the full mobility tensor and computed the averaged drag force via the harmonic mean:Where the are the eigenvalues of the mobility tensor; we have reported typical values in Table 5 .

Table 5

Drag averaged over spatial orientation

The table shows orientations averaged results, presented for four geometries used in the paper for . The results agree well with the directional drag results presented Table 3.

N	A_⊥
(S)
4 (4)	13.290
6 (16)	13.410
12 (8)	13.352
16 (114)	14.835

We further note that the ion mobility, , iswhere the product term is the charge on the ion and is the mobility defined in eq (8). During the initial review of the paper, we learned Larriba-Andaluz and Hogan had computed for ions composed of aggregates of spheres, in an executable form of their program IMoS (Ion Mobility Software), Coots, et al. (2020), which is available for download and public use. (Coots et al., 2020; Larriba-Andaluz & Hogan, 2013a,2013b; Larriba-Andaluz & Hogan, 2013a,2013b; Larriba-Andaluz & Hogan, 2014; Schoeman & Fielding, 2019) They had considered two approaches for calculating (essentially the mobility, , using TPMC (following Mackowski, with provisions for both cuboid and spherical enclosures), and the Mason-Schamp first approximation of the two-temperature theory: (Viehland & Mason, 1978)

Hence,Where, is the gas (molecular) density, is the reduced mass (note that under the assumptions that is the mass of the background gas, and is the mass of the much heavier ionic target), is an effective temperature (assuming the ion and the gas are at different temperatures and an ionic drift speed ), and is a collisional integral (computing it requires assuming both an ion-molecule potential function and scattering as elastic or inelastic). In the simple case of rigid sphere and isothermal interactions, one has: (Chapman & Cowling, 1970; Viehland & Mason, 1978)

Note that eq. (22) is simply the physical cross-section. For the diffusion coefficient this gives:

But this is different from the result of eq. (9):

The two results (23) and (24) agree only when the accomodation coefficient vanishes, i.e., the case of specular reflection. The differences are likely due to different assumptions regarding molecule-molecule (ion) interactions (diffuse scattering vs elastic collisions), and have led to several options for elastic, inelastic scattering, etc. in IMoS). (Larriba-Andaluz & Hogan, 2013a,2013b) The issue is not the first order Chapman-Enskog approximation in this case (it does not differ much from the exact result for the case of rigid spheres. (Tompson & Loyalka, 1987) (A Mathematica program for evaluation of collision integrals for elastic scattering for Lennard Jones potentials is given in the book by Ivchenko et al. (2007) and can be changed to suit other potentials and/or inelastic scattering).

Note that, in terms of the collisional integral, one can express the nondimensionalized drag force as:

We report results from the IMoS program for the non-dimensional drag in Table 6 ; the results are comparable with the results we obtained from our program.

Table 6

IMoS Computed Drag on Biomimetic Particles

The table shows results from the IMoS program for the nondimensionalized drag for the biomimetic geometries presented in Fig. 1. The results show that for each simulation geometry an increase in the total surface area of the particle corresponds to an increase in drag, defined in eq. (25), experienced by the particle.

N	α
	0	¼	½	¾	1
1	9.48	10.41	11.34	12.27	13.21
2	9.52	10.43	11.31	12.18	13.03
3	9.50	10.43	11.36	12.29	13.22
4	9.53	10.46	11.39	12.31	13.24
5	9.59	10.52	11.44	12.35	13.28
6	9.70	10.62	11.53	12.44	13.35
7	9.89	10.80	11.69	12.58	13.47
8	9.56	10.49	11.41	12.33	13.26
9	9.68	10.60	11.51	12.42	13.34
10	9.90	10.81	11.70	12.59	13.48
11	10.28	11.16	12.02	12.87	13.72
12	9.59	10.52	11.44	12.35	13.28
13	9.77	10.68	11.59	12.49	13.40
14	10.09	10.99	11.87	12.74	13.61
15	10.66	11.52	12.35	13.16	13.96
16	11.02	11.85	12.66	13.43	14.19

Drag in the continuum limit

Applications of both TPMC and the DSMC to computations for arbitrary Knudsen number would be most useful, however, the Monte Carlo methods have limitations for infinite domain problems because of decay (of moments such as velocity) at large distances and the convergence for Knudsen numbers in the continuum can be extremely slow or deceptive. If one can truncate the problem geometry without affecting the accuracy of the simulation, and if particle scattering events are isotropic, then one can use the so-called Walk on Spheres (WoS) method(s) to perform a particularly efficient (grid free) form of Monte Carlo by recasting integral equations as stochastic processes. They are particularly suitable for graphics problems, e.g., see the work of Sawhney and Crane (2020) for a recent review. Due to the restrictions of WoS methods, stemming from the subtlety of the connection between the Itô calculus and PDEs, especially their inability to handle infinite domains, and Neumann- or Robin-type boundary conditions, an alternate method was desired.

For more background information about the techniques discussed here, we recommend (in addition to other works cited in this section) the reference texts by Wyld and Powell (2020), Leal (2007), Pozrikidis (1992), respectively.

In the continuum regime, the drag and torque can be computed through solutions of the Stokes Equation, and for a rigid body(s) of arbitrary shape by solving the boundary integral equation corresponding to it:Where the Green’s function (tensor) used in eq. (26) is as defined in eq. (27). In these equations u is the fluid is the fluid translational velocity, and is the angular velocity, respectively, far from the fixed rigid bodies. Details of the derivation of eq. (26) from the Stokes equation are well explained by Leal (we refer the reader to chapter eight of that text) (Leal, 2007). For the computations of interest here we will assume , i.e. that the rotational component vanishes. Note that in the literature, authors may label eq. (27) as the Oseen tensor (especially when it is derived from Oseen's equations; fundamental solutions may be further called Oseenlets) or the Stokes tensor (with fundamental solutions sometimes called Stokeslets). To be specific, the difference between the two is that technically within an Oseen analysis there is consideration of convective (pure fluid) or advective (fluid as a carrier) acceleration, while in a Stokes analysis these effects are unconsidered. Here, is the stress tensor, is the vector of stress components we need to compute the drag, is a point on the body, and is the unit vector normal to the surface at directed from the fluid (gas) to the body. We obtain the drag then by integrating over the surface(s) of the body. We use the subscript ‘cont.’ to show that this is the continuum result, and to distinguish it from the free molecular result, which we will show with a ‘fm.’ subscript.

Eq. (26) is a system of three singular integral equations; Youngren and Acrivos (1975) solved it numerically, by using a boundary element method for certain spheroids. We have found it far simpler to solve it by using a singularity subtraction technique (for details of the technique we refer the reader to previous work on condensation rate on aggregates of spheres, spheroids and ellipsoids in the continuum regime, and on the spinning rotor gauge in the slip regime). (Loyalka, 1996; Loyalka & Griffin, 1993; Smith & Loyalka, 2014).

We write eq. (26) as:Where we have defined the integrated Green's tensor:

We evaluate the tensor from eq. (30) numerically, and use Gaussian quadrature for the second integral, noting that integrand vanishes for , and remains unevaluated at that point. For Gaussian quadrature of order , we use points corresponding to the polar and azimuthal integrations over the intervals and , respectively. We use a collocation method (over all surfaces in the geometry) on eq. (29) to solve the resulting set of rank three equations (there are three equations corresponding to each surface point, one for each spatial direction, and the total number of equations is thus , where is the number of bodies. We have not assumed any spatial symmetry, and there is considerable flexibility in choosing the integration and quadrature schemes; for example, one could use separate quadrature meshes in the polar and azimuthal integrations).

In the cases (spheres, spheroids including needles, discs, linear chains of spheres), that we used for verification, we found very satisfactory results even with low order quadrature. Additionally, we have successfully used the program for ellipsoids also, but are still testing and improving it. We will be reporting separately the full details of the technique with analytical benchmarking for both translation and rotation of arbitrary bodies (see the monograph by Williams and Loyalka (1991) for a detailed description of cases and results).

Note that for a single sphere of radius , one has the well-known Stokes result:

And for a linear chain of two touching spheres each of radius , one has the exact result: (Williams and Loyalka, 1991)

We report numerical results that we obtained for the case of a linear chain of equal sized spheres (assuming fluid flow , i.e. flow solely along the x-axis) in Table 7 . Our numerical results agree with (31), (32) and (32), almost precisely, even with a low order quadrature. We have shown added results also in Table 8 .

Table 7

Nondimensionalized Drag on a Linear Chain in the Continuum Regime

We report the nondimensionalized drag, benchmarks for the Stokes (continuum) problem for the geometric case of linear chains of touching spheres. In the table, is the number of spheres in that linear chain, and corresponds to quadrature points taken per sphere. To ease comparison with other results presented in this paper, we took the radius of the spheres to be one (under our nondimensionalization). For two spheres, the exact result from eq. (32) is 24.3212. For three and four spheres, the results are in good agreement with those of Geller, who used a boundary element technique (they obtained in the present normalization, values of 24.74, 28.98 and 33.24 for two-, three- and four-sphere chains respectively (Gelleret al., 1993).

S	2^Q
	1	2	3	4
1	18.244	18.849	18.850	18.850
2	23.466	24.325	24.337	24.323
3	27.582	29.016	29.076	29.057
4	31.424	33.310	33.409	33.387

Table 8

Nondimensionalized Drag on Selected Virus Mimetic Geometries in the Continuum Regime

We report the nondimensionalized drag, , for a sample of the biomimetic geometries considered in this study. We selected geometries for inclusion in this table if the number of bodies protruding from the main bulk was eight or less. Note that it appears in the continuum regime, increasing the bumpiness of the particle system increases the drag, but not dramatically so. The symbols , and keep their meanings as the index of the geometry and the number of quadrature points, respectively.

G	2^Q
	1	2	3
1	18.245	18.849	18.849
2	18.245	18.849	18.849
3	18.245	18.849	18.849
4	18.295	18.890	18.886
5	18.380	18.959	18.940
8	18.326	18.923	18.881
12	18.368	18.921	18.955

We have reported drag results, for selected configurations of the virus mimetic particle of Table 3 and corresponding to in Table 8. These results are helpful in tandem with those for the respective free molecular cases, in the approximations discussed in the next section. We have chosen here configurations with fewer spheres as the computations are time consuming, despite parallelization. We plan to report more extensive results with spikes shaped as thin needles, in future papers.

We have also calculated the drag tensor through variations of velocity , and obtained the orientation averaged drag, but are not reporting those results here. We plan to discuss these, together with alternate methods of calculating (approximating) the average, separately.

Like the results in the last section, one can examine a VES for computations in the continuum range. We have shown the results of this approximation for all sixteen geometries in Table 9 . Note that for the geometries 1–5, 8 and 12 the results of Table 9 are very close to those of Table 6, and it seems that the VES approximation is quite good for all the sixteen geometries considered here. Thus, further numerical computations for other geometries via the numerical technique do not appear necessary.

Table 9

VES Drag for Virus Mimetic Particles in the Continuum Regime

We report the results for the calculation of the nondimensionalized drag for the VES approximation for the geometries shown in Fig. 1, corresponding to all sixteen geometries used in Table 3. This table is the product of an indirect computation; first, we computed the drag on a bare sphere of varying radius, in the vicinity of one (to get good resolution around ), the result of which we interpolated to generate the results seen here. Note that in the table, is the geometry index, is the number of body protrusions, is the total volume of the geometry, is the VES radius, and is the nondimensionalized drag. We note that the original computations used eight quadrature points per sphere. We completed the same calculations with a of sixteen and found the results identical to at least four decimal places. Admittedly, we could have obtained the results for the bare sphere from eq. (31) without performing a Stokes calculation, but this was done for benchmarking purposes.

G	N	V	R_VES	F
1	0	4.1888	1.00000	18.850
2	1	4.1898	1.00008	18.851
3	1	4.1898	1.00008	18.851
4	4	4.1929	1.00033	18.855
5	8	4.1970	1.00065	18.861
6	16	4.2052	1.00130	18.874
7	32	4.2215	1.00260	18.898
8	6	4.1943	1.00049	18.858
9	14	4.2031	1.00114	18.871
10	30	4.2195	1.00244	18.895
11	62	4.2522	1.00502	18.944
12	8	4.1970	1.00065	18.861
13	20	4.2092	1.00162	18.880
14	44	4.2338	1.00357	18.916
15	92	4.2829	1.00743	18.989
16	114	4.3054	1.00919	19.022

Approximate expressions for drag over the entire Knudsen number range

Computations over the entire Knudsen number range for particles of arbitrary shape are challenging. Monte Carlo simulations as well as other direct numerical methods have extremely poor convergence in the continuum limit. Therefore, it is practical to consider an approximate expression for the drag for all Knudsen numbers as derived from a combination of the free molecular and the continuum results. One such expression is based on the harmonic approximation, known as the Sherman's approximation in Rarefied Gas Dynamics (Buckley & Loyalka, 1989; Cercignani et al., 1968; Lea & Loyalka, 1982; Loyalka, 1992):

Note that appx. (33) appears similar to a harmonic mean, but is in fact different, as it is meant to be an interpolation, not a statistic computed over a set of values. The subscripts ‘fm.’ and ‘cont.’ once again denote the free molecular and the continuum values of the drag, respectively. Defining the so-called Cunningham slip (or correction) factor as the ratio:

We can express the drag force (34) as:

Note that using the free molecular and the continuum drag expression for a single sphere, as in the Sherman approximation, we have (we have displayed below the dependence of the Cunningham Slip Factor):

The Knudsen number is once again:

and the mean free path is

We note that different authors specify distinct definitions for the mean free path, differing by multiplicative factors (such as ), and one needs to take care to use consistent definitions when comparing the work of different authors.

It has been found that the Sherman approximation (eqs. (33) - (36), with equal to one) results for the drag and Cunningham's Slip Factor may differ from (numeric) results obtained from the BGK model or the Boltzmann Equation (again with an accommodation coefficient, , of one), and Millikan's data by at most ten percent in the transition regime (Buckley & Loyalka, 1989; Cercignani et al., 1968; Lea & Loyalka, 1982; Loyalka, 1992).

Since we have reported both the free molecular and the continuum results in Tables 3 and 8, respectively, it is straightforward to use eq. (33) for the entire size range by using the appropriate values for a particular configuration of the particle for which the values have been given here; i.e., with the free-molecular drag:

and the continuum drag:Where, and (see eq. (28)) are the corresponding results for a configuration as reported in Tables 3 and 8, respectively, one simply has:

and,

For any other configuration, the programs created here provide the needed free molecular and the continuum drag values, and eq. (42) may be used to interpolate drag values on the full range of Knudsen numbers. Since for the geometries considered here the VES approximation has given good results for both the free molecular and continuum limits (the worst case was a difference of 10% in the free molecular limit), a reasonable estimate of can be obtained by using in place of , and corresponding single bare sphere results in eq. (48).

We further note that Dahneke's Adjusted Sphere Approximation (ASA) has found considerable use in the literature (Dahneke, 1973a, 1973b, 1973c). In this approximation, one calculates an adjusted sphere Knudsen number by solving the equation:

The two Knudsen numbers in eq. (43) necessarily differ; is a Knudsen number characteristic of the body, and is the Knudsen number of the adjusted sphere. We typify Millikan's results:

in describing the drag on a sphere, where the symbol describes the Cunningham slip (or correction) factor evaluated for the adjusted sphere, defined as:

And the coefficients have the approximate values of , and . Thus, one simply has:That is, will depend on all the factors (shape, orientation) describing the body. We can now estimate the drag on the body using eq. (44), to be described below.

The ASA approximation is easy to implement, and we have found that both the Sherman and the ASA approximation yield close results. We show results for the Cunningham slip factor in Sherman's approximation and the ASA approximation, in Fig. 2 , for a needle like spheroid (we take in Dahneke's nomenclature the parameters of the needle as , and , and the accommodation coefficient as ).

Fig. 2

Approximations of Sherman and Dahneke for Cunningham Slip Factor

We present a comparison of Sherman's approximation and Dahneke's ASA for the Cunningham Slip Factor for a specific spheroid. Note that for vanishing Knudsen number, (the continuum limit) and for diverging Knudsen number, (the free molecular limit) the difference between the two approximations vanishes. It is only in the transition regime that there is a small, but noticeable difference.

Zhang et al. (2012) further verified that the ASA approximation gives results that are comparable to the DSMC results for particle aggregates over a range of transition-regime Knudsen numbers (Chan & Dahneke, 1981).

We believe both the Sherman and the ASA approximations are especially useful, but their foundation is heuristic, not theoretical. Neither of them yields correct slip coefficients when examined in the near continuum regime, and their derivatives in the near free molecular regimes do not have logarithmic terms as they should. However, the error induced by these limitations is typically small despite them. Efforts to improve upon these approximations would remain of interest (see Buckley and Loyalka (1989) for an exploration).

We have further explored the effect of varying the radius of the small protrusions, while keeping a constant radius of the larger particle main body (we use this radius as the basic length unit, so varying it at the same time is unnecessary). The results of this comparison are shown in Table 10 for free-molecular computations and Table 11 for continuum computations, respectively. The results show that for small protrusions from a large body, the drag is insensitive to these perturbations, but the effect becomes larger for geometries with more protrusions. We will explore these effects) more in a later paper where we will consider non-spherical protrusions (needle like, among others) of varying sizes and examine the usefulness of the VES approximation further.

Table 10

Protrusions of varying size in the free molecular regime

We examine the effect of varying the size of the protrusions (shown in the column headings), while keeping the main-body radius constant for three different geometries in the free-molecular regime. Within each case there are three lines, the first line shows the computed nondimensionalized drag, the second line shows the VES drag for the system, and the third line is the ratio of the first two. The results show a slight increase in the error of the VES approximation with the accommodation coefficient, and a variance with respect to the size and number of protrusions that is not simple to quantify.

α	G	R/R_2
		1/20	1/16	5/64
0	2, 3	9.449	9.432	9.426
		9.454	9.455	9.456
		0.999	0.998	0.997
	4	9.488	9.472	9.519
		9.456	9.459	9.465
		1.003	1.001	1.006
	8	9.499	9.553	9.617
		9.457	9.462	9.471
		1.004	1.010	1.015
½	2, 3	11.303	11.318	11.295
		11.310	11.311	11.313
		0.999	1.001	0.998
	4	11.337	11.312	11.379
		11.312	11.317	11.324
		1.002	1.000	1.005
	8	11.367	11.427	11.482
		11.314	11.320	11.331
		1.005	1.009	1.013
1	2, 3	13.165	13.160	13.174
		13.166	13.167	13.170
		1.000	0.999	1.000
	4	13.169	13.156	13.219
		13.169	13.174	13.182
		1.000	0.999	1.003
	8	13.205	13.281	13.322
		13.171	13.178	13.190
		1.003	1.008	1.010

Table 11

Protrusions of varying size in the continuum regime

We examine the effect of varying the size of the protrusions (shown in the column headings), while keeping a static main-body radius of one for three different geometries in the continuum regime. For each case there are two lines, the first line is the drag computed by the GF method, while the second line is the ratio of the first line and the VES drag. The results show the VES approximation is good for the cases considered, but the quality of the approximation decreases as the size and number of the protrusions from the main body increases.

G	1/20	1/16	5/64
2, 3	18.849	18.849	18.850
	0.997	0.997	0.997
4	18.866	18.886	18.917
	0.997	0.999	1.001
8	18.864	18.881	19.157
	0.997	0.998	0.964

The gas-surface interaction models

For the free-molecular computations, we have used the classic diffuse-specular reflection model of Maxwell (1879) to explore the range of the drag values, but other models based on based on accommodation coefficients are also available (Cercignani & Lampis, 1971). Williams (1973), had considered various models for simultaneous drag and heat transfer on a sphere, and Chernyak and Sograbi (2019) also reported results for a sphere for models of interest. The relationship between the inbound and outbound molecular distributions is given by:

In eq. (47), the symbol represents a suitable gas-surface scattering operator. Under these conditions, eq. (5) can be evaluated either analytically or numerically. Thus, for example, for the model of Cercignani and Lampis (1971), Chernyak and Sograbi (assuming and close to 1), have reported (Chernyak & Sograbi, 2019):Where and are, respectively, the normal energy and tangential momentum accommodation coefficients. For values of and not close to one, the integration of eq. (16) requires numerical techniques. We computed these integrals and found that the above result is correct even for values of and close to 0.1 (within 2% or so of values obtained from numerical integrations). We will be reporting these numerical results/comparisons separately.

We should note that all gas-surface models are phenomenological, and most useful values for such coefficients come from experiments with careful control of surface conditions, and when the inverse problem (for extracting coefficients and fitting the theoretical results to the data) is not ill-posed. But for most applications, one deals with natural (or engineered) surfaces and particulate surface control is often unavailable.

One can make the problem less ambiguous by assuming relationships between the coefficients. This normally requires (unavailable, or at least inconvenient) experimental data, but datasets do exist in the literature. For example, one may obtain the tangential momentum accommodation coefficient and the coefficient of diffuse specular refection from Millikan's data.

Additional measurements made by Dahneke with an aerosol beam (Dahneke, 1973b, 1973c), and measurements with levitated and rotating steel spheres in gases such as , , , , , (and some mixtures of these) also provide experimental data (Bentz et al., 1997; Bentz et al., 1999; Bentz et al., 2001; Tekasakul, 1996).

All these experiments suggest that assuming a value for the accommodation coefficient, α, in its upper range (0.8–0.95) is proper for commonly encountered gasses. For example, the rotating sphere measurements showed values in the range of 0.88–0.97 for , 0.85 to 0.93 for , 0.80 to 0.95 for , and finally 0.83 to 0.89 for nitrogen, all at STP.

Larriba-Andaluz and Carbone (2021) have recently reviewed values inferred from experiments including those with differential mobility analyzers and other means, noting a general range of 0.8–0.9, with the values increasing as the particle size increases.

There remain interesting possibilities for analyzing the existing experimental data of Millikan and others (See Rader (1990) for a compilation for several gases), further with the use of the Cercignani-Lampis model in that the two accommodation coefficients in there can be determined by using the free molecular and the slip limits (these will provide two equations). We have noted the free molecular limit result above in eq. (26). Accurate results for the velocity slip coefficient for the Cercignani-Lampis model and general intermolecular interaction laws can be obtained for example, from the variational and similar general expressions given by Loyalka (1971), and further discussed in Ivchenko et al. (2007).

Since in the near term at least, there will remain uncertainties with the gas-surface interaction models and parameters, the practical solution is to have theoretical results available for a range of values of these parameters and assess the sensitivity of particle mobility and overall aerosol dispersion/transport and ES-DMA analysis to these variations.

It is also pertinent to note here that in interpretations of experiments, analysts often ignore wall effects. While small in absolute size, the nature of the confined (chamber) geometries may lead to these effects becoming significant. Thus, extensions of the present computations to particle diffusion in confined geometries will be important.

Experimental Measurements

We have discussed measurements of gas-surface interaction parameters above. As we have noted in the introduction, one can measure particle mobility via ES-DMA, and such measurements with well-defined virus like particles would be of considerable value.

There have been measurements on non-spherical particles with improved Millikan cells also, as discussed by Cheng and colleagues, (Cheng et al., 1988) and it would again be of interest to measure drag and hence the mobility of virus like particles in such cells. Levitation experiments, with a rotor gauge or similar device, to measure frictional torque on such particles should also be of interest. One could also use electrodynamic balances as described by Davis and Schweiger (Davis and Schweiger, 2002).

We should also note that in the continuum regime, one can additionally measure drag on particles via settling experiments (in high viscosity fluids), and for common cases, computations can be both verified and confirmed.

Summary

Our aim in this paper was to elucidate the role of virus structure on virus drag, diffusion, and mobility. The methods and the techniques considered also potentially apply to interpretation of virus structures, as well as virus separations through use of ES-DMS techniques and devices. We have considered model virus mimetic geometries (composed of spheres) and obtained free- molecular and continuum results. In the free molecular limit, we have used the Maxwell diffuse-specular reflection model and a Monte Carlo technique, and our results are in good agreement with those obtained from the IMoS program. In the continuum limit, we have solved the boundary integral equation associated with the Stokes equation by using a singularity subtraction and collocation technique and obtained satisfactory results. In both the free molecular and the continuum limits we have found that a volume equivalent sphere approximation gives satisfactory results for the geometries considered in this paper. We have described construction of approximate expressions, using the free molecular and continuum results, which are applicable over the entire size range of the virus size as compared to molecular mean free path (the Knudsen number). We have also discussed gas-surface interaction models and values of the accommodation coefficients reported in the literature. We have further briefly noted the experiments that would be useful in supplying useful data and verification of the theoretical results.

We will be reporting results for more complicated virus structures in the future as suggested in Fig. 3, in aiding understandings of virus and other biological particle transport, both in gases and in liquids. We would be particularly attentive to the applications of this work, and extensions of it, to virus structure determination and separations via the ES-DMA technique.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.