1. Introduction
The fluids produced and used by living organisms are inherently complex, containing a suspended microstructure such as deformable cells or polymers molecules that are immersed in a viscous fluid. Prime examples of such fluids are blood and mucus. While blood can be viewed as a suspension of deformable cells, mucus is instead a mesh of fibre bundles formed from the protein mucin and connected through entanglement or hydrophobicity (Lai et al. Reference Lai, Wang, Cone, Wirtz and Hanes2009). In our own bodies, mucus plays several important roles as a lubricating fluid (Nordgård & Draget Reference Nordgård and Draget2011; Boegh & Nielsen Reference Boegh and Nielsen2015), a protective layer (Allen & Flemström Reference Allen and Flemström2005; Mantelli & Argüeso Reference Mantelli and Argüeso2008) and as a barrier to pathogens and infection (Johansson et al. Reference Johansson, Phillipson, Petersson, Velcich, Holm and Hansson2008; Ensign, Cone & Hanes Reference Ensign, Cone and Hanes2014). Mucus helps facilitate digestion, protecting the soft tissue of the stomach from the acids used during digestion, and from bacteria such as Helicobacter pylori that cause stomach ulcers (Celli et al. Reference Celli2009). In relation to reproduction, mucus is present in the cervix (Ceric, Silva & Vigil Reference Ceric, Silva and Vigil2005; Rutllant, López-Béjar & López-Gatius Reference Rutllant, López-Béjar and López-Gatius2005), again serving to prevent infections, but also providing the fluid environment through which sperm swim. Mucin fibre density is known to decrease when fertilisation is possible, and sperm–mucus interactions have been put forward as a possible mechanism for selection of viable sperm (Bianchi et al. Reference Bianchi, De Agostini, Fournier, Guidetti, Tarozzi, Bizzaro and Manicardi2004; Suarez & Pacey Reference Suarez and Pacey2006; Holt & Fazeli Reference Holt and Fazeli2015).
As one might expect, the presence of the microstructure, which produces changes in the rheological properties of the surrounding fluid, alters the motility of cells moving through the fluid (Sznitman & Arratia Reference Sznitman and Arratia2015). Sperm cells swim by undulating their flagellum, and in purely viscous fluids, increased viscosity is seen to be associated with lower wave propagation velocity and smaller wave amplitude in sperm, but not necessarily with any change in average swimming speed (Brokaw Reference Brokaw1966; Smith et al. Reference Smith, Gaffney, Gadêlha, Kapur and Kirkman-Brown2009). But human cervical mucus is viscoelastic and shear thinning (Wolf et al. Reference Wolf, Blasco, Khan and Litt1977), and increased viscosity in mucus is linked to enhanced swimming speed (López-Gatius et al. Reference López-Gatius, Rutllan, López-Béjar and Labèrnia1994; Rutllant et al. Reference Rutllant, López-Béjar and López-Gatius2005). Conversely, experiments with Caenorhabditis elegans, which also moves through undulation, show that increased fluid elasticity leads to slower swimming speeds (Shen & Arratia Reference Shen and Arratia2011). As increasing the microstructure in the fluid changes the fluid viscoelasticity, in general, the relationship between suspension concentration and microorganism motility is not simple. As another example, experiments show that Escherichia coli, which propels itself through rotation of its helical flagella, travels faster in low-concentration polymeric fluids (up to approximately 20 %) than in a purely viscous fluid, but at high concentrations of polymers, travels slower (Schneider & Doetsch Reference Schneider and Doetsch1974; Martinez et al. Reference Martinez, Schwarz-Linek, Reufer, Wilson, Morozov and Poon2014).
To gain a more fundamental understanding of microorganism motility in complex biological fluids, researchers have turned to mathematical models based on continuum-level descriptions of non-Newtonian fluids. Initial theoretical work with Oldroyd-B and FENE-P viscoelastic models, on infinitely long undulating sheets and filaments, typically predicted that additional viscoelasticity would impede swimming (Lauga Reference Lauga2007; Fu, Wolgemuth & Powers Reference Fu, Wolgemuth and Powers2009). Later, computational studies in the nonlinear regime for finite-length undulatory swimmers painted a more complex picture, with both hindrance and enhancement (up to 30 %) found for different mechanical properties of the swimmer in a two-dimensional (2-D) Oldroyd-B fluid (Teran, Fauci & Shelley Reference Teran, Fauci and Shelley2010; Thomases & Guy Reference Thomases and Guy2014). The same broad picture was found in 2-D asymptotic and numerical work modelling the background fluid with the Brinkman equations for porous media (Leiderman & Olson Reference Leiderman and Olson2016), with both hindrance and enhancement (up to 45 %) found for finite-length swimmers with low and high bending stiffness, respectively. The same model in three dimensions also showed hindrance and enhancement (up to 15 %) (Cortez et al. Reference Cortez, Cummins, Leiderman and Varela2010; Olson & Leiderman Reference Olson and Leiderman2015), depending on bending stiffness.
While these models provided a number of key insights into how fluid rheology affects microorganism motility, the continuum fluid model assumes a large separation in length scale between the swimming microorganism and the microstructure, describing just one limit in a much larger parameter space (‘type II’ in Spagnolie & Underhill Reference Spagnolie and Underhill2023). In many notable cases, the characteristic length scales of the swimming body and the immersed microstructure are comparable (the border between ‘types II and IV’). For sperm moving in cervical mucus, for example, the gap sizes in the mucin fibre bundle mesh (${\sim }1\ \mathrm {\mu }{\rm m}$) are of the same order as the head of the cell (${\sim }5\ \mathrm {\mu }{\rm m}$) (Sheehan, Oates & Carlstedt Reference Sheehan, Oates and Carlstedt1986; Lai et al. Reference Lai, Wang, Cone, Wirtz and Hanes2009). In this scenario, it is therefore important to consider the fluid microstructure as discrete objects comparable in size to the swimmer. One intermediate method is to represent these discrete objects as point forces which act on the swimmer, where the spacing between the forces (e.g. in a lattice or a mesh) is of the same order as the swimmer size (Gniewek & Kolinski Reference Gniewek and Kolinski2010). Simulations of finite-length swimmers through such a mesh network have found mostly hindrance (Schuech, Cortez & Fauci Reference Schuech, Cortez and Fauci2022) but also speed enhancement of up to 30 % (Wróbel et al. Reference Wróbel, Lynch, Barrett, Fauci and Cortez2016). Simulations of infinitely long swimmers through bead–spring dumbbell suspensions (Zhang, Li & Ardekani Reference Zhang, Li and Ardekani2018) have shown both hindrance and enhancement (up to 10 %) for sufficiently short-wavelength swimmers, depending on the relative viscosity of the fluid around the cell body. Experiments and simulations (Park et al. Reference Park, Hwang, Nam, Martinez, Austin and Ryu2008; Majmudar et al. Reference Majmudar, Keaveny, Zhang and Shelley2012) of undulatory locomotion in lattices of rigid obstacles (‘type I’) have shown that speeds of finite-length swimmers can increase by nearly a factor of ten when the lattice spacing allows for the swimmer to constantly push itself along as it undulates. After introducing randomness in the obstacle positions and environment compliance through tethering springs (Kamal & Keaveny Reference Kamal and Keaveny2018), enhanced locomotion continues to be observed, although by a more modest but still substantial factor of two and a half. In addition, modelling the microstructure discretely also captures swimmer velocity fluctuations (Jabbarzadeh, Hyon & Fu Reference Jabbarzadeh, Hyon and Fu2014; Kamal & Keaveny Reference Kamal and Keaveny2018) due to collisions, which in turn produces effective swimmer diffusion over long times. Non-trivial changes in motility are also observed for swimmers propelled by helical flagella in heterogeneous environments. A recent theory to explain increased swimming speeds at low suspension concentrations is that the swimming motion of the microorganism leads to phase separation near the surface of the swimmer. Multiparticle collision dynamics simulations of finite-length helically propelled swimmers (Zöttl & Yeomans Reference Zöttl and Yeomans2019) suggest this could lead to speed increases between 10 % and 100 %. This mechanism yields similar increases in swimming speed for infinitely long undulatory swimmers (Man & Lauga Reference Man and Lauga2015). For recent reviews of theoretical and computational approaches to modelling swimming in complex media, see Li, Lauga & Ardekani (Reference Li, Lauga and Ardekani2021), Spagnolie & Underhill (Reference Spagnolie and Underhill2023); for a broad overview of sperm modelling, see Gaffney, Ishimoto & Walker (Reference Gaffney, Ishimoto and Walker2021).
In this paper, we perform simulations of an undulatory swimmer moving through 3-D suspensions of elastic filaments. We resolve hydrodynamic interactions using the force coupling method and elastic deformation of the swimmer and filaments using the computational framework described in Schoeller et al. (Reference Schoeller, Townsend, Westwood and Keaveny2021). A summary of the model is presented in § 2. In performing the simulations, we consider cases where the filaments are freely suspended in the fluid (§ 3), or tethered to one another via Hookean springs to form a network (§ 4). In both types of environments, we examine how the swimming speed is affected by filament concentration, filament stiffness and tether stiffness. In sufficiently compliant environments, we find that the swimming speed increases to a value of 10 % greater than that in viscous fluid for volume fractions up to 20 %. For stiff environments, however, the speed is non-monotonic with volume fraction, with peak enhancement of 5 % and reductions in speed as the volume fraction approaches 20 %. By systematically removing physical effects from the simulation, we show that speed enhancement is due to the hydrodynamic interactions, whereas speed reductions in stiff environments are due to reduced waveform amplitude as the swimmer becomes more confined.
2. Mathematical model for the swimmer and environment
2.1. General model for all filaments
In our study, we simulate the movement of an active, undulating swimmer, modelled as an active filament, travelling through a fully 3-D suspension or network of passive filaments. To perform these simulations, we employ the filament model developed in Schoeller et al. (Reference Schoeller, Townsend, Westwood and Keaveny2021). The model is a 3-D development of the 2-D model of swimming filaments first introduced in Majmudar et al. (Reference Majmudar, Keaveny, Zhang and Shelley2012), adapted for sperm in Schoeller & Keaveny (Reference Schoeller and Keaveny2018) and used in studies of 2-D obstacle arrays in Kamal & Keaveny (Reference Kamal and Keaveny2018). For full discussion, we refer the reader to Schoeller et al. (Reference Schoeller, Townsend, Westwood and Keaveny2021, § 2); we will summarise the model here in the context of a single filament.
In the model, the filament, sketched in figure 1, has length $L$ and a circular cross-section of radius $a$. The filament is parametrised by its arclength, $0 \leq s \leq L$, such that the position of any point on the filament centreline at time $t$ is $\boldsymbol {Y}(s,t)$. Additionally, we introduce a right-handed orthonormal frame at any $s$ and $t$, $\{\hat {\boldsymbol {t}}(s,t),\hat {\boldsymbol {\mu }}(s,t),\hat {\boldsymbol {\nu }}(s,t)\}$, that we will employ to keep track of the deformation of the filament.
To derive the equations of motion, we begin by considering the force and moment balances along the filament (Landau & Lifshitz Reference Landau and Lifshitz1986, § 19; Powers Reference Powers2010)
where $\boldsymbol {\varLambda }(s,t)$ and $\boldsymbol {M}(s,t)$ are the internal force and moment on the filament cross-section. The internal moments, $\boldsymbol {M}(s,t)$, are linearly related to the filament bending and twist through the constitutive law
where $K_B$ is the bending modulus and $K_T$ is the twist modulus. We also have preferred curvatures $\kappa _{\mu }(s,t)$ and $\kappa _{\nu }(s,t)$ and preferred twist, $\gamma _0(s,t)$, that can incorporate non-trivial equilibrium filament shapes (Lim et al. Reference Lim, Ferent, Wang and Peskin2008; Lim Reference Lim2010; Olson, Lim & Cortez Reference Olson, Lim and Cortez2013) and time-dependent filament motion. For the passive filaments comprising the suspension or network, the preferred curvatures and twist are zero, but for the swimmer, there is a time-dependent $\kappa _\nu$, as defined in § 2.2. The internal force, $\boldsymbol {\varLambda }(s,t)$, enforces the kinematic constraint
unlike the internal moment, the internal force is a priori unknown and will be solved for in the final part of the numerical method.
Also appearing in (2.2) are $\boldsymbol {f}$ and $\boldsymbol {\tau }$ which are the external force and torque per unit length on the filament, respectively. These consist of short-ranged, repulsive forces when filaments collide, but also the force and torque per unit length exerted on the filament by the surrounding fluid.
Next, we discretise the filament into $N$ segments of length $\Delta L$. Each segment has position $\boldsymbol {Y}_n$ and orientation vector $\hat {\boldsymbol {t}}_n$. After applying central differencing to (2.2), we have
where $\boldsymbol {M}_{n+1/2}$ and $\boldsymbol {\varLambda }_{n+1/2}$ are the internal moment and force between segments $n$ and $n+1$, and $\boldsymbol {f}_n$ and $\boldsymbol {\tau }_n$ are the external force and torque, respectively, per unit length on segment $n$. For the discrete system, $\boldsymbol {\varLambda }_{n+1/2}$ can be viewed as the Lagrange multiplier that enforces
the discrete version of (2.4).
Multiplying (2.5), (2.6) through by $\Delta L$, we arrive at the force and torque balances on each segment $n$
where $\boldsymbol {F}^C_n =\boldsymbol {\varLambda }_{n+1/2} - \boldsymbol {\varLambda }_{n-1/2}$, $\boldsymbol {T}^E_n = \boldsymbol {M}_{n+1/2} - \boldsymbol {M}_{n-1/2}$ and $\boldsymbol {T}^C_n = ({\Delta L}/{2}) \hat {\boldsymbol {t}}_n\times (\boldsymbol {\varLambda }_{n+1/2} + \boldsymbol {\varLambda }_{n-1/2}).$ The external forces and torques acting on each segment are $\boldsymbol {F}_n = \Delta L \boldsymbol {f}_n$ and $\boldsymbol {T}_n = \Delta L \boldsymbol {\tau }_n$ and are composed of the repulsive, or barrier, forces between segments and the hydrodynamic force and torque on each segment, $\boldsymbol {F}_n = \boldsymbol {F}^B_n - \boldsymbol {F}^H_n$ and $\boldsymbol {T}_n = - \boldsymbol {T}^H_n$. The barrier force, $\boldsymbol {F}^B_n$, prevents segments from overlapping and acts between pairs of segments of different filaments, or between non-neighbouring segments of the same filament. Specifically, expressing the centre-to-centre displacement between two segments as $\boldsymbol {r}_{nm} = \boldsymbol {Y}_n - \boldsymbol {Y}_m$, then $\boldsymbol {F}^B_n = \sum _m \boldsymbol {F}^B_{nm}$ where (Dance, Climent & Maxey Reference Dance, Climent and Maxey2004)
and $\boldsymbol {F}^B_{nm} = \boldsymbol {0}$ otherwise, where $r_{nm} = |\boldsymbol {r}_{nm}|$ and $\chi = 1.1$. The force strength $F^S$ is set to $51K_B/L^2$, where $K_B$ and $L$ are the bending modulus and length of the swimmer described in § 2.2. The steric forces serve to prevent the filament and swimmer segments from overlapping. The choice of $F^S = 51K_B/L^2$ is chosen to ensure that overlapping does not occur while also not introducing numerical stiffness that would lead to large increases in the number of Broyden iterations needed during a timestep. Additionally, this force captures steric interactions between the swimmer and filaments that are expected to occur in real swimmer–biopolymer suspensions where lubrication interactions break down due to polymer roughness and surface aberrations.
At this point we now have full expressions for the hydrodynamic forces and torques, $\boldsymbol {F}^H_n$ and $\boldsymbol {T}^H_n$, acting on all filament segments, in terms of the unknown Lagrange multipliers $\boldsymbol {\varLambda }_n$. As such, we now consider the hydrodynamic model which will allow us to close the system and obtain the motion of the filaments using an iterative method. Due to the negligible influence of fluid inertia, the hydrodynamic interactions are governed by the Stokes equations for a fluid with viscosity $\eta$. Due to the linearity of the Stokes equations, the velocities, $\boldsymbol {V}_n$, and angular velocities, $\boldsymbol {\varOmega }_n$, of the segments are linearly related to the hydrodynamic forces and torques on the segments, leading to the mobility problem
where $\boldsymbol {V}^\mathsf {T} = (\boldsymbol {V}_1^\mathsf {T},\ldots,\boldsymbol {V}_N^\mathsf {T})$ is the vector of all segment velocity components, and similarly for $\boldsymbol {\varOmega }$, $\boldsymbol {F}^H$ and $\boldsymbol {T}^H$.
The computation of the mobility matrix, $\mathcal {M}$, (or, more efficiently, the computation of its action on the force and torque vector) can be performed with a number of different hydrodynamic models and computational methodologies. In this paper, we use the matrix-free force coupling method (FCM) (Maxey & Patel Reference Maxey and Patel2001; Lomholt & Maxey Reference Lomholt and Maxey2003; Liu et al. Reference Liu, Keaveny, Maxey and Karniadakis2009). In FCM, the forces and torques exerted on the fluid by the segments are transferred through a truncated, regularised force multipole expansion in the Stokes equations. The delta functions in the multipole expansion are replaced by smooth Gaussians, and, using the ratios established for spherical particles, the envelope size of the Gaussians is based on the filament radius, $a$. After finding the fluid velocity, the same Gaussian functions are used to average over the fluid volume to obtain the velocity, $\boldsymbol {V}_n$, and angular velocity, $\boldsymbol {\varOmega }_n$, of each segment, $n$.
Having obtained the velocities and angular velocities of the segments, we now integrate in time. For the segment positions, we integrate
while for segment orientations, we have
where ${{\boldsymbol{\mathsf{q}}}}_n$ is the unit quaternion that provides the rotation matrix
To advance in time, we discretise (2.12), (2.13) using a geometric second-order backward differentiation scheme (Ascher & Petzold Reference Ascher and Petzold1998; Iserles et al. Reference Iserles, Munthe-Kaas, Nørsett and Zanna2000; Faltinsen, Marthinsen & Munthe-Kaas Reference Faltinsen, Marthinsen and Munthe-Kaas2001) that ensures that the quaternions remain unit length. The resulting equations, along with (2.7), yield a system of nonlinear equations that we solve using Broyden's method (Broyden Reference Broyden1965), a quasi-Newton iterative method, to obtain the updated segment positions, $\boldsymbol {Y}_n$, the quaternions associated with the segment orientations, ${{\boldsymbol{\mathsf{q}}}}_n$, and the Lagrange multipliers, $\boldsymbol {\varLambda }_n$.
For full details of the filament model, especially in regard to the quaternion representation and geometric integration scheme, we again refer the reader to Schoeller et al. (Reference Schoeller, Townsend, Westwood and Keaveny2021, § 2).
2.2. Swimmer
In all our simulations, the swimmer has length $L=33a$, is discretised into $N=15$ segments and its bending modulus and twisting modulus are equal, $K_T = K_B$. Swimmer motion is generated by planar undulations, driven by a preferred curvature
where $k = 2{\rm \pi} /L$. To connect with previous work in Kamal & Keaveny (Reference Kamal and Keaveny2018), we choose $K_0 = 8.25/L$ and set $\omega$ to yield a dimensionless sperm number $Sp = (4{\rm \pi} \omega \eta /K_B)^{1/4}L = 5.87$. The linear decay in amplitude for $s > L/2$ is introduced to limit otherwise very high curvatures toward the end of the swimmer (Majmudar et al. Reference Majmudar, Keaveny, Zhang and Shelley2012; Kamal & Keaveny Reference Kamal and Keaveny2018; Schoeller & Keaveny Reference Schoeller and Keaveny2018). The other preferred curvature, $\kappa _\mu$, and the preferred twist, $\gamma _0$, are set to zero. Under these conditions, in a quiescent fluid, a single swimmer swims its length $L$ in a time of $11.5T$, where $T = 2{\rm \pi} /\omega$ is the undulation period. The resultant periodic waveform has a wavelength equal to the length of the swimmer, broadly matching observations of sperm cells and C. elegans (Rikmenspoel Reference Rikmenspoel1965; Sznitman et al. Reference Sznitman, Purohit, Krajacic, Lamitina and Arratia2010); it is depicted in figure 11(a).
2.3. Domain size
All simulations are performed in a cubic periodic domain of size $87a\times 87a\times 87a = 2.64L \times 2.64L \times 2.64L$, where $L$ is the swimmer length, with a grid size of $288\times 288\times 288$ used to solve the Stokes equations as part of the FCM mobility problem. The domain size is chosen to ensure periodic images have a negligible effect on the swimmer's motion. Illustrations of the domain will be presented in the next section.
3. Unconnected filament suspensions
In our simulations, we will consider the motion of the swimmer through two types of filament environments: a suspension of unconnected filaments, and a spring-connected filament network. In this section, we look at the former, first summarising the approach that we use to generate this environment and then presenting the results of our simulations.
3.1. Filament environment
Producing the filament suspension follows a common and straightforward algorithm. We first place an undeformed swimmer (length $L=33a$) in the periodic domain. We then randomly seed (in position and orientation) $M$ undeformed filaments of length $L=33a$. This will inevitably lead to overlapping between filaments, which must be addressed before the simulation is run. For low filament volume fraction, $\phi < 5\,\%$, where $\phi = 4 {\rm \pi}a^3 L M / (3V\Delta L)$, overlap is resolved as part of the seeding process itself, where filaments are checked for overlap as they are seeded and the seed is rejected if overlap occurs. This is continued until all $M$ filaments are successfully introduced to the domain. For high volume fractions, $\phi > 5\,\%$, filaments are placed at random positions and orientations in the domain, regardless of any overlap. To remove overlap, the filaments are allowed to evolve under very short-ranged repulsive forces with segment mobility based on Stokes drag. The filaments are allowed to evolve until all overlap is removed and this final configuration is used as the initial condition for our simulations. An example of a suspension simulation initial condition for $\phi = 11.2\,\%$ is shown in figure 2(a).
Suspension simulations are conducted using two values of the non-dimensional relative bending modulus, $K_B' := K_B/K_B^{swim} = 10^{-3}$ and $1$, where $K_B^{swim}$ is the bending modulus of the swimmer. These values correspond to relatively flexible and stiff networks, respectively. We use non-dimensional parameters in our study, allowing us to transcend a link to any specific organism and environment, but we briefly present dimensional measurements from the literature for context. For C. elegans, estimates of $K_B^{swim}$ range widely, from $O(10^{-16})$ to $O(10^{-13})\ \text {N m}^2$ (Fang-Yen et al. Reference Fang-Yen, Wyart, Xie, Kawai, Kodger, Chen, Wen and Samuel2010; Sznitman et al. Reference Sznitman, Purohit, Krajacic, Lamitina and Arratia2010; Backholm, Ryu & Dalnoki-Veress Reference Backholm, Ryu and Dalnoki-Veress2013). For mammalian sperm cells, the estimated range is at least $O(10^{-21})$ to $O(10^{-19})\ \text {N m}^2$ (Rikmenspoel Reference Rikmenspoel1965; Lindemann, Macauley & Lesich Reference Lindemann, Macauley and Lesich2005; Gadêlha Reference Gadêlha2012). Estimates for mucin fibres are scarce in the literature, partly due to significant variation in their diameter, which ranges from 3 to 10 nm (Shogren, Gerken & Jentoft Reference Shogren, Gerken and Jentoft1989). Given that $K_B$ scales with the fourth power of the diameter for long, thin tubes, this results in a variation spanning four orders of magnitude. Assuming a typical biopolymer Young's modulus of $1\ \mathrm {GPa}$, the estimated $K_B$ ranges from $O(10^{-25})$ to $O(10^{-22})\ \mathrm {N\ m}^2$.
Simulations are also performed for six different filament volume fractions between $\phi = 1.9\,\%$ and $18.9\,\%$. Videos of simulations with $\phi =7.6\,\%$ for both $K'_B=10^{-3}$ and $1$ are included in the supplementary material available at https://doi.org/10.1017/jfm.2024.603.
3.2. Swimming speeds in filament suspensions
We start by recording and measuring the swimming speed of the undulating filament as it moves through the filament suspensions. For each value of $K'_B$ and $\phi$, we perform 10 independent simulations, each with different randomly generated initial filament configurations, run for 20 periods of undulation. As the active filament swims, its stroke adds periodic contributions to the swimmer's speed. To remove these contributions, the swimming speed is defined as the period average of the centre of mass velocity in the direction of the period-averaged swimmer orientation (Kamal & Keaveny Reference Kamal and Keaveny2018).
At any given time $t$, the instantaneous velocity of the centre of mass of the swimmer is given by
where $\boldsymbol {V}_n$ is the velocity of segment $n$. The instantaneous orientation vector of the swimmer is given by $\hat {\boldsymbol {p}} = \boldsymbol {p}/|\boldsymbol {p}|$ where
and where $\hat {\boldsymbol {t}}_n$ is the tangent vector of segment $n$. We can then calculate the period-averaged velocity of the swimmer, $\boldsymbol {V}_{period}(t)$, for time $t\geq T$, as
while the period-averaged orientation is $\hat {\boldsymbol {P}} = \boldsymbol {P}/|\boldsymbol {P}|$, where
The swimming speed, $V_{swim}(t)$, at any given time is then defined as
In figure 3(a), we present the normalised swimming speed $V_{swim}(t) / \langle V_{swim}^0 \rangle$, where $\langle V_{swim}^0 \rangle$ is the swimming speed in the absence of filaments, for each independent simulation, for increasing concentration and both values of the bending modulus. Along with this, we show in (b) of the same figure, $\langle V_{swim} \rangle /\langle V_{swim}^0 \rangle$, the average value of the swimming speed across both time and independent simulations. We see that in suspensions of relatively flexible filaments ($K'_B = 10^{-3}$), the speed of the swimmer is promoted up to 15 % due to the presence of the filaments, with more concentrated suspensions yielding higher speeds. For suspensions with stiffer filaments ($K'_B=1$), a small speed promotion of up to 5 % is observed at low concentrations, but at higher concentrations (above $\phi \sim 12\,\%$), the swimmer is on average hindered, displaying speed decreases of up to $20\,\%$.
The speed enhancement at low concentrations is consistent with the speed increases of up to 20 % seen in simulations of 2-D undulatory swimmers in Stokes–Oldroyd-B viscoelastic fluids (Teran et al. Reference Teran, Fauci and Shelley2010; Thomases & Guy Reference Thomases and Guy2014). In common with the planar undulatory motion in Kamal & Keaveny (Reference Kamal and Keaveny2018), the average swimmer slowdown for rigid, concentrated networks is not due to a uniform speed decrease across all independent runs; instead, it is due to large variations in the speed, including periods of time where the swimmer is trapped within the network in some simulations. But in contrast, our speed enhancement is significantly less than in simulations of interrupted planar swimming, where fixed (Majmudar et al. Reference Majmudar, Keaveny, Zhang and Shelley2012) or tethered (Kamal & Keaveny Reference Kamal and Keaveny2018) discrete obstacles lead to three- to tenfold speed increases. We also see less enhancement than in helical swimming with fixed obstacles (at 40 %: Leshansky Reference Leshansky2009), where swimming speed can increase monotonically with concentration (Klingner et al. Reference Klingner, Mahdy, Hanafi, Adel, Misra and Khalil2020). We discuss mechanisms for these distinctions in § 6.
The fluctuations in the swimmer speed seen in figure 3(a) increase with concentration, but this effect is much larger for the stiffer networks. At the highest concentration, the swimmer is sometimes sped up by a factor of 2, but also sometimes forced backwards in a number of simulations. This is quantified by the standard deviation in panel (c) of the same figure, where the mean standard deviation, $\langle \sigma \rangle$, increases almost linearly with concentration and reaches $0.4$ of the unhindered swimmer speed in the most concentrated suspensions of stiff filaments.
4. Spring-connected filament networks
4.1. Filament environment
A number of methods for generating linked filament networks exist in the literature. One approach, seen in 2-D (Head, Levine & MacKintosh Reference Head, Levine and MacKintosh2003a,Reference Head, Levine and MacKintoshb; DiDonna & Levine Reference DiDonna and Levine2006; Heussinger & Frey Reference Heussinger and Frey2006) and 3-D (Åström et al. Reference Åström, Kumar, Vattulainen and Karttunen2008) models of actin networks, is to seed straight filaments isotropically and add cross-links where filaments intersect, potentially removing loose ends (Wilhelm & Frey Reference Wilhelm and Frey2003; Huisman et al. Reference Huisman, van Dillen, Onck and Van der Giessen2007). Another option is to seed straight filaments, but let them grow into other shapes according to an algorithm, before forming cross-links (Buxton, Clarke & Hussey Reference Buxton, Clarke and Hussey2009).
Here, we use a probability-based algorithm to construct a phenomenological model of a generic network, based on Buxton & Clarke (Reference Buxton and Clarke2007). We place the swimmer, of length $L$, in the empty periodic domain and seed $N_{nodes}$ nodes randomly inside the domain at a distance at least $R_{min} = 2.2LN_{nodes}^{-1/3}$ from each other. The exponent of $-1/3$ corresponds to hard sphere packing inside a finite domain. For each pair of nodes, $i$ and $j$, with positions $\boldsymbol {r}_i$ and $\boldsymbol {r}_j$, respectively, we place a filament connecting these nodes with a probability $P_{conn}$ if the distance between nodes is $|\boldsymbol {r}_i - \boldsymbol {r}_j| = |\boldsymbol {r}_{ij}| < R_{max} = 1.2L$. The filament is placed along the vector $\boldsymbol {r}_{ij}$, with its centre at the midpoint of $\boldsymbol {r}_{ij}$ and formed of $N=\lfloor (r_{ij}-2r_{spacing})/\Delta L \rfloor$ segments, where $r_{spacing} = 2a$ is taken to ensure filaments do not overlap at the nodes. Linear spring forces $-k_s(x-\ell )\hat {\boldsymbol {x}}$, where $\boldsymbol {x}=x\hat {\boldsymbol {x}}$ is the end-to-end vector and $\ell$ is the natural spring length, are applied between all pairs of filament end segments that meet at a node. In doing this, we generate a network spanning the periodic domain with no external force holding it in place; in practice, the drag on the component filaments limits the network from moving as a rigid body. The seeding algorithm is summarised in figure 2(b).
To ensure that filaments do not overlap, this seeding procedure is carried out in a separate simulation, using friction-based mobility and an active collision barrier, until an equilibrium has been found. We once again define the volume concentration of the network by treating each segment as a sphere with its hydrodynamic radius $a$. That is, for $M$ filaments of length $L_i$ in a periodic domain with volume $V$, we have $\phi = 4 {\rm \pi}a^3 \sum _{i=1}^M{L_i} / (3V\Delta L)$. Under this definition, the concentration scales as the square of the number of nodes, $\phi \sim N_{nodes}^2$. Over a large number of independent random initial seeds, the concentration is normally distributed with a standard deviation that grows linearly with $N_{nodes}$.
Simulations are performed for different concentrations of filaments between $\phi = 3.4\,\%$ and $18.7\,\%$. We choose $N_{nodes}$ and $P_{conn}$ to produce networks of a given concentration with a common mean filament length of approximately $0.75L$, and common mean number of connections at each node (to within a reasonable tolerance). Our parameter choices are summarised in table 1.
The pair distribution function (PDF) of node separations under these configurations is displayed in figure 4(a). We see a bump near $R_{min}$ corresponding to close packing, but otherwise the PDF follows the shape for a uniform distribution, slightly enhanced. A render of one such network at 11.3 % volume concentration is shown in figure 2(c), with the nodes represented by black circles.
Simulations are performed for two different values of the network filament bending modulus, $K'_B = 10^{-3}$ and $1$ (as for the filament suspensions), and for three values of the spring modulus, $k'_s = 0,0.36,3.6$, where $k'_s := k_sL^3/K_B^{swim}$. The swimmer bending and twist modulus remains constant, and the twist modulus, $K_T$, for every filament, is set equal to its bending modulus.
Simulations of flagellar swimming through a viscoelastic network formed of virtually cross-linked nodes by Wróbel et al. (Reference Wróbel, Lynch, Barrett, Fauci and Cortez2016) show that the connectivity number matters: higher-connectivity networks are stiffer. Here, we keep this number constant so that we can explore the effect of the other parameters in the system: the network filament bending modulus, the filament concentration and the spring modulus of the networks. The distribution of connection numbers in each case is shown in figure 4(b), where the mean connectivity number in each of our simulations is between 16 and 19. For comparison, cubic lattices with the diagonals excluded and included have connectivity numbers of 6 and 26, respectively, and were used in Wróbel et al. (Reference Wróbel, Lynch, Barrett, Fauci and Cortez2016). Buxton & Clarke (Reference Buxton and Clarke2007) use connectivity numbers between 5 and 15.
4.2. Swimming speeds in filament networks
We now present results from simulations examining swimmer motion in the connected filament networks. We again explore the effects of filament bending stiffness and volume fraction, but now also include the spring constant of the connections. For each triplet of parameters $(K'_B,\phi,k'_s)$, five simulations are performed, each with a different initial random configuration, and each running for 20 undulation periods.
The mean speed and standard deviation over the simulation run, $\langle V_{swim} \rangle$ and $\langle \sigma \rangle$, are presented in figure 5. We perform simulations for three different values of the spring constant. In the case where $k'_s=0$, the simulations differ from the filament suspensions in § 3 only in the seeding algorithm. For comparison, $\langle V_{swim} \rangle$ and $\langle \sigma \rangle$ from the suspension simulations are reproduced in this figure. We see that, for $k'_s=0$ and $0.36$, the trend in the data follows that found for filament suspension simulations: swimmer speed is promoted with concentration in flexible networks ($K'_B=10^{-3}$) up to $\sim 15\,\%$, but rigid networks ($K'_B=1$) produce a small speed promotion followed by a hindrance as concentration increases.
The more striking differences with the suspensions simulations are observed when $k'_s=3.6$, where the filament tethers are stiffer than the swimmer. Here, the swimmer speed in the flexible network decreases monotonically with concentration, down to an approximate reduction of $12\,\%$. The swimmer is hindered even further in the stiff filament network, with speeds decreasing by approximately $35\,\%$ at the highest concentrations. Interestingly, at intermediate concentrations, there is little difference with the weaker-spring scenarios. The effects of the stronger spring constant also yield considerably higher fluctuations in swimming speed with a standard deviation of up to $0.2\langle V_{swim}^0 \rangle$, although notably not as high as in the unconnected suspensions from § 3 where the standard deviation can be as high as $0.4 \langle V_{swim}^0 \rangle$.
5. Investigation into mechanisms for motility changes
In the previous sections, we observed that the swimming speed is mildly promoted by the presence of filaments at low concentrations: the mere existence of the filaments (either flexible or stiff) in suspension produces this effect. The addition of weak connections between the filaments can provide an additional small speed increase across all concentrations and all rigidities. However, for cases with stiffer filaments or stiffer connecting springs, high concentrations of filaments tend to slow the swimmer down.
In this section, we aim to explore the mechanisms at play by which swimmer speed is promoted or hindered. Specifically, we seek to gain insight into these trends by examining the forces acting on the swimmer and changes in swimming gait due to the presence of the filaments. The repulsive barrier forces act on the swimmer as it passes through the suspension or network of filaments. These forces act to oppose contact with suspended filaments and can therefore act on the swimmer in all directions. We can consider these forces to act in two main ways: (i) by considering the total force acting on the swimmer, and (ii) by considering how the shape of the swimmer changes as a result of unequal forces along the swimmer body.
5.1. The role of steric forces directly
The barrier force can result in a non-trivial net force on the swimmer that can push or pull it through the environment. We compute the mean period-averaged barrier force in the swimming direction, $\langle F_{swim} \rangle$. This can be measured analogously to $\langle V_{swim}\rangle$ in (3.1)–(3.5), namely
noting that since the constraint forces on the swimmer sum to zero, the total barrier force is equal to the total hydrodynamic force, $\sum _n \boldsymbol {F}_n^B = \sum _n \boldsymbol {F}_n^H$.
The total force, $\langle F_{swim} \rangle$, in the suspensions is presented in figure 6(a). For all but the most concentrated case with the stiff filaments, the steric forces act against the direction of swimming. This is somewhat surprising as everywhere where this force is opposite the swimming direction, the velocity plot in figure 3 showed us that the swimmer is travelling faster than it does in the absence of filaments. The high-concentration, high-stiffness case appears an outlier, being pushed on average forwards, indicative of a different regime which will appear again in the spring networks.
Although the direction of the mean force appears at odds with the increase in swimming speed, an examination of the time series before averaging shows a direct correlation between force and speed fluctuations, as shown in figure 6(b). But, still, the region between the horizontal lines in this panel represents times when the forces oppose swimming, and yet the swimmer is still travelling faster than it would in the absence of filaments, suggestive of a uniform promotion in speed across the simulation.
The filament networks follow similar trends but paint a more nuanced picture, with the total forces presented in figure 7. As in the suspensions, the steric forces from interactions with flexible network filaments push against the direction of swimming, with larger concentrations leading to stronger opposing forces. Once again, this happens despite the swimmer moving faster in these networks than they would at 0 % concentration (figure 5), although this is not universal: at the highest spring constant the filament speed does drop below its filament-free value. Interestingly, increasing the spring constant from $k'_s=0$ to $0.36$ makes little difference to the force yet mildly promotes swimming speed. At low filament concentrations, the swimmer behaves similarly in stiff and flexible filament networks, however, we see that above a critical concentration the forces from stiff networks will push in the direction of swimming yet the swimmer slows down. The force is much weaker for the $(k'_s=0, K'_B=1)$ case than in the matching suspension case, which we attribute to the differences in filament seeding. Increasing the spring constant recaptures and enhances the behaviour seen in the suspension.
These results suggest that there are two force regimes for the swimmer: one where the swimmer is gently resisted, almost linearly with concentration; and one where the swimmer is pushed forwards by the network. What is surprising is that these two force regimes do not clearly translate to observable swimming speed regimes. In particular, the barrier forces cannot be responsible for the speed increase seen in suspensions at lower-to-mid concentrations. We conclude that something else is acting to mask the total applied force.
5.2. The role of the swimming gait
Along with modifying the swimmer's rigid body motion, steric forces can also act to change the swimming gait. Swimmer undulations produce collisions with the surrounding filaments, inducing lateral forces which can reduce the amplitude of the undulation and alter the beat plane, making the gait less effective. Sample gaits from the suspension simulations, in both the beating and out-of-beating plane, are illustrated at increasing concentration in figure 8(a). These stills show that in flexible suspensions, increased concentration leads to slightly reduced beating-plane amplitudes and minor out-of-plane movements. But in rigid suspensions, especially at the highest concentration (bottom row), sample gaits deviate significantly from beating behaviour in a filament-free environment, with out-of-plane amplitude sometimes surpassing in-plane amplitude.
We can quantify these observations by measuring the mean amplitude in the beating plane, presented in figure 8(b) for suspensions and in figure 9(a) for our networks. We find that increased filament stiffness, $K'_B$, and for the networks, increased spring constants, $k'_s$, lead to smaller amplitudes at higher concentrations.
We can assess further the efficacy of these reduced-amplitude gaits by examining how well they propel the swimming in the absence of filaments. To do this, we remove the rigid body motion from the swimmer in the filament environment and recompute it in a viscous fluid absent of any filaments using the remaining undulations. We describe this computation in detail in Appendix A. The resulting speed when the imposed-gait swimmer swims is shown in figures 8(c) and 9(b). Correlating strongly with the stroke amplitude in both the suspensions and the networks, we see that the imposed gait is less effective at swimming, and increasingly so for the higher $K'_B$ and (in the network) higher $k'_s$. Thus, any gait changes result in reduced swimming speeds rather than the increases that we observe in the full simulations.
5.3. The role of hydrodynamics
So far we have seen that steric interactions, at mid-to-low concentrations, hinder swimming speed directly through rigid body forces. Additionally, when we remove all suspended filaments and impose the swimming gait seen in the full simulations, the swimmer again moves with a reduced speed compared with its natural gait in an empty fluid.
To understand the speed increases noted in figure 3, having eliminated steric forces and swimming gait, we turn our attention to hydrodynamic interactions. To investigate this, we run our suspension simulations again, including all filaments but now without hydrodynamic interactions between the suspended bodies. We accomplish this by solving the fluid mobility problem in (2.11) using a resistive force theory (RFT) (Johnson & Brokaw Reference Johnson and Brokaw1979; Lauga & Powers Reference Lauga and Powers2009) with coefficients tuned to recover the correct filament dynamics as in Schoeller & Keaveny (Reference Schoeller and Keaveny2018). The details are provided in Appendix B.
The resulting force and speed of the swimmer in these simulations without hydrodynamics are shown in figure 10 and, similar to the gait analysis simulations in figure 8, show no modest increase in speed at low-to-mid concentrations. Unlike simulations with full hydrodynamics (reproduced from figure 3 in light grey in the background of figure 10), we see that the swimming speed in this concentration regime is strongly affected by the steric forces, as evidenced by the strong correlation between these forces (a) and the resulting speed (b). At these concentrations, also observe the reduction in swimmer amplitude present in the full simulations, figure 8.
Although the general shape of the total force graph is similar between RFT and the full simulations, see figures 6 and 10, it is noticeable that the backwards force response of $K'_B=10^{-3}$ is stronger, and of $K'_B=1$ is weaker, without the hydrodynamics. This can be attributed to the difference in hydrodynamic model since, when the hydrodynamics is included, both the swimmer and the suspended filaments can deform each other without contact; whereas in the drag-based model, the filaments only affect each other on contact through the repulsive forces. In short, the presence of hydrodynamics in figure 3 significantly changes the swimmer behaviour from a steric force driven regime seen in the absence of hydrodynamics, in figure 10.
5.4. Summary of investigation
In this examination of swimmer motion, we have analysed the swimming gait and the forces experienced by the swimmer. While we find that the fluctuations in the swimming speed are closely correlated with collisions between the swimmer and filaments, enhancement in the mean speed is not linked to the force experienced by the swimmer or changes in swimming gait. In fact, both effects would produce a monotonic decrease in the speed with filament concentration. This suggests that the enhancement is instead linked to hydrodynamic interactions between the swimmer and filaments. We confirm this by performing RFT simulations where hydrodynamic interactions between the swimmer and filaments are removed and for which we do not observe any enhanced speed.
6. Discussion and conclusions
In this paper, we performed simulations examining the changes in swimming speed experienced by an undulatory swimmer in filament suspensions and immersed filament networks. The simulations reveal a non-monotonic dependence of the speed on the filament concentration, with initial increases of the swimming speed to peak values of 5 %–15 % greater than the speed in a viscous fluid. These increases in speed are in line with predictions of swimming speeds using continuum models to describe the surrounding medium. Filament and connection flexibility are required to enhance the speed. Filaments and connections that are too stiff will instead hinder locomotion. Along with the speed, we examine speed fluctuations, which cannot be captured with continuum models. These fluctuations are found to be significant and increase monotonically with filament concentration.
To understand these changes in swimmer motion we examined closely the changes in swimming gait and force experienced by the swimmer. Although contact between the swimmer and its environment correlate with speed fluctuations, we cannot attribute the mean speed increase to these interactions, or indeed to changes in the swimming gait, since we see that independently these effects lead to a monotonic decrease in speed with filament concentration. Instead, we attribute the speed enhancement to hydrodynamic interactions between the swimmer and the filaments in the environment by showing that when we remove them, in RFT simulations, we no longer observe enhanced speed.
Previous work on planar undulatory motion between tethered obstacles showed a tenfold increase in the swimming speed in structured environments (Majmudar et al. Reference Majmudar, Keaveny, Zhang and Shelley2012), and an increase of a factor of two to three times in unstructured media (Kamal & Keaveny Reference Kamal and Keaveny2018). By contrast, our results here are significantly more modest. We attribute this more modest swimming speed enhancement to the freedom of movement afforded by three dimensions, and to a difference in the fixing of the obstacles. The previous work (Majmudar et al. Reference Majmudar, Keaveny, Zhang and Shelley2012; Kamal & Keaveny Reference Kamal and Keaveny2018) showed that planar motion was primarily governed by steric contact: motion could be hindered by trapping, and could be promoted by the provision of obstacles for the swimmer to push against. In contrast, 3-D motion allows for escape from constrained environments, diluting the ability for the structure to hinder or promote motion directly. In the earlier work, the obstacles are fixed to a point in space, allowing the forces on them to not be completely resolved by the fluid. Here, conversely, the forces on the suspended filaments are completely balanced by the surrounding fluid. Finally, we have concluded that speed enhancement in the current work is instead linked to hydrodynamic interactions with neighbouring filaments; an increased number of neighbours is another property of 3-D geometry.
We also see less evidence of frequent or sharp turns in the swimmer trajectories when compared with planar motion, resulting in swimming paths that are more linear and less diffusive. We note, however, that our simulations of $20T$ are well below the correlation time, $\tau _c$, which we estimate by fitting $\langle \hat {\boldsymbol {P}}(s)\boldsymbol {\cdot }\hat {\boldsymbol {P}}(s+t) \rangle$ with $\exp (-t/\tau _c)$. We find values of $\tau _c$ between $110T$ and $3500T$ in our parameter regime. Nonetheless, we find trends which agree with the planar study: increasing the concentration and stiffness of network filaments results in greater changes in swimmer orientation, reflecting that the swimmer is deflected more by larger numbers of suspended filaments.
The motivation for our work was to model a complex fluid, such as mucus, as an immersed 3-D filament network immersed in as viscous fluid, and elucidate the effects of physical interactions with the network on swimmer motion, especially with regard to the conditions under which this model recovers behaviour seen in viscoelastic experiments and in continuum models. Our results on swimming speeds, with mild speed enhancement at low concentrations of flexible networks but hindrance in stiff, concentrated networks, display the same trends as results obtained using 2-D and 3-D viscoelastic continuum models when the Deborah number or an analogous measure of resistance is increased. Interestingly, in these studies, either the waveform amplitude is relatively fixed (Teran et al. Reference Teran, Fauci and Shelley2010; Thomases & Guy Reference Thomases and Guy2014), or can correlate faster speeds with increased (free to vary) amplitudes (Olson & Leiderman Reference Olson and Leiderman2015; Leiderman & Olson Reference Leiderman and Olson2016); whereas the speed enhancement we observe happens despite a decrease in the swimming stroke amplitude. Along with these studies based on continuum descriptions, our results are in line with those obtained examining undulatory swimming through highly connected spring networks (Wróbel et al. Reference Wróbel, Lynch, Barrett, Fauci and Cortez2016). This study, where the swimmer waveform is prescribed and the network interact with the fluid only through the nodes, speed increases of approximately 20 % were also reported. We see a similar speed increases, especially in the more flexible networks, where the swimmer waveform is mainly unaffected by the suspended filaments, as in Wróbel et al. (Reference Wróbel, Lynch, Barrett, Fauci and Cortez2016).
The density and spacing of the network filaments in the simulations presented here were inspired by that found in mucin networks. Although we have considered a large range of concentrations, typical mucin concentration in mucus is between 0.2 % and 5.0 % v/w; water forms approximately 95 % (Gniewek & Kolinski Reference Gniewek and Kolinski2010; Leal, Smyth & Ghosh Reference Leal, Smyth and Ghosh2017). And in this regime, we have found mild speed promotion for an undulatory swimmer in both flexible and stiff filament environments, with a small further improvement with stiffer springs in the networks. This is consistent with positive sperm migration in human cervical mucus in the ovulatory phase of the menstrual cycle, when hydration increases to 98 % w/w and the mucin fibres can appear globular and disconnected (Bergman Reference Bergman1953; Morales, Roco & Vigil Reference Morales, Roco and Vigil1993; Brunelli et al. Reference Brunelli, Papi, Arcovito, Bompiani, Castagnola, Parasassi, Sampaolese, Vincenzoni and De Spirito2007; Curlin & Bursac Reference Curlin and Bursac2013). The luteal phase, correlating with highly reduced migration, sees hydration decreasing to 90 % w/w and the fibres appearing densely connected (Katz, Slade & Nakajima Reference Katz, Slade and Nakajima1997; Brunelli et al. Reference Brunelli, Papi, Arcovito, Bompiani, Castagnola, Parasassi, Sampaolese, Vincenzoni and De Spirito2007). In this less dilute regime, we begin to see speed reduction in stiff filament networks, which is enhanced by stiffer network connecting springs. It is important to note that while the chemical composition of mucins and polymers are well known (Leal et al. Reference Leal, Smyth and Ghosh2017), their mechanical properties at an individual fibre level, as well as the mechanical properties of their connections, are less well understood (Lai et al. Reference Lai, Wang, Cone, Wirtz and Hanes2009).
Classes of microswimmers are able to propel themselves through different locomotive styles: here we have found that low-to-mid concentrations of flexible networks are mildly favourable towards the undulatory motion of microorganisms such as sperm cells or C. elegans. Another significant locomotive style is to rotate helical flagella, seen in Pseudomonas aeruginosa or E. coli. Simulations of this motion, resolving the fluid using alternative methods such as multiparticle collision dynamics or regularised Stokeslets have seen speed promotion up to 60 % in confined tubes (LaGrone, Cortez & Fauci Reference LaGrone, Cortez and Fauci2019) and polymer solutions (Zöttl & Yeomans Reference Zöttl and Yeomans2019). In the latter case, speed up is attributed to chirality of the flagellum, which is absent in undulatory swimmers. This is interesting because one purpose of biological networks of the kind modelled here, such as mucus, is to block pathogens. We might expect that these networks are favourable to certain locomotive styles over others, and it remains as a future direction to determine this within the framework presented here.
Supplementary material
Supplementary material are available at https://doi.org/10.1017/jfm.2024.603.
Funding
The authors gratefully acknowledge support from EPSRC Grant EP/P013651/1.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Rigid body motion
In order to assess the effectiveness of the swimming gait that emerges in the filament environment, we determine the swimming speed produced by the gait in a viscous fluid absent of any filaments.
A.1. Extracting the swimming gait
Suppose that in the filament suspension, the swimmer's centre of mass has velocity $\boldsymbol {u}^{sus} = {\frac {{\rm d} \boldsymbol {Y}_{CoM}}{{\rm d} t}}$ and angular velocity $\boldsymbol {\omega }^{sus}$. Note that we use lowercase letters to represent the rigid body motion (vectors of length 3 in three dimensions) and uppercase letters for the segments of the swimming filament (indexed by $n$); the notation ‘${sus}$’ indicates these data are from the suspension simulation. The motion of segment $n$ can then be written as
where $\Delta \boldsymbol {Y}_n = \boldsymbol {Y}_n - \boldsymbol {Y}_{CoM}$, $\tilde {\boldsymbol {V}}_n$ and $\tilde {\boldsymbol {\varOmega }}_n$ are the segment velocities and angular velocities associated with swimmer deformation for swimming.
To extract $\tilde {\boldsymbol {V}}_n$ and $\tilde {\boldsymbol {\varOmega }}_n$, we need to find expressions for $\boldsymbol {u}^{sus}$ and $\boldsymbol {\omega }^{sus}$ in terms of the known $\boldsymbol {V}_n^{sus}$ and $\boldsymbol {\varOmega }_n^{sus}$. Summing (A1) over $n$, we find
since we must have $\sum _n\tilde {\boldsymbol {V}}_n = \boldsymbol {0}$ and $\sum _n\Delta \boldsymbol {Y}_n = \boldsymbol {0}$.
Now we consider taking the cross product of $\Delta \boldsymbol {Y}_n$ and (A1), and summing
where, on the right-hand side,
• term 1 is zero by the definition of $\tilde {\boldsymbol {V}}$;
• term 2 is zero since $\sum _n\Delta {\boldsymbol {Y}}_n = \boldsymbol {0}$;
• term 3 introduced the moment of inertia matrix, defined in terms of a skew-symmetric matrix $[\boldsymbol {a}]_\times$ such that $[\boldsymbol {a}]_\times \boldsymbol {b}=\boldsymbol {a}\times \boldsymbol {b}$
(A6a,b)\begin{equation} {{\boldsymbol{\mathsf{I}}}}_M ={-}\sum_n [\Delta\boldsymbol{Y}_n]_\times[\Delta\boldsymbol{Y}_n]_\times, \quad [\boldsymbol{a}]_\times := \begin{pmatrix}0 & -a_3 & a_2 \\ a_3 & 0 & -a_1 \\-a_2 & a_1 & 0\end{pmatrix}.\end{equation}
Hence
We can therefore now rearrange (A1), (A2) to find an expression for the local motion, $\tilde {\boldsymbol {V}}_n$ and $\tilde {\boldsymbol {\varOmega }}_n$, as a sum of fully evaluated terms
We can now take these expressions as representing the gait of the swimmer, and place them into an empty fluid.
A.2. Resulting swimming speed in a viscous fluid
Now consider the swimmer in a viscous fluid free from any suspended filaments. The mobility problem is, as in (2.11),
where we have dropped the $H$ for ‘hydrodynamic’ for notational convenience. Recall that $\mathcal {M}$, for a single swimming filament formed of $N$ segments, has size $6N \times 6N$ in three dimensions. Once again we can decompose the motion into a rigid body motion and a local deformation for each segment $n$
where $\tilde {\boldsymbol {V}}_n$ and $\tilde {\boldsymbol {\varOmega }}_n$ are specified from (A8) and the unknowns are $\boldsymbol {F}_n$, $\boldsymbol {T}_n$, $\boldsymbol {u}$ and $\boldsymbol {\omega }$. We solve this system for the unknowns, subject to the conditions that the total force and torque are zero
Equations (A10)–(A12) can be combined to form a system of equations with a symmetric matrix of size $(6N+6)\times (6N+6)$
recalling that $\tilde {\boldsymbol {V}}^\mathsf {T} = (\tilde {\boldsymbol {V}}_1^\mathsf {T},\ldots,\tilde {\boldsymbol {V}}_N^\mathsf {T})$ and similarly for $\tilde {\boldsymbol {\varOmega }}$, $\boldsymbol {F}$ and $\boldsymbol {T}$. Here, ${{\boldsymbol{\mathsf{I}}}}$ is the identity matrix and $[\Delta \boldsymbol {Y}]_\times$ is again the skew-symmetric matrix corresponding to the cross product from (A6a,b), in both cases suitably stacked for each segment $n$. Solving this system gives the swimming velocity in the viscous fluid.
Appendix B. Resistive force theory
Resistive force theory (Johnson & Brokaw Reference Johnson and Brokaw1979; Lauga & Powers Reference Lauga and Powers2009) can be used to provide the mobility matrix in (2.11). In RFT, the filament segments are hydrodynamically uncoupled, with the force and torque on a segment dependent only on the segment's orientation, velocity and angular velocity. In this way, it can be thought of as a drag-based model, but one where the force–velocity drag coefficients are different in the directions normal and tangent to the segment orientation.
In the mobility formulation of RFT, the velocity and angular velocity of segment $n$ are given by
where $\alpha _{||}$, $\alpha _{\perp }$, $\beta$ are mobility coefficients for the segments. The torque coefficient $\beta$ is simply that for a sphere in Stokes flow
while $\alpha _{||}$ and $\alpha _{\perp }$ are initially estimated using FCM simulations of straight filaments, following the approach of Schoeller & Keaveny (Reference Schoeller and Keaveny2018, supplementary material). These two parameters are then tuned so that the swimming speed and maximum beating amplitude match the full FCM simulation for a single swimmer in an empty fluid. Doing so, we choose the numerical values
A comparison of the resultant waveform for FCM and RFT is given in figure 11. The ratio $\alpha _{||}/\alpha _{\perp } \approx 1.44$ is similar to previous work (Schoeller & Keaveny Reference Schoeller and Keaveny2018). While for an infinitely thin filament, we expect a ratio of 2 (Johnson & Brokaw Reference Johnson and Brokaw1979), we observe that as the swimmer has an aspect ratio of $2a/L \approx 0.06$, we are not in this limiting regime. Note that RFT can also be presented in the resistance formulation, where this ratio is presented in terms of drag coefficients, $\xi _\perp /\xi _{||}$.