A single-layer based numerical method for the slender body boundary value problem
William H. Mitchell, Henry G. Bell, Yoichiro Mori, Laurel Ohm, Daniel Spirn
AA SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARYVALUE PROBLEM.
WILLIAM H. MITCHELL, HENRY G. BELL, YOICHIRO MORI, LAUREL OHM, AND DANIEL SPIRNA bstract . Fluid flows containing dilute or dense suspensions of thin fibers are widespread in biological and industrialprocesses. To describe the motion of a thin immersed fiber, or to describe the forces acting on it, it is convenient to workwith one-dimensional fiber centerlines and force densities rather than two-dimensional surfaces and surface tractions.Slender body theories o ff er ways to model and simulate the motion of immersed fibers using only one-dimensional data.However, standard formulations can break down when the fiber surface comes close to intersecting itself or other fibers.In this paper we introduce a numerical method for a recently derived three-dimensional slender body boundary valueproblem that can be stated entirely in terms of a one-dimensional distribution of forces on the centerline. The method isbased on a new completed single-layer potential formulation of fluid velocity which circumvents some of the traditionalconditioning issues associated with the unmodified single layer potential. We give numerical results demonstrating thegood conditioning and improved performance of the method in the presence of near-intersections.
1. I ntroduction
The use of small parameters to simplify di ffi cult modeling and simulation problems is one of the outstandingsuccesses of classical applied mathematics. In fluid mechanics, one important example of a small parameter isthe aspect ratio of an immersed structure. Biological examples of slender immersed objects include microtubulesinside cells [32, 25] and cilia and flagella external to cells [30, 17, 6, 19]; many industrial processes such aspapermaking also rely on the properties of sparse or dense fiber suspensions [11, 27]. It is very appealing to workwith the one-dimensional centerlines of these thin structures rather than with their two-dimensional surfaces orthree-dimensional volumes. In addition to the computational e ffi ciency that comes with lowering the dimension ofthe problem, there are also important theoretical advantages. For example, it is much simpler to formulate a modelfor the centerline density of forces on a fiber than it is to model fully two-dimensional surface tractions.An attempt to make physical sense of forces and velocities defined on the one-dimensional centerline instead ofthe two-dimensional surface is known as a slender body theory . The first generation of methods in this categorywas called resistive force theory [12]. These methods were based on treating the fiber as a succession of pro-late spheroids while ignoring nonlocal hydrodynamic interactions; that is, according to resistive force theory, thecenterline velocity at a given location depends on the local force applied there, but not on the forces applied else-where on the fiber. Subsequent improvements accounted for nonlocal hydrodynamic e ff ects [8, 2, 19]; a prominentexample is the Keller-Rubinow formulation [15], reformulated for a periodic fiber [33, 7]:(1) 8 πµ u KRC ( s ) = (cid:20)(cid:16) I − ee T (cid:17) − (cid:18) π(cid:15) (cid:19) (cid:16) I + ee T (cid:17)(cid:21) f ( s ) + (cid:90) (cid:40)(cid:32) I | r | + rr T | r | (cid:33) f ( t ) − I + e ( s ) e ( s ) T | sin( π ( s − t )) | /π f ( s ) (cid:41) dt . This equation assumes a unit-length closed fiber with arclength parameterization γ ( s ); from this we define r = γ ( s ) − γ ( t ) and e = γ (cid:48) ( s ). The fiber radius is (cid:15) , and µ is the fluid viscosity. The first term on the right-hand sideis a local contribution to the centerline velocity from the force f imposed there, while the integral term representsthe nonlocal hydrodynamic interactions with the rest of the fiber surface. This formulation, along with refinementsthat address the case where the fiber has free ends, has been widely used [13, 34, 25, 20, 17, 24, 10, 18]. However,the resistive force theory and its nonlocal successors all assume that the fiber does not closely approach itself orother fibers. In Keller and Rubinow’s derivation [15] this assumption is expressed in their choice of an expressionfor the velocity in the inner region using the method of matched asymptotic expansions. More generally, it has notbeen clear exactly which fully three-dimensional problem any of the prior slender body theories are approximating,that is, they are not derived from a boundary value problem (BVP) with explicit boundary conditions at the fibersurface. Date : February 4, 2021. a r X i v : . [ m a t h . NA ] F e b WILLIAM H. MITCHELL, HENRY G. BELL, YOICHIRO MORI, LAUREL OHM, AND DANIEL SPIRN
More recently, some of us presented a well-posed three-dimensional BVP [23] given only the one-dimensionalcenterline force density, which we summarize as follows. We assume that the length scales for the flow problemare small enough that the Stokes model is appropriate. Writing u for velocity, p for pressure, and µ for viscosity,the PDE for the fluid domain is:(2) = − ∇ p + µ ∇ u , = ∇ · u . This PDE must be augmented by boundary conditions. We assume that the fluid is infinite in extent and quiescent,that is, u ( x ) → as | x | → ∞ . Then the only boundary of the fluid domain is the surface of a single immersed fiberwhose centerline is a closed loop parameterized by a C function from the circle S into R . We assume a smallconstant fiber radius (cid:15) > X ( s , θ ) = γ ( s ) + (cid:15) cos( θ ) N ( s ) + (cid:15) sin θ B ( s ) . We call s the centerline coordinate and θ the circumferential coordinate . The Bishop frame vectors N and B have unit length and are perpendicular to each other and to the tangent vector γ (cid:48) ( s ); see Appendix A for anexplanation of why we use this frame instead of the simpler Frenet-Serret frame as well as the details of ournumerical implementation. We can now state the boundary conditions for the PDE: u ( X ( s , θ )) = c ( s ) , (4) (cid:90) π σ ( s , θ ) · ν ( s , θ ) J ( s , θ ) d θ = f ( s ) | γ (cid:48) ( s ) | . (5)The first equation (4) restricts the boundary values of u by requiring that the fluid velocity be constant on eachcross section of the fiber, a fiber integrity constraint that preserves the circular shape of the cross sections. Thecenterline velocity function c ( s ) is unknown and must be determined as part of the solution. The second condition(5), on the surface derivatives of u , states that the integral of surface traction around each cross section must equalthe imposed centerline force density f ( s ), which is given as part of the problem. On the left side of (5), the stresstensor σ is defined by the limiting value of σ i j = − p δ i j + µ ( ∂ u i /∂ x j + ∂ u j /∂ x i ) as we approach the boundary fromwithin the fluid, and the factor J ( s , θ ) is the surface Jacobian, J ( s , θ ) = | X s × X θ | . Note that integration of eitherside of (5) from s = s = π yields the net force exerted on the fiber by the fluid; the factor | γ (cid:48) ( s ) | appears onthe right side to account for the possibility of using a nonconstant parameterization speed.The boundary value problem defined by (2)-(4)-(5), first presented in [23], presents some interesting challengesfrom a numerical perspective. This paper presents a new computational procedure appropriate for this new ver-sion of slender body theory. The new formulation is based on a representation of the velocity as the sum of asingle-layer potential on the fiber surface and a distribution of point source singularities on its centerline; we callthis the completed single-layer potential . This formulation allows for convenient traction formulas while also cir-cumventing the traditional disadvantages of the unmodified single-layer potential, as we discuss in Section 2. Wethen present our discretization of the problem, wherein the unknowns are a finite number of Fourier coe ffi cientsof the single layer density, together with evenly spaced values of the centerline velocity c ( s ). The details of thisprocedure are covered in Section 3, while the related matter of quadrature for singular integrals on the surface of athin fiber appears in Appendix B. We study the convergence rate of our numerical method in Section 4, includinga comparison to exact solutions for a uniformly translating rigid torus. These ‘exact’ solutions were originallyreported only to four-digit accuracy and so we also recomputed them to 12-digit or higher accuracy followingthe method of [1] from 1982. Finally we compare our numerical method to the widely used and computationallye ffi cient slender-body theory of Keller and Rubinow. The two methods agree to approximately order O ( (cid:15) . ) whenthe centerline does not approach itself, a somewhat greater convergence rate than was rigorously proved in [23],which demonstrated agreement at order O (cid:16) (cid:15) | log (cid:15) | / (cid:17) . This O ( (cid:15) . ) agreement is close to the O ( (cid:15) | log (cid:15) | ) errorpredicted by the matched asymptotics of [10, 15, 13, 34]. We also present tests showing that the Keller-Rubinowslender body theory breaks down when the fiber surface comes close to self-intersection (Section 5). Thus, ourmethod may be more suitable for simulating densely packed fiber suspensions. We conclude with priorities forfuture extensions of the method, such as treatments of multiple fibers, free ends, dynamic problems, and inertialflows. SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 3
2. T he completed single - layer potential Two attractive and widespread representations of unbounded Stokes flows are the single- and double- layerpotentials: u i ( x ) = π (cid:90) D (cid:32) δ i j r + r i r j r (cid:33) ρ j ( y ) dS y , (6) u i ( x ) = − π (cid:90) D r i r j r k r ν k ( y ) µ j ( y ) dS y . (7)Here x is an observation point in the fluid domain, D is the surface of an immersed particle, r = x − y , the surfacenormal vector ν points out of the particle into the fluid, and ρ and µ are the single- and double-layer densities,respectively. We use Einstein’s implicit summation notation (sum over the repeated indices j and k ). The tensor δ i j / r + r i r j / r in (6) is the point force or Stokeslet on R so the single-layer potential is a sum of point forceslocated on the fiber surface. Similarly, − r i r j r k / r is the stresslet and the double-layer potential (7) is a sum ofpoint stresses located on the fiber surface. An arbitrary Stokes flow exterior to D and decaying at infinity can berepresented by a sum of these potentials; in this case the densities ρ and µ have physical meaning (the surfacetraction and surface velocity, respectively). When only one or the other of these potentials appears, there aresome flows which cannot be uniquely represented. For example, the single-layer potential is unable to representflows with a nonzero volume flux across D , and if the density ρ is taken as a multiple of the normal vector ν the resulting single-layer velocity field is zero. In the language of linear algebra, the single-layer potential as anoperator mapping ρ to u has a one-dimensional nullspace, and its range has codimension one within the spaceof all quiescent Stokes flows evaluated at the fiber surface. The double-layer potential is also deficient: it cannotrepresent any flow which exerts a net force or torque on a closed particle, and the velocity in the fluid domain is zeroif the density function µ matches a rigid-body motion. The double-layer potential therefore has a six-dimensionalnullspace and its range has codimension six. Thus the single-layer and the double-layer potential, used in isolation,are each inadequate to represent general particulate flows.Power and Miranda provided a completion of the double-layer potential in the form of a point force and apoint torque inside of an immersed solid object [28]. We sketch this formulation on the left side of Fig. 1. Thestrengths of these internal singularities can be defined as six independent integrals of the double-layer density µ .This modified flow representation operator is now of full rank: any surface velocity field, including one whichexerts a force and torque on the enclosed volume, can be uniquely represented. This procedure leads to a widelyused Fredholm integral equation of the second kind for solving Dirichlet problems. The completed double-layerpotential does have a disadvantage despite its wide use: if one needs the local surface traction (rather than the netforce), it generally requires evaluating hypersingular integrals of the density, although there are remedies for thisproblem in the case of rigid particle motions [14, 5, 22].The single-layer potential also has a long computational history. The first numerical implementation of a bound-ary integral method was based on an unmodified single-layer representation; see [36]. The single-layer represen-tation has the advantage that the local surface tractions are comprised of convergent integrals; however, thereare two problems with the unmodified single-layer formulation. The first is the existence of a nullspace (a zeroeigenvalue) discussed above. The second is that the operator carrying the density ρ to the surface velocity u hasarbitrarily small nonzero eigenvalues. In a numerical implementation, any surface integration procedure is sub-ject to numerical error, which can lead to invertible discrete linear operators even when the continuous operatoris singular. The condition number of the discrete operator increases as the integration procedure becomes moreaccurate, but the method can give acceptable results in an intermediate range where the discretization is neither toocoarse nor too fine. For problems where the surface velocity is prescribed, the completed double-layer methodslead to second-kind integral equations, which do not su ff er from this conditioning issue.In the more specific setting where the particle is a slender fiber, Koens and Lauga recently obtained versionsof slender-body theory by starting with either the single- or double-layer potentials and using matched asymptoticexpansions in the fiber radius [16]. They found that the use of a single-layer potential leads to a singular system ofequations; in particular, the first Fourier mode of the force density in the circumferential direction is not uniquelydetermined from the surface velocity, while no solution exists if the first mode of the surface velocity is nonzero.As they remark, this failure corresponds to the inability of the single-layer potential to represent volume changes. WILLIAM H. MITCHELL, HENRY G. BELL, YOICHIRO MORI, LAUREL OHM, AND DANIEL SPIRN F igure
1. The completed single- and double-layer potential representations of Stokes flow areconstructed by superposing fundamental solutions with singularities in the interior or distributedover the surface of an immersed particle. In Power and Miranda’s completion of the double layerpotential (at left), the flow consists of a distribution of stresslets on the surface together with apoint force and a point torque in the interior (shaded). Our proposed completion of the single-layer potential appears on the right; the flow consists of a distribution of point forces on thesurface together with one or more point sources in the interior (shaded). The total strength of theinternal point sources is proportional to the inner product of the surface normal with the single-layer density function.We now o ff er a completion procedure for single-layer potentials which resolves the nullspace issue and mayhave some advantages in comparison to the double-layer formulations. The procedure is very simple: in additionto the single-layer potential, we include the flow due to one or more point sources located internal to the particlewhose total strength is proportional to the inner product of the single-layer density with the surface normal vector.More concretely, let D be a closed particle surface, and let χ be a map carrying points on the surface to points inthe interior. For a convex particle, χ could be a constant map carrying surface points to the particle centroid. Inour setting with a closed slender fiber, we let χ carry surface points to the corresponding centerline point. We thendefine the velocity on the exterior to the particle by the modified single-layer equation(8) u i ( x ) = π (cid:90) D (cid:32) δ i j | x − y | + ( x i − y i )( x j − y j ) | x − y | (cid:33) ρ j ( y ) dS y + π (cid:90) D ν j ( y ) ρ j ( y )( x i − χ i ( y )) | x − χ ( y ) | dS y . The first term in this equation is the single-layer potential with density ρ . The second term is the sum of thevelocities at x due to point sources distributed at χ ( y ), each with strength ν ( y ) · ρ ( y ), where ν is the normalvector pointing out of the particle (into the fluid domain). We note that the first integrand becomes singular if x ison the particle surface, but the second integral remains regular even in that case. We are using the Green’s functionsfor unbounded flow, although the idea should also work for bounded flow domains. Although this procedure isexactly analogous to the widely used completion procedure for the double layer formulations, we have not foundany discussion of it in the literature.The main advantage of the completed single-layer representation is that the surface tractions are easy to compute.The single-layer potential has a known stress field, and so does the point source. At observation points x in thebulk fluid the stress tensor is σ ik = − π (cid:90) D ( x i − y i )( x j − y j )( x k − y k ) | x − y | ρ j ( y ) dS y + π (cid:90) D (cid:32) δ ik | x − χ ( y ) | − x i − χ i ( y ))( x k − χ k ( y )) | x − χ ( y ) | (cid:33) ν j ( y ) ρ j ( y ) dS y , (9)and the surface traction is the contraction of the stress with the surface normal in the limit as x approaches thefiber surface from the fluid domain, that is, from the side into which ν points. For the point source the integrands SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 5 F igure
2. We represent the fluid velocity using a distribution of point forces on the fiber surface(orange arrows) together with a distribution of point sources on the the centerline (green arrows).The strength of the point sources are taken equal to the fluxes of the single layer density througheach cross section, so that the volume creation rate at γ ( s ) is (cid:82) π ν ( s , θ ) · ρ ( s , θ ) J ( s , θ ) d θ . In thefigure, the orange arrows have a positive flux out of the cross section and accordingly the greenarrows indicate a point source rather than a sink.remain regular and for the single-layer part we use results from Pozrikidis [29] to arrive at t i ( x ) = − ρ i ( x ) − π (cid:90) D ( x i − y i )( x j − y j )( x k − y k ) | x − y | ν k ( x ) ρ j ( y ) dS + π (cid:90) D (cid:32) δ ik | x − χ ( y ) | − x i − χ i ( y ))( x k − χ k ( y )) | x − χ ( y ) | (cid:33) ρ j ( y ) ν k ( x ) ν j ( y ) dS y . (10)When χ is a constant function, we can prove that the velocity representation (8) is unique and can represent anarbitrary flow. We expect but do not prove that it also holds for nonconstant χ , which is a more computationallyappropriate choice for the slender fiber geometry. In the constant case the velocity can be more simply written as(11) u i ( x ) = π (cid:90) D (cid:32) δ i j | x − y | + ( x i − y i )( x j − y j ) | x − y | (cid:33) f j ( y ) dS y + π ( x i − χ i ( y )) | x − χ ( y ) | (cid:90) D ν j ( y ) f j ( y ) dS y so that the rate of volume creation at χ ( y ) is precisely the inner product of ν and f . Thus, to represent an arbitraryflow, one first determines the volume flux rate α and then takes an initial surface distribution f = α | D | ν . Thenthe velocity induced by f has the desired volume flux. Now the di ff erence between this flow and the desired oneis flux-free and can be represented in infinitely many ways by a single-layer potential, but only in one way by asingle-layer potential with density f satisfying (cid:82) D f · ν =
0. Then f = f + f is a density function inducing thedesired velocity. Uniqueness is also a consequence of the fact that a flow with zero flux can be uniquely representedby a single-layer potential whose density has zero inner product with the surface normal.We find that our discrete condition numbers do not increase as our quadrature is refined, although they appearto increase linearly with wavenumber as we resolve higher Fourier modes in the solution; see Table 1. Therefore,the completed single-layer potential is an e ff ective way to formulate our numerical method, as described below.In the remainder of this paper we use this completed single-layer velocity formulation to develop a computa-tional method suitable for simulating closed slender fibers. In a resonance with the finding by Koens and Laugathat the slender body theory based on a single-layer potential has a deficiency at mode one [16], we find that ourcorrection procedure modifies only the terms corresponding to modes zero, one, and two in the discrete version,and the greatest modification is to the first mode.3. D iscretization of the slender body BVPAs stated previously, we consider a single fiber in a quiescent fluid without boundary. The fiber centerlineis a closed loop. The fluid velocity is represented by the sum of a single-layer potential and a distribution ofpoint sources along the centerline; see Fig. 2. Let γ : [0 , π ) → R be a parameterization of the centerline, notnecessarily of constant speed. Let X ( s , θ ) parameterize the fiber surface. The circular cross sections are normal tothe centerline γ and have uniform radius (cid:15) . Let ν ( s , θ ) denote the unit surface normal vector pointing out of thefiber (into the fluid domain), and let J ( s , θ ) denote the Jacobian, J = | X s × X θ | . Then the completed single-layer WILLIAM H. MITCHELL, HENRY G. BELL, YOICHIRO MORI, LAUREL OHM, AND DANIEL SPIRN fluid velocity equation (8) reduces to:(12) u i ( x ) = π (cid:90) D (cid:32) δ i j r + r i r j r (cid:33) ρ j ( y ) dS y + π (cid:90) π R i R (cid:90) π ρ j ( s , θ ) ν j ( s , θ ) J ( s , θ ) d θ ds , where the vector r = x − y points from the surface integration point to the observation point, and the vector R = x − γ ( s ) points from the centerline integration point to the observation point. The first term is a single-layerpotential with density ρ , and the second term is a distribution of point sources over the centerline.Similarly, the equation (10) for the surface traction exerted by the exterior fluid on the fiber at a surface point x = X ( s ∗ , θ ∗ ) ∈ D becomes t i ( x ) = − ρ i ( x ) − π (cid:90) π (cid:90) π r i r j r k r ρ j ( y ) ν k ( x ) J ( s , θ ) d θ ds + π (cid:90) π (cid:32) ν i ( s ∗ , θ ∗ ) R − R i R j ν j ( s ∗ , θ ∗ ) R (cid:33) (cid:90) π ρ k ( s , θ ) ν k ( s , θ ) J ( s , θ ) d θ ds (13)with r and R now given by r = X ( s ∗ , θ ∗ ) − X ( s , θ ) and R = X ( s ∗ , θ ∗ ) − γ ( s ).We can now substitute these velocity and traction expressions into the boundary conditions (4)-(5). Let f ( s ) bea vector function of the centerline describing the force density on a cross section. Then we have0 = − c i ( s ∗ ) + π (cid:90) π (cid:90) π (cid:32) δ i j r + r i r j r (cid:33) ρ j ( s , θ ) J ( s , θ ) ds d θ + π (cid:90) π R i R (cid:90) π ρ j ( s , θ ) ν j ( s , θ ) J ( s , θ ) d θ ds (14) f i ( s ∗ ) (cid:12)(cid:12)(cid:12) γ (cid:48) ( s ∗ ) (cid:12)(cid:12)(cid:12) = (cid:90) π (cid:34) − ρ i ( s ∗ , θ ∗ ) − π (cid:90) π (cid:90) π r i r j r k r ρ j ( s , θ ) ν k ( s ∗ , θ ∗ ) J ( s , θ ) ds d θ + π (cid:90) π (cid:32) ν i ( s ∗ , θ ∗ ) R − R i R j ν j ( s ∗ , θ ∗ ) R (cid:33) (cid:90) π ρ k ( s , θ ) ν k ( s , θ ) J ( s , θ ) d θ ds (cid:35) J ( s ∗ , θ ∗ ) d θ ∗ , (15)where the unknowns are the velocity c ( s ), a function of the centerline only, and the single layer density ρ ( s , θ ). Thefirst equation says that the surface velocity is independent of the circumferential coordinate θ (the fiber integritycondition ). The second condition states that the integral of the surface traction around a circular cross section ofthe fiber matches f . The factor of | γ (cid:48) ( s ∗ ) | on the left of (15) accounts for the speed of the parameterization sothat upon integration from 0 to 2 π in s ∗ , we would find that the net force (right side) is equal to the integral of thecenterline force with respect to arclength (left side).To discretize this system of integral equations, we start by enforcing them only at finitely many points. For (14)we let ( s ∗ , θ ∗ ) range over a regular grid. Letting n s and n θ be odd integers giving the number of grid points in eachdirection, we take ( s ∗ , θ ∗ ) = (2 π j s / n s , π j θ / n θ ) for 0 ≤ j s < n s and 0 ≤ j θ < n θ . Similarly, we enforce (15) onlyfor s ∗ = π j s / n s with 0 ≤ j s < n s . This gives a total of 3( n s n θ + n s ) scalar equations. To obtain finitely manyunknowns, we seek the density ρ in a finite-dimensional space of complex exponentials: ρ (cid:96) ( s , θ ) = (cid:88) k s , k θ α (cid:96), k s , k θ exp (cid:16) √− k s s + k θ θ ) (cid:17) . Here the indices range over − ( n s − / ≤ k s ≤ ( n s − / − ( n θ − / ≤ k θ ≤ ( n θ − /
2. Consideringthe three possible values of the space dimension index (cid:96) , we have a total of 3 n s n θ unknown Fourier coe ffi cientsof ρ . The values of c i ( s ∗ ) at s ∗ = π j s / n s provide the remaining 3 n s unknowns, which we abbreviate by writing c i , j s = c i (cid:16) π j s n s (cid:17) for 0 ≤ j s < n s . Upon substituting this expression for ρ and moving the sums outside the integrals, SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 7 n θ = n θ = n θ = (cid:15) = − n s = . · . · . · n s =
21 1 . · . · . · n s =
63 1 . · . · . · n s =
189 1 . · . · . · (cid:15) = − n s = . · . · . · n s =
21 1 . · . · . · n s =
63 4 . · . · . · n s =
189 2 . · . · . · (cid:15) = − n s = . · . · . · n s =
21 1 . · . · . · n s =
63 3 . · . · . · n s =
189 1 . · . · . · T able
1. The condition numbers of the discrete linear systems increase as the fiber radius (cid:15) de-creases, and they do not increase as the quadrature is refined. The condition numbers generallyincrease with the centerline and circumferential discretization parameters n s and n θ , but in a com-plicated way which also depends on (cid:15) . The discrete systems reported in this table were generatedfor a trefoil knot centerline using a well resolved quadrature ( q n =
40 for all reported values; seeAppendix B for details). Further tests with a more refined quadrature q n =
50 (not reported here)gave relative changes of less than 1 / = − c i , j s + π (cid:88) j , k s , k θ α j , k s , k θ (cid:90) π (cid:90) π (cid:32) δ i j r + r i r j r + R i ν j ( s , θ ) R (cid:33) exp (cid:16) √− k s s + k θ θ ) (cid:17) J ( s , θ ) d θ ds (16) F i , j s = π n θ n θ − (cid:88) j θ = J ( s ∗ , θ ∗ ) (cid:40) (cid:88) j , k s , k θ α j , k s , k θ (cid:34) − δ i j exp (cid:16) √− (cid:0) k s s ∗ + k θ θ ∗ (cid:1)(cid:17) + π (cid:90) π (cid:90) π (cid:32) − r i r j ( r · ν ( s ∗ , θ ∗ )) r + ν i ( s ∗ , θ ∗ ) R ν j ( s , θ ) − R i ν j ( s , θ ) R · ν ( s ∗ , θ ∗ ) R (cid:33) exp (cid:16) √− k s s + k θ θ ) (cid:17) J ( s , θ ) ds d θ (cid:35)(cid:41) . (17)In these equations we are using the abbreviations s ∗ = π j s / n s , θ ∗ = π j θ / n θ , F i , j s = f i ( s ∗ ) | γ (cid:48) ( s ∗ ) | , r = X ( s ∗ , θ ∗ ) − X ( s , θ ), and R = X ( s ∗ , θ ∗ ) − γ ( s ). The discrete fiber integrity equation (16) holds for all 0 ≤ j s < n s and0 ≤ j θ < n θ , while the discrete averaged force equation (17) holds for 0 ≤ j s ≤ n s . Note that to arrive at (17) wehave replaced the outermost integral (in θ ∗ ) from (15) with a trapezoidal rule sum.To set up the linear algebra system, we have to evaluate the integrals in (16)-(17). The Stokeslet integrand in(16) and the stresslet integrand in (17) both have a 1 / r singularity as ( s , θ ) → ( s ∗ , θ ∗ ), so the numerical quadratureprocedure is a nontrivial problem. We give details of our method in Appendix B. An interesting feature of thisformulation is that the accuracy of the quadrature can be chosen independently of the matrix size. Once theintegrals have been computed, we have a dense and non-normal system of linear equations. The condition numbersfor the problems we considered range from approximately 10 to 10 . We chose to use the SVD for the linear solvebecause of its good performance with poorly conditioned systems; the computational expense of this method forour dense, non-normal system of linear equations is acceptable because the overall solution time is dominated bythe matrix assembly rather than the linear solve. WILLIAM H. MITCHELL, HENRY G. BELL, YOICHIRO MORI, LAUREL OHM, AND DANIEL SPIRN
Circumferential integrals of nonsingular terms.
The variable R = | X ( s ∗ , θ ∗ ) − γ ( s ) | has no dependence onthe circumferential integration variable θ ; moreover, its minimum as a function of s is (cid:15) rather than zero. Thereforethe terms with R in the denominators are more analytically tractable than those involving r = | X ( s ∗ , θ ∗ ) − X ( s , θ ) | .In particular, they can always be reduced to one-dimensional integrals and for many parameter values they simplyvanish. To see this, define the circumferential integrals(18) M ( k θ , s ) = (cid:90) π ν j ( s , θ ) exp (cid:16) √− k s s + k θ θ ) (cid:17) J ( s , θ ) d θ and use this notation to rewrite the R integrals as18 π (cid:90) π (cid:90) π R i ν j ( s , θ ) R exp (cid:16) √− k s s + k θ θ ) (cid:17) J ( s , θ ) d θ ds = π (cid:90) π R i R M ( s , k θ ) ds (19) 14 π (cid:90) π (cid:90) π (cid:32) ν i ( s ∗ , θ ∗ ) R ν j ( s , θ ) − R i ν j ( s , θ ) R · ν ( s ∗ , θ ∗ ) R (cid:33) exp (cid:16) √− k s s + k θ θ ) (cid:17) J ( s , θ ) ds d θ = π (cid:90) π (cid:32) ν i ( s ∗ , θ ∗ ) R − R i R · ν ( s ∗ , θ ∗ ) R (cid:33) M ( s , k θ ) ds . (20)Now the inner integrals can be solved explicitly by expanding J ( s , θ ) = (cid:15) | γ (cid:48) ( s ) | (1 − (cid:15) cos θκ ( s ) − (cid:15) sin θκ ( s )) and ν ( s , θ ) = cos θ N ( s ) + sin θ B ( s ); this leads to(21) M ( k θ , s ) = − π(cid:15) | γ (cid:48) ( s ) | (cid:16) N j ( s ) κ + B j ( s ) κ (cid:17) e √− k s s k θ = π(cid:15) | γ (cid:48) ( s ) | (cid:16) N j ( s ) + √− B j ( s ) (cid:17) e √− k s s k θ = − π(cid:15) | γ (cid:48) ( s ) | (cid:16) N j ( s ) + √− B j ( s ) (cid:17)(cid:16) κ + √− κ (cid:17) e √− k s s k θ = k θ > . The values for negative k θ are the complex conjugates of these. Thus, we could use one-dimensional quadrature for | k θ | < | k θ | ≥ . The integrals which reduce so conveniently to one-dimensionalquadrature problems all correspond to the point sources we introduced to complement the single-layer potential.As we mentioned in the previous section, a recent study [16] shows that a slender-body formulation based on anunmodified single-layer potential will be noninvertible precisely at mode k θ =
1. It is encouraging to see thatour correction procedure modifies the first mode most (order (cid:15) ) while leaving the | k θ | > (cid:15) .The reduction to one-dimensional quadrature for | k θ | ≤ R terms when | k θ | >
2, in accordance with(21). 4. E rror analysis of the numerical method
Vertically translating torus.
In order to compare our method to an exact solution, we consider a torus whosecenterline is a unit circle and whose cross sections have radius (cid:15) , as above. If the forcing is uniform and alignedwith the rotational symmetry axis of the body, the resulting velocity is also uniform (it does not vary along thecenterline) and in the same direction. Because of the simple geometry and because the computed velocity isa rigid-body motion, we can compare to Amarakoon’s analytical solution [1] which comes from separation ofvariables in toroidal coordinates. As in that work, we investigate the behavior of a nondimensionalized force as (cid:15) varies. The quantity we evaluate is F (cid:48)∞ = F (cid:48)∞ ( (cid:15) ) = F πµ U (1 + (cid:15) ) . Here U is the velocity of the torus and F is the net force, both measured in the direction of the symmetry axis. Thedenominator is the force on a sphere whose outer diameter (2 + (cid:15) ) coincides with the outer diameter of the torussurface, so we expect F (cid:48)∞ < F (cid:48)∞ → (cid:15) →
0. To be clear, the analytical solution [1] is posed as a resistanceproblem (set U = F ) while our method is essentially a mobility problem (set F = U ), but the nondimensionalization allows us to compare the results. SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 9 /(cid:15) . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . · − . . . . . . · − . . . · − T able
2. Drag coe ffi cients for a torus of centerline radius 1 and cross-sectional radius (cid:15) translat-ing along its symmetry axis; reported values are the net force scaled on 6 πµ U (1 + (cid:15) ). Previouslyreported values were obtained in the 1980s and only listed five digits [1], so we have recomputedthem. The accuracy estimates on the right refer to errors in these ‘exact’ solutions, not compar-isons to our numerical method.Before discussing the convergence rate of our numerical method, we will digress to comment further on thereference solutions. Although we refer to the results published in the 1982 paper [1] as “exact,” they were onlyreported to four or five digit accuracy. This is because the exact solution procedure yields an infinite set of equationsrelating the coe ffi cients of the toroidal harmonic expansion of the solution, and this system of equations has to betruncated and then solved with a numerical linear algebra procedure. Thus even the “exact” solutions are subject tosome numerical uncertainty. To get more digits, we reimplemented the procedure described by [1] and used largerlinear systems and double- rather than single-precision arithmetic. To assess the numerical error in these “exact”solutions, we computed the di ff erence between the two sides of the identity (cid:80) n ≥ nB n = (cid:80) ∞ n ≥ C n , which appearson the first page of [1] and report this quantity as a proxy for the accuracy of the results. Our recomputations of theexact solutions are given in Table 2 together with these accuracy estimates. These results, obtained through a strictreimplementation of the 1982 presentation, are su ffi cient to get as many digits as we require for verification of ourown numerical method. However, we note for completeness that O’Neill and Bhatt gave an improved version [26]of Goldman, Cox and Brenner’s classic work on the problem of a sphere moving near a plane wall [9], removing theneed to solve a linear system of equations. It is likely that similar methods could improve Amarakoon’s formulation[1]. However, we found that the unmodified 1982 algorithm already gives su ffi ciently accurate results, and we didnot attempt to improve it.In Fig. 3 we compare the accuracy of our numerical scheme to these recomputed exact solutions. For fixedgrid parameters ( n s , n θ ), the convergence rate is spectral in the number of quadrature nodes used to evaluate theintegrals, but it eventually reaches a plateau where other sources of error dominate. One source of error is from thepresence of unresolved Fourier modes in the single-layer density ρ ; this is clearly the issue for the (cid:15) = − curve E rr o r v s . e x a c t d i m e n s i o n l e ss n e t f o r c e n = 5 quadrature nodes n = 9 = 10 = 10 = 10 = 10 = 10 n = 13 F igure
3. A comparison of our numerical results to exact solutions for a simple torus with uni-form forcing along the symmetry axis. The horizontal coordinate gives the number of quadraturenodes used to evaluate each matrix entry. The vertical axis shows the relative error versus theexact solutions from Table 2. The method appears to be spectrally accurate with respect to thequadrature, subject to a ceiling on the accuracy imposed by the condition number (approximately100 /(cid:15) ) and the grid resolution parameters n s and n θ . To obtain a single scheme we should have thequadrature rule and the parameters n s and n θ increase together, but the optimal way to do this willdepend on the fiber geometry and forcing. We note that a small n θ is acceptable when (cid:15) becomessmall.in the left panel of Fig. 3, because the issue is resolved by increasing the circumferential discretization parameter n θ . For smaller values of the radius (cid:15) , there is apparently no error reduction at all when n θ increases from 5 to13; the higher circumferential modes are, unsurprisingly, irrelevant, which allows us to save some computationalexpense in the more complicated simulations discussed later. When the fiber centerline is more complicated thana simple circle, the centerline discretization parameter n s will also play an important role. The condition numberof our linear systems is approximately 100 /(cid:15) , and so there is also a loss of accuracy as (cid:15) → Trefoil knot.
We now study the convergence rate of our numerical method in the presence of complicatedcenterline geometry and higher-frequency forcing. Specifically, we consider the case of a trefoil knot where thecenterline and the applied force density are given respectively by(22) γ ( s ) = sin s + s cos s − s − sin 3 s f ( s ) = sin ks + ks − cos ks + ks . In the case where k =
1, the applied forcing is linear in the space variables and the resulting fluid flow resemblesan extensional flow (although it decays rather than grows with distance from the fiber). This flow field is illustratedin Fig. 4. We chose this trefoil curve because it has a nontrivial three-dimensional structure but does not comeclose to self-intersection, a more challenging problem that we will consider later in the paper. We found that theconvergence rate appears to be spectral (concave down on a log-log plot) until we reach a plateau in accuracywhich depends on (cid:15) but not the other parameters in the problem.5. C omparison to other slender body theories
We now compare the results of our BVP-based numerical method to the predictions of Keller-Rubinow (KR)slender-body theory as stated above (1). This formulation assumes that the length of the fiber has been scaled to1 and the parameterization has with constant speed. While our general numerical method does not require a unit-length fiber or even a constant-speed parameterization, in this section we numerically reparameterize all curves andthen use the constant-speed, unit length parameterizations for both the KR slender-body theory and our method.The evaluation of the one-dimensional integrals in the KR version deserves some comment. The two terms in
SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 11 ( n s , k )
4. To study the convergence of our numerical method, we consider a trefoil centerline γ ( s ) = (sin s + s , cos s − s , − sin 3 s ) and we vary the forcing frequency k , the fiberradius (cid:15) , the discretization parameter n s , and the number of quadrature nodes. The flow field for (cid:15) = . k = f ( s ) = ( x ( s ) , − y ( s ) , ff erence is that a true extensional flow would grow with distance from the center, whereasour velocity field decays with distance from the fiber. The fiber surface is colored according to theagreement of the surface velocity with Keller-Rubinow slender-body theory, with red indicatingthe maximum error. The error of our numerical method decreases exponentially in the number ofquadrature nodes until it reaches a floor which depends on the fiber radius (cid:15) . The many similarcurves illustrate that the accuracy is not sensitive to the parameter n s (number of Fourier modes inthe centerline direction) or the forcing wavenumber parameter k . The errors displayed in this fig-ure were computed by a comparison to a reference numerical solution with quadrature parameter39 (22496 nodes); the seven quadrature parameters displayed in the figure are evenly spaced from15 to 33. To compare results that were computed on di ff erent periodic grids, we use Fourier inter-polation to sample the centerline velocities on an evenly spaced 500-point grid, then subtract, takethe 2-norm of each of the resulting 500 vectors, and finally take the maximum ( ∞ -norm). We scaleby dividing by the norm of the reference solution so that the results we report are relative. The bluecolor corresponds to (cid:15) = .
05, while yellow and green correspond to (cid:15) = .
005 and (cid:15) = . k is equal to 1, 2, and 3, and we used n s = , , , . Among all of these choices, the only one that seems to matter is (cid:15) (evidently n s =
41 is alreadysu ffi cient to resolve the trefoil geometry). For all of these trials, the condition number of the dis-crete system is approximately 100 /(cid:15) , that is, increasing as the fiber radius shrinks but otherwiseinsensitive to the discretization parameters.the integrand both diverge like 1 / | s − t | , but they cancel to give a bounded integrand with a jump discontinuity. Inprinciple there is still a risk of machine arithmetic error when s ≈ t due to the subtraction of nearly equal quantities,but we found no problems of this kind in our experiments. For these integrals we used Gauss-Legendre quadraturewith 200 nodes on each of the subintervals [ s , s + / s + / , s + / s + / , s + t = s , but we also needed good resolution near s + / s + / Circular centerline, in-plane low-frequency forcing.
We begin with a fiber whose centerline is circular.Instead of a uniform force density directed along the circle’s symmetry axis, we consider an in-plane forcing with fiber radius 10 r e l a t i v e d i ff e r e n c e S B T v s . P D E mode 0mode 1mode 2reference ( )reference ( log( )) F igure
5. For a circular centerline and in-plane forcing at various wavenumbers, the discrepencybetween Keller-Rubinow slender body theory and our BVP-based method decays with fiber radiusat rate O ( (cid:15) . ). The error increases modestly with the wavenumber of the imposed centerline forcedensity.some sinusoidal variation:(23) γ ( s ) = π cos s sin s f ( s ) = cos ms . The results for the first three modes m = , , (cid:15) ranging from 10 − to 10 − are given in Fig. 5. Theconvergence rate is somewhat slower than (cid:15) .5.2. Tests with four curves.
We now compare Keller-Rubinow slender-body theory with our BVP-based methodusing several noncircular centerlines. We consider four closed curves: a planar ellipse of aspect ratio 2.5, a trefoilknot, the planar boundary of the unit ball of in the 4-norm, that is, the curve defined by x + y =
1, and finally afigure-eight loop. The initial parameterizations are respectively(24) γ ( s ) = cos s . s , sin s + s cos s − s − sin 3 s , (cos s + sin s ) − / cos s (cos s + sin s ) − / sin s , sin 2 s . s . s , although we then reparameterize for constant speed and unit length as mentioned above. The results are shown inFig. 6. The rate at which the discrepancy decays is between O (cid:16) (cid:15) . (cid:17) and O (cid:16) (cid:15) . (cid:17) across the four centerlines andthe three norms (1, 2, or ∞ ). That is, our method behaves similarly to the Keller-Rubinow formulation in thesetests.5.3. Exploring the limitations of SBT: self-intersection.
Keller and Rubinow’s derivation is based on the methodof matched asymptotic expansions; in the inner region, near the fiber surface, they assume that the fluid velocityis well approximated by the flow near a translating and rotating cylinder. This assumption is invalid when anothersection of the fiber is located within a distance of O ( (cid:15) ). Accordingly, we expect the Keller-Rubinow to break downat least locally in the presence of near self-intersections of the fiber surface. SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 13 R e l a t i v e e rr o r L convergence TrefoilEllipse4-BallEight 10 L convergence L convergence F igure
6. The discrepancy between our BVP-based numerical method and the Keller-Rubinowslender body theory decreases at a rate between O ( (cid:15) ) and O ( (cid:15) ) (dashed reference curves) as (cid:15) → n s =
201 for the trefoil, ellipse, and4-ball and n s =
301 for the figure-eight centerline. In all cases we took n θ = n q =
35 (resulting in about 12,000 quadrature nodes). In reporting the relative errorswe first subtract the centerline velocities of the full numerical and KR slender-body procedures,then take the 2-norm over the space dimension, then take either the 1 − , 2 − , or ∞ -norm of theresulting vector, and finally divide by the KR slender-body norm obtained in the same way. Wecomputed convergence rates using linear regression with the smallest three (cid:15) -values and found thatthe convergence rates are comparable for all four curves and all three norms; the fastest rate was1.78 (trefoil knot & 1-norm) and the slowest was 1.67 (figure eight & ∞ -norm).To be more quantitative about what it means for the surface to nearly intersect itself, consider the quantity [3, 21](25) σ ( γ ) = min s , t ∈ [0 , (cid:107) γ ( s ) − γ ( t ) (cid:107) sin | π ( s − t ) | . Here we are assuming a constant-speed, unit length parameterization so that σ is a geometric quantity. For a circle, σ = . σ = . σ = . σ = . σ = . (cid:15) should result in good agreement between our methodand the Keller-Rubinow formulation, just as in Fig. 6. To create a test where the two methods are likely todiverge, we consider a family of centerlines where we can move the near-intersection points closer together whilesimultaneously reducing the fiber radius so that the ratio between the gap size and the radius remains constant.These centerlines are initially defined by(26) γ ( s ) = cos( s )(1 + H cos 3 s )sin( s )(1 + H cos 3 s ) H sin 3 s , but then numerically reparameterized and scaled so that its speed is constant and the total length is 1. Here H < H → H →
1, the curve points γ ( π/ γ ( π ) and γ (5 π/
3) all approach the origin;we refer to the distance between any two of these points, after rescaling, as the gap size . Thus any value of the H gap σ able
3. Geometry data for the centerlines with near-self intersections defined by (26) after repa-rameterization. As H approaches 1, the centerline points γ ( π/
3) and γ ( π ) and γ (5 π/
3) all approachthe origin. The distance between any two of these is the gap ; this approaches zero along with theself-intersection quantity σ defined by (25).fiber radius (cid:15) exceeding half of the gap size would result in a self-intersection of the fiber surface. The relationshipbetween H , the gap size, and the self-intersection quantity σ is given in Table 3.As an example, Fig. 7 shows the centerline that results from setting H = .
8. The arrows indicate the centerlineforce function we impose, given by(27) f ( s ) = − cos (cid:18) s + π √ sin( s ) (cid:19) (cid:18) s + π √ sin( s ) (cid:19) . We chose this peculiar form for the forcing function because of its values at the three points of nearest self-intersetion, illustrated as black arrows in Fig. 7: f ( π/ = ˆ z , f (5 π/ = − ˆ z , and f ( π ) = ˆ x . That is, two of thethree branches that pass near the origin are being forced in opposite directions tangential to the fiber surface whilethe third branch is being forced in the plane normal to the centerline. This complicated scenario is designed to findbreakdowns in the Keller-Rubinow slender-body theory.Figures 8 and 9 illustrate the breakdown of the Keller-Rubinow formulation with self-intersections. In one set ofsimulations we fix H = . (cid:15) →
0. As expected, our method agrees with the Keller-Rubinow formulationat approximately order O ( (cid:15) . ) (gray polygon in Fig. 9) and this convergence is uniform in the centerline coordinate s (top panel of Fig. 8). Indeed, the fact that the convergence rate is essentially the same in the 1-, 2- and ∞ -normssuggests that the discrepancy between the two methods must decay uniformly in s . This situation changes, however,when we allow the parameter H to increase toward 1 while keeping the fiber radius at a constant fraction of thegap size. We did this with ratios of 1 /
10 (second row in Fig. 8, blue polygon in Fig. 9), 1 / . (third row inFig. 8, orange polygon in Fig. 9), and 1 /
100 (fourth row in Fig. 8, green polygon in Fig. 9). In all of these caseswe see that the local error does not decrease uniformly. Near the intersection points at s ∈ { π/ , π, π/ } the twomethods appear to be converging to di ff erent answers and the discrepancy is O (1). At fiber points away from theintersection regions, the discrepancy does decay at about the same rate. Thus, the overall error stagnates whenmeasured in the ∞ -norm (the upper boundary of the polygons in Fig. 9) but continues to decrease in the 1-norm(the lower boundary of the polygons in Fig. 9). The breakdown in the Keller-Rubinow formulation follows fromthe invalidity of their inner expansion of velocity when another part of the fiber surface lies within a distance of O ( (cid:15) ). It is interesting that we can detect this breakdown even when the other fiber surface is relatively far awayfrom the observation point (100 (cid:15) in the last set of tests).6. C onclusion and future work In this paper we presented evidence that the standard slender body theories based on matched asymptotic expan-sions can give inaccurate results when the fiber surface approaches itself, a situation that is common in biophysicaland industrial processes. The breakdown is local in the sense that these formulations still give good results onisolated sections of the fiber. As an alternative, we gave a numerical method for a fully three-dimensional slenderbody Stokes boundary value problem which can be stated with reference to a one-dimensional centerline forcedensity only. The numerical method is based on a completed single-layer potential and this leads to a discrete
SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 15 F igure
7. The centerline defined by (26) comes near to self-intersection when H ≈
1; here weillustrate the curve with H = .
8. The arrows depict the imposed centerline force density f ( s )given by (27). The three black arrows indicate the forcing imposed at the centerline points s = π/ s = π , s = π/ ˆ z , ˆ x , and − ˆ z . ✏ gap = 110
8. The breakdown of Keller-Rubinow slender body theory when the fiber surface ap-proaches itself is local. In the top panel, the centerline is fixed while the fiber radius decreases; thediscrepancy between the KR centerline velocity u KRC ( s ) and our centerline velocity function c ( s )decreases uniformly in s . In the lower three panels, the centerline shifts ( H → /
10, 1 / . , or 1 / . In these three cases we see that the pointwise discrepancy decreases alongmost of the fiber length as (cid:15) →
0. However, the errors are O (1) at the three locations s = π/ s = π , and s = π/
9. The gray polygon shows the convergence of the Keller-Rubinow slender body theoryand our BVP-based method when the centerline is fixed and the fiber radius shrinks, at about thesame order we observed in the other tests, O ( (cid:15) . ). The dashed reference curve is O ( (cid:15) ). Theupper boundary of the gray polygon gives the ∞ -norm while the lower boundary is the 1-norm;the 2-norm values appear as × markers inside the polygon. The similarity of the convergence in the1- and ∞ -norms suggests that the convergence must be uniform in the centerline coordinate s , as isthe case (top panel of Fig. 8). The blue, yellow, and green colored polygons describe tests whereinthe fiber centerline approaches itself while its radius decreases simultaneously. In these cases the ∞ -norm stagnates while the 1-norm continues to decrease, suggesting that the breakdown of theKeller-Rubinow formulation is local in s (lower three panels of Fig. 8).system whose condition number is stable under refinement of the numerical quadrature, although it does grow asthe fiber radius shrinks and as we resolve more Fourier modes in the solution.Some immediate next steps would be to account for multiple fibers [18, 11], fibers with free ends [24, 13], anddynamic problems [33, 20]. Another possible extension of this work would be to consider inertial flows, where theunderlying PDE is di ff erent but the fiber integrity and average force conditions at the boundary are still a reasonableway to make physical sense of one-dimensional forcing data.While our method is more accurate in the presence of near contacts, it is much more computationally intensivethan the Keller-Rubinow formulation. For example, the largest simulations reported here required assembling andsolving a dense system of 8712 linear equations, wherein the majority of the matrix entries are di ffi cult integralsrequiring two-dimensional quadrature. To simulate many fibers or to move from static to dynamic problems, wewill need some combination of a cheaper algorithm, larger machines and / or parallelism. At larger scales, we willalso need to reconsider our choice to use the full SVD for the numerical linear algebra. A fast multipole methodmay be an appropriate strategy for larger sized problems [31].Given the adequate performance of the inexpensive Keller-Rubinow formulation for isolated sections of thefiber, a hybridization of the two methods might provide another way to reduce the computational expense relativeto the method presented here. 7. A cknowledgements WM and HB were supported by NSF grant DMS-1907796. DS was supported in part by NSF grant DMS-2009352. LO was supported by NSF Postdoctoral Research Fellowship DMS-2001959. YM was supported by NSF(DMS-1907583) and the Math + X grant from the Simons Foundation. We acknowledge the generous hospitalityand computational resources of the IMA, where the project was initiated.
SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 17 R eferences [1] AMD Amarakoon, RG Hussey, Bill J Good, and Edward G Grimsal. Drag measurements for axisymmetric motion of a torus at lowreynolds number. The Physics of Fluids , 25(9):1495–1501, 1982.[2] GK Batchelor. Slender-body theory for particles of arbitrary cross-section in Stokes flow.
Journal of Fluid Mechanics , 44(3):419–440,1970.[3] Andrea Louise Bertozzi.
An Extension of the Smale-Birho ff Homoclinic Theorem, Melnikov’s Method, and Chaotic Dynamics inIncompressible Fluids . PhD thesis, Princeton University, 1987.[4] Michael Boyle. The integration of angular velocity.
Advances in Applied Cli ff ord Algebras , 27(3):2345–2374, 2017.[5] Eduardo Corona, Leslie Greengard, Manas Rachh, and Shravan Veerapaneni. An integral equation formulation for rigid bodies instokes flow in three dimensions. Journal of Computational Physics , 332:504–519, 2017.[6] Ricardo Cortez, Lisa Fauci, and Alexei Medovikov. The method of regularized Stokeslets in three dimensions: analysis, validation,and application to helical swimming.
Phys. Fluids (1994-present) , 17(3):031504, 2005.[7] Ricardo Cortez and Michael Nicholas. Slender body theory for stokes flows with regularized forces.
Communications in AppliedMathematics and Computational Science , 7(1):33–62, 2012.[8] RG Cox. The motion of long slender bodies in a viscous fluid part 1. general theory.
Journal of Fluid Mechanics , 44(4):791–810,1970.[9] Arthur Joseph Goldman, Raymond G Cox, and Howard Brenner. Slow viscous motion of a sphere parallel to a plane wall—i motionthrough a quiescent fluid.
Chemical engineering science , 22(4):637–651, 1967.[10] Thomas G¨otz.
Interactions of fibers and flow: asymptotics, theory and numerics . Doctoral dissertation, University of Kaiserslautern,2000.[11] J H¨am¨al¨ainen, Stefan B Lindstr¨om, T H¨am¨al¨ainen, and H Niskanen. Papermaking fibre-suspension flow simulations at multiple scales.
J. Engrg. Math. , 71(1):55–79, 2011.[12] GJ Hancock. The self-propulsion of microscopic organisms through liquids.
Proceedings of the Royal Society of London. Series A.Mathematical and Physical Sciences , 217(1128):96–121, 1953.[13] Robert E Johnson. An improved slender-body theory for Stokes flow.
Journal of Fluid Mechanics , 99(02):411–431, 1980.[14] Eric E Keaveny and Michael J Shelley. Applying a second-kind boundary integral equation for surface tractions in stokes flow.
Journalof Computational Physics , 230(5):2141–2159, 2011.[15] Joseph B Keller and SI Rubinow. Swimming of flagellated microorganisms.
Biophysical Journal , 16(2):151–170, 1976.[16] Lyndon Koens and Eric Lauga. The boundary integral formulation of stokes flows includes slender-body theory.
Journal of FluidMechanics , 850, 2018.[17] Eric Lauga and Thomas R Powers. The hydrodynamics of swimming microorganisms.
Reports on Progress in Physics , 72(9):096601,2009.[18] Lei Li, Harishankar Manikantan, David Saintillan, and Saverio E Spagnolie. The sedimentation of flexible filaments.
Journal of FluidMechanics , 735:705–736, 2013.[19] James Lighthill. Flagellar hydrodynamics.
SIAM review , 18(2):161–230, 1976.[20] Anke Lindner and Michael Shelley. Elastic fibers in flows.
Fluid-Structure Interactions in Low-Reynolds-Number Flows , 168, 2015.[21] Andrew J Majda and Andrea L Bertozzi.
Vorticity and incompressible flow. Cambridge texts in applied mathematics . World PublishingCorporation, 2002.[22] William H Mitchell and Saverio E Spagnolie. A generalized traction integral equation for stokes flow, with applications to near-wallparticle mobility and viscous erosion.
Journal of Computational Physics , 333:462–482, 2017.[23] Yoichiro Mori, Laurel Ohm, and Daniel Spirn. Theoretical justification and error analysis for slender body theory.
Communicationson Pure and Applied Mathematics , 73(6):1245–1314, 2020.[24] Yoichiro Mori, Laurel Ohm, and Daniel Spirn. Theoretical justification and error analysis for slender body theory with free ends.
Archive for Rational Mechanics and Analysis , 235(3):1905–1978, 2020.[25] Ehssan Nazockdast, Abtin Rahimian, Denis Zorin, and Michael Shelley. A fast platform for simulating semi-flexible fiber suspensionsapplied to cell mechanics.
Journal of Computational Physics , 329:173–209, 2017.[26] M.E. O’Neill and B.S. Bhatt. Slow motion of a solid sphere in the presence of a naturally permeable surface.
The Quarterly Journalof Mechanics and Applied Mathematics , 44(1):91–104, 1991.[27] Christopher JS Petrie. The rheology of fibre suspensions.
Journal of Non-Newtonian Fluid Mechanics , 87(2):369–402, 1999.[28] Henry Power and Guillermo Miranda. Second kind integral equation formulation of stokes’ flows past a particle of arbitrary shape.
SIAM Journal on Applied Mathematics , 47(4):689–698, 1987.[29] Constantine Pozrikidis et al.
Boundary integral and singularity methods for linearized viscous flow . Cambridge University Press,1992.[30] Bruce Rodenborn, Chih-Hung Chen, Harry L Swinney, Bin Liu, and HP Zhang. Propulsion of microorganisms by a helical flagellum.
Proceedings of the National Academy of Science , 110(5):E338–E347, 2013.[31] V Rokhlin. Rapid solution of integral equations of classical potential theory.
Journal of Computational Physics , 60(2):187 – 207,1985.[32] Michael J Shelley. The dynamics of microtubule / motor-protein assemblies in biology and physics. Annual Review of Fluid Mechanics ,48:487–506, 2016.[33] Michael J Shelley and Tetsuji Ueda. The stokesian hydrodynamics of flexing, stretching filaments.
Physica D: Nonlinear Phenomena ,146(1-4):221–245, 2000. [34] Anna-Karin Tornberg and Michael J Shelley. Simulating the dynamics and interactions of flexible fibers in stokes flows.
Journal ofComputational Physics , 196(1):8–40, 2004.[35] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peter-son, Warren Weckesser, Jonathan Bright, St´efan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov,Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, ˙Ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, DenisLaxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antˆonio H. Ribeiro,Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing inPython.
Nature Methods , 17:261–272, 2020.[36] G. K. Youngren and A. Acrivos. Stokes flow past a particle of arbitrary shape: a numerical method of solution.
Journal of FluidMechanics , 69(2):377–403, 1975. A ppendix A. T he B ishop frame and fiber surface Here we give more specific information on our surface geometry construction using a quaternion-based initialvalue problem. We first describe the continuous formulation and then give numerical details.A.1.
Defining the Bishop frame using an initial value problem.
Let γ : [0 , π ] → R be a twice-di ff erentiableand periodic curve parameterizing the centerline of a closed fiber. We do not require an arclength parameterization,that is, we require only | γ (cid:48) | > | γ (cid:48) | =
1. The unit tangent vector and its derivative with respect to s aredefined by T ( s ) = | γ (cid:48) | γ (cid:48) , (28) T (cid:48) ( s ) = | γ (cid:48) | γ (cid:48)(cid:48) − γ (cid:48) · γ (cid:48)(cid:48) | γ (cid:48) | γ (cid:48) . (29)Because the vector T is perpendicular to its derivative, it traces out a closed loop on the unit sphere as s ranges from0 to 2 π . We wish to find vectors N ( s ) and B ( s ) which together with T ( s ) form an orthonormal frame. The simplestmethod, due to Frenet, is to set N = T (cid:48) / | T (cid:48) | and B = T × N ; these definitions are local in the sense that we needto know γ and its derivatives only at a fixed s in order to determine the frame { T ( s ) , N ( s ) , B ( s ) } . However, theFrenet normal and binormal vectors are undefined wherever T (cid:48) = , or equivalently whenever the acceleration γ (cid:48)(cid:48) is a scalar multiple of the velocity γ (cid:48) . A far-from-pathological example where this occurs is the boundary of theunit ball in the (cid:96) norm in the plane, which has the polar parameterization r ( s ) = (cid:16) cos s + sin s (cid:17) / . We thereforeuse an alternative global frame which is defined as the solution to a certain initial value problem. Let T = T (0)and let N and B be any vectors completing the frame at s =
0. Then define ω = ω ( s ) by(30) ω = T × T (cid:48) + | γ (cid:48) | α L T where α is a constant scalar to be determined later and L is the total length of the path γ . Now consider thedi ff erential equation of rotation(31) v (cid:48) = ω × v . If the initial condition is v = T , we find that the solution is v ( s ) = T ( s ); to see this, use the triple cross productformula ( A × B ) × C = ( C · A ) B − ( C · B ) A to write(32) (cid:18) T × T (cid:48) + α π T (cid:19) × T = ( T · T ) T (cid:48) − ( T · T (cid:48) ) T + = T (cid:48) . This is a global definition for T which is identical to the local version (28). However, we can also solve (31) withthe initial conditions v = N and v = B to obtain unit vector functions N ( s ) and B ( s ). These complete theframe for each s ∈ [0 , π ] , for if v and v are solutions of (31), we have dds ( v · v ) = v · ( ω × v ) + ( ω × v ) · v = α determines the speed of rotation about T and therefore a ff ects the evolution of N and B butnot T . We require a value of α giving a periodic frame: N (2 π ) = N (0) and B (2 π ) = B (0). This value of α is notunique; one can always add an integer to obtain another periodic frame where N and B twist around T a di ff erentnumber of times on [0 , π ]. One method of choosing α is to solve (31) with α = { N (2 π ) , B (2 π ) } to { N , B } ; then we set α to be the angle by which the initial and final N and B di ff erand solve (31) again. SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 19
By defining the scalar functions κ = | γ (cid:48) | T (cid:48) · N , κ = | γ (cid:48) | T (cid:48) · B , κ = α L , we can make the di ff erential equation(31) equivalent to the usual formulation of the Bishop frame, which is to start with arbitrary κ i and evolve thedi ff erential equation(33) dds TNB = | γ (cid:48) | κ κ − κ κ − κ − κ TNB . There is a factor of | γ (cid:48) | on the right side of (33) because we have not assumed an arclength parameterization.We can now use the centerline curve γ and the orthonormal frame { T , N , B } to define the fiber surface: let(34) X ( s , θ ) = γ ( s ) + (cid:15) cos( θ ) N ( s ) + (cid:15) sin( θ ) B ( s ) . The Jacobian of this transformation and the normal vector pointing out of the fiber into the fluid can be calculateddirectly from X θ × X s by writing all vectors in the right-handed basis { T , N , B } . The surface normal is(35) ν ( s , θ ) = cos( θ ) N ( s ) + sin( θ ) B ( s )and the area element is(36) J ( s , θ ) = (cid:15) | γ (cid:48) | (cid:16) − (cid:15) cos( θ ) κ ( s ) − (cid:15) sin( θ ) κ ( s ) (cid:17) . A.2.
Quaternion-based numerical calculation of the Bishop frame.
We will focus on the case where the cen-terline γ is defined as the Fourier interpolant of some discrete set of points in R rather than by a symbolic formula.In this setting we wish to compute the frame vectors and their derivatives accurately. To reduce the number of inde-pendent quantities which must be numerically integrated from six (for N and B , since T is known independently)to four, we use a formulation in terms of quaternions. We write quaternions as vector-scalar pairs: q = ( z , r ) is thesame as q = z i + z j + z k + r . Then we employ a reformulation of (31) derived by [4]:(37) q (cid:48) =
12 ( ω , ∗ q . Here the symbol ∗ denotes the noncommutative multiplication of quaternions, and a reformulation of (37) usingthe vector dot and cross products is(38) q (cid:48) = (cid:0) z (cid:48) , r (cid:48) (cid:1) = ( r ω + ω × z , − ω · z ) . We then solve (37) with the initial condition q = ((0 , , ,
1) to obtain a solution as a path in R . From this wecan obtain T and N and B from ( T ( t ) , = q ( t ) ∗ ( T , ∗ q ( t ) − (39) ( N ( t ) , = q ( t ) ∗ ( N , ∗ q ( t ) − (40) ( B ( t ) , = q ( t ) ∗ ( B , ∗ q ( t ) − . (41)This formulation has the advantage that the computed frame { T , N , B } is exactly a rotated version of { T , N , B } ,even in the presence of numerical errors in the computation of q ( t ). The disadvantage is that numerical error in q ( t )may cause the globally computed T to di ff er from the local one, i.e. (39) may not be identical to (28). One mayassess the accuracy of the integration by comparing the two versions. Furthermore, one can take advantage of thelocal method of computing T by adding an auxiliary rotation spinning the global version of T towards the localone. This amounts to adding a term to (30):(42) ω ( t ) = T ( t ) × T (cid:48) ( t ) + α π T ( t ) + T num ( t ) × T ( t ) . In (42) the vector T is computed from the local definition (28) while the vector T num is the global version (39).Typically these are nearly identical and so their cross product is small, i.e. (30) and (42) are very similar. Theconstant factor of 10 is a heuristic which improves the accuracy of the computation without significantly increasingthe sti ff ness of the problem. We use the SciPy integration routine solve ivp [35] to solve (37) with both absoluteand relative error tolerances set to 5 · − . singularity at s = s ⇤ = ⇡
10. The quadrature is designed to integrate a singularity of order 1 / r at a source point.Above, we show the rule generated for a fiber surface which approaches the origin three times,specifically for the centerline (26) with H = . (cid:15) equal to one fifth of the gap size. Thecolor on the fiber surface indicates distance from the source point, or more precisely | γ ( s ) − γ ( s ∗ ) | .The source point ( s ∗ = π, θ ∗ =
0) is visible at the center of the figure on the rear fiber branch.The quadrature nodes cluster in nested squares around the singularity, as the right-hand figureindicates. Below this we plot the same quadrature nodes in the flat s θ -plane. The innermostregion is decomposed into six triangles and the outer region is broken into rectangles where weuse rescaled Gauss-Legendre-trapezoid grids. If the singularity moves in the circumferential ( θ )direction, the quadrature rule can be periodically shifted and reused; for di ff erent centerline values s ∗ the rule must be regenerated.A ppendix B. Q uadrature for single - and double - layer potentials on a slender fiber We now construct quadrature rules for the single- and double-layer potentials on thin tubes. Let the centerline γ ( s ) and fiber surface X ( s , θ ) be constructed as above. If the integrand were smooth, the doubly periodic geometrywould be well suited to quadrature on a regular grid (double trapezoid rule). However, we have not found a way totake advantage of the double periodicity for the singular integrands required in the velocity and traction integrals(16)-(17). Instead, we decompose the fiber surface into subregions by s . The inner region consists of all s satisfying | γ ( s ) − γ ( s ∗ ) | < (cid:15) ; this is then further split into six triangles, each of which has a vertex at the source point. Theouter region is divided into subrectangles where we use an exponentially rescaled Gauss-Legendre integration inthe s -direction and the trapezoid rule in the θ -direction, thus taking advantage of periodicity in one of the twodirections. The procedure accounts for situations where the fiber centerline approaches itself and resolves theintegrand more finely as needed. An example of the resulting quadrature nodes is given in Fig. 10. Our quadraturegeneration procedure relies on knowledge of the fiber centerline geometry and the critical centerline value s ∗ ; thusthe rules must be regenerated for each new centerline and even for each new value of s ∗ . However, they canbe shifted periodically and reused when s ∗ is fixed but θ ∗ changes. The accuracy and size of the quadrature iscontrolled by a single parameter q n and the total number of nodes is O ( q n ). The results presented in this paper usedvalues of q n ranging from 30 to 45. We give a more detailed description in the following subsections. SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 21
B.1.
Inner region: integration of the singularity using triangles.
We begin by describing a quadrature rule forthe square [ − , which is well adapted for a 1 / r singularity with angular dependence at the origin. Consider firstthe triangle T with vertices at (0 , ,
1) and (1 , T can be transformed to an integral on thesquare [0 , via(43) (cid:90) T ψ dA = (cid:90) (cid:90) ψ ( uv , u ) u dv du . This mapping expands the region around the singularity, providing a regularizing factor which yields a smooth andbounded integrand on the square. We use an outer product of Gauss-Legendre integration rules in the uv -domain,which leads to clustering of the nodes around the singularity in the triangular domain.We use similar strategies for the other five subtriangles depicted on the right panel of Fig. 10. This leads to aquadrature on the square [ − , . Finally, we shift and scale the square to the inner (singular) integration regionin the s θ -plane, that is, with θ ∈ ( θ ∗ − π, θ ∗ + π ) and s satisfying | γ ( s ) − γ ( s ∗ ) | < (cid:15) . The right and left halves arescaled di ff erently in s if | γ (cid:48) | is not constant near s ∗ .The total number of quadrature nodes in this inner region is 6 q n . We have not carefully optimized this quadraturefor the inner region because the outer region generally requires more nodes.B.2. Outer region: integration over multiple scales using rectangles.
It remains to integrate over a regionwhere the integrand is smooth, that is, the part of the fiber surface where | γ ( s ) − γ ( s ∗ ) | > (cid:15) . We use an outerproduct of one-dimensional rules: in the circumferential coordinate θ we use the trapezoid rule, while for thecenterline coordinate s we first subdivide and then use an exponentially rescaled Gauss-Legendre integration oneach subinterval.The subdivision in s allows us to respond to the fact that the integration on the outer region becomes morechallenging as (cid:15) →
0. Indeed, the length of the outer domain in s is approximately 2 π − (cid:15) | γ (cid:48) ( s ∗ ) | , and thesingularity therefore lies at a distance of 5 (cid:15) | γ (cid:48) ( s ∗ ) | from either end of the integration interval; thus as (cid:15) → (cid:15) is small we want to put moresubintervals near the ends of the integration region.The possibility of near self-intersection of the fiber centerline presents another integration challenge. In thissituation we may have a large increase in the value of the integrand inside a very narrow subregion which couldappear anywhere on the subinterval. We would like to concentrate some small integration subpanels in such self-intersection zones.We propose a method which addresses these issues, at the cost of becoming somewhat complicated and special-izing to a specific fiber centerline and a specific value of s ∗ . We begin by placing subdivision markers at any valueof s where h ( s ) = | γ ( s ) − γ ( s ∗ ) | has a local extremum. Additionally, we place subdivision markers at any locationwhere log ( h ( s ) / (5 (cid:15) )) has a positive integer value. This results in a concentration of fiber subintervals in any lo-cation where the centerline approaches γ ( s ∗ ), including the ends of the outer integration region; the concentrationbecomes more extreme as the fiber radius shrinks. However, if there are many local extrema this procedure mayproduce too many subdivisions in places where the integrand is relatively smooth. Thus, we remove the subdivi-sion markers at the extrema unless a centered-di ff erence estimation of d ds log ( h ( s )) gives a relatively large value(greater than four times the average on a regular 30-point grid).Once the subdivisions are finalized, we create outer-product rules for each subrectangle. We always use n q trapezoidal nodes in the periodic θ -direction, while the number of nodes in the s -direction varies with thelength of the subinterval and the logarithmic change in the distance to the singularity over that subinterval. Con-cretely, the number of quadrature nodes in s is the maximum of q n / q n log ( h max / h min ) / + q n ( s l − s r ) / (2 π ). The integration in s is done via Gauss-Legendre quadrature after anexponential transformation.The exponential transformation is as follows. Let the integration interval be [ s l , s r ] × [ − π, π ]. We already knowthe values h ( s l ) and h ( s r ) from the subdivision procedure referenced earlier. If these are equal, we use ordinaryGauss-Legendre integration in s . If the endpoint values of h are unequal, that is, if one end of the correspondingcenterline section is closer to the singularity than the other, we carry out the transformation(44) (cid:90) s r s l ψ ( s ) ds = (cid:90) ψ (cid:16) A + B exp( Ct ) (cid:17) · BC exp( Ct ) dt s R e s c a l e d d i s t a n c e t o s i n g u l a r i t y f ( s ) F igure
11. If the centerline is a figure-eight shape with a near self-intersection, γ (0) ≈ γ ( π ),we need a quadrature rule that integrates the region near the singularity at s = s = π carefully. We construct one by splitting the arclength s -domain into subregionswith boundaries located where the distance to the singularity takes geometrically spaced values,or at extrema where the second derivative is especially large (for example, there are three extremawithin the subinterval with green dots, but they are ignored; in contrast, the minimum at s = π becomes a subinterval boundary). The number of quadrature points on each subinterval is equalto the quadrature parameter ( q n = s or if f ( s ) takes on very di ff erent values inside the interval.The logic is complicated but the algorithm runs quickly and allows the use of smaller quadraturerules than would be necessary without this adaptive procedure.where(45) A = s r h ( s l ) − s l h ( s r ) h ( s r ) − h ( s l ) , B = ( s l − s r ) h ( s l ) h ( s r ) − h ( s l ) , C = log( h ( s r ) / h ( s l )) . This transformation comes from the fact that s ( t ) = A + Be Ct is the solution of the first-order boundary valueproblem(46) s (0) = s (cid:96) , s (1) = s r , dsdt = C (cid:32) s − s l s r − s l h ( s r ) + s r − ss r − s l h ( s l ) (cid:33) where the expression in large parentheses is the linear interpolation of h ( s ) on [ s l , s r ], and C is a proportionalityconstant whose value is determined as part of the solution of the BVP. The consequence is that the product of ds / dt and 1 / h ( s ) should be approximately constant, so the transformed problem can be integrated with fewer quadraturenodes. We apply Gauss-Legendre integration on the right-hand side of (44).B.3. Convergence of the quadrature procedure.
As the quadrature is an important ingredient in our overallnumerical method, and also an independent problem that may have value beyond the present paper, we presentsome convergence results for three model integrands. These problems have similar behavior to those actuallyneeded for our matrix assembly but we have simplified them somewhat to make the results in this section morereplicable by others.The integration domains are the tubes whose centerlines are given by (26) using H = .
9, and whose circularcross sections have radius (cid:15) ∈ { . , . , . , . } . We do not reparameterize or rescale the centerline.On each of these domains we consider three integrands which are singular at the point(47) X ( s ∗ , θ ∗ ) = / √ / + (cid:15) / √ / . SINGLE-LAYER BASED NUMERICAL METHOD FOR THE SLENDER BODY BOUNDARY VALUE PROBLEM. 23
This corresponds to taking s ∗ = π/ θ ∗ = N ( π/ = (0 . , . √ , r = r ( s , θ ) = X ( s , θ ) − X ( s ∗ , θ ∗ ), thedefinitions of our three model integrals are: I / r = (cid:90) π (cid:90) π r J ( s , θ ) d θ ds (48) I S L = π (cid:90) π (cid:90) π ( r · ˆ x )( r · ˆ y ) r cos(50 s + θ ) J ( s , θ ) d θ ds (49) I DL = − π (cid:90) π (cid:90) π ( r · ˆ x )( r · ˆ y )( r · ν ) r cos(50 s + θ ) J ( s , θ ) d θ ds . (50)The first model integrand is the simplest: it is just the reciprocal of the distance to the source point. This firstintegrand has a 1 / r singularity at the source point, but none of the angular dependence on the divergence ratethat the single- and double-layer integrands display. The second integrand is one component of the single-layerpotential multiplied by a Fourier basis function. The third is one component of the double-layer potential multipliedby a Fourier basis function. We have omitted the point-source terms for simplicity and because they are regular.As the results of Table 4 indicate, we generally get thirteen to fifteen digits for the first integrand using between10000 and 20000 quadrature nodes. The single- and double-layer integrals I S L and I DL approach zero more quicklyas (cid:15) →
0, but if we count the leading zeros as correct digits then we can claim a similarly good performance. Toimprove on the last block of the table, that is, to get more than five significant figures in I DL for the smallest value of (cid:15) , presents an interesting numerical challenge but is not a likely source of significant error in the overall numericalmethod presented here.The quadrature procedure presented here produces good results in the tests we considered, but we hope in thefuture to find a simpler and more lightweight method. (cid:15) q n total nodes 1 / r SL DL5 · − . e +
00 8 . e − − . e − . e +
00 8 . e − − . e − . e +
00 8 . e − − . e − . e +
00 8 . e − − . e − . e +
00 8 . e − − . e − . e +
00 8 . e − − . e − . e +
00 8 . e − − . e − . e +
00 8 . e − − . e − · − . e −
01 6 . e −
06 5 . e − . e −
01 6 . e −
06 6 . e − . e −
01 6 . e −
06 6 . e − . e −
01 6 . e −
06 6 . e − . e −
01 6 . e −
06 6 . e − . e −
01 6 . e −
06 6 . e − . e −
01 6 . e −
06 6 . e − . e −
01 6 . e −
06 6 . e − · − . e −
01 2 . e −
08 5 . e − . e −
01 2 . e −
08 5 . e − . e −
01 2 . e −
08 5 . e − . e −
01 2 . e −
08 5 . e − . e −
01 2 . e −
08 5 . e − . e −
01 2 . e −
08 5 . e − . e −
01 2 . e −
08 5 . e − . e −
01 2 . e −
08 5 . e − · − . e −
02 1 . e −
09 5 . e − . e −
02 1 . e −
09 5 . e − . e −
02 1 . e −
09 5 . e − . e −
02 1 . e −
09 5 . e − . e −
02 1 . e −
09 5 . e − . e −
02 1 . e −
09 5 . e − . e −
02 1 . e −
09 5 . e − . e −
02 1 . e −
09 5 . e − able
4. The numerical integration is challenging because of the singularity and the large aspectratio of the fiber surface. Here we report convergence results for three model integrands over fourorders of magnitude in the fiber radius. We get good results for the simplest problem I / r usingfewer than 10000 function evaluations; for the modified velocity and traction integrals I S L and I DL the accuracy is also good for large (cid:15) but decays when (cid:15) →
0. The results for I DL with (cid:15) = · −5