2.1. Tangential squirmers

The particle is modelled as a time independent tangential squirmer. The fluid velocity at the squirmer's surface is (Ishikawa et al. Reference Ishikawa, Simmonds and Pedley2006, p. 156)

(2.2)\begin{equation} \boldsymbol{u}^{sl}(\boldsymbol{e}_{s} ) = B^{sq}_{1} (1+ \beta (\boldsymbol{e}_{s}\boldsymbol{\cdot} \hat{\boldsymbol{r}}^{sq}) ) \left[ (\boldsymbol{e}_{s} \boldsymbol{\cdot} \hat{\boldsymbol{r}}^{sq})\hat{\boldsymbol{r}}^{sq} - \boldsymbol{e}_{s} \right], \end{equation}

and may be obtained by keeping the first two expansion coefficients from the full representation of the squirmer's tangential surface deformation field. The squirmer moves or swims in the direction of its orientation, $\boldsymbol {e}_{s}$. The ratio, $\beta = B^{sq}_{2} / B^{sq}_{1}$, describes the type of squirmer, the coefficients, $B^{sq}_{1}$ and $B^{sq}_{2}$, are referred to as squirming modes, and $\hat {\boldsymbol {r}}^{sq} = \boldsymbol {r}^{sq}/a$ is the radial unit vector pointing from the centre of the squirmer to the squirmer's surface. Far away from the squirmer, the velocity field created by the first mode, $B^{sq}_{1}$, is that of a source dipole and decays as $O(r^{-3})$. The velocity field created by the second mode, $B^{sq}_{2}$, is that of a force dipole and decays as $O(r^{-2})$. These two squirming modes have been shown to sufficiently capture the far field dynamics of a variety of biological microswimmers (Spagnolie & Lauga Reference Spagnolie and Lauga2012). In the literature, the sign of $B^{sq}_{2}$ is associated with the names pusher and puller for $B^{sq}_{2} <0$ and $B^{sq}_{2} > 0$. When $B^{sq}_{2} = 0$, the squirmer is termed neutral and becomes a source dipole swimmer.

Equation (2.2) expresses the fluid velocity at the squirmer's surface in a body fixed frame of reference attached to and moving with the squirmer. This surface slip velocity is axisymmetric relative to the squirmer's swimming orientation. The reciprocal theorem (Stone & Samuel Reference Stone and Samuel1996) provides an elegant method to calculate the squirmer's free space translational velocity, $\boldsymbol {U}_{fs}^{sq}$, from which it may be found that $\boldsymbol {U}_{fs}^{sq} = 2/3 B^{sq}_{1} \boldsymbol {e}_{s}$, implying that the first tangential squirming mode provides the only non-zero contribution to the squirmer's velocity.

3.1. The rigid non-porous container

In the case of the non-porous rigid container, the fluid mechanics problems in regions $\varOmega _{i}$ and $\varOmega _{e}$ are fully decoupled. The container is at rest with velocity $\boldsymbol {U}^{m} = \boldsymbol {0}$. The kinematic and dynamic boundary conditions on the surface of the container are given by the vanishing of the normal and tangential fluid velocity components and may be written as

(3.8a)\begin{gather} \left.\psi^{i}(r,\zeta)\right\rvert_{r=b} = 0, \end{gather}

(3.8b)\begin{gather}\left.\frac{\partial \psi^{i}(r,\zeta)}{\partial r}\right\rvert_{r=b} = 0. \end{gather}

Under a Cartesian world space basis (see figure 1) with $\boldsymbol {e}_{z}$ pointing upwards and $\boldsymbol {e}_{y}$ pointing into the page, the squirmer's swimming orientation, $\boldsymbol {e}_{s}$, points in the direction of $\boldsymbol {e}_{z}$. The translational velocity of the squirmer is $\boldsymbol {U}^{sq} = U_{z}^{sq}\boldsymbol {e}_{z}$. The kinematic and dynamic boundary conditions on the squirmer's surface are given by matching the squirmer's surface velocity to the fluid velocity at $r=a$ and may be written as

(3.9a)\begin{gather} \left. \left[ \boldsymbol{u}^{i}(r,\theta) - \left( \boldsymbol{U}^{sq} + \boldsymbol{\varOmega}^{sq} \times \boldsymbol{r} + \boldsymbol{u}^{sl}(r,\theta) \right) \right] \boldsymbol{\cdot} \boldsymbol{e}_{r}\right\rvert_{r=a} = 0, \end{gather}

(3.9b)\begin{gather}\left. \left[ \boldsymbol{u}^{i}(r,\theta) - \left( \boldsymbol{U}^{sq} + \boldsymbol{\varOmega}^{sq} \times \boldsymbol{r} + \boldsymbol{u}^{sl}(r,\theta) \right) \right] \boldsymbol{\cdot} \boldsymbol{e}_{\theta}\right\rvert_{r=a} = 0, \end{gather}

where $\boldsymbol {U}^{sq}$ and $\boldsymbol {\varOmega }^{sq}$ denote translational and angular velocities of the squirmer, and $\boldsymbol {u}^{sl}(r,\theta )$ is the squirmer's surface slip velocity defined in (2.2). Since the squirmer's surface velocity has no azimuthal or radial component, no net angular motion can be generated along $\boldsymbol {e}_{z}$, and so it must hold that $\boldsymbol {\varOmega }^{sq} = \boldsymbol {0}$. Using intrinsic coordinates on the surface of the squirmer and recalling the relationships (Happel & Brenner Reference Happel and Brenner1983, 4-5.3, 4-5.4)

(3.10a)\begin{gather} u^{i}_{r}(a,\theta) ={-}\frac{1}{a\sin\theta} \frac{\partial \psi^{i}(a,\theta)}{\partial \theta}, \end{gather}

(3.10b)\begin{gather}u^{i}_{\theta}(a,\theta) = \left.\frac{1}{a\sin\theta} \frac{\partial \psi^{i}(r,\theta)}{\partial r}\right\vert_{r=a}, \end{gather}

(3.10c)\begin{gather}\frac{\partial (r\sin\theta)}{\partial \theta} = \boldsymbol{e}_{r} \boldsymbol{\cdot} \boldsymbol{e}_{z}, \end{gather}

(3.10d)\begin{gather}\frac{\partial (r\sin\theta)}{\partial r} ={-}\boldsymbol{e}_{\theta} \boldsymbol{\cdot} \boldsymbol{e}_{z}, \end{gather}

the velocity boundary condition in (3.9a) simplifies as

(3.11)\begin{equation} \frac{\partial }{\partial \theta}\left( \psi^{i}(a,\theta) + \frac{1}{2}a^{2}U^{sq}_{z}\sin^{2}\theta \right) = 0. \end{equation}

After integrating around the squirmer in $\theta$, and setting the constant of integration to zero to be consistent with the general convention that streamfunctions vanish along the axis $\theta = 0$, (3.11) reduces to the expression

(3.12)\begin{align} \psi^{i}(r,\theta)\rvert_{r=a} &={-}\tfrac{1}{2}U^{sq}_{z}a^2\sin^{2}\theta,\nonumber\\ &={-}\tfrac{1}{2}U^{sq}_{z}a^2(1-\zeta^{2}). \end{align}

Similarly, using (3.10), (3.9b) simplifies as

(3.13)\begin{align} \left.\frac{\partial \psi^{i}(r,\theta)}{\partial r}\right\rvert_{r=a} &= a\sin^{2}\theta\left( B^{sq}_{1} (1+\beta\cos\theta) - U^{sq}_{z}) \right),\nonumber\\ &= a(1-\zeta^{2})\left( B^{sq}_{1} (1+\beta \zeta) - U^{sq}_{z}) \right). \end{align}

Based on the $\zeta$ dependence in the $\mathcal {J}_{2}$ and $\mathcal {J}_{3}$ eigenfunctions it is straightforward to see that the first two terms of (3.5) are sufficient to satisfy the boundary conditions on the squirmer's surface, (3.12), (3.13), and on the container's surface, (3.8a), (3.8b). Therefore, the streamfunction, $\psi ^{i}(r,\theta )$, defined over $\varOmega _{i}$, takes the form

(3.14) \begin{align} \psi^{i}(r,\theta) &= \left( A^{i}_{2}r^{2} + B^{i}_{2}r^{{-}1} + C^{i}_{2}r^{4} + D^{i}_{2}r\right)\mathcal{J}_{2}( \zeta)\nonumber\\ &\quad + \left(A^{i}_{3}r^{3} + B^{i}_{3}r^{{-}2} + C^{i}_{3}r^{5} + D^{i}_{3}\right)\mathcal{J}_{3}( \zeta). \end{align}

The procedure used to obtain the expansion coefficients, $\{A^{i}_{n},B^{i}_{n},C^{i}_{n},D^{i}_{n}\}$ for $n \in \{2,3\}$, in (3.14) is explained in appendix A.1. Solutions for these coefficients are provided in (A7). The expansion coefficients may be determined by using the orthogonality relations of the $\mathcal {J}_{n}$ polynomials from (3.7) on the four boundary conditions (3.8a), (3.8b), (3.12) and (3.13). Applying these orthogonality conditions over $\mathcal {J}_{2}( \zeta )$ and $\mathcal {J}_{3}( \zeta )$ yields eight linear equations in the eight unknown expansion coefficients. The solution of this system of equations yields expressions for the expansion coefficients in terms of $U^{sq}_{z}$. A force balance on the squirmer provides the final equation for determining $U^{sq}_{z}$.

In the absence of fluid inertia, a force balance on the neutrally buoyant, net force- and torque-free squirming particle reveals that the net hydrodynamic force and torque are zero

(3.15a)\begin{gather} \boldsymbol{F}^{H} = \boldsymbol{F}^{P} + \boldsymbol{F}^{D} = \boldsymbol{0}, \end{gather}

(3.15b)\begin{gather}\boldsymbol{T}^{H} = \boldsymbol{0}. \end{gather}

The net hydrodynamic force in (3.15a) is expressed as a sum of two separate contributions: a propulsive force $\boldsymbol {F}^{P}$ induced by the squirmer's slip velocity, and a drag force $\boldsymbol {F}^{D}$ induced by the rigid body motion of the squirmer through the fluid. Fundamentally, the squirmer generates directed motion by balancing the drag force with an internally generated equal and opposite propulsive force. Once the expansion coefficients of $\psi ^{i}(r,\theta )$ are determined, (3.15) may be used to solve for the translational velocity of the squirmer, $\boldsymbol {U}^{sq}$, subject to the formula (Happel & Brenner Reference Happel and Brenner1983, 4-14.18)

(3.16)\begin{equation} F^{H}_{z} = \mu {\rm \pi}\int_{\varGamma_{sq}}\left.\left[ (r \sin(\theta))^{3} \frac{\partial }{\partial r}\left( \frac{E^{2}\psi(r,\theta)}{(r \sin(\theta))^{2}} \right)r \right]\right\rvert_{r=a} \,\textrm{d}\theta. \end{equation}

Alternatively, since the squirmer is a sphere, the expression (Happel & Brenner Reference Happel and Brenner1983, 4-23.35)

(3.17)\begin{equation} F^{H}_{z} = 4 {\rm \pi}\mu D_{2}, \end{equation}

can be used to directly find $F^{H}_{z}$ from the $D_{2}^{i}$ coefficient.

Using the expansion coefficients from (A7) and solving either (3.16) or (3.17) together with net force- and torque-free constraints of (3.15), the squirmer's translational velocity is found to be

(3.18)\begin{equation} \boldsymbol{U}^{sq} = \frac{B^{sq}_{1} (b-a) \left(3 a^3+6 a^2 b+4 a b^2+2 b^3\right)}{3 \left(a^4+a^3 b+a^2 b^2+a b^3+b^4\right)} \boldsymbol{e}_{z}. \end{equation}

Equation (3.18) and the solution to the unit forced concentric sphere problem (Happel & Brenner Reference Happel and Brenner1983, pp. 130–133) are shown in figure 2. When compared with the forced particle–container solution, the squirmer–container solution shows much faster decay to its free space velocity, $\boldsymbol {U}_{fs}^{sq} = 2/3 B^{sq}_{1} \boldsymbol {e}_{z}$, as the container to particle size ratio increases. This faster decay is partially explained by the faster far field decay in the squirmer's disturbance velocity field when compared with the forced particle solution which decays as $O(r^{-1})$. The squirmer's slip velocity also provides an explicit mechanism by which fluid can move anti-parallel to the swimming orientation making the squirmer more efficient at moving in confinement than a forced particle.

Figure 2.

The non-dimensional centre of mass translational speeds of a forced particle, $U^{fp}_{z}$, and a squirmer, $U^{sq}_{z}$, are plotted as a function of the container–particle size ratio, $b/a$. The squirmer translational velocity approaches its free space squirming velocity rapidly in comparison to the forced particle as $b/a \to \infty$.

3.2. The rigid porous container

If the container is now made into a rigid porous membrane, the exterior and interior fluid regions become coupled across the membrane through the region $\varOmega _{p}$. The container now moves with velocity $\boldsymbol {U}^{m} = U^{m}\boldsymbol {e}_{z}$. Fluid transport in this region is modelled using a macroscopic approach, similar to Darcy's law; however, the porous region, $\varOmega _{p}$, is ultimately modelled as a thin permeable interface. The model for the porous container is partially motivated by the structure of Darcy's law which reads

(3.19)\begin{equation} \boldsymbol{u}^{p} ={-}\frac{\boldsymbol{\mathsf{K}}}{\mu} \boldsymbol{\nabla} p^{p}, \end{equation}

where $\boldsymbol{\mathsf{K}}$ is a second-order permeability tensor, and $p^{p}$, $\mu$ and $u^{p}$ are respectively the volume averaged pressure, viscosity and fluid filtration velocity. For a thin membrane, the gradient in pressure from (3.19) may be approximated using a first-order finite difference. For a finite second-order derivative in the pressure, this approximation is accurate to $O(l_{c})$ for membrane thickness $l_{c}$ under the assumption that $l_{c} \ll \{a,b,b-a\}$.

Instead of using $\boldsymbol{\mathsf{K}}$ it is more convenient to define a membrane resistance tensor, $\boldsymbol{\mathsf{R}}_{m} := l_{c} \boldsymbol{\mathsf{K}}^{-1}$, formulated as an inverse to permeability whose magnitude describes resistance to fluid transport across the membrane. Much like permeability, the membrane's resistance can be understood as an effective parameter, translating volume averaged fluid–solid quantities from a scale proportional to the membrane pore size, $a^{p}$, to a scale that may be several orders of magnitude larger.

Under a local spherical basis, $\{\boldsymbol {e}_{r}, \boldsymbol {e}_{\theta },\boldsymbol {e}_{\phi } \}$, the resistance tensor is diagonal with normal and tangential resistance parameters $R_{\perp } := R_{r}$ and $R_{\parallel } := R_{\theta } = R_{\phi }$; therefore, pressure driven flow normal to the membrane is controlled by $R_{\perp }$. Performing a normal stress balance on the membrane interface implies the following generalized condition:

(3.20)\begin{equation} -p^{e}+p^{i} + \mu ( \boldsymbol{e}_{r} \boldsymbol{\cdot} \boldsymbol{\tau}^{e} \boldsymbol{e}_{r} - \boldsymbol{e}_{r} \boldsymbol{\cdot} \boldsymbol{\tau}^{i} \boldsymbol{e}_{r}) = R_{{\perp}}\mu(\boldsymbol{u}^{p} - \boldsymbol{U}^{m}) \boldsymbol{\cdot} \boldsymbol{e}_{r}, \end{equation}

for viscous stress $\boldsymbol {\tau }$ and translational membrane velocity $\boldsymbol {U}^{m}$. Similar to (3.19), (3.20) states that a discontinuity in the normal stress drives the relative filtration velocity, $\boldsymbol {u}^{p} - \boldsymbol {U}^{m}$, normal to the membrane.

Since the order of the derivative in $\boldsymbol {u}^{p}$ is one less than in the Stokes equations, coupling Darcy's law to the Stokes equations requires making an additional assumption on a tangential boundary condition. If the vesicle were modelled as a Brinkman medium, an assumption on the tangential boundary condition would not be necessary. Rather, boundary conditions that express the continuity of the stress components between the free fluid in $\varOmega _{i}$ and $\varOmega _{e}$ and the Brinkman medium would be sufficient to create a well-posed system (Keh & Lu Reference Keh and Lu2005). When coupling to Darcy's law, no slip in the tangential direction is often the simplest condition to enforce, however, research on these boundary conditions has led to alternative tangential stress discontinuity conditions (Beavers & Joseph Reference Beavers and Joseph1967; Saffman Reference Saffman1971; Jones Reference Jones1973). Recently, Zampogna & Gallaire (Reference Zampogna and Gallaire2020) have developed a macroscopic model for a stress jump boundary condition across thin membranes. This model appears promising, however, to keep our analysis simple yet representative of a porous membrane, we have taken inspiration from the empirical model proposed by Beavers & Joseph (Reference Beavers and Joseph1967) which reads

(3.21)\begin{equation} \frac{\partial u}{\partial y} = \frac{\alpha}{\sqrt{k}}(u-q), \end{equation}

where $u$ is the mean velocity parallel to the surface, $q$ is the tangential volume flow rate inside the porous medium, $y$ is the coordinate normal to the surface, $\alpha$ is a dimensionless constant depending on the properties of the porous surface and $k$ is an isotropic permeability constant. Equation (3.21) implies that the fluid velocity at the porous interface does not have to be continuous and that the magnitude of the slip velocity is directly proportional to the shear stress. Jones (Reference Jones1973) solved for the flow past a spherical shell and proposed a generalization of (3.21) given by

(3.22)\begin{equation} \tau_{r\theta} = \frac{\alpha}{\sqrt{k}}(u_{\theta} - q_{\theta}). \end{equation}

The stress jump condition in (3.22) may become relevant when describing engineered membranes or vesicles featuring large pores, as tangential slip between the free fluid velocity and the velocity in the porous medium becomes much more physically plausible.

For a thin membrane, a tangential stress balance across the interface, taken together with (3.22), motivates the generalized condition

(3.23)\begin{equation} \mu(\boldsymbol{e}_{\theta} \boldsymbol{\cdot} \boldsymbol{\tau}^{e} \boldsymbol{e}_{r} - \boldsymbol{e}_{\theta} \boldsymbol{\cdot} \boldsymbol{\tau}^{i} \boldsymbol{e}_{r}) = R_{{\parallel}}\mu(\boldsymbol{u}^{p}-\boldsymbol{U}^{m})\boldsymbol{\cdot} \boldsymbol{e}_{\theta}, \end{equation}

where slip in the tangential component of the relative filtration velocity is driven by a discontinuity in the tangential stress at the interface. Equation (3.23) results in discontinuous partial derivatives

(3.24)\begin{equation} \left.\frac{\partial u^{i}_{\theta}}{\partial r}\right\rvert_{r=b} \ne \left.\frac{\partial u^{e}_{\theta}}{\partial r}\right\rvert_{r=b}, \end{equation}

and may be understood as one mechanism by which the squirmer's tangential slip velocity is transmitted to the container ultimately resulting in directed container motion.

The normal and tangential stress balances of (3.20) and (3.23) are supplemented with the conditions

(3.25a)\begin{gather} \boldsymbol{u}^{i}(r,\theta)|_{r=b}\boldsymbol{\cdot} \boldsymbol{e}_{r} = \boldsymbol{u}^{p}(r,\theta)|_{r=b}\boldsymbol{\cdot} \boldsymbol{e}_{r} = \boldsymbol{u}^{e}(r,\theta)|_{r=b}\boldsymbol{\cdot} \boldsymbol{e}_{r}, \end{gather}

(3.25b)\begin{gather}\boldsymbol{u}^{i}(r,\theta)|_{r=b}\boldsymbol{\cdot} \boldsymbol{e}_{\theta} = \boldsymbol{u}^{p}(r,\theta)|_{r=b}\boldsymbol{\cdot} \boldsymbol{e}_{\theta} = \boldsymbol{u}^{e}(r,\theta)|_{r=b}\boldsymbol{\cdot} \boldsymbol{e}_{\theta}, \end{gather}

which state that the fluid velocity is continuous at $r=b$.

The streamfunction, $\psi ^{i}(r,\theta )$, takes the same form as in (3.14) as the squirmer's boundary conditions have not changed. The conditions of (3.25) require that $\psi ^{e}(r,\theta )$ take the same form as $\psi ^{i}(r,\theta )$. The fluid velocity in $\varOmega _{e}$ decays to quiescence as $r \to \infty$ and requires $\lim _{r \to \infty }\psi ^{e}/r^{2} = 0$, implying that $A^{e}_{2}=C_{2}^{e}=A_{3}^{e}=C_{3}^{e}=0$. Therefore, the streamfunction, $\psi ^{e}(r,\theta )$, takes the form

(3.26)\begin{equation} \psi^{e}(r,\theta) = \left( B^{e}_{2}r^{{-}1} + D^{e}_{2}r\right)\mathcal{J}_{2}( \zeta) + \left( B^{e}_{3}r^{{-}2} + D^{e}_{3}\right)\mathcal{J}_{3}( \zeta). \end{equation}

Equations (3.14) and (3.26) collectively have 12 unknown expansion coefficients. The six boundary conditions, (3.12), (3.13), (3.20), (3.23), (3.25a) and (3.25b) are sufficient to determine $\psi ^{i,e}(r,\zeta )$ in terms of $U^{sq}_{z}$ and $U^{m}_{z}$ since they yield 12 equations after applying the orthogonality relation (3.7) over the $\mathcal {J}_{2}( \zeta )$ and $\mathcal {J}_{3}( \zeta )$ eigenfunctions. Force balances on the squirmer and container surfaces in combination with (3.16) or (3.17) are sufficient to determine $U^{sq}_{z}$ and $U^{m}_{z}$. Additional details pertaining to the structure and solution of these 12 linear equations are provided in appendix A.2.

Solutions for the 12 expansion coefficients from (3.14) and (3.26) are provided in (A23). The squirmer and membrane translational velocities are found to be

(3.27)\begin{gather}\boldsymbol{U}^{sq} = \frac{B^{sq}_{1} \left(3 a^5 R_{{\perp}} R_{{\parallel}}-5 a^3 b^2 R_{{\perp}} R_{{\parallel}}+2 b^5 R_{{\perp}} R_{{\parallel}}+10 b^4 (R_{{\perp}}+2 R_{{\parallel}})\right)}{3 \left((b^5-a^5) R_{{\perp}} R_{{\parallel}}+5 b^4 (R_{{\perp}}+2 R_{{\parallel}})\right)} \boldsymbol{e}_{z}, \end{gather}

(3.28)\begin{gather}\boldsymbol{U}^{m} = \frac{10 a^3 b B^{sq}_{1} (R_{{\perp}}-R_{{\parallel}})}{3 \left((b^5-a^5) R_{{\perp}} R_{{\parallel}}+5 b^4 (R_{{\perp}}+2 R_{{\parallel}})\right)}\boldsymbol{e}_{z}. \end{gather}

Both (3.27) and (3.28) are independent of $\beta$ and therefore $B^{sq}_{2}$ which indicates that all types of squirmers (pushers $\beta < 0$, neutral $\beta = 0$, and pullers $\beta >0$) move with the same translational speed.

The reason that (3.27) and (3.28) are independent of $B^{sq}_{2}$ is largely of mathematical consequence. As seen in (3.17), the $D_{2}^{i,e}$ coefficients determine the total hydrodynamic force on the squirmer and container, and are paired with the $\mathcal {J}_{2}( \zeta )$ eigenfunctions in (3.14). The streamfunction, $\psi ^{i}(r,\theta )$, must satisfy (3.13) on the surface of the squirmer. Since the $D^{i}_{2}$ coefficient is paired with the $\mathcal {J}_{2}( \zeta )$ eigenfunction, and the eigenfunctions are orthogonal, the $D^{i}_{2}$ expansion coefficient can only match the portion of the squirmer's velocity boundary condition that is proportional to $\sin ^{2}\theta$. The portion of the squirmer's tangential boundary condition that is proportional to $\sin ^{2}\theta$ is precisely independent of the $B^{sq}_{2}$ parameter. Therefore, the $B^{sq}_{2}$ will never be involved in the calculation of the hydrodynamic force on the squirmer. Since this force calculation is ultimately used to determine the squirmer's translational velocity, the $B^{sq}_{2}$ parameter will not contribute to the squirmer's velocity. In the same way, the streamfunctions $\psi ^{i,e}(r,\theta )$, have to match at $r=b$ which means that the $\mathcal {J}_{2}( \zeta )$ portions of $\psi ^{i,e}(b,\theta )$ must match. The $\mathcal {J}_{2}( \zeta )$ eigenfunction in $\psi ^{i}$ propagates the functional $B^{sq}_{1}$ dependence and $B^{sq}_{2}$ independence to the container's velocity.

3.3. Limiting behaviours

Several important limiting behaviours may be recovered from (3.27) and (3.28). Fixing the squirmer's size as $a$, if the size ratio $b/a \to \infty$, the squirmer's speed is affected less and less by confinement, and the free space squirming velocity, $2/3 B^{sq}_{1}$, is recovered with a container velocity that approaches zero. If $b/a \to 1$, then the squirmer perfectly transmits its boundary conditions to the container and the squirmer moves with the free space velocity, $2/3 B^{sq}_{1}$, with the container translating at a speed reduced by $(R_{\perp } - R_{\parallel })/(R_{\perp }+2R_{\parallel })$. Summarized more formally, for these two simple cases,

(3.29a,b)\begin{gather} \lim_{b\to \infty} \boldsymbol{U}^{sq} = \tfrac{2}{3}B^{sq}_{1}\boldsymbol{e}_{z},\quad \lim_{b\to \infty} \boldsymbol{U}^{m} =\boldsymbol{0}, \end{gather}

(3.30a,b)\begin{gather} \lim_{b\to a} \boldsymbol{U}^{sq} = \frac{2}{3}B^{sq}_{1}\boldsymbol{e}_{z},\qquad \lim_{b\to a} \boldsymbol{U}^{m} = \frac{2}{3}B^{sq}_{1}\left(\frac{R_{{\perp}} - R_{{\parallel}}}{R_{{\perp}}+2R_{{\parallel}}}\right)\boldsymbol{e}_{z}. \end{gather}

In the case where fluid travels through the membrane isotropically ($R_{\perp } = R_{\parallel }=R$) the container does not translate, but the squirmer translates with a speed

(3.31)\begin{equation} \lim_{R_{{\perp}}\to R_{{\parallel}}}\boldsymbol{U}^{sq} = \frac{B^{sq}_{1} \left(3 a^{5} R -5 a^{3} b^{2} R + 2 b^{5} R +30 b^{4}\right) }{ 3 \left(R(b^{5}-a^{5}) +15 b^{4}\right) } \boldsymbol{e}_{z}. \end{equation}

The cases in which the fluid is restricted to either leak normally or tangentially through the container results in two different container and squirmer velocities which implies that mechanisms for normal and tangential fluid leakage are not equivalent in the sense that only discontinuities in the normal stress may be driven by a finite pressure jump. There is no equivalent mechanism for driving fluid tangentially across the membrane. These two limiting cases correspond to where either $R_{\perp }$ or $R_{\parallel }$ diverges. In these cases, container and squirmer translational velocities become

(3.32)\begin{gather} \left. \begin{gathered} \lim_{R_{{\perp}} \to \infty}\boldsymbol{U}^{sq} = \frac{B^{sq}_{1} \left(3 a^5 R_{{\parallel}}-5 a^3 b^2 R_{{\parallel}}+2 b^5 R_{{\parallel}}+10 b^4\right) }{ 3 \left(R_{{\parallel}} (b^5-a^5) +5 b^4\right) } \boldsymbol{e}_{z},\\ \lim_{R_{{\perp}}\to \infty}\boldsymbol{U}^{m}= \frac{10 a^3 b B^{sq}_{1}}{3 \left(R_{{\parallel}} (b^5-a^5) +5 b^4\right)}\boldsymbol{e}_{z}, \end{gathered} \right\} \end{gather}

(3.33)\begin{gather} \left. \begin{gathered} \lim_{R_{{\parallel}} \to \infty}\boldsymbol{U}^{sq} = \frac{B^{sq}_{1} \left(3 a^5 R_{{\perp}}-5 a^3 b^2 R_{{\perp}}+2 b^5 R_{{\perp}}+20 b^4\right)}{3R_{{\perp}}(b^5-a^5) +30 b^4}\boldsymbol{e}_{z},\\ \lim_{R_{{\parallel}}\to \infty}\boldsymbol{U}^{m}={-}\frac{10 a^3 b B^{sq}_{1}}{3R_{{\perp}} (b^5-a^5) +30 b^4}\boldsymbol{e}_{z}. \end{gathered} \right\} \end{gather}

The limit in which both $R_{\perp } \to \infty$ and $R_{\parallel } \to \infty$ recovers the solution for the rigid non-porous container, (3.18). The case in which $\{R_{\perp },R_{\parallel }\} \to \{\infty , 0\}$ recovers the solution for a viscous drop, in which there is no viscosity discontinuity, yet a zero shear stress jump at the fluid–membrane interface. For a viscous drop, $\boldsymbol {U}^{sq}$ becomes (3.18) and $\boldsymbol {U}^{m} = 2/3B^{sq}_{1} a^{3}/b^{3}\boldsymbol {e}_{z}$. This solution agrees with the more general viscous drop solution presented by Reigh et al. (Reference Reigh, Zhu, Gallaire and Lauga2017). In the limits that $R_{\perp } \to 0$ or $R_{\parallel } \to 0$ the squirmer velocity goes to the free space solution, $\boldsymbol {U}^{sq} \to 2/3 B^{sq}_{1}\boldsymbol {e}_{z}$, yet the container moves at a reduced speed given by

(3.34a,b)\begin{gather} \lim_{R_{{\perp}} \to 0 }\boldsymbol{U}^{sq} = \frac{2}{3}B^{sq}_{1}\boldsymbol{e}_{z},\quad \lim_{R_{{\perp}} \to 0 }\boldsymbol{U}^{m} ={-}\frac{a^3 B^{sq}_{1}}{3 b^3}\boldsymbol{e}_{z}, \end{gather}

(3.35a,b)\begin{gather} \lim_{R_{{\parallel}} \to 0 }\boldsymbol{U}^{sq} = \frac{2}{3}B^{sq}_{1}\boldsymbol{e}_{z},\quad \lim_{R_{{\parallel}} \to 0 }\boldsymbol{U}^{m} = \frac{2 a^3 B^{sq}_{1}}{3 b^3}\boldsymbol{e}_{z}, \end{gather}

owing again to the difference in whether a normal or shear stress discontinuity drives fluid across the membrane.

3.4. Velocity fields

Velocity fields with streamlines have been constructed over the set of squirmer types $\beta \in \{-5,0,+5\}$ and resistances $R_{m}:= \{R_{\parallel },R_{\perp }\}$, for a fixed container to squirmer size ratio $b/a = 5$. Velocity fields for a $\beta = -5$ pusher, $\beta = 5$ puller, and $\beta = 0$ neutral squirmer are shown respectively in panel pairs 3(a,d), 3(b,e) and 3(c,f). The aforementioned pairs of panels are grouped by resistance parameterizations $R_{\parallel } < R_{\perp }$ for $R/a = \{10,100\}$ and $R_{\parallel } > R_{\perp }$ for $R/a = \{100,10\}$. The velocity fields are constructed in Cartesian coordinates by sweeping over the sets $\{\psi ^{i}(x,z,\phi ): b^{2} \ge x^{2}+z^{2} \ge a^{2}, \phi \in [0,{\rm \pi} ] \}$ and $\{\psi ^{e}(x,z,\phi ): x^{2}+z^{2}>b^{2}, \phi \in [0,{\rm \pi} ] \}$ and applying (3.2). Streamlines are constructed using standard formalisms.

Figure 3.

Velocity fields are shown for a container to particle size ratio $b/a = 5$, squirmer types $\beta \in \{-5,0,5\}$ and membrane resistance parameters $R_{m}=\{R_{\parallel },R_{\perp }\} = \{10,100\}$ and $\{100,10\}$. The centre arrow shows the squirmer's swim direction, $\boldsymbol {e}_{s} := \boldsymbol {e}_{z}$. The fluid velocity scales have normalized units $|\boldsymbol {u}(\boldsymbol {x},\boldsymbol {z})|/|\boldsymbol {U}^{sq}_{fs}|$.

Velocity fields in which $R_{\parallel } < R_{\perp }$ show that the fluid exits preferentially in tangential directions relative to the membrane surface. Conversely the condition, $R_{\parallel } > R_{\perp }$ forces fluid to enter and exit the membrane with a larger normal velocity component. Although translational velocities of the squirmer and container, in the concentric geometry, are independent of the squirmer's type, $\beta$, the flow fields are very much dependent on $\beta$. Neutral squirmers are seen to produce flow fields that are similar to the source dipole flow field. Pushers tend to draw fluid in radially (relative to $\boldsymbol {e}_{s}$) and expel it out axially. Pullers draw fluid in axially and expel it out radially. Both pushers and pullers expel fluid asymmetrically relative to the $xy$-plane, and thus generate net motion along their orientation $\boldsymbol {e}_{s}$ (in this case in $+\boldsymbol {e}_{z}$). Pushers and pullers should functionally behave like the point force dipole solution of Stokes flow; however, the container induces vortical flows in both the $z$-anterior and $z$-posterior portions of the container. This phenomenon is purely of hydrodynamic origin, owing to the confinement effects of the container.

4.1. Squirmer velocity contribution

Since the Stokes double-layer potential is able to fully represent flow fields that correspond to force- and torque-free surfaces, the double-layer potential can be used to elegantly describe the squirmer's disturbance velocity field. The squirmer's velocity contribution is formulated based on (4.5) as an indirect Fredholm integral equation of the second kind in an aphysical velocity density acting on the Stokes double-layer bivariate potential. Observing the form of (4.5) and the aforementioned jump properties of the double layer potential in (4.3), the squirmer's disturbance velocity field in (4.7) is found to be

(4.9)\begin{align} \boldsymbol{U}^{sq} + \boldsymbol{\varOmega}^{sq} \times ( \boldsymbol{x} - \boldsymbol{x}^{sq}_{c}) + \boldsymbol{u}^{sl}(\boldsymbol{x}) - \bar{\boldsymbol{u}}^{sq}(\boldsymbol{x}) = \boldsymbol{\phi}(\boldsymbol{x})\nonumber\\ + \int_{\varGamma_{sq}}\boldsymbol{\mathsf{K}}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\phi}(\kern1.5pt\boldsymbol{y})\,\textrm{d}S(\kern1.5pt\boldsymbol{y}), \quad \boldsymbol{x} \in \varGamma_{sq}, \end{align}

where the velocity boundary condition at the squirmer's surface is used in conjunction with (4.7) to couple the fluid velocity on the squirmer's surface to the squirmer's rigid body motion. The slip velocity $\boldsymbol {u}^{sl}(\boldsymbol {x})$, defined in (2.2), is solely responsible for generating the propulsive force, $\boldsymbol {F}^{P}$, that drives the squirmer whereas the rigid body motion terms, $\boldsymbol {U}^{sq}$ and $\boldsymbol {\varOmega }^{sq}$, give rise to the drag force $\boldsymbol {F}^{D}$ from (3.15a). Although the singularity present in the kernel $\boldsymbol{\mathsf{K}}$ is proportional to $1/r^{2}$, (4.9) is a second kind integral equation and admits rigorous analysis under the Fredholm–Riesz–Schauder theory (Hsiao & Wendland Reference Hsiao and Wendland2008; Kress Reference Kress2014) so long as the surface is sufficiently smooth.

Using the CDLBIEM completion relationships (Karrila & Kim Reference Karrila and Kim1989; Kim & Karrila Reference Kim and Karrila2005)

(4.10)\begin{gather} \boldsymbol{U}^{sq} ={-} \sum_{j=1}^{3}\boldsymbol{\varphi}^{j,RBM}(\boldsymbol{x})\langle\boldsymbol{\varphi}^{j,RBM}, \boldsymbol{\phi} \rangle, \end{gather}

(4.11)\begin{gather}\boldsymbol{\varOmega}^{sq} \times( \boldsymbol{x} - \boldsymbol{x}_{c}^{sq}) ={-} \sum_{j=4}^{6}\boldsymbol{\varphi}^{j,RBM}(\boldsymbol{x})\langle \boldsymbol{\varphi}^{j,RBM},\boldsymbol{\phi} \rangle, \end{gather}

where $\boldsymbol {x}_{c}^{sq}$ is the squirmer's centroid, (4.9) may be written solely in terms of the velocity potential $\boldsymbol {\phi }$. The, $\boldsymbol {\varphi }^{i,RBM}, i \in \{1,\ldots ,6\}$, are the rigid body motion (RBM) eigenfunctions of the homogeneous form of (4.3) with eigenvalue –1, which are proportional to $\boldsymbol {e}_{x}, \boldsymbol {e}_{y}, \boldsymbol {e}_{z}$, $\boldsymbol {e}_{x} \times (\boldsymbol {x} - \boldsymbol {x}_{c}^{sq})$, $\boldsymbol {e}_{y} \times (\boldsymbol {x} - \boldsymbol {x}_{c}^{sq})$ and $\boldsymbol {e}_{z} \times (\boldsymbol {x} - \boldsymbol {x}_{c}^{sq})$. In practice, when dealing with a mesh approximation to $\varGamma _{sq}$, these functions are numerically determined by making them orthonormal with respect to the natural inner product norm on the surface

(4.12)\begin{equation} \langle \boldsymbol{a}, \boldsymbol{b} \rangle = \int_{\varGamma_{sq}}a_{i}(\boldsymbol{x})b_{i}(\boldsymbol{x})\,\textrm{d}S(\boldsymbol{x}), \end{equation}

using the modified Gram–Schmidt process. For a perfectly spherical surface they take the forms (Phan-Thien & Kim Reference Phan-Thien and Kim1994)

(4.13a)\begin{gather} \boldsymbol{\varphi}^{i,RBM}(\boldsymbol{x}) = \frac{1}{\sqrt{|\varGamma_{sq}|}}\boldsymbol{e}_{i}, \end{gather}

(4.13b)\begin{gather}\boldsymbol{\varphi}^{i+3,RBM}(\boldsymbol{x}) = \frac{1}{\sqrt{I_{i}}}\boldsymbol{e}_{i} \times (\boldsymbol{x} - \boldsymbol{x}^{sq}_{c}), \end{gather}

where $|\varGamma _{sq}|$ is the surface area of the squirmer and

(4.14)\begin{equation} I_{i} = \int_{\varGamma_{sq}}\left[ (\kern1.5pt\boldsymbol{y}-\boldsymbol{x}^{sq}_{c})\boldsymbol{\cdot} (\kern1.5pt\boldsymbol{y}-\boldsymbol{x}^{sq}_{c}) - (\kern1.5pt\boldsymbol{y}-\boldsymbol{x}^{sq}_{c})_{i}^{2} \right] \,\textrm{d}S(\kern1.5pt\boldsymbol{y} ), \quad i \in \{1,\ldots,3\}, \end{equation}

are the surface moments.

Substituting (4.10) and (4.11) into (4.9) results in

(4.15)\begin{equation} \boldsymbol{u}^{sl}(\boldsymbol{x}) - \bar{\boldsymbol{u}}^{sq}(\boldsymbol{x}) = \boldsymbol{\phi}(\boldsymbol{x}) + \int_{\varGamma_{sq}}\boldsymbol{\mathsf{K}}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\phi}(\kern1.5pt\boldsymbol{y})\,\textrm{d}S(\kern1.5pt\boldsymbol{y})\nonumber\\ + \sum_{j=1}^{6}\boldsymbol{\varphi}^{j,RBM}(\boldsymbol{x})\langle \boldsymbol{\varphi}^{j,RBM}, \boldsymbol{\phi} \rangle, \end{equation}

where the integral equation's eigenvalue at $-1$ of rank six, corresponding to the operator $(\mathbb {I} + \mathcal {K})$, is mapped to zero. After solving for $\boldsymbol {\phi }$, post processing with the completion relations, (4.10) and (4.11), allows for recovery of $\boldsymbol {U}^{sq}$ and $\boldsymbol {\varOmega }^{sq}$.

4.2. Container velocity contribution

In formulating the container contribution, explicit indications of where $\boldsymbol {x}$ and $\boldsymbol {y}$ are located are made relative to the defined container normal vector, $\boldsymbol {n}$. Since the container normal points exterior to the container, quantities that depend on points $\boldsymbol {x}$ and $\boldsymbol {y}$ in $\varOmega _{e}$ and in the direction of $\boldsymbol {n}$, are denoted with a superscript $+$. Similarly, quantities that depend on $\boldsymbol {x}$ and $\boldsymbol {y}$ in volume $\varOmega _{i}$ bounded by $\varGamma _{c-}$, are denoted with a superscript $-$.

Under zero viscosity contrast inside and outside of the container, the Stokes boundary integral equations in the primary variables (surface velocities $\boldsymbol {u}(\boldsymbol {x})$ and tractions $\boldsymbol {f}(\boldsymbol {x})$) may be applied to both sides of the container to find that

(4.16) \begin{align} \boldsymbol{u}^{c}(\boldsymbol{x}) - \bar{\boldsymbol{u}}^{c}(\boldsymbol{x}) &= \left[ \boldsymbol{u}^{+}(\boldsymbol{x}) - \boldsymbol{u}^{-}(\boldsymbol{x}) \right] + \int_{\varGamma_{c}}\boldsymbol{\mathsf{K}}(\boldsymbol{x},\boldsymbol{y}) \left[ \boldsymbol{u}^{+}(\kern1.5pt\boldsymbol{y}) - \boldsymbol{u}^{-}(\kern1.5pt\boldsymbol{y})\right]\,\textrm{d}S(\kern1.5pt\boldsymbol{y})\nonumber\\ &\quad - \tfrac{1}{2}\int_{\varGamma_{c}}\boldsymbol{\mathsf{G}}(\boldsymbol{x},\boldsymbol{y}) \left[ \boldsymbol{f}^{+}(\kern1.5pt\boldsymbol{y}) - \boldsymbol{f}^{-}(\kern1.5pt\boldsymbol{y}) \right]\,\textrm{d}S(\kern1.5pt\boldsymbol{y}), \end{align}

where $\boldsymbol {u}^{c}(\boldsymbol {x})$ is the fluid velocity evaluated on the container/membrane surface. This fluid velocity is different from the membrane's physical point-wise velocity, $\boldsymbol {u}^{m}(\boldsymbol {x})$, unless no slip is used in conjunction with the standard kinematic boundary condition. The fluid velocity is assumed to be continuous across the membrane, as would be the case for pressure driven flow through a cylindrical pore (i.e. fully developed Hagen–Poiseuille flow). Imposing this condition as $\boldsymbol {u}^{+}(\boldsymbol {x}) = \boldsymbol {u}^{-}(\boldsymbol {x})$, (4.16) reduces to

(4.17)\begin{equation} \boldsymbol{u}^{c}(\boldsymbol{x}) - \bar{\boldsymbol{u}}^{c}(\boldsymbol{x}) ={-}\tfrac{1}{2}\int_{\varGamma_{c}}\boldsymbol{\mathsf{G}}(\boldsymbol{x},\boldsymbol{y}) \left[ \boldsymbol{f}^{+}(\kern1.5pt\boldsymbol{y}) - \boldsymbol{f}^{-}(\kern1.5pt\boldsymbol{y}) \right]\,\textrm{d}S(\kern1.5pt\boldsymbol{y}). \end{equation}

To proceed with modelling the squirmer–container system, a constitutive law, analogous to (3.23) and (3.20), is needed to describe the stress jump across the membrane. The jump in stress is modelled as

(4.18)\begin{align} \Delta \boldsymbol{f}(\kern1.5pt\boldsymbol{y}) &:= \left[ \boldsymbol{f}^{+}(\kern1.5pt\boldsymbol{y}) - \boldsymbol{f}^{-}(\kern1.5pt\boldsymbol{y}) \right],\nonumber\\ &={-}\left[ \boldsymbol{f}^{b}(\kern1.5pt\boldsymbol{y}) + \boldsymbol{f}^{\gamma}(\kern1.5pt\boldsymbol{y}) + \boldsymbol{f}^{p}(\kern1.5pt\boldsymbol{y})\right], \end{align}

which includes forces related to membrane bending $\boldsymbol {f}^{b}$, tension $\boldsymbol {f}^\gamma$ and resistance to permeable flow $\boldsymbol {f}^{p}$. In the present model the membrane is rigid and can only support a permeable force density expressed as

(4.19)\begin{align} \boldsymbol{f}^{p}(\kern1.5pt\boldsymbol{y}) &={-}\left[ \boldsymbol{u}^{c}(\kern1.5pt\boldsymbol{y})-\boldsymbol{u}^{m}(\kern1.5pt\boldsymbol{y}) \right ] \boldsymbol{\cdot} \left[ R_{n}\boldsymbol{n}(\kern1.5pt\boldsymbol{y})\boldsymbol{n}(\kern1.5pt\boldsymbol{y}) + R_{b}\boldsymbol{b}(\kern1.5pt\boldsymbol{y})\boldsymbol{b}(\kern1.5pt\boldsymbol{y}) + R_{t}\boldsymbol{t}(\kern1.5pt\boldsymbol{y})\boldsymbol{t}(\kern1.5pt\boldsymbol{y}) \right],\nonumber\\ &={-}\left\{ \boldsymbol{u}^{c}(\kern1.5pt\boldsymbol{y})-\left[ \boldsymbol{U}^{m} + \boldsymbol{\varOmega}^{m} \times (\kern1.5pt\boldsymbol{y} - \boldsymbol{x}^{c}_{c})\right] \right \}\nonumber\\ &\quad \boldsymbol{\cdot} \left\{ R_{n}\boldsymbol{n}(\kern1.5pt\boldsymbol{y})\boldsymbol{n}(\kern1.5pt\boldsymbol{y}) + R_{b}\boldsymbol{b}(\kern1.5pt\boldsymbol{y})\boldsymbol{b}(\kern1.5pt\boldsymbol{y}) + R_{t}\boldsymbol{t}(\kern1.5pt\boldsymbol{y})\boldsymbol{t}(\kern1.5pt\boldsymbol{y}) \right\}. \end{align}

This model is equivalent to the analytical formulation under the identifications $R_{\perp } := R_{n}$, $R_{\parallel } := R_{b} = R_{t}$. The vectors $\boldsymbol {t},\boldsymbol {b},\boldsymbol {n}$ are respectively the tangent, bitangent and normal vectors and provide a localized point-wise tangent space on the container's surface.

In addition, integral constraints on the total hydrodynamic force and torque on the membrane container must be made to fix the system

(4.20)\begin{gather} \int_{\varGamma_{c}}\Delta \boldsymbol{f}(\kern1.5pt\boldsymbol{y})\,\textrm{d}S(\kern1.5pt\boldsymbol{y}) = \boldsymbol{F}^{H} = \boldsymbol{0}, \end{gather}

(4.21)\begin{gather}\int_{\varGamma_{c}} (\kern1.5pt\boldsymbol{y} - \boldsymbol{x}_{c}^{c}) \times \Delta \boldsymbol{f}(\kern1.5pt\boldsymbol{y})\,\textrm{d}S(\kern1.5pt\boldsymbol{y}) = \boldsymbol{T}^{H} = \boldsymbol{0}, \end{gather}

where $\boldsymbol {x}_{c}^{c}$ is the container's centroid, and $\boldsymbol {F}^{H}$ and $\boldsymbol {T}^{H}$ are the hydrodynamic force and torque acting on the container.