Shape, confinement and inertia effects on the dynamics of a driven spheroid in a viscous fluid

Abstract

The dynamics of anisotropic particles in viscous flows underpin a wide range of processes in soft matter, microfluidics, and targeted drug delivery. Here, we investigate the motion of externally driven prolate and oblate spheroids suspended in a Newtonian fluid and confined within a square microchannel. Using lattice Boltzmann simulations, complemented by far-field hydrodynamic theory based on superposition of wall interactions, we systematically quantify how particle aspect ratio, strength of confinement, and fluid inertia influence the dynamics of a spheroid. For unconfined spheroids, we show that the translational velocity is maximised not for a sphere but for a prolate (end-on) or oblate (broadside-on) spheroid of specific aspect ratio. Under confinement, the optimal aspect ratio shifts toward oblate shapes due to the dominant contribution of wall-induced frictional resistance. Off-centre positioning introduces strong translation–rotation coupling, giving rise to two families of oscillatory trajectories—glancing and reversing—whose existence and structure are captured as closed orbits in the phase space. Weak fluid inertia breaks these closed loops: glancing trajectories spiral outward and merge with reversing trajectories, and new stable fixed points emerge. Together, these results reveal how modest deviations from sphericity or creeping-flow conditions profoundly alter the dynamics of driven particles in confined geometries. The predictions offer guidelines for optimising particle shape in microfluidic transport and highlight the rich nonlinear behaviour accessible in confined suspensions of nonspherical colloids.

---

1 Introduction

Dynamics of anisotropic particles, whether driven or self-propelled, is a problem of fundamental importance in fluid mechanics and statistical physics, and relevant to various engineering applications. Because of its close connection to fluctuation dissipation theorem,¹ tracking colloidal particles forms the foundation of microrheological measurements in soft and biological matter. Recent investigations have revealed the hydrodynamic consequences emerging from heterogeneous environments such as solid walls and interfaces on the Brownian dynamics of colloids. However, these studies have predominantly focused on spherical particles.^2–5 In contrast, the role of particle shape in these problems remains comparatively underexplored.

Targeted drug delivery represents one of the most promising applications of driven colloids,^6–10 yet most work to date has concentrated on spherical particles, largely due to the relative simplicity of their fabrication and their suitability as model systems.^8,11 Advances in synthesis methods, however, now enable the fabrication of nonspherical colloids, thereby opening new avenues of research.^12–14 Recent studies have demonstrated therapeutic advantages of nonspherical particles in cancer treatment,^15–17 but their broader potential as versatile drug carriers remains insufficiently investigated. Since such particles must move in biological fluids and complex, heterogeneous environments within the body,^18–21 it is essential to determine whether, and in what ways, particle shape plays a role.

Analysis of translational dynamics of a driven spheroid immersed in an unbounded, Newtonian, viscous fluid has a long history with exact solutions available in the low Reynolds number limit (Stokes flow).^22–24 The Reynolds number describes the ratio of inertial to viscous forces in the fluid, Re = Ua/ν where U and a are the characteristic velocity and size (the semi-major axis) of the spheroid respectively, and ν is the kinematic viscosity of the surrounding fluid. Limiting cases of the spheroid geometry – slender bodies, spheres and discs – are most often studied, and are amenable to analytical calculations.²⁵ Despite the long history, researchers continue to find the dynamics of spheroidal particles fascinating.^26–28

In an otherwise quiescent, infinite fluid at low Reynolds number, the dynamics of spheroids moving under an external force has two interesting characteristics:²⁹ (i) they do not undergo rotation due to reversibility constraints, and (ii) as they translate in the direction of the applied force, they also migrate transversely (perpendicular to the applied force). However, these characteristics are not preserved in the vicinity of a solid boundary, where hydrodynamic interactions with the wall generate additional torques and transverse forces.³⁰ The result of the consequent coupled translational and rotational motion is oscillatory dynamics of the spheroid. Depending on the initial conditions, spheroids may undergo ‘glancing’, ‘reversing’ or ‘periodic tumbling’ dynamics.³¹ Near a slanting wall additional modes such as sliding or wobbling motion may also be observed.³¹

Similarly, normal stresses arising from fluid inertia can exert torque on spheroids. For an unconfined spheroid, inertial torque tends to reorient its long axis perpendicular to the direction of the applied force.^32–36 Stronger inertial effects can induce wake instabilities and unsteady modes near the translating spheroid. In the case of a confined spheroid, both inertial forces and hydrodynamic interactions with the confining boundaries act simultaneously. Such conditions commonly arise in driven colloids and particles used for targeted drug delivery, where both the degree of confinement and the Reynolds number can vary significantly. In the present work, we investigate the effect of confinement and weak inertia ()—separately and in combination—on the translational dynamics of spheroids.

Based on the analysis near a single wall, Russel et al.³⁰ suggested that spheroids would exhibit oscillatory trajectories between two parallel walls. The reversibility constraints of Stokes flow imply that such glancing-reversing trajectories³¹ must be dependent on the initial conditions. For Re > 0, investigations of spheroids translating in narrow tubes^37–39 have revealed several distinct modes of motion such as translation at constant orientation, planar oscillations and spiralling trajectories. These modes are found to be independent of the initial conditions. This suggests the existence of a Reynolds-number-dependent transition in the translational modes of spheroids confined within a channel. Limited numerical studies in cylindrical channels have shown Reynolds-number-dependent radial positions in the cylinder at which spheroids translate without reorienting,⁴⁰ although it remains unclear whether these stable solutions arise specifically from the symmetry of the cylindrical geometry.

In confinement and other complex geometries, numerical techniques are required to analyse the dynamics of spheroids. Several techniques have been developed, which include the constrained-force technique,⁴¹ arbitrary Lagrangian–Eulerian based finite-element method,^33,40 direct numerical simulation (DNS) with Lagrange multiplier based fictitious domain methodology,^34,42 Stokesian dynamics,⁴³ boundary element method,³¹ immersed boundary method with DNS,⁴⁴ and the lattice Boltzmann method (LBM).^37–39,45

In this work, we employ a recently developed lattice Boltzmann framework⁴⁵ to study a model system: the dynamics of prolate and oblate spheroids immersed in a viscous fluid and confined within a channel of square cross section. The particle is driven by an external magnetic²¹ or gravitational field, which applies a force on the particle. We systematically investigate the roles of confinement and inertia on the translational dynamics of the spheroid. We consider only driven particles (by an external force) in this work which are experimentally easy to realise, say, compared to active particles for which the leading force description is a force dipole. We find that for a given volume, the maximum translational velocity of an unconfined spheroid does not occur for a sphere but rather for a prolate or an oblate shape depending on its orientation relative to the applied force. In confinement, the optimal aspect ratio shifts toward oblate shapes, highlighting a promising direction for future experimental studies in drug delivery applications. Furthermore, analysis of trajectories reveals that spheroids display glancing–reversing motions depending on the initial conditions, aspect ratio, and degree of confinement. At finite Reynolds numbers, inertial nonlinearities alter this behaviour, and by systematically increasing the Reynolds number from the Stokes regime to , we uncover the emergence of rich nonlinear dynamics in the motion of confined spheroids.

The paper is organised as follows. Section 2 introduces the model system and briefly reviews the existing literature of the analytical framework for translating spheroids in unconfined systems. We will discuss an approximate calculation based on the method of images in Section 2.2. This will be followed in Section 3 by the description of the numerical method based on lattice Boltzmann method used in this work to analyse the dynamics of confined spheroids. In Section 4 we will discuss the results in the following order: optimal aspect ratio for the unconfined and confined spheroids, effect of confinement on the glancing-reversing trajectories, and their phase-space description. We will end by highlighting the nonlinear features on the translational dynamics of the spheroid that emerge due to fluid inertia and then conclude with a discussion on the scope of the work in Section 5.

2 Theory

Consider a spheroid translating under the action of an external force F_e in an otherwise quiescent fluid as shown in Fig. 1a. The spheroid and the fluid are confined in a channel with a square cross section of side length L. The instantaneous orientation of the major axis of the ellipsoid is indicated by a unit vector ê. Let the spheroid be of semi-major axis a and semi-minor axes b = c. For prolate spheroids, a > b = c, and its eccentricity is defined as ; . Similarly, for oblate spheroids, a < b = c, and the eccentricity is ; .

Fig. 1 (a) Schematic representation of a translating spheroid with semi-major axis a and semi-minor axes b = c, located at (X, Y = h, Z = L/2), oriented at θ, experiencing an external force F_e, confined in a channel of square cross section of side length L. The side view of the square channel (Y–Z cross section) is shown in the inset. The walls are numbered as 1, 2, 3 & 4. (b) Fluid flow generated by the translating spheroid in the square channel, in an otherwise quiescent Newtonian fluid. The arrows indicate the velocity field, and are coloured according to the magnitude of velocity with red being the highest velocity.

The fluid motion is described by the incompressible Navier–Stokes equations:

∇·u = 0 (1)

ρ(∂_tu + u·∇u) = −∇P + η∇²u (2)

Here, P is the pressure field, u is the fluid velocity field, ρ is the density and η = νρ is the dynamic viscosity and ν is the kinematic viscosity of the fluid.

The force F_e may be any body force applied externally and acts on the particle to drive the particle to translational motion. For example, the external force field can be considered to arise from a magnetic field when the spheroid is a microrobot, or gravitational field when gravitational sedimentation is considered.

We consider the translation of unconfined and confined spheroids. Here, ‘unconfined’ refers to the translation of the spheroids in the bulk fluid where wall effects are not present. ‘Confined’ spheroids refer to particles that are enclosed within the four walls of the square channel, as shown in Fig. 1a, and thus the hydrodynamic interactions with the walls play a dominant role in their dynamics. The strength of the confinement is described by defining a confinement ratio, CR = 2b/L where L is the height (equal to the depth) of the channel.

2.1 Translation of the unconfined spheroid

We briefly discuss the translational dynamics of an unconfined spheroid first. At zero Reynolds number, the steady state velocity of a single, externally driven spheroid can be exactly determined^22,24,46,47 as

F_e = 6πηa(C_‖U_‖ê + C_⊥U_⊥ê_⊥) (3)

where ê and ê_⊥ are the unit vectors along and perpendicular to the long axis of the spheroid, respectively, U_‖ and U_⊥ are the translational velocities along ê and ê_⊥, respectively, i.e., U = U_‖ê + U_⊥ê_⊥, and C_‖ and C_⊥ are the corresponding drag coefficients given by ref. 48,for the prolate spheroids andfor the oblate spheroids. Under the slender-body approximation, i.e., for b/a → 0, eqn (4) and (5) reduce to and respectively. Similarly, in the limit of a disc, i.e., a/b → 0 eqn (6) and (7) reduce to and , respectively.²⁴

When the external force F_e is along the major or minor axis of the spheroid (along ê or ê_⊥), eqn (3) reduces to a single scalar equation relating the applied force and the corresponding translational velocity. For any other orientations of the spheroid, eqn (3) shows that the anisotropy in the drag coefficients, i.e., C_‖ ≠ C_⊥, results in transverse forces and the migration of the ellipsoid in a direction perpendicular to the direction of its translation.²⁹

2.2 Translation of the confined spheroid

Confining a force-driven spheroid in a channel can induce complex fluid flows and interesting dynamics. An example of the velocity field generated by a translating spheroid (produced using lattice Boltzmann simulations as mentioned in Section 3) is shown in Fig. 1b. In other words, hydrodynamic interactions of the spheroid with confining walls alter its translational velocity and induce an angular velocity, thus leading to coupled translational–rotational dynamics of the spheroid.

In this work, we only consider the case when F_e is along the direction of the channel length x̂ i.e., the externally applied force on the spheroid F_e is parallel to the confining walls. For simplicity, we place the particle at the midplane of the channel in the z direction, i.e., at z = L/2 and with long axis lying entirely in this midplane, i.e., the orientation vector ê in the x–y plane. Owing to the symmetry of this initial condition, and despite the coupled translational–rotational dynamics, the spheroid stays and moves only in this midplane. The complete solution of this problem is obtained numerically as described in Section 3. However, an approximate solution of the dynamics of the spheroid can be developed as follows.

Using the method of images expressions can be derived for the translational and rotational velocities of prolate and oblate spheroids placed near a wall and driven by an external force.³¹ These are far-field approximations at zero Reynolds number and near a single wall, series expressions in terms of inverse power of the distance from the wall and are valid for arbitrary orientations. Assuming a superposition principle,⁴⁹ namely that the hydrodynamic interactions of the spheroid with each wall are additive, but neglecting the corrections due to the simultaneous presence of multiple walls, and multiple reflections between the walls, we obtain

where U_bulk = U_‖ê + U_⊥ê_⊥ given by eqn (3). There are four confining walls for the channel as shown in Fig. 1a, and the translational and rotational velocity of the spheroid due to the hydrodynamic interaction with each wall is indicated by U^i,wall and Ω^i,wall. They are calculated as U^i,wall = F_e/(6πμR)Ũ^i,wall, Ω^i,wall = F_e/(6πμR²)[capital Omega, Greek, tilde]^i,wall, where R = a for prolate and R = b for oblate spheroids, F_e = |F_e|, and:The corresponding expressions for the wall on the right (i = 2) can be obtained by appropriate substitutions:

Ũ^2,wall_x = Ũ^1,wall_x(h → L − h, θ → −θ) (13)

Ũ^2,wall_y = −Ũ^1,wall_y(h → L − h, θ → −θ) (14)

[capital Omega, Greek, tilde]^2,wall_z = −[capital Omega, Greek, tilde]^1,wall_z(h → L − h, θ → −θ). (15)

The hydrodynamic interaction with the remaining pair of walls (i = 3, 4) are given by,

Ũ^3,wall_y = Ũ^4,wall_y = 0 (17)

[capital Omega, Greek, tilde]^3,wall_z = [capital Omega, Greek, tilde]^4,wall_z = 0. (18)

Analysing eqn (10)–(12) we can see that the far-field prediction is for translational and for angular velocity. Eqn (8)–(18) can be used to determine approximately the coupled translational–rotational dynamics of the spheroid. The complete dynamics can be captured using numerical simulations as described in the next section.

3 Simulation method

We use the lattice Boltzmann method to determine the dynamics of the fluid around the translating spheroid confined in the channel. Through the boundary conditions the motion of the spheroid and the fluid are completely coupled. The algorithm is implemented in our parallel lattice Boltzmann code Ludwig.^50,51 We briefly discuss the method below, and refer the reader to Thampi et al.⁴⁵ for a detailed discussion of the implementation.

3.1 Lattice Boltzmann method

The lattice Boltzmann method (LBM) defines a discrete velocity distribution function f_i(x, t) at each grid point x and time t that travels along the ith discrete direction with velocity c_i. The hydrodynamic observables, namely the density ρ(x, t) and momentum ρ(x, t)u(x, t) are obtained as the zeroth and first moment of the discrete distribution function:

In this work, the D3Q19 model is used, which discretises configuration and momentum space with a 3D cubic lattice with 19 discrete lattice vectors c_i.⁵² A lattice spacing of Δx = 1 and time step Δt = 1 are used for the spatial and temporal discretisation, respectively.

The discrete distribution functions f_i undergo successive collision and propagation operations:

f_i(x + c_iΔt, t + Δt) − f_i(x, t) = ΔtC_i(t). (21)

where in its simplest Bhatnagar–Gross–Krook form ref. 53 the collision operator contains a single relaxation time τ and is given as

C_i(t) = [f_i(x, t) − f^eq_i(x, t)]/τ (22)

The kinematic viscosity ν of the fluid is related to the relaxation time via . The equilibrium distribution:is constrained by the governing equations to be recovered, namely the Navier–Stokes equations. Here, c_s is the speed of sound and w^c_i are the weights associated with the velocity set c_i of the D3Q19 model. We use a generalised, multiple relaxation time (MRT) collision operator in our simulations.⁵⁴

3.2 Boundary conditions

The surface of the spheroid is described as

(x − x_c)^TA(x − x_c) = 1, (24)

where x_c is the centre of mass, x is any point on the surface of the spheroid and A is a 3 × 3 matrix whose eigenvalues are related to the inverse of the semi-axes of the spheroid.⁵⁵ Bounce-back boundary conditions are applied on the boundary nodes , which is an approximate surface of the spheroid.⁴⁵

The bounce back scheme results in momentum exchange between the fluid and the solid nodes. The total force F and torque T on the spheroid due to momentum exchange can be calculated by adding up the contributions from all boundary nodes as and , where⁵⁶

Here, is the post-collision distribution function, ρ₀ is the mean density of the fluid, and u_b is the local velocity of the boundary node calculated as

u_b = U + Ω × (x_b − x_c). (26)

where, U and Ω are the translational and angular velocities of the spheroid.

3.3 Dynamics of the spheroid

An implicit numerical scheme is used to update the translational and angular velocities of the spheroid from the force and torque determined in the previous section, accounting for particle inertia. In this scheme,⁵⁶ the total force and torque on the particle are written as,

F = F₀ − ζ^FU·U − ζ^FΩ·Ω (27)

T = T₀ − ζ^TU·U − ζ^TΩ·Ω (28)

where F₀ and T₀ are ‘velocity-independent’ forces and torques, determined from the post-collision distributions and ζ** are the drag coefficient matrices.⁴⁵

Then, the discretised conservation equations of linear and angular momentum of the translating spheroid read as

where M and I are respectively the mass and moment of inertia tensor of the spheroid. A Gaussian elimination method is used to solve eqn (29) and (30) simultaneously, thus ensuring the stability of the algorithm.

The position of the particle is then updated as

The orientation of the particle is updated as^57–59

q(t + Δt) = q̃(t)q(t), (33)

where q is the unit quaternion describing the orientation of the spheroid based on Euler angles,^60,61 is the mean angular velocity and ||·|| indicates the norm of the vector. This procedure avoids accumulated errors from successive matrix operations and singular matrix operations, which could arise if Euler angles were used to describe the orientation of the spheroid. It also avoids any possibility of a gimbal lock and requirement of renormalisation of quaternions to account for errors from numerical integration since the magnitude of the quaternion is exactly unity.

Finally, we note that in the implicit solution scheme in eqn (30) the rate of change of the moment of inertia of the spheroid is not zero in the laboratory coordinate system, which is why this calculation requires careful consideration. The time-dependent moment of inertia tensor I(t) can be determined from the quaternion q(t),^57,59,61

I(t) = (q(qIq⁻¹)^Tq⁻¹)^T, (34)

where I(t) = [0, I(t)] is a pure quaternion moment of inertia tensor based on the inertia tensor in the principal coordinate system. In the above, ⁻¹ and ^T represent the inverse and transpose operations respectively. Then the time derivative of the moment of inertia tensor in the second term in eqn (30) can be exactly determined.

3.4 Summary of the simulation method

The numerical procedure described in the previous section proceeds as follows. The lattice Boltzmann streaming and collision steps advance the fluid field around the moving spheroid. The mid-grid bounce-back condition is imposed on the particle surface, enabling momentum transfer between the spheroid and the surrounding fluid. This exchange enforces the no-slip boundary condition on both the spheroid and the channel walls while simultaneously providing the hydrodynamic force and torque acting on the particle. The resulting translational and rotational dynamics of the spheroid are then obtained by solving the linear and angular momentum balance equations, with the particle orientation updated using quaternions and standard quaternion operations. More details on the implementation of the algorithm and validation of the method have been reported previously.⁴⁵

The simulation parameters used in this study are summarised below. For the unconfined cases, a cubic computational domain of size 256³ with periodic boundary conditions applied on all boundaries is employed. The major and minor axes of the spheroid vary between 6–60 and 2–6, respectively. We have used particles with density ratio ρ_p/ρ_f = 1 in the simulations with an external body force specified using F_e. To investigate the role of geometric confinement, simulations are carried out in a domain of size 256 × N_y × N_z, where N_y = N_z is varied from 5–154 to vary the strength of confinement. When a spheroid is in close proximity to a wall, the hydrodynamic lubrication between the two can be under-represented on the discrete lattice. In this case a correction to the lubrication is made based on the surface separation at the point of closest approach. This prevents the spheroid touching the wall. Our results are independent of the specific implementation of this contact force. Periodic boundary conditions are applied on the open sides of the domain. To examine the influence of inertia, the kinematic viscosity of the fluid is varied 0.007–0.5. In all cases, simulations are initialised with a quiescent fluid. Each run is continued until a steady state is reached, and the corresponding results are presented below. In the Stokes flow regime, quasi-static simulations were employed wherever feasible to improve computational efficiency. In this approach, the fluid field is evolved alongside the particle velocities, but the particle is kept fixed in space. Holding the particle stationary eliminates unsteady effects and allows the hydrodynamic calculations to converge rapidly generating the instantaneous kinematic state for that configuration.

4 Results and discussion

We consider the dynamics of both prolate and oblate spheroids by varying the aspect ratio b/a → 0 to b/a → ∞. Below, we first consider the case of unconfined spheroids and demonstrate the existence of an optimal aspect ratio that yields the maximum translation velocity. We will then present the results of confined spheroids, and analyse the shift in the optimum aspect ratio as a function of strength of confinement. Furthermore, we present results that demonstrate the effects of confinement and inertia on the trajectories.

4.1 Optimum aspect ratio of an unconfined spheroid

Under the action of the external force F_e, the spheroid accelerates, but the drag force generated by the fluid opposes this motion. Balancing the two opposing forces, the particle undergoes a steady translation with a constant velocity. The variation in the translational velocity of unconfined prolate and oblate spheroids in an otherwise quiescent fluid as a function of aspect ratio b/a is shown in Fig. 2.

Fig. 2 Steady state translational velocity of an unconfined spheroid as a function of aspect ratio for a fixed particle volume. Aspect ratios b/a < 1 correspond to prolate spheroids, while b/a > 1 represent oblate shapes. The continuous curves show the analytical predictions (Section 2.1), the symbols denote the lattice Boltzmann simulation results, and the dashed lines indicate the asymptotic expressions for rod-like and disc-like limits. In the simulations, the spheroids were initially placed at the centre of the domain, with their symmetry axis either parallel (‖) or perpendicular (⊥) to the applied force (F_e = 0.05 lattice units). The spheroids retained this orientation through out the simulations. Re, calculated based on the steady state translational velocity of the spheroid, ≪1. For the sphere (b/a = 1), Re ≈ 0.023. Data plotted in blue (red) correspond to end-on (broadside-on) configurations, respectively. The aspect ratio at which the translational velocity is maximum in each case is marked by a vertical line.

The aspect ratio b/a < 1 for prolate and >1 for oblate spheroids. Two different cases are considered in Fig. 2: (i) the end-on translation, i.e., ê‖F_e, the spheroid translates with its axis of revolution oriented parallel to the direction of the external force and (ii) the broadside-on translation, i.e., ê⊥F_e, the spheroid translates with its axis of revolution oriented perpendicular to the direction of the external force. The corresponding velocities are denoted as U_‖ and U_⊥ respectively. The symbols are obtained from lattice Boltzmann simulations and the continuous lines are obtained from the analytical expression, eqn (3). The dashed lines on the extremes of both sides of the aspect ratio are analytical approximations – based on slender bodies for prolate and discs for oblate spheroids.

In typical applications driven by external fields such as magnetic or gravitational, the force density (force per unit volume) on the particle will be constant. Hence, in Fig. 2 the aspect ratio of the spheroid is varied by keeping the volume of the spheroid constant so that the total force acting on the translating particle is fixed. Furthermore, the translational velocities, U_‖ and U_⊥ are non-dimensionalised with , the Stokes velocity of an unconfined sphere of same volume V. Thus, U_‖/U_⊥ = U_⊥/U_sphere = 1 when b/a = 1.

The excellent match between the simulation data and the analytical solutions at all aspect ratios indicates the reliability of the numerical method used in this work. More importantly, Fig. 2 shows the existence of a maximum in the two cases considered. For the case of end-on translation (blue curve) the maximum in the translational velocity occurs at an aspect ratio b/a ≈ 0.512, that of a prolate shape. For the case of broadside-on translation (red curve) the maximum occurs at an aspect ratio of b/a ≈ 1.514, that of an oblate shape. The curves corresponding to two cases intersect at b/a = 1, that of the sphere, and the maxima in two cases occur on either side of the sphere.

The maximum in the translational velocity (or minimum in the drag force) has been observed for a variety of anisotropic particles in the context of gravitational sedimentation.^24,63–66 It exists because the viscous drag force acting on the translating particle has two contributions (i) friction drag and (ii) pressure drag, and these two contributions scale differently with change in the dimensions of the particle.^67–69 The skin friction arises from the shear stress on the particle's surface due to the viscosity of the fluid, and this contribution depends on the surface area of the particle. The pressure drag, also known as the form drag, arises as the particle pushes the fluid in front of it to move ahead. Thus, there is a pressure gradient across the translating particle, and the pressure drag contribution primarily depends on the projected area of the particle in the direction of translation.⁶⁸ The difference in the dependence of friction and pressure drag on the spheroid surface area results in the presence of maxima in Fig. 2, and can be rationalised as follows.

Consider a sphere of radius r translating under an external force. Both the total surface area, and the projected surface area scale as ∼r². Exact calculations show that the friction drag contributes 2/3rd and the pressure drag contributes 1/3rd of the total drag force.^67,68

To examine a small deviation from sphericity, let us consider a prolate spheroid with a > r. For end-on translation, the projected surface area is ∼b² while the total (or lateral) surface area is ∼ab. Under the constraint of fixed volume V ∼ ab², the projected area scales as ∼V/a and the total area as . Thus, changing the shape from a sphere to a prolate spheroid decreases the projected surface area as ∼1/a but increases the total surface area as . As a result, the reduction in pressure drag outweighs the increase in friction drag, leading to a net decrease in the total drag force; consequently, an end-on translating prolate spheroid attains a higher velocity than a sphere of equal volume. Exact computations confirm this scaling argument.^69,70

This trend however does not hold for highly elongated prolate spheroids (a ≫ b), where the lateral area far exceeds the projected area, i.e. ab ≫ b². In this regime, friction drag dominates, as seen in analytical results for needle-like bodies where pressure drag vanishes⁶⁷ and in numerical studies of spheroids.⁶⁹ Thus, for highly elongated particles, the total drag increases with a.

In other words, prolate spheroids that are closer to a sphere in shape (b/a ≲ 1) translate faster on increasing the length of the major axis, a, but a similar change of increasing a slows down highly elongated, rod-like spheroids (b/a ≪ 1). This contrast in the dependence of drag force at large (b/a ≲ 1) and small (b/a ≪ 1) aspect ratios gives rise to a maximum in the translational velocity at an intermediate value of aspect ratio for the translating prolate spheroid in the end-on configuration.

By contrast, when a sphere is deformed into an oblate spheroid translating in the end-on configuration, the increase in projected surface area (∼b²) exceeds the decrease in lateral area (∼ab), leading to a reduction in translational velocity relative to that of a sphere. This trend persists at larger aspect ratios (b/a ≫ 1), so that the translational velocity of an oblate spheroid decreases monotonically from the spherical limit (b/a ≳ 1). Similar arguments of anisotropic scaling of the relevant surface areas hold for broadside on translation as well, and a maximum velocity appears for an oblate spheroid in this configuration.

4.2 Confined spheroid: centreline of the channel

In this section, we investigate the influence of confinement on the translational velocity of a spheroid by placing it inside a channel with a square cross-section. As illustrated in Fig. 1a, the external force is applied along the open direction of the channel (the x-axis). For simplicity, the spheroid is positioned on the channel centreline in an end-on configuration, i.e., its axis ê is aligned with the channel axis. As in the previous section, we apply a constant body force density to the spheroid by keeping the volume of the spheroid and the applied body force F_e constant while changing the aspect ratio and the confinement ratio.

The effect of confinement is quantified using the confinement ratio, CR = 2b/L where L is the height (equal to the width) of the channel. As shown in Fig. 2a, the translational velocity of both prolate and oblate spheroids decreases with increasing confinement. Prolate spheroids exhibit lower translational velocities than oblate spheroids because they share a larger surface area with the channel walls and therefore experience greater hydrodynamic resistance.

Fig. 3a also includes translational-velocity data for spheroids (i) in a cylindrical channel and (ii) obtained from the approximate analytical method for comparison. Wakiya⁶² analysed the axisymmetric translation of spheroids in cylindrical channels using a combined analytical–numerical approach, and reported the translational velocity accurate up to in a tabulated form. Data extracted from his results for comparable aspect ratios—although available for only a few confinement ratios (CRs)—are shown in Fig. 3a as dash-dotted lines. The slightly larger translational velocities observed in the simulations arise from geometric differences; specifically, a square channel is marginally “less confining” than a cylindrical one due to the presence of corners. The results from the approximate analytical calculations, accurate up to and based on the superposition arguments discussed in Section 2.2, are also shown in Fig. 3a as dotted lines. These results systematically over-predict the translational velocity and exhibit only qualitative agreement with the simulations.

Fig. 3 (a) Translational velocity of spheroids (prolate, b/a ≈ 0.51 and oblate, b/a ≈ 1.57) moving in an end-on configuration along the centreline of a square channel, shown as a function of increasing confinement ratio, CR = 2b/L. The symbols ● are from the LBM simulations, the dash-dotted lines are spheroids in cylindrical channels (prolate, b/a = 0.5 and oblate b/a = 2)⁶² and dotted lines are theoretical predictions based on far-field analysis (Section 2.2). The ordinate is normalised with the translational velocity of the unconfined spheroid. Translational velocities of spheroids for (b) F‖ê and (c) F⊥ê configuration as a function of aspect ratio for fixed confinement ratios – ■: CR = 0, ▲: CR = 0.2, ♦: CR = 0.5 and ●: CR = 0.9. The ordinate is normalised with the translational velocity of an unconfined sphere and the optimum aspect ratio at each confinement ratio is indicated on top of each peak in (b) and (c).

We now examine how the translational velocity of confined spheroids varies with their aspect ratio. Fig. 3b shows the normalised translational velocity of spheroids moving along the channel centreline in the end-on configuration, with aspect ratios ranging from b/a < 1 (prolate) to b/a > 1 (oblate). As in the unconfined case, we find that for each confinement ratio there is an aspect ratio at which the translational velocity reaches a maximum. This optimum shifts to higher aspect ratios (toward the oblate side) as the degree of confinement increases.

The shift of the maximum translational velocity toward the oblate side arises from the nature of the dominant drag mechanisms. In this geometry, the primary contribution to hydrodynamic resistance comes from friction along the lateral surface of the particle, generated by the shear flow between the spheroid and the channel walls. For a given confinement ratio and in the end-on configuration, an oblate spheroid presents a smaller lateral surface area to the walls than a prolate one, leading to lower frictional drag, and therefore a higher translational velocity. However, as the aspect ratio becomes very large, the gap between the particle and the walls becomes quite small. This leads to strong shear and a substantial increase in frictional drag. Consequently, the total drag becomes very large as b/a ≫ 1. This implies that the translational velocity must peak at a finite aspect ratio, and that this optimum lies on the oblate side.

We also consider the broadside-on configuration of spheroids translating along the channel centreline. The corresponding results are presented in Fig. 3c. As in the previous cases, an optimal aspect ratio exists; it occurs for an oblate particle and shifts to larger values as the confinement increases. In this configuration, the location of the optimum is again determined by the frictional resistance between the particle and the channel walls, but for a different geometric reason. As the aspect ratio of an oblate spheroid increases, its flatter faces move farther from the walls while its thinner edges move closer. Consequently, the hydrodynamic resistance initially decreases at first because the dominant flat surfaces experience lower frictional forces, before increasing again at very large aspect ratios due to the progressively narrower gaps near the thin edges. The curves tend to flatten at very high aspect ratios, as the contribution of these thin edges to the overall frictional drag becomes relatively constant.

Overall, these results indicate that an optimal aspect ratio exists for maximising translational velocity in both confined and unconfined settings. For most configurations, this optimum lies on the oblate side, where oblate spheroids achieve higher translational velocities than spheres and prolate counterparts.

4.3 Confined spheroid: eccentric positions in the channel

In this section, we relax the assumption of symmetric positioning and orientation of the spheroid and examine its dynamics for all possible positions and orientations within the channel. We define h as the distance from wall 1 of the channel and θ as the orientation angle of the spheroid, as illustrated in Fig. 1a. To simplify the analysis, we restrict our investigation to the xy-plane at z = L/2. The translational velocity U_x generally depends on both h and θ. In addition, the spheroid acquires a transverse velocity U_y and an angular velocity Ω_z. The computed values of U_x, U_y and Ω_z as function of h and θ are presented as contour plots in the h–θ plane in Fig. 4. The simulation data along with theoretical predictions for U_x, U_y and Ω_z obtained from the far-field analysis are also shown as line plots in second row of Fig. 4.

Fig. 4 Contour plots of (a) translational velocity U_x, (b) transverse velocity U_y and (c) angular velocity Ω_z in the h–θ space for a channel confined spheroid of aspect ratio b/a ≈ 0.51 at a confinement ratio CR ≈ 0.15. The colour bars are shown at the top of each plot. h is normalised by the channel width L and θ is given in degrees (°). Ω_z > 0 corresponds to anti-clockwise rotation. In each figure, the inaccessible areas due to the finite size of the spheroid is coloured in grey. (d)–(f) Same data as line plots (★ symbols) along with results from the far field analysis (continuous lines) as a function of h at various θ. Results in (d) and (e) have been normalised with the translational velocity of the spheroid with same orientation (θ) in an unbounded fluid. Plots (a)–(c) correspond to a Re ≈ 0.25 where as plots (d)–(f) correspond to a Re ≈ 0.01.

Fig. 4a shows the translational velocity U_x of the spheroid along the direction of the applied force (along the x-axis). The translational velocity is maximal when h/L = 0.5 and θ = 0°, 180°, corresponding to the end-on orientation of the spheroid positioned along the channel centreline. At other positions and orientations within the channel the translational velocity decreases because (1) the pressure drag increases as the projected surface area grows with θ, and (2) the viscous drag increases due to hydrodynamic interactions with both the wall and the spheroid. The far-field predictions in Fig. 4d show qualitative agreement with the full numerical results; however, the quantitative agreement deteriorates near the channel walls.

Fig. 4b illustrates the variation in the transverse velocity U_y, i.e., the velocity of the spheroid along the y-axis. The transverse velocity relative to the reference wall can be either positive or negative, and it changes sign every 90° rotation of the spheroid, consistent with the symmetry of the system. The orientations at which the transverse velocity vanishes correspond to θ = 0°, 90°, 180°, and 270°. At these angles, the external force F̂ is either fully parallel or fully perpendicular to ê, and therefore the colloid does not migrate across the channel. In other words, hydrodynamic interactions with the wall in these configurations do not generate any transverse motion, in accordance with the reversibility constraints of Stokes flow. As in the previous case for U_x, the comparison between the full numerical results and the far-field predictions (Fig. 4e) is qualitative, with the agreement deteriorating near the walls.

The angular velocity of the spheroid, Ω_z in the h–θ phase space is shown in Fig. 4c. The rotation arises solely from the hydrodynamic torque on the particle, generated by the interplay between form drag, viscous drag, and the anisotropic shape of the spheroid. Depending on h and θ, the spheroid rotates clockwise (Ω_z < 0) or counterclockwise (Ω_z > 0) with the sign changing at the channel centreline. More importantly, the angular velocity also changes sign near the walls at θ = 45°, 135°, etc., indicating that the direction of rotation as the spheroid approaches a wall depends strongly on its incident orientation. This orientation-dependent reversal of rotation leads to the glancing–reversing trajectories of spheroids discussed in the next section. The far-field predictions of the angular velocity shown in Fig. 4f matches numerical simulations qualitatively, the asymmetry in the data in the numerical simulations is a finite Reynolds number effect (see Section 4.5).

4.4 Glancing-reversing trajectories of the confined spheroid

We now analyse the trajectories of the confined spheroid as it translates through the channel under the action of the external force. Two types of trajectories are observed—glancing and reversing—similar to those reported for motion near a single wall.³¹ These trajectories are illustrated schematically in Fig. 5a and b.

Fig. 5 Schematic representation of (a) glancing and (b) reversing spheroid. The dotted line represents the channel centreline. The instantaneous orientation, indicated by the arrow ê, is provided as a visual guide to illustrate the angular rotation of the spheroid. Glancing (reversing) trajectories span the full (half) channel width, with the spheroid approaching each wall almost parallel (perpendicular) to it. (c) and (d) Velocity vectors in the h–θ phase space for a spheroid of aspect ratio b/a ≈ 0.36 and b/a ≈ 0.7. h is normalised by the channel width L and θ is given in degrees (°). The background is coloured with the angular velocity; Ω_z > 0 corresponds to anti-clockwise rotation. Black vector arrows indicate the direction of the velocity of the spheroid in the phase space. The fixed points are highlighted with different markers (see text). In each figure, the inaccessible areas due to the finite size of the spheroid is coloured in grey.

In a glancing trajectory (Fig. 5a), the spheroid grazes the channel wall with its long axis oriented nearly parallel to the wall. Due to the repulsive hydrodynamic interaction (i.e., the transverse velocity induced by the wall), the spheroid begins to move away from the wall while undergoing anticlockwise rotation. This drift causes the spheroid to approach the opposite wall, where it is again repelled while rotating clockwise, eventually returning to the first wall. Thus, over one full cycle, the spheroid spans the entire channel width, with its orientation oscillating between angles less than and greater than 180°. In contrast, during a reversing trajectory (Fig. 5b), the spheroid remains near one wall and never crosses the channel centreline. In this case, the spheroid approaches the wall with its long axis nearly perpendicular to it, and hydrodynamic interactions cause its orientation to oscillate between angles less than and greater than 90° during each cycle. Thus, although both types of trajectories are oscillatory, the resulting motions are qualitatively different.

The origin of these two types of channel-confined spheroid trajectories can be understood by examining the corresponding phase-space dynamics. The velocity vectors in the h–θ phase space for spheroids of two different aspect ratios are shown in Fig. 5c and d. The phase space is coloured by the angular velocity Ω_z. As expected for a low–Reynolds-number system, the phase-space trajectories form closed loops, consistent with the oscillatory motion of the spheroid in the channel. With these phase portraits, we can identify the origin of the two distinct classes of trajectories.

The first type of trajectory are closed loops centred about h/L = 0.5 and θ = 180° (marked ● in the figures); this set corresponds to the glancing trajectories. Glancing trajectories of opposite sense of rotation are centred at 0° (or 360° equivalently). These trajectories occupy the entire width of the channel, thus crossing the channel centreline where the angular velocity changes sign. Glancing occurs when a spheroid approaching the bottom (or top) wall with its long axis nearly parallel to the wall undergoes anticlockwise (clockwise) rotation and subsequently drifts toward the opposite wall. There exist infinitely many glancing trajectories, each individual trajectory is uniquely determined by the initial position and orientation of the spheroid.

The second type, the reversing trajectories are sandwiched between the glancing trajectories. They are closed loops, centred around values of h/L close to either the top or bottom wall at θ = 90° (marked ★ in the figures). They occupy at most one half of the channel width. A pair of reversing trajectories with opposite sense of rotation exists symmetrically about the channel centreline (on either side of the symbol ■). An equivalent pair is centred about θ = 270°.

In general, the reversing trajectories are more skewed in phase space compared to the glancing trajectories. Importantly, these qualitative features are independent of the spheroid aspect ratio (see Fig. 5c and d). Thus, glancing and reversing trajectories are intrinsic characteristics of channel-confined, force-driven spheroidal particles.

The end-on configuration of a prolate spheroid translating along the centreline of the channel (h/L = 0.5, θ = 0°), discussed in Section 4.2, forms the centre of the glancing trajectories. This configuration corresponds to a neutrally stable point in the phase space. In addition to this, the phase space contains a pair of additional, neutrally stable centres—these are the centre points of the reversing trajectories—located on either side of the channel centreline and closer to the top and bottom walls. On the other hand, the broadside-on configuration at the channel centreline (h/L = 0.5, θ = 90°) is a saddle point formed between the pair of glancing trajectories (of opposite sense of rotation) and the pair of reversing trajectories (also of opposite sense of rotation), marked ■ in the figures. Consequently, the broadside-on state of a channel-confined spheroid at the centreline, is an unstable configuration; any small perturbation in position or orientation drives the spheroid away from this state and into either a glancing or a reversing trajectory.

4.5 Effect of fluid inertia and resulting bifurcations

In the results discussed so far, we have neglected the effects of fluid inertia and assumed the Stokes regime. We now relax this assumption and examine the influence of fluid inertia on the dynamics of a channel-confined spheroid. We fix the aspect ratio of the prolate spheroid at b/a ≈ 0.51 and vary the Reynolds number Re. The reported Re is based on the steady state velocity of the translating spheroid. Note that the quasi-static simulations cannot be used to study the effect of inertia; we obtain the complete trajectory of the spheroid using LBM simulations.

The phase-space trajectories in the h–θ plane for a channel-confined spheroid at Re ≈ 0.2 are shown in Fig. 6a. These trajectories are generated by initialising the spheroid at h ≈ 0.47, close to the channel centreline, and varying the initial orientations θ in the range 90°–180°. The initial and final points are marked with ★ and ◊ symbols respectively. Unlike the dynamics in the Stokes regime, the paths that the spheroid takes in phase-space are no longer closed loops. The glancing trajectories spiral outward and open up, and as a result, each glancing trajectory ultimately merges smoothly into a reversing trajectory. The reversing trajectories, in turn, spiral inward toward their centre point. Consequently, a spheroid that initially follows a glancing path will eventually transition into a reversing trajectory and then relax toward a broadside-on configuration (θ = 90°) near either the top or bottom wall.

Fig. 6 (a) Phase space trajectories of the spheroid at Re ≈ 0.2. The trajectories are obtained by initialising the spheroid at h/L ≈ 0.47 and θ varied in 90°–180°. Phase space trajectories of the spheroid at various Re (b) for initial orientation θ = 150° and (c) initial orientation θ = 170° for initial h/L ≈ 0.47. In each figure, the initial and final points are highlighted with ★ and ◊ respectively. The final point corresponds to the last data point obtained in the simulation lasting 6 × 10⁶ iteration steps, and may not necessarily correspond to a fixed point.

These qualitative changes arise from the introduction of fluid inertia. At finite Re, the flow no longer adjusts instantaneously to the particle configuration, and inertial forces develop that break the neutrally stable closed loops characteristic of Stokes flow. Glancing trajectories experience a net outward drift away from the centreline, while reversing trajectories experience an inward drift toward orientations near θ = 90°. As a result, the end-on configuration at the centre point, h/L = 0.5, θ = 180° loses neutral stability, and the system develops four stable broadside-on fixed points – two near the upper wall and two near the lower wall. The saddle node at h/L = 0.5, θ = 90° remains an unstable fixed point at this Reynolds number (Re ≈ 0.2). However, the phase-space structure undergoes significant changes as Re is increased further.

To systematically examine the influence of fluid inertia, we consider two initial orientations, θ = 150° and θ = 170° at h/L ≈ 0.47, and vary the Reynolds number Re. The resulting trajectories in the h–θ phase space are shown in Fig. 6b and c respectively. The dynamics exhibit pronounced nonlinearity: the outward spiralling characteristic of glancing trajectories and the inward spiralling associated with reversing trajectories are visible only at relatively small Re. As Re increases beyond , additional fixed points emerge, and the trajectories bear little resemblance to the aforementioned glancing and reversing paths observed in the Stokes regime. At the highest Reynolds number considered in this study (Re ≈ 5), the spheroid rapidly settles into a stable broadside-on configuration at the channel centreline (h/L = 0.5, θ = 90°), similar to that of an unconfined spheroid.^32–36 Thus, even modest levels of inertia can fundamentally alter the topology of the phase space and promote stable broadside-on alignment.

To map the dependence of the stable fixed points on the Reynolds number, we construct the bifurcation diagram shown in Fig. 7. The plot corresponds to a prolate spheroid of aspect ratio b/a ≈ 0.51, and in this diagram the final equilibrium position, denoted by (h_f/L) is plotted as a function of the Reynolds number Re. The diagram is obtained by initialising the spheroid at various positions (h/L) across the channel width and allowing it to evolve to its steady state. For a single stable solution exists at h_f/L = 0.5, corresponding to broadside-on translation at the channel centreline. In contrast, for small Reynolds numbers , this centreline state is unstable, and stable broadside-on motion occurs near the channel walls. This transition indicates the presence of a bifurcation at . A detailed characterisation of the mathematical nature of this bifurcation requires further analytical and numerical investigation and is beyond the scope of the present work.

Fig. 7 The bifurcation diagram shows effect of Re on the final position along the channel width (h_f/L). The stable (unstable) points are indicated by filled (open) symbols.

The presence of spiralling trajectories in the phase space indicates that fluid inertia also breaks the symmetry in the translational and angular kinematics of the spheroid with respect to the channel centreline. The extent of asymmetry is proportional to the Reynolds number Re. We find that both components of translational velocity and the angular velocity, U_x, U_y and Ω_z are visibly not symmetric with respect to the channel centreline when quasi-static simulations are performed at . With regard to the final configuration of the particle, previous studies have shown that fluid-inertia-induced torque on non-spherical particles cannot be neglected even at Reynolds numbers much smaller than unity.^36,71 The preferred orientation of a translating spheroid due to finite inertia is at θ = 90°, i.e., the broadside-on configuration.^{2,33–36,44,72} Our results show that, even in confinement this result holds; but the location of the equilibrium positions can be either at the channel centreline or near the walls depending upon Re as shown in Fig. 7. Hence, the approximate results obtained by Mitchell and Spagnolie³¹ which account for the wall induced hydrodynamic interaction can be substantially altered due to the presence of a small but finite Reynolds number.

5 Conclusions

Anisotropic particles appear widely in soft-matter systems and complex fluids, and even in simple Newtonian fluids they exhibit remarkably rich dynamics. In this work, we used lattice Boltzmann simulations to examine the motion of a force-driven spheroid and complemented these simulations with an analytical approach based on the superposition principle and far-field hydrodynamics. By systematically varying the particle shape, degree of confinement, and fluid inertia, we demonstrated how each of these factors influences the particle dynamics.

In unconfined conditions, we found that the shape corresponding to the maximum translational velocity for a fixed particle volume is not a sphere. Instead, the highest velocity occurs for a prolate spheroid of aspect ratio b/a ≈ 0.512 in an end-on configuration and for an oblate spheroid of aspect ratio b/a ≈ 1.514 in a broadside-on configuration. A similar trend persists in confined geometries, where the optimal shape typically corresponds to oblate spheroids. When a spheroid is placed at an off-centre position or orientation within the channel, transverse and angular velocities are generated, and the coupling between translation and rotation produces oscillatory trajectories. Two distinct types of oscillatory dynamics—glancing and reversing—were identified. In a glancing trajectory, the spheroid oscillates between the two channel walls with its long axis nearly parallel to the wall near each encounter. In contrast, reversing trajectories are confined to one half of the channel width, with the spheroid approaching the wall in a broadside-on configuration.

These oscillatory motions appear as initial-condition–dependent closed loops in phase space under Stokes flow. However, the inclusion of fluid inertia breaks these closed loops. At small but finite Reynolds numbers, glancing trajectories spiral outward and merge into reversing trajectories, which spiral inward toward a fixed point corresponding to a broadside-on orientation near the walls. As the Reynolds number increases further, the phase-space structure undergoes dramatic changes: new fixed points emerge, existing ones change stability, and the overall dynamical landscape becomes more intricate. We have characterised these transitions using a bifurcation diagram. While our study is primarily computational, it highlights the need for future investigations from a nonlinear dynamical perspective to fully unravel the complexity of this system.

In summary, our work reveals the rich and intricate dynamics associated with the classical problem of a force-driven spheroid in a confined channel. Further analytical and numerical studies are needed to fully understand the nonlinear behaviour uncovered here. Nevertheless, our systematic analysis offers predictions that can be readily tested in microfluidic experiments and provides guidance for optimising the shape of drug-delivery particles, illustrating the broad possibilities that nonspherical particles bring to applications in transport and design.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Data is generated through the open-source program LUDWIG (https://github.com/ludwig-cf/ludwig) and is available upon request.

Acknowledgements

Computational facilities at IIT Madras are duly acknowledged. SPT acknowledges the project support by I-Hub Foundation for Cobotics (IHFC), IITD and the travel support by IoE, IIT Madras. This work used the ARCHIE-WeSt High-Performance Computer (https://www.archie-west.ac.uk/) based at the University of Strathclyde. For the purpose of complying with UKRI's open access policy, the authors have applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising from this submission.

Notes and references

1. Y. Han, A. M. Alsayed, M. Nobili, J. Zhang, T. C. Lubensky and A. G. Yodh, Science, 2006, 314, 626–630.
2. N. Fares, M. Lavaud, Z. Zhang, A. Jha, Y. Amarouchene and T. Salez, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2411956121.
3. Y. Avni, S. Komura and D. Andelman, Phys. Rev. E, 2021, 103, 042607.
4. G. Boniello, C. Blanc, D. Fedorenko, M. Medfai, N. B. Mbarek, M. In, M. Gross, A. Stocco and M. Nobili, Nat. Mater., 2015, 14, 908–911.
5. M. Molaei, N. G. Chisholm, J. Deng, J. C. Crocker and K. J. Stebe, Phys. Rev. Lett., 2021, 126, 228003.
6. H. Zhou, C. C. Mayorga-Martinez, S. Pané, L. Zhang and M. Pumera, Chem. Rev., 2021, 121, 4999–5041.
7. V. A. Baulin, A. Giacometti, D. A. Fedosov, S. Ebbens, N. R. Varela-Rosales, N. Feliu, M. Chowdhury, M. Hu, R. Füchslin and M. Dijkstra, et al., Soft Matter, 2025, 21, 4129–4145.
8. G. Gompper, H. A. Stone, C. Kurzthaler, D. Saintillan, F. Peruani, D. A. Fedosov, T. Auth, C. Cottin-Bizonne, C. Ybert and E. Clément, et al., J. Phys.: Condens. Matter, 2025, 37, 143501.
9. T. Patiño Padial, S. Chen, A. C. Hortelão, A. Sen and S. Sánchez, Nat. Rev. Mater., 2025, 1–17.
10. X. Ju, C. Chen, C. M. Oral, S. Sevim, R. Golestanian, M. Sun, N. Bouzari, X. Lin, M. Urso and J. S. Nam, et al., ACS Nano, 2025, 19, 24174–24334.
11. T. He, Y. Yang and X.-B. Chen, Nanoscale, 2024, 16, 12696–12734.
12. T. Hueckel, G. M. Hocky and S. Sacanna, Nat. Rev. Mater., 2021, 6, 1053–1069.
13. J. Chen, N. E. Clay, N.-H. Park and H. Kong, Chem. Eng. Sci., 2015, 125, 20–24.
14. V. Lotito and T. Zambelli, Adv. Colloid Interface Sci., 2022, 304, 102642.
15. X. Zhu, C. Vo, M. Taylor and B. R. Smith, Mater. Horiz., 2019, 6, 1094–1121.
16. E. Ben-Akiva, J. W. Hickey, R. A. Meyer, A. Isser, S. R. Shannon, N. K. Livingston, K. R. Rhodes, A. K. Kosmides, T. R. Warren and S. Y. Tzeng, et al., Acta Biomater., 2023, 160, 187–197.
17. N. P. Truong, M. R. Whittaker, C. W. Mak and T. P. Davis, Expert Opin. Drug Delivery, 2015, 12, 129–142.
18. D. Esporrn-Ubieto, N. Ruiz-González, V. Di Carlo, D. Sánchez-deAlcázar, F. Lezcano, A. Pushkareva Fazullina and S. Sánchez, Adv. Funct. Mater., 2025, e10203.
19. X. Chen, B. Ran, H. Wu and G. Zeng, Adv. Funct. Mater., 2025, e04234.
20. M. Caraglio, H. Kaur, L. J. Fiderer, A. López-Incera, H. J. Briegel, T. Franosch and G. Muñoz-Gil, Soft Matter, 2024, 20, 2008–2016.
21. Z. Zhang, A. Klingner, S. Misra and I. S. Khalil, Sci. Rep., 2025, 15, 31041.
22. A. Oberbeck, J. Reine Angew. Math., 1876, 81, 62–80.
23. G. B. Jeffery, Proc. R. Soc. London, Ser. A, 1922, 102, 161–179.
24. J. Happel and H. Brenner, Low Reynolds number hydrodynamics: with special applications to particulate media, Springer Science & Business Media, 2012, vol. 1.
25. G. K. Batchelor, An introduction to fluid dynamics, Cambridge University Press, 1967.
26. K. Nissanka, X. Ma and J. C. Burton, J. Fluid Mech., 2023, 956, A28.
27. J. Gong, V. A. Shaik and G. J. Elfring, J. Fluid Mech., 2024, 984, A26.
28. A. Sharma, P. A. Bosler, R. Govindarajan and D. L. Koch, J. Fluid Mech., 2025, 1007, A53.
29. E. Guazzelli and J. F. Morris, A physical introduction to suspension dynamics, Cambridge University Press, 2011, vol. 45.
30. W. Russel, E. Hinch, L. G. Leal and G. Tieffenbruck, J. Fluid Mech., 1977, 83, 273–287.
31. W. H. Mitchell and S. E. Spagnolie, J. Fluid Mech., 2015, 772, 600–629.
32. R. Khayat and R. Cox, J. Fluid Mech., 1989, 209, 435–462.
33. P. Huang, H. H. Hu and D. D. Joseph, J. Fluid Mech., 1998, 362, 297–326.
34. T.-W. Pan, R. Glowinski and G. P. Galdi, J. Comput. Appl. Math., 2002, 149, 71–82.
35. J. Feng, H. H. Hu and D. D. Joseph, J. Fluid Mech., 1994, 261, 95–134.
36. V. Dabade, N. K. Marath and G. Subramanian, J. Fluid Mech., 2015, 778, 133–188.
37. H. Huang, X. Yang and X.-Y. Lu, Phys. Fluids, 2014, 26, 053302.
38. X. Yang, H. Huang and X. Lu, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 063009.
39. Z. Xia, K. W. Connington, S. Rapaka, P. Yue, J. J. Feng and S. Chen, J. Fluid Mech., 2009, 625, 249–272.
40. T. Swaminathan, K. Mukundakrishnan and H. H. Hu, J. Fluid Mech., 2006, 551, 357–385.
41. F. Fonseca and H. Herrmann, Phys. A, 2005, 345, 341–355.
42. T.-W. Pan, D. D. Joseph and R. Glowinski, Comput. Struct., 2005, 83, 463–478.
43. I. L. Claeys and J. F. Brady, J. Fluid Mech., 1993, 251, 411–442.
44. M. N. Ardekani, P. Costa, W. P. Breugem and L. Brandt, Int. J. Multiphase Flow, 2016, 87, 16–34.
45. S. P. Thampi, K. Stratford and O. Henrich, Phys. Rev. E, 2024, 109, 065302.
46. A. T. Chwang and T. Y. Wu, J. Fluid Mech., 1976, 75, 677–689.
47. L. G. Leal, Advanced transport phenomena: fluid mechanics and convective transport processes, Cambridge University Press, 2007, vol. 7.
48. X. Jiang, C. Xu and L. Zhao, J. Fluid Mech., 2024, 984, A40.
49. B. N. Radhakrishnan, A. Purushothaman, R. Dey and S. P. Thampi, Phys. Rev. Fluids, 2024, 9, 083302.
50. J.-C. Desplat, I. Pagonabarraga and P. Bladon, Comput. Phys. Commun., 2001, 134, 273–290.
51. Ludwig GitHub repository, https://github.com/ludwig-cf/ludwig.
52. T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer, 2017.
53. P. L. Bhatnagar, E. P. Gross and M. Krook, Phys. Rev., 1954, 94, 511–525.
54. R. Adhikari, K. Stratford, M. E. Cates and A. J. Wagner, Europhys. Lett., 2005, 71, 473–479.
55. G. Strang, Introduction to linear algebra, SIAM, 2022.
56. N.-Q. Nguyen and A. Ladd, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 046708.
57. F. Zhao and B. Van Wachem, Acta Mech., 2013, 224, 2331–2358.
58. F. Zhao and B. Van Wachem, Acta Mech., 2013, 224, 3091–3109.
59. B. van Wachem, M. Zastawny, F. Zhao and G. Mallouppas, Int. J. Multiphase Flow, 2015, 68, 80–92.
60. G. A. Voth and A. Soldati, Annu. Rev. Fluid Mech., 2017, 49, 249–276.
61. F.-G. Fan and G. Ahmadi, J. Aerosol Sci., 1995, 26, 813–840.
62. S. Wakiya, J. Phys. Soc. Jpn., 1957, 12, 1130–1141.
63. G. R. Carmichael, Ind. Eng. Chem. Process Des. Dev., 1982, 21, 401–403.
64. M. Pourali, N. O. Jaensson and M. Kröger, J. Fluid Mech., 2021, 924, A30.
65. E. Loth, Powder Technol., 2008, 182, 342–353.
66. S. K. Sanjeevi, J. F. Dietiker and J. T. Padding, Chem. Eng. J., 2022, 444, 136325.
67. T. Aoi, J. Phys. Soc. Jpn., 1955, 10, 119–129.
68. D. Leith, Aerosol Sci. Technol., 1987, 6, 153–161.
69. R. Ouchene, Phys. Fluids, 2024, 36, 053603.
70. R. Ouchene, M. Khalij, A. Tanière and B. Arcen, Comput. Fluids, 2015, 113, 53–64.
71. M. Z. Sheikh, K. Gustavsson, D. Lopez, E. Lévêque, B. Mehlig, A. Pumir and A. Naso, J. Fluid Mech., 2020, 886, A9.
72. X. Jiang, C. Xu and L. Zhao, J. Fluid Mech., 2024, 1000, A49.

---