Gaofeng Chen , Xikai Jiang, *

a State Key Laboratory of Nonlinear Mechanics, Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China

b School of Engineering Science, University of Chinese Academy of Sciences, Beijing 10 0 049, China

Keywords: Particle dynamics Low-Reynolds-number fluid Confinement Mobility Drift

ABSTRACT Dynamics of a spherical particle and the suspending low-Reynolds-number fluid confined by a cubic cav- ity were studied numerically. We calculated the particle’s hydrodynamic mobilities along x -, y -, and z- directions at various locations in the cavity. The mobility is largest in the cavity center and decays as the particle becomes closer to no-slip walls. It was found that mobilities in the entire cubic cavity can be determined by a minimal set in a unit tetrahedron therein. Fluid vortices in the cavity induced by the particle motion were observed and analyzed. We also found that the particle can exhibit a drift motion perpendicular to the external force. Magnitude of the drift velocity normalized by the velocity along the direction of the external force depends on particle location and particle-to-cavity sizes ratio. This work forms the basis to understand more complex dynamics in microfluidic applications such as intracellular transport and encapsulation technologies.

Particulate transport in the low-Reynolds-number fluid under confinement finds a wide variety of applications ranging from particle manipulation in lab-on-a-chip devices to macromolecular transport in biological tissues [1–3] . Recently, particle dynamics under total confinement have drawn much attention because it is important to microfluidic encapsulation techonologies and bio- logical functions in living cells [4–9] . In cells, biological particles in the cytoplasmic fluid are confined by the cell membrane, and their confined motion can affect intracellular activities such as cel- lular homeostasis, signaling, and metabolism [10,11] . Experimen- tal studies have revealed intriguing transport phenomena in the cytoplasm such as cytoplasmic streaming and hindered diffusion of proteins [9,12] . To understand mechanisms behind these phe- nomena, numerical and theoretical studies have been carried out over the past decades. In these modeling works, the cytoplasm was treated as a continuum fluid confined in a cavity [13,14] . Through Brownian and Stokesian dynamics simulations, effects of crowding, hydrodynamic and steric interactions, and confinement were found to play vital roles in active and passive motion of intracellular par- ticles [4,5,7,8,15] .

Single-particle dynamics in the low-Reynolds-number fluid un- der total confinement is also important, as it forms the basis to understand more complex dynamics in concentrated suspensions. In most prior studies, particles were considered to be confined in a spherical cavity, mainly because the shapes of many cells and droplet-based microreactors are spherical and this geometry makes it simple for analyzing results while revealing key features of in- tracellular transport [8,16–19] . However, less attention have been paid on cavities with other shapes. One of the shapes is the cube, which can be found in plant cells, cuboidal epithelium, and micro- capsules for drug delivery systems [20–22] . In a prior study, dy- namics of a spherical and an ellipsoidal particle in a low-Reynolds- number fluid confined by a cubic cavity were examined by bound- ary element method. The particles’ friction coefficients at the cav- ity center and their settling translational and angular velocities at different locations in the cavity were calculated [23] . Despite these progress, hydrodynamic mobilities of the particle at other locations in the cubic cavity, the particle’s drift motion perpendicular to the external force, and distributions of fluid velocities in the cavity re- main to be explored. In this work, we address these issues by nu- merical simulations.

We consider a rigid spherical particle in a quiescent low- Reynolds-number fluid confined by a cubic cavity. The schematic of the simulation system is shown in Fig. 1 a. The radius of the spherical particle isR, and the side length of the cubic cavity is 2h. Cartesian coordinate system is used in this study, and the origin (0, 0, 0) is the cavity center. Neglecting inertial effects and Brownian motion, the balance equation for force and torque on the particles is

Fig. 1. (a) Schematic of the simulation system for studying dynamics of a spherical particle and the suspending low-Reynolds-number fluid confined by a cubic cavity; (b) Discretization of the spherical particle.

Fig. 2. Normalized mobility of the spherical particle along x -, y -, or z-direction at the cavity center plotted against particle-to-cavity sizes ratio R/h . Blue line is our numerical result. Squares (circles) are results from boundary (finite) element method in Ref. [23] .

whereFHis the hydrodynamic force/torque vector andFextcon- tains external forces/torques. In our numerical method, the particle surface is discretized into a set ofNnodes. Each node is connected with its neighbouring nodes by elastic springs with a large spring constant, in order to model the rigid particle and maintain particle shape. The nodes are also connected with the particle center-of- mass by stiffsprings to avoid particle deformation [24] . Discretiza- tion of the spherical particle is shown in Fig. 1 b. Equation (1) is then translated into theNsurface nodes as

Fig. 3. (a - c) Two-dimensional distributions of M x , M y , and M z for a particle with R/h = 0 . 1 in the plane defined by 0 ≤x/h ≤1 , 0 ≤y/h ≤1 , and z = 0 . (d - f) M x , M y , and M z under different R/h versus scaled particle x -position along the line defined by 0 ≤x/h ≤1 and y = z = 0 . (g - i) M x , M y , and M z under different R/h versus scaled particle x -position along the line defined by 0 ≤x/h ≤1 , x = y , and z = 0 . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

here,i= 1,2,...,N,is the hydrodynamic force,is the spring force, andare external forces. The spring force acting on theith node by thejth node is

wherekis the spring constant,r0is the equilibrium spring length for each spring,rij=ri?rj, andriandrjare coordinates of theith andjth nodes. The external force on theith node is=Fext/N.

The equation of motion for all nodes on the particle surface is

whereR=(r1,r2,...,rN)denotes a 3Nvector containing nodal co- ordinates,U0denotes a 3Nvector with the undisturbed fluid ve- locity at nodal positions, andMis the 3N×3Nmobility tensor.U=(u1,u2,...,uN)=M·Fcontains 3Ndisturbed velocities from the hydrodynamic interaction (HI), andF=(f1,f2,...,fN)is a 3Nvector including non-HI forces on the surface nodes. The transla- tional and rotational motion of the rigid particle are realized by integrating Eq. (4) for all surface nodes, satisfying the balance equation for force and torque. The velocity fieldUdriven by nodal forcesFcan be obtained by solving the Stokes equation

wherepis the fluid pressure,μisthe fluidviscosity,uis thefluid velocity, andf(r)=fi(ri)δ(r?ri)is the force density exert- ing on the fluid. To avoid singularity due to point forces placed at nodes on the particle surface, the smoothing functionδ(r)is used to regularize the point forces, and it takes the form of a modified Gaussian function

The regularization parameterξis related to the characteristic node spacing on the particle surfaceh, i.e.,ξ～h?1. This is to ensure that the regularized force density is spread over the length scale of as- sociated surface elements on the particle, preventing the fluid from penetrating the particle surface [5,25,26] . For the undisturbed fluid velocity in Eq. (4) , we setU0 = 0 to model the quiescent fluid. No-slip boundary condition is applied on all walls of the cavity.

The aforementioned equations for particle dynamics and fluid motion are solved using the General Geometry Ewald-like Method. Details of the method can be found in Ref. [24] . The method has been validated for the sedimentation of a spherical particle be- tween two parallel walls [5] . It has also been used to study col- lision and segregation behavior of fluid-filled elastic capsules in confined simple shear flows [25,27] . In the remainder of this work, we will first validate the method against numerical solutions for the hydrodynamic mobility of a spherical particle at the center of the cubic cavity; we will then apply it to study dynamics of the particle and the suspending fluid under cubic confinement. Dimen- sionless variables are used in this study, and they are based on a set of characteristic scales. The characteristic length scale isl(hy- drodynamic radius of each node), the energy scale iskBTwherekBis Boltzmann constant andTis temperature, and the force scale iskBT/l.

We first calculate hydrodynamic mobility of the spherical parti- cle at the cavity center normalized by that of the same particle in unbounded fluid. Due to the symmetry of the cubic cavity, normal- ized mobilities alongx-,y-, andz-directions (Mx,My,Mz) are the same at the cavity center. Results for one of the mobilities under various particle-to-cavity sizes ratios (R/h) are shown in Fig. 2 . AsR/hincreases (particle size increases or cavity size decreases), nor- malized mobility of the particle decreases due to increased level of confinement and greater influence of no-slip walls on the particle. Results from Ref. [23] are also included in Fig. 2 , and our numerical results agree well with those reported in the literature.

We then calculateMx,My, andMzof a particle withR/h= 0.1 in a quarter of a midplane of the cubic cavity, which is defined by 0 ≤x/h≤1 , 0 ≤y/h≤1 , andz= 0 . Due to the symmetry of the cubic geometry, mobilities at the other three quarters of the midplane can be determined by those in the selected one. Two- dimensional distributions of the mobilities in the chosen plane are shown in Fig. 3 a–3 c. It can be seen that, in general, the mobil- ity is largest in the cavity center and decays as the particle be- comes closer to no-slip walls.MxandMyare found to be sym- metric about the diagonal of the chosen plane, whileMzitself is symmetric about the diagonal. We further select two lines in the plane and analyze variations of mobilities along them. The first line is defined by 0 ≤x/h≤1 andy=z= 0 , and relevant results are shown in Fig. 3 d–3 f; the second line is defined by 0 ≤x/h≤1 ,x=y, andz= 0 , and relevant results are shown in Fig. 3 g–3 i. In general, the mobility at a certain location decreases asR/hin- creases because the confinement level increases and the no-slip walls exert greater influence on the particle. On the first line,Mxis perpendicular to the wall atx=h, whileMyandMzare parallel to this wall andMy=Mz. Under a specificR/h,Mxis smaller thanMy(Mz) when the particle is between cavity center and cavity wall atx=h. As the particlex-position increases, mobilities decrease slowly (rapidly) near cavity center (wall). The decreasing trends forMxandMy(Mz) withx/hare different: the slowly decaying re- gions forMyandMzare wider than that forMx. It is also observed that, near the wall, mobilities decrease more abruptly asR/hin- creases. On the second line,Mx(My) is perpendicular to the wall atx=h(y=h)andMx=My, whileMzis parallel to these walls. Ob- servations on the first line for mobilities perpendicular and parallel to the wall also apply to those on the second line. One difference is that, as the particle on the second line is affected more greatly by the wall aty=h,Mzon the second line is smaller than that on the first line at the samex-position.

Fig. 4. Two-dimensional distributions of fluid velocity in x - and y -directions (u and v) and streamlines in the midplane of the cubic cavity at z = 0 . Red (blue) color denotes positive (negative) value. (a - c) Particle center is at the origin and the external force (F) is in positive x -direction. (d - f) Particle center is at (0 . 5 h , 0, 0) and F is in positive x -direction. (g - i) Particle center is at (0 . 5 h , 0, 0) and F is in positive y -direction. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Based on symmetric properties of the cubic geometry, we spec- ulate that particle mobilities in the entire cavity can be reduced to a minimal set in a subregion of the cavity. To find this subregion, we start by dividing the cubic cavity into eight subcubes using three midplanes atx= 0 ,y= 0 , andz= 0 . Mobilities in the entire cavity can then be determined by those in only one of the sub- cubes. We choose a representative subcube defined by the origin and its diagonal vertex (h,h,h) as shown in Fig. 1 a. Next, we con- sider the diagonal plane of symmetry of the cubic cavity, and the chosen subcube can be further divided into two triangular prisms. A representative prism defined by (0, 0, 0), (h,0,0), (h,h,0), (0, 0,h), (h,0,h), and (h,h,h) is chosen, and mobilities in the other prism can be calculated based on those in the chosen prism. Fi- nally, by taking advantage of the commutative property of mo- bilities in the cubic cavity, the chosen prism can be divided into three tetrahedra. The commutative property is related to diagonal planes of symmetry, and one example isMx(a,b,c)=Mx(a,c,b)=My(b,a,c)=My(c,a,b)=Mz(b,c,a)=Mz(c,b,a)wherea,b, andcare components of the particle center’s coordinate. One represen- tative tetrahedron defined by (0, 0, 0), (h, 0, 0), (h, 0,h), and (h,h,h) is marked out by yellow lines in Fig. 1 a. We arrive at the min- imal set formed by mobilities in one of the tetrahedra, which can be used to determine mobilities in the entire cavity.

We also analyze fluid velocities in the cavity induced by the particle motion. Three cases are considered: in the first, the parti- cle is at the cavity center and the external force is along positivex-direction; in the second, the particle is at (0.5h, 0, 0) and the ex- ternal force is the same as that in the first case; in the third, the particle location is the same as that in the second case while the external force is along positivey-direction. We choose a midplane of the cubic cavity atz= 0 , and analyze fluid velocities inx- andy-directions (uandv) and streamlines in this plane. As shown in Fig. 4 a,uis positive in the vicinity of the particle, but is negative near walls perpendicular toy-axis. From Fig. 4 b, we can see thatvis diagonally symmetric about the particle, withvbeing positive in the upper right and lower left regions of the plane and negative in the upper left and lower right regions. Fig. 4 c shows streamlines, and two vortices are observed near the particle. The line connect- ing two stagnation points in the centers of vortices is found to be perpendicular to the applied force. In the second case, the overall pattern of fluid velocity distributions is similar to those in the first case as shown in Fig. 4 d–f. Nonetheless, the distribution of nega- tiveubecomes more blunt near the wall atx=hin the second case;vis no longer diagonally symmetric about the particle, and its distribution near the wall becomes elongated alongy-direction compared to that in the first case. Fig. 4 g–i show results in the third case. As the particle is near the wall and the force direction changes, the overall flow pattern becomes different from those in previous cases. The maximal magnitude ofuon the left of the par- ticle is larger than that on the right; the region of nonzero fluid velocity on the left is wider than that on the right.vis positive near the particle, but is negative near the wall on the left of the particle. From the streamlines, only one vortex is observed on the left of the particle.

Fig. 5. Critical particle position x c /h for vortex transition as a function of R/h . To determine x c /h , the particle is placed on multiple locations along x -axis and an ex- ternal force in y -direction is applied. When 0 ≤x/h < x c /h , two vortices exist in the xy plane; when x c /h < x/h < 1 , only one vortex exits.

To determine the critical particle positionxc/hfor the transition from two vortices to one vortex, we place the particle on multiple locations alongx-axis and apply an external force iny-direction. By counting the number of vortices in thexyplane for each simula- tion, it is found thatxc/hdecreases with increasingR/has shown in Fig. 5 . This is because the no-slip wall nearest to the particle exerts greater influence on the particle with higherR/h, making the vortex between the particle and its nearest wall disappear at a smallerxc/h.

Fig. 6. (a) Two-dimensional distribution of normalized drift velocity along y -direction (U y /U x) for a particle with R/h = 0 . 2 in the plane defined by 0 ≤x/h ≤1 , y/h = 0 . 75 , and 0 ≤z/h ≤1 . The external force is along x -direction. (b) Normalized drift velocity of a particle with y/h = ?0 . 75 and z = 0 plotted against x/h for different R/h . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Prior studies found that the particle will drift perpendicular to the external force because of the anisotropy of the mobility tensor induced by no-slip walls [28] . Here we study the drift motion of the spherical particle in a cubic cavity. The particle is placed in a plane defined by 0 ≤x/h≤1,y/h= 0.75 , and 0 ≤z/h≤1 . The ex- ternal force alongx-direction is applied on the particle, and the main velocity alongx-direction (Ux) and drift velocity alongy- direction (Uy) are measured. Fig. 6 a shows two-dimensional distri- bution of the scaled drift velocity (Uy/Ux) in the plane. In general, the drift velocity near the wall atx=his larger than that near the cavity center. To examine the effect of confinement level on the drift velocity, we calculate variations ofUy/Ux, under variousR/h, along a line defined by 0 ≤x/h≤1,y/h= ?0.75 , andz= 0 , and results are presented in Fig. 6 b. We find that the scaled drift velocity at a certainx/hincreases asR/hincreases. For allR/hcon- sidered here, the increasing rate of the scaled drift velocity withx/hincreases in the interior of the cavity, but decreases when the particle gets very close to the wall. ForR/h= 0.1 , the scaled drift velocity increases (decreases) whenx/h<0.825 (x/h>0.825). Asx/hincreases and the particle moves further away from the sym- metry plane (midplane atz= 0), the flow around the particle be- comes more asymmetric with respect to the plane aty/h= ?0.75 . As a result, forces that the fluid exerts on the particle’s two hemi- spheres alongy-direction become more imbalanced, driving the magnitude of drift velocity up. On the other hand, as the particle moves closer to the wall, particle velocity tends to decreases due to the greater influence of the no-slip wall. Thus, the competition between effects of flow asymmetry and no-slip wall leads to the nonlinear variation of the drift velocity.

To conclude, we have studied hydrodynamic mobility and drift motion of a spherical particle in a low-Reynolds-number fluid con- fined by a cubic cavity using numerical simulations. Fluid veloci- ties induced by the particle motion were also analyzed for differ- ent particle positions and external forces. We found that mobili- ties in the entire cavity can be computed based on a minimal set in a unit tetrahedron in the cavity. Fluid vortices in the cavity in- duced by the particle motion were analyzed. It was also found that the particle’s drift velocity perpendicular to the external force de- pends on particle position and particle-to-cavity sizes ratio. This work forms the basis for understanding more complex particle dy- namics under total confinement, which could benefit microfluidic applications such as intracellular transport and encapsulation tech- nologies.

Declaration of Competing Interest

The authors declare that they have no known competing finan- cial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the Young Elite Scientists Spon- sorship Program by the Chinese Society of Theoretical and Applied Mechanics (CSTAM).