Swimming Mode of Two Interacting Squirmers under Gravity in a Narrow Vertical Channel

Abstract

The swimming mode of two interacting squirmers under gravity in a narrow vertical channel is simulated numerically using the lattice Boltzmann method (LBM) in the range of self-propelling strength 0.1 ≤ α ≤ 1.1 and swimming type −5 ≤ β ≤ 5. The results showed that there exist five typical swimming patterns for individual squirmers, i.e., steady upward rising (SUR), oscillation across the channel (OAC), oscillation near the wall (ONW), steady upward rising with small-amplitude oscillation (SURO), and vertical motion along the sidewall (VMS). The parametric space (α, β) illustrated the interactions on each pattern. In particular, the range of oscillation angle for ONW is from 19.8° to 32.4° as α varies from 0.3 to 0.7. Moreover, the swimming modes of two interacting squirmers combine the two squirmers’ independent swimming patterns. On the other hand, the pullers (β < 0) attract with each other at the initial stage, resulting in a low-pressure region between them and making the two pullers gradually move closer and finally make contact, while the result for the pushers (β > 0) is the opposite. After the squirmers’ interaction, the squirmer orientation and pressure distribution determine subsequent squirmer swimming patterns. Two pushers separate quickly, while there will be a more extended interaction period before the two pullers are entirely separated.

1. Introduction

In recent years, researchers have attracted considerable interest from self-propelled microorganisms because of their wide range of applications in nature and industrial processes. Microorganisms in nature swim under their energy in different forms [1,2,3,4], and the swimming Reynolds number (Re) of microorganisms is usually minimal because the scale of microorganisms is tiny [5]. The squirmer has been studied as a typical driving model of microorganisms and self-propelled particles [6,7]. Based on their driving model, the squirmer could be classified into the pusher, puller, and neutral squirmer [8]. The pusher and puller have different swimming characteristics, e.g., the swimming speed of the pusher increased significantly with increasing Re, but the swimming speed of the puller was not related to Re [9].

The behavior of squirmers swimming near walls is also of great interest. The swimming pattern of squirmers depended on the combined action of self-propulsion, gravity, and hydrodynamic interactions with walls [10]. Several distinct swimming behaviors exist, such as touching the wall and then moving away from the wall, swimming along the wall with periodic oscillations, and swimming parallel to the wall [11]. Fadda et al. [12] studied the dynamics of a chiral swimmer sedimenting on a flat plate. They revealed that different states of gravitational intensity produced different dynamics, but only the wall perpendicular to the direction of gravity affected the swimmer’s dynamics. Recently, Ouyang and Lin [13] investigated the behavior of a settling micro-swimmer in a narrow vertical channel and found four typical locomotive patterns, i.e., the vertical motion, attracted oscillatory, oscillatory, and upward motion. The patterns were related to Re and the density ratio of a swimmer to fluid. Moreover, the swimming pattern near the boundary of axially symmetric squirmers in an inertia-less Newtonian fluid with a no-slip interface was theoretically proposed by Ishimoto and Gaffney [14], indicating that the tangential squirmer is a relatively simple framework that enables predictions and classifications of the complexities associated with axially symmetric squirmers swimming near the boundary.

The study of microbial interactions provides a better understanding of microbial motility behavior. Several investigations had empirically indicated that the hydrodynamic interactions between microorganisms played a vital role in defining their dynamics, such as the dancing alga Volvox [15], changes in direction between two swimming [16], dynamic clustering in suspensions of motile bacteria [17], and a new form of collective swimming bacteria dynamics [18]. Saintillan and Shelley [19] studied the orientational order and dynamics numerically in suspensions of self-locomoting slender rods, considering the far-field interactions. They found that nematic suspensions of swimming particles were unstable at long wavelengths. Lushi et al. [20] used a spiral vortex state in confined suspensions of Bacillis subtilis to study the dynamics of swimmer suspensions. For near-field hydrodynamics between two interacting swimming microorganisms, Papavassiliou and Alexander [21] provided exact solutions of the Stokes equations for squirmers near a no-slip surface, which allowed the hydrodynamic interactions of swimming microscopic organisms to be determined in the near-field. Li et al. [11] numerically investigated the hydrodynamic interactions of swimming organisms at small to intermediate Re. They found that the inertial effects changed the contact time and dispersion dynamics of a pair of pusher swimmers and triggered the hydrodynamic attraction of two pullers. Their findings had significant potential implications for studying predator–prey interactions. Therefore, it is essential to understand the effect of inertia on the hydrodynamic interaction between two swimmers. Recently, More and Ardekani [22] studied numerically the interactions between a pair of model swimming organisms in configurations of approaching each other and moving side by side with finite inertia in a linearly density-stratified fluid and observed four types of interactions. Qi et al. [23] studied the hydrodynamic properties of a squirmer type of self-propelled particle in a simple shear flow. They found that increasing the swimming intensity of the squirmer would induce the squirmer to make a periodic and stable motion at a specific distance from the wall. However, no reports on the interaction between the two rising squirmers have been published. Not only do microorganisms have common sedimentation motions, but rising motions are also indispensable because of their physical characteristics [24] and survival needs [25].

Despite the significant progress made in the past, the study of swimming microorganisms in the presence of external forces is still insufficient, especially under the influence of gravity. A suspension of photosynthetic algae, denser than water, may have average swimming upward in still water. Based on this phenomenon, Hilla et al. [26] found that when the density of the upper regions of the suspension becomes more significant than the lower region it may lead to biological connections. Pedley [27] also reported that the upward swimming behavior of algae is influenced by gravity. In other words, the bottoms of algae are heavy, so they are affected by gravitational torque when they are titled. In addition, gravity can also determine the swimming behavior of some microorganisms, which are usually immersed in water and can sense gravity. For example, the weighted leverage within Loxodes can control the swimming direction [28], and Euglena gracilis (a photosynthetic flagellate) uses gravity to adjust its vertical position in the water to avoid excessive solar radiation [29]. Studying how microorganisms respond to external gravitational fields and how they interact under gravity can provide qualitative assistance in designing artificially intelligent swimming robots. With this in mind, our objective is to comprehensively analyze the motion behavior and interactions of two rising squirmers under the influence of gravity. The swimming patterns of two interacting squirmers under gravity in a narrow vertical channel are simulated numerically using the LBM for different self-propelling strengths and swimming types to understand the collective swimming behavior of microorganisms better. Unlike before, our work focuses on two light squirmers. We hope our results will provide new insights into the dynamical characteristics of squirmers.

1.1. Lattice Boltzmann Model

In this work, we use the single relaxation time LBM [30,31,32] to deal with the particle motion in the flow because of its streamlined programming and sufficient accuracy. The discrete lattice Boltzmann equation is expressed as follows:

(1)$f_i(\mathit{x}+{\mathit{e}}_i\Delta t,t+\Delta t)-f_i(\mathit{x},t)=-\frac{1}{\tau }[f_i(\mathit{x},t)-f_i{}^{(0)}(\mathit{x},t)],$

where f_i(x, t) is the distribution function for the microscopic velocity; e_i is in the ith direction; f_i⁽⁰⁾(x, t) is the equilibrium distribution function; Δt is the time step of the simulation, and τ is the relaxation time related to the fluid viscosity ν.

The density ρ and velocity u of the fluid are computed using the following formulations:

(2)$\rho ={\displaystyle \sum _i{f}_i},\rho \mathit{u}={\displaystyle \sum _i{f}_i{\mathit{e}}_i},$

The equilibrium distribution function, f_i⁽⁰⁾(x, t), is defined as follows:

(3)$f_i{}^{(0)}(\mathit{x},t)={\omega }_i\rho [1+\frac{3{\mathit{e}}_i\cdot \mathit{u}}{c^2}+\frac{9{({\mathit{e}}_i\cdot \mathit{u})}^2}{2c^4}-\frac{3{\mathit{u}}^2}{2c^2}].$

where c = Δx/Δt and Δx is the lattice spacing; for simplicity, Δx and Δt are both fixed at 1, which is common in using LBM. In addition, the weights are ω₀ = 4/9, ω_1–4 = 1/9, and ω_5–8 = 1/36.

For two-dimensional (2-D) simulations, the D2Q9 lattice model with nine velocities proposed by Qian et al. [30] is used. Note that the speed of sound has the relation c_s = c/3^1/2, and the discrete velocity vectors are given:

(4)${\mathit{e}}_i=\{\begin{array}{l}(0, 0) \mathrm{for} i=0 \\ (\pm 1, 0)c,(0, \pm 1)c \mathrm{for} i=1 \mathrm{to} 4 \\ (\pm 1, \pm 1)c \mathrm{for} i=5 \mathrm{to} 8\end{array} .$

By performing a Chapman–Enskog expansion, the macroscopic mass and momentum equations in the low Mach number limit can be recovered,

(5)$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathit{u})=0,$

(6)$\frac{\partial (\rho \mathit{u})}{\partial t}+\nabla \cdot (\rho \mathit{u}\mathit{u})=-\nabla p+\rho \nu {\nabla }^2\mathit{u},$

The kinematic viscosity of the fluid is determined using the equation ν = c_s²(τ − 0.5)Δt.

1.2. Squirmer Model

To enable the self-propelling behavior of particles, we used the squirmer model, which Blake [7] obtained and modified based on the results from Lighthill [6] to simulate the motion of self-propelled particles. According to Blake [7], the slip velocity at one point on the squirmer’s surface consists of radial and tangential surface velocity components:

(7)${\mathit{u}}^s(\mathit{r})={\displaystyle \sum _{n=0}^{\infty }{A}_n}\cos (n\theta )\widehat{\mathit{r}}+{\displaystyle \sum _{n=1}^{\infty }{B}_n}\sin (n\theta )\widehat{\theta },$

where $\widehat{\mathit{r}}$ and $\widehat{\theta }$ are the radial and tangential unit vectors at one point on the squirmer surface. For a more visual illustration, the unit vectors are shown in Figure 1.

For simplicity, the radial component was usually ignored, according to Blake [7]. Under the same assumptions, we adopted the simplified squirmer model represented by the second-order truncated tangential velocity,

(8)${\mathit{u}}^s(\theta )=(B_1\sin \theta +B_2\sin \theta \cos \theta )\widehat{\theta }.$

where B₁ represents the magnitude of the self-propelled force, which is reflected in the final swimming velocity of the squirmer; for example, in the Stokes flow limit, the velocity of a squirmer is Us = B₁/2. Note that in all cases, it is necessary to ensure B₁ > 0. Coefficient B₂ it determines the intensity of the vorticity around the squirmer. Here, a swimming type parameter β = B₂/B₁ is introduced; the squirmer can be divided into three different types: pullers (β > 0), pushers (β < 0), and neutrals (β = 0), based on the value of β. As the name implies, pullers generate thrust from their front (Figure 2a). The pushers push themselves using the cilia at the back of their body, thus generating thrust from the rear (Figure 2b).

1.3. Boundary Conditions

In LBM, special treatment of the moving boundary is usually required to ensure no-slip boundary conditions on the surfaces of squirmers. In this work, the improved bounce-back scheme proposed by Lallemand et al. [33] was used and, after improvement, this interpolation-based bounce-back scheme is stable, robust, and easy to implement in simulations.

The improved bounce-back scheme is shown in Figure 3, where solid squares indicate solid nodes inside the curved boundary, solid circles indicate fluid nodes outside the curved boundary, and hollow triangles indicate boundary nodes where the bounce-back occurs. Lallemand et al. [33] proposed that all fluid and solid nodes may change after the collision when the distribution function from the solid nodes to the fluid nodes needs to be recalculated based on the exact location of the boundary (i.e., boundary nodes). For example, the distribution function f₆ of the fluid node A originally needs to be calculated by the solid node B. However, node B is inside the curved boundary, so the distribution function f6 cannot be determined. A simple bounce-back scheme can be chosen to satisfy the no-slip boundary condition, such as making f₆ = f₈. However, this method assumes that the boundary node F is located precisely at the midpoint of the line between two lattice nodes, which is not universal. Therefore, an additional calculation of the distribution function f₆ is required. Lallemand et al. [33] introduced the parameter q, q = |AF|/|AB|, to obtain the location of the boundary nodes more accurately. Then, f₆ can be recalculated based on the near nodes using the interpolation method.

(9)$f_6(\mathrm{A})=\{\begin{array}{l}q(1+2q)f_8(\mathrm{B})+(1-4q^2)f_8(\mathrm{A})-q(1-2q)f_8(\mathrm{C})-2w_8\rho \frac{{\mathit{e}}_8\cdot {\mathit{u}}_{\mathrm{F}}}{c_s{}^2}\begin{array}{c}\end{array}\begin{array}{c}\end{array}q<\frac{1}{2}\hfill \\ \frac{1}{q(2q+1)}{f}_8(\mathrm{B})+\frac{2q-1}{q}{f}_6(\mathrm{C})-\frac{2q-1}{2q+1}{f}_6(\mathrm{D})-\frac{2w_8}{q(2q+1)}\frac{{\mathit{e}}_8\cdot {\mathit{u}}_{\mathrm{F}}}{c_s{}^2}\begin{array}{c}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\end{array}q\ge \frac{1}{2}\hfill \end{array}.$

where u_F is the moving velocity of the surface at the boundary node F, as shown in Figure 3.

The self-propelled motion of the squirrel can be achieved by updating the velocity of the boundary nodes by adding the slip velocity before interpolation, as shown in the following equation,

(10)${\mathit{u}}_{\mathrm{F}}={\mathit{u}}_b+{\mathit{u}}^s=\mathit{U}+\mathit{\Omega }\times \mathit{R}+{\mathit{u}}^{{}_s}.$

where U and Ω represent the translational and rotational speed of the squirmer and R is the position vector starting from the squirmer mass center to boundary node F.

During the calculation, all fluid-solid boundary links will be scanned at each time step, and the force and torque applied to the squirmer can be calculated by the momentum exchange method [34]. In addition, the increased force and torque when the squirmers swim through the fluid and cover the fluid nodes are also addressed [35]. In summary, the motion of the squirmer is defined by solving Newton’s equations using the net force and torque.

(11)$m\frac{d^2\mathit{x}}{d{t}^2}=\mathit{F}$

(12)$\frac{d(J\cdot \mathit{\Omega })}{dt}=\mathit{T}$

where x is the center of the squirmer’s mess, Ω and J represent the angular velocity and moment of inertia of squirmers, respectively. F and T represent the force and torque, respectively, when ρ_p (the density of squirmer) = ρ (the density of fluid), the squirmer was neutrally suspended.

1.4. Repulsive Force Model

The distance between squirmers or between the squirmer and wall surface in the simulation may be minimal. In order to prevent particles from penetrating each other or penetrating walls, the collision model proposed by Glowinski et al. [36] was used. B_i and B_j are assumed to be two circular squirmers with radii and centers of mass of R_i, R_j, and G_i, G_j, respectively. Then, the repulsive force, F_ij^p, between particles B_i and B_j must satisfy the following equation:

(13)$F_{ij}^p=\{\begin{array}{l}0\begin{array}{cc}& \end{array}\begin{array}{cc}& \end{array}\begin{array}{cc}& \end{array}\begin{array}{cc}& \end{array}\\ \frac{1}{{\epsilon }_p}{(\frac{R_i+R_j+h_c-d_{ij}}{h_c})}^2\frac{{\mathit{G}}_i{\mathit{G}}_j}{d_{ij}}\end{array}\begin{array}{c}\\ \end{array}\begin{array}{cc}& \mathrm{if}\\ & \mathrm{if}\end{array}\begin{array}{c}\\ \end{array}\begin{array}{c}\\ \end{array}\begin{array}{c}{d}_{ij}>R_i+R_j+h_c\\ d_{ij}\le R_i+R_j+h_c\end{array},$

where d_ij = |G_iG_j|, h_c is the range of the repulsive force (i.e., the repulsive force occurs at less than h_c), and ε_p is a “small” positive number and is set to 10⁻³. A similar approach is used to deal with the particle-wall interaction by introducing a virtual particle (G’_ij), as shown in Figure 4, where d_ij’ = |G_iG’_ij|. ε_w is similar to ε_p, which is used for the squirmer-wall collision, and the squirmer-wall repulsive force F_ij^w is given by,

(14)$F_{ij}^w=\{\begin{array}{l}0\begin{array}{cc}& \end{array}\begin{array}{cc}& \end{array}\begin{array}{cc}& \end{array}\begin{array}{cc}& \end{array}\\ \frac{1}{{\epsilon }_w}{(\frac{R_i+R_j+h_c-{d^{\prime }}_{ij}}{h_c})}^2\frac{{\mathit{G}}_i{\mathit{G}}_j{}^{\prime }}{d_{ij}{}^{\prime }}\end{array}\begin{array}{c}\\ \end{array}\begin{array}{cc}& \mathrm{if}\\ & \mathrm{if}\end{array}\begin{array}{c}\\ \end{array}\begin{array}{c}\\ \end{array}\begin{array}{c}{d}_{ij}>R_i+R_j+h_c\\ d_{ij}\le R_i+R_j+h_c\end{array},$

2. Flow and Parameters

The swimming pattern of two interacting squirmers under gravity in a narrow vertical channel (Tow-dimensional), as shown in Figure 5, is simulated numerically. Two squirmers with the same diameter d and density ρ_s swim in the channel filled with a fluid with density ρ and kinematic ν. In the simulation, unless otherwise specified, the parameters are fixed as follows: d = 24, L= 5d, H =40d, initial distance between particles L₁ = 2d, H_u = 15d, and H_d = 25d. The squirmers swim upward (in the opposite direction of gravity) because of ρ_s < ρ. The parameters above are all based on lattice units common in the LBM. This study used a moving computational domain to simulate an infinite channel. The Reynolds number is defined as:

(15)$R{e}_g=\frac{U_{g}d}{\nu } ,$

(16)$U_g=\frac{\left|{\rho }_s/\rho -1\right|g{r}^2}{4\pi K\nu },$

where g is the gravitational acceleration, r is the squirmer radius, ν is the kinematic viscosity of the fluid, and K is a constant related to the effect of the channel width and is related to the ratio of the channel width to squirmer diameter (W^* = W/d):

(17)$K=\frac{1}{In{W}^{*}-0.9157+1.7244/{(W^{*})}^2-1.7302/{(W^{*})}^4+2.4056/{(W^{*})}^6-4.5913/{(W^{*})}^8}.$

As mentioned above, a squirmer’s final steady swimming velocity in a Stokes flow is B₁/2, i.e., U_s = B₁/2. Therefore, the swimming Reynolds number is defined as Re_s = U_sd₁/ν. Additionally, parameter α = U_s/U_g indicates the intensity of the squirmer settlement, which describes the competition between the self-propelled force and gravity. The gravity-driven motion of squirmers is dominated by self-propulsion for α ≥ 1, called the cruising regime. The motion of the squirmers is dominated by gravity for α ≤ 1, referred to as the strong-gravity regime [10]. Parameter β = B₂/B₁ is used to indicate the swimming type of the squirmers. Therefore, to observe different dynamic states at low Re, Re_g = 0.1, α = [0.1, 1.1], and β = [−5, 5] are selected as shown in Figure 6a.

3. Validation

To check the accuracy of the LBM, we performed simulations for different β at Re_g = 0.1. As shown in Figure 6a, a squirmer with a diameter d = 40 moves to the right in a square area of size L × L = 800 × 800. The densities of the squirmer and fluid are ρ_f = ρ_s = 1, the initial position of the squirmer is in the center of the region at l × l = 400 × 400, and Re_s = 0.0005. As shown in Figure 6b, the squirmer starts swimming from rest, and all velocities eventually reach the theoretical value of B₁/2, represented by the black dashed line. The numerical simulation results are in good agreement with the theoretical solutions.

In addition, to further verify the accuracy of our method in dealing with the interaction between particles. The classical Drifting, Kissing, and Tumbling (DKT) motion of two circular particles were reproduced and compared with the results of Feng et al. [37]. The physical parameters selected are the same as those of Feng et al. [37]. The vertical channel is 20d wide (X-direction) and 80d high (Y-direction). The channel is filled with fluid, and the density and kinematic of the fluid are ρ = 1 and ν = 0.001, respectively. The density and radius of particles are ρ = 1.01 and d = 50, respectively. The initial velocities of the particles and fluid are zero. At the beginning of the release, the first particle is located 0.01d from the center of the channel at a height of 72d from the bottom, and the second particle is at the center of the channel at 68d. Due to the gravity drive, these two particles start to settle in the Y-direction. The repulsive force is added when the gap between particles is less than or equal to 2d + ξ and when the gap between the particle and wall is less than 2d + ξ. ξ is the range in which repulsion occurs; here, ξ = 1.5Δx. Figure 7a,b shows each particle center’s transverse and longitudinal coordinates, respectively. The comparison shows that our results are in great agreement with the simulation results obtained by previous authors, which proves the correctness of our method.

4. Results and Discussion

4.1. Swimming Mode

When squirmers swim in the channel under different α and β, their movement modes are diverse, and the collective movement modes depend on the swimming patterns of individual squirmers. Figure 8 shows the typical swimming patterns of individual squirmers. There are five patterns as shown in Figure 8a–e: steady upward rising (SUR) (Pattern 1), oscillation across the channel (Pattern 2), oscillation near the wall (Pattern 3), vertical motion along the center line with small-amplitude oscillation (Pattern 4), and vertical motion with invariant and variable orientation near the wall (Pattern 5). As α increases, the swimming velocity of the squirmer dominates, causing the squirmers in patterns 3 and 5 to swim in opposite directions.

Figure 9 shows the swimming modes of a pair of squirmers with different α (0.3≤ α ≤ 1.1) and β(−5 ≤ β≤ 5). It is observed that Mode 3, composed of patterns 1 and 5, is located in the region of β > 0. As the value of α and β increases, pattern 5 in Mode 3 changes from vertical up to vertical down, and the triangles change from black-filled to unfilled. The effect of changes in α and β on Mode 4 is the same as that of Mode 3, but Mode 4 is located in the region of β < 0. The swimming mode of all neutral squirmers is Mode 2, which is located in the regions with low β (|β| ≤ 3) but high α(|α|≥ 0.7). Mode 1, composed of pattern 1, is located at α = 0.3 and β = 5, and the trajectory is shown in Figure 8a. It can be seen that the pullers initially repel each other and then attract. After a period of interaction, the final vertical motion occurs in the center of the channel. The streamlines in Figure 10b show that the pullers have reached a delicate balance.

As previously described, the patterns of OAC are observed for all types of squirmers. Mode 2 at α = 0.9 and β = −3 is shown in Figure 11. The trajectory (Figure 11a) shows that two pushers attract each other until they lean together, after which there is a period of interaction. Eventually, they swim in a pattern from the leftmost end of the channel to the rightmost, i.e., OAC. As shown in Figure 11b–e, when the squirmers swim toward the wall, the swimming direction of the squirmers always points to the wall. When the squirmers come into contact with the wall, the swimming direction turns to the other side wall until making contact with the wall. The main difference between squirmers and conventional particles is the self-propelled direction, which can also be called the head. The trajectories and streamlines of Mode 4 are shown in Figure 12. One pusher swims near the wall while the other swims oscillate across the channel. Therefore, they are subjected to different swimming resistance, eventually causing the latter to leave the former far behind. When the former squirmer is closest to and farthest from the wall, the streamlines are presented in Figure 12b,c, where the wall pulls the squirmer back when the pusher tries to swim away from the wall.

4.2. Reasons for Different Modes

The trajectories of the two squirmers at α = 0.9 and different β values are presented in Figure 13. The squirmers exhibit an utterly different swimming behavior in the initial swimming phase because of the different β. For pullers (β = 5), as shown in Figure 13a, two pullers appear to repel at the beginning of the rising movement and then approach each other until fitting together. When the pullers fit together, the dynamic equilibrium between the two pullers is broken as the pullers continue to swim, which causes a large deflection in the horizontal position and results in various subsequent swimming modes. For the pushers (β = −5), as shown in Figure 13b, two pushers attract and move closer to each other, and then there is a period of equilibrium. Its duration is longer for pushers than that for pullers. The equilibrium is broken, resulting in a large deflection in the horizontal position, leading to different subsequent swimming modes.

To better understand the differences in the initial motions of the two pullers and pushers and the reasons for the loss of equilibrium of the two squirmers, the streamlines and pressure distribution are shown in Figure 14 and Figure 15, respectively. For pullers, the two sides of the pullers are high-pressure regions because a puller pulls the fluid to the sides of its body. In contrast, the upper and lower ends are low-pressure regions, which also causes the two pullers to repel each other initially (Figure 14a). However, at this time, the angular velocity of puller 1 on the left is negative (counterclockwise direction is positive), as shown in Figure 16a, thus making puller 1 rotate in a clockwise direction. In contrast, puller 2 on the right rotates in a counterclockwise direction, and thus the pullers’ swimming directions gradually point toward each other (Figure 14b). At this time, the high-pressure region between the pullers becomes a low-pressure region, while the repulsive force between them also becomes an attractive force. When the two pullers are close together, a low-pressure region is formed between them (Figure 14c), and equilibrium is maintained. However, as the pullers continued to swim, puller 2 suddenly acquired a greater negative angular velocity (Figure 16a), thereby breaking the equilibrium.

The heads and tails of the pushers are located in high-pressure regions because a pusher pushes the fluid toward the head and tail of its body, whereas the sides of their bodies are low-pressure regions. The opposite is true for pullers. The low-pressure region causes two pushers to be attracted to each other initially. Unlike the pullers, the angular velocity of pusher 1 on the left is favorable at this time (Figure 15b), i.e., it rotates counterclockwise. In contrast, pusher 2 on the right rotates clockwise. Then, two pushers swim backward, but the low-pressure region between them continues to pull them toward each other, with an equilibrium similar to that of the pullers. This back-to-back dynamic equilibrium is more stable than the head-to-head equilibrium, resulting in a more prolonged equilibrium for the pushers, as shown in Figure 16. However, as the pushers continue to swim, pusher 2 will gain a more positive angular velocity because there is a slight fluid disturbance, thus breaking the equilibrium.

When the dynamic equilibrium between the two squirmers is broken, their interaction continues until they separate from each other and swim separately. Moreover, the squirmer orientation and pressure distribution after the squirmers’ interaction determine their subsequent swimming patterns. Here, the parameter d’ = 2d is defined, i.e., when the distance between the centers of the two squirmers is larger than d’, the interaction between the two squirmers ends. In Figure 17, the points for d’ = 2d are marked. For α = 0.3, 0.5, and 0.7 (Figure 17a), the swimming pattern of the pushers is pattern 3 (Figure 8c), and the swimming angles all change in a homeopathic manner before contacting the wall, indicating the feasibility of d’. The instantaneous pressure distributions and squirmer orientation for α = 0.3, 0.5, and 0.7 are shown in Figure 18a–c, where the distribution and graph are similar. It can be seen that there is a low-pressure region in the lower right of the swimming direction of the pushers, which is the cause of the pusher’s counterclockwise rotation. This low-pressure region is attracted to the wall, thus pulling the pushers towards the wall and eventually forming the swimming pattern 3 (Figure 8c). When d’ = 2d for α = 0.9, the swimming direction of the pusher is nearly horizontal, the head and tail of the pusher are located in the high-pressure regions, and the low-pressure region is located below the body (Figure 18d), which allows the pusher to remain stable until it comes in contact with the wall. Eventually, the pusher interacts with the wall and forms the swimming pattern 2. When α increases to 1.1, the pressure distribution around the pusher is similar to that for α= 0.3, 0.5, and 0.7, as shown in Figure 18e. However, at this point, the pusher has a greater swimming intensity, thus generating the same swimming pattern but in the opposite direction.

For the pullers, the situation is different. The low-pressure region at the head of the below puller is attracted to the low-pressure region of the tail of the above puller, as shown in Figure 14e when the dynamic equilibrium is broken, which is opposite to the situation for the pushers shown in Figure 15e. Therefore, for pullers, there will be a more extended interaction period before the two pullers are entirely separated. Taking α= 0.5 and β= 5 as an example, when the two pullers swim near the wall at d’ = 2d, the low-pressure region at the head of the above puller and the tail of the below puller are both attracted by the wall, while there is a high-pressure region between the two pullers as shown in Figure 19a. These different forms of interaction with the wall cause their subsequent different swimming patterns, as shown in Figure 19b, where the puller swims vertically upward near the wall, as shown in Figure 8f.

5. Conclusions

The interaction between a pair of squirmers subject to gravity and the swimming mode of squirmers in a narrow vertical channel (2-D) are investigated in the range of 0.1 ≤ α ≤ 1.1 and −5 ≤ β ≤ 5. The main conclusions are summarized below.

The swimming modes of the two interacting squirmers are the combination of the two squirmers’ independent swimming patterns. The main swimming patterns are different for different types of squirmers. Neutral squirmers exhibit the pattern of OAC. The main swimming pattern for pullers is a combination of vertical motion along the center and near the wall. The main swimming pattern for pushers is a combination of OAC and ONW.

A pair of pullers swimming upward side-by-side will repel each other because of the high-pressure region between them. The pressure distribution changes with the deflection of the squirmer orientations, and the high-pressure region becomes a low-pressure one, resulting in subsequent attraction and close contact between the squirmers. The initial low-pressure region between the pushers causes them to be attracted and come into contact.

Lastly, the swimming modes in the range between 0.1 ≤ α ≤ 1.1 and −5 ≤ β ≤ 5 were summarized, and is was found that the different swimming modes are related to the rotation angles and the swimming velocities (after separation) of the two squirmers. Comparing the pushers at α =0.3 and α = 0.5, when β = −3, it can be seen that α = 0.3 has a smaller rotation angle (in the horizontal direction) and thus has a more significant horizontal velocity; while α = 0.5 is the opposite, leading to a transition between oscillation across channel and oscillation near the wall. A similar outcome was observed at α =0.5 and α = 0.7.

Author Contributions

Conceptualization, J.L. and G.G.; methodology, G.G. and D.N.; software, G.G.; validation, D.N.; formal analysis, G.G. and J.L.; investigation, G.G. and J.L.; resources, G.G.; data curation, G.G.; writing—original draft preparation, G.G.; writing—review and editing, J.L.; visualization, G.G.; supervision, D.N.; project administration, J.L.; funding acquisition, J.L. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Schematic diagram of the unit vectors.

Figure 2. Flow generated by (a) a puller and (b) a pusher.

Figure 3. Schematic of the bounce-back scheme in the LBM.

Figure 4. Schematic of the virtual particle.

Figure 5. Schematic of two squirmers swimming in a narrow channel.

Figure 6. Schematic of the validation of LBM: (a) schematic of a squirmer swimming in a square area; (b) evolution of the squirmer swimming velocity with time.

Figure 7. The position coordinates of the centers of the two particles: (a) Transverse coordinates; (b) longitudinal coordinates, ■: Feng et al. (2004) [37]; ●: Feng et al. (2004) [37].

Figure 8. Typical swimming patterns for individual squirmers at different α and β purple: pullers, red: pushers, green: neutrals). (a) Pattern 1: steady upward rising (SUR); (b) Pattern 2: oscillation across the channel; (c) Pattern 3: oscillation near the wall; (d) Pattern 4: vertical motion along the center line with small-amplitude oscillation; (e) Pattern 5: vertical motion with invariant orientation near the wall; (f) vertical motion with variable orientation near the wall.

Figure 8. Typical swimming patterns for individual squirmers at different α and β purple: pullers, red: pushers, green: neutrals). (a) Pattern 1: steady upward rising (SUR); (b) Pattern 2: oscillation across the channel; (c) Pattern 3: oscillation near the wall; (d) Pattern 4: vertical motion along the center line with small-amplitude oscillation; (e) Pattern 5: vertical motion with invariant orientation near the wall; (f) vertical motion with variable orientation near the wall.

Figure 9. Distributions of swimming modes with different β and α. Mode 1: composed of Pattern 1; Mode 2: composed of Pattern 2; Mode 3: composed of Patterns 1 and 5; Mode 4: composed of Patterns 2 and 3; Mode 5: composed of Pattern 4; Mode 6: composed of Pattern 5; Mode 7: composed of Pattern 3. The filled and unfilled shapes in Mode 3 and Mode 4 represent the same swimming patterns with different directions.

Figure 10. Trajectory and streamlines (solid black lines with arrows) of a pair of pullers (α = 0.3 and β = 5): (a) trajectory; (b) streamlines at t* = 224.23.

Figure 11. Trajectories and streamlines of a pair of pushers (α = 0.9 and β = −3): (a) trajectories; (b) streamlines for Point 1 at t* = 43.05; (c) streamlines for Point 2 at t* = 46.40; (d) streamlines for Point 3 at t* = 49.60; (e) streamlines for Point 4 at t* = 52.94.

Figure 12. Trajectories and streamlines of a pair of pushers (α = 0.5 and β = −3). (a) trajectories; (b) streamlines for Point 1 at t* = 60.01; (c) streamlines for Point 2 at t* = 62.79.

Figure 13. Trajectories of pairs of squirmers: (a) pullers at α = 0.9 and β = 5; (b) pushers at α = 0.9 and β= −5.

Figure 14. Instantaneous distribution of pressure and streamlines for pullers at α = 0.9 and β = 5. (a) t* = 0.16; (b) t* = 1.56; (c) t* = 2.49; (d) t* = 3.58; (e) t* = 4.05.

Figure 15. Instantaneous distribution of pressure and streamlines for pushers at α = 0.9 and β = −5. (a) t* = 0.31; (b) t* = 0.78; (c) t* = 1.25; (d) t* = 11.74; (e) t* = 12.36.

Figure 16. Evolution of the squirmer angular velocities (a) pushers at α = 0.9 and β= −5; (b) pullers at α = 0.9 and β = 5 (asterisk represents the normalized quantity).

Figure 17. Evolution of orientation angles of pushers with β = −3 at different α: (a) α = 0.3–0.7; (b) α = 0.9–1.1. The breakdown of equilibrium will result in different distributions of two pushers in the vertical position. For simplicity, the pusher located at the latter is chosen as the object of study, and the initial θ* is 0.5.

Figure 18. Instantaneous pressure distributions and squirmer orientation for β = −3 at different α at distance d’ = 2d where d’ is the distance between the centers of the two squirmers, (a) α= 0.3; (b) α = 0.5; (c) α = 0.7; (d) α = 0.9; (e) α = 1.1.

Figure 19. Instantaneous pressure distribution and squirmer orientation for α = 0.5 and β = 5: (a) point at which d’ = 2d; (b) pressure distribution and streamlines for a single puller.

References

1. Schuster, S.C.; Khan, S. The bacterial flagellar motor. Annu. Rev. Biophys. Biomol. Struct.; 1994; 23, pp. 509-539. [DOI: https://dx.doi.org/10.1146/annurev.bb.23.060194.002453]

2. Brennen, C.; Winet, H. Fluid mechanics of propulsion by cilia and flagella. Annu. Rev. Fluid Mech.; 1977; 9, pp. 339-398. [DOI: https://dx.doi.org/10.1146/annurev.fl.09.010177.002011]

3. Ajdari, A.; Stone, H.A. A note on swimming using internally generated traveling waves. Phys. Fluids; 1999; 11, pp. 1275-1277. [DOI: https://dx.doi.org/10.1063/1.869991]

4. Paxton, W.F.; Kistler, K.C.; Olmeda, C.C.; Sen, A.; Angelo, S.K.S.; Cao, Y.; Mallouk, T.E.; Lammert, A.P.E.; Crespi, V.H. Catalytic nanomotors: Autonomous movement of striped nanorods. J. Am. Chem. Soc.; 2004; 126, pp. 13424-13431. [DOI: https://dx.doi.org/10.1021/ja047697z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15479099]

5. Lauga, E.; Powers, T.R. The hydrodynamics of swimming microorganisms. Rep. Prog. Phys.; 2009; 72, 096601. [DOI: https://dx.doi.org/10.1088/0034-4885/72/9/096601]

6. Lighthill, M.J. On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers. Commun. Pure Appl. Math.; 1952; 5, pp. 109-118. [DOI: https://dx.doi.org/10.1002/cpa.3160050201]

7. Blake, J.R. A spherical envelope approach to ciliary propulsion. J. Fluid Mech.; 1971; 46, pp. 199-208. [DOI: https://dx.doi.org/10.1017/S002211207100048X]

8. Chwang, A.T.; Wu, T.Y.T. Hydromechanics of low-Reynolds-number flow. Part 2. Singularity method for Stokes flows. J. Fluid Mech.; 1975; 67, pp. 787-815. [DOI: https://dx.doi.org/10.1017/S0022112075000614]

9. Chisholm, N.G.; Legendre, D.; Lauga, E. A squirmer across Reynolds numbers. J. Fluid Mech.; 2016; 796, pp. 233-256. [DOI: https://dx.doi.org/10.1017/jfm.2016.239]

10. Rühle, F.; Blaschke, J.; Kuhr, J.T. Gravity-induced dynamics of a squirmer microswimmer in wall proximity. New J. Phys.; 2018; 20, 025003. [DOI: https://dx.doi.org/10.1088/1367-2630/aa9ed3]

11. Li, G.; Ostace, A.; Ardekani, A.M. Hydrodynamic interaction of swimming organisms in an inertial regime. Phys. Rev. E; 2016; 94, 053104. [DOI: https://dx.doi.org/10.1103/PhysRevE.94.053104] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27967048]

12. Fadda, F.; Molina, J.J.; Yamamoto, R. Dynamics of a chiral swimmer sedimenting on a flat plate. Phys. Rev. E; 2020; 101, 052608. [DOI: https://dx.doi.org/10.1103/PhysRevE.101.052608] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32575256]

13. Ouyang, Z.; Lin, J.Z. Behaviors of a settling microswimmer in a narrow vertical channel. Powder Technol.; 2022; 398, 117042. [DOI: https://dx.doi.org/10.1016/j.powtec.2021.117042]

14. Ishimoto, K.; Gaffney, E.A. Squirmer dynamics near a boundary. Phys. Rev. E; 2013; 88, 062702. [DOI: https://dx.doi.org/10.1103/PhysRevE.88.062702] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24483481]

15. Drescher, K.; Leptos, K.C.; Tuval, I. Dancing volvox: Hydrodynamic bound states of swimming algae. Phys. Rev. Lett.; 2009; 102, 168101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.102.168101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19518757]

16. Ishikawa, T.; Hota, M. Interaction of two swimming Paramecia. J. Exp. Biol.; 2006; 209, pp. 4452-4463. [DOI: https://dx.doi.org/10.1242/jeb.02537]

17. Chen, X.; Yang, X.; Yang, M. Dynamic clustering in suspension of motile bacteria. Europhys. Lett.; 2015; 111, 54002. [DOI: https://dx.doi.org/10.1209/0295-5075/111/54002]

18. Petroff, A.P.; Wu, X.L.; Libchaber, A. Fast-moving bacteria self-organize into active two-dimensional crystals of rotating cells. Phys. Rev. Lett.; 2015; 114, 158102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.114.158102]

19. Saintillan, D.; Shelley, M.J. Orientational order and instabilities in suspensions of self-locomoting rods. Phys. Rev. Lett.; 2007; 99, 058102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.99.058102]

20. Lushi, E.; Wioland, H.; Goldstein, R.E. Fluid flows created by swimming bacteria drive self-organization in confined suspensions. Proc. Natl. Acad. Sci. USA; 2014; 111, pp. 9733-9738. [DOI: https://dx.doi.org/10.1073/pnas.1405698111]

21. Papavassiliou, D.; Alexander, G.P. Exact solutions for hydrodynamic interactions of two squirming spheres. J. Fluid Mech.; 2017; 813, pp. 618-646. [DOI: https://dx.doi.org/10.1017/jfm.2016.837]

22. More, R.V.; Ardekani, A.M. Hydrodynamic interactions between swimming microorganisms in a linearly density stratified fluid. Phys. Rev. E; 2021; 103, 013109. [DOI: https://dx.doi.org/10.1103/PhysRevE.103.013109] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33601564]

23. Qi, T.T.; Lin, J.Z.; Ouyang, Z. Hydrodynamic behavior of self-propelled particles in a simple shear flow. Entropy; 2022; 24, 854. [DOI: https://dx.doi.org/10.3390/e24070854]

24. Guasto, J.S.; Rusconi, R.; Stocker, R. Fluid mechanics of planktonic microorganisms. Annu. Rev. Fluid Mech.; 2012; 44, pp. 373-400. [DOI: https://dx.doi.org/10.1146/annurev-fluid-120710-101156]

25. Ouillon, R.; Houghton, I.A.; Dabiri, J.O. Active swimmers interacting with stratified fluids during collective vertical migration. J. Fluid Mech.; 2020; 902, A23. [DOI: https://dx.doi.org/10.1017/jfm.2020.618]

26. Hill, N.A.; Pedley, T.J. Bioconvection. Fluid Dyn. Res.; 2005; 37, 1. [DOI: https://dx.doi.org/10.1016/j.fluiddyn.2005.03.002]

27. Pedley, T.J.; Kessler, J.O. Hydrodynamic phenomena in suspensions of swimming microorganisms. Annu. Rev. Fluid Mech.; 1992; 24, pp. 313-358. [DOI: https://dx.doi.org/10.1146/annurev.fl.24.010192.001525]

28. Fenchel, T.O.M.; Finlay, B.J. Photobehavior of the Ciliated Protozoon Loxodes: Taxic, Transient, and Kinetic Responses in the Presence and Absence of Oxygen 1. J. Protozool.; 1986; 33, pp. 139-145. [DOI: https://dx.doi.org/10.1111/j.1550-7408.1986.tb05578.x]

29. Ntefidou, M.; Iseki, M.; Watanabe, M. Photoactivated adenylyl cyclase controls phototaxis in the flagellate Euglena gracilis. Plant Physiol.; 2003; 133, pp. 1517-1521. [DOI: https://dx.doi.org/10.1104/pp.103.034223]

30. Qian, Y.H.; D’Humières, D.; Lallemand, P. Lattice BGK models for Navier-Stokes equation. Europhys. Lett.; 1992; 17, pp. 479-484. [DOI: https://dx.doi.org/10.1209/0295-5075/17/6/001]

31. Benzi, R.; Succi, S.; Vergassola, M.R. The lattice Boltzmann equation: Theory and applications. Phys. Rep.; 1992; 222, pp. 145-197. [DOI: https://dx.doi.org/10.1016/0370-1573(92)90090-M]

32. Chen, S.Y.; Doolen, G.D. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech.; 1998; 30, pp. 329-364. [DOI: https://dx.doi.org/10.1146/annurev.fluid.30.1.329]

33. Lallemand, P.; Luo, L.S. Lattice Boltzmann method for moving boundaries. J. Comput. Phys.; 2003; 184, pp. 406-421. [DOI: https://dx.doi.org/10.1016/S0021-9991(02)00022-0]

34. Mei, R.; Yu, D.; Shyy, W. Force evaluation in the lattice Boltzmann method involving curved geometry. Phys. Rev. E; 2002; 65, 041203. [DOI: https://dx.doi.org/10.1103/PhysRevE.65.041203]

35. Aidun, C.K.; Lu, Y.; Ding, E.J. Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation. J. Fluid Mech.; 1998; 373, pp. 287-311. [DOI: https://dx.doi.org/10.1017/S0022112098002493]

36. Glowinski, R.; Pan, T.; Hesla, T.I. A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: Application to particulate flow. J. Comput. Phys.; 2001; 169, pp. 363-426. [DOI: https://dx.doi.org/10.1006/jcph.2000.6542]

37. Feng, Z.G.; Michaelides, E.E. The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems. J. Comput. Phys.; 2004; 195, pp. 602-628. [DOI: https://dx.doi.org/10.1016/j.jcp.2003.10.013]