Three Simple Stokeslet Trajectories¹¹1Dedicated to William B. Russel, and especially his early work on suspension hydrodynamics

First we consider three identical particles in a vertical line, with center-to-center displacements parallel to gravity, which exerts identical forces to the particles. Let the positions of the particles be $x_{1}$, $x_{2}$, and $x_{3}$, counting in the direction of decreasing gravitational potential (i.e., towards the floor). Assume that the particles are in an infinitely-extended fluid. For simplicity here and in what follows, we work in distance units of particle radii $a$ and time units of Stokes times $6\pi\mu a^{2}/|\mathbf{f}|$, where $\mathbf{f}=4\pi a^{3}\Delta\rho\,\mathbf{g}/3$ is the buoyant force. The latter is the time it takes for a lone particle to translate a distance unit. This system of units absorbs the force magnitude into the time scale and uses force vectors of unit magnitude in the equations. Under these circumstances, the self-mobility of each particle $i$ is $\mathbf{M}_{ii}=\mathbf{I}$, and the pair mobility of particles $i$ and $j$ uses the form of the Stokeslet in Eq. 1.

The equation of motion of this particle system in its single spatial dimension is the following.

Because the spatial coordinates of interest are in the same direction as the force, the unit force is omitted in the equation above. Lone particles translate at a velocity of one, but their velocity increases above that due to hydrodynamic effects of the other two particles. It is simpler to analyze the center-to-center displacements instead of the centers themselves, since the system is translationally invariant. Let $x_{23}=x_{3}-x_{2}$ and $x_{12}=x_{2}-x_{1}$.

This pair of ordinary differential equations (ODEs) can be rewritten as a single differential equation in $x_{12}$ and $x_{23}$, omitting $t$, by the standard trick of “solving for $\mathrm{d}t$” in each equation and equating the resultant expressions in each.

After some algebra, we obtain a separable first-order ODE, which we can solve, revealing a constant of integration, a conserved parameter in the dynamical system of three particles.

It is proportional to the interparticle contribution to particle 2’s velocity in Eq. 4, rendering that particle’s velocity constant for all time.

The constancy of the velocity of the central particle is a non-trivial result. It is not merely the result of symmetry arguments, but it depends on the particular decay power of the long-range hydrodynamic interactions. One may verify this by considering particle 2’s velocity with an arbitrary decay power $\alpha$ for the interparticle mobility and checking if it is a constant of integration.

Taking its first temporal derivative, and combining terms, we arrive at the following.

Above, we introduced $k_{\alpha}$, a coefficient that weights the reciprocal powers in the interparticle mobility, which we may take to depend upon the power of the decay of the hydrodynamic interaction. For $\alpha=1$ we have $k_{1}=3/2$. Combining the two terms, we arrive at the following.

One may verify that the right-hand side vanishes for $\alpha=1$. Instead, for $\alpha=2$, we have the following.

The numerator is nonzero, so the quantity proportional to particle 2’s velocity, and therefore the velocity itself, is not conserved for other power-law decays in pair interaction.

Leichtberg et al. studied this three-particle problem with more accurate hydrodynamics through the discrete point boundary method ¹¹. They noted the following about the near constancy of the velocity of particle 2 (emphasis mine).

After the initial unsteady period (which is shrunk to zero at $\mathrm{Re}_{\infty}=0$, […]), spheres 1 and 2 possess essentially the same velocities, which are 30–40% greater than the velocity of sphere 3. The velocity of sphere 2 does not vary greatly from this point to the end of the experiment.

With the more realistic hydrodynamics used by these authors, the velocity of particle 2 is not conserved, but the results in the present work finally give rationale for the near constancy they noticed. It is due to the leading-order hydrodynamics and specifically the inverse-distance decay.

Second we consider the sedimentation of a pair of identical particles arranged horizontally above a flat, rigid wall. The Stokeslet pair mobility in this geometry can be constructed by a system of singularities at mirror-image points within the wall ^{2, 12}. It consists of a inifinitely-extended-fluid Stokeslet of Eq. 1 and several “reflections” across the wall’s surface: another Stokeslet, its gradient, and its Laplacian, with the latter two weighted by the distance from the forced particle to the wall. The velocity fields used to construct it enforce the no-slip condition at the surface: zero velocity for an assumed stationary wall. See Fig. 1 (c) and (d) for level sets of this sum.

We introduced a unit normal that points away from the wall and into the fluid $\mathbf{e}_{\text{w}}$, a tensor $\mathbf{P}\equiv\mathbf{I}-2\mathbf{e}_{\text{w}}\mathbf{e}_{\text{w}}$ that reflects across the wall, and a primed symbol $\mathbf{x}_{j}^{\prime}\equiv\mathbf{P}\cdot\mathbf{x}_{j}$ which is shorthand for the coordinate reflected from $\mathbf{x}_{j}$. As written, “target” particle index $i$ is the particle receiving the velocity disturbance produced by force on “source” particle $j$. Gradients are taken with respect to the reflected source particle’s position $\mathbf{x}_{j}^{\prime}$. The corresponding self mobility replaces the first term with an identity tensor.

When two identical particles with coordinates $\mathbf{x}_{1}(t)$ and $\mathbf{x}_{2}(t)$, with center-to-center separation vector $\mathbf{x}_{2}-\mathbf{x}_{1}$ perpendicular to $\mathbf{e}_{\text{w}}$, sediment towards a wall (i.e., gravity is parallel to $\mathbf{e}_{w}$), their distance from the wall $y(t)$ and their center-to-center separation distance $x(t)$ change simultaneously. The “un-reflected” infinitely-extended-fluid part of the Stokeslet interaction drives the two particles to the wall, it reducing $y$, but it does not affect $x$. The “reflections” embedded within the wall do, however, influence both $x$ and $y$, slowing the descent and pushing the particles apart. Our objective is to consider the extent of the latter effect. How far apart will particles spread apart from an initial center-to-center separation $L$?

For convenience, we consider the problem of sedimenting towards the wall in reverse, with particles falling from a wall instead, but the results apply when sedimenting towards the wall as long as certain quantities are negated. This means that particles are driven together rather than pushed apart. The linearity of the velocity with respect to the force means that our conclusions apply in the case of sedimentation towards the wall, as long as we suitably negate quantities.

After some algebra, the equations of motion for $x$ and $y$ are the following.

As-is, the set of equations is still difficult to analyze mathematically, and we seek to eliminate sub-dominant terms in the second equation to make things easier. To make progress, we assume two conditions. First, $x(t)$ does not change much from its large-enough initial value $L$. Second, $1/y$ can be considered small enough compared to $1$, the leading term in the second equation. We will justify our assumptions after our analysis by showing that our conclusions are consistent with them. If $L$ is indeed large enough, the $1/x$ and $1/\sqrt{x^{2}+y^{2}}$ terms are clearly $\ll 1$ for all $y$. So are the last two terms. At intermediate and large $y$, relative to $L$, the terms are comparable to reciprocal powers of $L$, which are also small by assumption. Thus, within our assumptions, at leading order we only need to retain the constant term in the second equation.

In other words, at leading order, the pair sediments towards the wall at the bare Stokes velocity, but hydrodynamic interaction from the wall still influences the center-to-center separation. Again, we can solve for $\mathrm{d}t$ in both equations, leaving us with a single equation for the horizontal separation, a non-separable first-order ODE.

Now, we utilize our assumption that $|x(t)-L|\ll L$ and linearize $x$ about $x=L$, introducing a variable $\Delta\equiv x-L$ and retaining only leading terms. This lets us integrate with respect to $y$ to obtain an estimate for $\Delta(y)$ in terms of $y$.

Above, we introduced the variables initial value $y_{0}$ to stand for the $y$ coordinate where $\Delta(y_{0})=0$, $u\equiv y/L$ a scaled height, and $u_{0}\equiv y_{0}/L$ its initial value. At $u\to\infty$, we have the “final” value of $\Delta$.

For convenience, we may select a $y_{0}$ large enough such that $1/y_{0}$ is smaller than some small error,³³3Choosing a large $y_{0}$ also avoids a known problem with the Stokeslet mobility near a wall, which causes particles very close to the wall to travel in the opposite direction of the force. but then we may also choose an $L\gg y_{0}$ larger than this, such that the $u_{0}$ terms contribute nothing to the limit above, which holds in the limit of large $L$. This justifies omission of the term proportional to $1/y$ in the dynamical equation for $y$. By this procedure, we obtain our final result, which is clearly $\ll L$, consistent with our assumptions.

This result, that the shrinking in horizontal separation of a widely-separated particle pair is both finite and largely independent of that separation, when given enough distance to travel away from the wall, is striking. One might instead expect the change in separation $|\Delta|$ to diminish as the initial separation $L$ of the particles increases. Despite this strangeness, the effect has a simple origin. The integral performed in Eq. 17 is “unitless.” In other words, the integral is independent of the choice of units for length, but the left-hand side of the equation is a ratio of a center-to-center separation and the unit $a$. Deeper than this, distances in the Green’s function for Stokes flow near a solid boundary are coupled together in groups with units of inverse length, such that when they are integrated over distance become unitless. Again, if the hydrodynamic interaction had a different spatial dependence, this result would not hold.

Now, we consider what happens when particles instead sediment towards rather than away from the wall. Due to the linearity of Stokes flow, negating the forces is equivalent to rewinding the trajectory we just analyzed. Therefore, the change in gap is just negated, and particles separate by three quarters of a radius.

It is worth considering how much headroom above a rigid wall a horizontal pair requires to realize most of this three quarters of a radius. To achieve 99% of this, assuming $u_{0}\ll 1$, we can solve $0.01=(1+6u^{2})/(1+4u^{2})^{3/2}\approx 3/(4u)$ for $u$, and we find that $u\approx 75$, or a height of $75L$.⁴⁴4Lengths much larger than a particle radius must be compared against length scales where fluid inertia becomes as important as viscosity. We may estimate the latter by $a/Re$, where $Re$ is the Reynolds number ⁵.

Finally, it must be mentioned that an infinite fluid above an infinitely-extending surface is an idealization. It is also atypical that a pair of particles be perfectly identical, isolated from other particles, and arranged horizontally. However, if the first two conditions can be guaranteed, deviation in horizontal orientation might not qualitatively change the result. Generalizing to larger particle counts in a line or a plane is not straightforward, as particles in the center of the arrangement would sediment faster than particles on the perimeter, ruining the linear or planar arrangement. If through an external potential the particles could be confined to a line or plane as they sediment, the configuration of the particles as they reach the wall might be interesting, especially if the line or plane of particles is initially disordered. The balance between infinite range (independence of $L$) and bounded influence (finite $\Delta$) may produce interesting effects. It should also be mentioned that in a real container that has walls on top and bottom, a pair of particles may contract together and then spread apart again as they fall from near the top wall to the bottom wall.

Third, we consider the sedimentation of another pair of identical particles, this time not restricted to a horizontal arrangement, through a specific ambient quadratic flow. This is a quadratic flow where the flow’s velocity is in the direction of gravity and its gradient is perpendicular to it.

Using the notation of Nadim and Stone ¹³, we write the equation of motion as follows, but absorbing constant factors into the flow curvature $\mathbf{K}$ for simplicity.⁵⁵5The flow curvature $\mathbf{K}$ here, in units convenient to sedimentation, is equal to $6\pi a^{3}\mathbf{K}/|\mathbf{f}|$ when instead using the $\mathbf{K}$ of Nadim and Stone.

The unit vector $\mathbf{e}_{\text{f}}$ points in the direction of the force, and the force magnitude is absorbed into the units. We may introduce another unit vector $\mathbf{e}_{\text{h}}$ that points “horizontally” and perpendicular to the force. For the problem at hand, the flow curvature $\mathbf{K}$ can be specified to $K\mathbf{e}_{\text{h}}\mathbf{e}_{\text{h}}\mathbf{e}_{\text{f}}$. Neither the Stokesian drift of the particles nor the particular flow curvature we consider changes the center-to-center displacement $\mathbf{x}_{12}$ in the direction of $\mathbf{e}_{\text{h}}$. Inserting the definition of the curvature and performing some algebra, we arrive at an equation of motion for the horizontal center $x_{\text{c}}\equiv(\mathbf{x}_{1}+\mathbf{x}_{2})\cdot\mathbf{e}_{\text{h}}/2$ and the vertical displacement $\Delta y\equiv\mathbf{x}_{12}\cdot\mathbf{e}_{\text{f}}$, written in terms of the flow curvature scalar $K$ and the constant horizontal displacement $\Delta x$.

We can also calculate the rate of change of the vertical center $y_{\text{c}}\equiv(\mathbf{x}_{1}+\mathbf{x}_{2})\cdot\mathbf{e}_{\text{f}}/2$, but it is extraneous and can be expressed using the closed set of equations above.

As we did previously, we can solve for $\mathrm{d}t$ in both parts of Eq. 21 and equate the results. The result is a constant of motion $H$ of the two-particle-in-quadratic-flow system.

The symbol for this constant is chosen because the dynamical system has a Hamiltonian formulation

where rates of change in “position” $q=\Delta y(t)$ and “momentum” $p=x_{\text{c}}(t)$ are written as suitable derivatives with respect to this Hamiltonian $H$.

The “position,” “momentum,” and “energy” in this formulation do not have units one would expect from mechanics. In this formulation, the “momentum” has units of length, and the “energy” has units of squared length per unit time. Hamiltonian formulations of suspension flows have appeared previously and have been used to explain the stability of horizontal polygonal arrangements of particles sedimenting through a quiescent fluid ¹⁴. Such formulations are exceptional in that they may be used to prove boundedness of phase-space trajectories and establish periodic solutions, as in the case of a harmonic oscillator. More on this later.

It is worth considering the linearized dynamical system, expanding in an assumed small parameter $\Delta y/\Delta x$, as this may reveal oscillation frequencies for certain values of the input parameters $x_{\text{c}}$, $\Delta x(0)$, $\Delta y(0)$, and $K$.

The linearized system can be written in matrix form.

It has a constant solution $x_{\text{c}}=\Delta y=0$, which is also true for the full nonlinear system. For $K=0$, the system has only a zero eigenvalue. Otherwise, it has a pair of eigenvalues with opposite signs.

If $K/\Delta x>0$, the eigenvalues are real, and the positive eigenvalue causes exponential departure from the neighborhood of the $x_{\text{c}}=\Delta y=0$ solution. If $K/\Delta x<0$, the eigenvalues are imaginary, permitting small-amplitude oscillations of period $2\pi/|\lambda|$ about the $x_{\text{c}}=\Delta y=0$ solution. The former, unstable case at is comparable to fluid in a vertical pipe being pumped upward while particles sediment downward. The latter, stable case at has fluid pumped in the direction of gravity.

Without loss of generality, we assume that $\Delta x>0$ from this point onwards. This makes $K>0$ map to unstable and $K<0$ to stable linearized scenarios. It also sets the ranges of the Hamiltonian in the two cases. For $K>0$, $H\in[0,+\infty)$, with the low end of the range pertaining to $x_{\text{c}}=0$ and $\Delta y=\pm\infty$. For $K<0$, $H\in(-\infty,3/4]$, with the high end of the range pertaining to $x_{\text{c}}=0$ and $\Delta y=0$.

Now we consider full nonlinear solutions, i.e., avoiding the linearization of Eq. 26. Example trajectories are sketched in Fig. 3. For $K<0$, like in the linearization, both $x_{\text{c}}$ and $\Delta y$ are bounded. The bounds are provided by extremizing $x_{\text{c}}(\Delta y)$ or $\Delta y(x_{\text{c}})$ at constant $H$ (i.e., along a trajectory) in Eq. 23. $x_{\text{c}}^{2}$ is maximized when $\Delta y=0$, when it takes a value of $[(3/4)-H]/(-K\Delta x)$. $\Delta y^{2}$ is maximized when $x_{\text{c}}=0$, where it takes a value of $\Delta x^{2}[9/(16H^{2})-1]$. Thus, the variables are always bounded along a trajectory. The symmetry of the phase portrait on negating $x_{\text{c}}$ and $\Delta y$ establishes the closedness of the trajectories in phase space, so the solutions for $K<0$ are periodic. Periodic $K<0$ and non-periodic $K>0$ trajectories in phase space are plotted in Fig. 4.

The oscillations found in this problem bear resemblance to those found on sedimenting prolate ellipsoids in vertical pipes ^{15, 16}, with obvious analogy between the Stokeslet pair and a prolate ellipsoid. In those studies, the oscillations are created by mobility gradients and other hydrodynamic effects near the pipe walls. It would be interesting to study the current oscillations in the presence of narrow tubes, and with more realistic hydrodynamics. For the former, the image system introduced by Liron may be a good starting point ¹⁷. For the latter, the approximation scheme used by Haber and Brenner, which models the quadratic flow in the vicinity of a solid particle as a velocity gradient, may be helpful ¹⁸.

Such quadratic flows may be produced naturally, rather than being forced by flow through a pipe. Suspension jets, or concentrated vertical assemblies of heavy particles, create an ensemble-averaged flow which is similar to the fluid-pumped-downwards stable case mentioned above ^{19, 20}. The suspension jets eventually break up due to a varicose instability into spherical suspension droplets, but it is worth considering the motion of pairs of particles inside those jets. Crosby and Lister provides an axial velocity field $v_{z}(r)$ for a theoretical suspension jet in cylindrical coordinates, written in the frame where there is no net flux in the interior of the jet $r\leq 1$.

Velocity components in $\theta$ and $r$ directions are assumed to be zero. The positive $z$ direction is the direction of sedimentation. The second derivative of this velocity field in the radial direction is relevant for a particle pair with a component of center-to-center separation across the diameter of the jet, as the flow curvature it generates can cause the particle pair to oscillate back and forth rather than depart the jet.

This analysis assumes that the other particles in the jet act in a smeared-out fashion, rather than disturbing the particles in the pair, an assumption that should be interrogated. The flow curvature then acts to localize the pair within the jet and eject the pair further out if they depart the jet, reflecting an affinity of particle pairs to the bulk of the suspension.

The author thanks Roseanna N. Zia for introducing him to suspension hydrodynamics and for bringing related work to his attention. He also thanks Christopher W. Macosko for helpful suggestions on figures.