Axisymmetric Stokes flow due to a point-force singularity acting between two coaxially positioned rigid no-slip disks

We investigate theoretically on the basis of the steady Stokes equations for a viscous incompressible fluid the flow induced by a Stokeslet located on the centre axis of two coaxially positioned rigid disks. The Stokeslet is directed along the centre axis. No-slip boundary conditions are assumed to hold at the surfaces of the disks. We perform the calculation of the associated Green's function in large parts analytically, reducing the spatial evaluation of the flow field to one-dimensional integrations amenable to numerical treatment. To this end, we formulate the solution of the hydrodynamic problem for the viscous flow surrounding the two disks as a mixed-boundary-value problem, which we then reduce into a system of four dual integral equations. We show the existence of viscous toroidal eddies arising in the fluid domain bounded by the two disks, manifested in the plane containing the centre axis through adjacent counterrotating eddies. Additionally, we probe the effect of the confining disks on the slow dynamics of a point-like particle by evaluating the hydrodynamic mobility function associated with axial motion. Thereupon, we assess the appropriateness of the commonly-employed superposition approximation and discuss its validity and applicability as a function of the geometrical properties of the system. Additionally, we complement our semi-analytical approach by finite-element computer simulations, which reveals a good agreement. Our results may find applications in guiding the design of microparticle-based sensing devices and electrokinetic transport in small scale capacitors.


Introduction
Manipulating colloidal particles suspended in viscous media is a challenging task and is of paramount importance in various fields of engineering and natural sciences. Frequently, taking into account the fluid-mediated hydrodynamic interactions between particles moving through a liquid is essential to predict the behaviour of colloidal suspensions and polymer solutions (Probstein 2005;Mewis & Wagner 2012). Recent advances in micro-and nanofluidic technologies have permitted the fabrication and manufacturing of channels with well-defined geometries and characteristic dimensions ranging from the micro-to the nanoscale. A deep understanding of the nature of the mutual interactions between particles and their confining interfaces is of crucial importance in guiding the design of devices and tools for an optimal nanoscale control of biological macromolecules. Notable examples include single-molecule manipulation (Turner et al. 1998;Campbell et al. 2004), DNA mapping for genomic applications Reisner et al. 2005;Persson & Tegenfeldt 2010), DNA separation and sorting (Doyle et al. 2002;Cross et al. 2007;Xia et al. 2012), and rheological probing of complex structures using Atomic Force Microscopy cantilevers (François et al. 2008(François et al. , 2009Dufour et al. 2012;Darwiche et al. 2013).
At these small scales, fluid flows are governed by low-Reynolds-number hydrodynamics, where viscous effects dominate over inertial effects (Kim & Karrila 2013). Solutions for fluid flows due to point forces, or Stokeslets, acting close to confining boundaries have been tabulated for various types of geometries, as summarised in the classic textbook by Happel & Brenner (1983). The study of the fluid-mediated hydrodynamic interactions in a channel confinement has received significant attention from many researchers over the past couple of years. In the following, we provide a survey of the current state of the art and summarise the relevant literature in this subject.
The first attempt to address the motion of a spherical particle confined between two infinitely extended no-slip walls dates back to Faxén (1921), who calculated in his PhD dissertation the hydrodynamic mobility parallel to the walls. These calculations have been performed when the particle is located in the quarter-plane or mid-plane between the two confining walls (Happel & Brenner 1983). Later, Oseen (1928) suggested that the hydrodynamic mobility between two walls could approximately be obtained by superposition of the contributions resulting from each single wall. A modified coherent superposition approximation has further been suggested by Benesch et al. (2003), providing the diffusion coefficients of a Brownian sphere in confining channels. These predictions were found to match more accurately the existing experimental data reported in the literature.
Exact solutions for a point-force singularity acting at an arbitrary position between two walls have first been obtained using the image technique in a seminal article by Liron & Mochon (1976). It has been noted that the effect of the second wall becomes important when the distance separating the particle from the closest wall is larger than approximately one tenth of the channel width (Brenner 1999). Using this solution, Liron (1978) further investigated the fluid transport problem of cilia between two parallel plates. A joint analytical-numerical approach (Ganatos et al. 1980a,b) as well as a multipole expansion technique (Swan & Brady 2010) were presented to address the motion of an extended particle confined between two hard walls. Meanwhile, Bhattacharya & Bławzdziewicz (2002) constructed the image system for the flow field produced by a force multipole in a space bounded by two parallel walls using the image representation for Stokes flow. In addition, compressibility effects were examined by Felderhof (2006Felderhof ( , 2010a. In this context, Hackborn (1990) investigated the asymmetric Stokes flow between two parallel planes due to a rotlet singularity, the axis of which is parallel to the boundary planes. Further, Ozarkar & Sangani (2008) prescribed an analytical approach using the image-system technique for determining the Stokes flow around particles in a thin film bounded by a wall and a gas-liquid interface. More recently, Daddi-  have provided the frequency-dependent hydrodynamic mobility functions between two planar elastic interfaces endowed with resistance toward shear and bending deformation modes.
Experimentally, Dufresne et al. (2001) reported direct imaging measurements of a colloidal particle diffusing between two parallel surfaces, finding a good agreement with the superposition approximation suggested by Oseen. In addition, video microscopy (Faucheux & Libchaber 1994) combined with optical tweezers (Lin et al. 2000;Tränkle et al. 2016) as well as dynamic light scattering (Lobry & Ostrowsky 1996) have also allowed for good agreement with available theoretical predictions. Further experimental investigations have focused on DNA conformation and diffusion in slit-like confinements (Stein et al. 2006;Balducci et al. 2006;Strychalski et al. 2008;Tang et al. 2010;Graham 2011;Dai et al. 2013;Jones et al. 2013).
Concerning collective properties, the behaviour of suspensions in a channel bounded by two planar walls has received a lot of attention. For instance, Bhattacharya et al. (2005) examined the fluid-mediated hydrodynamic interactions in a suspension of spherical particles confined between two parallel planar walls under creeping-flow conditions. In addition, Bhattacharya (2008) considered the collective motion of a two-dimensional periodic array of colloidal particles in a slit pore. Using a novel accelerated Stokesiandynamics algorithm, Baron et al. (2008) performed fully-resolved computer simulations to investigate the collective motion of linear trains and regular square arrays of particles suspended in a viscous fluid bounded by two parallel plates. Further,  analysed the far-field response to external forcing of a suspension of particles in a channel. Meanwhile, Swan & Brady (2011) presented a numerical method for computing the hydrodynamic forces exerted on particles in a suspension confined between two parallel walls. Furthermore, Saintillan et al. (2006) employed Brownian dynamics simulations to investigate the effect of chain flexibility on the cross-streamline migration of short polymers in a pressure-driven flow between two flat plates. The latter numerical study confirmed the existence of a shear-induced migration toward the channel centreline away from the confining solid boundaries.
The hydrodynamic problem of particles freely moving between plane-parallel walls in the presence of an incident flow has further been considered in still more details. Under an external flow, Uspal et al. (2013) showed how shape and geometric confinement of rigid microparticles can conveniently be tailored for self-steering. Jones (2004) made use of a two-dimensional Fourier-transform technique to obtain an analytic expression of the Green tensor for the Stokes equations with an incident Poiseuille flow. In addition, he provided the elements of the resistance and mobility tensors in this slit-like geometry. Besides, Bhattacharya et al. (2006) introduced a novel numerical algorithm based on transformations between Cartesian and spherical representations of Stokes flow to account for an incident Poiseuille flow. Meanwhile, Staben et al. (2003) presented a novel boundary-integral algorithm for the motion of a particle between two parallel planar walls in Poiseuille flow. The boundary-integral method formulated in their work allowed to directly incorporate the effects of the confining walls into the stress tensor, without requiring discretisation of the two walls. In this context, Griggs et al. (2007) and Janssen & Anderson (2007, 2008 employed boundary-integral methods to examine the motion of a deformable drop between two parallel walls in Poiseuille flow, where lateral migration towards the channel centre is observed. Geometric confinements significantly alter the behaviour of swimming microorganisms and can affect the motility of self-propelling active particles in a pronounced way (Lauga & Powers 2009;Menzel 2013Menzel , 2015Lauga 2016;Zöttl & Stark 2016;Bechinger et al. 2016;Ostapenko et al. 2018;Gompper et al. 2020;Shaebani et al. 2020). Surface-related effects on microswimmers can lead to crucial implications for biofilm formation and microbial activity. In a channel bounded by two walls, Bilbao et al. (2013) studied the locomotion of a model nematode, finding that the swimming organism tends to swim faster and navigate more effectively under confinement. Furthermore, Wu et al. (2015Wu et al. ( , 2016 investigated the effect of confinement on the swimming behaviour of a model eukaryotic cell undergoing amoeboid motion. There, the swimmer has been modeled as an inextensible membrane deploying local active force. It has been found that confinement can strongly alter the swimming gait. In addition, Brotto et al. (2013) described theoretically the dynamics of self-propelling active particles in rigidly confined thin liquid films. They demonstrated that, due to hydrodynamic friction with the nearby rigid walls, confined microswimmers do not only reorient themselves in response to flow gradients but they can also show reorientation in uniform flows. In this context, Mathijssen et al. (2016) investigated theoretically the hydrodynamics of self-propelling microswimmers in a thin film. Besides, Daddi-Moussa-Ider et al. (2018) examined the behaviour of a three-sphere microswimmer in a channel bounded by two walls, where different swimming states have been observed. More recently, amoeboid swimming in a compliant channel was numerically investigated (Dalal et al. 2020).
In all of the above-mentioned studies, the confining channel was assumed to be of infinite extent or periodically replicated along the lateral directions. Instead, we here consider the hydrodynamic problem for a point force acting near two coaxially positioned disks of finite radius. In many biologically and industrially relevant applications, finitesize effects become crucial for an accurate and reliable description of transport processes  (Daiguji et al. 2004), and the rheology of droplets, capsules, or cells in constricted/structured microchannels (Park & Dimitrakopoulos 2013;Le Goff et al. 2017;Trégouët et al. 2018Trégouët et al. , 2019, where boundary effects may play a pivotal role. In this contribution, we take a step toward addressing this context by presenting an analytical theory for the viscous flow resulting from a Stokeslet singularity acting along the centre axis of two coaxially positioned disks of no-slip surfaces. We formulate the hydrodynamic problem as a mixed-boundary-value problem, which we then transform into a system of dual integral equations. Along this path, we show that the solution of the flow field in the fluid region bounded by the two disks exhibits viscous toroidal eddies. In addition to that, we derive expressions for the hydrodynamics mobility functions and discuss the applicability and limitations of the superposition approximation. Moreover, we support our semi-analytical results by numerical simulations using a finite-element method (FEM), which leads to a good agreement.
The remainder of this article is organised as follows. In Sec. 2, we formulate the problem mathematically and derive the corresponding system of dual integral equations, from which the solution for the hydrodynamic flow fields can be obtained. We then make use of this solution in Sec. 3 to yield an integral expression of the mobility function of a point-like particle slowly translating along the axis of the disks. Concluding remarks and outlooks are contained in Sec. 4. In Appendix A, we detail the analytical derivation of the kernel functions arising in the resulting integral equations.

Mathematical formulation
We examine the axisymmetric flow induced by a Stokeslet singularity acting on the axis of two coaxially positioned circular disks of equal radius R. Moreover, we suppose that the disks are located within the planes z = −H/2 and z = H/2 with H denoting the separation distance between the disks. Their centres are positioned on the z axis. In addition, we assume that the surrounding viscous fluid is Newtonian, of constant dynamic viscosity η, and that the flow is incompressible.

Governing equations
In low-Reynolds-number hydrodynamics, the fluid dynamics is governed by the Stokes equations (Happel & Brenner 1983) where v and σ denote, respectively, the fluid velocity field and the hydrodynamic stress tensor. For a Newtonian fluid, the latter is given by σ = −pI + 2ηE, where p is the pressure field and E = 1 2 ∇v + (∇v) T is the rate-of-strain tensor, with the superscript T denoting a transpose. In addition, δ stands for the Dirac delta function, and F is the amplitude of a stationary point force acting on the fluid at position r 0 = hê z , where −H/2 < h < H/2, withê z denoting the unit vector along the z direction. See Fig. 1 for an illustration of the system setup. In the remainder of this article, we scale all the lengths involved in the problem by the separation H of the two disks.
We designate by the subscript 1 the variables and parameters in the fluid region underneath the plane containing the lower disk, for which z −1/2, by the subscript 2 the fluid domain bounded by the planes z = −1/2 and z = 1/2, and by the subscript 3 the region above the plane containing the upper disk, for which z 1/2. Since the system is axisymmetric, all field variables are thus functions of the radial and axial coordinates only. Accordingly, the Stokes equations (2.1) can be projected onto the cylindrical coordinate system as wherein v r and v z denote the radial and axial fluid velocities, respectively, and ∆ is the Laplace operator given by We note that the three-dimensional Dirac delta function is expressed in axisymmetric cylindrical coordinates as δ(r − r 0 ) = (πr) −1 δ(r)δ(z − h) (Bracewell 1999).
In an unbounded viscous fluid, i.e., in the absence of the disks, the solution of Eqs. (2.2) is given by the Oseen tensor, commonly denominated as the free-space Green function (Kim & Karrila 2013) with the distance from the position of the point force ρ = r 2 + (z − h) . The corresponding pressure field reads In the presence of the confining disks, the solution of the flow problem can be expressed as a superposition of the solution in an unbounded fluid, given above by Eqs. (2.4) and (2.5), and a complementary solution, the sum of the two solutions being required to satisfy the underlying regularity and boundary conditions. Then wherein v * and p * stand for the complementary solutions (also referred to as the image solution (Blake 1971)) for the velocity and pressure fields, respectively. For an axisymmetric Stokes flow, the general solution can be expressed in terms of two harmonic functions φ and ψ as (Imai 1973;Kim 1983 In each of the three fluid domains introduced above, the solution of Laplace's Eqs. (2.8) can be expressed in terms of Fourier-Bessel integrals as for i ∈ {1, 2, 3}, with λ denoting the wavenumber and J k the kth-order Bessel function of the first kind (Abramowitz & Stegun 1972). In addition, A ± i and B ± i are wavenumberdependent unknown coefficients, to be determined from the regularity and boundary conditions. Then, the components of the image velocity and pressure fields are given by

Boundary conditions and dual integral equations
As regularity conditions, we require for the image field vanishing velocity and pressure far away from the singularity location as ρ → ∞. This implies that In what follows, to simplify notations, we drop the plus sign in the fluid domain underneath the lower disk to denote A 1 = A + 1 and B 1 = B + 1 , and we drop the minus sign in the fluid domain above the upper disk to denote A 3 = A − 3 and B 3 = B − 3 . The boundary conditions consist of requiring (a) the natural continuity of the total fluid velocity field at the interfaces between the fluid domains, (b) vanishing total velocities at the surfaces of the disks (the no-slip and no-permeability boundary condition (Lauga et al. 2007)), and (c) continuity of the total viscous-stress vectors at the interfaces between the fluid domains outside the regions occupied by the disks. Mathematically, these conditions can be expressed as where the components of the stress vector are expressed in cylindrical coordinates for an axisymmetric flow field by Applying the continuity of the radial components of the fluid velocity at the surfaces occupied by the two disks yields the expressions of the wavenumber-dependent coefficients associated with the intermediate fluid domain bounded by the two disks as functions of those in the lower and upper fluid domains. Defining where the matrix Q is given by (2.14) Here, we have defined for convenience the abbreviations s = sinh(λ) and c = cosh(λ). In addition, φ ± = λ (λ ± 1) + se −λ . On the one hand, by addressing the no-slip velocity boundary conditions at the surfaces of the disks prescribed by Eqs. (2.11b) and projecting the resulting equations onto the radial and tangential directions, four integral equations on the inner domain are obtained, wherein the terms appearing on the right-hand sides in these equations are radial functions resulting from the evaluation of the terms associated with the flow velocity field induced by the free-space Stokeslet at the surfaces of the coaxially positioned disks. They are explicitly given by On the other hand, four integral equations on the outer domain are obtained by addressing the continuity of the hydrodynamic stress vector at z = ±1/2 prescribed by Eq. (2.11c). They can be cast in the form where we have defined the wavenumber-dependent quantities (2.13) and (2.14), Eqs. (2.15) through (2.18) form a system of four dual integral equations (Tricomi 1985) for the unknown wavenumber-dependent coefficients regrouped in X 13 . A solution of such types of dual integral equations with Bessel kernels can be obtained by the methods prescribed by Sneddon (1960Sneddon ( , 1966 and Copson (1961). A similar procedure has recently been employed by some of us to address the axisymmetric flow induced by a Stokeslet near a circular elastic membrane (Daddi-Moussa-Ider et al. 2019), and the asymmetric flow field near a finite-sized rigid disk (Daddi-Moussa-Ider et al. 2020). Once X 13 is determined from solving the dual integral equations derived above, the remaining wavenumber-dependent coefficients expressed by X 2 follow forthwith from Eqs. (2.13) and (2.14).
The core idea of our solution approach consists of expressing the solution of Eqs. (2.17) as definite integrals of the forms
(2.21d ) Next, by substituting Eqs. (2.19) into Eqs. (2.21) and interchanging the order of the integrations with respect to the variables t and λ, the equations associated with the inner problem can be expressed in the following final forms where the kernels L i : [0, R] 2 → R, for i ∈ {1, 2, 3, 4} are complex mathematical functions that are defined and provided in Appendix A. Equations (2.22) form a system of four Fredholm integral equations of the first kind (Smithies 1958;Polyanin & Manzhirov 1998) for the unknown functions f i (t), i ∈ {1, 2, 3, 4}. Due to the complicated nature of the kernel functions, we recourse to numerical solutions.

Numerical solution of the integral equations and comparison with FEM simulations
We now summarise the main steps involved in the numerical computations of the flow field. First, the integration over the intervals [0, R] in Eqs. (2.22) are partitioned into N subintervals and each integral is approximated by the standard middle Riemann sum (Davis & Rabinowitz 2007). The four resulting equations are evaluated at N values of t j that are uniformly distributed over the interval [0, R] such that t j = (j −1/2)(R/N ), with j = 1, . . . , N . Secondly, the discrete values of f i (t j ), with i ∈ {1, 2, 3, 4} are obtained by solving the resulting linear system of 4N equations. Thirdly, the four integrals in Eqs. (2.19) are converted into well-behaved definite integrals over [0, π/2] by using the change of variable λ = tan u and thus dλ = du/ cos 2 u. Thereupon, the resulting integrals are also approximated by the middle Riemann sum, and the wavenumber-dependent functions g i (λ k = tan u k ), k = 1, . . . , M , are evaluated at discrete values of u k such that u k = (k − 1/2)(π/2)/M . Fourthly, the values of X 2 at each discrete point λ k are readily obtained by inverting the linear system of four equations given by Eqs. (2.18). In addition, it follows from Eq. (2.13) that X 13 = Q −1 · X 2 . Finally, the image flow fields are obtained from Eqs. (2.10) by approximating, again, the integrals by the middle Riemann sum.
Even though the approach employed here may seem cumbersome at a first glance, it has the advantage of being amenable to straightforward implementation. Unlike many direct numerical simulation techniques which generally require discretisation of the entire three-dimensional fluid domain, or of at least an effectively two-dimenional domain when the axial symmetry is exploited, the integral formulation presented in this work reduces the solution of the flow problem to a set of one-dimensional integrals. Besides, the present semi-analytical approach might serve as a motivation for various theoretical investiga-tions of related problems that could possibly pave the way towards real engineering applications.
In Fig. 2, we present a log-log plot the variations of the discretisation error (Roy 2010) associated with the numerical computation of the amplitude of the image velocity field versus the number of discrete points used in the numerical integration of Eqs. (2.22) while keeping M = 10N in the discretisation of Eqs. (2.19) and (2.10). The error is estimated relative to the numerical solution on a finer gird size for N = 15000 and M = 150000 at three different points of the fluid domain. We observe that the error decays approximately algebraically as N −3/2 over the whole range of considered values of N and lies well below 10 −3 % for N 5000. We have checked that a similar behaviour is also found when varying the position of the Stokeslet or the evaluation point within the fluid domain.
To validate our semi-analytical solution, we perform direct numerical simulations for the same geometry as well. We use a piecewise quadratic finite-element discretisation of the Stokes problem stated by Eqs. (2.2) in cylindrical coordinates. Since such an equalorder discretisation does not satisfy the inf-sup condition, we add stabilisation terms of local projection type (Becker & Braack 2001). The numerical domain is artificially limited to (0, R) × (−Z, Z) with R, Z ∈ R being sufficiently large numbers so as to avoid spurious feedback to the region of interest close to the plates. In addition, the Dirac delta function forcing the flow is represented exactly in the variational formulation by means of where φ z is the test function corresponding to the vertical direction. Numerically, the singularity calls for very fine mesh resolution close to r 0 and in proximity to the coaxially positioned plates, which we accomplish by local mesh adaptivity (Braack & Richter 2006). Further details on the discretisation method and the solution of the resulting linear systems of equations can be found in Richter (2017). In Fig. 3, we represent the graphs of the resulting streamlines as well as contour plots of the total velocity field resulting from a Stokeslet singularity axisymmetrically acting at various positions along the axis of two coaxially disposed disks of unit radius. Here, we set the numbers of discrete points to N = 15000 and M = 150000 in our numerical evaluation of the analytical description. The magnitude of the scaled velocity field is shown on a logarithmic scale in order to better appreciate the difference in magnitude between the different fluid regions. In each panel, we depict on the left-hand side the results obtained via our semi-analytical approach derived in the present work. On the right-hand side in each panel, we include the corresponding flow fields determined via the FEM simulations. Good agreement between the two solution procedures is obtained over the whole fluid domain, demonstrating the robustness and applicability of our semi-analytical approach. Most noticeably, we observe the existence of adjacent counterrotating eddies, the axis of rotation of which is directed along the azimuthal direction. Accordingly, the resulting flow field in the inner region consists of toroidal eddies on account of the axisymmetric nature of the flow (Moffatt 1964). In contrast to that, descending streamlines are obtained in the outer region. For infinitely large disks, analogous toroidal structures have previously been identified and proven to decay exponentially with distance from the singularity position (Liron & Blake 1981). Moreover, we remark that the overall magnitude of the flow field becomes less important as the point force gets closer to a confining plate. This behaviour is accompanied by a notable increase of the asymmetry of the counterrotating eddies.
Having derived the solution of the flow problem due to an axisymmetric Stokeslet acting near two finite-sized coaxially positioned disks, we next employ our formalism to recover the solution earlier obtained by Liron & Mochon (1976) for a Stokeslet acting between two parallel planar walls of infinite extent along the transverse direction.

Solution for R → ∞
For infinitely large disks, the integral equations (2.21) in the inner domain become defined for the whole axis of positive real numbers. Accordingly, the solution for the unknown functions g i (λ), for i ∈ {1, 2, 3, 4} can be obtained using inverse Hankel transforms. By making use of the orthogonality property of Bessel functions (Abramowitz & Stegun 1972 (2.24) we readily obtain where we have defined the unknown vector g = (g 1 , g 2 , g 3 , g 4 ) T , the wavenumberdependent matrix (2.26) and whereψ = ψ + 1 ,ψ − 1 ,ψ + 2 ,ψ − 2 T gathers the inverse Hankel transforms of the previously introduced auxiliary functions defined by Eqs. (2.16). Specifically,

Hydrodynamic mobility
Our calculation of the flow field presented in the previous section can be employed in order to probe the effect of the two hard disks on the hydrodynamic drag acting on an enclosed point-like particle axially moving along the coaxially positioned axis. This effect is commonly quantified by the hydrodynamic mobility function, which relates the velocity of a particle to the net force exerted on its surface (Leal 1980;Swan & Brady 2007;, 2017Driscoll & Delmotte 2019). In a bulk Newtonian fluid of constant dynamic viscosity η, the mobility function µ of a spherical particle of radius a is given by the familiar Stokes law, which states that in this case the mobility is µ 0 = 1/(6πηa) (Stokes 1851). In the presence of the confining disks, the leading-order correction to the particle mobility for an axisymmetric motion along the axis is obtained by evaluating the image flow field at the particle position as (3.1) Evaluating the limit in the latter equation and scaling by the bulk mobility, the scaled correction to the particle mobility is obtained as (1 − λh) A + 2 − λB + 2 e λh + (1 + λh) is a positive dimensionless number commonly denominated as the correction factor of the Stokes steady mobility (Happel & Brenner 1983). Unfortunately, an analytical evaluation of this infinite integral is not auspicious. Therefore, we recourse to a numerical evaluation. For infinitely large disks, i.e., as R → ∞, the correction factor k in Eq. (3.2) can conveniently be cast into the simple integral form where we have defined the wavenumber-dependent function This result is found to be in full agreement with the expression obtained by Swan & Brady (2010), who used a two-dimensional Fourier transform technique. In Fig. 4, we present a linear-logarithmic plot of the correction factor of the mobility function versus the radius of the disks for various values of the singularity position. Results are obtained by integrating Eq. (3.3) numerically. We observe that the curves follow a sigmoid-logistic-like phenomenology, implying that the correction factor increases significantly in the range of small radii before it reaches a saturation value. The latter corresponds to the correction factor predicted near two infinitely large disks given by Eq. (3.4).
Next, in order to quantify the effect of finite disk size on the correction to the hydrodynamics mobility, we customarily define the radius R 99 for which the mobility near infinitely large disks is essentially reached, such that k(R 99 ) = 0.99k ∞ . In the inset of Fig. 4, we display the variations of R 99 versus h based on the data presented in the main plot. We observe that R 99 reaches a maximum value of about 0.62 at the mid-plane of the channel before it monotonically decreases with h. This observation suggests that, to a good approximation, the mobility near two infinitely large disks can adequately be used to estimate the mobility at arbitrary position along the axis provided that the radius-to-channel-height ratio is above 0.62. Hence, accounting for the finite-size effect here becomes crucial only for values below this threshold.
Finally, we comment on the applicability of the often-used approximation originally suggested by Oseen (1928) to predict the particle mobility between two boundaries by superimposing separately the leading-order effects of each boundary. Accordingly, where the leading-order correction to the mobility function for axisymmetric motion normal to one rigid circular disk has previously been obtained by Kim (1983) and is expressed by wherein ξ = b/R is a dimensionless parameter with b denoting the distance between the particle and the centre of the disk. This solution was obtained by formulating the flow problem in terms of a mixed-boundary-value problem and solving the resulting dual integral equations using an approach analogous to that employed in the present work. Notably, we recover for ξ → 0 the familiar correction to the hydrodynamic mobility near an infinitely extended plane solid wall of no-slip boundary condition at its surface, namely ∆µ Disk /µ 0 = −9a/(8b), as originally obtained by Lorentz using the reciprocal theorem more than a century ago (Lorentz 1907;Lee et al. 1979). We now assess the accuracy of the superposition approximation stated by Eq. (3.7) by direct comparison with the exact prediction given by Eq. (3.3). In Fig. 5, we plot the variations of the percentage relative error between the correction factors k Sup and k versus the radius of the disks R for various values of the particle position h. In the range of small values of R, the relative error amounts to small values, typically smaller than 10% for R < 0.1. Upon increasing R, the relative error gradually increases in a logistic-like manner, before it saturates on a plateau value as R gets larger. The maximum error is obtained for the particle located on the mid-plane between the two disks for h = 0 and is found to be of about 55% in the limit of infinite disk radius. Therefore, the superposition approximation cannot be applied properly in this case. Nonetheless, as the particle position gets closer to either disk, the maximum error notably decreases to amount to only about 12% for h = 0.4. Consequently, the superposition approximation can frequently be utilised in this range of values to predict the hydrodynamic mobility for axisymmetric motion along the axis of the disks.

Conclusions
To summarise, we have examined the axisymmetric Stokes flow resulting from a Stokeslet singularity acting on the axis of two coaxially positioned circular disks of equal radius. We have formulated the solution for the viscous incompressible flow field as a mixed-boundary-value problem, which we have then reduced into a system of dual integral equations for four unknown wavenumber-dependent functions. Most importantly, we have shown the existence of viscous toroidal eddies in the fluid region bounded by the two plates. In the limit of infinitely large disks, we have successfully recovered the classic solution by Liron & Mochon (1976) for a Stokeslet acting normal to two parallel planar walls.
Additionally, we have provided an integral expression of the hydrodynamic mobility function quantifying the effect of the confining plates on the motion of a point-like particle moving along the axis of the coaxially positioned disks. Furthermore, we have demonstrated that accounting for the finite-size effect of the disks becomes essential only below a threshold value of the radius-to-channel height. Beyond this value, the mobility near two infinitely large disks can appropriately be employed. Finally, we have tested the validity and robustness of Oseen's approximation that postulates that the particle mobility between two boundaries could approximately be predicted by superimposing the contributions from each boundary independently. We have found that this simplistic approximation works quite well as the particle gets closer to either boundary but severely breaks down when the particle is located in the mid-plane between the two disks.
The analytical approach in the present article is based on the assumption of flow axisymmetry. The Stokes flow induced by a Stokeslet directed along an arbitrary direction in the presence of two coaxially positioned disks would be worth being investigated in a future study. We conjecture that this solution might be obtained by making use of the Green and Neumann functions supplemented by the edge function following the approach by Miyazaki (1984). This solution can then be employed to evaluate the translational and rotational mobility functions of particles located at arbitrary positions between the two disks. Alternatively, the problem can possibly be approached differently by means of multipole expansion methods involving the expression of the relevant hydrodynamic fields using oblate spheroidal coordinates (Lee & Leal 1980). This approach has widely been employed in the context of micromechanics of heterogeneous composite materials and fracture analysis (Kushch & Sangani 2000;Kushch 2013). In principle, our calculations can be extended to account for higher-order correction factors in the aspect ratio between the radius of the disks and the distance between the particle and the bounding plates (Swan & Brady 2010), but this would require a very challenging effort.
For applications requiring the precise manipulation of single molecules at the nanoscale level, the no-slip boundary condition may need to be lifted. In this context, the effect of partial slip at the surfaces of the disks is commonly characterised by assuming that the velocity components of the fluid tangent to the surfaces of the disks is proportional to the rate of strain at the surfaces (Lauga & Squires 2005;Lasne et al. 2008). This is an interesting aspect that could be included in our formalism and represents a worthwhile extension of the problem for future studies. We hope that our study will prove useful to researchers as well as practitioners working on particulate flow problems involving finitely sized boundaries, and pave the way toward better design and control of various processes in micro-and nanofluidic systems.
Declaration of Interests. The authors report no conflict of interest.

Appendix A. Analytical expressions for the kernel functions
In this Appendix, we provide technical details regarding the analytical derivation of the kernel functions appearing in the system of Fredholm integral equations of the first kind given by Eqs. (2.22) of the main body of the article.
The kernel functions can be expressed as infinite integrals over the wavenumber λ as It can be shown that the first four integrals can conveniently be expressed in closed mathematical forms as where we have defined the abbreviations In addition, the integrals L 5 and L 6 have analytical forms and can be calculated directly form standard integration tables or software algebra systems such as Mathematica (Wolfram 1999) as where H(·) denotes the Heaviside step function.
In the following, we will show how the integrals given by Eqs. (A 1) can be evaluated analytically. The core idea of our approach consists of expressing these integrals in the form of Laplace transforms of Bessel functions of the first kind (Spiegel 1965;Widder 2015), and using the recurrence relation (Abramowitz & Stegun 1972) 2k In addition, we will employ the following identities providing closed-form expressions for the Bessel functions of the first kind of half-integer order in terms of the standard trigonometric functions, Evaluation of the integral L 1 By making use of the identity given by Eq. (A 7a), the integral L 1 stated by Eq. (A 1a) can be expressed as Using the change of variable x = λr and Euler's representation of the sine function, the latter integral can be expressed as This leads to Eq. (A 2a) after making use of the Laplace transform given by Eq. (A 5) for k = 1 and p = (1 − it)/r. We note that Im(z) = − Im(z) for z ∈ C, wherez denotes the complex conjugate of z.
Evaluation of the integral L 4 Finally, upon using Eq. (A 6) for k = 1/2 and the identities given by Eqs. (A 7), the integral L 4 can be decomposed into four parts In the following, we will make use when appropriate of the shorthand notations defined in Eq. (A 3a). By using the integral representation of the sine function given by Eq. (A 12), the first integral can be expressed as Next, the evaluation of the second integral is straightforward after expressing the firstorder Bessel function as a function of the zeroth and second order Bessel functions using the recurrence relation given by Eq. (A 6) for k = 1 to obtain L 4,2 (r, t) = r (2πt) Lastly, the fourth integral can be decomposed into two parts as L 4,4 (r, t) = L 4,4,1 (r, t) + L 4,4,2 (r, t) , (A 27) where L 4,4,1 (r, t) = (2α) −1 L 2,2 (r, t) and L 4,4,2 (r, t) = (2πt) − 1 2 r t ∞ 0 λ −1 e −λ J 2 (λr) sin (λt) dλ .