2.1 Piecewise-linear filaments

We consider a planar slender inextensible filament represented by $N$ piecewise-linear segments of constant length and described by arclength parameter $s\in [0,L]$ , where $L$ is the filament length. The segment endpoints are taken to correspond to material points, with their positions being denoted $\boldsymbol{x}_{1},\ldots ,\boldsymbol{x}_{N+1}$ and where $\boldsymbol{x}_{i}$ and $\boldsymbol{x}_{i+1}$ correspond to the $i$ th segment. We will refer to $\boldsymbol{x}_{1}$ as the base of the filament and correspondingly $\boldsymbol{x}_{N+1}$ as the tip, and denote the constant arclength associated with the material point $\boldsymbol{x}_{i}$ as $s_{i}$ . As the filament is assumed to be planar, valid in the absence of external symmetry-breaking features, without loss of generality we assume that it lies in a plane spanned by orthogonal unit vectors $\boldsymbol{e}_{x}$ and $\boldsymbol{e}_{y}$ , with $\boldsymbol{e}_{z}$ completing the orthonormal right-handed triad $\boldsymbol{e}_{x}\boldsymbol{e}_{y}\boldsymbol{e}_{z}$ , and we may write $\boldsymbol{x}_{i}=x_{i}\boldsymbol{e}_{x}+y_{i}\boldsymbol{e}_{y}$ . Following Moreau et al. (Reference Moreau, Giraldi and Gadêlha2018) we note that, once the segment lengths $v_{1},\ldots ,v_{N}$ are prescribed, the filament may be completely described by the $N+2$ scalars $x_{1},y_{1},\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{N}$ , where $\unicode[STIX]{x1D703}_{i}$ is defined as the angle between the $i$ th segment and the unit vector $\boldsymbol{e}_{x}$ as shown in figure 1. Explicitly, from these $N+2$ variables we may recover the filament endpoints as

(2.1) $$\begin{eqnarray}\displaystyle \boldsymbol{x}_{j}=\boldsymbol{x}_{1}+\mathop{\sum }_{i=1}^{j-1}(\cos \unicode[STIX]{x1D703}_{i}\boldsymbol{e}_{x}+\sin \unicode[STIX]{x1D703}_{i}\boldsymbol{e}_{y})v_{i}, & & \displaystyle\end{eqnarray}$$

and, using dots to denote derivatives with respect to time, the velocity of each material point is given by

(2.2) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{x}}_{j}=\dot{\boldsymbol{x}}_{1}+\mathop{\sum }_{i=1}^{j-1}(-\text{sin}\unicode[STIX]{x1D703}_{i}\boldsymbol{e}_{x}+\cos \unicode[STIX]{x1D703}_{i}\boldsymbol{e}_{y})\dot{\unicode[STIX]{x1D703}}_{i}v_{i}, & & \displaystyle\end{eqnarray}$$

again expressible using only the $N+2$ variables due to the imposed geometrical constraints.

Figure 1.

The two-dimensional piecewise-linear filament. We consider an inextensible filament composed of $N$ straight segments each of length $v_{i}$ , connected at material points $\boldsymbol{x}_{i}$ with each segment making angle $\unicode[STIX]{x1D703}_{i}$ with the global $\boldsymbol{e}_{x}$ axis. The endpoints $\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{N+1}$ shall be referred to as the base and tip respectively, and we note that, given segment lengths $v_{i}$ , the entire filament may be described by the location of the base and the angles $\unicode[STIX]{x1D703}_{i}$ .

2.2 Force distributions via regularised Stokeslet segments

To describe the low Reynolds number fluid dynamics pertinent to the filament we utilise the recent work of Cortez (Reference Cortez2018), namely the method of regularised Stokeslet segments (RSS). Here we briefly recapitulate the formulation of this methodology as presented by Cortez, relating the force density applied on the fluid by the filament to the fluid velocity via non-local hydrodynamics.

As in the popular method of regularised Stokeslets (Cortez Reference Cortez2001), we begin by considering solutions of the three-dimensional smoothly forced Stokes equations, which may be stated for viscosity $\unicode[STIX]{x1D707}$ and force $\boldsymbol{f}$ applied on the fluid at the origin as

(2.3a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}=\unicode[STIX]{x1D735}p-\boldsymbol{f}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D716}},\quad p,\boldsymbol{u}\rightarrow 0\quad \text{as }|\boldsymbol{x}|\rightarrow \infty , & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is the fluid velocity, $p$ is the pressure and $\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D716}}$ is a smooth approximation to the Dirac delta distribution dependent on the small parameter $\unicode[STIX]{x1D716}$ . Following § 2 of Cortez (Reference Cortez2018) we take

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{\unicode[STIX]{x1D716}}(\boldsymbol{x})=\frac{15\unicode[STIX]{x1D716}^{4}}{8\unicode[STIX]{x03C0}(|\boldsymbol{x}|^{2}+\unicode[STIX]{x1D716}^{2})^{7/2}}, & & \displaystyle\end{eqnarray}$$

for which the solution to (2.3) is known and given by

(2.5a,b ) $$\begin{eqnarray}\displaystyle 8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}\boldsymbol{u}(\boldsymbol{x})=\left[\left(\frac{1}{R}+\frac{\unicode[STIX]{x1D716}^{2}}{R^{3}}\right)\unicode[STIX]{x1D644}+\frac{\boldsymbol{x}\boldsymbol{x}^{\top }}{R^{3}}\right]\boldsymbol{f},\quad R(\boldsymbol{x})=\sqrt{|\boldsymbol{x}|^{2}+\unicode[STIX]{x1D716}^{2}}. & & \displaystyle\end{eqnarray}$$

By linearity of the Stokes equations the velocity at a point $\tilde{\boldsymbol{x}}$ due to a distribution of regularised Stokeslets along the filament with force density $\boldsymbol{f}(s)$ is therefore given by

(2.6a,b ) $$\begin{eqnarray}\displaystyle 8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}\boldsymbol{u}(\tilde{\boldsymbol{x}})=\int _{0}^{L}\left[\left(\frac{1}{R(\boldsymbol{r})}+\frac{\unicode[STIX]{x1D716}^{2}}{R(\boldsymbol{r})^{3}}\right)\unicode[STIX]{x1D644}+\frac{\boldsymbol{r}\boldsymbol{r}^{\top }}{R(\boldsymbol{r})^{3}}\right]\boldsymbol{f}(s)\,\text{d}s,\quad \boldsymbol{r}(s)=\tilde{\boldsymbol{x}}-\boldsymbol{x}(s).\quad & & \displaystyle\end{eqnarray}$$

Whilst the original method of regularised Stokeslets would entail approximating the force distribution by a finite number of smoothed point forces, the method of regularised Stokeslet segments instead considers a linear distribution of forces along the straight segments of the discretised filament. With the discretisation of the force distribution as a continuous piecewise-linear function along each of the segments, the fluid velocity is instead given by a sum of integrals over each of the segments, where each individual integral may be analytically evaluated. Parameterising the $j$ th segment by $\unicode[STIX]{x1D6FC}\in [0,1]$ so that $\boldsymbol{x}=\boldsymbol{x}_{j}-\unicode[STIX]{x1D6FC}\boldsymbol{v}$ for $\boldsymbol{v}=\boldsymbol{x}_{j}-\boldsymbol{x}_{j+1}$ , and additionally writing the force density as $\boldsymbol{f}=\boldsymbol{f}_{j}+\unicode[STIX]{x1D6FC}(\boldsymbol{f}_{j+1}-\boldsymbol{f}_{j})$ for force densities $\boldsymbol{f}_{j},\boldsymbol{f}_{j+1}$ at $\boldsymbol{x}_{j},\boldsymbol{x}_{j+1}$ respectively, the integral along the $j$ th segment is given by

(2.7) $$\begin{eqnarray}\displaystyle v_{j}\int _{0}^{1}\left[\left(\frac{1}{R}+\frac{\unicode[STIX]{x1D716}^{2}}{R^{3}}\right)\unicode[STIX]{x1D644}+\frac{(\tilde{\boldsymbol{x}}-\boldsymbol{x}_{j}+\unicode[STIX]{x1D6FC}\boldsymbol{v})(\tilde{\boldsymbol{x}}-\boldsymbol{x}_{j}+\unicode[STIX]{x1D6FC}\boldsymbol{v})^{\top }}{R^{3}}\right][\boldsymbol{f}_{j}+\unicode[STIX]{x1D6FC}(\boldsymbol{f}_{j+1}-\boldsymbol{f}_{j})]\,\text{d}\unicode[STIX]{x1D6FC}, & & \displaystyle\end{eqnarray}$$

where we have identified $|\boldsymbol{v}|=v_{j}$ and $R=\sqrt{|\tilde{\boldsymbol{x}}-\boldsymbol{ x}_{j}+\unicode[STIX]{x1D6FC}\boldsymbol{ v}|^{2}+\unicode[STIX]{x1D716}^{2}}$ . As noted by Cortez (Reference Cortez2018, § 2.2) this may be written as a sum of integrals of the form

(2.8) $$\begin{eqnarray}\displaystyle T_{m,q}=\int _{0}^{1}\unicode[STIX]{x1D6FC}^{m}R^{q}\,\text{d}\unicode[STIX]{x1D6FC} & & \displaystyle\end{eqnarray}$$

for $(m,q)\in \{(0,-1),(1,-1),(0,-3),(1,-3),(2,-3),(3,-3)\}$ , each of which is implicitly dependent on $j$ and may be evaluated by explicit computation of $T_{0,-1}$ and $T_{0,-3}$ and subsequent application of the recurrence relation

(2.9) $$\begin{eqnarray}\displaystyle v_{j}^{2}T_{m+1,q}=-\frac{m}{(q+2)}T_{m-1,q+2}-(\tilde{\boldsymbol{x}}-\boldsymbol{x}_{j})\boldsymbol{\cdot }\boldsymbol{v}T_{m,q}+\frac{1}{q+2}\unicode[STIX]{x1D6FC}^{m}R^{q+2}|_{0}^{1},\quad q\neq -2, & & \displaystyle\end{eqnarray}$$

which follows from the definition of $T_{m,q}$ and consideration of $(\text{d}/\text{d}\unicode[STIX]{x1D6FC})(\unicode[STIX]{x1D6FC}^{m}R^{q})$ .

The simple extension of this methodology to a half-space bounded by a no-slip planar wall via the image system of Ainley et al. (Reference Ainley, Durkin, Embid, Boindala and Cortez2008) (with the typographical error corrected Smith (Reference Smith2009)) is presented in detail by Cortez (Reference Cortez2018, § 3.3). Exemplified in figure 2, we will consider two configurations of half-space, one in which the infinite planar boundary is parallel to the plane containing the filament, given by $z=0$ , and another where the boundary is situated perpendicular to this plane and given by $y=0$ . In both cases the hydrodynamics may be described by the same regularised singularity representation, each requiring computation of $T_{m,q}$ for the additional values $(m,q)\in \{(0,-5),(1,-5),(2,-5),(3,-5),(4,-5),(5,-5)\}$ in order to include the additional necessary regularised singularities. Omitted from the previous work of Cortez, we compute

(2.10) $$\begin{eqnarray}\displaystyle T_{0,-5}=\left.\frac{(B+\unicode[STIX]{x1D6FC}C^{2})}{3R^{3}}\frac{(-B^{2}+C^{2}[3A^{2}+4\unicode[STIX]{x1D6FC}B+2\unicode[STIX]{x1D6FC}^{2}C^{2}+3\unicode[STIX]{x1D716}^{2}])}{(B^{2}-C^{2}[A^{2}+\unicode[STIX]{x1D716}^{2}])^{2}}\right|_{0}^{1}, & & \displaystyle\end{eqnarray}$$

where $A=|\tilde{\boldsymbol{x}}-\boldsymbol{x}_{j}|$ , $B=(\tilde{\boldsymbol{x}}-\boldsymbol{x}_{j})\boldsymbol{\cdot }\boldsymbol{v}$ and $C=|\boldsymbol{v}|=v_{j}$ . The remaining values of $T_{m,q}$ may be computed via the recurrence relation equation (2.9). Simple computation of the coefficients of the $T_{m,q}$ results in a matricial form of the velocity contribution at $\tilde{\boldsymbol{x}}$ from the $j$ th segment as a linear operator acting on $\boldsymbol{f}_{j}$ and $\boldsymbol{f}_{j+1}$ , the details of which are cumbersome and given in the supplementary material, available at https://doi.org/10.1017/jfm.2019.723. The overall velocity at the evaluation point $\tilde{\boldsymbol{x}}$ is then given by the sum of these contributions, resulting in a matricial equation of the form

(2.11) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(\tilde{\boldsymbol{x}})=\unicode[STIX]{x1D648}(\tilde{\boldsymbol{x}})\boldsymbol{F}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{F}=(f_{1,x},f_{1,y},\ldots ,f_{N+1,x},f_{N+1,y})^{\top }$ herein denotes the composite vector of force densities applied on the fluid at the segment endpoints $\boldsymbol{x}_{j}$ in the directions $\boldsymbol{e}_{x}$ and $\boldsymbol{e}_{y}$ , where we write $\boldsymbol{f}_{j}=f_{j,x}\boldsymbol{e}_{x}+f_{j,y}\boldsymbol{e}_{y}$ . Taking the evaluation point $\tilde{\boldsymbol{x}}$ as each $\boldsymbol{x}_{i}$ in turn yields a square system of linear equations in $2(N+1)$ variables, relating the velocities $\boldsymbol{u}(\boldsymbol{x}_{i})$ to the force distribution on the filament via non-local hydrodynamics. Application of the no-slip condition at the segment endpoints then gives a relation between the filament velocities and the forces applied on the fluid by the filament, which we write in brief as

(2.12) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{X}}=\unicode[STIX]{x1D63C}\boldsymbol{F}, & & \displaystyle\end{eqnarray}$$

for the square matrix $\unicode[STIX]{x1D63C}$ composed row-wise of blocks $\unicode[STIX]{x1D648}(\boldsymbol{x}_{i})$ for $i=1,\ldots ,N+1$ , and where $\dot{\boldsymbol{X}}=({\dot{x}}_{1},{\dot{y}}_{1},\ldots ,{\dot{x}}_{N+1},{\dot{y}}_{N+1})^{\top }$ . Given kinematic data, this linear system may be readily solved to give the force densities applied on the fluid by the filament.

Figure 2.

Schematic configurations of planar filaments near boundaries. (a) A filament moving parallel to the infinite planar boundary, at constant separation $h$ . The projection of the filament onto the boundary is shown as a thin black line. (b) A filament contained in a plane perpendicular to the boundary, with the plane of motion shown dashed.

2.3 Extension to efficient solution of elastohydrodynamic equations

We now proceed to combine the work of Cortez (Reference Cortez2018) and Moreau et al. (Reference Moreau, Giraldi and Gadêlha2018) to give an efficient scheme for solving non-local planar elastohydrodynamics, similar in concept to the piecewise-constant force density approach of Hall-McNair et al. (Reference Hall-McNair, Gallagher, Montenegro-Johnson, Gadêlha and Smith2019) but with continuous piecewise-linear force discretisation and the additional inclusion of an infinite planar boundary. As formulated by Moreau et al. (Reference Moreau, Giraldi and Gadêlha2018), we integrate the pointwise conditions of force and moment balance on the filament, given by

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{n}_{s}-\boldsymbol{f}=\mathbf{0}, & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{m}_{s}+\boldsymbol{x}_{s}\times \boldsymbol{n}=\mathbf{0} & \displaystyle\end{eqnarray}$$

for contact force and couple denoted $\boldsymbol{n},\boldsymbol{m}$ respectively and where a subscript of $s$ denotes differentiation with respect to arclength, noting that we have assumed slenderness of the filament. This yields the integrated equations

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle -\mathop{\sum }_{j=1}^{N}\int _{s_{j}}^{s_{j+1}}\boldsymbol{f}(s)\,\text{d}s=\boldsymbol{n}(0), & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle -\mathop{\sum }_{j=i}^{N}\int _{s_{j}}^{s_{j+1}}(\boldsymbol{x}(s)-\boldsymbol{x}_{i})\times \boldsymbol{f}(s)\,\text{d}s=\boldsymbol{m}(s_{i}),\quad i=1,\ldots ,N. & \displaystyle\end{eqnarray}$$

Here we have assumed conditions of zero contact force and couple at the tip of the filament, i.e. $\boldsymbol{n}(L)=\boldsymbol{m}(L)=\mathbf{0}$ , retaining generality at the base for the time being. With the assumption of piecewise-linear distributions of force density $\boldsymbol{f}$ over each segment, again with $\boldsymbol{f}=\boldsymbol{f}_{j}+\unicode[STIX]{x1D6FC}(\boldsymbol{f}_{j+1}-\boldsymbol{f}_{j})$ on the $j$ th segment for $\unicode[STIX]{x1D6FC}\in [0,1]$ , we may write these integrals as linear operators on $\boldsymbol{F}$ , yielding the system

(2.17) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D63D}\boldsymbol{F}=\boldsymbol{R}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{R}=(\boldsymbol{n}(s_{1})\boldsymbol{\cdot }\boldsymbol{e}_{x},\boldsymbol{n}(s_{1})\boldsymbol{\cdot }\boldsymbol{e}_{y},m(s_{1}),\ldots ,m(s_{N}))^{\top }$ and we are writing $\boldsymbol{m}(s_{i})=m(s_{i})\boldsymbol{e}_{z}$ by planarity of the filament. The matrix $\unicode[STIX]{x1D63D}$ is of dimension $(N+2)\times 2(N+1)$ , which under the assumption of equispaced segment endpoints has rows $\boldsymbol{B}_{k}$ given by

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}_{1}=\frac{\unicode[STIX]{x0394}s}{2}[1,0,2,0,2,\ldots ,2,0,1,0], & \displaystyle\end{eqnarray}$$

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}_{2}=\frac{\unicode[STIX]{x0394}s}{2}[0,1,0,2,0,2,\ldots ,2,0,1], & \displaystyle\end{eqnarray}$$

with the additional rows $\boldsymbol{B}_{i+2}$ each encoding the integrated expression

(2.20) $$\begin{eqnarray}\displaystyle \boldsymbol{e}_{z}\boldsymbol{\cdot }\mathop{\sum }_{j=i}^{N}\left\{\left[\frac{\unicode[STIX]{x0394}s}{2}(\boldsymbol{x}_{j}-\boldsymbol{x}_{i})+\frac{\unicode[STIX]{x0394}s^{2}}{6}\boldsymbol{t}_{j}\right]\times \boldsymbol{f}_{j}+\left[\frac{\unicode[STIX]{x0394}s}{2}(\boldsymbol{x}_{j}-\boldsymbol{x}_{i})+\frac{\unicode[STIX]{x0394}s^{2}}{3}\boldsymbol{t}_{j}\right]\times \boldsymbol{f}_{j+1}\right\}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x0394}s=L/N$ is the length of each segment, $\boldsymbol{t}_{j}=\cos \unicode[STIX]{x1D703}_{j}\boldsymbol{e}_{x}+\sin \unicode[STIX]{x1D703}_{j}\boldsymbol{e}_{y}$ is the local unit tangent and $i=1,\ldots ,N$ . Thus we may write the coupled elastohydrodynamic problem as

(2.21) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D63D}\unicode[STIX]{x1D63C}^{-1}\dot{\boldsymbol{X}}=\boldsymbol{R}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D63D}$ represents integration over segments and $\unicode[STIX]{x1D63C}$ encodes the relationship between the forces on the surrounding fluid and the velocities of material points on the filament, here non-local and assumed to be invertible.

As noted in § 2.1, the material velocities $\dot{\boldsymbol{x}}_{i}$ may be expressed in terms of the time derivatives of the base point $\boldsymbol{x}_{1}$ and the segment angles $\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{N}$ . Defining $\unicode[STIX]{x1D64C}$ to be the $2(N+1)\times (N+2)$ matrix such that

(2.22) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64C}\dot{\unicode[STIX]{x1D73D}}=\dot{\boldsymbol{X}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D73D}=(x_{1},y_{1},\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{N})^{\top }$ , we have $\unicode[STIX]{x1D64C}$ given explicitly by

(2.23)

where the subscript $P$ denotes that the $i$ th row of $\unicode[STIX]{x1D64C}$ is to be permuted to

(2.24) $$\begin{eqnarray}\displaystyle P(i)=\left\{\begin{array}{@{}ll@{}}2(i-1)+1, & i=1,\ldots ,N+1,\\ 2(i-N-1), & i=N+2,\ldots ,2N+2.\end{array}\right. & & \displaystyle\end{eqnarray}$$

This permutation of $\unicode[STIX]{x1D64C}$ allows us to define the blocks $\unicode[STIX]{x1D64C}_{1}$ and $\unicode[STIX]{x1D64C}_{2}$ as being strictly lower triangular matrices of dimension $(N+1)\times N$ with entries

(2.25) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{1}^{i,j}=\unicode[STIX]{x0394}s\left\{\begin{array}{@{}ll@{}}-\sin \unicode[STIX]{x1D703}_{j}, & j<i,\\ 0, & j\geqslant i,\end{array}\right. & \displaystyle\end{eqnarray}$$

(2.26) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{2}^{i,j}=\unicode[STIX]{x0394}s\left\{\begin{array}{@{}ll@{}}+\cos \unicode[STIX]{x1D703}_{j}, & j<i,\\ 0, & j\geqslant i.\end{array}\right. & \displaystyle\end{eqnarray}$$

Given the matrix $\unicode[STIX]{x1D64C}$ we may rewrite the full coupled elastohydrodynamic problem simply as

(2.27) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D63D}\unicode[STIX]{x1D63C}^{-1}\unicode[STIX]{x1D64C}\dot{\unicode[STIX]{x1D73D}}=\boldsymbol{R}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D63D}\unicode[STIX]{x1D63C}^{-1}\unicode[STIX]{x1D64C}$ being a square matrix of dimension $(N+2)\times (N+2)$ . Finally, and here for example assuming force and moment-free conditions at the filament base, the standard linear constitutive relation $m(s_{i})=EI\,\unicode[STIX]{x1D703}_{s}(s_{i})\approx EI\,(\unicode[STIX]{x1D703}_{i}-\unicode[STIX]{x1D703}_{i-1})/\unicode[STIX]{x0394}s$ for bending stiffness $EI$ , valid for $i=2,\ldots ,N$ , gives $\boldsymbol{R}$ explicitly in terms of $\unicode[STIX]{x1D73D}$ , yielding a linear system of low dimension that may be readily solved for $\dot{\unicode[STIX]{x1D73D}}$ , with the temporal dynamics then computable via existing ordinary differential equation (ODE) methods. As noted by Moreau et al., this formulation is readily extensible to the imposition of a variety of boundary conditions, in particular that of a cantilevered filament base, achieved by replacing the overall zero-force and zero-moment equations with ${\dot{x}}_{1}={\dot{y}}_{1}=\dot{\unicode[STIX]{x1D703}}_{1}=0$ .

We non-dimensionalise this system by scaling spatial coordinates with filament length $L$ , forces with $EI/L^{2}$ and time with some characteristic time scale $T$ , yielding the non-dimensional system

(2.28) $$\begin{eqnarray}\displaystyle -E_{h}\hat{\unicode[STIX]{x1D63D}}\hat{\unicode[STIX]{x1D63C}}^{-1}\hat{\unicode[STIX]{x1D64C}}\dot{\hat{\unicode[STIX]{x1D73D}}}=\hat{\boldsymbol{R}},\quad E_{h}=\frac{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}L^{4}}{EI\,T} & & \displaystyle\end{eqnarray}$$

for elastohydrodynamic number $E_{h}$ , where the notation $\hat{\cdot }$ denotes non-dimensional quantities. The dimensional and non-dimensional quantities are related by

(2.29a-d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63D}=L^{2}\hat{\unicode[STIX]{x1D63D}},\quad \unicode[STIX]{x1D63C}=\frac{1}{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}\hat{\unicode[STIX]{x1D63C}},\quad \unicode[STIX]{x1D64C}\dot{\unicode[STIX]{x1D73D}}=\frac{L}{T}\hat{\unicode[STIX]{x1D64C}}\dot{\hat{\unicode[STIX]{x1D73D}}},\quad \boldsymbol{R}=\frac{EI}{L}\hat{\boldsymbol{R}}, & & \displaystyle\end{eqnarray}$$

where we have multiplied the two force balance equations by $\unicode[STIX]{x0394}s=L/N$ and absorbed the dimensional scalings of $x_{1},y_{1}$ into $\unicode[STIX]{x1D64C}$ for convenience, with $\hat{\unicode[STIX]{x1D73D}}=(\hat{x}_{1},{\hat{y}}_{1},\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{N})^{\top }$ . In order to pose the non-dimensional equations in terms of the commonly used sperm number $Sp$ (Delmotte, Climent & Plouraboué Reference Delmotte, Climent and Plouraboué2015; Ishimoto & Gaffney Reference Ishimoto and Gaffney2018; Moreau et al. Reference Moreau, Giraldi and Gadêlha2018), we proceed following Ishimoto & Gaffney (Reference Ishimoto and Gaffney2018) to define

(2.30a,b ) $$\begin{eqnarray}\displaystyle Sp^{4}=\frac{\unicode[STIX]{x1D709}L^{4}}{EI\,T},\quad \unicode[STIX]{x1D709}=\frac{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\log (2L/\unicode[STIX]{x1D716})} & & \displaystyle\end{eqnarray}$$

for resistive force coefficient $\unicode[STIX]{x1D709}$ and $L/\unicode[STIX]{x1D716}=10^{3}$ , giving the approximate relation $E_{h}\approx 15.2\,Sp^{4}$ .

2.4 Implementation, verification and parameter choice

Both the calculation of force densities from kinematic data and the solution of full elastohydrodynamics were implemented in MATLAB^®. We utilise the inbuilt stiff ODE solver ode15s (Shampine & Reichelt Reference Shampine and Reichelt1997), making use of variable step sizes in order to satisfy configurable error tolerances that are typically set here at $10^{-6}$ , with the presented results insensitive to this choice. Prior to the recent work of Hall-McNair et al. (Reference Hall-McNair, Gallagher, Montenegro-Johnson, Gadêlha and Smith2019) the solution of non-local elastohydrodynamics has required significant computational work and minimal timestep for the solution of this stiff problem (Olson et al. Reference Olson, Lim and Cortez2013; Ishimoto & Gaffney Reference Ishimoto and Gaffney2018), and we replicate in our implementation the low computational cost associated with the integrated elasticity equations of Moreau et al. (Reference Moreau, Giraldi and Gadêlha2018). In particular, typical simulations of a cantilevered filament in background flow, explored in detail in §§ 3.4 and 3.5, have a typical runtime of 10 s, where $N=40$ and we simulate over 10 periods of oscillation of the background flow on modest hardware (Intel^® Core™ i7-6920HQ CPU).

Verification was first performed on the reduction to an unbounded domain, with the full elastohydrodynamics being verified by comparison with the resistive force theory results of Moreau et al. (Reference Moreau, Giraldi and Gadêlha2018) for filament relaxation. Further drag calculations were compared against the work of Cortez (Reference Cortez2018) for the case of uniform motion, with the degree to which the no-slip condition is numerically satisfied along the filament length having previously been investigated in detail (Cortez Reference Cortez2018, § 3.1). The regularised image system in a half-space was then verified by explicit evaluation of the fluid velocity on the no-slip boundary, with the numerical result being zero to machine precision, and further by noting that far-field results were seen to converge to those of the free-space system. Further qualitative checks were performed, an example being the successful reproduction of the intuitive result that relaxation time scales increase for filaments with reduced separation from the no-slip stationary boundary, along with numerical comparison against the previous works of Pozrikidis regarding the elastohydrodynamics of filaments in steady shear flow (Pozrikidis Reference Pozrikidis2010, Reference Pozrikidis2011), a brief example of which being presented in appendix B. Results of additional verification against the boundary element method of Ramia et al. (Reference Ramia, Tullock and Phan-Thien1993) are shown in figure 3(a). Sufficient accuracy is typically achieved with $N=20$ segments, similar to the discretisation used by Cortez (Reference Cortez2018), although we typically take $N=40$ to enable the capturing of high-curvature filaments.

As noted by Cortez (Reference Cortez2018), the use of regularised Stokeslet segments allows the regularisation parameter to represent the radius of the filament being modelled, verified here by comparison with the boundary element method of Ramia et al. (Reference Ramia, Tullock and Phan-Thien1993) in figure 3, subject to the phenomenological condition that $\unicode[STIX]{x1D716}$ be less than the length of the segments, which appears necessary for convergence. We will take $\unicode[STIX]{x1D716}/L=10^{-3}$ unless otherwise stated, with typical slender filament aspect ratios yielding $\unicode[STIX]{x1D716}/L$ in the range $(10^{-2},10^{-3})$ (Yonekura, Maki-Yonekura & Namba Reference Yonekura, Maki-Yonekura and Namba2003; Ishimoto & Gaffney Reference Ishimoto and Gaffney2016). Whilst the results that follow naturally depend on the size the of the regularisation parameter, it being used as a proxy for the filament radius, this dependence appears to be predominantly quantitative, with the same qualitative conclusions holding across the range of physically relevant values of $\unicode[STIX]{x1D716}/L$ .

2.5 Endpoint effects

Inherent to the method of regularised Stokeslet segments as presented here is the presence of apparent large variations in computed force densities at the tips of filaments. This was initially observed by Cortez (Reference Cortez2018, § 3.1, figure 4), who remarked that these endpoint effects may be reduced when $\unicode[STIX]{x1D716}/L$ is small. Indeed, for $\unicode[STIX]{x1D716}/L$ in our range of physical interest these endpoint effects contribute minimally to overall drag calculations, in particular having little effect on the computed total drag exerted on a filament and the computed force densities away from the filament tip, whose presentation we focus on herein. Exploration of the effects of non-uniform segment lengths yields the observation that these oscillatory endpoint effects remain limited to approximately the 3 segments proximal to the ends of the filament, thus the integral contribution of these oscillations may be reduced by the clustering of segments near the filament tips. For our typical parameters of $N=40$ and $\unicode[STIX]{x1D716}/L=10^{-3}$ we observe a difference in integrated drag along of filaments of less than 2 % between linear, quadratic and Chebyshev segment endpoints, hence we proceed with calculations utilising segments of uniform length.

2.6 Boundary-corrected resistive force theory

The leading-order approximation of slender body hydrodynamics that is resistive force theory (RFT) was first introduced by Gray and Hancock (Hancock Reference Hancock1953; Gray & Hancock Reference Gray and Hancock1955), relating the local drag on a body to its local tangential and normal velocities $u_{t},u_{n}$ via

(2.31a,b ) $$\begin{eqnarray}\displaystyle f_{t}=-C_{t}u_{t},\quad f_{n}=-C_{n}u_{n}. & & \displaystyle\end{eqnarray}$$

Here $f_{t},f_{n}$ are the local tangential and normal components of force applied on the fluid, with the constant resistive coefficients $C_{t},C_{n}$ typically being functions of filament aspect ratio. As RFT relates forces only to local velocities, calculations using this simple local drag theory are not able to account for the presence of a boundary. Brenner (Reference Brenner1962) and Katz et al. (Reference Katz, Blake and Paveri-Fontana1975) posed corrections to resistive force theory for straight filaments in asymptotic regimes far from or near to an infinite no-slip boundary, with the latter having been verified against high-accuracy boundary element methods by Ramia et al. (Reference Ramia, Tullock and Phan-Thien1993). In the case of filament motion confined to a plane parallel to the boundary, as exemplified in figure 2(a), by this wall-corrected resistive force theory (W-RFT) a straight filament parallel to the boundary has resistive coefficients given by

(2.32a,b ) $$\begin{eqnarray}\displaystyle C_{t}=\frac{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\log \left({\displaystyle \frac{2}{\unicode[STIX]{x1D716}}}\right)-0.807-{\displaystyle \frac{3L}{8h}}},\quad C_{n}=\frac{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\log \left({\displaystyle \frac{2}{\unicode[STIX]{x1D716}}}\right)+0.193-{\displaystyle \frac{3L}{4h}}}\quad \text{if }L\ll h,\quad \qquad & & \displaystyle\end{eqnarray}$$

(2.33a,b ) $$\begin{eqnarray}\displaystyle C_{t}=\frac{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\log \left({\displaystyle \frac{2h}{\unicode[STIX]{x1D716}}}\right)},\quad C_{n}=\frac{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\log \left({\displaystyle \frac{2h}{\unicode[STIX]{x1D716}}}\right)}\quad \text{if }L\gg h, & & \displaystyle\end{eqnarray}$$

where $h$ is the boundary separation of the filament and $\unicode[STIX]{x1D716}$ has been taken as the filament radius. Similarly, and as in figure 2(b), when aligned parallel to the boundary and moving in a plane perpendicular to the surface the coefficients are given by

(2.34a,b ) $$\begin{eqnarray}\displaystyle C_{t}=\frac{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\log \left({\displaystyle \frac{2}{\unicode[STIX]{x1D716}}}\right)-0.807-{\displaystyle \frac{3L}{8h}}},\quad C_{n}=\frac{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\log \left({\displaystyle \frac{2}{\unicode[STIX]{x1D716}}}\right)+0.193-{\displaystyle \frac{3L}{2h}}}\quad \text{if }L\ll h,\quad \qquad & & \displaystyle\end{eqnarray}$$

(2.35a,b ) $$\begin{eqnarray}\displaystyle C_{t}=\frac{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\log \left({\displaystyle \frac{2h}{\unicode[STIX]{x1D716}}}\right)},\quad C_{n}=\frac{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}}{\log \left({\displaystyle \frac{2h}{\unicode[STIX]{x1D716}}}\right)-1}\quad \text{if }L\gg h, & & \displaystyle\end{eqnarray}$$

where all coefficients are as summarised by Brennen & Winet (Reference Brennen and Winet1977) and (2.35) and (2.33) additionally require $\unicode[STIX]{x1D716}\ll h$ . In free space we will adopt the resistive force coefficients defined by the large- $h$ limit of (2.32) and (2.34), and refer to this simpler theory as free-space resistive force theory (F-RFT). As given by Brenner and Katz et al., the leading-order boundary corrections to $C_{t}$ and $C_{n}$ are $O(L/h)^{3}$ when $L\ll h$ , and $O(h/L)$ when $L\gg h$ , although the overall error in the resistive force approximation remains logarithmic in the filament aspect ratio.