Helical Swimming in Stokes Flow Using a Novel Boundary-element Method
We apply the boundary-element method to Stokes flows with helical symmetry, such as the flow driven by an immersed rotating helical flagellum. We show that the two-dimensional boundary integral method can be reduced to one dimension using the helical symmetry. The computational cost is thus much reduced while spatial resolution is maintained. We review the robustness of this method by comparing the simulation results with the experimental measurement of the motility of model helical flagella of various ratios of pitch to radius, along with predictions from resistive-force theory and slender-body theory. We also show that the modified boundary integral method provides reliable convergence if the singularities in the kernel of the integral are treated appropriately.
The boundary-element method (BEM) Pozrikidis:1992td () is a reliable and accurate tool for studying the zero-Reynolds number hydrodynamics of motile microorganisms Lauga:RPP2009 (), especially in situations where the detailed geometry of the structure of the microorganism plays a role Phan-Thien:1987 (). When properly formulated, the BEM leads to robust solutions with high orders of convergence Power:1987 (); Gonzalez:2009 (); Keaveny:JCP2011 (). However, the trade-off for high reliability is that the computational cost of this method is typically high, increasing drastically with the number of mesh elements used for spatial discretization. In practice, this heavy computational load can be reduced either by approximating the singular kernel in a regularized form, as in the method of regularized Stokeslets Cortez:2001 (); Cortez:PF2005 () or by replacing the force density with a piecewise constant function Smith:2009kh (). With these approximations, fewer nodes are required for a given boundary geometry. However, to achieve robust convergence, regularizing the singular kernel requires a delicate selection of the numerical parameters describing the width and spacing of the regularized Stokeslets Cortez:PF2005 (); Bouzarth:2011 (). Likewise, there is no natural prescription for assigning the patches of constant force, and clever choices of the constant-force patches are often necessary Smith:2009kh (). However, small motile organisms often have symmetric bodies or ciliary beat patterns, and these symmetries have not always been exploited in computation. For instance, bacteria such as Eshericia coli swim by rotating helical flagella Berg:1973td (), and the spirochete Leptonema illini has a helically shaped cell body Charon:1984 (). These structures typically have many pitches, and thus approximate perfect helical symmetry. Another example is the array of beating cilia on an actively swimming Paramecium, which coordinate to form metachronal wave pattens. This wavefront follows a counterclockwise gyration along the cell body Horridge:1969vt (), again leading to an approximate helical symmetry.
Symmetries in fluid-structure interaction problems have been used by others to reduce the number of unknowns and thus simplify numerical computation. For an axially symmetric low-Reynolds-number swimmer Purcell:1977 (), the Green’s function in a boundary-integral method is modified to include integration along the azimuthal direction Pozrikidis:1992td (), and the two-dimensional (2) boundary-integral problem can thus be reduced to a one-dimension (1) problem Thaokar:2007 (); Spagnolie:2010 (). Helical symmetry has been applied in Lighthill’s slender-body-theory calculation of the motility of rotating flagella Lighthill:1976 (); lighthill1996b (), where every point on the helical body-centerline is regarded as identical. More generally, it has been shown that the 3 Navier-Stokes equation for flow with helical symmetry can be reduced to a 2 problem Childress:1989 (). This technique was used to study flow within a helical pipe Childress:1989 (); Zabielski:1998 (), and the helical wake of a rotating propeller at large Reynolds numbers Delbende:2011 ().
In this study we show, in general, that symmetry in the surface domain leads to a reduction in dimension of the boundary-integral method. More specifically, for a helical symmetry, a 2 boundary integral can be reduced to 1. The idea of this dimension reduction is straightforward: if we know the force vectors on an arbitrary circumference around a helical body, we know immediately the force distribution on the entire surface of the body, since the whole surface can be reconstructed from these identical circumferences. The computational benefit from exploiting this symmetry is evident. Suppose the surface of the helix is approximated by a mesh with nodes around each circumference, and circumferences, for a total number of mesh points of . In a full boundary-element method, we are required to invert a matrix. By exploiting the helical symmetry we reduce the problem to that of inverting a much smaller matrix.
As an example application, we use this modified boundary-element method to model the motion of a bacterium propelled by a helical flagellum, and compute the swimming speed of a rotating helix subject to zero force. We demonstrate that our numerical method is highly reliable, and achieves third-order convergence with proper singularity reduction. We also compare quantitatively our simulation results with some experimental studies on the motility of model helical swimmers, along with those predicted from other theoretical tools, such as resistive force theory and slender-body theory.
Ii Helical Symmetry and the Boundary Integral Method
In the boundary-integral method, the flow field of a Stokes fluid can be represented in terms of two boundary distributions Pozrikidis:1992td (), involving the Green’s function , the associated stresslet , and the force density on the domain ,
The first distribution on the right-hand side of Eq. (1) is the single-layer potential, and the second distribution is the double-layer potential. The Green’s function , also known as the Stokeslet, is
and the stresslet is
We are interested in the motion of a rigid particle, in which case the velocity may be written solely in terms of the single-layer potential kimbook ():
The velocity field is linear in the force density, . For a prescribed velocity field on the meshed surface, (, 2, , N), the force density, , can be determined using the discretized version of (Eq. 4):
where the matrix is the matrix
and the ’s are matrices of Eq. (2). Here is the area occupied by the mesh element at position , and the indices = 1, 2, , . The matrix is typically not sparse, and we must solve linear equations to find the force density. The number of equations we are required to solve increases linearly with the number of nodes on the meshed surface. However, if the surface possesses symmetries, we may reduce the size of the linear system dramatically.
Let us consider a system with helical symmetry, such as a long helical filament that rotates and translates along its axial direction in a Stokes flow. A segment of the filament is shown in Fig. 1. The axis of the helix is along . Consider an arbitrary circumference of the filament, defined by intersecting the surface by the cross-section normal to the body centerline. Note that these cross-sections at different points of the helix are related by rotations about , since, as we review below, the normal and binormal vectors of the Serret-Frenet frame Powers2010 () lie in these cross-sections, and the Serret-Frenet frame of a helix rotates about as the arclength increases. Now suppose is the force-density vector located at an arbitrary position . Because of the helical symmetry, there exists a point on the circumference such that the force density satisfies
where is rotation about axis by . The angle is given by
where is the angle between the axis and the projection of the normal onto the - plane (Fig. 1). The dashed curve in Fig. 1 shows a contour that goes through and exhibits the helical symmetry. Force densities along this contour can always be computed from the force-density vector at , using equation Eq. (7). Moreover, the angle is the phase of the helical wave along the body-centerline, and varies by over one period of the filament. This more generalized notion of is useful for determining when the surface normal is parallel to the axial direction. By applying this strategy, the force densities over the entire helical surface can be expressed as those on the circumference , and the boundary-integral formulation can be written as
Since the flow field, , also satisfies the helical symmetry, the only flow velocities we need to consider are those distributed on the circumference. Thus the integral equation for is reduced to one dimension:
Here, the tensor is our modified Stokeslet:
where the integration is performed on the helical contour that contains the point . Recall that along this contour (such as the dashed curve in Fig. 1), the force densities are of the same magnitude but have varying orientation. The modified Stokeslet depends only on the geometry of the surface, and it can be computed once the mesh is given.
Equations (10) and (11) form a complete set of equations to solve. The solution to the boundary integral problem is decomposed into two steps: (1) compute by performing the integral in Eq. (11) for each pair of nodes along , and (2) solve the linear equations Eq. (10) for . The number of linear equations is reduced to the number of mesh nodes along a single circumference, for each of the three spatial dimensions. To solve Eq. (10), we must invert the matrix , which is . This matrix is much smaller than the matrix , which must be inverted in the full boundary-element-method case. The spatial resolution along the body centerline determines the accuracy of the modified Stokeslet . The computational cost can thus be reduced by orders of magnitudes for the same number of mesh elements. It should be noted that the integrals introduced by Eq. (11) add no additional computational cost when compared to the cost of computing the matrix , since each element of need only be computed once.
Iii Body-centerline coordinates
In Sec. II, we described the general idea of how the symmetry of a helical domain can reduce the dimensionality of the boundary-integral method from three to one. In this section, we implement the idea by introducing a coordinate frame that rotates about at the same rate the cross-sections of the helix rotate. In this coordinate system, the components of the force density are the same for every point on the body centerline.
Figure 2 shows the rotating frame . This frame is the Serret-Frenet frame. We briefly summarize the properties of this frame for a helical curve. The filament body-centerline follows the path , where is the arc length along the body centerline. If the pitch of the centerline is , then is the arc length of one pitch of the centerline. Defining the pitch angle to be the angle between the centerline tangent and the axis, we find and , and therefore . Thus the centerline tangent vector is given by . The normal to the curve is given by the direction of , or , and the binormal is . The moving frame is hence defined by .
Note that the Serret-Frenet frame is obtained by rotating the space-fixed frame by about the axis, and then rotating the resulting frame by about . The angles and are Euler angles, and we may relate coordinates in the two frames by
where describes the body-centerline in the new coordinate system. The operation in (12) is a change of basis from the space-fixed frame to the body-fixed frame that accounts properly for the different choices of origins, and . The rotation operators and are
Figure 2 (b) shows a cross-section of helical filament in this new coordinate . It should be noted that such coordinate system is not unique. For a round filament, this cross-section is circular if the selected coordinate systems orient along the body-centerline.
Now let , , and denote a scalar, vector, and tensor, respectively, in the Cartisian coordinates , and let , , and denote the same in the moving coordinate system . The relations between these quantities in the two coordinate systems are given by
where angle is constant and is a linear function of . In the rotating coordinate frame, the helical symmetry eliminates the dependence on :
By substituting the force densities and velocity fields in Eq. (4) by Eq. (15), By using Eq. (4) to express the force density and velocity in terms of the coordinates , the boundary integral formula Eq. (15) becomes
where is the 2-D Jacobian arising from the change of coordinates from to . To close the problem, and in Eq. (18) are both constrained to lie on the circumference , i.e., . The modified Stokeslet is
where is the total arc length of the filament, and the tensor is the rotation operator that takes the space-fixed frame to the moving frame,
where the vector is defined as . Likewise, since , we can choose the origin of the coordinate so that
where is the arc length of the filament within each period. Thus Eq. (19) can be written as
where is the number of helical pitches (or the wave number), and is phase of the body-centerline at . To assure that the helical symmetry is a valid approximation, we use the convention that is located in the middle of the filament and .
Iv Numerical simulations of a helical swimmer
In this section we implement the boundary-integral method for a rotating helix with many turns. Now that we have the boundary-integral equations (21) and (23), the next step is to discretize them. We use an rectangular grid with points per helical pitch, so that the total number of grid points is . Using the trapezoidal rule for the integrals, we find
where , , , } is the weight function for the trapezoidal rule. Note that the sum omits the terms where the integrand is singular. These terms, with and , are accounted for in . To ensure good accuracy, the contributions to the integral from the singular parts of the integrand are computed analytically; see Appendix B. For a given rigid body with prescribed velocity, our task is to solve Eqs. (24–25) for the force per unit length by inverting the matrix . Note that without exploiting the helical symmetry, we would have to invert a matrix.
We consider two cases: (1) a tethered helix, which rotates but is prevented from translating, and (2) a swimming helix, which is subject to zero net force. Note that in both cases an external torque drives the motion. For the tethered helix, the velocity is given by . In the helical coordinate system , the helix velocity takes the form .
Figure 3 shows a typical simulation result for the force density on a tethered helix with parameters and . To ensure helical symmetry, the number of the pitches is set at . Given the spatial grids, matrices are computed using Eq. (25). Discretized force densities , , are thus obtained by inverting the linear equation (Eq. (24)). As shown in Fig. 3, converges as increases. Among the three components of the force density, the component converges fastest. The details of the convergence analysis are described in Appendix C. Using the helical symmetry, we reconstruct the force density on the entire helical surface (Fig. 3(b)). The figure shows three components of the stress tensor distribution, (, , ), with components along the normal, tangential, and bi-normal directions, respectively (inset of Fig. 3(b)). Such high spatial resolution of the stress distribution provides an accurate method in computing many mechanical features on helical propulsion, such as its net power consumption.
Using the same algorithm, we also calculate the swimming speed of a rotating helix. For a free swimmer with swimming speed , the velocity of points on the helix is
The swimming speed is determined by the condition of vanishing axial force,
The free-swimming speed exhibits similar convergence properties as those of the force densities (see Appendix C). As also demonstrated in Appendix C, by properly removing the lowest order of the numerical errors, a robust third-order convergence can be obtained.
V Comparison with experiments and theories
In this section we compare the results of our boundary-integral method for a swimming helix with experiments and the predictions of resistive-force theory and two different slender-body theories.
v.1 Experimental system
The experimental system is described in our previous work Liu:PNAS2011 (): a helical filament a few centimeters long ( cm) and millimeters wide ( mm) simultaneously rotates and translates along its axial direction in a viscous fluid. To achieve force-free swimming, we fix its rotation rate , and vary the translation speed until the net hydrodynamic force on the helix vanishes. The translation speed corresponding to is the free-swimming speed. The fluid is a high molecular weight silicone oil with kinematic viscosity St. In the regime of our experiment, the fluid is Newtonian, with the viscosity almost independent of the shear rate. For the typical rotation rate, rad/s, the Reynolds number , and the inertia of the fluid is thus negligible. As the helix is inserted in the fluid, we find that the force-free swimming speed, , does not vary with time once about one helical turn has been immersed in the fluid. Figure 5 shows the results of our measurements for two different filament thicknesses and a few different pitch angles Liu:PNAS2011 ().
v.2 Resistive force theory
Resistive-force theory, or local drag theory, is the simplest approximation for computing the force on a thin body at low Reynolds number Gray:1955 (). The force per unit length, , that a small segment of the body exerts on the fluid is taken to be proportional to the body’s local velocity, u, relative to the fluid far away, with different proportionality constants and for motion perpendicular and parallel to the body centerline (see Fig. 4):
where is fluid viscosity. The linearity of Stokes equations implies that the total force per wavelength that the helix exerts along the axis is , where and are constants proportional to viscosity . To compute , note that the force per wavelength required to pull a helix along its axis but with is
Likewise, is determined by the force required to prevent a rotating helix from translating,
The swimming speed is determined by demanding that the total force vanish, . As proposed by Gray and Hancock Gray:1955 (), the drag coefficients are and . Adopting for simplicity the coefficients in the limit , we find . Lighthill also proposed a set of drag coefficients Lighthill:1976 (), and , which were optimized for helical filaments. The predictions of resistive-force theory for as a function of are shown in Fig. 5. Note that the curves are not symmetric about , even in the limit that . The maximum speed is at a pitch angle less than because although the thrust is maximized at , the drag is minimized at , since .
v.3 Slender-body theory
We now turn to slender-body theory in which the filament is modeled by a one-dimensional distribution of singular solutions to Stokes equations. The flow at any point on the surface of the filament is given by integrating the contributions to the flow from all the other parts of the filament. Since the filament is represented by a one-dimensional distribution of singular solutions, the accuracy of slender-body theory is controlled by the aspect ratio, , with the error vanishing as .
v.3.1 Lighthill’s slender-body theory
Lighthill gave a simple physical derivation for an integral relating the velocity of a point on a filament to the distribution of forces per unit length acting on the filament Lighthill:1976 (); lighthill1996b (). The singular solutions are point forces, or Stokeslets, and source doublets. Let denote the centerline of the filament, where is arclength. For our helix, , where is the polar angle.
Since Stokes flow depends only on the instantaneous velocity of the filament, and since we consider rigid body motion, it is sufficient to consider the position of the helix at only one instant of time. Denoting the vector from one point on the helix to another by , Lighthill’s slender-body theory formula for the velocity of a point on the centerline of the filament is
where is the part of that is normal to the filament centerline, and the short distance cutoff, Lighthill:1976 (); lighthill1996b () , where ‘e’ is the natural exponent. Lighthill argued that the errors in his formula can be as small as , and Childress showed that the errors are no worse than Childress1981 ().
It is convenient to parametrize the helix by the angle . The short-distance cutoff corresponds to a cutoff in the angle, defined by . In our experiments, the pitch angle is changed for a given helix by stretching the wire, which changes the pitch and radius , but keeps the contour length fixed. As before we denote the contour length of one helical pitch by . Since , the cutoff angle is approximately . This expression is accurate to three significant figures over the range of that we measure.
Just as in our resistive-force theory calculation, the force per unit length at a given point has only a component. Since our experiments show that the zero-force swimming speed is independent of immersed length once the immersed length is greater than one or two wavelengths, we may consider an infinite helix for which the magnitude of is uniform. Thus, we find
where . This integral is readily evaluated numerically. The predictions of Lighthill’s formula for are shown in Fig. 5 for two different values of .
v.3.2 Johnson’s slender-body theory
Johnson gave a more rigorous derivation of slender-body theory JOHNSON:1980p98 (), building on ideas of Keller and Rubinow keller_rubinow1976 (). By assuming that the filament has tapered ends, with a radius , Johnson derived slender-body theory formulas with accuracy. The velocity of a point on the filament is broken into local and nonlocal terms, , where
with is the local tangent vector, and
For the helix, , where
The expressions and are readily evaluated numerically. The predictions of Johnson’s slender-body theory for force-free swimming helices are also plotted in Fig. 5.
The simulation results from our modified boundary-element method agrees best with the experimental measurements. Although the resistive-force theory with captures the qualitative trend of the dependence of on , it is not very accurate even for filaments with [Fig 5]. The results are much improved when taking into account the finite , by using Gray and Hancock’s or Lighthill’s drag coefficients. The agreement with Lighthill’s version is noticeably better at small pitch angles for the two groups of filament thicknesses we test. However, these resistive force theories give qualitatively incorrect behavior at larger pitch angles (e.g., ), since the distance between adjacent pitches becomes shorter. The inaccuracy arises because resistive-force theory does not properly account for the hydrodynamic interactions between different parts of the helix Rodenborn:2013 ().
For the boundary-element method, the number of helix pitches being simulated is large enough, e.g., , so that the free-swimming speed is already independent on the detailed selection of . Note also that there are walls in the experiment, but the simulations and theories do not account for the wall; apparently the wall is far enough from the helix to have no effect on the speed. Despite the different error estimates, the results from both slender body theories (Lighthill’s and Johnson’s) are virtually identical. Note that the slope of the vs. [Fig. 5] curve is different for the asymptotic resistive force theory and the slender-body theories (or boundary-element technique) near . The difference arises because the asymptotic resistive force theory does not account for the thickness of the filament; sufficiently close to , a helix of nonzero thickness and many pitches long will intersect itself. For this reason, and because this regime is not physically relevant, we do not carry out our slender-body calculations and boundary-element simulations very close to .
Vi Other applications
Our technique is useful not only for infinite rigid helical filaments with circular cross section. It can be applied to many other geometries, such as helices with non-circular cross section, helices confined to circular tubes, and non-rigid bodies with helical symmetry. With minor modification, we can even apply our technique to situations that do not exhibit perfect helical symmetry, such as finite-length helices or bodies with non-uniform geometry. Two examples of these extensions are presented here.
vi.1 Confined geometry and non-rigid body
First we consider swimming in confined geometry. Noting that a long cylinder is a special case of a helix with pitch angle , we apply helical symmetry to the study of a helical swimmer in a co-axial cylindrical tube [inset of Fig. 6(a)]. The force densities distributed on the helical filament and the cylindrical wall are mapped to those on two circumferences along the above two surfaces, respectively. The dimensional reduction inherent in our technique allows us to study situations when the confining wall is very close to the filament surface without significantly increasing the number of grid points. For example, when the radius, , of the confining tube is comparable to the minimum tube radius, , the number of grid points required is comparable to that required for an unconfined helix and much smaller than the number of grid points required to accurately simulate a helix in a tight-fitting tube using a conventional boundary-element technique. Figure 6(a) shows that, at a given rotation rate , the free-swimming speed almost increases monotonically with decreasing .
Besides rigid bodies, our technique can also be extended to deforming structures by keeping the double-layer potential in Eq. 1. For instance, a transverse helical wave can propagate along the body and lead to motility. In rigid-body motion, all points of the helical filament rotate about the axis of the filament. For a helical wave, the filament deforms at every point, and the cross-sections of the filament do not rotate about the filament centerline as the deformation progresses. However, the centerline of a filament carrying a helical traveling wave of frequency rotates about the helix axis with rotation frequency . An example of swimming motility by such a helical wave is shown in Fig. 6(a), where we see that vs. is almost identical to the rigid-rotation case. This is reasonable since the filament here is still extremely thin, i.e., .
vi.2 Non-uniform geometry
As shown previously, we map the force densities on the entire helical surface to those on a single circumference either around the filament or around the confining structure. This mapping is no longer valid for short, finite-length filaments, where end effects are important, e.g., when . Nevertheless, the method can still be of value, and in such situations, instead of mapping the force densities to a single circumferences, it is natural to map them to a few circumferences spaced along the length of the filament. For this case, we introduce interpolation operators, and the orthogonal mapping in Eq. (7) becomes
where is the number of circumferences onto which the force densities are mapped, and is the weight function due to interpolation. The size of the linear equations to solve now becomes . When , the formula above (Eq. (37)) reduces to the case with full helical symmetry (Eq. (7)). Our preliminary studies show that even if is still much less than (e.g., in Fig. 6), we obtain reliable results due to the approximate helical symmetry. The computational cost is thus much reduced. Figure 6(b) shows an example of such an application. We study the length-dependence of the motility of a helical filament , and compare the results with that obtained from the experiments as described in our previous work Liu:PNAS2011 (). Here, the filament has a finite length and is modeled as an elongated spheroid shape with the radius of its cross-section . The aspect ratio is selected such that the mean radius satisfies . In both the experiments and the numerical simulation with , the free-swimming speed saturates at . We also validate here the use of helical symmetry for finite-length helices by letting . The error due to the enforced helical symmetry oscillates but eventually vanishes at . Detailed analysis and further applications of this extension to non-uniform body shape will be reported in a separate work.
We have shown that by exploiting helical symmetry, we are able to reduce a two-dimensional boundary-integral method to a one-dimensional method, without throwing away the detailed small-scale structure, such as the finite thickness of the flagellum and its distance to the surrounding structures. Our strategy is ideal for situations with symmetry that persists in the entire fluid flow and structure such as the flow generated by an infinitely long rotating helical filament. Nevertheless, in cases such as our experiment, with , the results for infinitely long systems with perfect helical symmetry can be applied to finite-length systems. Meanwhile, as demonstrated above, this boundary-integral method with reduced dimension can be applied to studying other more complicated fluid effects and can be extended to non-uniform structures. More results with application to these systems will be reported in future.
Acknowledgements.This work was supported by National Science Foundation Grant No. CBET-0854108.
Appendix A Integral kernel in the body-centerline coordinates
In Cartesian coordinates, a point on a surface with helical symmetry can be written as
where and are the radius and the phase angle of the helix, , are the radius and the polar angle of the surface intercepted by a plane (see Fig. 2). According to the definition of the Stokeslet (Eq. 2) and using the convention that , the integral kernel for modified Stokeslets (expressed in Eq. (23)) can be expressed as
Here, is a rotation operator defined in the body-centerline coordinates :
and vectors , are projections of vector () in the body-centerline coordinates. More specifically,
Appendix B Singularity reduction
In order to avoid the numerical divergence of the singular kernel , we evaluate the integral separately for the modified Stokeslet when is below grid size. As shown in Fig. 7, a small patch of surface (defined by the phase angel and polar angle ) is cropped out. The integral over the rest of the surface is performed numerically, using the trapezoidal rule. We evaluate the integral over that small patch in an analytical form, to the order of accuracy that is no worse than .
To formulate an infinitesimal expansion of such integral, we consider an arbitrary point on the patch , with a small displacement from its center . The displacement can be written as
where the second term on the right-hand side is due to surface curvature, and can be expanded in infinitesimal form as
where tensor is defined as
To the lowest order of infinitesimal expansion, vectors and in Eq. (39) become
with and , the rescaled distance (or ), can be expressed quadratically as
This quadratic form facilitates the analytic integration of the singular kernel,