Choice of computational method for swimming and pumping with nonslender helical filaments at low Reynolds number Free

Low Reynolds number flows induced by filaments are commonly encountered in microscale scenarios. They occur in the swimming of microorganisms such as bacteria and spermatozoa,¹ as well as artificial flagella.^2–6 Furthermore, cilia- and flagella-induced microscale pumping is used in applications such as particle transport in biological systems^7,8 and the regulation and transport of fluids by fabricated systems utilizing bacterial flagella.⁹ Filament-induced flows have been modeled using Resistive Force Theory (RFT), slender body theory (SBT), and boundary element methods (BEMs). For the slender filaments encountered in flagellar (bacterial and eukaryotic) swimming, there is some guidance in the literature as to which of these methods to use. Resistive force theory has been found to be inaccurate relative to slender body theory results,¹⁰ while slender body theory and boundary element methods have compared favorably.^11,12 Recently, it has been shown that for slender helical filaments, RFT also gives misleading results in comparison to experiments.¹³ However, not all helical filaments of interest are as slender as flagella, and there is less guidance available for the choice of method for thicker filaments. Thicker helical filaments can be realized in the helical cell bodies of many motile bacteria and are also important for magnetically actuated microscale swimming used for microrobotics, which could be used in medical applications such as in vivo surgery.^2,3,6,14,15

While accurate and efficient computational models are needed for these filament-induced flows, it has not been established that which methods are appropriate for a given geometry given the wide range of filament thicknesses and helical pitches encountered, making it difficult for practitioners to make an informed choice without extensive numerical studies. Here, we compare the popular RFT, SBT, and two other methods based on the use of regularized Stokeslets. Each has different advantages and disadvantages in terms of ease of use, accuracy, and computational cost. Our objective is to evaluate the accuracy and computational requirements of these methods across the range of geometries encountered in microscale applications in order to provide a guideline which may be used to make an informed choice of which numerical method to use for a given filament geometry.

What helical geometries are observed in the study of biological and artificial pumping and swimming applications? Consider a rigid helical filament whose geometry depends on three dimensionless parameters (see Fig. 1(a)): the ratio of the filament radius to helical radius $ã=a/R$⁠, the ratio of helical pitch to helical radius $p̃=p/R$⁠, and the number of turns. In this study, we vary the parameters $ã$ and $p̃$ to encompass the range of geometries encountered across most applications. Because Rodenborn et al.¹³ find that the accuracy of SBT and the method of regularized Stokeslets are not dependent on total arclength, in this study we keep the arclength fixed at 3 pitch lengths.

Bacterial flagella are the thinnest filaments in our study with a radius of a = 0.012 μm.¹⁶ For bacteria such as Salmonella and E. coli, the flagella may undergo polymorphic transformations between helical forms categorized by the number (n) of “R”-type protofilaments in the flagellum.^16–18 These polymorphic forms have a wide variation in both $ã$ and $p̃$⁠. For example, the normal (n = 2) form has values of $ã=0.023$ and $p̃=1.611$⁠, while the curly I (n = 5) form, seen in flagellar unbundling, has values of $ã=0.029$ and $p̃=8.724$⁠.¹⁸ In Fig. 1(b), we have noted $ã$ and $p̃$ for the n = 1, 2, and 9 polymorphic forms. These different bacterial flagellar forms are important not only for swimming and chemotaxis but also in micropumping applications. While we do not take into account the fluid structure interactions in this study, our results are relevant to any approach which would add deformations to the hydrodynamic interactions, such as for understanding the hydrodynamics associated with transformations between polymorphic configurations of these bacterial flagella.

At the slender end of the spectrum, we also find the flagella of various mammalian spermatozoa, whose waveforms were presented by Brennen and Winet.¹⁹ Human spermatozoa have a filament radius of a ≈ 0.25 μm, around 20 times larger than that of bacterial flagella. However, due to their larger helical radius, their values of $ã$ and $p̃$ are on the same order as those of bacterial flagella, with $ã$ ranging from 0.01 for bull spermatozoa to 0.05 for human spermatozoa, and $p̃≈10$⁠.¹⁹ 

For the slender regime of bacterial and eukaryotic flagella, Brennen and Winet note that the slender body theory originally proposed by Lighthill²⁰ or the resistive force theory of Gray and Hancock²¹ should be applicable. Although there are many well-known variants of slender body theory,^22–26 in this study, as a representative SBT, we use the implementation proposed by Higdon,²⁷ which is equivalent to Lighthill’s and which provides a relatively easy to implement numerical framework.

Rodenborn et al. focus on the slender regime up to $ã=0.0625$⁠, but thicker helical filaments are common. For example, spirochetes have helical or wavy cell bodies, such as L. illini, with $ã=0.835$ and $p̃=7.977$⁠, and B. burgdorferi, with $ã=0.868$ and $p̃=14.895$⁠.²⁸ Along with these spirochetes, thick filaments include artificial swimmers such as the magnetically actuated nanopropellers fabricated and studied by Ghosh et al.^2,3,29 with $ã=1.188$ and $p̃=5.688$⁠, and the artificial flagella studied by Peters and others^4,5 with $ã=0.483$ and $p̃=5.517$⁠. In this regime, we would expect the previously mentioned theories to fail as the filament is no longer slender.²⁷ Instead, to obtain accurate results, typically one must employ a BEM.³⁰ As a simple implementation of a BEM, we utilize the method of regularized Stokeslets, which covers the surface of an object with continuous analogues of the classical Stokes flow singularities^14,31–34 and has gained popularity due to its ease of use. Recent work on swimming sperm³⁵ and the beating patterns of cilia⁸ have also distributed regularized Stokeslets along filament centerlines rather than at the surface, which yields significant savings in computational effort, so we also examine such centerline distributions of regularized Stokeslets. Optimal parameters have been found for centerline distributions for slender ellipsoids,³⁶ which may be used to model straight rods, but not for helical filaments.

To summarize the range of biological and artificial helical filaments encountered in microscale applications, in Fig. 1(b), we plot a number of the examples mentioned above in $ã$ and $p̃$ space. In the following investigation, we will examine geometries with $ã$ between 0.05 and 0.95 in increments of 0.1, and $p̃$ between 10 and 50 in increments of 10, which encompass most of the physical and artificial microswimmers seen in Fig. 1(b).

Our investigation evaluates the accuracy of resistive force theory, the centerline distribution of regularized Stokeslets, and the slender body theory for two physical scenarios by comparing to benchmark results obtained from surface distributions of regularized Stokeslets, which are verified by convergence studies. The first scenario we consider is that of a swimmer with an ellipsoidal head and helical flagellum (seen in Fig. 2). The swimmer is force-free (F = 0) and is driven by setting a constant relative angular velocity ω between the body and flagellum. The outputs we use to evaluate accuracy in this scenario are the translational velocity V and the total torque N_f on the flagellum. The second scenario we consider is a pump created by a rotating helical filament attached to a no-slip plane (seen in Fig. 2). Since the filament is fixed to a plane, the translational velocity is V = 0, and to drive the pump, we apply a constant torque N to the filament. The outputs used to evaluate accuracy are the total force on the filament F and the angular velocity Ω. These two scenarios differ enough that they may be used as representatives for a broad range of applications.

Finally, we combine evaluation of the accuracy of resistive force theory, slender body theory, and the centerline distribution of regularized Stokeslets with an analysis of the relative computational cost for each method. The end result is a framework that can guide the choice of an optimal numerical method for a given filament geometry, which we hope will be broadly useful to researchers encountering flows caused by helical filaments in microscale applications.

In practice, SBT and methods involving regularized Stokeslets involve a similar procedure: the flow is created by an unknown force distribution, either on a surface for the method of regularized Stokeslets or on the filament centerline for slender body theory and centerline regularized Stokeslets. No-slip boundary conditions are imposed by setting the flow velocities at a set of collocation points equal to the velocity specified from the kinematics of the object. The result is the formation of a linear system of equations which may be solved to find the force distribution and hence velocity field everywhere. In equations, if r⁰ is the origin of the body frame, V is the translational velocity, and Ω is the angular velocity about r⁰, then the velocity at point α is given by

where Δr^α = r^α − r⁰. Since this is a linear equation, we may write it in terms of the 3 × 6 matrix k^α,

This allows us to represent the velocity at the collocation points of the system as

where K ∈ ℝ^3N×6. Note that we use a capital letter U to represent the 3N × 1 vector containing the N velocities of all collocation points and similarly a capital K which is composed of all N k^α. To satisfy no-slip boundary conditions, the velocity specified by the kinematics is set equal to the fluid velocity at the collocation points. Thus, the fluid velocity vector U can also be written as

where G ∈ ℝ^3N×3N and P ∈ ℝ^3N×1 contains the 3 × 1 force p^α at each collocation point. We call the submatrices g^α,β ∈ ℝ^3×3 the kernel. The kernel represents the velocity at point α due to a force at point β. Each numerical method in our study can be specified by choosing a different kernel.

For the case of a rotating helical filament attached to a no-slip plane, we choose r⁰ = 0 at the attachment point at the plane surface (see Fig. 2) and represent the total force and torque on the body as

which can be written as

Enforcing the no-slip boundary condition on the surface of the filament lets us to form a relationship between the velocities V and Ω, and the total force F and torque N on the filament. Solving Eq. (4) for P = G⁻¹U and using Eqs. (3) and (7),

where R = K^TG⁻¹K is the resistance matrix with 3 × 3 translational, rotational, and coupling components denoted by R_T, R_R, and R_C, respectively. We apply a constant torque N tangential to the axis of the filament and set the translational velocity V = 0 since the filament remains anchored above the x = 0 plane. These constraints yield the system

with solution

For the swimmer with an ellipsoidal head, we choose the origin of the body-fixed frame r⁰ as the flagellum attachment point on one end of the head (see Fig. 2). A surface distribution of regularized Stokeslets is always used to model the head, and these collocation points are denoted as $rhα1$ for 1 ≤ α₁ ≤ N₁. The collocation points on the flagellum are denoted as $rfα2$ for 1 ≤ α₂ ≤ N₂. In this scenario, V and Ω represent the translation and rotation about the origin. Here, the kinematics specify

where ω = (ω, 0, 0)^T is the specified relative rotation rate between the head and flagellum and

As in Section II, this is a linear system of equations and can be written as

where K ∈ ℝ^{3(N₁+N₂)×6}, U_h ∈ ℝ^3N₁×1 is the velocity at each collocation point on the head, U_f ∈ ℝ^3N₂×1 is the velocity at each collocation point on the flagellum, and $L̇∈R3(N1+N2)×1$ is

On the other hand, the hydrodynamics specify

where G_hh ∈ ℝ^3N₁×3N₁ and G_fh ∈ ℝ^3N₂×3N₁ give the velocity on the head and flagellum, respectively, due to the point forces $phα1$ located on the head, and G_hf ∈ ℝ^3N₁×3N₂ and G_ff ∈ ℝ^3N₂×3N₂ give the velocity on the head and flagellum, respectively, due to the point forces $pfα2$ located on the flagellum. Enforcing no-slip boundary conditions requires that the kinematic and hydrodynamic velocities to be equal, which yields the equation

In this scenario, we also require that the net force and net torque be zero, so Eq. (7) implies

This equation can be solved for V and Ω to obtain

where $M=KTG−1K−1$ is the inverse of the resistance matrix, called the mobility matrix. In this scenario, the translational velocity V and the torque on the flagellum N_f are the outputs. To solve for the torque, we find P_f through Eq. (16) since V and Ω are now known, giving

We also examine the accuracy of resistive force theory for the case of the swimmer. In resistive force theory, long-range hydrodynamic interactions are ignored, so we treat the hydrodynamics of the flagellum and head separately through their resistance matrices. The resistance matrix of the ellipsoidal head R_h may be calculated from analytical results, e.g., as presented in the textbook of Kim and Karrila.³⁷ Resistive force theory is used to calculate the resistance matrix of the helix as follows. The local force per unit length on the filament is

where v_rel = V + Ω × r, $tˆ$ is the local unit tangent vector, and the coefficients ξ_∥ and ξ_⊥ are drag coefficients obtained by slender body analysis. For the theories of Gray and Hancock²¹ and Lighthill,²⁰ these coefficients are

Note that in Lighthill’s formula, we use the pitch as suggested by analysis of helical filaments.³⁸ The Lighthill coefficients might be expected to be more applicable since in a comparison of RFT and SBT for planar thin filaments, Johnson and Brokaw¹⁰ found that the Lighthill coefficients adequately predicted local forces and torques on the filament, while the Gray and Hancock coefficients required ∼35% adjustments to produce comparable accuracy. Integrating the force and torque r × f_RFT along the filament for various translational and rotational velocities yields the resistance matrix R_f of the flagellum. The origin of rotation for both R_h and R_f is the flagellum attachment point.

Given the resistance matrices, we require that the total force and torque be zero,

and using Eq. (8), we may express the above equation as

This equation is solved for V and Ω to give

which gives us one of our outputs, the translational velocity V. Using the relationship

we find our second output, the torque on the flagellum N_f.

In each of the physical scenarios described in Section I, we consider a helical filament oriented in the positive x-direction of non-dimensional radius $ã$ and non-dimensional helical pitch $p̃$ whose arclength is fixed at $L̃=L/R=3p̃$⁠. The filament centerline is described by

The exponential term describes a tapering region near the attachment point, where $kE=p̃/(2π)$ is the characteristic length of the tapering region.¹⁴ 

For the helical swimmer, we use an ellipsoidal cell body with semi-major and semi-minor axes which are $ãh=ah/R=10.204$ and $b̃h=bh/R=5.102$⁠, respectively. The major axis is placed along the x-axis, and a helical flagellum whose geometry is prescribed in Eq. (25) is attached at x = 0. The helix is translated away from the ellipsoid to create a gap between the head and flagellum in order to ensure accurate results for nearby moving surfaces.¹⁴ An example of this geometry is seen in Fig. 2.

Here, we examine the rotation of a helical filament with its end at a distance of $h̃=10ã$ above a no-slip plane (Fig. 2). To model this scenario, we translate the centerline prescribed above by $h̃$ in the x-direction. This setup is inspired by pumping applications where pumping flows may be caused by the movement of filaments.

For the benchmark cases, we employ the method of regularized Stokeslets by distributing smoothed analogues of the classical Stokeslet on the filament surface or body.^31–33 In the following calculations, we make use of a standard smoothing function which depends on a “blob” parameter δ that describes the spread of the point forces and allows us to obtain non-singular numerical values for all calculations. A full description of the regularized Stokeslet is found in the literature,^31–33 and the kernel g^α,β (defined in Section II) that we use for the method of regularized Stokeslets is given by Eq. (10b) in the work of Cortez et al.³³ 

For discretization, the filament is divided into cross sections perpendicular to the centerline of the filament. We place n_c evenly spaced Stokeslets on the perimeter of each cross section, similar to the discretization described previously by Hyon et al.¹⁴ The cross sections are spaced so that the distance between Stokeslets (Δs) on a given cross section is the same as the distance between successive cross sections. Therefore, $Δs=2ãsin(π/nc)$⁠. This spacing means that the total number of Stokeslets increases rapidly as we increase the number of Stokeslets per cross section and as the filament radius decreases. For a given spacing, the blob parameter δ is set equal to the distance between Stokeslets, δ = Δs.

For the discretization of the swimmer with an ellipsoidal head, each successive filament cross section is rotated by an angle π/n_c with respect to its normal axis in order to stagger the Stokeslets. Also, the ends of the filament are capped with hemispheres that are discretized using the grid method given by Cortez et al.³³ The ellipsoidal head is discretized by cutting into cross sections perpendicular to the major (fore-aft) axis. If we define Δs_h as the spacing between Stokeslets on the cell body and successive cross sections, the number of Stokeslets on each cross section depends on the radius of the cross section and is determined by n_c = 2πr_h/Δs_h, where $rh=b̃h1−(x/ãh)2$⁠. Similarly to the flagellum, each successive cross section of the head is rotated by an angle Δs_h/2r_h to stagger the Stokeslets on the cell body. For our investigation, Δs_h = 0.1 is the non-dimensional Stokeslet spacing on the head, and N₁ = 750 is the total number of regularized Stokeslets on the head.

The case of the filament rotating above a no-slip plane uses a slightly different discretization. The cross sections are not rotated about the normal axis of the filament, and there are no hemispherical caps added to the end of the filament. Most importantly, the kernel is modified to satisfy the no-slip boundary conditions at the planar boundary (x = 0). For the method of regularized Stokeslets, we employ a full regularized Blakeslet distribution which consists of the regularized Stokeslet and its corresponding image system mirrored about the x = 0 plane:³¹ for a collocation point on the filament surface at r^α = (h^α, y^α, z^α), its image is located at a point $r*α=(−hα,yα,zα)$⁠. In this scenario, the kernel g^α,β is derived from Eq. (21) in the work of Ainley et al.³¹ 

The flows around filaments have also been computed using a centerline distribution along the center of the filament instead of at the filament surface. For this method, we use the same kernels as in Section II E 1. For the swimmer, the head is always modelled using a surface distribution of regularized Stokeslets with spacing and blob size as defined in Sec. II E 1.

For both the swimmer and the pump, we discretize the filament centerline by placing regularized Stokeslets with blob size δ at equally spaced points along the filament centerline with spacing Δs. There is no a priori choice of Stokeslet spacing Δs or blob size δ, although physically the blob size is expected to be of the same order of magnitude as the filament radius. These parameters must be tuned for each geometry to obtain the most accurate result. We report the optimal parameter values for each geometry we study later; these provide starting points for studies determining the optimal parameter values for other geometries in future studies. Despite this, the positive aspects of this discretization are twofold. First, it still relies on the method of regularized Stokeslets, and so we may use the same kernels found in the works of Cortez et al.³³ and Ainley et al.,³¹ which are relatively easy to implement. Second, the number of Stokeslets needed for the discretization is much lower than that of the surface distribution, increasing the computational efficiency of the method.

In our investigation, we use the slender body theory outlined by Higdon²⁷ as a representative SBT. This method calculates the flow due to a distribution of classical Stokeslets and dipoles along the centerline of the filament. The force distribution f along the centerline is unknown, while the dipole strength $d=−(ã2/4μ)f⊥$ is specified to preserve no-slip boundary conditions on the surface of the filament. In Higdon’s method, the filament centerline is discretized into small cylindrical segments of length 2q, and the Stokeslet (and hence dipole strength) is assumed to be constant within each segment. The integral of the Stokeslet and dipole distribution along each segment can be evaluated analytically to yield the flow produced by each segment, which defines the kernel. The expression for the kernel we use can be found in Eq. (27) in the work of Higdon.²⁷ Similarly to the centerline distribution of regularized Stokeslets discussed in Sec. II E 2, the parameter q must be tuned to ensure the most accurate results, but the advantage is far fewer discretization elements than a surface discretization, especially for thin filaments.

For the swimmer, we always model the ellipsoidal head using the surface distribution of regularized Stokeslets described in Section II E 1. Therefore, the G_hh and G_fh submatrices in Eq. (16) are always formulated in terms of the regularized Stokeslet kernel discussed in Section II E 1. The other submatrices G_hf and G_ff are formulated using the kernel for the slender body theory.

For the case of a helical filament above a plane, we must include an image system to satisfy the no-slip condition on the plane. The singularities distributed along the centerline of the filament are Stokeslets and dipoles, and so the system is not simply a Blakeslet distribution since we must include the image system of the dipole.²⁵ To evaluate the contribution of the images, we integrate the contribution from each image segment numerically using Simpson’s rule since they are non-singular at the collocation points on the filament centerline. The kernel for this method is defined by the sum of Eq. (27) in the work of Higdon²⁷ and the numerical integral of the Stokeslet and dipole image systems provided by Blake and Chwang.²⁵ 

In order to assess the accuracy of these methods, we examine an error metric which includes both output parameters, the swimming velocity V and torque on the filament N_f in the case of the swimmer, and the net force F and the angular velocity Ω in the case of the pump,

where we add the sum-squared error of each respective output. The “b” superscript denotes the output from the benchmark case using the surface distribution of regularized Stokeslets.

In previous work,¹⁴ we have shown that for geometries similar to those used here for swimming, the method of regularized Stokeslets agrees well with SBT and boundary element method results presented in the work of Phan-Thien et al.¹¹ For this study, we create benchmark results using the surface distribution of regularized Stokeslets and perform convergence studies for selected geometries for swimming and pumping. These sample cases correspond to $p̃=20$ and $ã=0.05,0.35,0.65,0.95$⁠. For the swimmer, we examine the x-component of the non-dimensionalized translational velocity against the Stokeslet spacing Δs. For the pump, we examine the x-component of the non-dimensionalized force against Δs.

Fig. 3 presents the convergence studies for the swimmer and pump, respectively. In both cases, the quantity examined (either velocity or force) is normalized by a Richardson extrapolant obtained as described below, as a function of number of regularized Stokeslets used in the discretization. The extrapolant d_e for each case is obtained from the expectation that the error for each data point d_i obeys a power law relation to the discretization spacing Δs,

where C and γ are constants. In each case, we examine the data for which γ is in the asymptotic range to obtain a fit for C and γ in the least-squared sense.

Ainley et al. find slightly greater than first order convergence when using a blob size δ ∼ Δs^0.9 for the problem of a sphere translating near a wall.³¹ Since we use the similar δ = Δs, we expect similar rates of convergence, and indeed in nearly all the cases, we find that γ ≈ 1. To provide estimates of error, in the insets of Fig. 3, we show the fits to d_i/|d_e| vs. Δs for γ = 1 over the asymptotic range. The swimming and pumping geometries with $ã=0.05$ (red crosses) display higher order behavior. In those cases, for consistency, we still find the extrapolant using γ = 1, which yields a conservative estimate of error. For the swimmer, the non-dimensionalized x-component of velocity for every geometry has converged to within 4% of the linear extrapolant. For the pump, the non-dimensionalized x-component of force has converged to within 5% of the linear extrapolant.

Based on these convergence studies, benchmark results for all the other geometries are generated using 8300 − 12 000 regularized Stokeslets for the swimmer and 10 000 − 12 500 regularized Stokeslets for the pump, for which we expect an error of less than 5%. These convergence studies also highlight the computational cost of this method and the motivation behind our investigation of other computationally cheaper numerical methods. The computational cost for all methods scales as N³, where N is the total number of collocation points. As we will see in Secs. III B and III C, the method of regularized Stokeslets distributed on the filament centerline as well as the slender body theory provides a significant reduction in the total number of collocation points. This reduction comes at the cost of reduced accuracy, and so we examine the error due to each of these methods to assess their utility in various physical regimes.

We now investigate the accuracy of the computationally cheaper methods outlined in Section II by comparing the results with the benchmarks obtained from calculations using the converged surface distribution of regularized Stokeslets. First, we present findings for the method of regularized Stokeslets distributed on the filament centerline. For each of the 50 filament geometries, we minimize the total percentage error metric E as defined in Section II F with respect to the blob parameter δ and the Stokeslet spacing Δs. Tables I and II report the value of these optimal parameters. From the tables, we see that for those geometries where the error can be made small (<10%), the optimal blob size is typically between $ã$ and $3ã$⁠, while the spacing is similar to the blob size. In this regime of accuracy, both the optimal spacing and blob size tend to increase as $p̃$ increases. Our result for optimal blob size is consistent with the results of Bouzarth and Minion,³⁶ who found an optimal blob size of 1.43a for slender prolate ellipsoids with aspect ratios >10, which might be used to model straight rods.

Fig. 4 shows the total error for each geometry for pumping and swimming, respectively. In these figures, the circle size is scaled according to the percent error, and interpolated contour lines are drawn at 5, 10, 20, and 30 percent error in order to give a visual representation of the regions of accuracy. From the results, we make the following observations. First, the centerline distribution of regularized Stokeslets is more accurate for the swimmer than the pump. Second, this numerical method is most accurate for thin filaments (where $ã$ is relatively small) and not accurate for the thickest filaments examined. For example, comparing to Fig. 1(b), this numerical method is quite accurate (with errors less than 5%) in the modeling of the flagella of bacteria such as Salmonella and E. coli as well as spermatozoa, but not accurate (with errors of approximately 25% in each scenario) for spirochetes such as L. illini and B. burgdorferi.

We perform a similar investigation by comparing results for the slender body theory outlined in Section II to the benchmark results. In this numerical method, the total error is minimized with respect to the parameter q (recall that the length of the individual cylindrical segments is 2q). Tables I and II report the optimal value of q for each geometry. From the tables, we see that accurate results are produced when both $q/ã>5$ (the segment length is at least 5 times the diameter) and N ≳ 45 (≳15 segments per turn). For slender body theory, we only report error values for the first 25 geometries since the error grows rapidly as a/R increases. This is consistent with expectations for this method: as the filament radius approaches the helical radius, the filament is less slender, and so it follows that a slender body theory approach would break down.

Fig. 5 shows the total error for each geometry for pumping and swimming, respectively. The total error in the slender body theory is within a tolerable range (<10%) only for $ã≲0.10$⁠. We note that the total error using the method of regularized Stokeslets on the filament centerline is also quite low for these particularly slender geometries, so in this range, we must make a judgement based on computational cost to compare the two methods.

Table III shows the results of RFT using the drag coefficients from the theories of Gray and Hancock²¹ and Lighthill.²⁰ Throughout the whole range, neither gives accurate results; for the coefficients derived by Lighthill, the total percent error is minimal at E_Total = 26% and for the coefficients derived by Gray and Hancock, the total percent error is minimal at E_Total = 33.6%, both for the geometry specified by $ã=0.05$ and $p̃=50$⁠.

Although the total error is large throughout the range tested, it is interesting to note that if one only considers the error in swimming speed (separately reported in Table III), the coefficients derived by Gray and Hancock produce accurate results for the swimming speed (E_V < 10%) for $ã≤0.65$ and $p̃≥30$⁠, a large range of geometries. In this regime, however, we find large errors in the filament torque. On the other hand, the theory of Lighthill produces accurate results for the filament torque for $ã<0.15$⁠, but always gives large errors for the translational velocity. Thus, RFT may give accurate results for individual quantities, as has previously been observed in the literature, not only for swimming speeds but also for quantities such as distance traveled per revolution,¹⁵ but does not accurately capture the full dynamics for any investigated geometry. Furthermore, without a more accurate benchmark comparison such as that performed here, it is difficult to know in advance if any particular quantity of interest will be accurately calculated using RFT.

To summarize the results up to this point, based on a convergence study, we created benchmark results for swimming and pumping flows caused by helical filaments using a surface distribution of regularized Stokeslets. By comparison to the benchmark results, we found that resistive force theory does not give accurate results (in terms of total error) for any of the geometries studied. We also used the benchmark case to optimize the relevant parameters of a centerline distribution of regularized Stokeslets and slender body theory with respect to a total error metric and obtained regimes of accuracy for each numerical method. Together, Figs. 1, 4, and 5 allow us to make a judgement about which numerical method is appropriate for a wide variety of biological and artificial filament geometries. The regimes of accuracy for more cost efficient methods (SBT and centerline distributed regularized Stokeslets) overlap, so it is of interest to examine which numerical method will be optimal in terms of computational efficiency as well as total error. The computational cost of the numerical methods scales as N³, where N is the number of collocation point. As noted before, a surface distribution of regularized Stokeslets (or any boundary element method) requires more and more collocation points as the filament becomes thinner since the spacing of the cross sections becomes smaller, but for slender filaments, SBT and centerline distributions of regularized Stokeslets are accurate and much cheaper computationally.

The overall choice for the cheapest numerical method which is still accurate can be summarized in Fig. 6. In this figure, we first choose an error threshold of 10% against the benchmarks obtained from the surface distribution of regularized Stokeslets. A contour is drawn for an interpolated 10% error threshold for each of the cheaper methods, which delineate three regions labeled near the bottom of each plot as “SBT” for the slender body theory, “CL” for the centerline distribution of regularized Stokeslets, and “SURFACE” for the surface distribution. These regions show the method which, given the 10% error threshold, is the most computationally efficient for a given geometry. Note that the regions are similar for the swimming and pumping scenarios, which suggests that these recommendations are broadly appropriate across many applications.

In order to provide information about the computational costs in each region, in Fig. 6, we also write the number of collocation points for each geometry tested. For the surface distribution region, we report the number of collocation points obtained by decreasing the number of collocation points for each geometry by lowering the number of regularized Stokeslets per cross section n_c until we reach the error threshold of 10%. Points marked with an asterisk denote geometries for which the error of the surface distribution against the benchmark case was always less than 10%, even for the lowest number of regularized Stokeslets per cross section (n_c = 3). In these cases, the value of N using n_c = 3 is reported. This number allows us to compare to the costs required for the other methods. For the centerline distribution and SBT regions, we report the number of collocation points needed for the optimal accuracy parameters.

Within its regime of accuracy, we see that slender body theory requires approximately half as many collocation points as the centerline distribution of regularized Stokeslets. This means that the slender body theory may be up to ∼10 times more computationally efficient in these regimes, which include more slender forms of bacterial flagella as well as spermatozoa. This result is consistent with the experimentally confirmed results of Rodenborn et al.¹³ It should be noted that due to the addition of the dipole singularity (as well as the corresponding image system in the case of the pump), this method requires a larger initial investment for implementation than the method of regularized Stokeslets. Taking this into account, it becomes a value judgement in this regime as to whether the computational savings will be worth the initial effort in the long run.

For $0.15<ã<0.40$⁠, the centerline distribution of regularized Stokeslets provides a very efficient model for a wider range of physical examples. This range may include some bacterial flagella which are enveloped by a sheath, which makes them thicker than a bare flagellum, such as those of H. pylori³⁹ or Vibrio species.⁴⁰ Finally, as expected for the thickest filaments, such as the cell bodies of spirochetes or artificial helical microswimmers, the surface distribution of regularized Stokeslets is required for accurate results but increases the required number of collocation points by at least an order of magnitude.

In conclusion, the numerical methods tested — resistive force theory, slender body theory, and the method of regularized Stokeslets — have gained popularity due to their ease of implementation. We investigated their use throughout the range of geometries likely to be encountered in applications. In Tables I and II, we present optimized parameters for each method we study. Although in other geometries these precise parameters may not be optimal, they should provide useful starting points for tuning the parameters in other applications and geometries. For the case of a 10% accuracy threshold, Fig. 6 shows how an appropriate choice of method can be made for a wide range of filament geometries. Thus, our results provide necessary guidance for selection of computational method and should be useful for many future investigation involving filament induced flows.

This work was supported by National Science Foundation Award Nos. DMR-1307497 and CBET-1252182.