Stokes flow singularities in a two-dimensional channel: a novel transform approach with application to microswimming

A transform method for determining the flow generated by the singularities of Stokes flow in a two-dimensional channel is presented. The analysis is based on a general approach to biharmonic boundary value problems in a simply connected polygon formulated by Crowdy & Fokas in this journal. The method differs from a traditional Fourier transform approach in entailing a simultaneous spectral analysis in the independent variables both along and across the channel. As an example application, we find the evolution equations for a circular treadmilling microswimmer in the channel correct to third order in the swimmer radius. Significantly, the new transform method is extendible to the analysis of Stokes flows in more complicated polygonal microchannel geometries.


Introduction
The study of low Reynolds number flows in confined geometries bounded by no-slip walls is of perennial interest in fluid dynamics [1]. It has received a resurgence of attention in recent years owing to new applications in emerging technologies in micro-and nanofluidics, biological fluid dynamics and colloid science where the small scales render the role of inertia negligible. The relevant boundary value problems are linear but finding solutions even for the fundamental singularities of Stokes flows in just two dimensions can be challenging if the geometries are complicated. The monograph by Pozrikidis [1] documents solutions for the fundamental Stokes singularities in various geometries and surveys the variety of mathematical techniques available for their analysis.
This paper focuses on Stokes flow in the relatively simple geometry of a two-dimensional channel, but it offers a new mathematical perspective. The geometry is an important one in applications. Davis [2] has presented the solutions to various problems involving distributions of point forces, or stokeslets, in a two-dimensional channel with a view to understanding the blocking properties of periodic arrays of wall-attached barriers. His analysis makes use of Fourier transforms with respect to the variable along the channel length. Liron & Mochon [3] have found the Stokes flow field due to a three-dimensional force singularity at an arbitrary location and arbitrary orientation between two parallel plates. Section 3.5 of the monograph by Pozrikidis [1] describes a standard Fourier transform analysis of a two-dimensional point force in a channel.
Although it is not the principal focus of this article, the present study was motivated by problems of computing the swimming dynamics of low Reynolds number organisms in microchannels. In a recent paper, the authors [4] have given a general prescription for computing the dynamical system governing the motion of a model two-dimensional circular treadmilling swimmer in an arbitrary domain confined by no-slip walls. By an asymptotic expansion, the governing ordinary differential equations are determined correct to third order in the swimmer radius: for a given fluid region, rather than solving a boundary value problem for the full finitearea swimmer, one solves boundary value problems for certain point singularities of Stokes flow-namely, a torque-free stresslet and an irrotational (source) quadrupole-in that domain. Even though this is usually simpler than solving the full boundary value problem for the finitearea swimmer, it remains mathematically challenging for most geometries and this challenge has led to the results herein.
Our principal aim here is to present a transform approach to finding analytical solutions for those singular Stokes flows in the channel geometry. We adapt a new general transform approach to biharmonic boundary value problems in a simply connected polygon formulated by Crowdy & Fokas [5], who originally showcased their method for canonical semi-strip problems in plane elastostatics. The method differs from the usual Fourier transform approach where transforms are taken only with respect to the independent variable along the channel length [1,2]; instead, a simultaneous spectral analysis with respect to the two independent variables both along and across the channel is performed. It offers various mathematical advantages, not least of which is the fact that the approach can be generalized to any simply connected polygonal geometry. This work appears to be the first to apply the general method of [5] to problems of slow viscous Stokes flows, where biharmonic boundary value problems also arise. It showcases the method in one of the simplest non-trivial situations (the channel geometry) and therefore constitutes a valuable case study. The simplicity of the geometry also allows an instructive comparison with traditional Fourier transform methods detailed in an appendix. As an illustrative example, in §5, we apply the method to our motivating problem of treadmilling microswimmers in a channel. The streamfunction then satisfies the biharmonic equation where ∇ 2 is the two-dimensional Laplacian. Crowdy & Fokas [5] present their approach to biharmonic boundary value problems in a polygon in terms of ψ and its various partial derivatives, but there is an alternative equivalent presentation of the same method in terms of the so-called Goursat functions [6]. With z = x + iy the Goursat functions, here denoted by f (z) and g(z), are analytic functions of z in the fluid region with the property that the solution of (2.2) is represented in the form

Stokes flows in two dimensions
The Goursat functions are familiar to fluid dynamicists using complex variable methods for Stokes flows [6] so we adopt this formulation and seek integral representations of the required physical quantities u, v, fluid pressure p and fluid vorticity ω by making use of the following relations of the latter quantities to the Goursat functions [6]:

Point stresslet in a channel
One important singularity of Stokes flow is the stresslet [1]. Consider a channel occupied by fluid at zero Reynolds number in the region −∞ < x < ∞, 0 < y < h. A torque-free stresslet of strength Ṽ e iα at a point z 0 requires that f (z) and g (z) have the local expansions [7] f (z) = Ṽ e 2iα z − z 0 + locally analytic function and g (z) = Ṽ e 2iα z 0 (z − z 0 ) 2 + locally analytic function.
Our objective is to find the flow induced in the channel by such a singularity. The velocity as |x| → ∞ is taken to vanish. We focus on this stresslet for two reasons: first, it is relevant to the study of swimming organisms considered later; second, we have been unable to find the solution reported elsewhere (although it can in principle be derived as an appropriate parametric derivative, with respect to the singularity position, of the solution for a point stokeslet in a channel given by Davis [2]). It should be pointed out that the quantity Ṽ e iα is just a complex constant but its designation has been chosen so as to be consistent with our later application of the results to a model treadmilling swimmer. It is convenient to introduce the meromorphic functions The functions f s (z) and g s (z) have precisely the required local form (3.1) at z 0 and no other singularities in the channel. Both functions are periodic with period 2 hi, i.e. 3) The derivative f s (z) is also 2 hi-periodic and this will be useful later. We now pose the ansatz where f R (z), g R (z) are analytic in the channel and vanish as |x| → ∞; note that the second term in the definition of g s (z) has been chosen to ensure that It should be emphasized that it is not necessary to introduce the hyperbolic functions (3.2) and one can proceed instead by defining f s (z) and g s (z) to have the simple rational function forms in (3.1). However, we have found that the periodicity properties (3.3) greatly simplify the algebraic manipulations to follow. The no-slip conditions on the two channel walls are On substitution of the ansatz (3.4), we find that, onz = z, and, onz = z − 2 hi, where known functions appear on the right-hand sides of (3.7) and (3.8). The Schwarz conjugate of an analytic function h(z) is defined by Following the general formulation in [5] the two functions f R (z) and g R (z), which are analytic in the channel and decaying in the far-field, will be represented by where ρ 1 (k) and ρ 2 (k) are defined by and, similarly, where ρ 3 (k) and ρ 4 (k) are defined by The set {ρ n (k)|n = 1, 2, 3, 4} is called the set of spectral functions [5,8], and there are relations between these spectral functions referred to as global relations. In this case, there are two global relations given by (3.14) In this case, these relations can be viewed as a consequence of Cauchy's theorem: for example, the condition ρ 2 (k) = −ρ 1 (k) is equivalent to the statement that the integral of f R (z) e −ikz around the channel is zero because, for k ∈ R, the latter function is analytic inside the channel.  (3.15) or, on making the substitution k → −k in the second integral, Similarly, we can write It is worth emphasizing that it is this spectral representation of the analytic functions f R (z) and g R (z) that can be extended to arbitrary simply connected polygons [5,8,9].

(a) Preliminary observations
The representations of f R (z) and g R (z) above involve only real values of k. Hence, for real k, complex conjugation implies The integral on the left-hand side will appear in the transform of the boundary conditions on the wall at y = 0. Similarly, on writing x = z − ih in the integral expression for ρ 2 (k), we can write On taking a complex conjugate for real k, along the real z-axis: The function R 1 (k) is known. On use of the result, available from an integration by parts, Similarly, if we multiply the second boundary condition (3.8) by e −ikz and integrate over the upper boundary, we find Equations (3.27) and (3.28) can be analysed, together with the global relations (3.14), to determine the spectral functions. Indeed, on use of (3.14) equation (3.28) becomes On subtraction of (3.27) from (3.30), we find After taking a complex conjugate and letting k → −k, we find Elimination of ρ 1 (−k) from (3.31) and (3.33) produces The function R(k) can be computed explicitly using residue calculus to give Given ρ 1 (k), the spectral function ρ 3 (k) follows from (3.27) together with which also follows from a residue calculation. All physical quantities are now available from the representations, dependent only on f s (z), g s (z) and the associated spectral functions ρ 1 (k) and R 1 (k), given by This is the required solution. Significantly, these integrals are readily computed numerically owing to the exponential decay of all integrands as k → ∞. As k → 0, Note that b is real-valued and c becomes purely imaginary if y 0 = h/2 or sin2α = 0. Expressions (3.40) can be used to verify that all integrands in (3.38) are regular at k = 0. Formula (3.39) is also useful in facilitating explicit computation of the limiting behaviour of the integrands in any numerical evaluation procedure. Finally, it should be pointed out that although the integrals defining all physical quantities in (3.38) are regular at k = 0 the integrand in the expression (3.16) for f R (z) is singular at k = 0 but this is rendered immaterial by a freedom to redefine f (z) → f (z) + pz + q and g (z) → g (z) + q for constants p ∈ R, q ∈ C without affecting any of the physical quantities. It is instructive, for purposes of comparison, to solve the problem of a point stresslet in a channel using a more traditional approach involving a Fourier transform only with respect to x. This separate analysis in appendix A also affords us a numerical check on (3.38).

Point quadrupole in a channel
It is straightforward to adapt the approach to other singularities. Indeed, the solution for a stokeslet singularity [1] in a channel can be derived by simple adaptation of the above approach and compared with the result of Davis [2], who employed standard Fourier transform methods akin to those used in appendix A for the stresslet problem. For use in §6, consider the separate problem in which an irrotational quadrupole is situated at z 0 with the corresponding function f (z) analytic at z 0 and g (z) behaving like g (z) = 2 Ṽ e 2iα (z − z 0 ) 3 + locally analytic function. use the same notation as in §3). But it turns out that we can bypass a full analysis by noticing that the functions f (z) and g (z) relevant for the point quadrupole problem can be derived from those of §3 by computing the mixed parametric derivatives This generation of the flow due to a quadrupole by taking parametric derivatives with respect to the singularity position is equivalent to use of the ∇ 2 0 operator by Davis & Crowdy [4] in their analysis of low Reynolds number swimmers in confined domains (the application to be considered in §5). It follows that, for the point quadrupole problem, with spectral functions in the integral representations for the corresponding f R (z) and g R (z) now given by and

Application: microswimmers in a channel
One application of the foregoing analysis is to modelling the low Reynolds number swimming of micro-organisms in a channel. In 1963, Lord Rothschild [10] conducted an experiment in which a sample of bull sperm in a very viscous fluid was placed between two glass plates; after time he measured the distribution of the sperm between the two plates and found a tendency for them to aggregate at the channel walls. The experiment was repeated recently by Berke et al. [11] and a hydrodynamic explanation of the wall attraction phenomenon proposed. But more general dynamical behaviour of low Reynolds number swimmers near walls has been observed, including steady translation and periodic motions along the walls [12,13]. To rationalize this behaviour using simple mathematical models, Crowdy & Or [7] have proposed the study of a simple circular 'treadmiller' comprising a cylindrical circular swimmer with an imposed tangential velocity profile; this profile actuates motion in the spirit of an 'envelope model' of surface ciliatory motion [14]. Crowdy & Or [7] showed that, in free space, such a swimmer has the singularity distribution comprising a torque-free stresslet and a superposed irrotational quadrupole; as an approximation, they studied the motion of such a singularity combination near a straight no-slip wall and found good qualitative agreement with numerical [12] and laboratory [15] experiments. Davis & Crowdy [4] have since found the general evolution equations for the same swimmer correct to O( 3 ), where is the swimmer radius. We now use the latter equations to study the dynamics of such a treadmilling swimmer in a channel. The swimmer, which is oriented at some angle α(t) to the real axis, has an imposed tangential surface velocity of the functional form 2 sin 2(φ − α(t)), where φ is an angular coordinate around the swimmer boundary. The interaction of the Stokes flow it generates with any nearby solid walls generally results in dynamical evolution of its centre z 0 (t). A schematic of such a swimmer located in a channel of width h is shown in figure 1.
The evolution equations given in Davis & Crowdy [4] have the following form. Suppose ψ R denotes the streamfunction of the reflected velocity generated by the wall interaction of the stresslet considered in §3, where  Figure 1. Circular treadmilling swimmer, of the type proposed in [7], centred at z 0 (t) at orientation α(t) in a channel of width h. The radius of the swimmer is 1. Its evolution equations, correct to O( 3 ), can be determined from the results in [4] using the transform results of § §3 and 4. (Online version in colour.) Suppose too that the streamfunctionψ R for the reflected velocity generated by the point quadrupole of §4 can similarly be determined in the form The evolution of the swimmer, correct to O( 3 ), is given [4] by the solution of the system and The main point here is that the analysis of § §3 and 4 allow all the quantities (5.6) appearing in the approximate evolution equations (5.5) to be calculated and, thus, the swimmer motion determined. Figure 2 shows some typical swimmer trajectories for α(0) = π/4 and x(0) = 0, y(0) = 0.1(0.1)0.9 for h = 1. Consistent with the experimental observations of wall aggregation [10,11], it is found that the swimmer is attracted to the nearest wall of the channel with only swimmers exactly positioned on the channel centre-line remaining stationary. A range of more general initial orientations were explored and all swimmers off the centre-line are similarly found to be attracted to the nearest wall. We found no initial conditions for which the swimmer propagates along the channel. It is important to recall that these evolution equations are approximate and only correct to O( 3 ); it would be of interest, for purposes of comparison, to compute the full unapproximated evolution equations but this would probably require a numerical treatment using, for example, boundary integral methods [1].

Discussion
Although our motivating application was the dynamics of low Reynolds number swimming organisms, the principal contribution of this paper is mathematical: a novel transform approach to finding the flows due to basic point singularities of the Stokes equations in a simple twodimensional channel has been given. Other types of Stokes flow in the channel can be studied using adaptations of the same approach. The solutions can be derived using alternative approaches, such as traditional Fourier methods, but the method here, based on a general formulation by Crowdy & Fokas [5], is significant in that it can be generalized to finding flows in more complicated simply connected polygonal domains (a similar solution scheme for the Laplace equation in a polygon has been given by Fokas & Kapaev [9]). In those cases, the determination of f (z) and g (z) will follow similar steps: a convenient form of the singular terms f s (z) and g s (z) must be chosen based on considerations from the geometry being analysed and the integral representations for the correction terms f R (z) and g R (z) will be similar except that every additional side of the polygonal domain will introduce a new spectral function. The set of spectral functions will again satisfy a global relation which must be analysed to determine any unknown spectral functions [5,8].
Here, the global relations are manifested in equations (3.14). In many cases, the analysis of these involves the solution of a Riemann-Hilbert problem for the spectral functions [5,8] but the problem in a channel geometry studied here does not require any such Riemann-Hilbert analysis and the spectral functions can be determined using only algebraic manipulations. Other boundary value problems for which a Riemann-Hilbert analysis can be avoided, and only algebraic manipulations used, have been reported by Fokas [8]. The latter monograph also provides an interesting survey of how the transform method for Stokes flow used here fits into a unified framework for boundary value problems developed in recent years by Fokas and collaborators.
Extensions of our approach to finding Stokes flow solutions, and studying microswimmer dynamics, in more complicated polygonal channel geometries are currently under investigation.