Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Ryoichi
Yamamoto
*^{a},
John J.
Molina
^{a} and
Yasuya
Nakayama
^{b}
^{a}Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan. E-mail: ryoichi@cheme.kyoto-u.ac.jp
^{b}Department of Chemical Engineering, Kyushu University, Fukuoka 819-0395, Japan

Received
16th December 2020
, Accepted 9th March 2021

First published on 17th March 2021

A general method is presented for computing the motions of hydrodynamically interacting particles in various kinds of host fluids for arbitrary Reynolds numbers. The method follows the standard procedure for performing direct numerical simulations (DNS) of particulate systems, where the Navier–Stokes equation must be solved consistently with the motion of the rigid particles, which defines the temporal boundary conditions to be satisfied by the Navier–Stokes equation. The smoothed profile (SP) method provides an efficient numerical scheme for coupling the continuum fluid mechanics with the dispersed moving particles, which are allowed to have arbitrary shapes. In this method, the sharp boundaries between solid particles and the host fluid are replaced with a smeared out thin shell (interfacial) region, which can be accurately resolved on a fixed Cartesian grid utilizing a SP function with a finite thickness. The accuracy of the SP method is illustrated by comparison with known exact results. In the present paper, the high degree of versatility of the SP method is demonstrated by considering several types of active and passive particle suspensions.

Theoretical approaches have been limited to simple idealized systems (i.e., spheroidal particles in dilute suspensions), while experiments can be difficult to analyze, leaving computer simulations as one of the preferred methods. While this still requires advanced simulation methods and considerable computational resources, significant progress has been made in the past few decades, and there are now several well-established numerical techniques to explicitly account for the hydrodynamic interactions in a suspension of particles.^{5–11} However, not all of these methods can be (easily) applied to dispersions with non-Newtonian host solvents, or solvents with internal degrees of freedom, which severely limits the physical and biological systems they can study. To address this issue, we have proposed a direct numerical simulation (DNS) method called the smoothed profile (SP) method,^{12,13} which simultaneously solves the coupled equations of motion for both the host fluid and the particles, and allows for complex non-Newtonian fluids at arbitrary Reynolds numbers. The basic idea behind our method, and the key to its simplicity and efficiency, is to replace the sharp particle interfaces with smoothed (diffuse) profiles. This allows us to directly couple the fluid and particle motions, effectively treating the system as a single component fluid. This leads to the main advantage of our model with respect to alternative approaches, which is the ability of using a fixed Cartesian grid to solve the fluid equations of motion, as this allows for improved accuracy and efficiency (e.g., by adopting pseudo-spectral solvers based on the Fast-Fourier transform). We note that our method is closely related to Araki and Tanaka's fluid particle dynamics (FPD) method,^{14} which models the particles as a highly viscous fluid. Furthermore, while a similar approach to deal with the solid–fluid interaction has been proposed in earlier work,^{15} our proposed method presents a more general and robust framework.

The SP method is our approach to include the hydrodynamic interactions in a colloidal suspension. As presented here, it uses a DNS in which the Navier–Stokes equation is solved numerically. However, it is also possible to achieve the same goal without solving this equation directly by using a proxy which yields the correct hydrodynamics. The three most popular of these alternative computational fluid dynamics (CFD) methods, often used to simulate soft matter, are the lattice Boltzmann method (LBM),^{10,16–20} the dissipative particle dynamics (DPD) method,^{21} and the multi-particle collision dynamics (MPC)^{22–27} method. While the implementation details and inherent approximations vary between each of these methods (LBM, DPD, MPC, and DNS), they all rely on the same physical principle to couple the particle and fluid motion: local conservation of momentum. In fact, the SP method has also been successfully used together with the LBM.^{28}

Originally, we developed the SP method to simulate colloidal particles in liquid crystals^{29,30} and electrolyte solutions.^{13.,31–34} It has now been extended to study diffusion,^{35–37} coagulation,^{38,39} sedimentation,^{40–45} and rheological processes^{46–51} of colloidal particles in incompressible fluids. Recently, the method has also been adapted to allow for arbitrarily shaped rigid bodies,^{52} active particles,^{53} and self-propelled micro-swimmers.^{54–60} Extensions for complex host fluids such as compressible,^{61–63} viscoelastic,^{64,65} and multi-phase^{66–69} fluids have also been developed. Important contributions for further extending the SP method have also been published independently by other groups.^{70–72} The purpose of this paper is to provide a comprehensive summary of the SP method and its application to several problems of particulate systems in which the hydrodynamic interactions play crucial rules.

∇·u_{f} = 0, | (1) |

(∂_{t} + u_{f}·∇)u_{f} = ρ_{f}^{−1}∇·σ, | (2) |

σ = −pI + η[∇u_{f} + (∇u_{f})^{t}] | (3) |

Ṙ_{I} = V_{I}, _{I} = skew(Ω)·Q_{I}, | (4) |

M_{I}_{I} = F^{H}_{I} + F^{C}_{I} + F^{E}_{I} + G^{V}_{I}, | (5) |

_{I} = N^{H}_{I} + N^{C}_{I} + N^{E}_{I} + G^{Ω}_{I}, | (6) |

(7) |

We will now explain how to solve the system of coupled equations shown above with the SP method. In this approach, the velocity field is defined over the whole domain, in order to include both the host fluid and the solid particles, as

u(x,t) = (1 − ϕ)u_{f}(x,t) + ϕu_{p}(x,t), | (8) |

(9) |

There are two main difficulties when describing the whole domain, including the solid and fluid regions, as a continuum, which is what we propose for the SP method. First, the continuum field must satisfy the boundary condition at the particle surface. Second, the hydrodynamic forces need to be evaluated to solve for the particle dynamics. Here, it is crucial that the particle boundary conditions and the hydrodynamic forces be evaluated consistently with each other. To solve this problem, several methods have been developed. In the distributed Lagrange multiplier (DLM) method,^{75} a Lagrange multiplier field is introduced in the continuum equation, within the particle domain, in order to simultaneously solve for both the fluid velocity field and the particle velocity. In the fictitious boundary method based on a multi-grid finite element approach,^{76} hydrodynamic forces are explicitly calculated by introducing a fictitious boundary indicator. In the immersed-boundary method (IBM),^{77–79} markers distributed over the particle surface are used to apply frictional forces between the particle and the fluid velocity field near the surface. It is also possible to have a direct forcing of the particle velocity field in the particle domain, rather just at the particle surface, by adopting a volume-of-fluid representation for the particle.^{15} In the force coupling method (FCM),^{80} all forces on the particle, other than the hydrodynamic ones, are used to define a force field (assuming a Gaussian profile) that is then applied to the particle domain in the continuum equations. In the fluid particle dynamics (FPD) method,^{14,81} in addition to the particle force field, a large artificial viscosity is applied within the inner particle region to penalize fluid-like behaviour and thus enforce the rigid-body constraints on the velocity field. In the SP method, the calculation of the hydrodynamic forces and the boundary conditions are separated by using a fractional-step approach, without introducing an artificially large viscosity in the particle domain (in contrast to FPD).

∇·u = 0, | (10) |

(∂_{t} + u·∇)u = ρ^{−1}∇·σ + ϕf_{p}, | (11) |

(12) |

(13) |

(14) |

(15) |

(16) |

(17) |

(18) |

(19) |

Fig. 2 Schematic representation of the cross section of a spherical SP particle. Δ is the lattice spacing, a is the particle radius, and ξ is the interfacial thickness. The particle surface now has a finite volume ∼πa^{2}ξ supported by several grid points on a fixed 3D Cartesian grid. Reproduced from Phys. Rev. E, 71, 036707,^{12} Copyright 2005, with permission from the American Physical Society. |

For non-spherical particles, Luo et al.^{70} have generalized the definition of the SP function to arbitrary geometries by replacing |r_{I}(x, t)| − a with a signed distance function d_{I}(x, t) to the surface of the I-th particle. However, to define d_{I}, this requires that we define the position of the particle surface. Alternatively, complex shapes can be constructed using a set of spherical beads,^{52} which is the approach we have adopted.

For this, consider the I-th particle as being composed of a collection of n_{I} spherical beads. Since the constituent beads can overlap, a large class of particle geometries can be easily represented as an assembly of spherical beads. The SP function for such an assembly of beads is constructed as

(20) |

(21) |

Once the SP function of a complex shape has been defined, the hydrodynamic impulse on a particle and the corresponding momentum impulse on the medium are evaluated using eqn (15), (16) and (18).

Since the fluid inertia is solved for in the above SP algorithm, the added mass effect is accounted for in the hydrodynamic impulse in eqn (15). To specify the inertial effect, the integrand of the hydrodynamic impulse (15) is decomposed as

(22) |

(23) |

The accuracy of the flow field obtained from the SP method will depend on both the time increment h and the interfacial thickness ξ. Luo et al.^{70} reported that the accuracy does not improve monotonically as h decreases, but instead reaches an optimal value when ξ^{2} = O(ηh/ρ). This relation reflects the fact that a sufficiently large interfacial scale is required to resolve the viscous Stokes layer generated by the rigidity constraint impulse (18).

Due to the fractional step approach used in the SP method, it is relatively easy to apply different hydrodynamic solvers. In fact, the SP method has already been applied in combination with different flow solvers: the pseudo-spectral method,^{12,13,31,35–40,46–49,54,61} spectral element method,^{70,82} finite difference method,^{83} finite volume method,^{71} and LBM.^{28}

3.1.1 Thermal fluctuations.
Let us consider the Brownian motion of spherical particles dispersed in an incompressible fluid. Following eqn (4)–(6), the motion of the i-th particle is now governed by a set of Langevin equations:

where the external forces/torques are set to be zero, and it is assumed that the inter-particle interactions are given by a central-force potential, such that N^{C} = 0.

where the square brackets denote an equilibrium ensemble average. Here, α^{n} represents the noise intensity associated with each of the translational (n = V) and rotational (n = Ω) degrees of freedom of the particles. The noise intensity is controlled in such a way that the variance in the linear and angular velocities is constant, that is, 〈V_{i}^{2}〉 = C_{1} and 〈Ω_{i}^{2}〉 = C_{2}, where C_{1} and C_{2} are constant parameters. The time evolution of the noise intensity is described by α^{V}(t + h) = α^{V}(t)e^{1−Vi2(t)/C1} and α^{Ω}(t + h) = α^{Ω}(t)e^{1−Ωi2(t)/C2}, which play the role of a thermostat. The temperature of the system is defined by the diffusive motion of the dispersed particles. The translational particle temperature k_{B}T^{V} is determined by the long-time diffusion coefficient D^{V} of a spherical particle, as given by the Stokes–Einstein relation, k_{B}T^{V} = 6πηaD^{V}, with D^{V} obtained from the simulations. Similarly, the rotational particle temperature k_{B}T^{Ω} can be determined by the rotational diffusion coefficient D^{Ω}.

Ṙ_{i} = V_{i}, | (24) |

M_{i}_{i} = F^{H}_{i} + F^{C}_{i} + G^{V}_{i}, | (25) |

(26) |

The stochastic forces G^{V}_{i} and torques G^{Ω}_{i} satisfy the following properties:

〈G^{n}_{i}(t)〉 = 0, 〈G^{n}_{i}(t)G^{n}_{j}(0)〉 = α^{n}Iδ(t)δ_{ij}, | (27) |

3.1.2 Velocity auto-correlation of a single particle.
To test the validity of the present approach, we first consider the translational and rotational motions of a single Brownian particle in an incompressible Newtonian fluid. The translational velocity auto-correlation function β_{V}〈V_{i}(t)·V_{i}(0)〉 is shown in Fig. 3. Our numerical data (+) clearly approach the analytical solution in the long-time limit, exhibiting the appropriate power law decay Bt^{−3/2}. A similar trend is found for the rotational velocity auto-correlation function , as shown in Fig. 4.

Fig. 3 Pluses (+): translational velocity auto-correlation function for Brownian particles at β_{V} = 1.2 and φ = 0.002. The solid line is obtained from the analytical solution for the translational particle velocity of the corresponding drag problem. The dotted line shows the power-law Bt^{−3/2}, where B = 1/12ρ_{f}(πν_{f})^{3/2}, whereas the dashed line shows the exponential decay obtained from the Markovian Langevin equation. Reproduced from J. Phys. Soc. Jpn., 77, 074007,^{35} Copyright 2008, with permission from the Physical Society of Japan. |

Fig. 4 Crosses (×): rotational velocity auto-correlation function for Brownian particles at β_{Ω} = 1 and φ = 0.002. Pluses (+): the relaxation response of the angular velocity of a particle under a constant external torque N_{0}. The solid line is obtained from the analytical solution for the angular velocity of the corresponding drag problem. The dashed line shows the power-law Ct^{−5/2}, where C = π/32ρ_{f}(πν_{f})^{5/2}, and the dashed-dotted line shows the exponential decay obtained from the Markovian Langevin equation. Reproduced from J. Phys. Soc. Jpn., 77, 074007,^{35} Copyright 2008, with permission from the Physical Society of Japan. |

3.1.3 Brownian coagulation with hydrodynamics.
For over 100 years, ever since Smoluchowski presented his theory of Brownian coagulation in colloidal dispersions, there have been numerous theoretical, experimental, and simulation studies performed to assess the validity of this theory. Based on balance equations and a simple kinetic theory of coagulation, but ignoring the effect of inter-particle interactions, Smoluchowski's theory predicts a coagulation rate that is twice that of actual colloidal dispersions.^{84–88} To address this discrepancy, it has been suggested to introduce a correction factor into the coagulation rate, so as to account for the contributions from these missing inter-particle interactions.^{89–91} Several experimental evaluations of this correction factor have been performed, and they suggest that the reduction in the aggregation rate is mostly caused by the (repulsive) hydrodynamic lubrication forces between the particles.^{84–88,92} However, it should be noted that these examinations assume the validity of the correction factor itself. To understand the role played by the hydrodynamic interactions it is necessary to directly measure them.

with n:m = 200:100, where r is the inter-particle distance and ε is the interaction strength. The large value of the exponent in eqn (28) results in a very short-range interaction, appropriate for the rapid coagulation of sticky particles. For the results presented here, we chose a strong enough attractive potential, so that coagulated particles were prevented from breaking apart due to thermal fluctuations. For this, the value of ε was fixed to be 7.03k_{B}T. To improve the statistical significance of the data, five independent simulations were performed for each set of parameter values, and the results were averaged over.

with α^{0} = 1. This equation by Heine et al. reduces to the Smoluchowski result K/K^{0} = 1 in the dilute volume fraction limit, when (−logφ)^{−2.7} ≪ 1. However, since this formula was obtained from simulations using Langevin dynamics, it ignores the hydrodynamic interactions between the particles. To account for these missing contributions, we can treat α^{0} as a scaling or correction factor. Using a least-squares fit of α^{0} to our simulation results, we obtained α^{0} = 0.326. This result agrees well with the form of eqn (29), despite the differences in the way that the hydrodynamic interactions are considered. From this, we can conclude that while hydrodynamic interactions affect the magnitude of the coagulation rate, they do not affect the particle concentration dependency, at least in the dilute limit.

We have performed simulations of monodisperse spherical particles of diameter σ = 2a, with a = 4Δ, ξ = 2Δ, in a Newtonian host fluid. All simulations have been performed in a 3D cubic box, of length L = 256, under periodic boundary conditions, starting from a random initial configuration of particles. The number of particles was set to N_{p} = 187, 281, 375, 625, 938, 1877, 3754 and 6258. These conditions correspond to a particle volume fraction of φ = 0.003, 0.0045, 0.006, 0.01, 0.015, 0.03, 0.06 and 0.1, respectively, where φ = πσ^{3}N_{p}/6L^{3}. The particles interact via a modified Lennard–Jones (LJ) potential, defined as

V_{LJ}(r) = 4ε[(σ/r)^{m} − (σ/r)^{n}] | (28) |

The units of all variables are based on the three fundamental constants Δ, μ, and ρ_{f}. Accordingly, the unit of time is τ = ρ_{f}Δ^{2}/μ, and the unit of energy is ε_{0} = Δμ^{2}/ρ_{f}. Unless otherwise stated, we set Δ = 1, τ = 1, μ = 1, ρ_{f} = 1, ρ_{p} = 1, a = 4, σ = 8, ξ = 2, the Lennard–Jones time unit τ_{M} = ((M_{i}σ^{2})/ε)^{1/2} ≈ 131. Assuming dispersions of neutrally buoyant particles of radius 1 μm in water at room temperature, our unit length Δ and time τ correspond to 0.25 μm and 0.0625 μs, respectively. The discretized time step is h = 0.07481. Although this value is chosen from the stability condition of the Navier–Stokes equation, it can be used safely for the particle's equations because h is smaller than the oscillation period of the particle in the potential well .

The volume fraction dependence of Smoluchowski's coagulation rate K, normalized by the early-stage coagulation rate K^{0}, is shown in Fig. 5. This ratio is seen to asymptotically approach a fixed value in the dilute volume fraction limit, as φ → 0, whereas it gradually increases with increasing volume fraction. The coagulation rate we have obtained in the dilute limit is approximately 0.3–0.5 times smaller than that predicted by Smoluchowski's theory, in good agreement with previous experimental results (shown as open symbols in the figure), which reported values of K/K^{0} between 0.4 and 0.6.^{84–88} Likewise, the increase at higher volume fractions is in qualitative agreement with experiments and simulation results.^{85,93–95} The solid line in Fig. 5 shows the ratio of the coagulation rates K/K^{0} in the high-concentration limit, as given by the regression formula of Heine et al:^{93}

(29) |

Fig. 5 The volume fraction dependence of the normalized coagulation rate K/K^{0}. The filled symbols show our simulation results, and the open symbols the experimental results of Higashitani et al.^{86} The arrows indicate Smoluchowski's prediction, while the solid line shows the fit to the data using eqn (29). Reproduced from Phys. Rev. E, 86, 051403,^{39} Copyright 2012, with permission from the American Physical Society. |

3.2.1 Settling velocity at small Re.
We now consider the sedimentation of non-Brownian (T = 0) spherical particles interacting via a truncated LJ or Weeks–Chandler–Anderson potential

under a gravitational force F^{E}_{I} = (0, 0, −F_{g}). Simulations are performed using cubic simulation boxes, of length L = 64Δ, 128Δ or 256Δ, under periodic boundary conditions, with Δ being the grid spacing and our unit of length. The particle radius a is set to 4Δ and the particle volume fraction ranges from 0.01 to 0.5. To avoid any inertial effects, and remain within the Stokes regime, we perform all simulations at a particle Reynolds number of Re = ρ_{f}aV_{s}/η = 0.14, with being the corresponding Stokes velocity for a single particle.

(30) |

To account for the sedimentation of colloidal dispersions, we need to consider the fact that a sedimenting particle will experience a hindered settling due to the presence and the motion of the other particles. This includes the drag force induced by the back flow of the fluid, as well as any direct particle–particle interactions. The net effect of these interactions is to reduce the average settling velocity of the suspension, in comparison with the terminal velocity of a single isolated particle sedimenting under the same gravitational force. If a simulation is to reproduce this hindered settling, which has been amply studied, both theoretically^{96–101} and experimentally,^{102} it must accurately characterize the many-body hydrodynamic interactions at their origin. Thus, this hindered settling presents itself as a suitable parameter to assess whether or not the hydrodynamic interactions are properly described. We verify our SP method by comparing its predictions for the average settling velocity of a suspension, with theoretical,^{96–101} experimental^{102} and alternative simulation^{103,104} results, as shown in Fig. 6. The simplest semi-empirical relation to describe this behaviour is the power-law function proposed by Richardson and Zaki,^{96}V_{sed}/V_{s} = (1 − φ)^{n}, where V_{sed} = 〈V_{iz}〉 is the average settling velocity of the particle and n is the power law exponent. Different studies^{97,101–105} have obtained different values for this exponent, ranging from 4.7 to 6.55 in the Stokes flow regime. From our results, we find that a value of 5.3 provides the best fit to the data, which is well within the range of previous values. Our results are also in good agreement with the theoretical work of Hayakawa et al.,^{99} which also considers the effect of the hydrodynamic interactions. The agreement of our simulation results with previous experimental, theoretical, and simulation studies further validates our method.

Fig. 6 Average particle settling velocity V_{sed}, normalized by the Stokes velocity V_{s}, as a function of the particle volume fraction. This observed reduction with volume fraction highlights the hindered settling of these particles. Solid lines show the various theoretical predictions,^{96,98–101} and the symbols represent the experimental^{102} and simulation^{103} data. Reproduced from Soft Matter, 9, 10056–10068,^{40} Copyright 2013, with permission from the Royal Society of Chemistry. |

3.2.2 Settling velocity at Re > 1.
Let us now consider the inertial effects on the sedimentation of non-Brownian particles. For this, we vary the Reynolds number Re between 0.05 and 10, by increasing the strength of the gravitational force. In principle, our method is also applicable for Re > 10, but this requires that we reduce the time step and decrease the grid spacing, dramatically increasing the computational cost of our simulations. Our simulation results, presented in Fig. 7(a–c), show that the theoretical predictions obtained for the Stokes regime are no longer applicable due to the presence of inertial effects. To capture this deviation, one can use a modified expression for the sedimentation velocity, of the form:

where the terminal velocity, as a function of Re, is computed according to V_{t} = Reη/ρ_{f}σ. Here, both the scaling factor k and the exponent n are allowed to vary with Re. Koch et al. have reported k values of 0.92 ± 0.03 and 0.86 ± 0.04 for Re = 1 and 10, respectively, and suggested to use the following quadratic polynomial for n:

Previous work by Koch and Di Felice details the mechanism responsible for the initial decay in the hindered settling function at a low volume fraction, but their proposed mathematical expressions cannot be applied at φ ≃ 0. In order to include this missing information, we suggest a modified hindered settling expression:

where the values of k and n are, respectively, determined by

and

These simple fitting functions capture the behaviour of sedimenting particle dispersions as a function of Re and volume fraction. We note that these expressions are reduced to the Richardson and Zaki^{96} law at Re ≤ 0.1, whereas at 0.1 < Re ≤ 10 the second term (0.1 < Re ≤ 10) provides for the fast decay in the sedimentation velocity, as shown in Fig. 7.

V_{sed}/V_{t} = k(1 − φ)^{n}, | (31) |

n = 4.23 − 0.0526Re + 0.00111Re^{2}. | (32) |

V_{sed}/V_{t} = k(1 − φ)^{n} + (1 − k)exp(−φ/0.008), | (33) |

(34) |

(35) |

Fig. 7 Average sedimentation velocity V_{sed}, normalized by the terminal velocity V_{t} of an isolated sphere, as a function of the volume fraction for Re = 0.5, 1 and 10. Simulation results are fitted to the hindered settling function of eqn (33) and compared with other theoretical predictions^{101,106} and simulation results.^{104} Legends apply to all three figures. Reproduced from RSC Adv., 4, 53681–53693,^{43} Copyright 2014, with permission from the Royal Society of Chemistry. |

3.3.1 Implementation of sheared boundary conditions.
Understanding the macroscopic rheological properties of colloidal dispersions remains a challenging problem in colloidal science, particularly when considering non-Newtonian host fluids and non-spherical particles. Of particular interest are the non-Newtonian behaviours that appear whenever the suspension is subjected to shear, such as shear thinning or shear thickening, both of which are strongly related to changes in the microstructure of the suspension. In addition to experiments, DNSs are left as the only viable tool to measure the rheology. The SP method is an ideal tool for these studies, as it allows us to directly solve the Navier–Stokes equation for (Brownian and non-Brownian) arbitrarily shaped rigid particles under shear. We can account for the fluid and particle inertial effects and do not require the separation of time scales between fluid and particle relaxation times.

where a caret () is used here to denote tensorial components in the oblique frame. After the flow equations are solved, we transform back to the fixed laboratory frame to compute the particle–fluid interactions and solve for the particle dynamics.

and the particle contribution to the stress can be obtained by subtracting the average fluid contribution, s = Σ − 〈σ〉. In the case of zig-zag shear flow, a similar expression is obtained but with the external force field ρf_{s} appearing instead of the particle constraint force ρϕf_{p}. The apparent and intrinsic viscosities of the dispersions, η and [η], respectively, are defined as

with η_{f} = 〈σ^{12}〉/.

The easiest way to introduce a simple shear flow U = ye_{x} in the system is to add an external forcing term ρf_{s} into the modified Navier–Stokes equation (11), the role of which is to maintain the target shear flow (). When performing bulk 3D simulations under periodic boundary conditions, a linear zig-zag (possibly time-dependent) profile with two kinks (where the force is applied) can be used to satisfy the boundary constraints.^{46,47} The results obtained using this external forcing are in good agreement with alternative DNS methods and theoretical predictions, but the method presents some considerable drawbacks: first, the shear rate is not a precisely controlled quantity; second, there are non-physical kinks where the velocity gradient changes sign (at which point particles translate with no net rotation); and third, the zero-wavevector transport coefficients cannot be obtained (as the technique uses a finite-wavelength perturbation). To overcome these difficulties, one can adopt Lees–Edwards boundary conditions, solving the continuum equations on an (oblique) time-dependent coordinate system whose points are advected by the shear flow. Fig. 8 presents a schematic illustration of this coordinate transform. At time t = t_{0}, a solid spherical particle is shown dispersed in a solvent, which is discretized on a rectangular grid (Fig. 8(a)). After the shear-flow is applied (t > t_{0}), both the particle and the grid points are convected by the flow, with the grid undergoing a shear-induced deformation, in contrast to the particle whose shape remains unchanged. This situation can be represented in a fixed coordinate system (Fig. 8(b)), as well as a time-dependent oblique coordinate system which tracks the shear-flow (Fig. 8(c)), and in which the particle shape is deformed. This formulation is equivalent to the SLLOD algorithm used in non-equilibrium MD simulations to model (arbitrary) homogeneous and divergence-free flows, but the implementation of this algorithm within a grid-based system requires a careful treatment.^{49,108,109} The appropriate contravariant form of the Navier–Stokes equation for the peculiar velocity ξ = u − U within this time-dependent oblique coordinate system, with basis vectors (Ê_{1} = e_{x} + γ(t)e_{y}, Ê_{2} = e_{2}, Ê_{3} = e_{3}), is given as^{51}

(36) |

(37) |

Fig. 8 A schematic representation of the coordinate transformation used to model the effect of shear-flow in the DNS simulations. (a) A solid spherical particle is discretized on a rectangular lattice at time t = t_{0}. Both the particle and the solvent, represented by the grid points, are convected by the applied shear-flow (t > t_{0}). While the grid is deformed, the shape of the solid particle remains unchanged. This is depicted in (b) the fixed or laboratory reference frame and (c) the transformed or oblique reference frame. The transformation between (b) and (c) is defined in eqn (7) of our earlier paper.^{49} Reproduced from J. Chem. Phys., 134, 064110,^{49} Copyright 2011, with the permission of AIP Publishing. |

Within this formulation, the instantaneous volume-averaged stress of the flowing dispersion is defined in terms of the particle constraint force as

(38) |

(39) |

(40) |

3.3.2 Rheology of dilute suspensions.
The suspension viscosity for dilute systems of spherical and non-spherical particles is calculated and compared with analytical results at zero Reynolds numbers, as well as the results obtained using alternative (high-precision) methods. SP method simulation results for a single particle under simple shear can reproduce the known analytical solutions for the fluid velocity, stress distribution and particle angular velocity.^{110} Furthermore, the results at finite Reynolds numbers (Re ≲ 10) are also in agreement with the detailed finite element simulations of Mikulencak and Morris^{107} (see Fig. 9). Simulations for two-particle systems give similarly accurate stress predictions, in agreement with Stokesian dynamics calculations, provided that the surface-to-surface distance between the particles is larger than the grid spacing. When this condition is violated, we cannot expect the SP method to correctly account for the lubrication forces. However, if required, such lubrication corrections can be directly incorporated into the inter-particle interactions.^{111} The results for rigid linear bead chains are also in good agreement with theoretical predictions and SD simulations, allowing us to reproduce Jeffery's orbits, as well as the time-varying particle contribution to the stress.^{112} Simulations of Brownian flexible and rigid bead chains have revealed a tumbling motion, with a frequency F that depends only on Pe, following F ∝ Pe^{α}, and which exhibits a crossover between Brownian α = 2/3 and shear-flow α = 1 dominated behaviour.^{48} Furthermore, this transition in the tumbling frequency is shown to correspond to the transition between shear thinning and the second Newtonian regime,^{50} as first predicted by Hinch and Leal.^{113,114}

Fig. 9 (left) Angular velocity ω_{z} and intrinsic viscosity [η] of a single spherical particle as a function of Re. Filled symbols are from SP method simulations using Lees–Edwards boundary conditions, with a particle radius-to-system size ratio of ε/L = 0.0625; open symbols are from the finite element simulations of Mikulencak and Morris,^{107} for ε = 0.01 (squares) and ε = 0.1 (circles). (right) SP method results, with ε = 0.0625 and 10^{−3} ≲ Re ≲ 1, under zig-zag shear for the low-shear η_{0}- and high-shear η_{∞}-limiting viscosity of concentrated particle dispersions as a function of the particle volume fraction. Solid and dashed lines show the fits to the data using the Krieger–Dougherty relations, with Φ_{M} being the packing volume fraction at which the viscosity diverges. The left figure is reproduced from J. Fluid Mech., 792, 590–619,^{51} Copyright 2016, with permission from Cambridge University Press. The right figure is reproduced from Phys. Rev. E, 80, 061402,^{46} Copyright 2009, with permission from the American Physical Society. |

3.3.3 Rheology of non-dilute suspensions.
To understand the relationship between the non-Newtonian viscosity of spherical particle dispersions and the underlying microstructural changes, we report simulation results over a wide range of volume fractions φ ≲ 0.56 and Peclet numbers Pe = 6πη_{f}a^{3}/K_{B}T (10^{−2} < Pe < 10). For low to intermediate particle volume fractions φ ≲ 0.3, the viscosities are constant over the examined range of Peclet numbers, signalling Newtonian behaviour. However, for higher volume fractions (φ ≳ 0.4), shear-thinning behaviour is obtained, with a sharp decrease in the viscosity at high Peclet numbers. The low- and high-shear limiting viscosities (shown in Fig. 9), both monotonically increasing functions of volume fraction, are in good agreement with the semi-empirical Krieger–Dougherty relationship. Analysis of the structural relaxation time τ_{p}, defined as the time required for the Brownian particles to diffuse a distance equal to their radius, shows that shear thinning starts at approximately τ_{p} ∼ 1, reminiscent of the non-Newtonian behaviour of supercooled liquids.^{115} Furthermore, to study the viscoelastic properties, we also consider an oscillatory shear flow and measure the storage and loss moduli, which are shown to be in excellent agreement with experimental data.^{47}

In contrast to incompressible fluids, for which the speed of sound is infinite and the hydrodynamic interactions propagate instantaneously, followed by a temporal evolution due to viscous diffusion, in compressible fluids the sound propagates at finite speeds. This difference strongly affects the time evolution of the hydrodynamic interactions. To quantify the relative importance of the sound propagation, compared to the viscous diffusion, we define the compressibility factor ε as

(41) |

The fluid dynamics is now governed by the following hydrodynamic equations:

(42) |

(43) |

(44) |

(45) |

We perform 3D simulations for two spherical particles (a = 4), labelled 1 and 2, in a cubic simulation box (L = 256), under periodic boundary conditions, for various compressibility factors. At t = 0 an impulsive force is exerted on particle 2, and their correlated motions are analyzed in order to characterize the hydrodynamic interactions between them. The cross-relaxation tensor γ_{12} is defined as the normalized change in the velocity of particle 1:

(46) |

γ_{12}(R, t) = γ^{‖}_{12}(R,t) + γ^{⊥}_{12}(R,t)(I − ). | (47) |

(48) |

Fig. 10 Velocity field around particles in a compressible fluid after the application of an impulsive force on the right particle in the parallel (a and c) or perpendicular (b and d) direction. Simulation results at t/τ_{ν} = 6.25 × 10^{−1} for compressibility factors with values (a and b) ε = 0.1 and (c and d) ε = 0.6. The divergence of the velocity field ∇·v is described by a colour scale extending from negative (darker) to positive (lighter) divergence. The value of divergence is normalized by τ_{ν}/Re_{p}. Reproduced from Phys. Fluids, 25, 046101,^{62} Copyright 2013, with the permission of AIP Publishing. |

3.5.1 A model for ternary systems.
We now consider the adsorption of solid particles at the fluidic interface of moving droplets, which is relevant to many industrial processes, including the stabilization of emulsions and foams,^{116} the armouring of droplets in capillary tubes,^{117} and the recovery of mineral particles by rising gas bubbles.^{118,119} The encapsulation of air bubbles in seawater presents another striking example in which the addition of solid particles dramatically changes the dynamics of these multi-phase systems.^{120} Performing DNS of such a ternary, fluid–fluid–particle system is challenging due to the complexity of capturing all the relevant physical mechanisms, including capillary effects, three-phase flow hydrodynamics, and inter-particle interactions. Most studies have focused on the dynamics of single particles trapped at planar fluidic interfaces^{121–123} or at the fluidic interface of an immobile spherical droplet.^{124,125}

where denotes the spatial extent of our ternary system, k_{B} is the Boltzmann constant, T_{0} is the reference temperature, v_{0} is the reference unit volume, and f is the non-dimensional free energy density. The variables (ϕ, ϕ_{A}, and ϕ_{B}) in the free energy density f are usually referred to as “order parameters”.^{126} Because the sum total of the volume fractions is constrained to be unity ϕ_{A} + ϕ_{B} + ϕ = 1, the free energy density can be conveniently recast as a function of only two order parameters, which we take to be ψ = ϕ_{A} − ϕ_{B} and ϕ. For simplicity, We refer to ψ as “the” order parameter and ϕ as the particle or fibre (i.e., bead-chain) volume fraction. Consistent with previous studies,^{127,128} we adopt the following form for the free energy densities of the mixture:

where f_{0}(ψ) is a double well function with minima at ψ = −1 and ψ = +1, f_{1}(ψ, ϕ) is an energy penalty introduced to impose ψ = ψ_{0} inside the particles, ψ* is a control parameter that sets the affinity of the beads to the fluid phases, and ξ is the interfacial thickness defined in eqn (19).

As suggested by Araki and Tanaka,^{129} imposing ψ = ψ_{0} inside the solid domain determines the bead affinity. A positive value ψ* > 0 results in a stronger affinity to fluid A, whereas ψ* < 0 results in a stronger affinity to fluid B, and ψ* = 0 results in an equal affinity to the two fluid phases. In fact, for the special case when ψ* = 0, it can be shown that the bulk driving force term δf_{b}/δψ coincides with that suggested by Kim,^{130} resulting in a factor of 3/2 for the f_{1} term in eqn (52). We stress the fact that ψ* is an input parameter, while ψ_{0} is an output parameter. In all the simulations we have performed, the value of ψ_{0} coincides exactly with the theoretical prediction of eqn (53). Because it is straightforward to verify that ψ = ψ_{0} from a contour of the order parameter, we only specify the value of ψ_{0} in what follows.

where M is the mobility of the binary fluid, u is the velocity field, and μ_{ψ} = δ/δψ is the chemical potential associated with the order parameter ψ.

where ϕf_{p} is again the penalty term that enforces the rigid-body constraints on the beads and μ_{ϕ} = δ/δϕ is the chemical potential with respect to the fibre volume fraction ϕ. The last source term ϕ∇μ_{ϕ} is motivated by the ternary fluid studies of Boyer et al.,^{131} but our tests show that it is only necessary to improve the stability of the solution in 2D. For 3D simulations, remarkable agreement with theory can be achieved without it.^{68}

Consider a ternary system in which solid particles are dispersed in a binary fluid mixture. A schematic representation of this setup is given in Fig. 11. The field ϕ(x, t) will denote the volume fraction of the solid constituent at position x and time t. The binary fluid mixture will separate into its two immiscible fluid constituents, “A” and “B”. Constituent A represents the host fluid, whereas constituent B represents the fluid inside the droplet phase. The volume fractions of the two fluid constituents are denoted as ϕ_{A}(x, t) and ϕ_{B}(x, t). The separation of the binary fluid mixture into A and B phases is driven by the minimisation of the following free energy functional:

(49) |

(50) |

(51) |

(52) |

Fig. 11 Schematic representation of a reference ternary system, showing a fluid (B) droplet rising in the host fluid (A). The solid particles (S) adsorb at the fluidic interface of the binary fluid. Reproduced from Phys. Rev. Fluids, 3, 094002,^{68} Copyright 2018, with permission from the American Physical Society. |

The bulk term f_{b} = f_{0} + f_{1} is a function with a single minimum located at ψ = ψ_{0} inside the particle domain (see Fig. 12). This bulk term ensures that ψ = +1 in fluid A, ψ = −1 in fluid B, and ψ = ψ_{0} in each solid bead. The last term in eqn (50), ξ^{2}|∇ψ|^{2}/2, gives an excess energy contribution across each diffuse interface. We note that the control parameter ψ* and the position of the single well minimum ψ_{0} are dependent on each other. At equilibrium, δf_{b}/δψ = 0, which results in the following dependency:

(53) |

Fig. 12 Bulk free energy density f_{b} = f_{0} + f_{1}. Reprinted from J. Comp. Phys., 409, G. Lecrivain, T. B. P. Grein, R. Yamamoto, U. Hampel, and T. Taniguchi, Eulerian/Lagrangian formulation for the elasto-capillary deformation of a flexible fibre,^{69} 109324, Copyright 2020, with permission from Elsevier. |

The order parameter ψ(x, t) is updated in time following a modified Cahn–Hilliard equation:

(54) |

We use the SP method to solve for the coupled hydrodynamics and particle dynamics. For the type of system we are interested in, namely the deformation of flexible fibres, capillary effects provide the dominant contributions. Since a mismatch in the fluid densities or viscosities would only alter the transient dynamical behavior, but not the equilibrium steady-states, we will only consider the case of low-Reynolds number flows with equal densities and viscosities for all constituents, such that ρ_{A} = ρ_{B} = ρ and η_{A} = η_{B} = η. As such, the two volume fractions of the pure fluid phases, ϕ_{A} and ϕ_{B}, are not directly needed for our implementation. The total fluid velocity, satisfying the incompressibility condition ∇·u = 0, is advanced in time by solving the following modified Navier–Stokes equation:

(55) |

3.5.2 Elasto-capillary deformation of a fibre.
For simplicity, we have solved the equations in their non-dimensional form, using the Reynolds number Re, the Weber number We, the Peclet number Pe, and the elasto-capillary number Ec, defined as

where ρ_{0} = ρ, η_{0} = η, L_{0} = Δ, and U_{0} are the reference density, viscosity, length and velocity scales, respectively. Here U_{0} can be the terminal velocity of a single particle in the reference/host fluid. The diffusion coefficient is defined as D_{0} = e_{0}M, where e_{0} = k_{B}T_{0}/v_{0} is the reference free energy density of eqn (49). The reference surface tension is given by σ_{0} = e_{0}L_{0}. The term k_{0} is a measure of the fibre stiffness, with the axial beam force between adjacent beads at a separation r given by |F^{b}(r)| = k|(r − 2a)|. The two limiting cases of Ec = 0 and Ec → ∞ correspond to independent (non-attached) beads and a fully rigid fibre, respectively.^{69} We find that the dynamics of the ternary system is particularly sensitive to the ratio of the Weber number to the Reynolds number We/Re. For We/Re ≳ 1, we observe large amplitude oscillations in the motion of the fibre at the A/B fluidic interface, which can require long simulation times to reach equilibrium. For this reason, in what follows we set this ratio to be We/Re ≃ 10.

(56) |

(57) |

(58) |

(59) |

We report on simulations to study the capillary-induced tumbling of a fibre in a stratified fluid, which allow us to assess the validity of our ternary fluid model by checking the equilibrium orientation. We consider a stagnant fluid with no externally imposed flow, such that the fibre rotation is driven solely by the capillary effects at the three-phase contact line. Fig. 13 shows simulation snapshots illustrating the capillary-induced tumbling of a flexible fibre with equal affinity to the two fluid constituents in a stratified 3D system. The A/B fluidic interface is represented by the ψ = 0 isosurface, drawn here using a wire-frame representation. These results are obtained for a fiber with the number of beads N_{b} = 8, bead radius r_{b} = 6, particle interfacial thickness ξ = 2Δ, in a rectangular system of size 120 × 120 × 32, under periodic boundary conditions. The Reynolds, Weber, Peclet, and elasto-capillary numbers are Re = 0.02, We = 0.1, Pe = 1, and Ec = 10, respectively. Further simulations can be found in ref. 69.

Fig. 13 Capillary-induced tumbling of a flexible fibre in a 3D stratified system. Reprinted from J. Comput. Phys., 409, G. Lecrivain, T. B. P. Grein, R. Yamamoto, U. Hampel, and T. Taniguchi, Eulerian/Lagrangian formulation for the elasto-capillary deformation of a flexible fibre,^{69} 109324, Copyright 2020, with permission from Elsevier. |

Next, we have considered the bending of a fibre at the interface of a fluid droplet. At time t = 0, the droplet is initialised to a sphere filled with the fluid constituent A, that is, ψ = 1. Fig. 14 shows the simulation snapshots of the equilibrium solutions, for a fibre with an equal affinity to both fluid constituents, ψ_{0} = 0. The droplet is represented by the ψ = 0 isosurface. Details regarding the unsteady dynamics leading up to the encapsulation process can be found in ref. 69. These results are obtained with N_{b} = 10, r_{b} = 6, ξ = 2Δ, Re = 0.02, We = 0.1, Pe = 1, for a system size of 128 × 120 × 120.

Fig. 14 3D equilibrium solutions of droplet encapsulation for high values of E_{c}. Reprinted from J. Comput. Phys., 409, G. Lecrivain, T. B. P. Grein, R. Yamamoto, U. Hampel, and T. Taniguchi, Eulerian/Lagrangian formulation for the elasto-capillary deformation of a flexible fibre,^{69} 109324, Copyright 2020, with permission from Elsevier. |

3.6.1 A model for charged colloidal dispersion.
The dynamics of charged colloidal dispersions is crucial to many physical, chemical, and biological systems, as well as of great importance for the engineering field.^{1} Among them, the electrophoresis of charged colloidal particles is one of the most salient examples. When a charged particle is placed in an electrolyte solution, it will develop an electric double layer (EDL), as a cloud of counter-ions is attracted to the surface of the particle. Upon the application of an external electric field, this charged particle will move with an electrophoretic mobility determined by the balance between the electrostatic driving force and the hydrodynamic frictional force acting on the particle. This mobility will depend strongly on the EDL whose distribution becomes anisotropic due to the external field and the friction between the ions and the fluid. The time evolution of the colloidal particles, the ions, and the host fluid are described by the coupled equations for the rigid-body dynamics, hydrodynamics (Navier–Stokes), and electrostatics (Poisson) with the appropriate boundary conditions on the surface of the colloidal particles. We briefly outline our numerical model for such charged colloidal dispersions and how it differs from the standard model discussed in the introduction, and then demonstrate the reliability of the SP method by comparing our numerical results with well-established approximate theories.^{132–135} Finally, the electrophoretic mobilities of dense dispersions, where simulation results deviate from mean-field type theories,^{136,137} as well as the field-induced pearl formation, are presented.

with the incompressibility condition ∇·v = 0, where Ψ_{ex} = −E·r is the electrostatic potential due to the imposed external electric field E. The electrostatic potential Ψ(r) is determined by solving the Poisson equation ε∇^{2}Ψ = −ρ_{e} with ε being the dielectric constant of the host fluid. For simplicity we consider only dielectrically matched particles, such that ε_{f} = ε_{p} = ε.

We include all the electrostatic interactions, which consist of the driving from the external electric field E, as well as the forces due to the distribution of charges in the system, in F^{H}_{i}.

where f_{α} is the ionic friction coefficient, I is the unit tensor, and n is a unit-vector field defined by n = −∇ϕ/|∇ϕ|. In our method, the physical ionic density fields are actually C_{α}(r,t) = (1 − ϕ(r, t))C_{α}*(r,t), which are defined using the auxiliary density field C_{α}*(r,t). This definition explicitly avoids the penetration of ions into the colloidal domains without using any artificial potentials, which would require smaller time increments. The operator (I − nn) in eqn (62) ensures the conservation of C_{α} by enforcing the no-penetration condition, as n·∇μ_{α} = 0 is directly assigned at the diffuse particle interface. In this way the total charge neutrality of the system is guaranteed. The chemical potential is defined as^{138,139}

If we set v = 0 in eqn (62), the equilibrium (t → ∞) ionic density corresponds to the Boltzmann equation:

The bulk salt concentration _{α} is related to the Debye screening length , where λ_{B} = e^{2}/4πεk_{B}T is the Bjerrum length, which is approximately 0.72 nm in water at 25 °C.

Let us consider N charged spherical particles of total charge Ze (e being the unit charge) distributed uniformly on its surface. We define the spatial distribution of the surface charge to be eq(r) = Ze|∇ϕ|/4πa^{2}, with ϕ being the SP function. The particles are dispersed in a host fluid containing ions of species α, with charges Z_{α}e. The local number density of the α-th ionic species is C_{α}(r, t) at time t, and the total charge density is smoothly defined over the entire domain by . The complete dynamics of the system is obtained by solving the following time evolution equations.^{12,13}

(i) The Navier–Stokes equation:

ρ(∂_{t} + v·∇)v = −∇p + η∇^{2}v − ρ_{e}∇(Ψ + Ψ_{ex}) + ϕf_{p}, | (60) |

(ii) Newton's and Euler's equations of motions:

(61) |

(iii) Advection–diffusion equation:

(62) |

μ_{α} = k_{B}TlnC_{α}* + Z_{α}e(Ψ + Ψ_{ex}), | (63) |

C_{α}* = _{α}exp[−Z_{α}e(Ψ + Ψ_{ex})/k_{B}T]. | (64) |

Simulations are performed in a 3D cubic box, of length L = 64Δ (Δ = 4πλ_{B} is the grid spacing and our unit of length), under periodic boundary conditions. We use spherical particles of radius a = 5 and interfacial thickness ξ = 2. The host fluid contains a 1:1 electrolyte composed of monovalent counter-ions (α = +) and co-ions (α = −). The units of energy and electrostatic potential are k_{B}T and k_{B}T/e, respectively. The non-dimensional mobility parameter m_{α} = 2εk_{B}Tf_{α}/3ηe^{2} is set to be m_{+} = m_{−} = 0.184, which corresponds to that of KCl solution at 25 °C. Finally, our unit of time τ = Δ^{2}f_{+}/k_{B}T would correspond to 0.44 μs.

3.6.2 Electrophoresis of charged colloids.
We first consider a single charged particle moving with velocity V = (−V, 0, 0) under a constant external electric field E = (E, 0, 0). The electrophoretic mobility V/E as a function of the zeta potential ζ, defined as the electrostatic potential on a slipping plane, is given by

in the limit when ζ is small.^{1} The Smoluchowski equation, corresponding to f = 1, is valid in the limiting case when κa → ∞,^{132} while the Hückel equation, corresponding to f = 2/3, is valid in the opposite limit of κa → 0.^{133} An expression for f = f_{H}(κa) that is valid for an arbitrary κa has been derived by Henry.^{134} These equations all predict a mobility that is proportional to ζ; however, this is not true for large values of ζ, when the relaxation effects due to deformation of the EDL can no longer be neglected. O’Brien and White have proposed an approximate theory that is also valid for larger values of ζ.^{135}

V/E = fεζ/η | (65) |

We perform simulations for the electrophoresis of a single particle in the linear response regime E ≲ 0.15 and compare their results with the O’Brien–White theory. A constant uniform electric field of magnitude E = 0.1 is applied. The terminal velocity V is calculated as a function of the charge on the particle, for 50 ≤ −Z ≤ 750, with κ^{−1} = 10, which would correspond to a salt concentration of 11 μmol l^{−1} at 25 °C in water. In particular, we set the surface charges to be Z = −50, −100, −200, −300, −400, −500, and −750, corresponding to a dimensionless zeta potential of y ≡ eζ/k_{B}T = 0.525, 1.044, 2.035, 2.927, 3.692, 4.332, and 5.510, respectively. We choose ν = η/ρ = 5 to maintain a small Reynolds number Re = aV/ν. In our simulations, as well as the O’Brien–White theory, the zeta potential is assumed to be the electrostatic potential at the particle surface, ζ = Ψ|_{surface}. Our simulation results are in excellent agreement with the O’Brien–White theory, even in the non-linear regime y ≥ 2, with an error of only a few percent.^{31} This agreement is unprecedented when compared to alternative simulation methods at this level of description.

Our simulation method applies just as easily to dense dispersions, allowing us to examine the effect of particle concentration on the electrophoretic mobility. The linearized theory for the mobility of a single particle, as described by eqn (65), is still valid for dense dispersions, as long as E is small; however, f now depends on both κa and φ. Simulations for Z = −100 and E = 0.1, at various particle volume fractions φ ≡ 4πa^{3}N/3L^{3}, are performed to measure the value of f(κa, φ) = ηV/εζE. We consider Debye lengths κ^{−1} of 5 and 10, which correspond to κa = 1 and 0.5, and salt concentrations of 44 μmol l^{−1} and 11 μmol l^{−1}, respectively. Fig. 15 shows typical simulation snapshots for systems with κa = 1 and (a) FCC, (b) BCC, and (c) random particle configurations.^{140} The colour map indicates the charge density, shown over a cross-sectional slab perpendicular to the z-axis. In the case of FCC and BCC, E is applied perpendicular to the (1, 0, 0) and (1, 1, 1) faces, but we obtain very small differences in the electrophoretic mobility (within only 1%).

Fig. 15 Simulation snapshots for the electrophoresis of dense dispersions with (a) FCC, (b) BCC, and (c) random particle configurations. The colour map represents the total ionic charge density in a plane perpendicular to the z-axis. The electric field is applied in the +x direction, normal to the (1, 0, 0) face for the FCC and BCC lattices. Reproduced from Phys. Rev. Lett., 96, 208302,^{31} Copyright 2006, with permission from the American Physical Society. |

The mobility coefficient f(κa, φ) for κa = 1 and 0.5 is plotted as a function of φ in Fig. 16(a and b), respectively. It is seen that f decreases rapidly with increasing φ. Furthermore, we see no significant difference in f with respect to the particle configurations (FCC, BCC, or random). A theoretical model for the electrophoretic mobility based on the cell model has been proposed by Levine and Neale.^{136} For this, they have studied a spherical particle (radius a) at the center of a spherical container or cell (radius b), and calculated the velocity V as a function of κa and φ = (a/b)^{3}. A similar expression for the mobility coefficient f, also based on this cell model, has been proposed by Ohshima.^{137} Ohshima's prediction is shown in Fig. 16, together with our numerical results. The agreement between our simulations and Ohshima's theory is generally good for small Debye lengths (κ^{−1} = 5 = a), but shows considerable deviations at large Debye lengths (κ^{−1} = 10). In both cases, the simulation results predict larger values of f than Ohshima's theory as φ increases. We consider that this deviation arises because of the overlapping of the EDLs for larger φ, as this effect is completely neglected in the theory. To test this, we estimate the effective particle radius a + κ^{−1} of the ionically dressed particles and define an effective volume fraction φ_{eff} ≡ 4π(a + κ^{−1})^{3}N/3L^{3} = (1 + (κa)^{−1})^{3}φ. Fig. 16 clearly shows that our results agree well with Ohshima's theory for φ_{eff} ≤ 1, where the overlap effect is small. However, for φ_{eff} > 1, where the degree of EDL overlap is large, deviations between our simulations and the theory are notable. These are the first simulations to provide quantitative data necessary to examine the reliability of Ohshima's theory in dense colloidal dispersions. Our results are consistent with more recent studies that have taken EDL overlap into account.^{141–143}

Fig. 16 The mobility coefficient f(κa, φ) as a function of the volume fraction φ for (a) κa = 1 and (b) κa = 0.5. The solid lines show the results given by the approximate theory proposed by Ohshima.^{137} This theory is confirmed to be accurate for φ_{eff} ≤ 1 but tends to deviate from our numerical results for φ_{eff} > 1, where the overlap of the EDLs becomes notable. Reproduced from Phys. Rev. Lett., 96, 208302,^{31} Copyright 2006, with permission from the American Physical Society. |

Fig. 17 Simulation snapshots depicting the time evolution of 12 free particles, initially configured as two stacked linear chains aligned parallel to the oscillating electric field, for ζ = 0.5, κa = 2.5, E = 0.3, and ω* = 0.001. Assuming a colloidal radius of 1 μm and an ionic diffusion constant of D_{α} = 2 × 10^{−9} m^{2} s^{−1}, the simulation times t correspond to (a) 2.1 × 10^{−4} s, (b) 1.4 × 10^{−2} s and (c) 8.4 × 10^{−2} s. The colour map represents the charge density. Reproduced from Soft Matter, 14, 4520–4529,^{34} Copyright 2018, with permission from the Royal Society of Chemistry. |

3.6.3 Field-induced order formation.
The field-induced anisotropic interactions between like-charged colloidal particles is studied using the SP method. In particular, we are able to include the polarization of the EDL under external AC electric fields (at frequency ω*) explicitly. Consider a field that is parallel to the distance vector between two particles, such that the induced dipoles (arising from the polarization of the EDL) are aligned in an end-to-end fashion, with the partial positive charge on one particle attracted to the partial negative charge of its neighbour. This situation is depicted in Fig. 17, where particles in the same row feel an attractive inter-particle force (whereas particles on opposite rows are repelled). If the induced dipoles are strong enough to overcome their Coulombic repulsion, i.e., if the field strength is large enough, the particles will experience an effective attraction. Since this attractive interaction depends on the polarization of the EDL, it will also depend strongly on ζ and κa. For a fixed value of the field strength E_{0}, the interaction is attractive for small ζ but repulsive for large ζ. This is due to the fact that at large ζ (high surface potential) more counter-ions are attracted to the charged particle surface, leading to an increase in the charge density around the particle. In turn, this makes the particle more insensitive to external perturbations and more difficult to polarize. In the end, this results in smaller induced dipole–dipole interactions between the particles. Similarly, for thin double layers (large κa) the EDL is more difficult to polarize, compared to thick double layers (small κa), because the counter-ions are more strongly bound to the particle surface. Furthermore, we find that the dynamic response of the colloids and ions is frequency dependent. In particular, we observe two dynamically distinct regimes, corresponding to whether ω* is faster or slower than the ionic diffusion rate ω_{κ}* over the size of the EDL. In the low-frequency regime ω* < ω_{κ}*, the polarization of the EDL has enough time to develop in response to the electric field E(t) and can therefore induce a strong inter-particle force. In contrast, in the high-frequency regime ω* > ω_{κ}* the ions cannot follow the field fast enough to fully polarize the EDL. This results in an inter-particle force that decreases with increasing ω*, up until the point where the system reproduces the results in the absence of an external electric field E_{0} = 0.

To characterize a viscoelastic suspending medium using the SP method, the viscoelastic stress is added to the continuum eqn (11):

ρ(∂_{t} + u·∇)u = ∇·σ + ∇·σ_{p} + ρϕf_{p}, | (66) |

The shear rheology of a suspension in an Oldroyd-B fluid is investigated by using the SP method.^{64,65} The Oldroyd-B model exhibits a shear-rate-independent viscosity and finite elasticity, which serves as a model for Boger fluids. The polymer stress in the single-mode Oldroyd-B fluid is governed by

(67) |

(68) |

Since shear thickening in Boger fluid media occurs at dilute concentrations, this thickening is attributed mainly to the interaction between a particle and the viscoelastic stress around the particle. The flow around a particle has been previously examined by different DNS studies.^{64,154,163} Since the flow is forced to circumvent the particle, the direction of the shear flow necessarily changes around the particle. Depending on the conformation orientation, the elastic stress caused by the stretched conformation contributes to an increase in either the macroscopic shear stress or the macroscopic first normal stress difference. Variations in several quantities along a recirculating streamline around a particle at Wi = 2 and β = 0.5 are shown in Fig. 18. The location in Fig. 18 is specified by the angle ϕ_{s} from the velocity gradient direction. When σ_{p,xy}/η_{p} > 1, the elastic stress contributes to bulk shear thickening (Fig. 18(a)). For elastic thickening to occur, both the strong stretch of the conformation and the appropriate orientation of the conformation are necessary. For the conformation to be stretched, a large shear rate is not necessarily required (Fig. 18(b)), but the alignment of the major axes of the conformation and the strain-rate tensors, i.e., θ_{C} ≈ θ_{D} in Fig. 18(c), is important. The efficiency of the conformation stretch can be evaluated by the effective molecular elongation rate,^{174}_{p} = n_{1}n_{1}:D shown in Fig. 18(d), where n_{1} and D are the primary directions of C and the strain-rate tensor, respectively.^{64} Furthermore, the flow modulated by a particle is almost irrotational in front of and behind the particles,^{64,154,163} which promotes conformation stretching because the vorticity induced misalignment of C from D is suppressed. These effects of a suspended particle do not depend on the particle shape; therefore, elastic thickening is expected to occur even in suspensions of non-spherical particles.

Simulations of multi-particle systems in an Oldroyd-B fluid succeeds in quantitatively reproducing the experimentally measured shear-thickening behaviour with a Boger fluid matrix up to φ ≤ 1 for the first time.^{65}

3.8.1 A model for swimming systems.
Micro-swimmers, defined as particles that self-propel in a viscous fluid, among which we find spermatozoa, bacteria, and algae, are ubiquitous in biology and life sciences and are one of the main representatives of active matter systems. From a purely theoretical standpoint, the study of active matter presents itself as fertile ground on which to test and develop theories of non-equilibrium systems. From an applied perspective, particularly with the introduction of synthetic micro-swimmers, there is also the prospect of developing novel medical applications, such as cargo transport systems for targeted drug delivery. However, understanding the role played by hydrodynamic interactions in determining the physical properties of the system remains an incredibly challenging task.^{175}

where P_{n}′ is the derivative of the n-order Legendre polynomial and B_{n} and C_{n} are the coefficients for the n-th order polar and azimuthal modes, respectively. Most studies neglect the azimuthal modes (C_{n} = 0) and all polar modes beyond the second order B_{n} = 0 (n > 2), as this results in the simplest non-trivial model that can represent both pushers (generating propulsion at the back) and pullers (generating propulsion at the front). In this case, the surface velocity becomes (α = B_{2}/B_{1})

where B_{1} determines the steady-state swimming speed for a single particle in an unbounded fluid (U = 2/3B_{1}) and B_{2} is the strength of the stresslet. The ratio of these modes α (also denoted as β in the literature) determines the type of swimmer, with α > 0 and α < 0 corresponding to pullers and pushers, respectively, and α = 0 to the so-called neutral or potential swimmers. A schematic diagram of the squirmer representation for pushers/pullers, together with the different kinds of self-generated flows, is shown in Fig. 19. The nature of swimming strongly affects the hydrodynamics of the generated flow, with the flow generated by neutral swimmers undergoing 3D decay as r^{−3}, while pullers (pushers) generate contractile (extensile) flows with a weaker decay as r^{−2}.^{181} This has important consequences for collective behaviour, as shown below.

To solve this equation, we incorporate an additional step in the fractional-step approach. After updating for the advection and viscous terms and obtaining u*, but before enforcing the rigidity constraints through ϕf_{p}, we solve for the surface constraint force ϕf_{sq} to obtain an updated velocity field u** that satisfies the squirmer boundary conditions. For this, we consider an additional SP function ϕ^{sq}, which is non-zero only within the interface domain of the particle

The new velocity field can then be obtained as

where the first term is the one that actually fixes the slip boundary conditions, the second term is added to ensure local momentum conservation (i.e., the squirmer experiences a counteracting force if it pushes the fluid at the boundary), and the third term guarantees that the final fluid velocity field is incompressible. Since the updated particle velocities are not yet known at this point of the calculation, we adopt an iterative solution and enforce the slip velocity with respect to V_{I}^{†} (W_{I}^{†}) using the velocities at the previous time step at an initial guess. The hydrodynamic force and torque are computed as before, using u** instead of u*, and used to update the particle velocities. This procedure is iterated to ensure consistency between the velocities V^{†}_{I} used to enforce the slip boundary conditions and the particle velocities at the next time step V^{n=1}_{I}. The method has been validated for the motion of a single swimmer, where both the particle velocity and the fluid velocity show good agreement with theoretical predictions.^{54}

We focus here on a generic model for micro-swimmers called squirmers. Originally developed by Lighthill and Blake as a model for the ciliary propulsion of paramecia,^{176,177} this model has become standard for studying the hydrodynamic interactions of rigid micro-swimmers.^{178} Instead of modelling the synchronized beating of cilia at the surface, which would result in a computationally prohibitive description, one can reproduce the same effects (at the scale of the particle size) by imposing a prescribed boundary condition at the surface of the rigid particle. Here, we consider only spherical swimmers (radius a), but extensions of the squirmer model to spheroidal particles also exist.^{179,180} Let {ẽ_{i}} denote the ortho-normal basis vectors defining the body-axis reference frame, with ẽ_{3} = ẽ being the fixed swimming axis of the particle. A given point on the particle surface is parametrized in terms of the polar (colatitude) θ and azimuthal λ angles, with θ and λ being the corresponding unit tangent vectors, ϱ the unit vector in the radial direction and θ = arccosϱ·ẽ. The imposed “slip” velocity for the squirmer (assuming zero radial velocity) is^{181}

(69) |

(70) |

Fig. 19 Schematic representation of the propulsion mechanism and flow profiles of a pusher and a puller (a and b), respectively. These swimmers can be represented using Blake's squirming model, in which the detailed propulsion mechanism is replaced by a specified slip velocity at the surface of the particles for pushers and pullers (c and d), respectively. Reproduced from Soft Matter, 9, 4923–4936,^{54} Copyright 2013, with permission from the Royal Society of Chemistry. |

To solve for the hydrodynamics of these micro-swimmer suspensions using the SP method, we introduce an additional constraint force ϕf_{s}q to the continuum eqn (11) to enforce the slip velocity at the surface^{54,58}

ρ(∂_{t} + u·∇)u = ∇·σ + ρϕf_{p} + ρϕf_{sq} | (71) |

(72) |

(73) |

(74) |

3.8.2 Self-diffusion of a micro-swimmer.
Experiments of tracers in swimming suspensions show that hydrodynamic interactions with swimmers, which continuously stir the fluid, are responsible for enhancing the tracer diffusion (beyond that due to thermal fluctuations).^{182,183} Theoretical predictions based on a kinetic theory approach, as well as simulations and experiments, show that this induced diffusion scales with the flux of swimmers (∝UΦ_{swimmer}).^{184,185} Early simulations of squirmer suspensions indicated a similar diffusive behaviour for swimmers, but with a diffusion coefficient D_{T} that decreases with increasing volume fraction (∝Φ_{swimmer}^{−1}).^{186} This difference is due to the distinct dynamics of tracers and swimmers, where the former can move only as a result of interactions with other particles, while the latter propel themselves along straight paths, deviating only when interacting with other particles. These interactions can be either hydrodynamic or steric in nature. The hydrodynamic interactions, which give rise to rapidly decaying velocity fluctuations as the particles respond to the evolving configuration of the system, are expected to similarly affect both tracers and swimmers. This short time scale is given approximately by the average time it takes for a swimmer to move a distance equal to its diameter. In contrast, the velocity fluctuations caused by the steric interactions decay over a time scale determined by the average time between collisions, which can be estimated to be ∝Φ_{swimmer}^{−1} from kinetic theory arguments. Clearly, the contributions to the dynamics of the swimmers coming from the hydrodynamic interactions are completely masked by their intrinsic motion. This is the reason why the diffusion coefficients of tracers and swimmers show opposite tendencies.^{187} However, with a careful choice of reference system, it is possible to define an effective diffusion coefficient for the swimmers _{T}, in which the contributions from the particle's self-propulsion have been removed. This effective diffusion is measuring only the contribution from the hydrodynamic interactions. This is done by computing the particle displacements in the rotating body frame and subtracting the average contribution of the swimming motion along ẽ. In this way, we can show the correspondence between tracer diffusion and (effective) swimmer diffusion (see Fig. 20).^{54,56} At high swimmer concentrations, both factors are essentially the same, while at low concentrations the asymmetry of the flow field generated by the swimmers gives rise to a negative correlation in the tracer velocity (as the swimmer moves past) and thus reduces the diffusion coefficient.

Fig. 20 (top) Mean-squared displacements for a suspension of pullers (α = +2) at various volume fractions. Dashed lines show the full mean-squared displacements, as measured in the laboratory frame, while solid lines show the effective mean-squared displacements (after removing the self-propulsion term) in the body frame of the swimmers. (bottom) Effective translational _{T} and rotational D_{R} diffusion coefficients as a function of volume fraction for pullers with α = +1, 2. Adapted from Soft Matter, 9, 4923–4936,^{54} Copyright 2013, with permission from the Royal Society of Chemistry. |

3.8.3 Collective dynamics of micro-swimmers.
It is well known that squirmers in bulk fluid tend to exhibit polar order^{188} provided |α| is small, with the highest degree of order obtained for α = 0. As there are no direct alignment interactions, this collective behaviour is due exclusively to the hydrodynamic interactions among particles. Furthermore, this non-zero polar order is stronger at lower volume fractions, which suggests that it could arise from repeated binary interactions. Indeed, using collision data from dilute simulations, we have computed the transition probabilities for the incoming and outgoing angles of the particles and used them to develop a coarse-grained binary collision model.^{59} Using this simplified representation, we are able to predict the bulk polar order of swimmer suspensions arising from repeated binary collisions. Our results are in good agreement for dilute suspensions of pushers and neutral swimmers but show large discrepancies for pullers, particularly for intermediate α. This can be explained as a breakdown of the binary collision hypothesis, as these systems are known to form transient dynamic clusters in bulk fluid.^{189}

In practice, swimmers usually move in complex environments near boundaries or interfaces, which strongly affects their dynamics and the nature of the hydrodynamic interactions.^{190} Studies indicate that varying the strength and type of confinement can give rise to dramatically different collective behaviours.^{191,192} Simulations of swimmers between flat parallel walls, with a narrow confinement ∼2σ that results in quasi-2D systems, show that neutral squirmers exhibit dense/dilute phase separation, whereas pushers and pullers tend to exhibit a hexagonal phase.^{193} Even though the hydrodynamics are strongly screened by the walls, the near-field hydrodynamic interactions play a crucial role in the particle alignment, and polar ordering is observed only for neutral particles.^{55} Upon increasing the separation of the wall, squirmers tend to accumulate at the boundary in a manner that depends on the swimmer type.^{191}

We have performed simulations of squirmers under confinement between two flat parallel walls and shown that when the wall separation is large enough (≳40σ, i.e., larger than the characteristic size of the propagating flock), weak pullers develop travelling density waves that bounce back and forth between the walls.^{58} A similar travelling wave-like motion has been reported in systems with explicit alignment interactions between the particles,^{194,195} but our results show that hydrodynamic interactions alone are enough to yield this type of collective behaviour. We observe that the strength and type of swimming (pusher vs. puller) are crucial to observe this behaviour but not the degree of global alignment. Furthermore, this behaviour is a manifestation of the bulk behaviour of the suspension and is not due to the presence of walls. This is evidenced by the fact that the velocity of these waves is found to correspond to the sound velocity measured in the bulk fluid, which means that this phenomenon can be interpreted as a pseudo-sound property of the system.

The travelling wave-like motion is illustrated in Fig. 21, where we show simulation snapshots, together with the time evolution of the (local) density ρ, perpendicular polar order parameter Q_{l} and velocity v_{l}, as a function of the height in the channel (computed within parallel slabs of width 2Δ), for a dilute suspension of squirmers (volume fraction of 13%). These simulations are performed in a rectangular simulation box, with periodic boundaries in the x and z directions, but confined between two flat parallel walls along the y direction (at y = 1Δ and 159Δ). Here, for simplicity, the walls are constructed out of (inert) spherical particles of equal radius to the swimmers, which are pinned in place and unable to move or rotate from their initial configuration. Simulations for various α values are performed, with the most interesting results obtained for weak pullers satisfying 0 < α < 1. Fig. 21 shows that the pullers with α = ±0.5 show strong spatiotemporal density fluctuations, first developing two waves travelling in opposite directions, before they merge into a single persistent propagating wave, which bounces back and forth between the top and bottom walls. No such behaviour is observed for the pushers. The onset of this nontrivial collective behaviour is due purely to the differences in the hydrodynamic flow fields generated by the swimmers and their corresponding hydrodynamic interactions, as no direct alignment potential is used.

Fig. 21 Simulation results for pullers and pushers (α = ±0.5) at a volume fraction of 13% and wall separation W/a = 40. (a) Simulation snapshots at T = 11T_{0}, (b) local order of pullers defined as Q_{l} = 〈ẽ_{y}〉_{plane}, (c) local velocity of pullers in the y direction (perpendicular to the walls) defined as _{l} = 〈V_{y}/U_{0}〉_{plane}, (d and e) plane averaged density for pushers and pullers at each height y/W. The density is normalized by the global average value, ρ_{0}, and time and velocity are normalized by T_{0} = W/U_{0} and U_{0}, where W represents the separation of parallel flat walls. The snapshots are at approximately t = 11T_{0}, as marked by dashed lines. Reproduced from Phys. Rev. E, 93, 043114,^{57} Copyright 2016, with permission from the American Physical Society. |

Recently, we have extended this model to include the rotlet term, in addition to the pusher/puller term, and DNS were performed to study the single particle motion of a swimmer near a flat wall.^{60}

For Brownian particles, we have developed a numerical method for consistently implementing the thermal fluctuation in an incompressible host fluid.^{35} The validity of the method was confirmed by comparing our numerical results against the fluctuation–dissipation theorem for a single spherical system. Simulations for concentrated dispersions, as well as for a polymeric chain, were then performed. For the former, we found that the hydrodynamic long-time tail in the velocity auto-correlation function, which is clearly observed for a single-particle dispersion, becomes weaker as the volume fraction is increased. For the latter case, we found very good agreement with the theoretical predictions of the Zimm model.^{35} We have also investigated the rate of rapid Brownian coagulation over a range of volume fractions from φ = 0.003 to φ = 0.1. In the dilute regime the rate of rapid Brownian coagulation is reduced to approximately 0.3–0.5 times the theoretical value predicted by Smoluchowski. This is consistent with previous experimental and theoretical results taking the hydrodynamic effects into account. Because we explicitly include the exact hydrodynamic interactions in our simulations, the fact that we can reproduce this reduction in the coagulation rate is clear evidence for their hydrodynamic origin.^{39}

We have considered non-Brownian sedimenting particles at small Reynolds number (Re ≈ 0.14) over a wide range of volume fractions (0.01 ≤ φ ≤ 0.5). The good agreement between our simulations and other theoretical, experimental and simulation results further illustrates the validity of our method.^{40} Upon increasing the Re, the inertial forces give rise to a drafting–kissing–tumbling (DKT) mechanism that results in a deficit of neighboring particles. This entails significant phenomenological changes and strongly affects the transport properties of the suspension. At small Reynolds number, Re = 0.1, the fluid exhibits a fore-aft symmetry around the particle, whereas at high Reynolds number, Re > 1, an asymmetric behaviour is clearly observed. The theoretical predictions for the settling velocity in the Stokes regime are no longer applicable when these inertial effects are relevant. To extend the predictive capabilities of previously proposed relations, we proposed an empirical formula valid for Re ≤ 10, which captures the rapid decay of the average sedimentation velocity at low volume fractions and the strong influence of the anisotropic local structure.^{43}

To study the rheological behaviour of colloidal dispersions we have extended the SP method to systems under steady or oscillatory shear flow, by using the Lees–Edwards boundary conditions.^{49,51} In this way, all the relevant physical quantities, such as the local and total shear stress, can be directly computed during the simulation. We simulated three systems under simple-shear to test the validity of our model: a spherical particle, a rigid-bead chain, and the collision of two spherical particles. We compared the viscosity predicted by the simulations to theoretical predictions and Stokesian dynamics calculations and found good overall agreement. We are also able to capture the shear-thinning behaviour of concentrated colloidal dispersions.^{51}

Hydrodynamic interactions are propagated by both viscous diffusion and sound propagation. To study the time evolution of the hydrodynamic interactions, and the relative importance of these two mechanisms, we have extended the SP method to perform DNS of particles in compressible media. In an incompressible fluid, with an infinite speed of sound, hydrodynamic interactions propagate instantaneously, followed by a temporal evolution due to viscous diffusion. On the other hand, in a compressible fluid, sound propagates at a finite speed, which affects the time evolution of the hydrodynamic interactions. This effect can be quantified by measuring the ratio of the time scales associated with viscous diffusion and sound propagation,^{62} which defines a compressibility factor. We quantify hydrodynamic interactions for a two-particle system in a compressible fluid through the velocity correlation of the particles. We observe excellent agreement with theoretical predictions.

We have considered multi-phase fluids, consisting of particles dispersed in an A/B fluid mixture, in order to study the attachments/detachments of particles to/from the fluid–fluid interfaces.^{66} We extended the SP method, within a Cahn–Hilliard framework, in order to simulate such multi-phase flows. The method was applied to perform DNS of a flexible fibre, represented as a chain of spherical beads, interacting with a fluidic interface. The fibre can undergo stretching, bending, and twisting deformations. The capillary force, acting at the three-phase contact line, is calculated using a ternary diffuse-interface model. First we validated the fibre deformation and ternary diffuse-interface model against theoretical solutions. We then performed 2D and 3D simulations for the elasto-capillary bending of the fibre by an immersed droplet. We were able to simulate both partial wrapping and complete encapsulation. We found that the fibre curvature increases linearly with the square of the elasto-capillary length, regardless of the degree of structural deformation.^{69}

We have simulated electrohydrodynamic phenomena in charged colloidal dispersions through a suitable extension of the SP method. The time evolution of the colloidal particles, ions, and host fluid were simultaneously computed by solving the Newton–Euler, advection–diffusion, and Navier–Stokes equations, allowing us to fully capture the electrohydrodynamic coupling of the system.^{13} We have calculated the electrophoretic mobilities of charged spherical particles under several conditions. Comparisons with approximate theories showed quantitative agreement for dilute dispersions, without requiring any fitting parameters; however, our predictions exhibit considerable deviations for dense dispersions.^{31} Next, we considered the field-induced anisotropic interactions between like-charged colloidal particles, where the polarization of the electric double layer was explicitly computed under external AC electric fields.^{34} These interactions were found to depend on the magnitude E_{0} and frequency ω of the applied field, as well as the zeta potential, the Debye length, and the relative orientation of the particles. We also performed simulations for systems of six and twelve colloidal particles to study the stability of pearl-chain-like configurations.

We have simulated particle dispersions in viscoelastic media. The shear-thickening behaviour in dilute suspensions was investigated utilizing an Oldroyd-B model.^{64} Here, the fluid elasticity is determined by the Weissenberg number Wi and the relative contribution of the polymer stress that is measured by the viscosity ratio 1 − β = η_{p}/(η_{s} + η_{p}), where η_{p} and η_{s} are the polymer and solvent viscosities, respectively. As 1 − β increases, while the suspension viscosity increases, the growth rate of the normalized polymer stress with Wi is suppressed. An analysis of the flow and the conformation dynamics around a particle reveals that, for large values of 1 − β, the conformation stretching is suppressed because of the modulated flow by the polymer stress. This effect of β on the polymer stress is indicative of a complex coupling between the fluid elasticity and the flow. The experimental shear thickening in a Boger fluid matrix was quantitatively reproduced by multi-particle calculations using an Oldroyd-B fluid matrix up to φ ≤ 1.^{65}

The hydrodynamic interactions of a suspension of self-propelled particles were studied using a modified version of the SP method that represents micro-swimmers as squirmers.^{54} The diffusive behaviour that arises due to the swimming motion of the particles reveals that there are two basic mechanisms at play: the hydrodynamic interactions caused by the squirming motion of the particles and the particle–particle collisions. This dual nature gives rise to two distinct time and length scales, and thus to two distinct diffusion coefficients, which we obtained by a suitable analysis of the swimming motion within the co-rotating reference frame of the particles. These results have allowed us to gain a better understanding of the complex hydrodynamic interactions of self-propelled swimmers. In addition, simulations were performed to study the collective motion of model swimmers in confinement.^{57} We showed that certain swimming mechanisms can lead to travelling wave-like collective motion, even without any direct alignment interactions between the particles. It was also shown that, by varying the swimming mechanism, this collective motion can be suppressed, contrary to the common perception that hydrodynamic effects are completely screened at high volume fractions. From an analysis of bulk systems, it was shown that this travelling wave-like motion, which can be characterized as a pseudo-acoustic mode, is due mainly to the intrinsic swimming property of the particles.

(75) |

(76) |

(77) |

Furthermore, to reduce the memory requirements of the code, we can take advantage of the fact that the three components of the curl are not linearly independent and solve for only two of them (for every given wave vector). Let â = _{k}(∇ × A) = ik × Â be the Fourier transform of the curl of some arbitrary vector field A. We define the following invertible mapping (†: C^{3} → C^{2}) such that

(78) |

The form of the integrator depends on the particular stress tensor s that is used. For the case of a Newtonian fluid, for which ∇·s = η∇^{2}u, the vorticity equation becomes

(79) |

∂_{t}a = a − (t,a(t)) | (80) |

(81) |

a(t_{n} + h) = e^{h}[a(t_{n}) − ^{−1}(1 − e^{−h})(t_{n},a(t_{n}))] | (82) |

= a(t_{n}) + (e^{−h} − 1)(a(t_{n}) + ^{−1}(t_{n},a(t_{n}))) | (83) |

(84) |

(85) |

(86) |

F = Z·U | (87) |

(88) |

(89) |

Simulations at Re ≃ 0 for a single non-spherical particle, constructed as a rigid assembly of non-overlapping beads of equal radius a, are performed with a fixed particle velocity U = (V^{x}, V^{y}, V^{z}, 0, 0, 0). The hydrodynamic forces at steady-state F are measured to compute the friction coefficients. We consider six families of axisymmetric regularly shaped particle agglomerates: (even/odd) v-shaped, h-shaped, hexagonal, and rectangular arrays (see Fig. 22(1–f)). Except for the v- and w-shaped particles, for which there is a weak coupling between translation and rotation, the friction matrices are all diagonal. In all cases, however, the translational friction matrix K^{TT} is diagonal.

Fig. 22 (a–f) Non-spherical regular-shaped agglomerates used in the simulations. (g) Correlation between the friction coefficients obtained from the SP method (DNS) simulations and the exact solutions to the SE (as given by HYDROLIB). Adapted from J. Chem. Phys., 139, 234105,^{52} Copyright 2013, with the permission of AIP publishing. |

We vary the number of beads and relative bead configuration while keeping the same geometrical shape for a total of 84 different geometric configurations. For simplicity, we report the kinetic form factors K, defined in terms of the friction coefficients by

(90) |

We now consider the lubrication interactions between nearby particles moving relative to each other. These are one of the most important and difficult to capture hydrodynamic effects. We computed the squeeze lubrication interaction between two approaching spheres in a finite system. The velocity and pressure distributions are shown in Fig. 23(a). Fig. 23(b) shows the normalized approach velocity for a pair of particles of equal radius (a), versus the normalized gap distance (h/a) between them, obtained using the SP method (black circles) as well as three different theoretical predictions. Our simulation data accurately reproduce the expected crossover from the near-field (h/a<1)^{204} to the far-field (h/a > 1)^{2,203} behaviour around h/a∼ 0.7. It is found that SD^{5} tends to underestimate the mobility in the crossover regime. The validity of the SP method has been further confirmed using other critical tests, such as the friction coefficient of chiral objects^{52} where the coupling between translational and rotational motions becomes important, the volume fraction^{13} and Reynolds number^{43} dependence of the drag coefficient, the motion and stresslet of a sphere under shear flow,^{51} the electrokinetics in charged colloidal systems,^{13,31} and the pair interaction between two squirmers,^{54} among others.

Fig. 23 (a) Snapshot of two approaching colloidal particles of radius a = 4, under a constant attractive force F. The black arrows indicate the velocity field and density map represents the pressure field (a change in color from red to blue corresponds to a change from high to low pressure). (b) The normalised relative velocity of approaching spheres versus the normalised gap distance between them (h/a). A small oscillating behaviour seen in the present simulation data (black circles) is due to the numerical artifacts of using a finite grid size (Δ). The solid curves correspond to three different theoretical predictions: far-field asymptotics computed using the Ewald summed Rotne–Prager–Yamakawa tensor (dotted line),^{2,203} near-field exact solution (dashed line),^{204} and the Stokesian Dynamics predictions (solid line).^{5} Reproduced from Eur. Phys. J. E Soft Matter, 26, 361–368,^{13} Copyright 2008, with kind permission of The European Physical Journal (EPJ). |

- W. B. Russel, W. Russel, D. A. Saville and W. R. Schowalter, Colloidal dispersions, Cambridge University Press, 1991 Search PubMed.
- H. Yamakawa, Modern theory of polymer solutions, Harper & Row, 1971 Search PubMed.
- C. Y. Son and S. Lee, J. Chem. Phys., 2011, 135, 224512 CrossRef PubMed.
- A. Sokolov and I. S. Aranson, Phys. Rev. Lett., 2012, 109, 248109 CrossRef PubMed.
- J. F. Brady and G. Bossis, Annu. Rev. Fluid Mech., 1988, 20, 111–157 CrossRef.
- B. Cichocki, B. Felderhof, K. Hinsen, E. Wajnryb and J. Bl/awzdziewicz, J. Chem. Phys., 1994, 100, 3780–3790 CrossRef CAS.
- A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605–8613 CrossRef CAS.
- H. H. Hu, N. A. Patankar and M. Zhu, J. Comput. Phys., 2001, 169, 427–462 CrossRef CAS.
- R. Glowinski, T.-W. Pan, T. I. Hesla, D. D. Joseph and J. Periaux, J. Comput. Phys., 2001, 169, 363–426 CrossRef CAS.
- A. Ladd and R. Verberg, J. Stat. Phys., 2001, 104, 1191–1251 CrossRef CAS.
- M. Cates, K. Stratford, R. Adhikari, P. Stansell, J. Desplat, I. Pagonabarraga and A. Wagner, J. Phys.: Condens. Matter, 2004, 16, S3903 CrossRef CAS.
- Y. Nakayama and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 036707 CrossRef PubMed.
- Y. Nakayama, K. Kim and R. Yamamoto, Eur. Phys. J. E: Soft Matter Biol. Phys., 2008, 26, 361–368 CrossRef CAS PubMed.
- H. Tanaka and T. Araki, Phys. Rev. Lett., 2000, 85, 1338–1341 CrossRef CAS PubMed.
- T. Kajishima and S. Takiguchi, Int. J. Heat Fluid Flow, 2002, 23, 639–646 CrossRef CAS.
- A. J. Ladd, Phys. Fluids A, 1993, 5, 299–310 CrossRef CAS.
- A. J. Ladd, Phys. Rev. Lett., 1996, 76, 1392 CrossRef CAS PubMed.
- A. J. Ladd, Phys. Fluids, 1997, 9, 491–499 CrossRef CAS.
- A. Ladd, Phys. Rev. Lett., 2002, 88, 048301 CrossRef CAS PubMed.
- N.-Q. Nguyen and A. J. Ladd, J. Fluid Mech., 2005, 525, 73 CrossRef.
- P. Hoogerbrugge and J. Koelman, Europhys. Lett., 1992, 19, 155 CrossRef.
- J. Padding and A. Louis, Phys. Rev. Lett., 2004, 93, 220601 CrossRef CAS PubMed.
- S. H. Lee and R. Kapral, J. Chem. Phys., 2004, 121, 11163–11169 CrossRef CAS PubMed.
- M. Hecht, J. Harting, T. Ihle and H. J. Herrmann, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 011408 CrossRef PubMed.
- I. O. Götze, H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 046705 CrossRef PubMed.
- S. Poblete, A. Wysocki, G. Gompper and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 033314 CrossRef PubMed.
- P. D. Godfrin, N. E. Valadez-Pérez, R. Castaneda-Priego, N. J. Wagner and Y. Liu, Soft Matter, 2014, 10, 5061–5071 RSC.
- S. Jafari, R. Yamamoto and M. Rahnama, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 026702 CrossRef PubMed.
- R. Yamamoto, Phys. Rev. Lett., 2001, 87, 075502 CrossRef CAS PubMed.
- R. Yamamoto, Y. Nakayama and K. Kim, J. Phys.: Condens. Matter, 2004, 16, S1945–S1955 CrossRef CAS.
- K. Kim, Y. Nakayama and R. Yamamoto, Phys. Rev. Lett., 2006, 96, 208302 CrossRef PubMed.
- C. Shih and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 062317 CrossRef PubMed.
- C. Shih, J. J. Molina and R. Yamamoto, Mol. Phys., 2015, 113, 2511–2522 CrossRef CAS.
- C. Shih, J. J. Molina and R. Yamamoto, Soft Matter, 2018, 14, 4520–4529 RSC.
- T. Iwashita, Y. Nakayama and R. Yamamoto, J. Phys. Soc. Jpn., 2008, 77, 074007 CrossRef.
- T. Iwashita and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 031401 CrossRef PubMed.
- T. Iwashita, Y. Nakayama and R. Yamamoto, Prog. Theor. Phys. Suppl., 2009, 178, 86–91 CrossRef CAS.
- R. Yamamoto, K. Kim, Y. Nakayama, K. Miyazaki and D. R. Reichman, J. Phys. Soc. Jpn., 2008, 77, 084804 CrossRef.
- Y. Matsuoka, T. Fukasawa, K. Higashitani and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 051403 CrossRef PubMed.
- A. Hamid, J. J. Molina and R. Yamamoto, Soft Matter, 2013, 9, 10056 RSC.
- A. Hamid and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 022310 CrossRef PubMed.
- A. Hamid and R. Yamamoto, J. Phys. Soc. Jpn., 2013, 82, 024004 CrossRef.
- A. Hamid, J. J. Molina and R. Yamamoto, RSC Adv., 2014, 4, 53681–53693 RSC.
- A. Hamid, J. J. Molina and R. Yamamoto, Mol. Simul., 2015, 41, 968–973 CrossRef CAS.
- M. Shakeel, A. Hamid, A. Ullah, J. J. Molina and R. Yamamoto, J. Phys. Soc. Jpn., 2018, 87, 064402 CrossRef.
- T. Iwashita and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 061402 CrossRef PubMed.
- T. Iwashita, T. Kumagai and R. Yamamoto, Eur. Phys. J. E: Soft Matter Biol. Phys., 2010, 32, 357–363 CrossRef CAS PubMed.
- H. Kobayashi and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 041807 CrossRef PubMed.
- H. Kobayashi and R. Yamamoto, J. Chem. Phys., 2011, 134, 064110 CrossRef PubMed.
- H. Kobayashi and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 051404 CrossRef PubMed.
- J. J. Molina, K. Otomura, H. Shiba, H. Kobayashi, M. Sano and R. Yamamoto, J. Fluid Mech., 2016, 792, 590–619 CrossRef CAS.
- J. J. Molina and R. Yamamoto, J. Chem. Phys., 2013, 139, 234105 CrossRef PubMed.
- N. Oyama, K. Teshigawara, J. J. Molina, R. Yamamoto and T. Taniguchi, Phys. Rev. E, 2018, 97, 032611 CrossRef CAS PubMed.
- J. J. Molina, Y. Nakayama and R. Yamamoto, Soft Matter, 2013, 9, 4923 RSC.
- J. B. Delfau, J. Molina and M. Sano, EPL, 2016, 114, 24001 CrossRef.
- J. J. Molina and R. Yamamoto, Mol. Phys., 2014, 112, 1389–1397 CrossRef CAS.
- N. Oyama, J. J. Molina and R. Yamamoto, Phys. Rev. E, 2016, 93, 043114 CrossRef PubMed.
- N. Oyama, J. J. Molina and R. Yamamoto, J. Phys. Soc. Jpn., 2017, 86, 101008 CrossRef.
- N. Oyama, J. J. Molina and R. Yamamoto, Eur. Phys. J. E: Soft Matter Biol. Phys., 2017, 40, 95 CrossRef PubMed.
- F. Fadda, J. J. Molina and R. Yamamoto, Phys. Rev. E, 2020, 101, 052608 CrossRef CAS PubMed.
- R. Tatsumi and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 066704 CrossRef PubMed.
- R. Tatsumi and R. Yamamoto, Phys. Fluids, 2013, 25, 046101 CrossRef.
- R. Tatsumi and R. Yamamoto, J. Chem. Phys., 2013, 138, 184905 CrossRef PubMed.
- Y. Matsuoka, Y. Nakayama and T. Kajiwara, Soft Matter, 2020, 16, 728–737 RSC.
- Y. Matsuoka, Y. Nakayama and T. Kajiwara, J. Fluid Mech., 2021, 913, A38 CrossRef CAS.
- G. Lecrivain, R. Yamamoto, U. Hampel and T. Taniguchi, Phys. Fluids, 2016, 28, 83301–123305 CrossRef.
- G. Lecrivain, R. Yamamoto, U. Hampel and T. Taniguchi, Phys. Rev. E, 2017, 95, 063107 CrossRef PubMed.
- G. Lecrivain, Y. Kotani, R. Yamamoto, U. Hampel and T. Taniguchi, Phys. Rev. Fluids, 2018, 3, 094002 CrossRef.
- G. Lecrivain, T. B. P. Grein, R. Yamamoto, U. Hampel and T. Taniguchi, J. Comput. Phys., 2020, 409, 109324 CrossRef.
- X. Luo, M. R. Maxey and G. E. Karniadakis, J. Comput. Phys., 2009, 228, 1750–1769 CrossRef.
- M. Fujita and Y. Yamaguchi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 026706 CrossRef PubMed.
- F. Balboa Usabiaga, I. Pagonabarraga and R. Delgado-Buscalioni, J. Comput. Phys., 2012, 235, 701–722 CrossRef.
- H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen and J. M. Yeomans, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 14308–14313 CrossRef CAS PubMed.
- J. V. José and E. J. Saletan, Classical Dynamics: A Contemporary Approach, Cambridge University Press, 1998 Search PubMed.
- R. Glowinski, T. W. Pan, T. I. Hesla and D. D. Joseph, Int. J. Multiphase Flow, 1999, 25, 755–794 CrossRef CAS.
- D. Wan and S. Turek, Int. J. Numer. Methods Fluids, 2006, 51, 531–566 CrossRef.
- C. S. Peskin, Acta Numer., 2002, 11, 479–517 CrossRef.
- M. Uhlmann, J. Comput. Phys., 2005, 209, 448–476 CrossRef.
- P. Atzberger, P. Kramer and C. Peskin, J. Comput. Phys., 2007, 224, 1255–1292 CrossRef.
- S. Lomholt and M. R. Maxey, J. Comput. Phys., 2003, 184, 381–405 CrossRef.
- S. Vincent, J. C. B. De Motta, A. Sarthou, J.-L. Estivalezes, O. Simonin and E. Climent, J. Comput. Phys., 2014, 256, 582–614 CrossRef.
- X. Luo, A. Beskok and G. E. Karniadakis, J. Comput. Phys., 2010, 229, 3828–3847 CrossRef CAS PubMed.
- S. Gallier, E. Lemaire, L. Lobry and F. Peters, J. Comput. Phys., 2014, 256, 367–387 CrossRef.
- J. Lichtenbelt, C. Pathmamanoharan and P. Wiersema, J. Colloid Interface Sci., 1974, 49, 281–285 CrossRef CAS.
- W. Hatton, P. Mcfadyen and A. Smith, J. Chem. Soc., Faraday Trans. 1, 1974, 655–660 RSC.
- K. Higashitani and Y. Matsuno, J. Chem. Eng. Jpn., 1979, 12, 460–465 CrossRef CAS.
- H. Sonntag and K. Strenge, Coagulation kinetics and structure formation, Plenum, New York, 1987 Search PubMed.
- H. Holthoff, S. U. Egelhaaf, M. Borkovec, P. Schurtenberger and H. Sticher, Langmuir, 1996, 12, 5541–5549 CrossRef CAS.
- N. Fuchs, Z. Phys., 1934, 87, 736–743 CrossRef.
- L. Spielman, J. Colloid Interface Sci., 1970, 33, 562–571 CrossRef CAS.
- E. Honig, G. Roebersen and P. Wiersema, J. Colloid Interface Sci., 1971, 36, 97–109 CrossRef CAS.
- K. Higashitani, T. Tanaka and Y. Matsuno, J. Colloid Interface Sci., 1978, 63, 551–560 CrossRef CAS.
- M. Heine and S. Pratsinis, Langmuir, 2007, 23, 9882–9890 CrossRef CAS PubMed.
- G. Urbina-Villalba, J. Toro-Mendoza, A. Lozsán and M. García-Sucre, J. Phys. Chem. B, 2004, 108, 5416–5423 CrossRef CAS.
- T. Trzeciak, A. Podgóski and J. Marijnissen, In 5th World Congress on Particle Technology, Orlando, FL, 2006, p. 202d Search PubMed.
- J. Richardson and W. Zaki, Chem. Eng. Sci., 1954, 3, 65–73 CrossRef CAS.
- J. Garside and M. R. Al-Dibouni, Ind. Eng. Chem. Proc. Des. Dev., 1977, 16, 206–214 CrossRef CAS.
- G. Batchelor, J. Fluid Mech., 1972, 52, 245–268 CrossRef.
- H. Hayakawa and K. Ichiki, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, R3815 CrossRef CAS PubMed.
- J. F. Brady and L. J. Durlofsky, Phys. Fluids, 1988, 31, 717–727 CrossRef CAS.
- R. Di Felice, Int. J. Multiphase Flow, 1999, 25, 559–574 CrossRef CAS.
- H. Nicolai, B. Herzhaft, E. Hinch, L. Oger and E. Guazzelli, Phys. Fluids, 1995, 7, 12–23 CrossRef CAS.
- J. Padding and A. Louis, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 011402 CrossRef CAS.
- E. Climent and M. Maxey, Int. J. Multiphase Flow, 2003, 29, 579–601 CrossRef CAS.
- F. Cunha, G. Abade, A. Sousa and E. Hinch, J. Fluids Eng., 2002, 124, 957–968 CrossRef.
- X. Yin and D. L. Koch, Phys. Fluids, 2007, 19, 093302 CrossRef.
- D. R. Mikulencak and J. F. Morris, J. Fluid Mech., 2004, 520, 215–242 CrossRef.
- A. Onuki, J. Phys.: Condens. Matter, 1997, 9, 6119 CrossRef CAS.
- S. Toh, K. Ohkitani and M. Yamada, Phys. D, 1991, 51, 569–578 CrossRef.
- C. J. Lin, J. H. Peery and W. R. Schowalter, J. Fluid Mech., 1970, 44, 1–17 CrossRef.
- N. Q. Nguyen and A. J. Ladd, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 12 Search PubMed.
- G. B. Jeffery, R. Soc., 1922, 102, 161–179 Search PubMed.
- L. G. Leal and E. J. Hinch, J. Fluid Mech., 1971, 46, 685–703 CrossRef.
- E. J. Hinch and L. G. Leal, J. Fluid Mech., 1972, 52, 683–712 CrossRef.
- R. Yamamoto and A. Onuki, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 58, 3515–3529 CrossRef CAS.
- N. Taccoen, F. m. c. Lequeux, D. Z. Gunes and C. N. Baroud, Phys. Rev. X, 2016, 6, 011010 Search PubMed.
- Y. E. Yu, S. Khodaparast and H. A. Stone, Soft Matter, 2017, 13, 2857–2865 RSC.
- S.-Y. Tan, S. Ata and E. J. Wanless, J. Phys. Chem. B, 2013, 117, 8579–8588 CrossRef CAS PubMed.
- G. Lecrivain, G. Petrucci, M. Rudolph, U. Hampel and R. Yamamoto, Int. J. Multiphase Flow, 2015, 71, 83–93 CrossRef CAS.
- B. D. Johnson and R. C. Cooke, Science, 1981, 213, 209–211 CrossRef CAS.
- K. Stratford, R. Adhikari, I. Pagonabarraga and J. C. Desplat, J. Stat. Phys., 2005, 121, 163–178 CrossRef.
- G. B. Davies, T. Krüger, P. V. Coveney and J. Harting, J. Chem. Phys., 2014, 141, 154902 CrossRef PubMed.
- H. Mehrabian, J. Harting and J. H. Snoeijer, Soft Matter, 2016, 12, 1062–1073 RSC.
- T. Krüger, S. Frijters, F. Günther, B. Kaoui and J. Harting, Eur. Phys. J.: Spec. Top., 2013, 222, 177–198 Search PubMed.
- X.-C. Luu and A. Striolo, J. Phys. Chem. B, 2014, 118, 13737–13743 CrossRef CAS PubMed.
- D. M. Anderson, G. B. McFadden and A. A. Wheeler, Annu. Rev. Fluid Mech., 1998, 30, 139–165 CrossRef.
- H. Shinto, Adv. Powder Technol., 2012, 23, 538–547 CrossRef CAS.
- T. Araki and S. Fukai, Soft Matter, 2015, 11, 3470–3479 RSC.
- T. Araki and H. Tanaka, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 061506 CrossRef PubMed.
- J. Kim, Comput. Methods Appl. Mech. Eng., 2007, 196, 4779–4788 CrossRef.
- F. Boyer, C. Lapuerta, S. Minjeaud, B. Piar and M. Quintard, Transp. Porous Media, 2010, 82, 463–483 CrossRef.
- M. v. Smoluchowski, Z. Phys. Chem., 1918, 93, 129 CrossRef.
- E. Huckel, Phys. Z., 1924, 25, 204–210 Search PubMed.
- D. Henry, Proc. R. Soc. London, Ser. A, 1931, 133, 106–129 CAS.
- R. W. O'Brien and L. R. White, J. Chem. Soc., Faraday Trans. 2, 1978, 1607–1626 RSC.
- S. Levine and G. H. Neale, J. Colloid Interface Sci., 1974, 47, 520–529 CrossRef.
- H. Ohshima, J. Colloid Interface Sci., 1997, 188, 481–485 CrossRef CAS.
- J.-P. Hansen and H. Löwen, Annu. Rev. Phys. Chem., 2000, 51, 209–242 CrossRef CAS PubMed.
- J.-L. Barrat and J.-P. Hansen, Basic concepts for simple and complex liquids, Cambridge University Press, 2003 Search PubMed.
- EPAPS Document No. E-PRLTAO-96-051621, http://www.aip.org/pubservs/epaps.html.
- F. Carrique, F. J. Arroyo, M. L. Jiménez and Á. V. Delgado, J. Phys. Chem. B, 2003, 107, 3199–3206 CrossRef CAS.
- V. Lobaskin, B. Dünweg, M. Medebach, T. Palberg and C. Holm, Phys. Rev. Lett., 2007, 98, 176105 CrossRef PubMed.
- T. Palberg, M. Medebach, N. Garbow, M. Evers, A. B. Fontecha, H. Reiber and E. Bartsch, J. Phys.: Condens. Matter, 2004, 16, S4039 CrossRef CAS.
- D. Boger, J. Non-Newtonian Fluid Mech., 1977, 3, 87–91 CrossRef CAS.
- D. F. James, Annu. Rev. Fluid Mech., 2009, 41, 129–142 CrossRef.
- I. E. Zarraga, D. A. Hill and D. T. Leighton, J. Rheol., 2001, 45, 1065–1084 CrossRef CAS.
- R. Scirocco, J. Vermant and J. Mewis, J. Rheol., 2005, 49, 551–567 CrossRef CAS.
- S.-C. Dai, F. Qi and R. I. Tanner, J. Rheol., 2014, 58, 183–198 CrossRef CAS.
- R. I. Tanner, J. Non-Newtonian Fluid Mech., 2015, 222, 18–23 CrossRef CAS.
- M. Yang and E. S. G. Shaqfeh, J. Rheol., 2018, 62, 1379–1396 CrossRef CAS.
- S. Dai and R. I. Tanner, Rheol. Acta, 2020, 59, 477–486 CrossRef CAS.
- W. R. Hwang, M. A. Hulsen and H. E. Meijer, J. Non-Newtonian Fluid Mech., 2004, 121, 15–33 CrossRef CAS.
- G. D'Avino, F. Greco, M. A. Hulsen and P. L. Maffettone, J. Rheol., 2013, 57, 813–839 CrossRef.
- A. Vázquez-Quesada, P. Español, R. I. Tanner and M. Ellero, J. Fluid Mech., 2019, 880, 1070–1094 CrossRef.
- N. A. Patankar and H. H. Hu, J. Non-Newtonian Fluid Mech., 2001, 96, 427–443 CrossRef CAS.
- M. Ahamadi and O. G. Harlen, J. Comput. Phys., 2008, 227, 7543–7560 CrossRef.
- M. Ahamadi and O. G. Harlen, J. Non-Newtonian Fluid Mech., 2010, 165, 281–291 CrossRef CAS.
- Y. J. Choi, M. A. Hulsen and H. E. Meijer, J. Non-Newtonian Fluid Mech., 2010, 165, 607–624 CrossRef CAS.
- Y. J. Choi and M. A. Hulsen, J. Non-Newtonian Fluid Mech., 2012, 175-176, 89–103 CrossRef CAS.
- N. O. Jaensson, M. A. Hulsen and P. D. Anderson, J. Non-Newtonian Fluid Mech., 2016, 235, 125–142 CrossRef CAS.
- N. O. Jaensson, M. A. Hulsen and P. D. Anderson, J. Non-Newtonian Fluid Mech., 2015, 225, 70–85 CrossRef CAS.
- M. Yang, S. Krishnan and E. S. Shaqfeh, J. Non-Newtonian Fluid Mech., 2016, 233, 181–197 CrossRef CAS.
- M. Yang and E. S. G. Shaqfeh, J. Rheol., 2018, 62, 1363–1377 CrossRef CAS.
- W. R. Hwang and M. A. Hulsen, Macromol. Mater. Eng., 2011, 296, 321–330 CrossRef CAS.
- R. Pasquino, G. D'Avino, P. L. Maffettone, F. Greco and N. Grizzuti, J. Non-Newtonian Fluid Mech., 2014, 203, 1–8 CrossRef CAS.
- S. Krishnan, E. S. Shaqfeh and G. Iaccarino, J. Comput. Phys., 2017, 338, 313–338 CrossRef.
- Y. K. Lee and K. H. Ahn, J. Non-Newtonian Fluid Mech., 2017, 244, 75–84 CrossRef CAS.
- C. Fernandes, S. A. Faroughi, O. S. Carneiro, J. M. Nóbrega and G. H. McKinley, J. Non-Newtonian Fluid Mech., 2019, 266, 80–94 CrossRef CAS.
- R. B. Bird, C. F. Curtiss, R. C. Armstrong and O. Hassager, Dynamics of Polymeric Liquids, Volume 2: Kinetic Theory, Wiley, 1987 Search PubMed.
- R. G. Larson, Constitutive equations for polymer melts and solutions: Butterworths series in chemical engineering, Butterworth-Heinemann, Boston, 1988 Search PubMed.
- C. W. Macosko, Rheology: Principles, Measurements, and Applications, VCH Publishes, 1994 Search PubMed.
- R. I. Tanner, Engineering rheology, OUP, Oxford, 2000, vol. 52 Search PubMed.
- R. G. Owens and T. N. Phillips, Computational rheology, Imperial College Press, London, 2002, vol. 14 Search PubMed.
- M. Pasquali and L. Scriven, J. Non-Newtonian Fluid Mech., 2002, 108, 363–409 CrossRef CAS.
- M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189 CrossRef CAS.
- M. J. Lighthill, Annu. Rev. Fluid Mech., 1969, 1, 413 CrossRef.
- J. R. Blake, J. Fluid Mech., 1971, 46, 199–208 CrossRef.
- M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper and H. Rieger, Nat. Rev. Phys., 2020, 2, 181–199 CrossRef.
- S. R. Keller and T. Y. Wu, J. Fluid Mech., 1977, 80, 259 CrossRef.
- K. Kyoya, D. Matsunaga, Y. Imai, T. Omori and T. Ishikawa, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 1–6 CrossRef PubMed.
- O. S. Pak and E. Lauga, J. Eng. Math., 2014, 88, 1–28 CrossRef.
- X. L. Wu and A. Libchaber, Phys. Rev. Lett., 2000, 84, 3017–3020 CrossRef CAS PubMed.
- K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci and R. E. Goldstein, Phys. Rev. Lett., 2009, 103, 1–4 CrossRef PubMed.
- G. L. Miño, J. Dunstan, A. Rousselet, E. Clément and R. Soto, J. Fluid Mech., 2013, 729, 423–444 CrossRef.
- A. Jepson, V. A. Martinez, J. Schwarz-Linek, A. Morozov and W. C. Poon, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 3–7 CrossRef PubMed.
- T. Ishikawa and T. J. Pedley, J. Fluid Mech., 2007, 588, 437–462 CrossRef.
- T. Ishikawa, J. T. Locsei and T. J. Pedley, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 1–15 CrossRef PubMed.
- A. A. Evans, T. Ishikawa, T. Yamaguchi and E. Lauga, Phys. Fluids, 2011, 23, 111702 CrossRef.
- F. Alarcón and I. Pagonabarraga, J. Mol. Liq., 2013, 185, 56–61 CrossRef.
- C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006 CrossRef.
- G. J. Li and A. M. Ardekani, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 1–12 Search PubMed.
- L. Zhu, E. Lauga and L. Brandt, J. Fluid Mech., 2013, 726, 285 CrossRef.
- A. Zöttl and H. Stark, Phys. Rev. Lett., 2014, 112, 1–5 CrossRef PubMed.
- H. Chaté, F. Ginelli, G. Grégoire, F. Peruani and F. Raynaud, Eur. Phys. J. B, 2008, 64, 451–456 CrossRef.
- T. Ohta and S. Yamanaka, Prog. Theor. Exp. Phys., 2014, 2014, 1–7 Search PubMed.
- KAPSEL Home Page, https://kapsel-dns.com.
- S. M. Cox and P. C. Matthews, J. Comput. Phys., 2002, 176, 430–455 CrossRef CAS.
- M. Hochbruck and A. Ostermann, Acta Numer., 2010, 19, 209–286 CrossRef.
- M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford Science Publications, 1987 Search PubMed.
- S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications, Dover Publications, 2005 Search PubMed.
- K. Hinsen, Comput. Phys. Commun., 1995, 88, 327–340 CrossRef CAS.
- B. Cichocki and K. Hinsen, Phys. Fluids, 1995, 7, 285–291 CrossRef CAS.
- C. Beenakker, J. Chem. Phys., 1986, 85, 1581–1582 CrossRef CAS.
- D. J. Jeffrey and Y. Onishi, J. Fluid Mech., 1984, 139, 261–290 CrossRef.
- OCTA Home Page, http://octa.jp.

This journal is © The Royal Society of Chemistry 2021 |