3.1 Reciprocal relations

The boundary integrals associated with the boundary element method are constructed via the reciprocal relation, which imposes a constraint on two sets of solutions $(\boldsymbol{u},\unicode[STIX]{x1D748})$ and $(\boldsymbol{u}^{\prime },\unicode[STIX]{x1D748}^{\prime })$ for the fluid equations within the same domain but with different boundary conditions. For Stokes flow, the reciprocal relation is (Pozrikidis Reference Pozrikidis1992)

with $\unicode[STIX]{x1D749}=2\unicode[STIX]{x1D707}\boldsymbol{D}$ . However, this no longer holds for linear viscoelastic media such as the Maxwell fluid we are considering, where ${\mathcal{L}}\unicode[STIX]{x1D749}=2\unicode[STIX]{x1D707}\boldsymbol{D}$ .

With the latter deviatoric stress, ${\mathcal{L}}\unicode[STIX]{x1D749}=2\unicode[STIX]{x1D707}\boldsymbol{D}$ , we consider ${\mathcal{L}}\unicode[STIX]{x1D70E}_{ij}$ instead of $\unicode[STIX]{x1D70E}_{ij}$ , and derive a similar reciprocal relation for a linearly viscoelastic fluid, given below by (3.9). The derivation is directly analogous to that presented in Pozrikidis (Reference Pozrikidis1992) for Newtonian flows, where further details may be found. To proceed, first consider the identity

Focusing on the second term of the right-hand side above, note that the pressure term in the stress tensor does not contribute due to fluid incompressibility, so that

where the last equality holds from the constitutive relation (2.3).

Hence we have

By analogy, we have the symmetric relation

Here, noting the definition of $D_{ij}$ and relabelling by $i\rightarrow j$ and $j\rightarrow i$ , we have

yielding the relation

The left-hand side of this relation, however, vanishes from the momentum-balance equations as

Hence we finally obtain the reciprocal relation for linearly viscoelastic flow,

As with the Newtonian case (Pozrikidis Reference Pozrikidis1992), this reciprocal relation provides the keystone for the boundary element numerical methods, as detailed below.

3.2 Viscoelastic Stokeslet

In this subsection, we consider a Stokeslet for the linear viscoelastic fluid $G_{ij}$ , which is the Green function of this linear system and thus the velocity field when a point force at $\boldsymbol{x}_{0}$ is applied, defined as ${\mathcal{L}}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{\prime })_{ij}=-\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{0})$ , or

where $P_{j}$ is the associated pressure and the singularity location, $\boldsymbol{x}_{0}$ , in general may be a function of time. We can therefore express the Stokeslet for the linear viscoelastic fluid by the solutions for the Newtonian Stokes flow,

with the superscript $S$ referring to Stokes flow.

Let the constant vector $\boldsymbol{q}$ be a force strength, and introduce a set of solutions for the linear viscoelastic fluid flow,

whereby we have

This is equivalent to

where $\unicode[STIX]{x1D70E}_{ij}^{\prime S}$ is the stress tensor for the Newtonian Stokes flow, given by $\unicode[STIX]{x1D70E}_{ij}^{\prime S}=(1/8\unicode[STIX]{x03C0})T_{ijp}^{S}q_{p}$ . Thus (3.14) can be rephrased as

In principle there is an undetermined freedom regarding the initial condition. We thus assume the flow is stationary at initial time, $t_{0}$ , and then in turn $\unicode[STIX]{x1D70E}_{ij}^{\prime }$ is uniquely specified as the operator ${\mathcal{L}}$ is linear.

More generally, note that, in developing a boundary element method, one is constructing the flow field from the linear superposition of singularity solutions (Pozrikidis Reference Pozrikidis1992), and the choice of singularity solutions is the same as choosing a basis for a representation of the velocity field. Thus here we are using Newtonian velocity singularities, since $G_{ij}=G_{ij}^{S}$ , but the associated singularity pressure and stress fields are now more complex, and given by ${\mathcal{L}}P_{j}=P_{j}^{S}$ , ${\mathcal{L}}T_{ijp}=T_{ijp}^{S}$ , to ensure the reciprocal theorem holds. The latter in turn enables the construction of the boundary element method to find the decomposition of the flow field into singularity solutions even with linear viscoelasticity. However, the pressure and stress fields do not enter the final formalism and thus the net influence of linear viscoelasticity only significantly arises in the force and torque balances for the swimmer, which are, at most, modified in a simple manner. In turn, this enables a rapid adaptation of standard boundary element codes to study swimming in linearly viscoelastic media, as we shall derive and detail below.

3.3 Boundary integral equations

Based on the relations derived above, we are in a position to formulate the viscoelastic boundary integral equations, starting with the reciprocal relation (3.9). Substitution of the velocity and the stress of the point force flow $(\boldsymbol{u}^{\prime },\unicode[STIX]{x1D748}^{\prime })$ , obtained in (3.11) and (3.15), into the reciprocal relation gives us the relation for the velocity and the stress of our swimmer problem $(\boldsymbol{u},\unicode[STIX]{x1D748})$ as

This relation is integrated over the fluid region except for a small sphere of radius $\unicode[STIX]{x1D716}$ around the singular point $\boldsymbol{x}_{0}$ , $\unicode[STIX]{x1D6FA}\backslash \unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D716}}$ , as depicted in figure 1. Noting the zero-velocity constraint on the wall boundary and at infinity, we have

where the vector $\boldsymbol{n}_{j}$ denotes the outward normal of the surface region, as in figure 1. The time derivative implicit within the operator ${\mathcal{L}}$ on the swimmer surface $S$ should be clarified. We take a larger region that covers the surface $S$ , with separation $\unicode[STIX]{x1D6FF}>0$ from it, and consider the time derivatives for the variables outside $S$ but in the larger region, where the derivative is unambiguously defined, before taking the limit of $\unicode[STIX]{x1D6FF}\rightarrow 0$ .

Exactly as in the Newtonian case (Pozrikidis Reference Pozrikidis1992, p. 26), the integral over $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{\unicode[STIX]{x1D716}}$ converges to $8\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}u_{p}(\boldsymbol{x}_{0})$ in the limit of $\unicode[STIX]{x1D716}\rightarrow 0$ , yielding the viscoelastic boundary integral equation,

The first of these integrals is known as a single-layer boundary integral, while the second is known as a double-layer boundary integral. The latter requires care in its treatment, as the integrand is singular, and not just weakly singular. The evaluation of these integrals in the Newtonian case is fully documented (e.g. Pozrikidis Reference Pozrikidis1992) and immediately generalises to (3.18), as the double-layer singularity, $T_{ijk}^{S}$ , is the Newtonian singularity. However, we do not follow this approach here. Instead we explicitly assume the swimmer volume is invariant. Then it is much simpler to eliminate the second integral of (3.18), and write $u_{p}(\boldsymbol{x}_{0})$ in terms of only a single-layer boundary integral, as is well known for the Newtonian boundary element method (Pozrikidis Reference Pozrikidis1992) and as we now demonstrate for linear viscoelasticity.

To proceed, let $^{\ast }$ denote variables associated with a fictitious virtual flow within the interior of the swimmer $\unicode[STIX]{x1D6FA}_{s}(t)$ on making the supposition that the swimmer interior is constituted from the same linear viscoelastic medium as external to the swimmer. For example, let $\boldsymbol{u}^{\ast }$ denote the velocity field of this internal swimmer flow; its incompressibility immediately requires the restriction that swimmer volume is invariant. Given this, the same reciprocal relation (3.9) then holds in the interior region $\unicode[STIX]{x1D6FA}_{s}$ , and integrating over this region yields

after using the divergence theorem and continuity of velocity at the cell surface. The vector $\boldsymbol{n}_{j}$ again denotes the outward normal to the swimmer surface. We use the solutions $(\boldsymbol{u}^{\prime },\unicode[STIX]{x1D748}^{\prime })$ and $(\boldsymbol{u}^{\prime \ast },\unicode[STIX]{x1D748}^{\prime \ast })$ for a point force at $\boldsymbol{x}_{0}$ outside $\unicode[STIX]{x1D6FA}_{s}$ , and substituting (3.19) into (3.18) we have the single-layer integral equation,

We therefore find the same expression for the single-layer integral equation as in the Newtonian case with

giving the linear viscoelastic generalisation of the jump in the traction between the exterior and interior flows. When we take the limit of $\boldsymbol{x}_{0}$ approaching the swimmer surface $S$ , we obtain the single-layer boundary integral equation for the surface velocity,

3.4 Force and torque balance equations

For the swimmer problem, we seek the linear and angular velocities $\boldsymbol{U}$ and $\unicode[STIX]{x1D734}$ of the swimmer when its deformation is given. In addition to the no-slip boundary condition on the swimmer surface, equation (3.22), the equations of motion for the swimming body are required to close the system. At negligible Reynolds number, these reduce to the conditions of zero force and zero torque on the swimmer, given there are no external forces and torques. The zero-force condition,

needs to be expressed in terms of the viscoelastic traction jump, $\boldsymbol{f}^{VE}$ . From the momentum balance (2.1) and the divergence theorem, we have

with the temporal operator in ${\mathcal{L}}$ commuting with integration on the right-hand side, as the integration domain is time-invariant. Using the momentum-balance relation and the divergence theorem once more gives the general commutation of the surface integral and the linear operator  ${\mathcal{L}}$ :

Furthermore, note that for the virtual interior flow

and hence we obtain the relation

The zero-force condition (3.23) therefore entails

for all time.

The zero-torque condition

is rewritten in a similar manner. Applying the operator ${\mathcal{L}}$ to this integral together with the divergence theorem gives us

The integrand of the volume integral vanishes due to the momentum balance and the symmetry of the stress tensor, with the same reasoning applying to the internal flow relationship:

Thus the integral (3.30) reduces to

which yields the zero-torque condition

Equations (3.22), (3.28) and (3.33) complete the swimmer problem, which gives us the swimmer velocities, $\boldsymbol{U}$ and $\unicode[STIX]{x1D734}$ . This set of equations has the same form as the Newtonian case, though the corresponding tractions are not identical. Thus the swimmer velocities, and therefore the swimmer trajectories, are identical to those in the Newtonian Stokes flow for any finite prescribed shape deformation, regardless of the linear viscoelasticity of the Maxwell medium. Hence Purcell’s scallop theorem holds for finite-deformation swimming in a linear Maxwell fluid.

4.1 Wall interaction

When a swimmer is less than submicrometre length scales from a boundary, we can no longer neglect molecular-based swimmer–wall interactions. While there are no comprehensive data on such interactions, it has been reported that the wall interaction of Escherichia coli can be attractive or repulsive according to the solutes present (Klein et al. Reference Klein, Clapp and Dickinson2003) and on a length scale of approximately 100 nm. One also regularly observes that sperm easily stick to naked glass, requiring the use of dilute albumin to enable swimming adjacent to coverslips (Smith et al. Reference Smith, Gaffney, Gadêlha, Kapur and Kirkman-Brown2009b ). Also, particular proteins distributed on the sperm head can bind to biological surfaces such as the epithelium of the mammalian oviduct and the zona pellucida (Okabe Reference Okabe2013; Suarez Reference Suarez2016).

Hence, as used in previous studies (Ishimoto et al. Reference Ishimoto, Cosson and Gaffney2016; Ishimoto & Gaffney Reference Ishimoto and Gaffney2016) and assuming conditions that the swimmer does not bind to the surface, we introduce a repulsive interaction potential for a repulsive force per unit area between the wall and the swimmer surface at the position $\boldsymbol{x}$ . This is given via

where $r$ is the distance from the surface and $\boldsymbol{e}_{r}$ denotes the outward normal of the wall, which we consider hereafter as a no-slip rigid infinite plane boundary, and the constants $g$ and $d$ represent the force strength and the force length scale, respectively.

The force and torque balance equations, (3.28) and (3.33), thus need to be modified to accommodate this wall interaction. The force balance is given by

Applying the operator ${\mathcal{L}}$ to both sides, and using (3.25) and arguments analogous to those presented in the previous section, we have

The torque balance equation is similarly obtained, to give

Note that the singularity solution used in the boundary integral equation (3.18) or (3.22) satisfies appropriate boundary conditions of no slip on $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{w}$ and zero for $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{\infty }$ . Hence the integral kernel in the presence of a wall, $G_{ip}^{S}$ , is that of the Blakelet (Blake Reference Blake1971), which we employ for numerical computation below. Finally, observe that the structure of the equations is that of the standard, Stokes flow, boundary element method, enabling rapid generalisation of existing boundary element numerical libraries for linear viscoelastic calculations.

4.2 Numerical implementation

Here we use the single-layer formulation for microswimming problems (Ishimoto & Gaffney Reference Ishimoto and Gaffney2014a ,Reference Ishimoto and Gaffney b ) to consider deformation swimming in a linear Maxwell fluid with a constitutive relation given by ${\mathcal{L}}\unicode[STIX]{x1D749}=2\unicode[STIX]{x1D707}\boldsymbol{D}$ . This requires solving the linear set of equations (3.22), (4.3) and (4.4), in which the traction $\boldsymbol{f}^{VE}(\boldsymbol{x})$ and the velocities $\boldsymbol{U}$ and $\unicode[STIX]{x1D734}$ are unknown but the swimmer deformation $\boldsymbol{v}$ is given. The only difference from the Newtonian case is the additional operator ${\mathcal{L}}$ in the force and torque balance equations (4.3) and (4.4). The swimmer surface $S$ is discretised by $N$ triangular meshes, and the surface integrals are numerically performed by Gaussian quadratures supported by the BEMLIB library (Pozrikidis Reference Pozrikidis2002). The resulting $3N+6$ coupled linear equations are then solved via an LU decomposition. The time marching of the system is implemented by the Heun method (Smith et al. Reference Smith, Gaffney, Blake and Kirkman-Brown2009a ). As exemplars, we consider a spherical tangential squirmer, a bacterial swimmer and a sperm-like swimmer, as detailed in the following sections.

4.3 Spherical tangential squirmer

The first example swimmer is a spherical squirmer, which is driven via a tangential surface velocity profile on its surface. In particular, we consider a sphere of radius $a$ with its centre a height $h$ above an infinite plane wall at $z=0$ , and introduce a swimmer orientation vector $\boldsymbol{d}$ within the $xz$ plane, making an angle $\unicode[STIX]{x1D711}$ with the plane wall; see figure 2. The surface velocity profile is given by

where the angle $\unicode[STIX]{x1D703}$ is measured from the orientation vector $\boldsymbol{d}$ , as also depicted in figure 2. This swimmer is also known as a neutral squirmer (Ishikawa, Simmons & Pedley Reference Ishikawa, Simmons and Pedley2006), and is repelled by the wall in the absence of wall interactions (Ishimoto & Gaffney Reference Ishimoto and Gaffney2013).

Figure 2.

(a) A schematic of the spherical squirmer. (b) Example swimmer trajectories near a wall with an enlarged view shown in the inset. (c) Relation of the initial and final scattering angles for wall scattering, including Newtonian theoretical predictions for shallow scattering.

The sphere radius and the surface velocity are non-dimensionalised to a unit radius and a non-dimensional velocity $V=3/2$ . Owing to this scaling, the non-dimensional swimming velocity in the absence of a wall becomes unity. The number of mesh elements on the swimmer surface is $N=2048$ , and the meshes are iteratively generated by the BEMLIB library accompanying Pozrikidis (Reference Pozrikidis2002). The discretisation for each time step is given by $\unicode[STIX]{x0394}t=1/400$ , and the constants of the wall interactions are set as $(g\unicode[STIX]{x1D707}/T,d/a)=(400,1/50)$ . The large scale of the interaction strength entails the repulsive force is very large once the swimmer is in the vicinity of the wall, while the small interaction length scale entails that the swimmer has to be very close to the wall to experience the repulsive force.

We performed numerical computations, changing the initial angle of attack and the Deborah number, which is the non-dimensional relaxation time, given by $De=\unicode[STIX]{x1D706}/T$ , where $T=3a/2V$ is the characteristic time scale of this system. In figure 2(b), the computed trajectories of the squirmer are plotted for an initial angle of attack $\unicode[STIX]{x1D711}_{init}=0.2\unicode[STIX]{x03C0}$ with the Deborah numbers $De\in \{0,1,10\}$ . It is found that the viscoelasticity weakly increases the surface repulsion and weakly influences the scattering angle and the duration that the swimmer remains near the wall. In particular, on varying the initial angles of attack for these swimmers, the relation between the initial and final scattering angles $(\unicode[STIX]{x1D711}_{init},\unicode[STIX]{x1D711}_{fin})$ is plotted in figure 2(c). In these simulations, the swimmer is initially located at $(x,z)=(0,5)$ and the final angles are obtained from the trajectory only once the simulation times reaches $50T$ , so that the swimmer is in the very far field of the wall. In the plots of figure 2(c), for the cases with shallow initial angles, the swimmer does not approach the area with wall interactions, and thus the scattering angles are well explained by the far-field point singularity theory,

provided by Spagnolie & Lauga (Reference Spagnolie and Lauga2012). Nonetheless, as the initial angle steepens, the final angles deviate from the shallow-angle theoretical prediction and vary non-monotonically with the Deborah number.

4.4 Bacterial swimmer

Figure 3.

(a) A schematic of the monotrichous bacterial swimmer with its flagellum, together with the axes of the inertial and body-fixed reference frames and the underlying no-slip wall. (b) The time evolution of the separation distance from the proximal flagellum to the wall, with differences of $z$ from the $De=0$ case depicted for $De=0.1$ and $De=1$ in the inset.

The second example is a monotrichous bacterial swimmer, which is known to accumulate near a no-slip wall (Berke et al. Reference Berke, Turner, Berg and Lauga2008) and numerical simulations have shown that hydrodynamic interactions are sufficient to explain even the subtle aspects of this accumulation behaviour (Giacché, Ishikawa & Yamaguchi Reference Giacché, Ishikawa and Yamaguchi2010; Goto, Nakai & Aoki Reference Goto, Nakai and Aoki2010; Shum, Gaffney & Smith Reference Shum, Gaffney and Smith2010). Hence this example examines the extent to which viscoelastic effects may alter stable swimming near a surface, as opposed to the scattering of the previous example.

The virtual bacterium consists of a spheroidal cell body, whose polar radius and equatorial radius are, respectively, $a_{1}$ and $a_{2}$ , together with a single cylindrical flagellum of length $L$ and radius $r$ (figure 3 a). Following Shum et al. (Reference Shum, Gaffney and Smith2010), the flagellum waveform is taken to be a helical wave, and the centreline of the flagellum is given in terms of the body-fixed coordinates $(\hat{x}(\hat{z},t),{\hat{y}}(\hat{z},t),\hat{z})$ by

where $A$ is the amplitude of the helix, $k$ is the wavenumber and $\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}/T$ is the angular velocity in terms of the time period $T$ . The parameter $k_{E}$ governs the scale of the initial increase in the flagellar envelope as the flagellum leaves its junction with the cell body at $\hat{\boldsymbol{x}}=\mathbf{0}$ . The functions (4.7) and (4.8) are defined for $\hat{z}\in [\hat{z}_{end},0]$ , where the coordinate of the flagellar distal end $\hat{z}_{end}<0$ is given by

In the simulation, the parameters of the model bacterium are those of optimal power efficiency, as obtained by Shum et al. (Reference Shum, Gaffney and Smith2010) and given by $a_{1}/L=0.154$ , $a_{2}/L=0.0922$ , $k_{E}L=10.0$ , $A/L=0.079$ , $kL/2\unicode[STIX]{x03C0}=1.953$ and $r/L=0.003$ . The system is non-dimensionalised using the flagellum length for the length scale $L$ , the flagellum rotation period for the time scale $T$ and the medium viscosity for the force scale, such that the non-dimensional viscosity is unity. Note however that, for the virtual bacterial model discussed by Shum et al. (Reference Shum, Gaffney and Smith2010), the flagellum rotates with a constant applied torque at its proximal end. In this study, in contrast, the flagellum motion is fully prescribed in the body-fixed frame by the functions (4.7) and (4.8), and the bacterial motion is therefore subject to the same kinematical problem as the squirmer case.

The surface meshes are generated as in Shum et al. (Reference Shum, Gaffney and Smith2010), with the head mesh constructed in the same way as the squirmer, with a subsequent stretching to generate the ellipsoidal head shape. The flagellum is meshed as a hexagonal tube with hemispherical caps at the ends, noting a small cell–flagellum gap is introduced in practice to numerically accommodate the relative motion of these surfaces as in previous studies (Shum et al. Reference Shum, Gaffney and Smith2010; Ishimoto & Gaffney Reference Ishimoto and Gaffney2014a ). The total number of mesh elements on the swimmer surface is set to $N=494$ , and the discretisation for each time step is given by $\unicode[STIX]{x0394}t=1/400$ . The constants of the wall interactions are once more set to $(g\unicode[STIX]{x1D707}/T,d/L)=(400,1/50)$ .

We conducted numerical simulations with different Deborah numbers $De=\unicode[STIX]{x1D706}/T\in \{0,0.1,1\}$ . Note that at high Deborah numbers the forces and torques acting on the flagellum can differ substantially relative to the Newtonian case, and therefore the flagellar motion may change even if the internal torque generation is invariant. In such a highly elastic medium, the prescribed motion is in fact likely to be altered. In contrast to the previous idealised example of a squirmer, which we assumed can always maintain the prescribed deformation even with large elastic forces, here we focus on moderate values of the Deborah number. In figure 3(b), the time evolution of the distance from the wall boundary and the flagellar proximal end is plotted, and the differences are relatively minor, as particularly emphasised by the inset. The swimmer is initially located so that the flagellar proximal end is at $\boldsymbol{x}=(0,0,0.5)$ and its initial angle of attack relative to the wall is $0.2\unicode[STIX]{x03C0}$ . After the swimmer reaches the region where the wall interaction is significant, the cell body rotates away from the surface and the cell subsequently relaxes to stable swimming adjacent to the wall.

4.5 Sperm-like swimmer

The third example swimmer is based on the morphology of a sperm, as illustrated in figure 4(a). Sperm in particular can be significantly affected by steric interactions when a sudden change in surface topography induces scattering (Kantsler et al. Reference Kantsler, Dunkel and Goldstein2013). In addition, it is theoretically predicted that steric interactions significantly influence the behaviour of a sperm cell near a boundary (Montenegro-Johnson, Gadelha & Smith Reference Montenegro-Johnson, Gadelha and Smith2015; Ishimoto & Gaffney Reference Ishimoto and Gaffney2015, Reference Ishimoto and Gaffney2016) where, in the latter citation, we examine the impact of viscoelasticity on beating cells adhered to a wall. Hence, here we consider how viscoelasticity impacts sperm swimming trajectories near a wall, which necessitates a viscoelastic boundary element method in the absence of the simplification associated with cell attachment.

Figure 4.

(a) A schematic of the sperm-like swimmer with its flagellum above a no-slip plane boundary, together with body-fixed and inertial reference frame axes. (b) The time evolution of the separation distance of the cell body–flagellum junction from the wall, with differences of $z$ from the Newtonian case shown for $De=0.1$ and $De=0.5$ in the inset.

The virtual sperm-like swimmer has a flattened ellipsoidal head of human sperm and a cylindrical tail flagellum of length $L$ and radius $r$ . The head and tail morphologies are motivated by the human sperm geometrical model provided in Ishimoto & Gaffney (Reference Ishimoto and Gaffney2016), with a prescribed flagellar waveform that is an ellipsoidal helical wave (Smith et al. Reference Smith, Gaffney, Blake and Kirkman-Brown2009a ), with an ellipsoidal cone envelope. The centreline of the flagellum is again defined as a function of $\hat{z}\in [\hat{z}_{end},0]$ , via

where the parameter $\unicode[STIX]{x1D6FC}$ represents the major to minor aspect ratio of the elliptical cross-section through the enveloping cone of the flagellar beat, and $\hat{z}_{end}<0$ satisfies the fixed length constraint given by (4.9). The parameters used in the simulation are taken from Ishimoto & Gaffney (Reference Ishimoto and Gaffney2016), where morphological justification is provided; in particular we have $\unicode[STIX]{x1D6FC}=0.2$ , $A/L=0.2$ , $kL/2\unicode[STIX]{x03C0}=1.5$ and $r/L=0.0223$ . As with the bacterial swimmer, the system is non-dimensionalised using the flagellum length $L$ for the length scale, the beat period $T$ for the time scale and the medium viscosity for the force scale, such that the non-dimensional viscosity is unity.

The surface meshes are generated in a similar manner to the bacterial swimmer, details of which are provided in our previous studies (Ishimoto & Gaffney Reference Ishimoto and Gaffney2014a , Reference Ishimoto and Gaffney2016). The total number of mesh elements on the swimmer surface is again $N=494$ , and the discretisation for each time step is given by $\unicode[STIX]{x0394}t=1/1600$ . The constants of the wall interactions are set to $(g\unicode[STIX]{x1D707}/T,d/L)=(400,1/50)$ while the swimmer is initially located so that the head–flagellum junction is at $\boldsymbol{x}=(0,0,0.5)$ and its angle of attack to the wall is  $0.2\unicode[STIX]{x03C0}$ .

The results of similar numerical simulations with different Deborah numbers $De=\unicode[STIX]{x1D706}/T\in \{0,0.1,0.5\}$ are shown in figure 4(b) at integer times, or every beat cycle, with enlarged views inset. As with the bacterium, we do not consider large elasticity, because the prescribed waveform is unlikely to be conserved; instead, we consider moderate Deborah numbers once more. In figure 4(b), the time evolution of the distance between the wall and the head–flagellum junction is plotted, and the differences are relatively minor here too, as emphasised by the inset. Similarly to the bacterial case, the swimmer rotates on encountering the wall interactions, but the sperm head is slightly directed towards the wall. Furthermore, when the cell trajectory has asymptoted to boundary swimming, the conical envelope of the flagellum has rotated sufficiently downwards so as to induce intermittent wall interactions during the beat period. However, despite the fact that these interactions differ in the contribution to the momentum balances for the cell according to medium rheology, the actual trajectory of the sperm is not significantly altered on comparing swimming in a Newtonian fluid to that in a linear Maxwell fluid.

5.1 Model

The magnetic swimmer is taken to be a rigid helix of length $L$ and radius $r$ , as considered in the section discussing the bacterial swimmer and depicted in figure 3. The centreline of the helix is given in terms of the body-fixed coordinates $(\hat{x}(\hat{z}),{\hat{y}}(\hat{z}))$ by

where $A$ is the radius of the helix and $k$ is the wavenumber, with the latter giving the pitch of the helix. The $\hat{z}$ coordinate ( $\hat{z}\in [\hat{z}_{end},0]$ ) is related to the arclength via the fixed length constraint given by (4.9).

We assume the swimmer is a hard magnet with a fixed magnetisation $\boldsymbol{M}=M_{0}\hat{\boldsymbol{e}}_{x}$ along the body-fixed $\hat{x}$ axis. Therefore the helix experiences a force and torque in the presence of a magnetic flux density $\boldsymbol{B}$ given by (Qiu & Nelson Reference Qiu and Nelson2015)

where $C_{0}$ is the constant magnetic volume of the object. We further assume the external magnetic flux density is that of a uniform rotating magnetic field given with respect to the inertial-frame basis via

where $B_{0}$ is constant and $T$ is the period of rotation. Note that such fields are commonly using to guide magnetic particles (Qiu & Nelson Reference Qiu and Nelson2015) and also that an inspection of (5.2) reveals that the magnetic force on the helix vanishes due to the spatial uniformity of the magnetic flux density.

5.2 Computational results

The effects of linear viscoelasticity can be treated by considering the torque balance equation. Similarly to the wall interaction case (4.4), we have

where we consider the linear Maxwell model, ${\mathcal{L}}=1+\unicode[STIX]{x1D706}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t$ , and the Deborah number is given by $De=\unicode[STIX]{x1D706}/T$ .

For the shape parameters, we set $A/L=0.05$ , $r/L=0.003$ and $kL/2\unicode[STIX]{x03C0}=10$ , which gives the pitch of the helix to be $0.2L$ . We proceed by non-dimensionalising the length, time and force scales such that in the non-dimensional model the swimmer length, time period and viscosity are given by $L=1$ , $T=1$ and $\unicode[STIX]{x1D707}=1$ , which results in a non-dimensionalised constant representing the strength of the magnetic torque, given by $F_{mag}=C_{0}M_{0}B_{0}(\unicode[STIX]{x1D707}L^{3}/T)$ . We also set the initial condition of the magnetic helix so that the body-fixed frame and the inertial reference frame coincide, and we denote the origin of the body-fixed frame relative to the inertial frame by $\boldsymbol{X}(t)$ . With these initial conditions, the swimmer experiences no torque at the initial time $t=0$ . The resulting boundary integral equations are analogous to those in figure 3 and are computed with $N=366$ mesh elements and the time discretisation $\unicode[STIX]{x0394}t/T=1/500$ .

Figure 5.

The dynamics of permanently magnetised swimmers. (a) The movement of the helical magnetic swimmer within a Newtonian fluid for a rotating magnetic flux density and varying non-dimensional magnetic torque, $F_{mag}$ . The inset indicates the time evolution of the orientation of the swimmer magnetisation vector, as measured by $\unicode[STIX]{x1D703}_{mag}$ , the angle between $\boldsymbol{e}_{x}$ and $\hat{\boldsymbol{e}}_{x}$ , measured from the $x$ axis. (b) Movement of the helical magnetic swimmer with a non-dimensional magnetic torque $F_{mag}=0.03$ in fluids with different Deborah numbers. The inset indicates the corresponding time evolution of the orientation of the magnetisation vector as in (a). The red curves in (a) and (b) are therefore the same.

We first consider the motion in a Newtonian fluid, and the time evolution of the helix is plotted in figure 5(a) with varying non-dimensional magnetic torque, $F_{mag}$ . The plots shows the $z$ component of the swimmer, $\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{e}_{z}$ , as a function of time, where $\boldsymbol{e}_{z}$ denotes the unit inertial-frame basis vector in the $z$ direction. In the inset, the orientation of the magnetisation vector is shown via a plot of $\unicode[STIX]{x1D703}_{mag}$ , the angle between $\boldsymbol{e}_{x}$ and $\hat{\boldsymbol{e}}_{x}$ , as measured from the $x$ axis. For a sufficiently strong magnetic torque, the helix rotation is slave to that of the external field (figure 5 a inset), and the helix moves along the $z$ axis at constant speed. However, for weaker magnetic torques, there is a phase difference between the magnetisation vector and the external magnetic flux density. Once this phase difference exceeds  $\unicode[STIX]{x03C0}$ , the magnetisation vector then counter-rotates to relax to the external magnetic flux density. In turn, the helix also moves in the opposite direction, resulting in trajectories with oscillations (figure 5 a).

We then proceed to consider non-zero Deborah number, with the same initial conditions, and focus on swimming with $F_{mag}=0.03$ , which is insufficient to induce a rotation of the helix in a constant direction within a Newtonian fluid. Analogous plots of the position and the orientation of the magnetic swimmer are shown in figure 5(b) with increasing Deborah numbers. When $De=0.1$ , the phase difference between the magnetisation vector and the external magnetic vector is slightly suppressed, but nonetheless is still large enough to induce counter-rotation. Further increases of the Deborah number, however, amplify the effective torque that the swimmer experiences, and result in the magnetisation vector following the external magnetic field, and thus the swimmer rotates uniformly in time, with a dramatic increase in swimming propagation.