Axisymmetric flow due to a Stokeslet near a finite-sized elastic membrane

Axisymmetric flow due to a Stokeslet near a finite-sized elastic membrane

Abdallah Daddi-Moussa-Ider Institut für Theoretische Physik II: Weiche Materie, Heinrich-Heine-Universität Düsseldorf, Universitätsstraße 1, 40225 Düsseldorf, Germany    Badr Kaoui Biomechanics and Bioengineering Laboratory (UMR 7338), CNRS, Université de Technologie de Compiègne, 60200 Compiègne, France    Hartmut Löwen Institut für Theoretische Physik II: Weiche Materie, Heinrich-Heine-Universität Düsseldorf, Universitätsstraße 1, 40225 Düsseldorf, Germany
July 20, 2019

Elastic confinements play an important role in many soft matter systems and affect the transport properties of suspended particles in viscous flow. On the basis of low-Reynolds-number hydrodynamics, we present an analytical theory of the axisymmetric flow induced by a point-force singularity (Stokeslet) directed along the symmetry axis of a finite-sized circular elastic membrane endowed with resistance toward shear and bending. The solution for the viscous incompressible flow surrounding the membrane is formulated as a mixed boundary value problem, which is then reduced into a system of dual integral equations on the inner and outer sides of the domain boundary. We show that the solution of the elastohydrodynamic problem can conveniently be expressed in terms of a set of inhomogeneous Fredholm integral equations of the second kind with logarithmic kernel. Basing on the hydrodynamic flow field, we obtain semi-analytical expressions of the hydrodynamic mobility function for the translational motion perpendicular to a circular membrane. The results are valid to leading-order in the ratio of particle radius to the distance separating the particle from the membrane. In the quasi-steady limit, we find that the particle mobility near a finite-sized membrane is always larger than that predicted near a no-slip disk of the same size. We further show that the bending-related contribution to the hydrodynamic mobility increases monotonically upon decreasing the membrane size, whereas the shear-related contribution displays a minimum value when the particle-membrane distance is equal to the membrane radius. Accordingly, the system behavior may be shear or bending dominated, depending on the geometric and elastic properties of the system. Our results may find applications in the field of nanoparticle-based sensing and drug delivery systems near elastic cell membranes.


I Introduction

Geometric confinements are commonly encountered in soft matter systems and in particular affect the behavior and transport properties of colloidal suspensions in a viscous medium Waigh (2016); Graham (2011); Yamamoto et al. (2008); Nakayama and Yamamoto (2005); Diamant (2009); Molina, Nakayama, and Yamamoto (2013). Hydrodynamic interactions between nanoparticles and elastic cell membranes play a key role in a large variety of biological and technological applications. A notable example being drug delivery and targeting using nanoparticle carrier systems, which navigate through blood vessels to reach disease sites such as tumors and inflammation areas Veiseh, Gunn, and Zhang (2010); Naahidi et al. (2013); Al-Obaidi and Florence (2015); Liu et al. (2016). During their uptake by living cells via endocytosis Doherty and McMahon (2009); Meinel et al. (2014); Agudo-Canalejo and Lipowsky (2015), the behavior of nanoparticles is strongly affected by hydrodynamic interactions with living cells composing of lipid and protein membranes.

In low-Reynolds-number hydrodynamics, the presence of nearby interfaces is known to drastically modify the flow field around immersed objects because of the long-range nature of the fluid-mediated hydrodynamic interactions. Over the last few decades, considerable research efforts have been devoted to the theoretical study of the slow (creeping) motion of a small particle moving close to a rigid boundary MacKay and Mason (1961); Gotoh and Kaneda (1982); Cichocki and Jones (1998); Lauga and Squires (2005); Swan and Brady (2007); Franosch and Jeney (2009); Swan and Brady (2010); Felderhof (2012); De Corato et al. (2015); Huang and Szlufarska (2015); Rallabandi, Hilgenfeldt, and Stone (2017), a fluid-fluid interface separating two immiscible viscous fluids Lee, Chadwick, and Leal (1979); Berdan II and Leal (1982); Bickel (2006, 2007); Bławzdziewicz, Ekiel-Jeżewska, and Wajnryb (2010a, b), or a soft membrane endowed with surface elasticity Felderhof (2006); Shlomovitz et al. (2014); Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a, b); Daddi-Moussa-Ider and Gekle (2016); Daddi-Moussa-Ider et al. (2018a), finding that particle diffusion perpendicular to the interface is significantly hindered compared to that along the direction parallel to the interface.

Unlike a solid-liquid or a liquid-liquid interface, an elastic interface stands apart as it introduces a memory effect in the system that causes a long-lasting anomalous subdiffusive behavior on nearby particles Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a). In addition, a rigid object that is dragged tangent to an elastic wall in a viscous fluid experiences a lift force that is directed normal to its direction of motion Skotheim and Mahadevan (2004, 2005); Yin and Kumar (2005); Snoeijer, Eggers, and Venner (2013); Wang, Dhong, and Frechette (2015); Wang et al. (2017); Wang, Tan, and Frechette (2017). This lift mechanism is a direct consequence of the elastic deformation of the interface, which breaks the time-reversal symmetry of Stokes flows. Experimentally, particle motion near confining interfaces has been investigated using optical tweezers Faucheux and Libchaber (1994); Dufresne, Altman, and Grier (2001); Schäffer, Nørrelykke, and Howard (2007), video microscopy Dufresne et al. (2000); Eral et al. (2010); Cervantes-Martínez et al. (2011); Dettmer et al. (2014), or evanescent wave dynamic light scattering Holmqvist, Dhont, and Lang (2007); Michailidou et al. (2009); Lisicki et al. (2012); Rogers et al. (2012); Michailidou et al. (2013); Lisicki et al. (2014). Meanwhile, the influence of a nearby elastic interface has been investigated using magnetic particle actuation Irmscher et al. (2012), optical traps Shlomovitz et al. (2013); Boatwright et al. (2014); Jünger et al. (2015); Tränkle, Ruh, and Rohrbach (2016), or rotating coherent scattering microscopy Jünger, Olshausen, and Rohrbach (2016).

In this paper, we calculate theoretically the axisymmetric flow field induced by a Stokeslet situated along the axis of a finite-sized circular membrane possessing resistance toward shear and bending. The solution of the flow equations is formulated as a mixed boundary value problem, which is then reduced to a set of convergent Fredholm integral equations of the second kind. We then derive semi-analytical expressions of the frequency-dependent mobility function that relates the velocity of a particle moving near a membrane to the hydrodynamic force exerted on its surface. More prominently, we show that the system behavior may be dominated by shear or bending, and that depending on the membrane size and the distance separating the particle from the membrane.

The remaining of the paper is organized as follows. In Sec. II, we state the elastohydrodynamic problem and introduce a relevant model for the membrane that incorporates both shear and bending deformation modes. We then formulate in Sec. III the solution of the flow equations and derive the corresponding mixed boundary value problem which we solve in Sec. IV for idealized membrane with pure shear or pure bending. In Sec. V, we provide semi-analytical expressions of the hydrodynamic mobility functions. Concluding remarks are contained in Sec. VI. Further technical details which are not essential to the understanding of results are relegated to appendices.

Ii Mathematical formulation

Figure 1: Cross section of the system setup. A solid spherical particle of radius  translating perpendicular to a finite-sized elastic membrane of radius , under the action of an external force . The surrounding fluid is Newtonian characterized by a constant dynamic viscosity . The membrane is composed of a hyperelastic material that possesses resistance toward shear elasticity and bending stiffness.

We consider the axisymmetric motion of a solid spherical particle of radius , initially located a distance  above a finite-sized membrane of radius  that is extended in the  plane. The unit vector  is directed normal to the undeformed membrane. The particle is moving under the action of an arbitrary time-dependent external force , as schematically illustrated in Fig. 1. Here, we investigate the system behavior in the far-field limit, where the particle size is small compared to the particle-membrane distance. The fluid is Newtonian of constant dynamic viscosity  and the flow is assumed to be incompressible. The membrane is modeled as a two-dimensional hyperelastic circular disk endowed with resistance toward shear and bending. The shear elasticity of the membrane is described by the Skalak model Skalak et al. (1973), which incorporates both the resistance toward shear and area dilation. The Skalak model is characterized by the elastic shear modulus  and the area expansion modulus , both of which are related by the dimensionless number . The latter is commonly known as the Skalak parameter Kaoui, Krüger, and Harting (2012); Freund (2014), and is typically very large for red blood cells, so as to express the area incompressibility of cell membranes. The resistance of the membrane toward bending is modeled by the well-established Helfrich model, which is commonly used as a relevant model for bilayer lipid vesicles or liposomes.

For a small displacement of the membrane relative to a horizontal plane of reference, the traction jump equations arising from these two deformation modes are expressed in a linearized form by Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a, b); Daddi-Moussa-Ider and Gekle (2016)


where denotes the membrane displacement field and is the traction jump vector across the membrane. Moreover, is the Laplace-Beltrami operator Deserno (2015), which, for a given function , is defined in Cartesian coordinates as . In addition, is the dilatation function, and denotes the position vector of the material points of the membrane in the planar configuration of reference. We note that a comma in indices represents a partial spatial derivative with respect to the corresponding variable.

Assuming low-Reynolds-number hydrodynamics Leal (1980); Happel and Brenner (2012); Kim and Karrila (2013), the fluid velocity and viscous stress fields, respectively denoted as  and , satisfy the stationary Stokes equations


wherein  is the force density acting on the surrounding fluid due to the presence of the suspended particle. The total force  is obtained by integrating over the particle surface. The viscous stress tensor is where denotes the pressure field and is the rate-of-strain tensor. Moreover, the components of the traction jump vector appearing on the right-hand side of Eqs. (1), are related to the fluid stress tensor by


It is worth mentioning that we have omitted the unsteady term in the Stokes equations since, in realistic situations, it yields a negligible contribution to the induced flow field near elastic interfaces Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a).

In the following, we approximate the hydrodynamic force by its first multipole moment such that . Accordingly, the elastohydrodynamic problem reduces to solving the governing equations of fluid motion for a point-force singularity positioned at the center of the particle and located above the planar membrane at position .

Since the flow is axisymmetric, the problem can more conveniently be solved using cylindrical coordinates. The stationary Stokes equations governing the fluid motion in the cylindrical coordinate system read Kim and Karrila (2013)


where  is the axisymmetric Laplace operator, which, for a function  expressed in cylindrical coordinates, is given by .

Iii Formulation of the solution of the flow problem

iii.1 Solution form

Due to the linearity of the Stokes equations, the solution of the flow problem can be written as a linear superposition of the free-space Green’s function and a complementary solution that is required to satisfy the boundary conditions prescribed at the elastic membrane. In an unbounded fluid, the radial and axial velocities due to a Stokeslet are Kim and Karrila (2013)


where and is the distance from the Stokeslet position. In the far-field limit, the flow velocity field decays as . This implies that hydrodynamic interactions are long ranged and can strongly be altered by geometric confinements. The corresponding Stokeslet solution for the pressure field reads


In the presence of a confining interface, the solution of the flow problem for the velocity and pressure fields can be presented in the form


where and are the complementary (image system) solution needed to satisfy the boundary conditions imposed at the membrane. For an axisymmetric flow, a convenient solution form of the stationary Stokes equations has been given by Imai Imai (1973), and can be expressed as Kim (1983)


wherein and are two harmonic functions, satisfying the axisymmetric Laplace equation, to be determined from the underlying boundary conditions. Far away from the singularity position, it can readily be checked that the solution form given above satisfies the governing equations stated by Eqs. (4).

We now denote by  and  the unknown harmonic functions in the upper-half space for , and by and  the corresponding functions in the lower half-space for . The general solution of the harmonic equations  and  in cylindrical coordinates for an axisymmetric problem can conveniently be expressed in terms of infinite integrals over the wavenumber , as Watson (1995)


where  denotes the th order Bessel function of the first kind Abramowitz and Stegun (1972). Moreover, the functions and are unknown wavenumber-dependent functions to be determined from the natural continuity of the velocity field across the membrane, together with the discontinuity of the hydrodynamic stresses, as derived from the shear and bending properties of the membrane.

By making use of Eqs. (8) and (9), and interchanging the derivative and integral operators, the solution for the velocity and pressure fields is given by


It worth mentioning that, since the problem is two dimensional, the equations of fluid motion can also be expressed in terms of the stream function . Accordingly, the solution of the flow problem could be reduced to the search of a single scalar function instead of simultaneously solving for the velocity and pressure fields. The stream functions in the upper- and lower-half spaces can be presented in the following form


where the first term corresponds to the Stokeslet solution, whereas the integral term corresponds to the image system solution. Accordingly, the radial and axial components of the flow velocity can, respectively, be calculated as


iii.2 Derivation of the mixed boundary value problem

Having presented the general form solution of the axisymmetric problem, we next determine the expressions of the unknown functions  and from the underlying boundary conditions. In the following, we will demonstrate that the solution of the elastohydrodynamic problem can be formulated as a mixed boundary value problem, which can then be reduced into a system of dual integral equations.

The continuity of the radial and axial velocity components at the plane of reference yields and . As a result, it follows that .

We now define, for convenience, the wavenumber-dependent functions and . We will show in the sequel that and are in fact functions associated with the shear and bending deformation modes, respectively.

On the one hand, the tangential and normal tractions across the membrane are continuous for . This results into the following integral equations for the outer boundary


On the other hand, membrane resistance toward shear and bending introduces jumps in the traction vector across the membrane for . The traction jump equations given by Eqs. (1), are expressed in cylindrical coordinates as


wherein and  are the radial and axial displacements of the material points of the membrane relative to their initial positions in the undeformed planar state. In an axisymmetric problem, the membrane displacement field is a function of the radial distance only.

In order to accomplish a closure of the present elastohydrodynamic problem, we require a relationship between the velocity and displacement fields. For that purpose, we assume, for simplicity, a no-slip boundary condition at the undisplaced membrane. Accordingly, the velocity of the fluid at is supposed to be identical to that of the displaced material points of the membrane. The effect of partial slip can be explored in future studies. Mathematically, the no-slip condition reads Bickel (2006, 2007)


If the membrane undergoes a large deformation, the no-slip boundary condition stated above should rather be applied at the displaced membrane, see for instance Refs. Sekimoto and Leibler, 1993; Weekley, Waters, and Jensen, 2006; Salez and Mahadevan, 2015; Saintyves et al., 2016; Rallabandi et al., 2017; Daddi-Moussa-Ider et al., 2018b; Rallabandi et al., 2018. Since we restrict our attention here to the system behavior in the small deformation regime such that , applying the no-slip condition at the undisplaced plane of the membrane should be sufficient for our investigations.

The resulting mathematical problem can conveniently be solved using Fourier transforms in time. Here, we employ the usual convention of a negative exponent for the forward Fourier transform. Transforming Eq. (15) into Fourier domain yields


By making use of Eq. (16), the traction jump equations across the elastic membrane given by Eqs. (14) can be expressed as




are characteristic length scales for shear and bending, respectively, as previously defined in earlier works Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a, b); Daddi-Moussa-Ider and Gekle (2016). Moreover, is a dimensionless number associated with the Skalak model.

Upon substitution of the general solution for the velocity and pressure fields given by Eqs. (10) into Eq. (17), the following integral equations for the inner boundary are obtained


Here, we have defined, for convenience, the frequency-dependent radial functions


where, the subscripts S and B denote shear and bending, respectively.

iii.3 Solution for an infinite membrane

Before proceeding with the solution of the flow problem near a finite-sized membrane, we first recall the solution for the axisymmetric Green’s function near an infinitely-extended elastic membrane. This solution has been derived earlier in closed form by some of us, by means of a two-dimensional Fourier transform technique Bickel (2007), see, for instance, Ref. Daddi-Moussa-Ider, Guckenberger, and Gekle, 2016a for more details regarding the derivation. In the present framework, the solution for an infinite membrane can be recovered by taking the limit in Eqs. (19) as goes to infinity. Accordingly, the shear- and bending-related functions  and  can readily be obtained by performing the inverse Hankel transforms Ditkin and Prudnikov (1965), to obtain


Using the relations and , it follows that


which are in full agreement with Ref. Daddi-Moussa-Ider, Guckenberger, and Gekle, 2016a. As  and are both taken to infinity (corresponding to a membrane with infinite shear and bending moduli, or to a vanishing actuation frequency), these expressions reduce to the well-known solution by Blake Blake (1971), for a Stokeslet acting normal to an infinitely-extended planar hard wall.

Iv Solution of the mixed boundary value problem

Due to the decoupled nature of the shear and bending deformation modes, it appears to be more convenient to solve the elastohydrodynamic problem by considering the effects of shear and bending independently. The overall flow field is the superposition of the individual flow fields resulting from these two deformation modes. Nevertheless, this decoupling behavior only occurs near a single planar elastic membrane. A nonlinear coupling between shear and bending has been observed for curved elastic interfaces Daddi-Moussa-Ider and Gekle (2017); Daddi-Moussa-Ider, Lisicki, and Gekle (2017a, b, 2018), or for two closely coupled Auth, Safran, and Gov (2007) or warped Košmrlj and Nelson (2014) fluctuating membranes.

iv.1 Shear contribution

For an idealized membrane with pure shear elasticity, the corresponding dual integral equations reads


where is the radially-symmetric function given explicitly by Eq. (20a). The goal is to solve for the unknown function  such that the dual integral equations on both the inner and outer boundaries are satisfied.

To reduce the order of the Bessel function and bring the equations into a more familiar form, we multiply both members by  and differentiate with respect to  afterwards. Consequently, the dual integral equation can be rewritten as


where we have defined

The solution of the resulting dual integral equations can be obtained following the resolution recipes given by Sneddon Sneddon (1966). The starting point consists of writing the solution in the following integral form


where  is an unknown function to be determined. This choice is motivated by the fact that the integral equation on the outer boundary is automatically satisfied. In fact, by substituting the form solution given above into Eq. (24b), and interchanging the order of integration, one obtains


after making use of the identity Abramowitz and Stegun (1972)


where denotes the Heaviside function.

Substituting the form solution given by Eq. (25) into Eq. (24a) associated with the inner boundary, and interchanging the order of integration, yields


On the one hand, it follows from Eq. (27) that


On the other hand, an integration by parts yields


upon making use of the identity Abramowitz and Stegun (1972)


By inserting Eqs. (29) and (30) into the integral equation given by Eq. (28), one gets


After some rearrangements and mathematical manipulations, the latter equation can be presented in the form of an inhomogeneous Fredholm integral equation of the second kind with a logarithmic kernel. Specifically,


Further derivation details are contained in Appendix A. Remarkably, the solution at the origin is trivial by noting that .

Due to the somewhat complex nature of the resulting integral equation at hand, an analytical solution is far from being trivial. Therefore, recourse to a numerical resolution approach is necessary. For that purpose, we express the solution as a finite series of terms in powers of  as


and sample the result at Chebyshev nodes. These correspond to the roots of the Chebyshev polynomial Mason and Handscomb (2002) of the first kind of degree . We then perform the integration analytically and solve the resulting finite linear system of equations for the unknown complex coefficients . Inserting the form series solution given by Eq. (34) into Eq. (33) and interchanging between the sum and the integral operators, the resulting linear system of equations reads


where we have defined the dimensionless parameter  representing the ratio of the particle-membrane distance to the radius of the membrane. In particular, corresponds to the infinite-size limit, whereas  holds for a bulk fluid, i.e., in the absence of the membrane. The Chebyshev nodes are mapped over the interval and are given by


Moreover, the series functions  are defined by


the expressions of which are explicitly given by

The latter result can readily be shown by recurrence over . We further mention that the series  can also be expressed for arbitrary parity of  in terms of the Lerch transcendent function Abkarian, Faivre, and Viallat (2007) (which is implemented in computer algebra systems such as Mathematica or Maple as LerchPhi). Specifically,

The expression of for are provided for convenience in Tab. 1.

Table 1: Analytical expressions of the first eleven terms of the series functions  defined in the main text by Eq. (37). Here, corresponds a Chebyshev node mapped over the interval as given by Eq. (36).

iv.2 Bending contribution

We next consider the bending-related contribution to the flow field and search for solution of the dual integral equations


where  is the unknown function, and is the known radially-symmetric function given by Eq. (20b). The dual integral equations have the familiar form encountered in various mixed boundary value problems where the kernel is expressed in terms of the zeroth-order Bessel function.

Similarly, we search a solution of the integral form


where we assume that . This form solution satisfies the integral equation on the outer boundary after making use of Eq. (27). Performing three successive integrations by parts yields


where, once again, we have made use of the identity given by Eq. (31). We have further assumed that vanishes at the domain boundary, so as to ensure convergence of the overall improper integral.

By making use of Eqs. (27) and (40), the integral equation on the inner boundary becomes


The latter can be written in the form of an inhomogeneous Fredholm equation of the second kind by using the same resolution procedure described in Appendix A. Integrating the resulting equation with respect to  twice yields


where , and are four unknown integration constants to be determined from the imposed boundary conditions, namely . The function appearing in the kernel of the integral is given by


Since an analytical solution is far from being intuitive, we attempt a numerical evaluation of the integral equation by considering a solution of the form


and sample the result at Chebyshev nodes, in the same way as previously done for the shear-related part. After some rearrangement, we obtain


where we have defined the series functions


the general term of which is given by




wherein . Moreover, denotes the ceiling function, which maps a variable  to the least integer greater than or equal to . Depending on the parity of , the general term of  and are, respectively, expressed by




Expressions for become rather lengthy and cumbersome as  gets larger, and thus have not been provided here. These can more advantageously be obtained using computer algebra systems.

V Hydrodynamic mobility

The calculation of the flow field induced by a point-force acting near the elastic membrane can be employed to assess the effect of the interface on the slow motion of a suspended particle moving in its vicinity. This membrane-induced effect is quantified by the hydrodynamic mobility functions, which are tensorial quantities that bridge between the velocity of a sedimenting particle and the force exerted on its surface Batchelor (1976). In a bulk fluid, the particle mobility is isotropic and may be obtained from by the Stokes resistance formula as , where and denotes the Kronecker symbol. The presence of the elastic membrane introduces a correction to the hydrodynamic mobility that is dependent not only on the particle size and the viscosity of the suspending fluid, but also on the distance separating the particle from the membrane, as well as on the forcing frequency in the system.

In the following, we present analytical expressions for the frequency-dependent mobility of a spherical particle translating perpendicular to a finite-sized membrane possessing either shear elasticity or bending rigidities. The mobility correction near a membrane endowed simultaneously with both shear and bending deformation modes can be obtained by linear superposition of the shear- and bending-induced corrections as obtained independently. Again, this is true only for a single planar membrane where a decoupling between shear and bending deformation modes exists Daddi-Moussa-Ider, Guckenberger, and Gekle (2016b).

v.1 Shear contribution

For an idealized elastic membrane with only-shear resistance, such as an artificial capsule designed for drug delivery Barthès-Biesel (2016), the bending-related part in the solution of the flow problem vanishes. Thus, and . In the point-particle approximation, the scaled correction to the particle mobility at leading order for the motion perpendicular to the finite-sized membrane can directly be calculated by evaluating the correction to the axial component of the flow field at the particle position. Specifically,


As shown in the previous section, an analytical solution for  is not available due to the complicated nature of the resulting Fredholm integral equation. Therefore, a numerical technique has been employed in order to overcome this difficulty. Substituting the form series solution given by Eq. (34) into Eq. (25) and interchanging between the sum and the integral operator, yields


By inserting the latter equation into Eq. (51) and performing integration with respect to the variable , the shear-related part in the correction to the frequency-dependent mobility can be obtained in a scaled form as


where we have defined the series functions


which, depending on the parity of , is explicitly given by