An integral model based on slender body theory, with applications to curved rigid fibers

Andersson, Helge I.; Celledoni, Elena; Ohm, Laurel; Owren, Brynjulf; Tapley, Benjamin K.

doi:10.1063/5.0041521

We propose a novel integral model describing the motion of both flexible and rigid slender fibers in viscous flow and develop a numerical method for simulating dynamics of curved rigid fibers. The model is derived from nonlocal slender body theory (SBT), which approximates flow near the fiber using singular solutions of the Stokes equations integrated along the fiber centerline. In contrast to other models based on (singular) SBT, our model yields a smooth integral kernel which incorporates the (possibly varying) fiber radius naturally. The integral operator is provably negative definite in a nonphysical idealized geometry, as expected from the partial differential equation theory. This is numerically verified in physically relevant geometries. We discuss the convergence and stability of a numerical method for solving the integral equation. The accuracy of the model and method is verified against known models for ellipsoids. Finally, we develop an algorithm for computing dynamics of rigid fibers with complex geometries in the case where the fiber density is much greater than that of the fluid, for example, in turbulent gas-fiber suspensions.

I. INTRODUCTION

The dynamics of thin fibers immersed in fluid play an important role in many biological and engineering processes, including micro-organism propulsion,^8,28,43,49 rheological properties of fiber suspensions used to create composite materials,^16,20,41 and deposition of microplastics in the ocean.³¹ Here the term “fiber” is used to refer to a particle with a very large aspect ratio. In many of the applications mentioned, the cross sectional radius of the fiber is small compared to the length scales of the surrounding fluid, which can be well approximated locally by Stokes flow. This allows for the development of computationally tractable mathematical models describing the interaction between the fiber and the surrounding fluid.

Due to the linearity of the Stokes equations, the three dimensional flow about a body can be fully described by an expression over only the two dimensional surface of the body;⁴² however, for flexible particles with complex shapes or for multiple interacting particles, this quickly becomes both analytically and computationally prohibitive. In the case of slender fibers, a more tractable option is to exploit the thinness of the fiber by approximating it as a one dimensional curve. This idea forms the basis for slender body theory (SBT). Models based on slender body theory in general are popular because they yield simple, efficient expressions for the velocity of filaments in fluid, allowing for the simulation of many interacting fibers with complex, semiflexible shapes. The most basic form of SBT (placing singular point forces known as Stokeslets along the fiber centerline) dates back to works by Hancock,²¹ Cox,¹⁴ and Batchelor.⁴ Later developments in singular SBT, due to Keller and Rubinow,²⁵ Lighthill,²⁹ and Johnson,²⁴ involved adding higher order corrections to the point force to account for the finite radius of the fiber. The most natural choice of higher order correction is often referred to as the doublet [see discussion following Eq. (9)]. We will refer to these methods based on distributing Stokeslets and doublets along the fiber as classical nonlocal SBT to distinguish from some more recent developments.

Classical SBT gives rise to an expression which exactly satisfies the unforced Stokes equations away from the fiber, and, to leading order (with respect to the fiber radius) satisfies the boundary conditions for a well-posed boundary value problem for the Stokes equations.^35,36 This expression has served as the basis for various numerical methods.^45,54,55 However, one issue with classical SBT is that the velocity expression is singular along the fiber centerline, and the usual methods for obtaining an expression for the velocity of the fiber itself—involving a nonstandard finite part integral—give rise to high wavenumber instabilities.^18,45,55 To address this, Tornberg–Shelley⁵⁵ regularize the integral kernel using an additional parameter proportional to the fiber radius.

To more generally avoid some of the difficulties of integrating a singular kernel, Cortez^10,12,13 developed the method of regularized Stokeslets. Here, instead of placing singular solutions of the Stokes equations along the fiber centerline, regularized Stokeslets are used. Regularized Stokeslets satisfy the Stokes equations with forcing given by a smooth approximation to the identity—or blob function—whose width is controlled by a parameter which can be chosen to be proportional to the fiber radius. Unlike classical SBT, this results in an expression for the fluid velocity that is nonsingular along the actual centerline of the fiber, allowing for a simpler representation of the fiber velocity. Many recent computational models for thin fibers rely on the method of regularized Stokeslets.^{5,11,48,58,59} However, many choices of blob function are possible and there is not a canonical procedure for choosing one. Additionally, many commonly used blob functions introduce an additional nonzero body force into the fluid away from the fiber surface.⁶¹

Most recently, Maxian et al.³² developed a fiber model that is asymptotically equivalent to SBT but based on the Rotne–Prager–Yamakawa (RPY) tensor^44,60 commonly used to model hydrodynamically interacting spheres. The model also places a curve of (singular) Stokeslets plus doublets along the fiber centerline, but replaces the region around the singular part of the Stokeslet/doublet kernel with the RPY regularization. The RPY kernel is divergence-free and known to be positive definite, making it a good choice for modeling particles in close proximity. The discontinuous kernel, however, makes the model more difficult to compare to the partial differential equation (PDE) solution of Refs. 35 and 36, which is one of the main goals of the model presented here.

We aim to make use of the fact that classical SBT closely approximates the solution to a well-posed boundary value problem^35,36 for the fluid velocity outside of the fiber, although the conventional way to obtain an expression for the velocity of the filament itself gives rise to instabilities which must later be corrected. Regularized Stokeslets yield a simpler expression for the fiber velocity, but can introduce errors outside of the filament and give rise to a fiber velocity which may fundamentally differ from the aforementioned PDE solution (see Remark III.2). Thus we consider a different approach to deriving a fiber velocity expression from classical SBT. Beginning with the fundamental premise of classical SBT—placing singular Stokeslets along the fiber centerline along with doublets to cancel the angular dependence across each fiber cross section—we aim to devise a model which is analytically and computationally attractive (in that it does not exhibit high wavenumber instabilities) with a physically meaningful derivation.

Our integral model is based on classical SBT but involves a smooth kernel which incorporates the (possibly varying) fiber radius in a natural way. Since the integral kernels are smooth, the model resembles the method of regularized Stokeslets with an arc length-dependent regularization similar to Ref. 58; however, we derive our model from usual (singular) Stokeslets and doublets. As such, we avoid introducing a nonzero body force throughout the fluid outside of the fiber,⁶¹ and avoid introducing additional parameters into the basic first-kind formulation of the model. The model relies on the asymptotic cancelation of angular-dependent terms along the fiber surface (see Sec. III for details), leaving an expression that retains a dependence on the fiber radius in a natural way.

Furthermore, we include a systematic way of comparing mapping properties among different fiber models based on Ref. 34, which involves calculating the spectra of the integral operators from various models in the toy scenario of a straight-but-periodic fiber with constant radius. In this model geometry, our integral operator is negative definite, as is the well-posed partial differential equation (PDE) operator of Refs. 35 and 36 which it is designed to approximate (see Ref. 34). This is in contrast to other models based on (non-regularized) slender body theory which give expressions for the fiber velocity involving further asymptotic expansion with respect to the fiber radius.^24,25,29 These models exhibit an instability as the eigenvalues of the operator cross zero at a high but finite wavenumber.

The model we derive initially yields a first-kind Fredholm integral equation for the force density along the fiber centerline. Such integral equations are known for being ill-posed (see Ref. 26, Chapter 15.1), as they do not necessarily have a bounded inverse at the continuous level. Numerical discretization alone can provide sufficient regularization to invert first-kind integral equations at the discrete level, but to make our model more suitable for inversion, we use an integral identity to regularize the expression into a second-kind equation. The second-kind regularization preserves the asymptotic accuracy of the model while improving the conditioning and invertibility of the corresponding numerical method. The regularization also serves to ensure that the discretized operator is negative definite, even in the presence of numerical errors, by bounding the spectrum away from zero. We distinguish this type of regularization from the method of regularized Stokeslets, since our regularization is not a key component of the model derivation. In particular, we can directly compare our model with regularization to our model without, which we will do repeatedly throughout the paper. We also distinguish this regularization from the procedure used by Tornberg and Shelley,⁵⁵ since we are not correcting for a high wavenumber instability. This allows us to compare the numerical behavior of our regularized and unregularized models at the discrete level even for very fine discretization. Moreover, the regularization used here affects all directions (both normal and tangent to the slender body centerline) in the same way.

The solution of the resulting second-kind Fredholm integral equation is a force density along the slender body centerline which we integrate to find the total force and torque on the rigid fiber. We implement a numerical method based on the Nyström method with Gauss–Legendre quadrature for solving the second-kind Fredholm integral equation (see Ref. 2, Chapter 12.4). Numerical tests confirm its convergence. Not surprisingly, we note significant improvements in the conditioning of the second-kind vs first-kind formulation of the model. We also numerically verify the spectral properties of the model in different geometries.

The model applies to both semiflexible and rigid fibers; however, the invertibility properties of the second kind model make it particularly well suited for simulating rigid filaments. We present an algorithm for dynamic simulations of rigid fibers where the fiber density is assumed to be much greater than that of the fluid, for example, in turbulent gas-fiber suspensions. The rigidity of the fiber can be exploited such that only matrix-vector products need to be performed within the time loop. We compare the dynamics of our model to the well-studied dynamics of a slender prolate spheroid.^6,9,23 We then apply our model to compare the dynamics of curved fibers whose centerlines deviate randomly from straight lines by varying magnitudes.

The structure of the paper is as follows. Section II presents the slender body model, which is derived in greater detail and justified via spectral comparisons with other slender body theories in Sec. III. In Sec. IV we discuss a method for numerically solving Fredholm integral equations and integrating the result, and demonstrate the convergence of the method for our model. Section V outlines a fast algorithm for computing the dynamics of a rigid slender fiber in viscous flow. We apply the dynamical algorithm to simulate the dynamics of fibers with complex shapes. Finally, we comment on conclusions and outlook for the model in Sec. VI.

A. Fiber geometry

We begin by introducing some notation for the slender geometries considered throughout the paper. Fix ε, L with $0 < ε ≪ L$ and let $X_{ext} : [- \sqrt{L^{2} + ε^{2}}, \sqrt{L^{2} + ε^{2}}] \to ℝ^{3}$ denote the coordinates of a C² curve in $ℝ^{3}$ ⁠, parameterized by arc length s. Defining $e_{s} (s) = \frac{d X_{ext}}{d s} / | \frac{d X_{ext}}{d s} |$ ⁠, the unit tangent vector to $X_{ext} (s)$ ⁠, we parameterize points near $X_{ext} (s)$ with respect to the orthonormal frame $(e_{s} (s), e_{n_{1}} (s), e_{n_{2}} (s))$ defined in Ref. 36. Letting

e_{r} (s, θ) : = cos θ e_{n_{1}} (s) + sin θ e_{n_{2}} (s),

we define the slender body $Σ_{ε}$ as

Σ_{ε} : = {x \in ℝ^{3} : x = X_{ext} (s) + ρ e_{r} (s, θ), ρ < ε r (s), s \in [- \sqrt{L^{2} + ε^{2}}, \sqrt{L^{2} + ε^{2}}]} .

(1)

Here the radius function $r \in C^{2} (- \sqrt{L^{2} + ε^{2}}, \sqrt{L^{2} + ε^{2}})$ is required to satisfy $0 < r (s) \leq 1$ for each $s \in (- \sqrt{L^{2} + ε^{2}}, \sqrt{L^{2} + ε^{2}})$ ⁠, and r(s) must decay smoothly to zero at the fiber endpoints $\pm \sqrt{L^{2} + ε^{2}}$ ⁠. There are many admissible radius functions r which can be considered. For the simulations in this paper, we will use a thin prolate spheroid as our geometrical model for a slender fiber. In this case, the radius function r(s) is given by

r (s) = \frac{1}{\sqrt{L^{2} + ε^{2}}} \sqrt{L^{2} + ε^{2} - s^{2}} .

(2)

We consider the subset

X : = {X_{ext} (s) : - L \leq s \leq L}

(3)

extending from focus to focus of the prolate spheroid (2), and define $X (s)$ to be the effective centerline of the slender body so that $r = O (ε)$ at the effective endpoints $s = \pm L$ ⁠.

The slender body model described in Sec. II may also be used in the case of a closed curve, in which case we take $X (L) = X (- L)$ and consider $s \in ℝ / 2 L$ ⁠. We may take the radius function $r \equiv 1$ in this case.

II. SLENDER BODY MODEL

To describe the motion of the thin fiber $Σ_{ε}$ (1) in Stokes flow, we will use an expression derived from the classical nonlocal slender body theory.^18,24,55 Letting $f (s, t)$ denote the force per unit length exerted by the fiber on the surrounding fluid at time t, we approximate the velocity $\frac{\partial X}{\partial t}$ of the fiber relative to a given background flow $u_{0}$ by

8 π μ (\frac{\partial X}{\partial t} - u_{0} (X (s, t), t)) = - 2 log (η) f (s, t) - \int_{- L}^{L} (S_{ε, η} + \frac{ε^{2} r^{2} (s')}{2} D_{ε}) f (s', t) d s',

(4)

S_{ε, η} (s, s', t) = \frac{I}{{(| \bar{X} |^{2} + η^{2} ε^{2} r^{2} (s))}^{1 / 2}} + \frac{\bar{X} {\bar{X}}^{T}}{{(| \bar{X} |^{2} + ε^{2} r^{2} (s))}^{3 / 2}},

(5)

D_{ε} (s, s', t) = \frac{I}{{(| \bar{X} |^{2} + ε^{2} r^{2} (s))}^{3 / 2}} - \frac{3 \bar{X} {\bar{X}}^{T}}{{(| \bar{X} |^{2} + ε^{2} r^{2} (s))}^{5 / 2}},

(6)

where $\bar{X} (s, s', t) = X (s, t) - X (s', t)$ ⁠. Here $η \geq 1$ is a parameter which can be chosen to yield either a first kind (η = 1) or a second-kind (⁠ $η > 1$ ⁠) Fredholm equation for f. Notice that η must also appear in the first term of $S_{ε, η}$ in order to retain the asymptotic consistency of the model (4). This is due to an integral identity (14) used to convert the integral model from a first-kind equation for f. The model accounts for a varying radius r(s) through the denominators of each term as well as the coefficient of $D_{ε}$ ⁠. Note that since r(s) is nonzero for $- L \leq s \leq L$ ⁠, the integral kernel is smooth for each $s \in [- L, L]$ ⁠. We provide a more detailed derivation of Eqs. (4)–(6) in Sec. III. We note that when η = 1, the expression (4) looks a bit similar to SBT using regularized Stokeslets, but—as we detail in the next section—the expression is derived through different means (namely, the near-cancelation of angular dependent terms in classical SBT along the surface of the filament) and, in particular, the appearance of $ε^{2} r^{2}$ in the denominator is not ad hoc but rather the best approximation of the slender body PDE of Refs. 35 and 36.

The model given by Eqs. (4)–(6) and the analysis in Sec. III can be used to describe both flexible and rigid fibers. In Sec. V we apply our model to the dynamics of a rigid fiber, since the invertibility properties of Eqs. (4)–(6) make the model especially suitable for simulating rigid filaments.

In the case of a rigid fiber, at each time t we additionally impose the constraint

\frac{\partial X}{\partial t} = v + ω \times X (s),

(7)

where $v$ ⁠, $ω \in ℝ^{3}$ are the linear and angular velocity of the fiber (see Refs. 19, 33, and 54). The total force $F (t)$ and torque $T (t)$ exerted on the slender body at time t are computed from the line force density $f (s, t)$ via

\int_{- L}^{L} f (s, t) d s = F (t), \int_{- L}^{L} X (s, t) \times f (s, t) = T (t) .

(8)

When $v$ and $ω$ are prescribed and one aims to solve for F and T, this is known as the resistance problem. Conversely, the case when F and T are given and the rigid fiber velocity is sought is known as the mobility problem. Note that for both the resistance and mobility problems along a thin fiber, using Eq. (4) to relate fiber velocity to force involves inverting the integral equation to solve for the force density f. Thus we are particularly concerned with the invertibility of Eq. (4). In Sec. V, we use Eqs. (4), (7), and (8) to solve the resistance problem, which is of interest when the density of the fiber is much larger than the density of the fluid.

III. DERIVATION AND JUSTIFICATION OF THE SLENDER BODY MODEL

Our model for the motion of the fiber is based on classical nonlocal slender body theory, where the fluid velocity $u^{SB} (x, t)$ at any point x away from the fiber centerline $X (s, t)$ is approximated by the integral expression

\begin{matrix} 8 π μ (u^{SB} (x, t) - u_{0} (x, t)) = - \int_{- L}^{L} (s (x - X (s', t)) + \frac{ε^{2} r^{2} (s')}{2} d (x - X (s', t))) f (s', t) d s' \\ s (x) = \frac{I}{| x |} + \frac{x x^{T}}{{| x |}^{3}}, d (x) = \frac{I}{{| x |}^{3}} - \frac{3 x x^{T}}{{| x |}^{5}}, \end{matrix}

(9)

where $u_{0} (x, t)$ is the fluid velocity in the absence of the fiber and μ is the fluid viscosity. The force-per-unit-length $f (s, t)$ exerted by the fluid on the body is distributed between the generalized foci of the slender body at $s = \pm L$ ⁠. The expression $\frac{1}{8 π μ} s (x)$ is the free space Green's function for the Stokes equations in $ℝ^{3}$ ⁠, commonly known as the Stokeslet, while $\frac{1}{8 π μ} d (x) = \frac{1}{16 π μ} Δ s (x)$ is a higher order correction to the velocity approximation, often known as a doublet. The doublet coefficient $\frac{ε^{2} r^{2}}{2}$ is chosen to cancel the leading order (in ε) angular dependence in the fluid velocity at the surface of the actual 3D filament. This coefficient can be obtained via matched asymptotics, or by the following heuristic. Since the purpose of the doublet is to cancel the angular dependence over each 2D cross section of the fiber, we consider Stokes flow in $ℝ^{2}$ due to a point force at the origin of strength f. In polar coordinates $x = {(ρ cos θ, ρ sin θ)}^{T}$ ⁠, the velocity due to the Stokeslet at $ρ > 0$ is given by

u^{s} (ρ, θ) = \frac{1}{4 π} (- log ρ I + \frac{1}{2} (\begin{matrix} 1 + cos 2 θ & sin 2 θ \\ sin 2 θ & 1 - cos 2 θ \end{matrix})) (\begin{matrix} f_{1} \\ f_{2} \end{matrix}),

where $I$ is the 2D identity matrix. To eliminate the θ-dependence on the circle $ρ = ε$ ⁠, we note that

\begin{matrix} Δ u^{s} (ρ, θ) = \frac{\partial^{2} u^{s}}{\partial ρ^{2}} + \frac{1}{ρ} \frac{\partial u^{s}}{\partial ρ} + \frac{1}{ρ^{2}} \frac{\partial^{2} u^{s}}{\partial θ^{2}} \\ = - \frac{1}{2 π ρ^{2}} (\begin{matrix} cos 2 θ & sin 2 θ \\ sin 2 θ & - cos 2 θ \end{matrix}) (\begin{matrix} f_{1} \\ f_{2} \end{matrix}) . \end{matrix}

Therefore, the θ-dependence in the velocity due to the Stokeslet at $r = ε$ can be canceled by adding a doublet term (⁠ $\frac{1}{2} Δ u^{s}$ ⁠) with coefficient $\frac{ε^{2}}{2}$ ⁠,

u^{SB} = u^{s} + \frac{ε^{2}}{4} Δ u^{s} .

The expression (9) is valid for describing flows around fibers which are not highly curved (i.e., with maximum centerline curvature $≪ 1 / ε$ ⁠) and do not come close to self-intersection (⁠ $| X (s) - X (s') | / | s - s' | \geq C$ for C independent of ε). The force density f must also be sufficiently regular. Given these constraints, in the stationary setting, the velocity field given by Eq. (9) is an asymptotically accurate approximation to the velocity field around a three-dimensional semi-flexible rod satisfying a well-posed slender body PDE, defined in Refs. 35 and 36 as the following boundary value problem for the Stokes equations,

\begin{matrix} - μ Δ u + \nabla p = 0, div u = 0 in ℝ^{3} \ \bar{Σ_{ε}} \\ \int_{0}^{2 π} (σ n) |_{(φ (s), θ)} j_{ε} (φ (s), θ) φ' (s) d θ = - f (s) on \partial Σ_{ε} \\ u |_{\partial Σ_{ε}} = u (s), unknown but independent of θ \\ | u | \to 0 as | x | \to \infty . \end{matrix}

(10)

Here $σ = μ (\nabla u + {(\nabla u)}^{T}) - p I$ is the fluid stress tensor, $n (x)$ denotes the unit normal vector pointing into $Σ_{ε}$ at $x \in \partial Σ_{ε}, j_{ε} (s, θ)$ is the Jacobian factor on $\partial Σ_{ε}$ ⁠, and $φ (s) : = \frac{s \sqrt{L^{2} + ε^{2}}}{L}$ is a stretch function to address the discrepancy between the extent of f and the extent of the actual slender body surface. Given a force density $f \in C^{1} (- L, L)$ which decays like r(s) at the fiber endpoints (⁠ $f (s) \sim r (φ (s))$ as $s \to \pm L$ ⁠), the difference between the slender body approximation $u^{SB}$ and the solution of Eq. (10) is bounded by an expression proportional to $ε | log ε |$ ⁠. Note that r(s) need not be spheroidal (2) for this error analysis to hold, but r(s) must decay smoothly to zero at the physical endpoints of the fiber at $s = \pm \sqrt{L^{2} + ε^{2}}$ ⁠.

In Sec. V we consider in greater detail the resistance problem for a rigid fiber, which is actually a case of the ‘inverse problem’ for Eq. (10): instead of prescribing the force density $f (s)$ along the filament, the fiber velocity $u (s)$ is given and we must solve for $f (s)$ ⁠. The slender body PDE (10) is then simply Dirichlet problem for the Stokes equations; however, it is unclear what type of decay in f (if any) is necessary for the SBT expression (9) to accurately approximate the PDE solution very near the fiber endpoints. Nevertheless, in Appendix D, we provide numerical evidence that the force density arising from inverting the expression (4) does exhibit decay at the fiber endpoints, both in the case of a prolate spheroid and a cylinder with hemispherical caps.

A key component of the well-posedness theory for the slender body PDE to which (9) is an approximation is the fiber integrity condition on $u |_{\partial Σ_{ε}}$ ⁠. The fiber integrity condition requires the velocity across each cross section s of the slender body to be constant; i.e., the velocity $u (x)$ at any point $x (s, θ) = X (s) + ε r (s) e_{r} (s, θ) \in \partial Σ_{ε}$ satisfies $\partial_{θ} u (x (s, θ)) = 0$ ⁠. This is to ensure that the cross-sectional shape of the fiber does not deform over time. An important aspect of the accuracy of slender body theory is that the expression (9) satisfies this fiber integrity condition to leading order in ε. Specifically, by Propositions 3.9 and 3.11 in Refs. 35 and 36, respectively, we have that for $x (s, θ) \in \partial Σ_{ε}$ ⁠,

| \partial_{θ} u^{SB} (x (s, θ)) | \leq C (ε | log ε | {| | f | |}_{C^{1} (- L, L)} + ε {‖ \frac{f}{r} ‖}_{C^{0} (- L, L)}),

(11)

i.e., the angular dependence in $u^{SB} (x)$ over each cross section s of the slender body is only $o (ε log ε)$ ⁠.

Another important general feature of the slender body PDE (10) is that the operator mapping the force data $f (s)$ to the θ-independent fiber velocity $u |_{\partial Σ_{ε}} (s)$ is negative definite [see (Ref. 34); note that the sign convention for f is opposite].

Now, the velocity expression (9) is singular at $x = X (s, t)$ and can be used only away from the fiber centerline; however, Eq. (9) presents a starting point for approximating the velocity of the slender body itself. Various methods can be used to obtain an expression for the relative velocity of the fiber centerline $\frac{\partial X (s, t)}{\partial t}$ which depends only on the arc length parameter s and time t. The most common way to go from Eq. (9) to an expression independent of θ is to perform an asymptotic expansion about ε = 0.^18,24,40,55 However, as alluded to in the Introduction, this leads to issues at high frequency modes along the fiber (we will come back to this point later). Here we consider a different approach to deriving a limiting centerline expression from Eq. (9) which evidently results in a negative definite integral operator mapping f to $u |_{\partial Σ_{ε}}$ ⁠. We then regularize this first-kind integral equation in an asymptotically consistent way to yield the second-kind integral equation (4). We detail our approach here and provide further justification in Sec. III A using a model geometry.

The first step in approximating $\frac{\partial X (s, t)}{\partial t}$ is to evaluate (9) on the surface of the slender body at $x = X (s, t) + ε r (s) e_{r} (s, θ, t)$ ⁠. Written out, the velocity field along the fiber surface is given by

\begin{matrix} 8 π μ (u^{SB} (x (s, θ, t), t) - u_{0} (X (s, t), t)) \\ = - \int_{- L}^{L} (\frac{I}{| R |} + \frac{\bar{X} {\bar{X}}^{T} + ε r (\bar{X} e_{r}^{T} + e_{r} {\bar{X}}^{T}) + ε^{2} r^{2} e_{r} e_{r}^{T}}{{| R |}^{3}} \\ + \frac{ε^{2} r^{2} (s')}{2} (\frac{I}{{| R |}^{3}} - 3 \frac{\bar{X} {\bar{X}}^{T} + ε r (\bar{X} e_{r}^{T} + e_{r} {\bar{X}}^{T}) + ε^{2} r^{2} e_{r} e_{r}^{T}}{{| R |}^{5}})) f (s', t) d s', \end{matrix}

(12)

where unless otherwise specified, we have $r = r (s), \bar{X} = \bar{X} (s, s', t) = X (s, t) - X (s', t)$ and $R = R (s, s', θ, t) = \bar{X} + ε r (s) e_{r} (s, θ, t)$ ⁠. Now, along the fiber surface, the expression (12) satisfies the fiber integrity condition to leading order in ε; i.e., the terms containing $e_{r} (s, θ, t)$ in Eq. (12) vanish to $o (ε log ε)$ ⁠, by Eq. (11). Because of this, to obtain an approximation to the velocity of the fiber itself which depends only on arc length, we could simply select a single curve along the length of the filament—i.e., fix $θ = θ^{*}$ or even $θ = θ^{*} (s)$ —and use the expression (12) evaluated along this curve as the approximate velocity of the fiber.⁴⁰ This yields an integral expression with a smooth, divergence-free kernel with clear physical meaning. However, this also involves a choice of $θ^{*}$ and subsequent computation of a normal vector at each point along the fiber, which is unnecessarily complicated given that we know from Eq. (11) that the terms containing θ are small.

In particular, both the Stokeslet and doublet include a θ-dependent term with $ε^{2} r^{2} e_{r} e_{r}^{T}$ in the numerator. Due to the form of R in the denominator, both of these terms are $o (1)$ at $s = s'$ ⁠; however, upon integrating in $s'$ ⁠, these terms cancel each other asymptotically to order $ε log ε$ ⁠. In particular, by Lemmas 3.5 and 3.7 in Refs. 35 and 36, respectively, we have

\begin{matrix} | \int_{- L}^{L} \frac{ε^{2} r^{2} e_{r} e_{r}^{T}}{| R |^{3}} f (s') d s' - 2 e_{r} e_{r} \cdot f (s) | \leq C ε ({| | f | |}_{C^{1} (- L, L)} + {‖ \frac{f}{r} ‖}_{C^{0} (- L, L)}), \\ | - \int_{- L}^{L} \frac{ε^{2} r^{2} (s')}{2} \frac{3 ε^{2} r^{2} e_{r} e_{r}^{T}}{| R |^{3}} f (s') d s' + 2 e_{r} e_{r} \cdot f (s) | \leq C ε ({| | f | |}_{C^{1} (- L, L)} + {‖ \frac{f}{r} ‖}_{C^{0} (- L, L)}) . \end{matrix}

As we can see, the O(1) contributions from both of these terms exactly cancel, leaving only higher order (in ε) contributions. Furthermore, the terms $ε r (\bar{X} e_{r}^{T} + e_{r} {\bar{X}}^{T})$ in both the Stokeslet and doublet approximately integrate to zero in $s'$ ⁠, since, by Lemmas 3.4 and 3.6 in Refs. 35 and 36, respectively, we have

| \int_{- L}^{L} ε^{m} r^{m} (s') \frac{ε r (\bar{X} e_{r}^{T} + e_{r} {\bar{X}}^{T})}{{| R |}^{m + 3}} f (s') d s' | \leq C ε (| log ε | {| | f | |}_{C^{1} (- L, L)} + {‖ \frac{f}{r} ‖}_{C^{0} (- L, L)}), m = 0, 2.

Finally, the $e_{r}$ term in each denominator from ${| R (s, θ, t) |}^{2} = {| \bar{X} |}^{2} + 2 ε r e_{r} \cdot \bar{X} + ε^{2} r^{2}$ is also only $O (ε log ε)$ ⁠, since, again using Lemmas 3.4 and 3.6 in Refs. 35 and 36,

| \int_{- L}^{L} (\frac{ε^{m} r^{m} (s') f (s')}{{| R |}^{m + 1}} - \frac{ε^{m} r^{m} (s') f (s')}{{(| \bar{X} |^{2} + ε^{2} r^{2})}^{\frac{m + 1}{2}}}) d s' | \leq C ε (| log ε | {| | f | |}_{C^{1} (- L, L)} + {| | \frac{f}{r} | |}_{C^{0} (- L, L)}), m = 0, 2.

Due to these cancelations and the fact that dropping these terms still approximates the slender body PDE solution of Refs. 35 and 36 to at least $O (ε log ε)$ ⁠, we may eliminate all terms containing $e_{r} (s, θ, t)$ in Eq. (12) to obtain a θ-independent expression which approximates the velocity of the fiber itself,

\begin{matrix} 8 π μ (\frac{\partial X}{\partial t} - u_{0} (X (s, t), t)) = - \int_{- L}^{L} (\frac{I}{{(| \bar{X} |^{2} + ε^{2} r^{2} (s))}^{1 / 2}} + \frac{\bar{X} {\bar{X}}^{T}}{{(| \bar{X} |^{2} + ε^{2} r^{2} (s))}^{3 / 2}} \\ + \frac{ε^{2} r^{2} (s')}{2} (\frac{I}{{(| \bar{X} |^{2} + ε^{2} r^{2} (s))}^{3 / 2}} - \frac{3 \bar{X} {\bar{X}}^{T}}{{(| \bar{X} |^{2} + ε^{2} r^{2} (s))}^{5 / 2}})) f (s', t) d s' . \end{matrix}

(13)

The expression (13) serves as the model underlying our final slender body velocity expression (4). Again, expression (13) looks somewhat similar to SBT using regularized Stokeslets, but is instead derived by the near-cancelation of angular-dependent terms in the classical SBT expression (9) along the fiber surface. In particular, $ε^{2} r^{2}$ remaining in the denominators here is simply what remains after eliminating these angular-dependent terms. In fact, due to the integral identity (14) which will introduced below as a means of converting (13) into a second-kind integral equation, we can see that altering this term severly affects the local behavior of the operator. For example, multiplying this $ε^{2} r^{2}$ denominator term by a constant other than 1 will introduce an O(1) dispartiy between the approximation and the slender body PDE unless corrected via an additional local term (see also Remark III.2 in Sec. III A 3 related to the method of regularized Stokeslets).

One further limitation to note about the centerline expressions (4) and (13) is that because the model is essentially 1D, in certain special cases (i.e., when the fiber is straight and its axis is perfectly aligned with the flow), the slender body approximation, in contrast to a truly 3D fiber, does not pickup on fluid gradients (see Sec. V B 1).

In Sec. III A, we show that in a simplified setting, (13) results in a negative definite operator mapping the force density f to the fiber velocity $\frac{\partial X}{\partial t}$ ⁠, whereas other models which rely on further asymptotic expansion of Eq. (12) about ε = 0 do not, and incur high wavenumber instabilities. This phenomenon is well known for the Keller–Rubinow model,^18,25 but for other possible centerline expressions, including models similar to Lighthill,²⁹ this high wavenumber instability has not been documented previously. It seems that our model (13) may be the simplest that can be obtained by expanding from Eq. (12) while still guaranteeing a negative definite operator.

Now, since the integral operator in Eq. (13) has a smooth kernel, the expression (13) yields a first-kind Fredholm integral equation for f when the fiber velocity $\frac{\partial X}{\partial t}$ is supplied. Describing the motion of a rigid fiber involves inverting this expression to solve for f, which in general is an ill-posed problem for a first-kind equation. Thus we want to regularize the integral operator (13) to create a second-kind integral equation while keeping the same order of accuracy in the map $f \mapsto \frac{\partial X}{\partial t}$ ⁠.

We first note that, for $η > 1$ ⁠, we have the following identity:

\int_{- L}^{L} (\frac{1}{{(| \bar{X} |^{2} + ε^{2} r^{2} (s))}^{1 / 2}} - \frac{1}{{(| \bar{X} |^{2} + η^{2} ε^{2} r^{2} (s))}^{1 / 2}}) g (s') d s' = 2 log (η) g (s) + o (η ε log (η ε)) .

(14)

Proof. By Lemma 3.8 in Ref. 36, for a > 0 sufficiently small, we have

\begin{matrix} \int_{- L}^{L} (\frac{g (s^{'})}{{(| \bar{X} |^{2} + a^{2} r^{2} (s))}^{1 / 2}} - \frac{g (s^{'})}{| \bar{X} |} + \frac{g (s)}{| s - s^{'} |}) d s^{'} \\ = log (\frac{2 (L^{2} - s^{2}) + 2 \sqrt{{(L^{2} - s^{2})}^{2} + a^{2} r^{2} (s)}}{a^{2} r^{2} (s)}) + o (a log a) . \end{matrix}

(15)

Subtracting (15) with $a = η ε$ from Eq. (15) with $a = ε$ and using that

| log (\frac{(L^{2} - s^{2}) + \sqrt{L^{2} + ε^{2} r^{2}}}{(L^{2} - s^{2}) + \sqrt{L^{2} + η^{2} ε^{2} r^{2}}}) | = | log (\frac{(L^{2} - s^{2}) + \sqrt{L^{2} + ε^{2} r^{2}}}{(L^{2} - s^{2}) + \sqrt{L^{2} + η^{2} ε^{2} r^{2}}}) - log (1) | \leq C ε^{2},

we obtain Eq. (14). □

Using Eq. (14), we replace the first term in the integrand of Eq. (13) to obtain (4). We can compare the expression (4) to that of Tornberg and Shelley,⁵⁵ where a regularization of the Keller–Rubinow model is used to obtain a second-kind integral equation for f. One thing to note is that, due to the form of the local term in our model (4), the effect of the regularization parameter η is the same in all directions (both tangent and normal to the fiber centerline). This is not necessarily the case for the Tornberg and Shelley model (see Sec. III A 3 for a spectral comparison given a simplified fiber geometry).

A. Spectral comparison of slender body integral operators

In this subsection we provide evidence that our model (4) is well suited for approximating the map $\frac{\partial X}{\partial t} \mapsto f$ needed to simulate the motion of a rigid fiber. Here we consider the spectrum of the integral operator taking the force density f to the fiber velocity $\frac{\partial X}{\partial t}$ in the nonphysical but nevertheless instructive case of a straight, periodic fiber with constant radius ε. In this scenario we can explicitly calculate the eigenvalues of both the slender body PDE operator (10) and the integral operator (13) and related models. This allows us to directly compare the properties of different models in the same simple setting and serves as a starting point for understanding more complicated geometries. In particular, we expect this analysis to roughly capture the high wavenumber behavior of these models in different geometries—on length scales much smaller than the variation in curvature and fiber radius. The high wavenumber behavior is of particular interest for the invertibility and stability of the slender body theory integral operator.

For comparison, we first recall the form of the eigenvalues of the slender body PDE (10), calculated in Ref. 34. In Sec. III A 2, we consider the model (13), before regularization, and show that the integral operator is negative definite. We compare the spectrum of Eq. (13) to three other possible models based on the slender body theory which do not result in negative definite operators. Then in Sec. III A 3, we consider the regularized version of our model (4) and compare its spectrum to the regularized model of Tornberg and Shelley.⁵⁵ We note that in our model, a uniform regularization parameter appears to give the best approximation of the slender body PDE spectrum in directions both normal and tangent to the slender body centerline, whereas in the Tornberg–Shelley model, the parameter required by the tangential direction may not be optimal in the normal direction.

1. Spectrum of the slender body PDE

Here we consider a straight, periodic fiber with constant radius ε. We take the fiber centerline to be two-periodic and lie along the z-axis, $X (z) = z e_{z}, z \in ℝ / 2 ℤ$ ⁠, and for simplicity take μ = 1 and zero background flow. We consider the stationary setting and omit the time dependence in our notation; in particular, we denote the fiber velocity by $\bar{u} (z)$ to distinguish from the fluid velocity away from the fiber.

We consider this scenario because we can explicitly calculate the eigenvalues of the slender body PDE (10) as well as various possible integral expressions for approximating the map $f \mapsto \bar{u}$ ⁠. In particular, the eigenvectors of this map can be decomposed into tangential (⁠ $e_{z}$ ⁠) and normal (⁠ $e_{x}, e_{y}$ ⁠) directions and are given by $f_{m} (z) = e^{iπkz} e_{m}, m = x, y, z$ ⁠. We may then explicitly solve for $λ_{k}^{m}$ satisfying

\bar{u} (z) = λ_{k}^{m} f_{m} (z), m = x, y, z,

(16)

for both the slender body PDE operator and various approximations based on the slender body theory. To avoid logarithmic growth of the corresponding bulk velocity field at spatial infinity, we will ignore translational modes (k = 0) in the following spectral analysis. Clearly these modes are important, especially for a rigid body; however, we are mainly interested in the high wavenumber behavior of these operators. High wavenumber instabilities are a known issue for nonlocal slender body theory,^18,45,55 and the following analysis likely captures the behavior of these models at high wavenumbers (small length scales) even in curved geometries.

To begin, the eigenvalues of the slender body PDE operator (10) mapping f to $\bar{u}$ were calculated in Ref. 34 (Proposition 1.4). Note that the sign convention in this paper is opposite, as we are considering f to be the hydrodynamic force exerted by rather than on the slender body. For the slender body PDE, the eigenvalues satisfying (16) in the tangential and normal directions, respectively, are given by

λ_{k}^{m} = {\begin{matrix} - \frac{2 K_{0} K_{1} + π ε | k | (K_{0}^{2} - K_{1}^{2})}{4 π^{2} ε | k | K_{1}^{2}}, & m = z \\ - \frac{2 K_{0} K_{1} K_{2} + π ε | k | (K_{1}^{2} (K_{0} + K_{2}) - 2 K_{0}^{2} K_{2})}{2 π^{2} ε | k | (4 K_{1}^{2} K_{2} + π ε | k | K_{1} (K_{1}^{2} - K_{0} K_{2}))}, & m = x, y, \end{matrix}

(17)

where each $K_{j} = K_{j} (π ε | k |)$ ⁠, j = 0, 1, 2, is a $j th$ order modified Bessel function of the second kind. Note that both sets of eigenvalues $λ_{k}^{z}$ and $λ_{k}^{x}, λ_{k}^{y}$ are strictly negative and decay to 0 at a rate proportional to $1 / | k |$ as $| k | \to \infty$ ⁠. We will compare our approximation and various other slender body approximations to Eq. (17).

2. Pre-regularization comparison

Before we consider the regularized version (4) of our model, we consider the base model (13) and compare its spectrum to other existing models based on slender body theory, before regularization. In the straight-but-periodic scenario, our model (13) becomes the periodization of the expression

\bar{u} (z) = - \frac{1}{8 π} \int_{- 1}^{1} (\frac{I}{{({\bar{z}}^{2} + ε^{2})}^{1 / 2}} + \frac{{\bar{z}}^{2} e_{z} e_{z}^{T}}{{({\bar{z}}^{2} + ε^{2})}^{3 / 2}} + \frac{ε^{2}}{2} (\frac{I}{{({\bar{z}}^{2} + ε^{2})}^{3 / 2}} - 3 \frac{{\bar{z}}^{2} e_{z} e_{z}^{T}}{{({\bar{z}}^{2} + ε^{2})}^{5 / 2}})) f (z - \bar{z}) d \bar{z} .

(18)

For this geometry, we may calculate the eigenvalues $λ_{k}^{m}$ satisfying (16), which are given by

λ_{k}^{m} = {\begin{matrix} - \frac{1}{8 π} \int_{- 1}^{1} \frac{2 {\bar{z}}^{4} + 2 ε^{2} {\bar{z}}^{2} + \frac{3}{2} ε^{4}}{{({\bar{z}}^{2} + ε^{2})}^{5 / 2}} e^{- i π k \bar{z}} d \bar{z}, & m = z \\ - \frac{1}{8 π} \int_{- 1}^{1} \frac{{\bar{z}}^{2} + \frac{3}{2} ε^{2}}{{({\bar{z}}^{2} + ε^{2})}^{3 / 2}} e^{- i π k \bar{z}} d \bar{z}, & m = x, y . \end{matrix}

(19)

These integrals may be computed explicitly to obtain

λ_{k}^{m} = {\begin{matrix} - \frac{1}{8 π} ((4 + π^{2} ε^{2} k^{2}) K_{0} (π ε | k |) - 2 π ε | k | K_{1} (π ε | k |)), & m = z \\ - \frac{1}{8 π} (2 K_{0} (π ε | k |) + π ε | k | K_{1} (π ε | k |)), & m = x, y . \end{matrix}

(20)

Here K₀ and K₁ are zero and first order modified Bessel functions of the second kind, respectively. The eigenvalues $λ_{k}^{m}$ lie along the curves plotted in Fig. 1. Importantly, these eigenvalues satisfy the following lemma.

FIG. 1.

View large Download slide

Log-scale plot of the tangential (a) and normal (b) eigenvalues $λ_{k}^{m}$ of the operator mapping $f \mapsto \bar{u}$ in various slender body models for a straight-but-periodic fiber. Our model (blue) results in strictly negative eigenvalues in both the tangential and normal directions, as does the slender body PDE (dotted). The Keller–Rubinow approximation (green) exhibits instabilities at wavenumbers $| k | \approx 0.2 / ε$ (tangential direction) and $| k | \approx 0.6 / ε$ (normal direction) as the eigenvalues of the operator mapping $f \mapsto \bar{u}$ become positive. For the modified Lighthill models, the normal direction eigenvalues $λ_{k}^{x}$ and $λ_{k}^{y}$ (red) remain negative at high wavenumber, but in the tangential direction, the eigenvalues of Modified Lighthill 1 (red) become positive when $| k | > 0.5 / ε$ ⁠. Furthermore, the tangential eigenvalues of Modified Lighthill 2 (magenta) do not agree with the slender body PDE at low wavenumber.

Lemma III.1. For all $| k | \geq 1$ and $m = x, y, z$ , the eigenvalues $λ_{k}^{m}$ given by Eq. (20) satisfy $λ_{k}^{m} < 0$ ⁠.

Proof. The case $m = x, y$ is immediate, since $K_{0} (t) > 0$ and $K_{1} (t) > 0$ for any t > 0.

For the tangential direction m = z, we first note that, by Lemma 1.16 in Ref. 34, we have

1 \leq \frac{K_{1} (t)}{K_{0} (t)} \leq 1 + \frac{1}{2 t},

for all t > 0. Letting $g (t) = (4 + t^{2}) K_{0} (t) - 2 t K_{1} (t)$ ⁠, it suffices to show that $g (t) / K_{0} (t) > 0$ ⁠. But

\frac{g (t)}{K_{0} (t)} = 4 + t^{2} - 2 t \frac{K_{1} (t)}{K_{0} (t)} \geq 3 + t^{2} - 2 t > {(t - \sqrt{3})}^{2} \geq 0.

□

Now, at a continuous level, regularization is necessary to make sense of inverting the integral operator (18), since K₀ and K₁ decay exponentially as $| k | \to \infty$ ⁠. However, at a discrete level, numerical approximation of Eq. (18) will be invertible, albeit with a large condition number, due to Lemma III.1. This negativity does not hold for other popular slender body approximations which rely on further asymptotic expansion of Eq. (13) with respect to ε to obtain a limiting centerline velocity expression. In particular, we consider the models of Keller and Rubinow²⁵ and of Lighthill.²⁹

The Keller–Rubinow model, proposed in Ref. 25 and further studied by Refs. 18, 24, 45, and 55, is equivalent to a full matched asymptotic expansion of Eq. (12) about ε = 0. In the straight-but-periodic setting, the Keller–Rubinow expression for the slender body velocity is given by

8 π \bar{u} (z) = - ((I - 3 e_{z} e_{z}^{T}) - 2 log (π ε / 8) (I + e_{z} e_{z}^{T})) f (z) - (I + e_{z} e_{z}^{T}) \frac{π}{2} \int_{- 1}^{1} \frac{f (z - \bar{z}) - f (z)}{| sin (π \bar{z} / 2) |} d \bar{z} .

(21)

The eigenvalues of the periodic Keller–Rubinow operator taking f to $\bar{u}$ have been calculated in Refs. 18, 45, and 55 and are given by

λ_{k}^{m} = {\begin{matrix} \frac{1}{4 π} (1 + 2 log (π ε | k | / 2) + 2 γ), & m = z \\ - \frac{1}{8 π} (1 - 2 log (π ε | k | / 2) - 2 γ), & m = x, y . \end{matrix}

(22)

Here $γ \approx 0.5772$ is the Euler gamma.

In both the tangent and normal directions, however, the Keller–Rubinow approximation runs into stability issues at moderately high wavenumbers, apparent in Fig. 1 at $| k | = \frac{2 e^{- γ - 1 / 2}}{π ε} \approx 0.217 / ε$ (tangent) and $| k | = \frac{2 e^{- γ + 1 / 2}}{π ε} \approx 0.589 / ε$ (normal). In particular, the curve containing the eigenvalues $λ_{k}^{m}$ crosses zero and becomes negative. This is an issue both because the slender body PDE eigenvalues (17) are strictly negative, and because, for arbitrary ε, there is no clear way to guarantee that $λ_{k}^{m} \neq 0$ ⁠, especially for more complicated fiber geometries. Thus some sort of regularization of Eq. (21) is necessary before approximating the inverse map $\bar{u} \mapsto f$ ⁠.

In addition to the Keller–Rubinow model, we consider what we will term the modified Lighthill approach to deriving a fiber velocity approximation. This approach, due to Lighthill,²⁹ also begins with the classical SBT expression (12) but uses asymptotic integration of the doublet term to arrive at an expression for the fiber velocity. We explore the Lighthill method in detail in Appendix A, but plot the resulting spectrum in Fig. 1.

The takeaway here is that, at least in the case of a straight, periodic fiber, our model (13), before regularization, captures the negative-definiteness of the slender body PDE and provides a better approximation than other models based on classical SBT.

3. Regularized comparison

To make our model truly suitable for inversion, we need to regularize the integral kernel as in Eq. (4). In the straight-but-periodic setting, the operator in Eq. (4) becomes the periodization of

\begin{matrix} 8 π \bar{u} (z) = - 2 log (η) f (z) - \int_{- 1}^{1} (\frac{I}{{({\bar{z}}^{2} + η^{2} ε^{2})}^{1 / 2}} + \frac{{\bar{z}}^{2} e_{z} e_{z}^{T}}{{({\bar{z}}^{2} + ε^{2})}^{3 / 2}} \\ + \frac{ε^{2}}{2} (\frac{I}{{({\bar{z}}^{2} + ε^{2})}^{3 / 2}} - 3 \frac{{\bar{z}}^{2} e_{z} e_{z}^{T}}{{({\bar{z}}^{2} + ε^{2})}^{5 / 2}})) f (z - \bar{z}) d \bar{z} . \end{matrix}

(23)

The eigenvalues of Eq. (23) are then given by

λ_{k}^{m} = {\begin{matrix} - \frac{1}{8 π} (2 log (η) + 2 K_{0} (η π ε | k |) + (2 + π^{2} ε^{2} k^{2}) K_{0} (π ε | k |) - 2 π ε | k | K_{1} (π ε | k |)), m = z \\ - \frac{1}{8 π} (2 log (η) + 2 K_{0} (η π ε | k |) + π ε | k | K_{1} (π ε | k |)), m = x, y . \end{matrix}

(24)

For $η > 1$ ⁠, the spectrum of our operator is bounded away from 0 and Eq. (23) is a second-kind integral equation for f.

We can compare the behavior of Eq. (23) with the Tornberg–Shelley regularization of the Keller–Rubinow model. In Refs. 45 and 55, the high wavenumber instability in Eq. (21) is removed by replacing the denominator of the integral term, which vanishes at $\bar{z} = 0$ ⁠, with an expression proportional to ε at $\bar{z} = 0$ ⁠. Using the relation

\int_{- 1}^{1} (\frac{π}{| 2 sin (π z / 2) |} - \frac{1}{| z |}) d z = - 2 log (π / 4),

(25)

to rewrite Eq. (21), a regularization $δ ε, δ > 0$ ⁠, is added to the denominator to obtain

8 π \bar{u} (z) = - ((I - 3 e_{z} e_{z}^{T}) + 2 log (δ) (I + e_{z} e_{z}^{T})) f (z) - (I + e_{z} e_{z}^{T}) \int_{- 1}^{1} \frac{f (z - \bar{z})}{{({\bar{z}}^{2} + δ^{2} ε^{2})}^{1 / 2}} d \bar{z} .

(26)

Here we have also used that the second term in the original Keller–Rubinow integral expression can now be integrated up to $O (ε^{2})$ errors to nearly cancel the logarithmic term in Eq. (21), leaving only $log (δ)$ ⁠. The idea is to then choose δ such that all eigenvalues of the operator taking $f \mapsto \bar{u}$ are negative. Since the integral kernel is now smooth, Eq. (26) is now a second-kind integral equation for f.

The eigenvalues of this δ-regularized Keller–Rubinow operator are given by

\begin{matrix} λ_{k}^{m} & = {\begin{matrix} - \frac{1}{4 π} (- 1 + 2 log δ + 2 K_{0} (δ π ε | k |)), & m = z \\ - \frac{1}{8 π} (1 + 2 log δ + 2 K_{0} (δ π ε | k |)), & m = x, y . \end{matrix} \end{matrix}

(27)

Since K₀ is positive, $λ_{k}^{z}$ is guaranteed to be negative and bounded away from 0 as long as $δ > \sqrt{e}$ (see Fig. 2).

FIG. 2.

View large Download slide

Log-scale plot of the tangential (a) and normal (b) eigenvalues $λ_{k}^{m}$ of our regularized model (23) (blue) with $η = 1.5$ and the Tornberg–Shelley δ-regularized model (26) (red) with $δ = \sqrt{e} + 0.5$ ⁠. Note that the regularization parameter η in our model affects the tangential and normal eigenvalues in a similar way; in particular, $η > 1$ is required in both cases to ensure that Eq. (23) is a second-kind integral equation. In the δ-regularized model, the tangential direction requires $δ > \sqrt{e}$ ⁠, but the normal direction does not, resulting at least visually in a greater disparity between the $λ_{k}^{x}, λ_{k}^{y}$ for the PDE (dotted) and the δ-regularized approximation.

Note that in our model (23), the regularization parameter η affects the spectrum of the operator mapping f to $\bar{u}$ in the same way in both the tangential and normal directions. In particular, in both directions, $η > 1$ is required to obtain the desired second-kind integral equation. In the Tornberg–Shelley model, the bound $δ > \sqrt{e} \approx 1.649$ is required to ensure negativity of the tangential eigenvalues, but this lower bound does not apply to the normal direction; in fact, $δ > e^{- 1} \approx 0.368$ is sufficient for ensuring strictly negative normal eigenvalues. This may mean that the η-regularization in our model is more physically reasonable; see Fig. 2.

In Ref. 34, it is shown that using the δ-regularized model (26) to approximate the map $\bar{u} \mapsto f$ yields $ε^{2}$ convergence to the slender body PDE for sufficiently smooth $\bar{u}$ ⁠. It is also shown that the constant in the resulting error estimate has the form $C_{1} δ^{2} (1 + log (δ)) + C_{2} / (- 1 + log (δ))$ for constants C₁ and C₂. We expect that a similar error estimate and analogous η dependence hold for our model (23); i.e., the constant should look like $C_{1} η^{2} + C_{2} / log (η)$ ⁠. If $C_{1} \approx C_{2}$ ⁠, this yields an optimal η of approximately 1.5. This should give a rough guideline for a good choice of η for more general curved geometries, at least in the periodic setting.

Remark III.2. We can also consider using the method of regularized Stokeslets to rederive the Keller–Rubinow model (see Ref. 13 ). Here the following choices of blob functions are used in place of Dirac deltas to derive the regularized Stokeslet and doublet, respectively:

ϕ_{S} (R) = \frac{15}{8 π} \frac{δ^{4} ε^{4}}{{(| R |^{2} + δ^{2} ε^{2})}^{7 / 2}}, ϕ_{D} (R) = \frac{3}{4 π} \frac{δ^{2} ε^{2}}{{(| R |^{2} + δ^{2} ε^{2})}^{5 / 2}} .

Note that we have modified the notation fromRef.13 to emphasize that the blob “width” will be taken to be proportional to the fiber radius ε, and to more easily compare with the δ-regularization of Tornberg–Shelley. For the straight-but-periodic fiber, this method yields a nearly identical expression to Eq. (26), but with a different logarithmic factor in front of the local terms: $- log (\sqrt{δ^{2} + 1} / δ)$ in place of $log (δ)$ . Due to the low wavenumber expansion (A6) of the Bessel function K₀, however, we note that the $log (δ)$ term in Eq. (27) exactly cancels the leading order dependence of $K_{0} (δ π ε | k |)$ on δ, yielding an expression consistent with the slender body PDE (17) when $| k |$ is small. When $δ ≪ 1$ , we have $- log (\sqrt{δ^{2} + 1} / δ) \approx log (δ)$ , but recall that $δ > \sqrt{e}$ is required for Eq. (27) to be negative for all k. Thus this particular choice of blob function in the method of regularized Stokeslets appears to yield an expression for the fiber velocity which fundamentally differs from the slender body PDE solution, although a different choice of blob function may yield closer agreement. Note that this low wavenumber descrepancy occurs whether we start from the non-periodic or periodic regularized expressions mentioned inRef.13 , due to the identity (25).

IV. NUMERICAL DISCRETIZATION OF THE SLENDER BODY MODEL

We turn now to numerically simulating thin rigid fibers in flows. We begin by generally discussing the numerical solution of Fredholm integral equations where the result must be integrated (i.e., to find the total force and torque on a rigid fiber). We apply these general methods to the slender body model (4) and perform convergence tests. We note improvements in conditioning and stability for the second kind (⁠ $η > 1$ ⁠) vs first kind (η = 1) integral equation. Finally, we look at the spectrum of the discretized integral operator in different geometries to verify the negative definite nature of the operator.

A. Solving the second-kind Fredholm integral equation

Denote by $K : L^{2} ([- L, L], ℝ^{3}) \to L^{2} ([- L, L], ℝ^{3})$ the integral operator

K [f] (s) : = \int_{- L}^{L} K (s, s') f (s') d s' .

(28)

Then a Fredholm integral equation of the first kind reads

y (s) = K [f] (s) .

(29)

It is well known that the inversion of such an integral operator is an ill-posed problem, meaning that the solution may not be unique or not even exist.^2,22,26 Furthermore, small perturbations to the left hand side of Eq. (29) can lead to relatively large perturbations of the solution $f (s)$ ⁠. The ill-posedness of this problem can be circumvented by regularizing the integral operator into a second-kind Fredholm integral equation, which takes the form

y (s) = (α I + K) [f] (s),

(30)

for some parameter α. Discretization of Eq. (30) yields a linear system with a far better condition number. The connection between Eq. (30) and our model is illustrated in Sec. IV B.

Numerical methods for solving Fredholm integral equations are well documented^22,57 and the approach we adopt is based on the Nyström method (see Ref. 2, Chapter 12.4). For rigid fibers, after numerically inverting a second-kind Fredholm integral equation, linear functionals (8) will also need to be applied to the resulting $f (s)$ to find the total force and torque.

We consider the numerical approximation of a general linear functional of $f (s)$ ⁠, given by

ϕ_{M} (f) = \int_{- L}^{L} M (s) f (s) d s .

(31)

Here $M (s) \in ℝ^{3 \times 3}$ is a bounded, smooth operator and $f (s)$ is found by numerically inverting a second-kind Fredholm integral equation of the form (30). The numerical method is obtained discretizing the Eq. (30) by replacing the integral with a convergent quadrature formula with nodes $- L = s_{1} < s_{2} < \dots < s_{n} = L$ and weights $w = {(w_{1}, w_{2}, \dots, w_{n})}^{T} \in ℝ^{n}$ ⁠, and requiring the numerical approximation $f_{i}^{[n]} \approx f (s_{i})$ to satisfy

y (s_{i}) = α f_{i}^{[n]} + \sum_{j = 1}^{n} w_{j} K (s_{i}, s_{j}) f_{j}^{[n]} for i = 1, \dots, n .

(32)

Introducing the vectors ${\underline{f}}^{[n]} = {({(f_{1}^{[n]})}^{T}, \dots, {(f_{n}^{[n]})}^{T})}^{T}$ and $\underline{y} = {(y {(s_{1})}^{T}, \dots, y {(s_{n})}^{T})}^{T}$ ⁠, Eq. (32) can be written compactly as

\underline{y} = (α I + \underline{K} \underline{W}) {\underline{f}}^{[n]} .

(33)

Here I denotes the $3 n \times 3 n$ identity matrix, and

\underline{W} = diag (w) \otimes I, \underline{K} = (\begin{matrix} K (s_{1}, s_{1}) & \dots & K (s_{1}, s_{n}) \\ ⋮ & ⋱ & ⋮ \\ K (s_{n}, s_{1}) & \dots & K (s_{n}, s_{n}) \end{matrix}) \in ℝ^{3 n \times 3 n}

(34)

with $\otimes : ℝ^{n_{1} \times m_{1}} \times ℝ^{n_{2} \times m_{2}} \to ℝ^{(n_{1} n_{2}) \times (m_{1} m_{2})}$ the Kronecker product of matrices and $I$ the 3 × 3 identity matrix. We then approximate (31) by the same quadrature formula

ϕ_{M} (f) \approx \sum_{i = 1}^{n} w_{i} M (s_{i}) f_{i}^{[n]} = ({\vec{1}}^{T} \otimes I) \underline{M} \underline{W} {\underline{f}}^{[n]} : = ϕ_{M}^{[n]},

(35)

where

\underline{M} = (\begin{matrix} M (s_{1}) & 0 \\ ⋱ \\ 0 & M (s_{n}) \end{matrix}) \in ℝ^{3 n \times 3 n}

(36)

and $\vec{1} = {(1, \dots, 1)}^{T} \in ℝ^{n} .$ Here we have used $ϕ_{M}^{[n]}$ to denote the approximation of $ϕ_{M} (f)$ obtained by quadrature. After inserting the solution of Eq. (33), we obtain

ϕ_{M}^{[n]} = ({\vec{1}}^{T} \otimes I) \underline{M} \underline{W} {(α I + \underline{K} \underline{W})}^{- 1} \underline{y} .

(37)

Remark IV.1.The numerical approximation $ϕ_{M}^{[n]}$ shares the same convergence as the underlying quadrature method. This is illustrated in Appendix B.

B. Application to the slender body model and convergence tests

We apply the numerical method from Sec. IV A to approximate the force and torque on a slender body. Note that the Eq. (8) is given by setting $M (s) = I$ and $M (s) = \hat{X} (s)$ in the functional (31). That is,

F = ϕ_{I} (f) and T = ϕ_{\hat{X}} (f) .

(38)

Letting $α = 2 log (η)$ and

K (s, s') = S_{ε, η} (s, s') + \frac{ε^{2} r^{2} (s')}{2} D_{ε} (s, s'),

(39)

y (s) = - 8 π μ (v - \hat{X} (s) ω - u_{0} (X (s, t), t)),

(40)

our model (4) is of the form (30), and we may write the discretization of Eq. (4) in the form (33). Here we have introduced the hat operator $\hat{.} : ℝ^{3} \to s o (3)$ which maps vectors in $ℝ^{3}$ to 3 × 3 skew symmetric matrices by

ω = (\begin{matrix} ω_{1} \\ ω_{2} \\ ω_{3} \end{matrix}) \mapsto \hat{ω} = (\begin{matrix} 0 & - ω_{3} & ω_{2} \\ ω_{3} & 0 & - ω_{1} \\ - ω_{2} & ω_{1} & 0 \end{matrix}) .

(41)

Here, $s o (3)$ is the Lie algebra of SO(3), and such that $ω \times v = \hat{ω} v$ for $ω, v \in ℝ^{3}$ ⁠.

Denote the numerical approximations to Eq. (38) by

F^{[n]} = ϕ_{I}^{[n]} and T^{[n]} = ϕ_{\hat{X}}^{[n]} .

(42)

Defining the matrices $Φ$ and $Ψ \in ℝ^{3 \times 3 n}$ as

Φ = ({\vec{1}}^{T} \otimes I) \underline{W} {(α I + \underline{K} \underline{W})}^{- 1},

(43)

Ψ = ({\vec{1}}^{T} \otimes I) \underline{W} \underline{X} {(α I + \underline{K} \underline{W})}^{- 1},

(44)

we may then write Eq. (42) as

F^{[n]} = Φ \underline{y} and T^{[n]} = Ψ \underline{y} .

(45)

In Secs. IV B 1 and IV B 2 we perform convergence tests for our discrete model (45) for both a thin ring and a prolate spheroid. With these geometries we are able to calculate accurate reference solutions against which we can compare the accuracy of our numerical solution. Furthermore, we will look at how the conditioning of the linear system associated with the discretized integral operator improves as the regularization parameter η is increased from η = 1 to $η > 1$ ⁠.

Remark IV.2. For very large aspect ratios, e.g., $L / ε \approx O (10^{3})$ or larger, the kernel becomes very nearly singular meaning one must take n very large to accurately resolve the $O (ε)$ length scales in the kernel. In this case, the quadrature can be improved by implementing special quadrature methods that take into account the near singular nature of the integral kernel.^1,53 For modest aspect ratios, e.g., $L / ε \approx O (10^{2})$ , this is not an issue as one can accurately resolve the kernel with a few hundred points. As noted in Ref. 46 , even local slender body theory, i.e., just the leading order fiber velocity approximation $8 π u (s) = 2 log (ε) (I - e_{s} e_{s}^{T}) f (s)$ , yields “reasonable predictions” for the behavior of particles with aspect ratio larger than 20. We expect that the integral model (4) should be more physically realistic than the local approximation, and in some of the following numerical tests we consider aspect ratios down to about 20.

1. Thin ring translating with unit velocity

As a convergence test, we use Eq. (45) to calculate the force on a thin ring of unit length in the xy-plane translating in the z direction with unit velocity in zero background flow. We will consider both the first- and second-kind formulations of the model. In this setting, the force on the ring can be calculated to arbitrary high precision by evaluating elliptic integrals, which can be used as a reference solution. For a circular centerline parametrized by

X (s) = {(\frac{cos (π s)}{2 π}, \frac{sin (π s)}{2 π}, 0)}^{T},

the z-component of our unregularized (η = 1) model becomes

8 π μ = - \int_{- \frac{1}{2}}^{\frac{1}{2}} \frac{\sqrt{2} π (3 ε^{2} π^{2} - cos (2 π (s - s')) + 1)}{{(2 ε^{2} π^{2} - cos (2 π (s - s')) + 1)}^{3 / 2}} f^{z} (s') d s' .

(46)

As in the straight-but-periodic geometry of Sec. III A, the eigenfunctions of this operator are the Fourier modes $f_{k}^{z} (s) = exp (i 2 π k s)$ ⁠. The force $F = {(F, 0, 0)}^{T}$ is therefore given by

F = \int_{- \frac{1}{2}}^{\frac{1}{2}} f^{z} (s) d s = \frac{8 π μ}{λ_{0}^{z}},

(47)

where $λ_{0}^{z}$ is the k = 0 eigenvalue. This can be found by evaluating the integral in Eq. (46) with $f^{z} (s) = f_{0}^{z} (s) = 1$ ⁠, which gives

λ_{0}^{z} = - c_{ε} (2 ϕ_{K} (c_{ε}) + ϕ_{E} (c_{ε})) .

(48)

Here $c_{ε} = \sqrt{{(ε^{2} π^{2} + 1)}^{- 1}}$ ⁠, and

ϕ_{K} (x) = \int_{0}^{1} \frac{1}{\sqrt{1 - θ^{2}} \sqrt{1 - x^{2} θ^{2}}} d θ, ϕ_{E} (x) = \int_{0}^{1} \frac{\sqrt{1 - x^{2} θ^{2}}}{\sqrt{1 - θ^{2}}} d θ

(49)

are the complete elliptic integrals of the first and second kind, respectively.

For $ε = 0.05, 0.025, 0.01$ and 0.005, Eq. (46) is discretized using trapezoidal quadrature, and we numerically approximate F by Eq. (45). Figure 3 plots the numerical error as a function of n for four different values of ε. We observe spectral convergence of the numerical error to machine precision, which is consistent with the error estimates (B21). We note that the condition number of the unregularized discrete integral operator grows exponentially as n increases, as shown in Fig. 4(a). However, because we are considering a rigid fiber with constant radius, computing F has a regularizing effect which lessens the impact of this ill-conditioning in the final force calculation. This may be contrasted with the prolate spheroid, where, as we will see in Sec. IV B 2, the conditioning does have a noticeable effect on the error. Nevertheless, we note that by setting $η > 1$ we can improve the condition number of the linear system [see Fig. 4(b)]. We also note that there is a $1 / ε$ dependence on n for a given accuracy. This can be circumvented by using a special quadrature method that takes into account the kernel (see Remark IV.2).

FIG. 3.

View large Download slide

The approximate drag force $F^{[n]}$ on a thin ring translating broadwise with unit velocity converges with spectral accuracy to the true force F.

FIG. 4.

View large Download slide

The condition numbers associated with the discretized versions of the (a) unregularized (η = 1) and (b) regularized (⁠ $η = 1.5$ ⁠) slender body models for calculating the force on a thin ring. Note the change in scale between the two figures.

2. Prolate spheroid with artificial fluid velocity field

We next use Eq. (45) to compute the drag force for a stationary prolate spheroid immersed in an artificial fluid velocity field. The particle centerline is aligned in the z-direction, parameterized by $X (s) = {(0, 0, s)}^{T}, s \in [- 1, 1]$ ⁠. The fluid velocity field $u (s) = {(u (s), 0, 0)}^{T}$ is designed such that $f (s) = {(f^{x} (s), 0, 0)}^{T}$ is a known analytic function. We choose this function to be a Gaussian $f^{x} (s) = exp (- \frac{s^{2}}{ε^{2}})$ such that the force decays to zero at the fiber endpoints and use high order Gauss–Lobatto quadrature for the discretization of the integral operator. Denote the set of n quadrature nodes by ${s_{i}}_{i = 1}^{n}$ ⁠. Inserting the above expression for $f^{x} (s)$ into our model (18), the fluid velocity at s_i is found by solving the integral

u (s_{i}) = \frac{- 1}{8 π} (2 log (η) exp (- \frac{s_{i}^{2}}{ε^{2}}) + \int_{- 1}^{1} \frac{ε^{2} r {(s_{i})}^{2} + \frac{1}{2} ε^{2} r {(s^{'})}^{2} + {(s_{i} - s^{'})}^{2}}{{(ε^{2} r {(s_{i})}^{2} + {(s_{i} - s^{'})}^{2})}^{3 / 2}} exp (- \frac{s^{' 2}}{ε^{2}}) d s'),

(50)

where the ellipsoidal radius function is given by Eq. (2). We also take the viscosity μ = 1. To solve for $u (s_{i})$ for $i = 1, \dots, n$ ⁠, the integral in Eq. (50) is evaluated to machine precision using MATLAB's built-in integral function, which uses adaptive quadrature. For this fluid velocity field, the total force $F = {(F, 0, 0)}^{T}$ on the ellipsoid is found by

F = \int_{- 1}^{1} exp (- \frac{s^{2}}{ε^{2}}) d s = \sqrt{π} ε erf (\frac{1}{ε}) .

(51)

We compute numerical approximations to F using Eq. (45) for four choices of ε. We initially set η = 1 and compute these numerical approximations for the non-regularized, first-kind equation. The numerical errors are presented in Fig. 5(a). We see that the error converges spectrally up to a certain point where the method begins to diverge due to numerical instabilities and poor conditioning of the discrete integral operator, which is plotted in Fig. 5(b).

FIG. 5.

View large Download slide

The errors (a) and condition numbers (b) associated with the unregularized (η = 1) numerical method for the calculation of the force on a prolate spheroid for different values of ε.

However, by choosing $η > 1$ ⁠, we can amend the condition number and therefore improve the accuracy of the numerical solution. In Fig. 6, we fix $ε = 0.025$ and calculate the numerical errors for four choices of η. We see from Fig. 6(a) that the numerical error converges spectrally to machine precision for all such choices of η. Furthermore, we observe from Fig. 6(b) that the condition number of the discrete integral operator is bounded by a value that becomes smaller for larger η. We note that in practice, the modeling error is much larger than machine precision as we will see in Sec. V B 2.

FIG. 6.

View large Download slide

The errors (a) and condition numbers (b) associated with the regularized numerical method for the calculation of the force on a prolate spheroid for $ε = 0.025$ ⁠. Similar results are observed for other values of ε.

C. Spectrum of the slender body operator in different geometries

One important unresolved question about the slender body model (4) is the effect of different geometries, including curvature, endpoints, and non-uniform fiber radius, on the spectrum of the integral operator. The main difficulty is that the integral kernel (5),(6) is only well defined along the centerline of the fiber. Since the kernel is so dependent on the shape of the fiber centerline, it is difficult to prove general properties for it. Although we cannot analytically determine the spectrum of the continuous operator in general, we can determine the eigenvalues of the discrete operator $(2 log (η) I + \underline{K} \underline{W})$ (33). We consider first the unregularized version η = 1, recalling that in the straight-but-periodic geometry of Sec. III A, the continuous operator was provably negative definite. Ideally we would like to see evidence that this negative definiteness persists in general geometries, as this would be the physically correct behavior and also would agree with the underlying slender body PDE operator (10).

We begin by calculating the eigenvalues ${λ_{i}}_{i = 1}^{3 n}$ of $\underline{K} \underline{W}$ for the thin ring. Letting $λ_{max} = {max}_{i} (λ_{i})$ ⁠, in Fig. 7(a) we plot $λ_{max}$ vs n for five different values of ε. Note that for very large n relative to $ε^{- 1}$ [roughly $n = O (ε^{- 2})$ ], we begin to see numerical error resulting in very small positive eigenvalues of $\underline{K} \underline{W}$ (denoted by red markers). However, the magnitude of these positive eigenvalues is on the order of machine precision and may be attributed to round-off errors.

FIG. 7.

View large Download slide

Magnitude of the maximum eigenvalue of the non-regularized discrete slender body operator $\underline{K} \underline{W}$ for the thin ring (a) and the prolate spheroid (b). Blue markers mean $λ_{max} < 0$ while red markers mean that $λ_{max} > 0$ ⁠.

We next consider the effects of endpoints and a non-uniform radius by calculating the eigenvalues of $\underline{K} \underline{W}$ for a slender prolate spheroid (2), keeping in mind the above level of numerical error. In Fig. 7(b) we again plot $λ_{max}$ vs n for four different values of ε. Again for $n = O (ε^{- 2})$ we begin to see small positive eigenvalues which are significantly larger than for the thin ring [around $O (10^{- 10})$ ]. However, the magnitude of the positive eigenvalues is still very small and bounded as n increases. It is not clear whether this is a numerical artifact or an actual eigenvalue crossing 0 for the continuous operator. At any rate, the non-regularized operator would never actually be used for simulations with such large n because the condition number of $\underline{K} \underline{W}$ is prohibitive [see Fig. 5(b)]. It appears that a very reasonable choice of regularization parameter η will ensure that none of these near-zero eigenvalues actually cross zero.

As a final test, we calculate the spectrum of $\underline{K} \underline{W}$ for randomly but systematically generated curvy fibers with complicated shapes (Fig. 8). Here the magnitude of the fiber's deviation from a straight line is controlled by a small parameter $δ \geq 0$ ⁠. The fiber shapes are generated by interpolating m points $(x_{i}, y_{i}, z_{i}) \in ℝ^{3}, i = 1, \dots, m$ ⁠, with cubic splines. Here $z_{i} = (i - 1) \frac{2 L}{m}$ while $x_{i}, y_{i} \in [- δ, δ]$ are given by a random number generator and are of size at most δ. Setting δ = 0 corresponds to a straight fiber. Examples of the fiber centerline for m = 10 and four different values of δ are given in Fig. 8.

FIG. 8.

View large Download slide

The centerlines of four curved fiber shapes with curviness parameterized by (a) $δ \approx 3 * 10^{- 4}$ ⁠, (b) $δ \approx 2 * 10^{- 3}$ ⁠, (c) $δ \approx 9 * 10^{- 3}$ ⁠, and (d) $δ = 5 * 10^{- 2}$ ⁠.

We fix $ε = 0.1$ and use the spheroidal radius function (2). Taking m = 10, we generate 6 different curvy fibers for different magnitudes $δ \in [0, \frac{1}{10}]$ ⁠. For each fiber we compute the spectrum ${λ_{i}^{δ}}_{i = 1}^{n}$ of its corresponding (non-regularized) integral operator $\underline{K} \underline{W}$ ⁠. We plot the most positive eigenvalue $λ_{max}^{δ} = {max}_{i} (λ_{i}^{δ})$ for each fiber in Fig. 9(a). For each value of δ we note that there is an eigenvalue crossing zero when $n = O (ε^{- 2})$ ⁠. As δ increases and the magnitude of the curviness of the fiber increases, we can note a slight increase in the magnitude of the largest positive eigenvalue, but $λ_{max}^{δ}$ is still small—roughly $O (10^{- 8})$ ⁠. Again, we can be assured to have a negative spectrum bounded away from 0 by a reasonable choice of regularization $η > 1$ ⁠. This effect is displayed in Fig. 9(b), which shows the maximum eigenvalue $λ_{max}^{δ, η}$ of the now regularized discrete integral operator $(2 log (η) I + \underline{K} \underline{W})$ for a fixed value of ε and δ and varying values of η. We see here that for all choices of $η > 1$ in this range, the spectrum of $(2 log (η) I + \underline{K} \underline{W})$ remains negative definite.

FIG. 9.

View large Download slide

The magnitude $λ_{max}^{δ, η}$ of the maximum eigenvalues for the unregularized (a) and regularized (b) discrete integral operators for the curved fibers. For (b) we fix $ε = 0.1$ and $δ = 0.001$ and consider different regularizations η. The color blue denotes a negative maximum eigenvalue and red denotes a positive maximum eigenvalue.

V. DYNAMICS OF CURVED RIGID FIBERS

We next use the slender body model (4) and the discretization procedure of Sec. IV to simulate the dynamics of curved rigid fibers in Stokes flow. After outlining the dynamical equations, we validate the model against known dynamical models for a slender prolate spheroid. Finally, we compare the rotational dynamics of randomly curved fibers as in Fig. 8 to straight fibers.

A. Dynamical equations

Assuming that the particle to fluid density ratio is large $ρ_{p} / ρ_{f} ≫ 1$ ⁠, such as in gas-solid fiber suspensions,^15,27,30,38 the dynamics of the slender body are governed by the following rigid body equations. The angular momentum m of a rigid particle with torque $T (t)$ is found by solving

\dot{m} = m \times ω + T,

(52)

where $ω = J^{- 1} m$ for moment of inertia tensor J. Each of these quantities is defined in a reference frame whose axes are co-rotating and co-translating with the fiber. The fiber orientation (with respect to a fixed inertial reference frame) is specified using Euler parameters $q \in ℝ^{4}$ which satisfy the constraint $| | q | |_{2} = 1$ and are determined by solving the ordinary differential equation (ODE)

\dot{q} = \frac{1}{2} q w,

(53)

where $w = {(0, ω^{T})}^{T} \in ℝ^{4}$ ⁠. Here, $q w$ is the Hamilton product of two quaternions.¹⁷ That is, by letting $q = (q_{0}, q)$ and $r = (r_{0}, r)$ denote quaternions for $q_{0}, r_{0} \in ℝ$ and $q, r \in ℝ^{3}$ ⁠, then their Hamilton product is given by

q r = (q_{0} r_{0} - q \cdot r, q_{0} r + r_{0} q + q \times r) .

(54)

The translational dynamics are given by Newton's second law

\dot{p} = F,

(55)

where $p = v m$ is the inertial frame linear momentum for a fiber of mass m. The position of the fiber center of mass is found by solving

\dot{x} = v .

(56)

The ODEs (52)–(56) are integrated using the second order Strang splitting method of Ref. 50.

Recall Eq. (45) for $F^{[n]}$ and $T^{[n]}$ ⁠. Since $F^{[n]}$ and $T^{[n]}$ depend linearly on the linear and angular momenta p and m, we may update them according to the linear equation

(\begin{matrix} F^{[n]} \\ T^{[n]} \end{matrix}) = A (\begin{matrix} p \\ m \end{matrix}) + b,

(57)

where A is a negative definite dissipation matrix and $b$ is due to the background fluid velocity and is independent of p and m. We have that

A = (\begin{matrix} Φ (\vec{1} \otimes (I / m)), & Φ (- \underline{X} (\vec{1} \otimes J^{- 1})) \\ Ψ (\vec{1} \otimes (I / m)), & Ψ (- \underline{X} (\vec{1} \otimes J^{- 1})) \end{matrix}) and b = - (\begin{matrix} Φ \underline{u} \\ Ψ \underline{u} \end{matrix}),

(58)

where m and J are the filament mass and moment of inertia tensor, respectively. We have also introduced the vector $\underline{u} = {(u_{0} {(X (s_{1}))}^{T}, \dots, u_{0} {(X (s_{n}))}^{T})}^{T}$ containing the background fluid velocities at the location of the quadrature nodes along the centerline.

1. Overview and cost of algorithm

The algorithm used to compute the dynamics of a slender fiber is as follows:

Define particle geometry $X (s)$ ⁠, ε, regularization parameter η, and discretization n.
Choose a quadrature rule and compute the matrices $\underline{W}$ and $\underline{K}$ ⁠.
Compute the matrices $Φ, Ψ$ ⁠, and A from Eqs. (43), (44), and (58).
Time loop: for $t = 0, Δ t, \dots, m Δ t$

Compute $F^{[n]}$ and $T^{[n]}$ using Eq. (57)
Numerically integrate the ODEs (52)–(56).

For step (2), we use the trapezoidal quadrature rule for closed fibers (i.e., a periodic integration interval) or Gauss–Lobatto quadrature rule for fibers with open ends. For step (4b), we use a splitting method.⁵⁰ We note that for simulations where the fluid velocity field is calculated from a direct numerical simulation of the Navier–Stokes equations, the fluid field needs to be approximated onto the centerline of the particle using an interpolation method.⁵¹

The above algorithm exploits the rigidity of the fiber by using the fact that A, $Φ$ ⁠, and $Ψ$ are constant in time and therefore can be computed outside of the time loop. The calculation of these matrices, which involves solving a linear system, is the most costly operation in the algorithm but only needs to be done once. If, for example, Gaussian elimination is used, this step has complexity of $O (n^{3})$ ⁠. Within the time loop, however, the most costly operation is the calculation of $F^{[n]}$ and $T^{[n]}$ ⁠, which involves only $3 \times 3 n$ by $3 n \times 1$ matrix-vector products, which has O(n) complexity. We assume that the cost of numerically integrating the ODEs is negligible compared to this. For a single fiber, the total complexity of the algorithm is therefore $O (n^{3} + n m)$ ⁠, where m is the total number of time steps used in the simulation. Hence, for simulations where many time steps are needed, the algorithm scales by O(n). We remark that for problems where the background flow is zero, the cost of computing $F^{[n]}$ and $T^{[n]}$ is independent of n (after A has been computed) and therefore is O(1). This is relevant, for example, when simulating fibers sedimenting in a still fluid under the influence of gravity.³⁷

B. Numerical validation of model dynamics

1. Dissipation matrix of a prolate spheroid

Here we compare our model and numerical method with accurate closed form expressions for the force and torque given by Brenner⁶ and Jeffery.²³ These expressions are valid for an ellipsoid when the fluid Jacobian is approximately constant throughout the volume of the particle. When the flow is linear, these terms are essentially exact and therefore serve as a good reference model against which to validate our model.

The purpose of this numerical experiment is therefore twofold. First, we aim to show that our model converges to the reference model as $ε \to 0$ ⁠. This is primarily to validate the accuracy of the model. However, the numerical approximation of the force and torques also introduces a numerical error that is related to the discretization parameter n. Clearly, taking n too small means that we will not exploit the accuracy of the model to its entirety. On the other hand, it is unwise to take n as large as possible as this will incur unnecessary computational costs that go to minimizing numerical error beyond the accuracy of the model. So the second question we address here is what is an ideal choice of discretization parameter to use such that the numerical error is roughly the same as the modeling error.

Using $η = 1 + ε^{2}$ ⁠, the dissipation matrix for our slender body model A is numerically approximated by Eq. (58). The reference dissipation matrix A_sph is found using the closed form expressions from Jeffery and Brenner, which are given in Appendix C. Denote the six eigenvalues of A and A_sph by λ_i and $λ_{i}^{sph}$ ⁠, respectively. Note that due to symmetry of the spheroid, $λ_{1} = λ_{2}$ and $λ_{4} = λ_{5}$ and similarly for the eigenvalues of A_sph. Furthermore, the slender body model is essentially a one dimensional filament and therefore $λ_{6} = 0$ meaning that spinning motion about the centerline does not dissipate. This is in contrast to the Jeffrey term, which does dissipate spinning motion. We remark that this phenomenon only occurs in the case where the centerline is perfectly straight. Hence for curved fiber geometries where the application of the slender body is most useful, this nonphysical phenomenon is not observed. Note that for this geometry the dissipation matrices are diagonal and therefore the eigenvalues are directly proportional to the calculation of $F^{[n]}$ and $T^{[n]}$ in zero background flow.

The eigenvalues of A are calculated using Eq. (58) after discretizing Eq. (30) on the Gauss–Lobatto nodes. The values $| λ_{i} - λ_{i}^{sph} |$ for i = 1, 3, 4 are plotted in Fig. 10 as a function of the discretization parameter n. We see that λ_i converges exponentially to a point near $λ_{i}^{sph}$ ⁠, which is likely due to the slender body modeling error. As ε decreases, we make two observations. First, for large n the rate at which λ_i converges to $λ_{i}^{sph}$ is approximately $- ε^{2} η^{2} log (ε η)$ ⁠, as seen by the horizontal dash-dot lines. Second, as ε decreases, the convergence rate slows down and one must use a larger value of n to reach the most accurate solution. This means that one must pay careful attention to the choice of n when taking ε to be very small. In fact, we observe empirically that the convergence rate is approximately bounded by $e^{- 4 ε n}$ ⁠. Motivated by this, we will take n in future experiments to be approximately the intersection of these two lines, that is,

n \approx - \frac{log (- ε^{2} η^{2} log (ε η))}{4 ε} .

(59)

FIG. 10.

View large Download slide

The difference in the dissipation matrix eigenvalues $| λ_{i} - λ_{i}^{sph} |$ ⁠, i = 1, 3, 4 as a function of n for four different values of ε: (a) $ε = 0.02$ ⁠, (b) $ε = 0.01$ ⁠, (c) $ε = 0.005$ ⁠, and (d) $ε = 0.0025$ ⁠. The black dashed lines are $e^{- 4 ε n}$ and the horizontal dashed-dotted lines are $- ε^{2} η^{2} log (ε η)$ ⁠.

2. Prolate spheroids rotating in shear flow

Now we calculate the dynamics of a prolate spheroid in shear flow $u = {(z, 0, 0)}^{T}$ using our model and compare it with that of the accurate Jeffrey model. The fiber is initially aligned at rest in the z-direction and its rotational dynamics are calculated by integrating Eq. (52) on the interval $t \in [0, 100]$ using the splitting method of Ref. 50 with a small step size of h = 0.01. The simulation was repeated with h = 0.05 with no significant changes to the results and it is therefore concluded that time integration errors are negligible. We repeat the experiment for 20 values of ε logarithmically spaced in the interval $[0.1, 0.001]$ and choose n using Eq. (59) and $η = 1 + ε^{2}$ ⁠. As the spheroids are axisymmetric, they only experience a torque about their y axis, hence all of other angular momentum components are zero (to machine precision). Three examples of the rotational dynamics are shown in Fig. 11. It is seen here that as ε becomes smaller, the dynamics more closely resemble the Jeffery model.

FIG. 11.

View large Download slide

The y component of a spheroid rotating in shear flow for three different values of ε: (a) $ε = 0.1$ ⁠, (b) $ε = 0.048329$ ⁠, and (c) $ε = 0.01833$ ⁠. The solid line is our slender body model and the dashed line is due to Jeffery.

The relative difference between the angular momenta of the Jeffery and slender body solutions is calculated and averaged over the simulation. This average relative error is then plotted against the corresponding value of ε in Fig. 12. We see that the average relative error decreases with ε. It is observed that in the region $0.01 < ε < 0.1$ the error converges at a faster rate than in the region $0.001 < ε < 0.01$ ⁠. This could be partially explained by the fact that wider particles (larger ε) experience a greater resistive force as seen by the regions where m_y nearly reaches zero. This means that the particle spends more time in the shear plane where the fluid velocity is zero and hence the slender body model does not experience a large torque. However, the fluid gradient is non-zero in this orientation and therefore the Jeffery model, which depends only on the fluid gradient, still experiences a constant torque. This means that compared to the Jeffery model, thicker fibers will see a greater difference in the torque term when the fiber is aligned in the shear plane than thinner fibers.

FIG. 12.

View large Download slide

The relative difference in m_y between the slender body and Jeffery solutions averaged over the interval $[0, 100]$ ⁠.

C. Dynamics of randomly curvy fibers

Understanding how different shaped particles rotate in shear flow is an important step in understanding their dynamics in more complex flows.⁵² Here we simulate the dynamics of the randomly curvy fibers of Fig. 8 as they rotate in shear flow. In particular, we show how the rotational variables deviate from a straight fiber as δ becomes larger.

We generate 100 different fiber shapes with m = 10 using 10 different values of δ logarithmically spaced in the interval $[5 \times 10^{- 5}, 5 \times 10^{- 2}]$ ⁠. The 100 fibers are placed in shear flow $u = {(z, 0, 0)}^{T}$ and their rotational dynamics are calculated on the interval $t \in [0, 100]$ ⁠. The moment of inertia tensor is approximated by placing point masses along the centerline and using the formula

J_{i, i} = \sum_{j = 1}^{k} m_{j} {(X_{i} (s_{j}) - c_{i})}^{2}, for i = 1, \dots, 3,

(60)

where $X_{i} (s_{j})$ is the ith component of the centerline function at the point s_j on the centerline and c_i is the ith component of the fiber center of mass. We weight m_j by the cross-sectional radius and use a very large value for k, e.g., $k = 10^{4}$ ⁠. Here we take $ε = 0.01$ and use the spheroidal radius function (2) along with $η = 1 + ε^{2}$ ⁠.

Figure 13(a) shows the angular momentum m of three fibers compared to the δ = 0 case. As the δ = 0 fiber is perfectly straight, it does not exhibit spinning motion and its angular momentum is purely in the m_y component. This is in contrast to the fibers with a non-zero value of δ, in which case some of the momentum is transferred to m_x. We therefore compare the value $\sqrt{m_{x}^{2} + m_{y}^{2}}$ between the fibers to account for this. We see here that the $δ = 0.017783$ solution is visually very similar to the δ = 0 solution. We notice a significant difference between the other two solutions. Figure 13(b) shows the angle θ between the z-axis of the particle reference frame (that is, a frame that is rotating with the fiber) and the x-axis of a fixed inertial reference frame. As the $δ \neq 0$ fibers are not symmetric, they slowly rotate out of the xz-plane and therefore after a long time, we see much more significant discrepancies in θ.

FIG. 13.

View large Download slide

The rotational variables of four fibers with different values of δ. (a) shows the angular momentum and (b) is the angle between the fiber's long axis and the x-axis of the inertial frame.

To quantify the effect that δ has on the angular momentum, we calculate the difference in the angular momentum $Δ m$ by subtracting off the δ = 0 solution and averaging over the time interval $t \in [92, 100]$ ⁠, which corresponds to roughly one period of rotation. This value is averaged over all the fibers with similar values of δ and is expressed as a percentage of the δ = 0 solution, which we denote by $% Δ m$ ⁠. The results are plotted in Fig. 14(a). We notice that $% Δ m$ is linearly proportional to δ. We observe that at the end of the simulation the $δ = 0.0003$ fibers correspond to roughly 1% discrepancy in angular momentum and $δ = 0.0015$ corresponds to roughly 7.5% discrepancy.

FIG. 14.

View large Download slide

(a) shows the difference in angular momentum $Δ m$ between the curved fibers and the δ = 0 solution after 100 time units and averaged over all the fibers with similar δ. The black dashed line is $O (δ)$ ⁠. (b) shows the discrepancy $Δ θ$ in the angle between the centerline and the x-axis after roughly one rotation.

The difference in θ after one rotation as a function of δ is displayed in Fig. 14(b). The $δ = 0.0003$ solution corresponds to about a $3 °$ difference in θ and the $δ = 0.0015$ solution corresponds to about an $8 °$ difference.

VI. CONCLUSIONS

We have developed an integral model for the motion of a thin filament in a viscous fluid based on the nonlocal slender body theory. The model relies on standard singular Stokeslets and doublets but makes use of the fiber integrity condition—the near-cancelation of angular-dependent terms along the fiber surface—in a novel way to yield an integral expression for the fiber velocity with a smooth kernel which retains dependence on the (possibly varying) fiber radius in a natural way. We include a systematic way of comparing mapping properties of different models using the simplified geometry of a straight-but-periodic filament. In this simple geometry, we can show that our integral operator is negative definite and compares favorably to other models, and we expect similar high wavenumber behavior for curved filaments with constant radius. It is less clear how a non-constant radius affects the spectrum; however, numerical tests indicate that the discretized integral operator is very close to negative definite. Nevertheless, to ensure invertibility, we develop an asymptotically consistent regularization to convert the first-kind Fredholm integral equation for the force density along the fiber into a second-kind equation and show that this second-kind regularization improves the stability and conditioning of the discretized equation. We numerically solve the integral equation using the Nyström method² and show how constraining the fiber motion to be rigid can be exploited for fast computation of fiber dynamics. We validate the method and model against the prolate spheroid model of Jeffery,²³ and apply the method to study the rotational deviation of randomly curved rigid fibers from straight fibers.

While the fibers considered here are rigid, the model can also be used to simulate the dynamics of semiflexible filaments. The invertibility properties of the integral equation make it particularly well suited for handling simulations involving inextensible fibers, where an additional line tension equation must be solved at each time step.^32,55 We may also consider the effects of different choices of radius functions on the model properties, similar to what is done in Ref. 58, although we note the necessity of smooth decay in our radius function near the fiber endpoints.

To build on the dynamic simulations for rigid fibers, we aim to consider the effects of fiber shape on particle deposition and aggregation. We are especially interested in more complicated background flows, including suspensions of rigid fibers in turbulence. The novel modeling approach advocated herein will enable earlier explorations based on the point-particle approach⁷ to be extended to curved fibers particles.

ACKNOWLEDGMENTS

This work has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement (No. 691070) as well as the SPIRIT project (No. 231632) under the Research Council of Norway FRIPRO funding scheme. E. Celledoni, B. Owren, and B. K. Tapley would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Geometry, compatibility, and structure preservation in computational differential equations (2019) where part of the work on this paper was undertaken. E. Celledoni and B. Owren also acknowledge funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skodowska-Curie grant agreement No. 860124. L. Ohm was supported by a University of Minnesota Doctoral Dissertation Fellowship and NSF Postdoctoral Research Fellowship Grant No. 2001959. The authors are listed in alphabetical order.

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

Appendix A: Modified Lighthill model

Here we consider the modified Lighthill approach to deriving a fiber velocity approximation from classical SBT (12). This approach takes advantage of the fact that the doublet term of Eq. (13) only has an O(1) contribution to the fiber velocity very close to $s' = s$ ⁠, and thus can be integrated asymptotically to leave only a local term. This results in a model similar to that of Lighthill,²⁹ which was derived via different reasoning but also includes a local doublet term and a nonlocal Stokeslet contribution (see Remark A.1).

There are two ways to consider the nonlocal Stokeslet contribution. The first expression, which we will term Modified Lighthill 1, is given by the periodization of

\bar{u} (z) = - \frac{1}{8 π} ((I - e_{z} e_{z}^{T}) f (s) + \int_{- 1}^{1} (\frac{I}{{({\bar{z}}^{2} + ε^{2})}^{1 / 2}} + \frac{{\bar{z}}^{2} e_{z} e_{z}^{T}}{{({\bar{z}}^{2} + ε^{2})}^{3 / 2}}) f (z - \bar{z}) d \bar{z}) .

(A1)

Here the local term $(I - e_{z} e_{z}^{T})$ comes from asymptotically integrating the doublet term of Eq. (12) (see estimate 3.65 of Ref. 35 for more details). Note that in Eq. (A1), the Stokeslet term inside the integral is equal to $f / ε$ when $\bar{z} = 0$ ⁠.

For the second expression, which we will call Modified Lighthill 2, the $e_{z} e_{z}^{T}$ component of the Stokeslet term is normalized to give the same order contribution at $\bar{z} = 0$ as in Eq. (12); namely, $(I + e_{z} e_{z}^{T}) f / ε$ ⁠. This yields the periodization of the expression

\bar{u} (z) = - \frac{1}{8 π} ((I - e_{z} e_{z}^{T}) f (s) + \int_{- 1}^{1} \frac{I + e_{z} e_{z}^{T}}{{({\bar{z}}^{2} + ε^{2})}^{1 / 2}} f (z - \bar{z}) d \bar{z}) .

(A2)

Remark A.1. The actual model proposed by Lighthill in Ref. 29 , written in the periodic, straight setting, has the form

\bar{u} (z) = - \frac{1}{8 π} (2 (I - e_{z} e_{z}^{T}) f (z) + \int_{| \bar{z} | > q} \frac{I + e_{z} e_{z}^{T}}{| \bar{z} |} f (z - \bar{z}) d \bar{z}); q = ε \sqrt{e} / 2.

(A3)

At first glance, this looks like a slightly different model from Eqs. (A1) and (A2), due to the 2 in front of the $(I - e_{z} e_{z}^{T}) f (z)$ term. However, the extra factor here is precisely due to the removal of the section $| \bar{z} | \leq q$ from the integral term. Indeed, if we consider the integrand of Eq. (A1), we note that

\begin{matrix} \int_{- q}^{q} (\frac{I}{{({\bar{z}}^{2} + ε^{2})}^{1 / 2}} + \frac{{\bar{z}}^{2} e_{z} e_{z}^{T}}{{({\bar{z}}^{2} + ε^{2})}^{3 / 2}}) f (z - \bar{z}) d \bar{z} = (2 log (2 q / ε) (I + e_{z} e_{z}^{T}) - 2 e_{z} e_{z}^{T}) f (z) + O (ε^{2} / q^{2}) \\ = (I - e_{z} e_{z}^{T}) f (z) + O (ε^{2} / q^{2}) \end{matrix}

for q as in Eq. (A3). Now, this particular choice of q is not large relative to ε, so the $O (ε^{2} / q^{2})$ error term is not small asymptotically. However, this is merely a heuristic and we will not be considering the expression (A3) in greater depth here. Furthermore, the expressions (A1) and (A2) are more amenable to calculating eigenvalues.

The eigenvalues of Eq. (A1) are given by

\begin{matrix} λ_{k}^{m} & = {\begin{matrix} - \frac{1}{4 π} (2 K_{0} (π ε | k |) - π ε | k | K_{1} (π ε | k |)), & m = z \\ - \frac{1}{8 π} (1 + 2 K_{0} (π ε | k |)), & m = x, y . \end{matrix} \end{matrix}

(A4)

Now the normal eigenvalues $λ_{k}^{x}$ and $λ_{k}^{y}$ are always negative. However, there is still a high wavenumber instability in the tangent direction. In particular, $λ_{k}^{z} = 0$ when $π ε | k | \approx 1.55265$ ⁠, and becomes positive at higher wavenumbers (see Fig. 1). Thus the instability issue is not fully resolved by expanding only the doublet term of Eq. (12).

For Modified Lighthill 2, the eigenvalues of Eq. (A2) are given by

\begin{matrix} λ_{k}^{m} & = {\begin{matrix} - \frac{1}{2 π} K_{0} (π ε | k |), & m = z \\ - \frac{1}{8 π} (1 + 2 K_{0} (π ε | k |)), & m = x, y . \end{matrix} \end{matrix}

(A5)

Here the eigenvalues $λ_{k}^{x}$ and $λ_{k}^{y}$ in the normal directions are identical to Eq. (A4), but the tangential eigenvalues $λ_{k}^{z}$ are very different. In fact, they are too different: Recall that near t = 0, the modified Bessel functions $K_{0} (t)$ and $K_{1} (t)$ satisfy

K_{0} (t) = - log (t / 2) - γ + O (t^{2}); t K_{1} (t) = 1 + O (t^{2}) .

(A6)

Therefore, at low wavenumber (⁠ $k = O (1)$ ⁠), the tangential eigenvalues of Modified Lighthill 2 (A2) look like

λ_{k}^{z} = \frac{1}{2 π} (log (π ε | k | / 2) + γ) + O (ε^{2} k^{2}) .

This does not agree with the low wavenumber behavior of the slender body PDE (17) (see Fig. 1). It appears that the normalization in Modified Lighthill 2 (A2) results in the wrong model.

For the sake of completeness, we also consider a modification of our model (13) in which the $\bar{X} {\bar{X}}^{T}$ terms are normalized as in Modified Lighthill 2 (A2) to yield a nonzero contribution to the fiber velocity when $s = s'$ ⁠. In the case of the periodic straight centerline, the modified version of our model becomes the periodization of

\bar{u} (z) = - \frac{1}{8 π} \int_{- 1}^{1} (\frac{I + e_{z} e_{z}^{T}}{{({\bar{z}}^{2} + ε^{2})}^{1 / 2}} + \frac{ε^{2}}{2} \frac{I - 3 e_{z} e_{z}^{T}}{{({\bar{z}}^{2} + ε^{2})}^{3 / 2}}) f (z - \bar{z}) d \bar{z} .

(A7)

The eigenvalues of Eq. (A7) are given by

\begin{matrix} λ_{k}^{m} & = {\begin{matrix} - \frac{1}{4 π} (2 K_{0} (π ε | k |) - π ε | k | K_{1} (π ε | k |)), & m = z \\ - \frac{1}{8 π} (2 K_{0} (π ε | k |) + π ε | k | K_{1} (π ε | k |)), & m = x, y . \end{matrix} \end{matrix}

(A8)

Now, the eigenvalues $λ_{k}^{x}$ and $λ_{k}^{y}$ in the directions normal to the fiber are unchanged from our original expression (20). However, the tangent eigenvalues $λ_{k}^{z}$ are now given by the same expression as Modified Lighthill 1 (A4), which we recall exhibits a high wavenumber instability (Fig. 1).

Appendix B: Convergence and error bounds of numerical method

We are interested in obtaining an estimate for the error when approximating (31) by its discrete approximation (37), which we denote by

d^{[n]} = ϕ_{M} (f) - ϕ_{M}^{[n]} = \int_{- L}^{L} M (s) f (s) d s - \sum_{j = 1}^{n} w_{j} M (s_{j}) f_{j}^{[n]} .

(B1)

This error will depend on the error committed in the numerical approximation of Eq. (30) by the solution ${\underline{f}}^{[n]}$ of Eq. (33). For this reason, we first analyze the convergence of Nyström's method (see Ref. 2, Chapter 12.4) in using Eq. (33) to approximate the solution of Eq. (30). At each quadrature node, we define the error of this approximation as

e_{i}^{[n]} : = f (s_{i}) - f_{i}^{[n]}, for i = 1, \dots, n,

(B2)

and let ${\underline{e}}^{[n]} : = {({(e_{1}^{[n]})}^{T}, \dots, {(e_{n}^{[n]})}^{T})}^{T}$ denote the error vector. We want to show that $| | {\underline{e}}^{[n]} | |_{\infty} \to 0$ as $n \to \infty$ ⁠. Let $\underline{f} : = {(f {(s_{1})}^{T}, \dots, f {(s_{n})}^{T})}^{T}$ and define ${\underline{τ}}^{[n]} : = {(τ_{1}^{T}, \dots, τ_{n}^{T})}^{T}$ with components

τ_{i} : = y (s_{i}) - α f (s_{i}) - \sum_{j}^{n} K_{i, j} w_{j} f (s_{j}),

(B3)

the truncation error for the discrete second kind Eq. (33)—i.e., the residual obtained replacing ${\underline{f}}^{[n]}$ by $\underline{f}$ in Eq. (33). We obtain

(α I + \underline{K} \underline{W}) \underline{f} = \underline{y} - {\underline{τ}}^{[n]} .

(B4)

It is easily seen using Eq. (30) that

τ_{i} = \int_{- L}^{L} K (s_{i}, s') f (s') d s' - \sum_{j}^{n} K_{i, j} w_{j} f (s_{j}),

(B5)

which is simply quadrature error, and for any convergent quadrature formula we have

\lim_{n \to \infty} | | τ^{[n]} | |_{\infty} = 0.

(B6)

We next bound the norm of the error ${\underline{e}}^{[n]}$ by the norm of ${\underline{τ}}^{[n]}$ to prove the convergence of the method. Subtracting Eq. (33) from Eq. (B4) we obtain a linear system satisfied by ${\underline{e}}^{[n]}$ ⁠,

(α I + \underline{K} \underline{W}) {\underline{e}}^{[n]} = - {\underline{τ}}^{[n]} .

(B7)

From Ref. 2 [Chapter 12.4, Theorem 12.4.4 and Equation (12.4.51)], we have that for sufficiently large n, say $n \geq n^{*}$ ⁠, the matrix $(α I + \underline{K} \underline{W})$ is invertible and

| | {(α I + \underline{K} \underline{W})}^{- 1} | |_{\infty} \leq C_{1} \forall n \geq n^{*} .

(B8)

Thus we can conclude that

| | {\underline{e}}^{[n]} | |_{\infty} \leq | | {(α I + \underline{K} \underline{W})}^{- 1} | |_{\infty} | | {\underline{τ}}^{[n]} | |_{\infty} \leq C_{1} | | {\underline{τ}}^{[n]} | |_{\infty} .

(B9)

Since C₁ is independent of n for $n \geq n^{*}$ and $| | {\underline{τ}}^{[n]} | |_{\infty} \to 0$ as $n \to \infty$ ⁠, this implies that

\lim_{n \to \infty} | | {\underline{e}}^{[n]} | |_{\infty} = 0.

Consider now the quadrature error

δ^{[n]} : = \int_{- L}^{L} M (s) f (s) d s - \sum_{j = 1}^{n} w_{j} M (s_{j}) f (s_{j}) .

(B10)

From Eq. (B1) we obtain

d^{[n]} = δ^{[n]} - \sum_{j = 1}^{n} w_{j} M (s_{j}) e_{j},

(B11)

and using Eq. (B7) the total discretization error for our methods is given by

d^{[n]} = ({\vec{1}}^{T} \otimes I) \underline{W} \underline{M} {(α I + \underline{K} \underline{W})}^{- 1} {\underline{τ}}^{[n]} + δ^{[n]} .

(B12)

Since both $δ^{[n]}$ and $τ^{[n]}$ are quadrature errors, $| | {(α I + \underline{K} \underline{W})}^{- 1} | | \leq C_{1}$ for all $n \geq n^{*}$ ⁠, and M is bounded, the method converges at the same rate as the underlying quadrature.

1. Convergence of numerical method for closed loop geometry

By applying the formula (B12), we now show how one can achieve spectral convergence in the case of a closed fiber geometry with constant radius ε and periodic integration domain. In this setting, we will use trapezoidal quadrature. We begin by bounding the norms of the integration kernels to which we apply the trapezoidal quadrature rules to, namely the integrals (28) and (31). Using this, and some smoothness assumptions, we are able bound the quadrature errors $τ_{i}^{[n]}$ and $δ^{[n]}$ using classical error estimates. This leads to a bound on the total error $d^{[n]}$ for both the force and torque calculation.

Let C₂ be a constant such that

| | f (s') | |_{\infty} \leq C_{2} for s \in [- L, L] .

(B13)

From the definition of $K (s, s')$ [Eqs. (5), (6), and (39)] in the constant radius case, we observe that

| | K (s, s') | |_{\infty} \leq \frac{3}{2 ε},

(B14)

with equality when $s = s'$ ⁠. From Eq. (38) we have $| | M (s) | |_{\infty} = 1$ for the force calculation, while for the torque calculation, $M (s) = \hat{X} (s)$ and therefore

| | M (s) | |_{\infty} \leq max_{s \in [- L, L]} | | X (s) | |_{1} .

(B15)

Therefore, we can bound the integration kernels of Eqs. (28) and (31) by

| | K (s, s') f (s') | |_{\infty} \leq \frac{3}{2 ε} C_{2}

(B16)

and

| | M (s) f (s) | |_{\infty} \leq | | M (s) | |_{\infty} C_{2} .

(B17)

Note that in the constant radius case, $K (s, s')$ has the same regularity as $X (s)$ ⁠. If we assume that $X (s), f (s)$ ⁠, and M(s) are analytic, then using [56, Theorem 3.2] we can bound the trapezoidal rule quadrature error from Eq. (B3) by

| | τ_{i}^{[n]} | |_{\infty} \leq \frac{6 L C_{2}}{ε (e^{a n} - 1)} for i = 1, \dots, n .

(B18)

Similarly, we can bound Eq. (B10) by

| | δ^{[n]} | |_{\infty} \leq \frac{4 L | | M (s) | |_{\infty} C_{2}}{e^{a n} - 1} .

(B19)

Here a is some constant. Using Eq. (B12), the total discretization error is therefore bounded as

| | d^{[n]} | |_{\infty} \leq (| | ({\vec{1}}^{T} \otimes I) \underline{W} \underline{M} {(α I + \underline{K} \underline{W})}^{- 1} | |_{\infty} \frac{3}{2 ε} + | | M (s) | |_{\infty}) \frac{4 L C_{2}}{e^{a n} - 1} .

(B20)

Using that $| | \underline{M} | |_{\infty} \leq | | M (s) | |_{\infty}, | | \underline{W} | |_{\infty} = \frac{2 L}{n}$ and C₁ is given by Eq. (B8), this simplifies to

| | d^{[n]} | |_{\infty} \leq (\frac{6 C_{1} L}{2 ε} + 1) \frac{4 L | | M (s) | |_{\infty} C_{2}}{e^{a n} - 1} .

(B21)

Hence, the method shares the same exponential convergence as the underlying trapezoidal rule. We remark that one could perform an analogous analysis for open ended fiber geometries with, e.g., Gauss–Lobatto quadrature, and derive similar results. Furthermore, we also remark that one could require less stringent regularity assumptions on the integration on the kernels or the fiber centerline $X (s)$ ⁠, e.g., $M (s) f (s) \in C^{2 m + 2} [- L, L]$ ⁠. Then Ref. 3 (Theorem 5.5) can be used to derive asymptotic error estimates for $τ_{i}^{[n]}$ and $δ^{[n]}$ of order $O (h^{2 m + 2})$ ⁠. Nonetheless, we do observe spectral convergence in numerical experiments in the following sections, as predicted by the bound (B21).

Appendix C: Dissipation matrix of a prolate spheroid

The non-dimensionalized body frame resistance tensor R₁ for a spheroid with aspect ratio λ was derived by Oberbeck³⁹ and is given by

R_{1} = 16 π λ diag (\frac{1}{χ_{0} + α_{0}}, \frac{1}{χ_{0} + β_{0}}, \frac{1}{χ_{0} + λ^{2} γ_{0}}) .

(C1)

The constants χ₀, α₀, β₀, and γ₀ were calculated by Siewert⁴⁷ and are presented for a prolate (⁠ $λ > 1$ ⁠) spheroid

χ_{0} = \frac{- κ_{0} λ}{\sqrt{λ^{2} - 1}},

(C2)

α_{0} = β_{0} = \frac{λ^{2}}{λ^{2} - 1} + \frac{λ κ_{0}}{2 {(λ^{2} - 1)}^{3 / 2}},

(C3)

γ_{0} = \frac{- 2}{λ^{2} - 1} - \frac{λ κ_{0}}{{(λ^{2} - 1)}^{3 / 2}},

(C4)

κ_{0} = ln (\frac{λ - \sqrt{λ^{2} - 1}}{λ + \sqrt{λ^{2} - 1}}) .

(C5)

The torques $N = {(N_{x}, N_{y}, N_{z})}^{T}$ that describe the rotational forces acting on an ellipsoid in creeping Stokes flow in the body frame were calculated by Jeffery²³ and are presented in their non-dimensional form with zero background flow

N_{x} = - \frac{16 π λ}{3 (β_{0} + λ^{2} γ_{0})} [(1 + λ^{2}) ω_{x}],

(C6)

N_{y} = - \frac{16 π λ}{3 (α_{0} + λ^{2} γ_{0})} [(1 + λ^{2}) ω_{y}],

(C7)

N_{z} = - \frac{32 π λ}{3 (α_{0} + β_{0})} ω_{z} .

(C8)

Here $ω = {(ω_{x}, ω_{y}, ω_{z})}^{T}$ is the body frame angular velocity, which is related to body frame angular momentum by $m = J ω$ ⁠. Taking derivatives of N with respect to m gives for the rotational dissipation matrix

R_{2} = - \frac{16 λ}{3} diag (\frac{(1 + λ^{2})}{(β_{0} + λ^{2} γ_{0})}, \frac{(1 + λ^{2})}{(α_{0} + λ^{2} γ_{0})}, \frac{2}{(α_{0} + β_{0})}) J^{- 1} .

(C9)

The full dissipation matrix used for the calculation in Fig. 10 is given by

A_{sph} = (\begin{matrix} R_{1} & 0 \\ 0 & R_{2} \end{matrix}) .

(C10)

Appendix D: Endpoint behavior of model

Here we numerically determine the behavior at the fiber endpoints of the force density $f (s)$ that results from inverting the model (4). Although the end point behavior of the corresponding slender body PDE for the “inverse problem” is unknown, it is possible that decay in $f (s)$ is required to accurately approximate the PDE solution up to the fiber endpoints. We consider two different free-end geometries—the slender prolate spheroid and a cylinder with hemispherical caps, both translating with uniform unit velocity—and note that some decay in f is indeed observed at the endpoints of the filament (Fig. 15). Note that the prolate spheroid force density appears to be better behaved than the cylindrical fiber with hemispherical caps, which exhibits erroneous-looking oscillations toward the fiber ends. This may mean that the model (4) is better suited for modeling fibers whose radii decay more gradually toward the fiber endpoints.

FIG. 15.

View large Download slide

The computed force-per-unit-length $f (s)$ for the prolate spheroid (a), (b) and cylinder with hemispherical caps (c), (d) with centerline aligned with the x-axis. The left figures show the x-component of the force density for the spheroid (a) and cylinder (c) translating with unit speed in the x-direction, while the right figures show the y-component of the force density for the spheroid (b) and cylinder (d) translating with unit speed in the y-direction. Here we use the regularization $η = 1.1$ and take $n = 1 / ε$ discretization points.

References

1.

L.

af Klinteberg

and

A. H.

Barnett

, “

Accurate quadrature of nearly singular line integrals in two and three dimensions by singularity swapping

,”

BIT Numer. Math.

61

,

83

–

36

(

2021

).

https://doi.org/10.1007/s10543-020-00820-5

Google Scholar

Crossref

2.

K.

Atkinson

and

W.

Han

,

Theoretical Numerical Analysis

(

Springer

,

2005

), Vol.

39

.

Google Scholar

Crossref

3.

K. E.

Atkinson

,

An Introduction to Numerical Analysis

(Wiley,

1978

).

Google Scholar

4.

G.

Batchelor

, “

Slender-body theory for particles of arbitrary cross-section in Stokes flow

,”

J. Fluid Mech.

44

(

3

),

419

–

440

(

1970

).

https://doi.org/10.1017/S002211207000191X

Google Scholar

Crossref

5.

E. L.

Bouzarth

and

M. L.

Minion

, “

Modeling slender bodies with the method of regularized stokeslets

,”

J. Comput. Phys.

230

(

10

),

3929

–

3947

(

2011

).

https://doi.org/10.1016/j.jcp.2011.02.017

Google Scholar

Crossref

6.

H.

Brenner

, “

The stokes resistance of an arbitrary particle–iv arbitrary fields of flow

,”

Chem. Eng. Sci.

19

(

10

),

703

–

727

(

1964

).

https://doi.org/10.1016/0009-2509(64)85084-3

Google Scholar

Crossref

7.

N. R.

Challabotla

,

L.

Zhao

, and

H. I.

Andersson

, “

On fiber behavior in turbulent vertical channel flow

,”

Chem. Eng. Sci.

153

,

75

–

86

(

2016

).

https://doi.org/10.1016/j.ces.2016.07.002

Google Scholar

Crossref

8.

S.

Chattopadhyay

and

X.-L.

Wu

, “

The effect of long-range hydrodynamic interaction on the swimming of a single bacterium

,”

Biophys. J.

96

(

5

),

2023

–

2028

(

2009

).

https://doi.org/10.1016/j.bpj.2008.11.046

Google Scholar

Crossref

PubMed

9.

A. T.

Chwang

and

T. Y.-T.

Wu

, “

Hydromechanics of low-Reynolds-number flow. Part 2: Singularity method for Stokes flows

,”

J. Fluid Mech.

67

(

4

),

787

–

815

(

1975

).

https://doi.org/10.1017/S0022112075000614

Google Scholar

Crossref

10.

R.

Cortez

, “

The method of regularized stokeslets

,”

SIAM J. Sci. Comput.

23

(

4

),

1204

–

1225

(

2001

).

https://doi.org/10.1137/S106482750038146X

Google Scholar

Crossref

11.

R.

Cortez

, “

Regularized Stokeslet segments

,”

J. Comp. Phys.

375

,

783

–

796

(

2018

).

https://doi.org/10.1016/j.jcp.2018.08.055

Google Scholar

Crossref

12.

R.

Cortez

,

L.

Fauci

, and

A.

Medovikov

, “

The method of regularized Stokeslets in three dimensions: Analysis, validation, and application to helical swimming

,”

Phys. Fluids

17

(

3

),

031504

(

2005

).

https://doi.org/10.1063/1.1830486

Google Scholar

Crossref

13.

R.

Cortez

and

M.

Nicholas

, “

Slender body theory for Stokes flows with regularized forces

,”

Commun. Appl. Math. Comput. Sci.

7

(

1

),

33

–

62

(

2012

).

https://doi.org/10.2140/camcos.2012.7.33

Google Scholar

Crossref

14.

R.

Cox

, “

The motion of long slender bodies in a viscous fluid part 1. general theory

,”

J. Fluid Mech.

44

(

4

),

791

–

810

(

1970

).

https://doi.org/10.1017/S002211207000215X

Google Scholar

Crossref

15.

A. R.

Esmaeili

,

B.

Sajadi

, and

M.

Akbarzadeh

, “

Numerical simulation of ellipsoidal particles deposition in the human nasal cavity under cyclic inspiratory flow

,”

J. Braz. Soc. Mech. Sci. Eng.

42

(

5

),

1

–

13

(

2020

).

https://doi.org/10.1007/s40430-020-02345-y

Google Scholar

Crossref

16.

X.

Fan

,

N.

Phan-Thien

, and

R.

Zheng

, “

A direct simulation of fibre suspensions

,”

J. Non-Newton. Fluid Mech

74

(

1

),

113

–

135

(

1998

).

https://doi.org/10.1016/S0377-0257(97)00050-5

Google Scholar

Crossref

17.

H.

Goldstein

,

C.

Poole

, and

J.

Safko

,

Classical Mechanics

(Addison-Wesley, Reading, MA,

2002

).

Google Scholar

18.

T.

Götz

, “

Interactions of fibers and flow: Asymptotics, theory and numerics

,” Doctoral dissertation (

University of Kaiserslautern

,

2000

).

Google Scholar

19.

K.

Gustavsson

and

A.-K.

Tornberg

, “

Gravity induced sedimentation of slender fibers

,”

Phys. Fluids

21

(

12

),

123301

(

2009

).

https://doi.org/10.1063/1.3273091

Google Scholar

Crossref

20.

J.

Hämäläinen

,

S. B.

Lindström

,

T.

Hämäläinen

, and

H.

Niskanen

, “

Papermaking fibre-suspension flow simulations at multiple scales

,”

J. Eng. Math.

71

(

1

),

55

–

79

(

2011

).

https://doi.org/10.1007/s10665-010-9433-5

Google Scholar

Crossref

21.

G.

Hancock

, “

The self-propulsion of microscopic organisms through liquids

,”

Proc. R. Soc. Lond. A

217

(

1128

),

96

–

121

(

1953

).

https://doi.org/10.1098/rspa.1953.0048

Google Scholar

Crossref

22.

P. C.

Hansen

, “

Numerical tools for analysis and solution of fredholm integral equations of the first kind

,”

Inverse Probl.

8

(

6

),

849

(

1992

).

https://doi.org/10.1088/0266-5611/8/6/005

Google Scholar

Crossref

23.

G. B.

Jeffery

, “

The motion of ellipsoidal particles immersed in a viscous fluid

,”

Proc. R. Soc. Lond. A

102

(

715

),

161

–

179

(

1922

).

https://doi.org/10.1098/rspa.1922.0078

Google Scholar

Crossref

24.

R. E.

Johnson

, “

An improved slender-body theory for Stokes flow

,”

J. Fluid Mech.

99

(

02

),

411

–

431

(

1980

).

https://doi.org/10.1017/S0022112080000687

Google Scholar

Crossref

25.

J. B.

Keller

and

S. I.

Rubinow

, “

Slender-body theory for slow viscous flow

,”

J. Fluid Mech.

75

(

4

),

705

–

714

(

1976

).

https://doi.org/10.1017/S0022112076000475

Google Scholar

Crossref

26.

R.

Kress

,

V.

Maz'ya

, and

V.

Kozlov

,

Linear Integral Equations

(

Springer

,

1989

), Vol.

82

.

Google Scholar

Crossref

27.

S.

Kuperman

,

L.

Sabban

, and

R.

van Hout

, “

Inertial effects on the dynamics of rigid heavy fibers in isotropic turbulence

,”

Phys. Rev. Fluids

4

(

6

),

064301

(

2019

).

https://doi.org/10.1103/PhysRevFluids.4.064301

Google Scholar

Crossref

28.

E.

Lauga

and

T. R.

Powers

, “

The hydrodynamics of swimming microorganisms

,”

Rep. Progr. Phys.

72

(

9

),

096601

(

2009

).

https://doi.org/10.1088/0034-4885/72/9/096601

Google Scholar

Crossref

29.

J.

Lighthill

, “

Flagellar hydrodynamics

,”

SIAM Rev.

18

(

2

),

161

–

230

(

1976

).

https://doi.org/10.1137/1018040

Google Scholar

Crossref

30.

C.

Marchioli

,

M.

Fantoni

, and

A.

Soldati

, “

Orientation, distribution, and deposition of elongated, inertial fibers in turbulent channel flow

,”

Phys. Fluids

22

(

3

),

033301

(

2010

).

https://doi.org/10.1063/1.3328874

Google Scholar

Crossref

31.

J.

Martin

,

A.

Lusher

,

R. C.

Thompson

, and

A.

Morley

, “

The deposition and accumulation of microplastics in marine sediments and bottom water from the irish continental shelf

,”

Sci. Rep.

7

(

1

),

10772

(

2017

).

https://doi.org/10.1038/s41598-017-11079-2

Google Scholar

Crossref

PubMed

32.

O.

Maxian

,

A.

Mogilner

, and

A.

Donev

, “

Integral-based spectral method for inextensible slender fibers in stokes flow

,”

Phys. Rev. Fluids

6

(

1

),

014102

(

2021

).

https://doi.org/10.1103/PhysRevFluids.6.014102

Google Scholar

Crossref

33.

Y.

Mori

and

L.

Ohm

, “

An error bound for the slender body approximation of a thin, rigid fiber sedimenting in Stokes flow

,”

Res. Math. Sci.

7

(

8

),

8

(

2020

).

https://doi.org/10.1007/s40687-020-00206-7

Google Scholar

Crossref

34.

Y.

Mori

and

L.

Ohm

, “

Accuracy of slender body theory in approximating force exerted by thin fiber on viscous fluid

,”

Stud. Appl. Math.

(published online) (

2021

).

https://doi.org/10.1111/sapm.12380

Google Scholar

35.

Y.

Mori

,

L.

Ohm

, and

D.

Spirn

, “

Theoretical justification and error analysis for slender body theory

,”

Comm. Pure Appl. Math.

73

(

6

),

1245

–

1314

(

2020

).

https://doi.org/10.1002/cpa.21872

Google Scholar

Crossref

36.

Y.

Mori

,

L.

Ohm

, and

D.

Spirn

, “

Theoretical justification and error analysis for slender body theory with free ends

,”

Arch. Ration. Mech. Anal.

235

,

1905

–

1978

(

2020

).

https://doi.org/10.1007/s00205-019-01458-6

Google Scholar

Crossref

37.

R.

Newsom

and

C.

Bruce

, “

The dynamics of fibrous aerosols in a quiescent atmosphere

,”

Phys. Fluids

6

(

2

),

521

–

530

(

1994

).

https://doi.org/10.1063/1.868347

Google Scholar

Crossref

38.

D. O.

Njobuenwu

and

M.

Fairweather

, “

Simulation of inertial fibre orientation in turbulent flow

,”

Phys. Fluids

28

(

6

),

063307

(

2016

).

https://doi.org/10.1063/1.4954214

Google Scholar

Crossref

39.

A.

Oberbeck

, “

Uber stationare flussigkeitsbewegungen mit berucksichtigung der inner reibung

,”

J. Reine Angew. Math.

81

,

62

–

80

(

1876

).

https://doi.org/10.1515/crll.1876.81.62

Google Scholar

Crossref

40.

L.

Ohm

,

B. K.

Tapley

,

H. I.

Andersson

,

E.

Celledoni

, and

B.

Owren

, “

A slender body model for thin rigid fibers: Validation and comparisons

,” arXiv:1906.00253 (

2019

).

Google Scholar

41.

C. J.

Petrie

, “

The rheology of fibre suspensions

,”

J. Non-Newtonian Fluid Mech.

87

(

2

),

369

–

402

(

1999

).

https://doi.org/10.1016/S0377-0257(99)00069-5

Google Scholar

Crossref

42.

C.

Pozrikidis

,

Boundary Integral and Singularity Methods for Linearized Viscous Flow

(

Cambridge University Press

,

1992

).

Google Scholar

Crossref

43.

B.

Rodenborn

,

C.-H.

Chen

,

H. L.

Swinney

,

B.

Liu

, and

H.

Zhang

, “

Propulsion of microorganisms by a helical flagellum

,”

Proc. Natl. Acad. Sci.

110

(

5

),

E338

–

E347

(

2013

).

https://doi.org/10.1073/pnas.1219831110

Google Scholar

Crossref

44.

J.

Rotne

and

S.

Prager

, “

Variational treatment of hydrodynamic interaction in polymers

,”

J. Chem. Phys.

50

(

11

),

4831

–

4837

(

1969

).

https://doi.org/10.1063/1.1670977

Google Scholar

Crossref

45.

M. J.

Shelley

and

T.

Ueda

, “

The Stokesian hydrodynamics of flexing, stretching filaments

,”

Phys. D

146

(

1

),

221

–

245

(

2000

).

https://doi.org/10.1016/S0167-2789(00)00131-7

Google Scholar

Crossref

46.

M.

Shin

and

D. L.

Koch

, “

Rotational and translational dispersion of fibres in isotropic turbulent flows

,”

J. Fluid Mech.

540

,

143

(

2005

).

https://doi.org/10.1017/S0022112005005690

Google Scholar

Crossref

47.

C.

Siewert

,

R.

Kunnen

,

M.

Meinke

, and

W.

Schröder

, “

Orientation statistics and settling velocity of ellipsoids in decaying turbulence

,”

Atmos. Res.

142

,

45

–

56

(

2014

).

https://doi.org/10.1016/j.atmosres.2013.08.011

Google Scholar

Crossref

48.

D. J.

Smith

, “

A boundary element regularized stokeslet method applied to cilia-and flagella-driven flow

,”

Proc. R. Soc. A

465

(

2112

),

3605

–

3626

(

2009

).

https://doi.org/10.1098/rspa.2009.0295

Google Scholar

Crossref

49.

S. E.

Spagnolie

and

E.

Lauga

, “

Comparative hydrodynamics of bacterial polymorphism

,”

Phys. Rev. Lett.

106

(

5

),

058103

(

2011

).

https://doi.org/10.1103/PhysRevLett.106.058103

Google Scholar

Crossref

PubMed

50.

B.

Tapley

,

E.

Celledoni

,

B.

Owren

, and

H. I.

Andersson

, “

A novel approach to rigid spheroid models in viscous flows using operator splitting methods

,”

Numer. Algorithms

81

,

1423

–

1441

(

2019

).

https://doi.org/10.1007/s11075-019-00666-1

Google Scholar

Crossref

51.

B. K.

Tapley

,

H. I.

Andersson

,

E.

Celledoni

, and

B.

Owren

, “

Computational methods for tracking inertial particles in discrete incompressible flows

,” arXiv:1907.11936 (

2019

).

Google Scholar

52.

I.

Thorp

and

J.

Lister

, “

Motion of a non-axisymmetric particle in viscous shear flow

,”

J. Fluid Mech.

872

,

532

–

559

(

2019

).

https://doi.org/10.1017/jfm.2019.367

Google Scholar

Crossref

53.

A.-K.

Tornberg

, “

Accurate evaluation of integrals in slender-body formulations for fibers in viscous flow

,” arXiv:2012.12585 (

2020

).

Google Scholar

54.

A.-K.

Tornberg

and

K.

Gustavsson

, “

A numerical method for simulations of rigid fiber suspensions

,”

J. Comput. Phys.

215

(

1

),

172

–

196

(

2006

).

https://doi.org/10.1016/j.jcp.2005.10.028

Google Scholar

Crossref

55.

A.-K.

Tornberg

and

M. J.

Shelley

, “

Simulating the dynamics and interactions of flexible fibers in Stokes flows

,”

J. Comput. Phys.

196

(

1

),

8

–

40

(

2004

).

https://doi.org/10.1016/j.jcp.2003.10.017

Google Scholar

Crossref

56.

L. N.

Trefethen

and

J.

Weideman

, “

The exponentially convergent trapezoidal rule

,”

SIAM Rev.

56

(

3

),

385

–

458

(

2014

).

https://doi.org/10.1137/130932132

Google Scholar

Crossref

57.

S.

Twomey

, “

On the numerical solution of fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature

,”

J. ACM (JACM)

10

(

1

),

97

–

101

(

1963

).

https://doi.org/10.1145/321150.321157

Google Scholar

Crossref

58.

B. J.

Walker

,

M. P.

Curtis

,

K.

Ishimoto

, and

E. A.

Gaffney

, “

A regularised slender-body theory of non-uniform filaments

,”

J. Fluid Mech.

899

,

A3

(

2020

).

https://doi.org/10.1017/jfm.2020.434

Google Scholar

Crossref

59.

B. J.

Walker

,

K.

Ishimoto

,

H.

Gadêlha

, and

E. A.

Gaffney

, “

Filament mechanics in a half-space via regularised stokeslet segments

,”

J. Fluid Mech.

879

,

808

–

833

(

2019

).

https://doi.org/10.1017/jfm.2019.723

Google Scholar

Crossref

60.

H.

Yamakawa

, “

Transport properties of polymer chains in dilute solution: Hydrodynamic interaction

,”

J. Chem. Phys.

53

(

1

),

436

–

443

(

1970

).

https://doi.org/10.1063/1.1673799

Google Scholar

Crossref

61.

B.

Zhao

,

E.

Lauga

, and

L.

Koens

, “

Method of regularized stokeslets: Flow analysis and improvement of convergence

,”

Phys. Rev. Fluids

4

(

8

),

084104

(

2019

).

https://doi.org/10.1103/PhysRevFluids.4.084104

Google Scholar

Crossref

2021

Author(s)

An integral model based on slender body theory, with applications to curved rigid fibers

I. INTRODUCTION

A. Fiber geometry

II. SLENDER BODY MODEL

III. DERIVATION AND JUSTIFICATION OF THE SLENDER BODY MODEL

A. Spectral comparison of slender body integral operators

1. Spectrum of the slender body PDE

2. Pre-regularization comparison

3. Regularized comparison

IV. NUMERICAL DISCRETIZATION OF THE SLENDER BODY MODEL

A. Solving the second-kind Fredholm integral equation

B. Application to the slender body model and convergence tests

1. Thin ring translating with unit velocity

2. Prolate spheroid with artificial fluid velocity field

C. Spectrum of the slender body operator in different geometries

V. DYNAMICS OF CURVED RIGID FIBERS

A. Dynamical equations

1. Overview and cost of algorithm

B. Numerical validation of model dynamics

1. Dissipation matrix of a prolate spheroid

2. Prolate spheroids rotating in shear flow

C. Dynamics of randomly curvy fibers

VI. CONCLUSIONS

ACKNOWLEDGMENTS

DATA AVAILABILITY

Appendix A: Modified Lighthill model

Appendix B: Convergence and error bounds of numerical method

1. Convergence of numerical method for closed loop geometry

Appendix C: Dissipation matrix of a prolate spheroid

Appendix D: Endpoint behavior of model

References

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

An integral model based on slender body theory, with applications to curved rigid fibers

I. INTRODUCTION

A. Fiber geometry

II. SLENDER BODY MODEL

III. DERIVATION AND JUSTIFICATION OF THE SLENDER BODY MODEL

A. Spectral comparison of slender body integral operators

1. Spectrum of the slender body PDE

2. Pre-regularization comparison

3. Regularized comparison

IV. NUMERICAL DISCRETIZATION OF THE SLENDER BODY MODEL

A. Solving the second-kind Fredholm integral equation

B. Application to the slender body model and convergence tests

1. Thin ring translating with unit velocity

2. Prolate spheroid with artificial fluid velocity field

C. Spectrum of the slender body operator in different geometries

V. DYNAMICS OF CURVED RIGID FIBERS

A. Dynamical equations

1. Overview and cost of algorithm

B. Numerical validation of model dynamics

1. Dissipation matrix of a prolate spheroid

2. Prolate spheroids rotating in shear flow

C. Dynamics of randomly curvy fibers

VI. CONCLUSIONS

ACKNOWLEDGMENTS

DATA AVAILABILITY

Appendix A: Modified Lighthill model

Appendix B: Convergence and error bounds of numerical method

1. Convergence of numerical method for closed loop geometry

Appendix C: Dissipation matrix of a prolate spheroid

Appendix D: Endpoint behavior of model

References

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only