1 Introduction
Colloidal dispersions exhibit a rich variety of phenomena and have been providing time-honoured problems in the field of physical chemistry for the past century (Wagner & Brady Reference Wagner and Brady2009; Mewis & Wagner Reference Mewis and Wagner2012), starting with the pioneering work of von Smoluchowski (Reference von Smoluchowski1906) and Einstein (Reference Einstein1911), among others. In particular, the rheometry or viscosity measurement has proved to be one of the most challenging issues. Theoretical descriptions based on hydrodynamic theories have been continuously improved upon, but they are still valid only at low concentrations. In addition, new phenomena have continued to be discovered, such as the effect of frictional interactions on the discontinuous shear thickening (Seto et al.
Reference Seto, Mari, Morris and Denn2013; Mari et al.
Reference Mari, Seto, Morris and Denn2014) and the emergence of irreversibility in the configurations of a sheared suspension at low Reynolds (
$Re$
) number (Pine et al.
Reference Pine, Gollub, Brady and Leshansky2005; Corte et al.
Reference Corte, Gerbode, Man and Pine2009). Microrheology, which has found extensive applications outside colloidal science (Mason, Bibette & Weitz Reference Mason, Bibette and Weitz1995; Mizuno et al.
Reference Mizuno, Head, MacKintosh and Schmidt2008; Guo et al.
Reference Guo, Ehrlicher, Jensen, Renz, Moore, Goldman, Lippincott-Schwartz, Mackintosh and Weitz2014; Head et al.
Reference Head, Ikebe, Nakamasu, Zhang, Villaruz, Kinoshita, Ando and Mizuno2014), is now being used to probe the dynamics of active matter systems (Marchetti et al.
Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013), which show anomalous viscosities due to their inherent swimming motion (Sokolov & Aranson Reference Sokolov and Aranson2009; Mussler et al.
Reference Mussler, Rafaï, Peyla and Wagner2013).
When dealing with many-body interactions between colloidal particles, particularly when hydrodynamic effects are significant, theoretical predictions are limited to simple systems at low volume fractions. In such cases, numerical simulations have proved to be a powerful tool to probe the dynamics of colloidal particle suspensions. However, traditional molecular dynamics simulations are insufficient, as they are unable to probe the necessary length and time scales. Thus, a simplified approach is necessary to simulate such systems. The most widely accepted method for solving this issue is Stokesian dynamics (SD) (Brady & Bossis Reference Brady and Bossis1988), along with recently developed approximations such as fast lubrication dynamics (Bybee Reference Bybee2009). In this framework, the Stokes equation is solved with boundary conditions determined by the shape of the rigid particles. Although such methods provide a high degree of accuracy, and valuable information has been obtained using them (Foss & Brady Reference Foss and Brady2000; Sierou & Brady Reference Sierou and Brady2004; Seto, Botet & Briesen Reference Seto, Botet and Briesen2011; Kumar & Higdon Reference Kumar and Higdon2012), they are restricted to simple host fluids at zero
$Re$
. To consider complex host solvents at arbitrary
$Re$
, various mesoscopic approaches have been invented in order to take the hydrodynamic effects into account: lattice Boltzmann method (LBM) (Ladd Reference Ladd1993; Succi Reference Succi2001; Cates et al.
Reference Cates, Stratford, Adhikari, Stansell, Desplat, Pagonabarraga and Wagner2004), fluid–particle dynamics (Tanaka & Araki Reference Tanaka and Araki2000), smoothed profile method (SPM) (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005; Nakayama, Kim & Yamamoto Reference Nakayama, Kim and Yamamoto2008), distributed-Lagrange-multiplier/fictitious domain method (Glowinski et al.
Reference Glowinski, Pan, Hesla and Joseph1998, Reference Glowinski, Pan, Hesla, Joseph and Périaux2000), force coupling method (FCM) (Maxey & Patel Reference Maxey and Patel2001; Yeo & Maxey Reference Yeo and Maxey2010), and particle-based hydrodynamic simulation approaches, such as the multi-particle collision method (MPC) (Ji et al.
Reference Ji, Jiang, Winkler and Gompper2011; Poblete et al.
Reference Poblete, Wysocki, Gompper and Winkler2014), among others. The boundary integral or boundary element method, originally developed by Youngren & Acrivos (Reference Youngren and Acrivos1974), Kim & Karrila (Reference Kim and Karrila2005) and Muldowney & Higdon (Reference Muldowney and Higdon2006), has also been used extensively within the fluid dynamics community, though mostly for low-Reynolds-number flow problems. In this paper, the methodology and the quantitative validity of the SPM for rheological calculations will be considered. In particular, we calculate the viscosity of the system by measuring its response to an applied external shear flow. The standard methodologies for implementing such a flow in a simulation include the moving-wall boundary condition and the Lees–Edwards boundary condition. In practice, the former, which explicitly includes the walls driving the flow, can be employed for all the simulation methods (this seems to be the preferred approach when using LBM). However, this approach introduces significant boundary effects, which can be problematic if one intends to study the dynamics of bulk systems. To remove these effects, one must use suitable periodic boundaries. For systems under simple-shear flow, these are the Lees–Edwards boundary conditions (Lees & Edwards Reference Lees and Edwards1972).
In this paper, an improved formulation for shear-flow simulations using the SPM, along with demonstrative simulation results intended to establish the validity of the method, is presented. The organization of the paper is as follows. In § 2, the implementation of the SPM using the Lees–Edwards boundary conditions is presented. In particular, we show how to obtain all components of the local stress fluctuation. In § 3, the quantitative validity of the formulation will be examined using rheological simulations (under steady shear flow) for three example cases, which have been solved theoretically (with various degrees of approximation) in the literature: (1) a single spherical particle, (2) a rigid bead chain, and (3) a collision between two spherical particles. For these three cases, the validity of the viscosity obtained by SPM will be quantitatively examined, by comparing them with the available theoretical descriptions and also with the SD results. In addition, we have also considered the shear-thinning behaviour of concentrated spherical particle dispersions. In § 4, we summarize our work and present our conclusions.
2 Model and methods
2.1 The smoothed profile method for shear-flow simulations
To simulate diverse types of colloidal dispersions under shear flow, we further extend the simulation method developed by two of the authors (Kobayashi & Yamamoto Reference Kobayashi and Yamamoto2011). In this method, the shear flow is driven by an imposed coordinate flow (i.e. a moving grid). The dispersion is simulated by using the SPM, which is a direct numerical simulation method for solving the motion of the host fluid and the particles simultaneously. This is accomplished by the introduction of a diffuse interface between the two phases (solid/fluid), which replaces the original sharp interface, and allows us to define all field variables over the entire computational domain. When solving the Navier–Stokes equation the system is then considered as a single-component fluid, without any boundaries. The effect of particles is introduced through an additional forcing term in the Navier–Stokes equation. In this sense our method is similar to the distributed-Lagrange-multiplier method (Glowinski et al. Reference Glowinski, Pan, Hesla and Joseph1998), the fluid particle dynamics method (Tanaka & Araki Reference Tanaka and Araki2000), and the force coupling method (Maxey & Patel Reference Maxey and Patel2001); the only difference lies in how the coupling between the particles and the fluid is achieved. Since the system lacks periodicity under the Lees–Edwards conditions, the solution method employed for previous SPM studies (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005; Nakayama et al. Reference Nakayama, Kim and Yamamoto2008; Molina & Yamamoto Reference Molina and Yamamoto2013), which used the fast Fourier transform (FFT) to efficiently solve the equations of motion on a discrete mesh, cannot be directly applied. Thus, one must construct a new generalized (time-dependent) coordinate system – ‘oblique coordinate’ – in which the system is periodic in all directions (Onuki Reference Onuki1997; Kobayashi & Yamamoto Reference Kobayashi and Yamamoto2011). Alternative methods have been used to study such systems, we note, for example, the work by Hwang, Hulsen & Meijer (Reference Hwang, Hulsen and Meijer2004a ,Reference Hwang, Hulsen and Meijer b ), who used Glowinski’s distributed Lagrangian method to study the dynamics of inertialess particles under simple shear in 2D, and the study of Villone et al. (Reference Villone, D’Avino, Hulsen, Greco and Maffettone2011) into the particle migration in a viscoelastic fluid under Poiseuille flow. In contrast to our method, these approaches require a complex treatment of the solid–fluid interface (by the use of a finite element method). The virtue of the SPM is in the introduction of a diffuse interface between the two phases, which greatly simplifies the calculations.
With the introduction of a smooth interface of thickness
${\it\zeta}$
between the host fluid and the particles, the incompressible Navier–Stokes equation is modified to properly consider the dynamics of the whole system, consistent with the particle equations of motion. In the absence of shear flow, this modified Navier–Stokes equation takes the form (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005)
$$\begin{eqnarray}\displaystyle (\partial _{t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}={\it\rho}^{-1}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\sigma}+{\it\phi}\boldsymbol{f}_{p}, & & \displaystyle\end{eqnarray}$$
with
$\boldsymbol{u}$
the total velocity field,
${\it\rho}$
the dispersion density,
${\bf\sigma}$
the viscous stress tensor, and
${\it\phi}\boldsymbol{f}_{p}$
the body force required to satisfy the rigid body constraints. The incompressibility condition is applied to the total velocity
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}$
; this also guarantees the solid–fluid impermeability (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005). The total velocity field is defined in terms of the particle and fluid velocity fields,
$\boldsymbol{u}_{p}$
and
$\boldsymbol{u}_{f}$
, respectively, as
$\boldsymbol{u}(\boldsymbol{x})=(1-{\it\phi}(\boldsymbol{x}))\boldsymbol{u}_{f}(\boldsymbol{x})+{\it\phi}(\boldsymbol{x})\boldsymbol{u}_{p}(\boldsymbol{x})$
, with
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\phi}(\boldsymbol{x})=\mathop{\sum }_{i}^{N}{\it\phi}_{i}(\boldsymbol{x};\boldsymbol{R}_{i},\unicode[STIX]{x1D618}_{i}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\phi}(\boldsymbol{x})\boldsymbol{u}_{p}(\boldsymbol{x})=\mathop{\sum }_{i}^{N}{\it\phi}_{i}(\boldsymbol{x};\boldsymbol{R}_{i},\unicode[STIX]{x1D618}_{i})(\boldsymbol{V}_{i}+{\it\bf\Omega}_{i}\times \boldsymbol{r}_{i}), & \displaystyle\end{eqnarray}$$
$N$
is the total number of particles,
${\it\phi}_{i}$
is the smooth profile function of particle
$i$
, with position, orientation matrix, velocity and angular velocity given by
$\boldsymbol{R}_{i}$
,
$\unicode[STIX]{x1D618}_{i}$
,
$\boldsymbol{V}_{i}$
and
${\it\bf\Omega}_{i}$
, respectively, and
$\boldsymbol{r}_{i}=\boldsymbol{x}-\boldsymbol{R}_{i}$
. The dispersion density is defined as
${\it\rho}=(1-{\it\phi}){\it\rho}_{f}+{\it\phi}{\it\rho}_{p}$
, with
${\it\rho}_{f}$
and
${\it\rho}_{p}$
the fluid and particle densities, respectively. Throughout this work, we will assume buoyancy-free particles with the same density as the fluid, such that
${\it\rho}={\it\rho}_{p}={\it\rho}_{f}$
. The particle phase field
${\it\phi}_{i}(\boldsymbol{x})$
is defined such that
${\it\phi}=1$
inside the particle domain,
${\it\phi}=0$
in the fluid domain, and
$0\leqslant {\it\phi}\leqslant 1$
within the interfacial regions. For a spherical particle of radius
$a$
, we thus require 
$$\begin{eqnarray}\displaystyle & & \displaystyle {\it\phi}_{i}(\boldsymbol{x}):\left\{\begin{array}{@{}ll@{}}1.0,\quad & |\boldsymbol{x}-\boldsymbol{R}_{i}|<a-{\it\zeta}/2,\\ 0.0<{\it\phi}_{i}<1.0,\quad & a-{\it\zeta}/2\leqslant |\boldsymbol{x}-\boldsymbol{R}_{i}|\leqslant a+{\it\zeta}/2,\\ 0.0,\quad & |\boldsymbol{x}-\boldsymbol{R}_{i}|>a+{\it\zeta}/2.\end{array}\right.\end{eqnarray}$$
Explicit expressions for
${\it\phi}_{i}$
of spherical particles can be found in the literature (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005; Nakayama et al.
Reference Nakayama, Kim and Yamamoto2008), as well as an extension to handle arbitrary rigid bodies (Molina & Yamamoto Reference Molina and Yamamoto2013). For simplicity, we assume a Newtonian fluid of viscosity
${\it\eta}_{f}$
, so that the Newtonian stress tensor takes the form
$$\begin{eqnarray}\displaystyle {\bf\sigma}=-p\unicode[STIX]{x1D644}+{\it\eta}_{f}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}+^{t}(\boldsymbol{{\rm\nabla}}\boldsymbol{u})), & & \displaystyle\end{eqnarray}$$
where
$p$
is the pressure (
$\unicode[STIX]{x1D644}$
the unit tensor). By letting the viscous stress tensor act on the entire domain, we indirectly enforce the no-slip boundary condition at the particle surface, as it will tend to reduce any difference between
$\boldsymbol{u}_{f}$
and
$\boldsymbol{u}_{p}$
over the interface region (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005; Luo, Maxey & Karniadakis Reference Luo, Maxey and Karniadakis2009).
For simulations under shear, using an oblique coordinate system with Lees–Edwards boundary conditions, care must be taken when writing down the Navier–Stokes equation (2.1), since the basis vectors, as well as the grid points on which the field variables are defined, are time-dependent. Under the assumption that the fluid reacts instantaneously to the imposed flow, the (incompressible) Navier–Stokes equation (2.1) can be written as (see appendix A)
$$\begin{eqnarray}\displaystyle & \displaystyle (\partial _{t}+\hat{{\it\xi}}^{{\it\nu}}\hat{{\rm\nabla}}_{{\it\nu}})\hat{{\it\xi}}^{{\it\mu}}={\it\rho}^{-1}\hat{{\rm\nabla}}_{{\it\nu}}\hat{{\it\sigma}}^{{\it\nu}{\it\mu}}+\hat{{\it\phi}}\hat{f}_{p}^{{\it\mu}}-2\dot{{\it\gamma}}(t)\hat{{\it\xi}}^{2}{\it\delta}^{{\it\mu},1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{{\rm\nabla}}_{{\it\mu}}\hat{u} ^{{\it\mu}}=\hat{{\rm\nabla}}_{{\it\mu}}\hat{{\it\xi}}^{{\it\mu}}=0, & \displaystyle\end{eqnarray}$$
$\hat{{\it\xi}}^{{\it\mu}}$
is the contravariant component of the peculiar fluid velocity, i.e. the fluid velocity relative to the shear-flow velocity
$\boldsymbol{U}=\dot{{\it\gamma}}y\boldsymbol{e}_{x}$
$$\begin{eqnarray}\displaystyle {\bf\xi}=\boldsymbol{u}-\boldsymbol{U} & & \displaystyle\end{eqnarray}$$
with carets (
$\hat{\cdot }$
) used to specify tensor components in the oblique frame. The Newtonian stress tensor can be expressed in oblique coordinates as
${\bf\sigma}=\hat{{\it\sigma}}^{{\it\mu}{\it\nu}}\hat{\boldsymbol{E}}_{{\it\mu}}\hat{\boldsymbol{E}}_{{\it\nu}}$
$$\begin{eqnarray}\displaystyle \hat{{\it\sigma}}^{{\it\mu}{\it\nu}}=-\unicode[STIX]{x1D60E}^{{\it\mu}{\it\nu}}\hat{p}+{\it\eta}_{f}(\unicode[STIX]{x1D60E}^{{\it\mu}{\it\gamma}}\hat{{\rm\nabla}}_{{\it\gamma}}\hat{u} ^{{\it\nu}}+\unicode[STIX]{x1D60E}^{{\it\nu}{\it\gamma}}\hat{{\rm\nabla}}_{{\it\gamma}}\hat{u} ^{{\it\mu}}), & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D60E}^{{\it\mu}{\it\gamma}}$
is the metric tensor (A 3). The force due to the viscous stresses in the fluid can be evaluated solely in terms of the peculiar velocity, since
$$\begin{eqnarray}\displaystyle \hat{{\rm\nabla}}_{{\it\nu}}\hat{{\it\sigma}}^{{\it\nu}{\it\mu}}=-\unicode[STIX]{x1D60E}^{{\it\nu}{\it\mu}}\hat{{\rm\nabla}}_{{\it\nu}}\hat{p}+{\it\eta}_{f}\,\unicode[STIX]{x1D60E}^{{\it\nu}{\it\gamma}}\hat{{\rm\nabla}}_{{\it\nu}}\hat{{\rm\nabla}}_{{\it\gamma}}\hat{{\it\xi}}^{{\it\mu}}, & & \displaystyle\end{eqnarray}$$
where we have used the fact that the velocity field is solenoidal (2.5b
), and the Laplacian of the coordinate flow is zero (
${\rm\nabla}^{2}U^{{\it\mu}}=0$
). The fluid equations of motion, (2.5)–(2.8), are discretized using a Fourier spectral scheme in space and a first-order Euler scheme in time. To simultaneously solve the fluid equations of motion alongside the particle equations of motion, we use the same fractional-step algorithm introduced by Nakayama et al. (Reference Nakayama, Kim and Yamamoto2008) to handle spherical particles, and later extended by Molina & Yamamoto (Reference Molina and Yamamoto2013) for arbitrary rigid bodies.
The solution procedure can be summarized as follows (see appendix B): (1) we solve for the viscous advection and diffusion terms, by integrating (2.5) in oblique space over the whole domain (without the body-force term
${\it\phi}\boldsymbol{f}_{p}$
). We simultaneously update the particle positions and orientations. (2) The updated fluid velocity is then transformed to rectangular coordinates, and the momentum exchange over the particle domain is computed, in order to determine the hydrodynamic force. (3) The hydrodynamic force (together with any external or inter-particle force) is used to update the particle velocities. (4) Finally, the particle rigidity (using the new particle velocities) is enforced by applying the body-force term
${\it\phi}\boldsymbol{f}_{p}$
, to obtain the updated velocity field over the whole computational domain. We note that the incompressibility condition is enforced twice, during the first and final steps, when the fluid velocity is being updated. The last step, which enforces the rigidity constraint, is crucial to achieve an accurate coupling between the fluid and the particles.
The body-force term acts as a penalty function that makes the velocity field match the rigid body velocities inside the particle domain. Naively, one thinks that a smaller time step leads to smaller error. However, the impulse that is added to the particle region as the constraint is imposed will in turn create an impulse in the fluid. The Stokes layer over which this impulse appears scales as
$\sqrt{{\it\nu}{\rm\Delta}t}$
(
${\it\nu}$
the kinematic viscosity), i.e. it will become thinner as the time step is decreased. If the thickness of this Stokes layer falls below the resolution of the grid, the error will increase (Luo et al.
Reference Luo, Maxey and Karniadakis2009). This situation, in which a smaller time step gives a better control of the rigidity constraint, but a thinner and less resolved Stokes layer, results in a non-monotonic error as a function of the time step. A detailed error analysis of the SPM has been performed by Luo et al., they find that the optimum value for the time step
${\rm\Delta}t$
depends on the specific time-stepping scheme that is used, but in general it will scale as
${\rm\Delta}t\sim {\it\zeta}^{2}/{\it\nu}$
. As expected, a smaller value of
${\it\zeta}$
will provide a more accurate representation, but results in a smaller optimum time step. For simplicity, we do not use the precise time step which minimizes this time-stepping error, instead, we use the stability condition given by the momentum diffusion term, which gives
${\rm\Delta}t=1/{\it\nu}K_{max}^{2}$
(
$K_{max}$
the largest Fourier mode in our pseudo-spectral scheme), which is independent of
${\it\zeta}$
. For typical values of
${\it\zeta}=1,2$
, there is no significant difference between the two time steps (the error is a function of
$\sqrt{{\it\nu}{\rm\Delta}t}/{\it\zeta}$
), and a similar degree of accuracy is obtained (Nakayama & Yamamoto Reference Nakayama and Yamamoto2005).
Among the alternative simulation methods that have been proposed to model homogeneous shear simulations, the SPM with the Lees–Edwards boundary conditions (Kobayashi & Yamamoto Reference Kobayashi and Yamamoto2011) closely resembles that used by Rogallo (Reference Rogallo1981) and Brucker et al. (Reference Brucker, Isaza, Vaithianathan and Collins2007) to study turbulent flows, and Onuki’s investigation of phase separating systems (Onuki Reference Onuki1997). In these cases, the fluid equations are solved on a grid that is moving with the shear flow, but in our case we have added an additional forcing term
${\it\phi}\boldsymbol{f}_{p}$
to treat particles dispersed in the fluid. The standard approach requires the remeshing of the grid at some specified strain rate, in order to maintain numerical stability (Onuki Reference Onuki1997). In our case, we reset the grid (
${\it\gamma}=0$
) whenever
${\it\gamma}(t)=\dot{{\it\gamma}}t=1$
, since, at this point in time, the stretched oblique grid overlaps with the ortho-normal Cartesian grid. Unfortunately, such remeshing introduces discontinuities which lead to a loss of turbulent kinetic energy and rate of dissipation (Brucker et al.
Reference Brucker, Isaza, Vaithianathan and Collins2007). These discontinuities are not an issue for the shear rates (Reynolds numbers) considered in this work, but they can become drastic at very high shear rates, such as those used by Rogallo to study turbulent flows. Fortunately, Brucker et al. have provided a clever fix to this problem. They noticed that the effect of the moving frame could be accounted for by introducing a phase shift into the Fourier transforms themselves, while still maintaining a fixed orthogonal mesh. The drawback to this method is the fact that standard FFT routines cannot be used, but it is entirely compatible with the current SPM framework.
Finally, we note that there is an alternative method to perform shear-flow simulations with the SPM, which uses an external force to maintain a linear zig-zag profile in the fluid (
$v_{x}(y=\pm L_{y}/4)=\pm \dot{{\it\gamma}}L_{y}/4$
, with
$L_{y}$
the length of the box along the shear-gradient direction) (Iwashita & Yamamoto Reference Iwashita and Yamamoto2009). As the zig-zag flow matches the periodic boundary conditions, there is no need to use a time-dependent oblique coordinate frame. However, the use of this zig-zag profile presents several difficulties. First, the shear rate cannot be controlled precisely and, more importantly, there are artificial kinks in the flow field, where the sign of the velocity gradient changes sign. This can lead to unwanted situations when the particles cross the kinks in the flow. In particular, if a particle is placed at the position of the kinks
$y=\pm L_{y}/4$
, it will exhibit translation with no net rotation. For non-spherical particles this can result in a considerable increase of the shear stresses in the fluid. At moderate concentrations, we can expect the relative effect of the particles near the kink regions to be reduced; however, for non-spherical particles or very dense suspensions, it is not immediately clear what the effect is on the dynamics and structure of the system. In addition, one cannot directly obtain the zero-wavevector transport coefficients by applying this finite-wavelength perturbation to the system. The present approach, based on the Lees–Edwards formulation, is free of such problems.
2.2 Stress calculation
The rheometry of suspensions is one of the most crucial parts of colloidal science (Mewis & Wagner Reference Mewis and Wagner2012). Because of the temporally slow and spatially large degrees of freedom characterizing colloidal suspensions, the shear stress depends strongly on the imposed shear-flow velocity. For the suspension rheology of non-deformable colloidal particles, both the inter-particle and hydrodynamic interactions contribute to the macroscopic stress behaviour (Batchelor Reference Batchelor1977; Brady & Bossis Reference Brady and Bossis1988; Brady Reference Brady1993; Bender & Wagner Reference Bender and Wagner1995). The average stress tensor of the dispersion
$\langle {\it\bf\Sigma}\rangle$
is given as
$$\begin{eqnarray}\displaystyle \langle {\it\bf\Sigma}\rangle =\langle {\bf\sigma}_{f}\rangle +\langle \unicode[STIX]{x1D668}\rangle , & & \displaystyle\end{eqnarray}$$
where, by assuming ergodicity, ensemble averages
$\langle \cdots \rangle$
are averages over volume and time, with
$\langle {\bf\sigma}_{f}\rangle$
the bulk contribution due to the fluid
$$\begin{eqnarray}\displaystyle \langle {\bf\sigma}_{f}\rangle \equiv -\langle p\rangle \unicode[STIX]{x1D644}+2{\it\eta}_{f}\langle \unicode[STIX]{x1D65A}\rangle , & & \displaystyle\end{eqnarray}$$
$\langle \unicode[STIX]{x1D65A}\rangle$
the bulk rate-of-strain tensor, and
$\langle \unicode[STIX]{x1D668}\rangle$
the particle contribution to the stress. The dispersion viscosity
${\it\eta}$
can then be obtained from the shear component of the stress as
$$\begin{eqnarray}\displaystyle {\it\eta}=\frac{\langle {\it\Sigma}^{12}\rangle }{\dot{{\it\gamma}}}={\it\eta}_{f}+\frac{\langle \unicode[STIX]{x1D634}^{12}\rangle }{\dot{{\it\gamma}}}. & & \displaystyle\end{eqnarray}$$
For what follows, we define the intrinsic viscosity for a particle dispersion at solid volume fraction
${\it\Phi}=4/3{\rm\pi}(a/L)^{3}N$
, as
$$\begin{eqnarray}\displaystyle [{\it\eta}]=\frac{{\it\eta}-{\it\eta}_{f}}{{\it\eta}_{f}{\it\Phi}}=\frac{\langle \unicode[STIX]{x1D634}^{12}\rangle }{{\it\eta}_{f}{\it\Phi}\dot{{\it\gamma}}}. & & \displaystyle\end{eqnarray}$$
To evaluate the heterogeneity of the stress in a realization, we define the deviatoric shear stress as
$$\begin{eqnarray}\displaystyle {\it\delta}{\bf\sigma}={\bf\sigma}-\langle {\bf\sigma}_{f}\rangle . & & \displaystyle\end{eqnarray}$$
Within the SPM formalism, the instantaneous volume-averaged stress of the flowing dispersion is given by Iwashita & Yamamoto (Reference Iwashita and Yamamoto2009) and Kobayashi & Yamamoto (Reference Kobayashi and Yamamoto2011)
$$\begin{eqnarray}\displaystyle {\it\bf\Sigma}=\frac{1}{V}\int \text{d}\boldsymbol{r}\,[{\bf\sigma}-\boldsymbol{r}{\it\rho}{\it\phi}\boldsymbol{f}_{p}+\boldsymbol{r}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}({\it\rho}\boldsymbol{u})], & & \displaystyle\end{eqnarray}$$
where the last term in square brackets comes from the convective momentum-flux tensor required to define the full stress tensor. The particle contribution to the stress is then
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D668}={\it\bf\Sigma}-\langle {\bf\sigma}_{f}\rangle . & & \displaystyle\end{eqnarray}$$
This instantaneous particle stress
$\unicode[STIX]{x1D668}$
can be used to define the corresponding values of the instantaneous viscosity and intrinsic viscosity, as in (2.11) and (2.12). To further evaluate the role of the inter-particle forces, and the accuracy of the lubrication interactions within the SPM, we divide the body force into hydrodynamic and direct particle–particle Lennard-Jones (LJ) force terms
$\boldsymbol{f}_{p}=\boldsymbol{f}_{p}^{H}+\boldsymbol{f}_{p}^{C}$
(see appendix C), such that
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\bf\Sigma}=\langle {\bf\sigma}_{f}\rangle +\unicode[STIX]{x1D668}_{H}+\unicode[STIX]{x1D668}_{C}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D668}_{H}=\frac{1}{V}\int \text{d}\boldsymbol{r}\,[{\it\delta}{\bf\sigma}-\boldsymbol{r}{\it\rho}{\it\phi}\boldsymbol{f}_{p}^{H}+\boldsymbol{r}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}({\it\rho}\boldsymbol{u})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D668}_{C}=-\frac{1}{V}\int \text{d}\boldsymbol{r}\,\boldsymbol{r}{\it\rho}{\it\phi}\boldsymbol{f}_{p}^{C}, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D668}=\unicode[STIX]{x1D668}_{H}+\unicode[STIX]{x1D668}_{C}$
. In particular, for spherical particles,
$\unicode[STIX]{x1D668}_{C}$
can be evaluated analytically and reduces to the virial expression for the stress (Evans & Morriss Reference Evans and Morriss2008)
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D668}_{C}=-\frac{1}{V}\mathop{\sum }_{i<j}\boldsymbol{R}_{ij}\boldsymbol{F}_{ij}^{C}, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{R}_{ij}=\boldsymbol{R}_{j}-\boldsymbol{R}_{i}$
is the distance vector from particle
$i$
to particle
$j$
, and
$\boldsymbol{F}_{ij}$
is the force on particle
$i$
due to particle
$j$
. We note that, in SD simulations, the sharp interface of the particle is taken into account, and no inter-particle interactions are required. Particle overlap is prevented by including the lubrication forces (Brady & Bossis Reference Brady and Bossis1988; Brady Reference Brady1993). In the current formulation of the SPM, we are not able to accurately recover the lubrication forces when the surface-to-surface distance between particles becomes smaller than the grid spacing. Therefore, we include a truncated LJ-type potential to prevent particle overlap. This inter-particle interaction, which depends only on the centre-to-centre distance, gives rise to a repulsive force when two particles are within contact distance
$r\simeq 2a$
.
3 Results
When solid or viscous objects are immersed in a host fluid under a linear shear flow, they induce an inhomogeneity in the flow field, and a decoupling of the rotation and strain flows takes place, due to the imposed boundary conditions. For the case of uniform shear, the velocity field can be decomposed into pure rotation and pure strain. Measured in the shear plane, with respect to the mean flow direction, the strain field is extended in the
${\rm\pi}/4$
angle and compressed in the
$3{\rm\pi}/4$
angle (Batchelor Reference Batchelor1967). Such anisotropic behaviour in the flow and strain fields is ubiquitous, for example in the flow behaviour of polymeric liquids (Milner Reference Milner1993) or supercooled liquids (Furukawa & Tanaka Reference Furukawa and Tanaka2009). When a particle of non-negligible size is immersed in the flow, there is a pure rotational flow around the particle due to the no-slip boundary condition, and the flow is modulated so that the strain tensor becomes anisotropic around the particle (note that this property is reproduced using MPC, a particle-based hydrodynamic simulation method (Ji et al.
Reference Ji, Jiang, Winkler and Gompper2011)).
We consider here the simulation results for four different cases of particles in simple-shear flow in the
$x$
–
$y$
plane: (1) a single particle, (2) a single rigid linear bead chain, (3) two colliding spherical particles, and (4) a dense dispersion of spherical particles. To test the validity of our method, we compare the results obtained in the first three cases with available analytical theories and SD calculations. We take as canonical simulation units the grid spacing
${\it\Delta}$
, the fluid density
${\it\rho}_{f}={\it\rho}$
and the viscosity
${\it\eta}_{f}$
. All remaining units can be obtained from these three; the units of mass, time and pressure are
${\it\rho}{\it\Delta}^{3}$
,
${\it\rho}{\it\Delta}^{2}/{\it\eta}_{f}$
and
${\it\eta}_{f}/{\it\rho}{\it\Delta}^{2}$
, respectively. All the simulations presented here were performed for periodic cubic systems of lateral size
$L=64$
(volume
$V=L^{3}$
), under steady shear. For the low-
$Re$
simulations, the shear rate
$\dot{{\it\gamma}}$
was chosen small enough (
$Re\simeq 10^{-3}-10^{-2}$
) that inertial effects are negligible and we can compare our results with those obtained at
$Re=0$
. We take the
$x$
–
$y$
plane to be the shear plane, with
$x$
the shear-flow direction,
$y$
the shear-gradient direction, and
$-z$
the vorticity axis. Unless stated otherwise, the thickness of the smooth interface is fixed at
${\it\zeta}/{\it\Delta}=2$
; the particle resolution was varied depending on the system under study, we have used
$a/{\it\Delta}=2,4,8$
. Of particular relevance to our study is the ratio between the particle radius and the system size
${\it\epsilon}=a/L$
. For large values of
${\it\epsilon}$
, finite size effects and the use of the periodic boundaries conditions become important. However, for
${\it\epsilon}\lesssim 0.1$
, these effects can be safely ignored, at least with regards to the stress and viscosity calculations (this is not true for the calculation of diffusion coefficients for example) (Iwashita & Yamamoto Reference Iwashita and Yamamoto2009).
3.1 Single particle under shear flow
For the case of a single spherical particle in shear flow, the equations governing the Stokes flow (
$Re=0$
) can be solved analytically (Lin, Peery & Schowalter Reference Lin, Peery and Schowalter1970; Mikulencak & Morris Reference Mikulencak and Morris2004). The expression for the peculiar velocity field is
$$\begin{eqnarray}\displaystyle {\bf\xi}=\dot{{\it\gamma}}\left(\begin{array}{@{}c@{}}\displaystyle -\frac{a^{5}y}{2r^{5}}-\frac{5a^{3}x^{2}y}{2r^{5}}\left(1-\frac{a^{2}}{r^{2}}\right)\\ \displaystyle -\frac{a^{5}x}{2r^{5}}-\frac{5a^{3}xy^{2}}{2r^{5}}\left(1-\frac{a^{2}}{r^{2}}\right)\\ \displaystyle -\frac{5a^{3}xyz}{2r^{5}}\left(1-\frac{a^{2}}{r^{2}}\right)\end{array}\right), & & \displaystyle\end{eqnarray}$$
where
$a$
is the particle radius. Note that in a region close to the surface of a large enough particle
$(r-a\ll a,~a\rightarrow \infty )$
, the flow field represents an extensional flow, because
${\bf\xi}=\boldsymbol{u}-\dot{{\it\gamma}}y\boldsymbol{e}_{x}\simeq -\dot{{\it\gamma}}a^{5}/2r^{4}(y/r,x/r,0).$
The resulting deviatoric Newtonian stress is given by
$$\begin{eqnarray}\displaystyle {\it\delta}{\it\sigma}^{12}=\dot{{\it\gamma}}{\it\eta}\left[-\frac{a^{3}}{2r^{3}}\left(5-\frac{8a^{2}}{r^{2}}\right)+\frac{5a^{3}x^{2}y^{2}}{r^{7}}\left(5-\frac{7a^{2}}{r^{2}}\right)\right], & & \displaystyle\end{eqnarray}$$
where
${\it\delta}{\it\sigma}^{12}$
is anisotropic with a periodicity of
${\rm\pi}/2$
, because
$4x^{2}y^{2}=(r\cos {\it\theta})^{4}\sin 2{\it\varphi}$
in spherical polar coordinates, where
${\it\theta}=\arccos (z/\sqrt{x^{2}+y^{2}+z^{2}})$
and
${\it\varphi}=\arctan (y/x)$
are the polar and azimuthal angles. On the
$z=0$
plane,
${\it\delta}{\it\sigma}^{12}$
exhibits a sign inversion slightly outside the spherical surface in the region
$\sqrt{19/15}a\leqslant r\leqslant \sqrt{8/5}a$
.
Figure 1. Comparison of the peculiar fluid velocity
$(1-{\it\phi}){\bf\xi}$
and normalized deviatoric Newtonian shear stress field
${\it\delta}{\it\sigma}^{12}/{\it\eta}_{f}\dot{{\it\gamma}}$
for a single spherical particle, under simple shear (
$\dot{{\it\gamma}}=10^{-3}$
), placed at the origin: (a,d) analytical solution (unbounded flow) and (b,e) numerical solution using the SPM (with periodic boundary conditions), (c,f) normalized shear stress difference
$|{\it\delta}{\it\sigma}^{12}-{\it\delta}{\it\sigma}^{12}|/{\it\eta}_{f}\dot{{\it\gamma}}$
. Solid lines show the particle surface assuming a sphere of radius
$a$
, dashed lines show the surface of the outer boundary
$a^{\prime }=a+{\it\zeta}/2$
(the radius of the particle including the interfacial boundary layer). (a–c) Shows results for
$a/{\it\Delta}=4$
(
${\it\epsilon}=0.0625$
), (d–f) for
$a/{\it\Delta}=8$
(
${\it\epsilon}=0.125$
).
To show that the SPM provides the same velocity profile and stress distribution, a simulation of one spherical particle in simple-shear flow is conducted (
$\dot{{\it\gamma}}=10^{-3}$
). The centre of mass of the spherical particle, of radius
$a/{\it\Delta}=4$
(
${\it\epsilon}=0.0625$
) is initially set at the centre of the box
$\boldsymbol{r}=\mathbf{0}$
, where the shear-flow velocity vanishes. The particle Reynolds number for this system is
$Re={\it\rho}\dot{{\it\gamma}}a^{2}/{\it\eta}_{f}=1.6\times 10^{-2}$
. As expected, the particle position remains constant during the simulation, with the particle rotating about its centre with a mean angular velocity of
$\dot{{\it\gamma}}/2$
. Figure 1(a–c) shows the comparison between the (a) analytical and (b) numerical solutions for the peculiar fluid velocity field
$(1-{\it\phi}){\bf\xi}$
and the deviatoric Newtonian stress field
${\it\delta}{\it\sigma}^{12}$
(with the density map representing the magnitude). Simulation results are presented as time averages in the steady state (over a time interval
$\dot{{\it\gamma}}t\sim 10^{2}$
). Both the velocity and stress fields exhibit very good agreement, which is also supported by the plot in figure 1(c), which shows the distribution of the difference in the stress
${\it\delta}{\it\sigma}^{12}$
between the analytical and numerical calculations. Since, in our numerical simulations, the particle/fluid interface is smeared out, the fluid can enter the particle interface region (
$a-{\it\zeta}/2<r<a$
), giving rise to relatively high stress differences. To check the effect of the boundaries, we have also performed simulations for a larger value of
${\it\epsilon}=0.125$
(figure 1
d–f). We again obtain good qualitative agreement with the analytical results, although notable differences in the velocity at the system boundaries start to appear, accompanied by an increase in the stress differences. This is expected, since the analytical solution we are comparing against is for a spherical particle in an unbounded flow, and our simulations assume periodic boundary conditions.
For this system, we expect that Einstein’s law for the viscosity
${\it\eta}={\it\eta}_{f}(1+5/2{\it\phi})$
should be applicable, since the dispersion is dilute, with a volume fraction of
${\it\Phi}=4/3{\rm\pi}(a/L)^{3}=1.0\times 10^{-3}$
(
${\it\epsilon}=0.0625$
), such that
${\it\eta}=1.00255$
. Indeed, from the numerical calculations using the SPM, the shear viscosity is calculated to be 1.00244. For the system at
${\it\epsilon}=1.020$
, we obtain
${\it\eta}=1.019$
, compared with the theoretical value of
${\it\eta}=1.020$
. This agreement between the simulation and Einstein’s law was already reported by Kobayashi & Yamamoto (Reference Kobayashi and Yamamoto2011). In addition, we have also checked the dependence of the viscosity on
$Re$
, by varying the shear rate
$10^{-4}\leqslant \dot{{\it\gamma}}\leqslant 1.6$
(
$a/{\it\Delta}=4,{\it\epsilon}=0.0625$
). Our results are given in figure 2, together with previous simulation data (Mikulencak & Morris Reference Mikulencak and Morris2004) obtained using a finite element method for
${\it\epsilon}=0.01,0.1$
. Excellent agreement is obtained up to
$Re\simeq 10^{1}$
, however, differences start to appear as
$Re$
is increased further. At high
$Re$
, we should decrease the time step, to properly account for the advection, as well as increase the resolution, to resolve the small-scale variations in the flow. The latter is straightforward, but the former requires special care due to the non-monotonic dependence of the error in the SPM method with respect to the time step. Such considerations are outside the scope of this work.
Figure 2. (a) Angular velocity
${\it\omega}_{z}$
and (b) intrinsic viscosity
$[{\it\eta}]$
as a function of
$Re$
for a single spherical particle under simple-shear flow, for
${\it\epsilon}=0.0625$
(filled circles). Results obtained by Mikulencak & Morris for two different system sizes
${\it\epsilon}=0.01$
(open squares) and
${\it\epsilon}=0.1$
(open circles), are also shown.
3.2 Single bead chain
In practice, most industrial and technological applications which involve the motion of particles under an imposed flow deal with non-spherical particles. Even for the simple case of a single particle in shear flow, a wide variety of dynamic modes can be found for high
$Re$
. Recent work by Huang et al. (Reference Huang, Yang, Krafczyk and Lu2012) for spheroidal particles in Couette flow (
$Re\leqslant 700$
), using a multi-relaxation-time lattice Boltzmann method, has shown the existence of distinct tumbling and log rolling modes, which are sensitive to both the initial orientation and the precise value of
$Re$
. Here, in an attempt to validate our model, we focus on the case of an isolated axisymmetric particle under simple shear in the low-
$Re$
regime, for which analytic solutions are known (
$Re=0$
). These were first obtained by Jeffery for ellipsoidal particles (Jeffery Reference Jeffery1922), but they are applicable to any rigid axisymmetric body once one defines its equivalent aspect ratio
$r_{e}$
(which need not coincide with the geometric aspect ratio) (Bretherton Reference Bretherton1962; Cox Reference Cox2006; Zhang et al.
Reference Zhang, Jack, Smith and Montgomery-Smith2011). The orientation of the body is parametrized using the three Euler angles
${\it\theta},{\it\phi},{\it\psi}$
(see figure 3). Here,
${\it\theta}$
is the polar angle between the vorticity axis
$z$
and the symmetry axis of the body,
${\it\phi}$
is the angle of rotation in the shear plane, and
${\it\psi}$
the angle of rotation around the symmetry axis of the body (not shown). Without loss of generality, we can consider the fluid velocity at the centre of mass to be zero, hence the particle will experience no net translation (this is accomplished by placing the centre of mass of the chain at the centre of the simulation box). The particle will therefore experience a pure rotational motion, completely specified by the three Euler angles. As discovered by Jeffery, the particles will follow closed periodic trajectories, specified by an orbit constant
$C$
, which depends only on the initial orientation. These are the so-called Jeffery’s orbits and they take the following form
Figure 3. Schematic representation of a (rigid) linear chain with
$n=5$
beads in simple-shear flow.
$$\begin{eqnarray}\displaystyle \tan {\it\phi}\! & = & \displaystyle \!r\tan \left(\frac{\dot{{\it\gamma}}t}{r+1/r}\right)+\tan {\it\phi}_{0},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \tan {\it\theta}\! & = & \displaystyle \!\left(\frac{Cr}{r^{2}\cos ^{2}{\it\phi}+\sin ^{2}{\it\phi}}\right)^{1/2},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\it\psi}\! & = & \displaystyle \!\int _{0}^{t}\text{d}t(\dot{{\it\gamma}}/2-\dot{{\it\phi}})\cos {\it\theta},\end{eqnarray}$$
$C$
is given by 
$$\begin{eqnarray}\displaystyle C=\frac{\tan {\it\theta}}{r}(r^{2}\cos ^{2}{\it\phi}+\sin ^{2}{\it\phi})^{1/2}. & & \displaystyle\end{eqnarray}$$
We recall that
$r$
in (3.3) and (3.4) should be replaced with the equivalent aspect ratio
$r_{e}$
when dealing with general axisymmetric bodies (instead of the ellipsoids originally considered by Jeffery). This equivalent aspect ratio can be obtained if one knows the period of oscillation, since
$$\begin{eqnarray}\displaystyle T=\frac{2{\rm\pi}}{\dot{{\it\gamma}}}(r+r^{-1}). & & \displaystyle\end{eqnarray}$$
An alternative, proposed by Cox (Reference Cox2006), is to measure the ratio of the torques on the body in a fixed horizontal and vertical orientation, as
$$\begin{eqnarray}\displaystyle r_{e}=\sqrt{\frac{[{\it\omega}]_{V}}{[{\it\omega}]_{H}}}=\sqrt{\frac{[{\it\tau}]_{V}}{[{\it\tau}]_{H}}}, & & \displaystyle\end{eqnarray}$$
where
${\it\tau}_{V}$
(
${\it\tau}_{H}$
) is the hydrodynamic torque exerted on a particle fixed in the vertical
${\it\phi}=0$
(horizontal
${\it\phi}={\rm\pi}/2$
) position within the shear plane
${\it\theta}={\rm\pi}/2$
. However, we have obtained more accurate results by using the relation between the oscillation period
$T$
and the aspect ratio
$r$
: we use the oscillation period obtained from the simulation results of a free particle and invert (3.5) to obtain the equivalent aspect ratio. Finally, the viscosity of a dilute suspension of such axisymmetric rigid bodies (in which the particle–particle interactions are ignored) is also known (Batchelor Reference Batchelor1970; Hinch & Leal Reference Hinch and Leal1972), and in the absence of Brownian diffusion, takes the form
$$\begin{eqnarray}\displaystyle [{\it\eta}]=A\sin ^{4}{\it\theta}\sin ^{2}2{\it\phi}+2B\sin ^{2}{\it\theta}+\frac{2}{I_{\Vert }}, & & \displaystyle\end{eqnarray}$$
where the intrinsic viscosity
$[{\it\eta}]$
is defined such that
${\it\eta}={\it\eta}_{0}(1+[{\it\eta}]{\it\Phi})$
, with
${\it\Phi}$
the particle number density, and
$A$
and
$B$
given by
$$\begin{eqnarray}\displaystyle & \displaystyle A=\frac{J_{\Vert }}{I_{\Vert }J_{\bot }}+\frac{1}{I_{\Vert }}-\frac{2}{I_{\bot }}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle B=\frac{1}{I_{\bot }}-\frac{1}{I_{\Vert }}, & \displaystyle\end{eqnarray}$$
where the
$I$
and
$J$
are form factors which depend only on the shape (aspect ratio) of the particle; assuming
$r>1$
they are given by
$$\begin{eqnarray}\displaystyle & \displaystyle I_{\Vert }=2r\int _{0}^{\infty }\frac{\text{d}x\,}{(r^{2}+x)^{1/2}(1+x)^{3}}=(3s+2r^{2}-5)r^{2}\frac{q}{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle J_{\Vert }=r\int _{0}^{\infty }\frac{\text{d}x\,x}{(r^{2}+x)^{1/2}(1+x)^{3}}=((1-4r^{2})s+2r^{2}+1)r^{2}\frac{q}{4}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle I_{\bot }=r(r^{2}+1)\int _{0}^{\infty }\frac{\text{d}x\,}{(r^{2}+x)^{3/2}(1+x)^{2}}=(r^{2}(1-3s)+2)(r^{2}+1)q, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle J_{\bot }=r\int _{0}^{\infty }\frac{\text{d}x\,x}{(r^{2}+x)^{3/2}(1+x)^{2}}=((2r^{2}+1)s-3)r^{2}q, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle s=\frac{a\text{cosh}\,r}{r(r^{2}-1)^{1/2}}\quad q=\frac{1}{(r^{2}-1)^{2}}. & \displaystyle\end{eqnarray}$$
$n=3$
spherical beads, particle radius
$a/{\it\Delta}=2$
, at the centre the cubic simulation box of length
$L=64$
(
${\it\epsilon}=an/L\simeq 9\times 10^{-2}$
), under simple-shear flow with shear rate
$\dot{{\it\gamma}}=5\times 10^{-4}$
(
$Re={\it\rho}\dot{{\it\gamma}}(an)^{2}/{\it\eta}_{f}=1.8\times 10^{-2}$
). We consider three different initial orientations,
${\it\theta}={\rm\pi}/10,{\rm\pi}/4,{\rm\pi}/2$
with
${\it\phi}=0$
, which corresponds to orbit constants of
$C=0.32,1,\infty$
. As explained above, the effective aspect ratio of this rigid body
$r_{e}=2.79$
can be obtained after measuring the period of oscillation. This value of
$r_{e}$
can then be used to obtain the corresponding Jeffery orbits, which are compared to our simulation results in figure 4. Excellent agreement is obtained for the orbit of the spherical chain. Results for larger system size
$L$
, smaller
${\it\epsilon}$
, exhibit no system size dependence (not shown). However, when comparing the results for the intrinsic viscosity, we observe a clear difference between both sets of data. To validate our simulation procedure, we have also computed the viscosity given by SD, which provides the exact solution (at
$Re=0$
). The SD results are obtained by taking the SPM simulation data for the trajectories and velocities and using them to solve the corresponding resistance problem. Good agreement with our data, particularly for chain configurations parallel or perpendicular to the shear flow (
${\it\phi}=0,{\rm\pi}/2$
), is obtained. For configurations with
${\it\phi}\sim {\rm\pi}/4$
, for high values of the orbit constant, there is a clear overestimation of the intrinsic viscosity in the SPM measurements. We believe this discrepancy, which can be of the order of
$\simeq 10\,\%$
, is due to the diffuse interface used to represent the fluid–particle boundary. Specifically, due to the strong fluid flow along this diffuse particle surface. The reason this effect is highly orientation-dependent is simple: for both
${\it\phi}=0$
and
${\it\phi}={\rm\pi}/2$
the fluid velocity parallel to the symmetry axis of the body is small, and does not vary strongly along the surface, in contrast to the case for
${\it\phi}={\rm\pi}/4$
(where the velocity along the body of the particle increases with height). The effective size of the particles can be considered to increase slightly depending on the relative orientation with respect to the shear flow. To test this, we have validated against SD calculations, using the same effective aspect ratio as before, but varying only the radius of the spherical beads. The results for the particle shear stress
$\unicode[STIX]{x1D634}^{12}$
, summarized in figure 5, indeed show a strong dependence on the hard-sphere radius of the particles. A difference of only
$5\,\%$
, which is well within the relative size of our diffuse interface, can account for the variations between our data and the reference SD data. Finally, we see that even though an effective aspect ratio can be used to transfer Jeffery’s solution to arbitrary axisymmetric rigid bodies, the stress on the fluid (not surprisingly) depends strongly on the shape of the particle, and it cannot be accurately represented using the same aspect ratio (as shown by the poor comparison between Jeffery’s solution and both the SD and SPM simulations).
Figure 4. Jeffery’s orbits for a single spherical bead chain, of geometric aspect ratio
$r=3$
and effective aspect ratio
$r_{e}\simeq 2.8$
. The Euler angles (
${\it\psi},{\it\phi},{\it\theta}$
) characterizing the orientation of the chain (see figure 3), are given in (a), (b) and (c), respectively, for one full period. The corresponding intrinsic viscosity is given in plot (d). Simulation data using our SPM model is represented using the filled symbols and Jeffery’s solution (using the effective aspect ratio
$r_{e}$
) is given by the dashed lines. For the intrinsic viscosity
$[{\it\eta}]$
, data obtained from SD calculations are also given (solid lines).
Figure 5. Particle shear stress
$\unicode[STIX]{x1D634}^{12}$
for a single rigid bead chain of aspect ratio
$r=3$
(
$r_{e}\simeq 2.8$
), for different orbit parameters: (a)
$C=\infty$
, (b)
$C=1.0$
, (c)
$C=0.32$
. Simulation data obtained using our SPM model, with a diffuse fluid/particle interface, is given by the solid symbols. Results from SD calculations using the same hard-sphere radius
$a$
as the SPM simulations is given by the solid line, results obtained using a modified hard-sphere particle radius
$a^{\prime }/a=0.95,1.05$
, are given as dashed lines.
3.3 Collision of a pair of particles
Interactions between colloidal particles strongly affect the rheology and flocculation behaviour under shear flow. Due to the hydrodynamic interactions from the host fluid, a pair of immersed particles, which are initially separated, will never come into contact, thanks to the lubrication forces between them. Accurately representing these interactions, which diverge as the surface-to-surface separation between the particles goes to zero, poses significant problems when dealing with many-particle systems. However, simulation frameworks for dealing with hydrodynamically interacting particles have already been developed (Durlofsky, Brady & Bossis Reference Durlofsky, Brady and Bossis1987; Brady & Bossis Reference Brady and Bossis1988; Ladd Reference Ladd1988), which take these lubrication effects into account, at least at the level of SD. Numerical studies based on this formulation have shown that hydrodynamic effects can account for the concentration dependence of the viscosity of semi-dilute suspensions (Batchelor & Green Reference Batchelor and Green1972; Chwang & Wu Reference Chwang and Wu1974). There have been many discussions regarding the rheology of more highly concentrated suspensions, where the viscosity dependence on the concentration
${\it\phi}$
is compared between various simulations (Ladd Reference Ladd1990; Boek et al.
Reference Boek, Coveney, Lekkerkerker and van der Schoot1997; Sierou & Brady Reference Sierou and Brady2001) and experiments (van der Werff & de Kruif Reference van der Werff and de Kruif1989; Shikata & Pearson Reference Shikata and Pearson1994; Singh & Nott Reference Singh and Nott2003), and the inter-particle forces have been shown to have a dramatic effect on the viscosity.
The simplest example for the rheology of interacting colloidal particles is that of two colliding spheres under shear flow, which we turn to now. This problem was first considered by Batchelor and Green some 40 years ago, when they analytically solved for the interaction between a pair of colloids in an infinite fluid, as a function of the relative particle separation distance
$\overline{\boldsymbol{r}}=\boldsymbol{r}_{2}-\boldsymbol{r}_{1}$
, under the Stokes approximation (Batchelor & Green Reference Batchelor and Green1972; Russel, Saville & Schowalter Reference Russel, Saville and Schowalter1992). In their theory, the interaction force dipoles between a pair of particles is solved, and gives rise to the following particle contribution to the shear stress
$$\begin{eqnarray}\displaystyle V\unicode[STIX]{x1D634}^{12}=\frac{20}{3}{\rm\pi}a^{3}{\it\eta}_{f}\dot{{\it\gamma}}\left[(1+K(\overline{r}))+\displaystyle \frac{\overline{x}^{2}+\overline{y}^{2}}{\overline{r}^{2}}L(\overline{r})+\frac{2\overline{x}^{2}\overline{y}^{2}}{\overline{r}^{4}}M(\overline{r})\right], & & \displaystyle\end{eqnarray}$$
where
$\overline{\boldsymbol{r}}=(\overline{x},\overline{y},\overline{z})$
and the coefficients
$K,L,M$
are defined as
$$\begin{eqnarray}\displaystyle K(\overline{r})\! & = & \displaystyle \!-2\tilde{r}^{-5}+o(\tilde{r}^{-6}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle L(\overline{r})\! & = & \displaystyle \!-{\textstyle \frac{5}{2}}\tilde{r}^{-3}+10\tilde{r}^{-5}+{\textstyle \frac{25}{4}}\tilde{r}^{-6}+o(\tilde{r}^{-6}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle M(\overline{r})\! & = & \displaystyle \!{\textstyle \frac{25}{2}}\tilde{r}^{-3}-35\tilde{r}^{-5}+25\tilde{r}^{-6}+o(\tilde{r}^{-6}).\end{eqnarray}$$
$a$
) and
$\tilde{r}=(\overline{r}/a)$
is the relative distance normalized by this radius
$a$
.
Figure 6. A pair of particles in a host fluid under shear flow is shown for three consecutive times
$\dot{{\it\gamma}}t=-2.53,~-1.17,$
and
$0.87$
, for a shear rate
$\dot{{\it\gamma}}=10^{-3}$
. The peculiar velocity field
${\bf\xi}$
and deviatoric Newtonian shear stress
${\it\delta}{\it\sigma}^{12}/{\it\eta}_{f}\dot{{\it\gamma}}$
are displayed using the arrows and colour map, respectively.
We performed a simulation for a pair of spherical particles, radius
$a/{\it\Delta}=8$
(
${\it\zeta}/{\it\Delta}=1$
), to confirm whether our rheological calculation can accurately reproduce the behaviour described above, for two approaching particles. As before, we impose a constant shear flow (
$\dot{{\it\gamma}}=10^{-3}$
). The two particles are initially placed at positions
$\boldsymbol{r}_{\pm }=(\pm 20{\it\Delta},\pm {\rm\Delta}y,0)$
, such that they are displaced from the origin by a distance of
${\rm\Delta}y=a$
in the
$y$
-direction (shear-gradient direction). The imposed shear flow will pull the particles together, ensuring that they will collide. Therefore,
${\rm\Delta}y$
represents the impact parameter of this particle collision. We define the time origin (
$t=0$
) when the particles have the same position along the
$x$
-axis (shear-flow direction):
$x_{+}(t=0)=x_{-}(t=0)$
. Thus, negative times correspond to instants before the collision, when particles are being squeezed together, positive times to instants after the collision, when the flow is pulling the particles apart.
In figure 6, a cross-section view of the spatial distribution of the deviatoric Newtonian shear stress
${\it\delta}{\it\sigma}^{12}$
is plotted as a colour map, together with a
$2$
D projection of the peculiar velocity field
${\bf\xi}$
(represented by the arrows) at
$z=0$
. As shown in figure 6(a),
${\it\delta}{\it\sigma}^{12}$
and the velocity
${\bf\xi}$
have anisotropic structures with
${\rm\pi}/2$
rotational periodicity, as in the one-particle case (see figure 1). When the particles are in near contact with each other (b) (at
$\dot{{\it\gamma}}t=-1.17$
), the direct inter-particle LJ interaction will give a significant contribution to the stress. The surrounding fluid then starts flowing into the interstitial space between the particles, to allow subsequent particle separation. Because the particles are nearly at contact and are pushing out the interstitial fluid between them,
${\it\delta}{\it\sigma}^{12}$
is largely negative, as observed around the contact point. At
$\dot{{\it\gamma}}t=0.87$
, the particles are being pulled apart, with the fluid coming into the region between them,
${\it\delta}{\it\sigma}^{12}$
is then strongly positive, as shown in (c). Around the particle surface in particular, the deviatoric Newtonian stress
${\it\delta}{\it\sigma}^{12}$
represents the stress due to the particle pushing or dragging the fluid.
The shear stress
$\unicode[STIX]{x1D634}^{12}$
, can be decomposed into the contributions given in (2.9). In particular, within the SPM framework, the effect of the hydrodynamic interactions is represented by
$\unicode[STIX]{x1D668}_{H}$
, and it is of interest to see how large it is compared to the pure inter-particle interactions (
$\unicode[STIX]{x1D668}_{C}$
). To validate the accuracy of the stress calculations within the SPM, we compare against Batchelor’s analytical solution and SD results. For this, we use the trajectory of the particles
$\boldsymbol{r}_{1}(t)$
and
$\boldsymbol{r}_{2}(t)$
obtained from our SPM simulations. Using (3.11), we obtain an estimate for the shear stress
$\unicode[STIX]{x1D634}^{12}$
within Batchelor’s force dipole approximation. For the comparisons with SD, the particle trajectories are used to solve the corresponding mobility problem, both with and without lubrication corrections. For both cases (Batchelor and SD), we must define the hard-sphere radius of the particles. However, given the diffuse interface used in our simulations, it is not clear what value we should use. For this, we have fitted an effective radius
$a_{eff}$
, in such a way that the stress at large separations
$\dot{{\it\gamma}}t\simeq 2$
obtained from our SPM simulations matches the stresses obtained from Batchelor’s solution or SD. We obtained an effective radius
$a_{eff}/a\simeq 0.975$
(
$a_{eff}/{\it\Delta}=7.8$
), which is consistent with the SPM parameters, since the interface region is smeared out over distances
$7.5\leqslant r/{\it\Delta}\leqslant 8.5$
(
${\it\xi}/{\it\Delta}=1.0$
). Furthermore, in this large-separation limit, where lubrication forces are negligible, both Batchelor’s solution and SD predict the same values for the stress (using equivalent particle radii). This shows that any effects due to the use of periodic boundary conditions can be ignored. However, for the results obtained during the collision event, when
$\overline{r}\simeq 2a$
and the interface of both particles can overlap, strong differences are seen between Batchelor’s results and the SD calculations. In addition, there is also a strong dependence on the precise value used for the effective radius
$a_{eff}$
, which makes a direct comparison with our SPM results difficult. These differences are considerably reduced at large particle separations, and can be made arbitrarily small by increasing the particle resolution (as discussed below). Finally, we note that when we compare our results with SD, we do not included any additional inter-particle repulsive forces (
$\unicode[STIX]{x1D668}_{C}^{SD}=0$
); we thus compare the stress obtained from the SPM calculations (with the short-range LJ-type repulsion) with the purely hydrodynamic SD stress. Figure 7 shows the time evolution of the total deviatoric shear stress
$\unicode[STIX]{x1D634}^{12}$
obtained from our simulations, as well as the pure hydrodynamic contribution to the shear stress
$\unicode[STIX]{x1D634}_{H}^{12}$
, together with the relative positions
$\overline{x}$
and
$\overline{y}$
, and the surface-to-surface distance between the particles (
$\overline{r}-a$
). Our SPM simulation results are compared to Batchelor’s theory (3.11) and to SD calculations (with and without the lubrication corrections), using two different particles sizes: (1) an equivalent hard-sphere radius
$a/{\it\Delta}=8.0$
, and (2) the effective hard-sphere radius
$a_{eff}/{\it\Delta}=7.8$
.
Figure 7. Collision process for a pair of particles in linear shear flow. (a) The total deviatoric shear stress
$\unicode[STIX]{x1D634}^{12}$
obtained from our simulations is shown with filled blue points, and the deviatoric shear stress
$\unicode[STIX]{x1D634}_{H}^{12}$
is plotted with blue open circles, for particles of radius
$a/{\it\Delta}=8$
(
${\it\epsilon}=0.125$
). SPM results for the shear stress are compared with Batchelor’s approximation (dashed–dotted lines), and SD calculations with (solid lines) and without (dashed lines) lubrication corrections. The pink and blue lines correspond to results obtained for an equivalent hard-sphere radius
$a$
and an effective hard-sphere radius
$a_{eff}/a=0.975$
, respectively. (b) The surface-to-surface distance between the particles
$\overline{r}-2a$
, as well as the
$x$
(
$y$
) components of their relative position
$\overline{x}$
(
$\overline{y}$
).
Figure 8. Relative deviation of the particle contribution to the shear stress
$\unicode[STIX]{x1D634}^{12}$
obtained by SPM with respect to Batchelor’s results is plotted as a function of particle resolution
$a/{\it\Delta}$
. This deviation is estimated in the time region
$\dot{{\it\gamma}}t\simeq 2$
where
$\unicode[STIX]{x1D668}^{SPM}=\unicode[STIX]{x1D668}_{H}^{DNS}$
. The impact parameter is set as
${\rm\Delta}y=a$
.
During the collision process, as seen in figure 7, when the particles come close to each other (
$\overline{r}\gtrsim 2$
;
$\dot{{\it\gamma}}t\sim -2$
), the deviatoric shear stress
$\unicode[STIX]{x1D634}^{12}$
obtained from the simulations
$\unicode[STIX]{x1D668}^{SPM}$
begins to increase due to the hydrodynamic interactions. When this distance reaches a value close to the particle diameter, the particles are subjected to an inter-particle repulsive interaction, which corresponds to a gap between
$\unicode[STIX]{x1D668}^{SPM}$
and
$\unicode[STIX]{x1D668}_{H}^{SPM}$
and partially represents the lubrication interaction in real dispersions (as evidenced by the comparison with the SD calculations that include the lubrication corrections
$\unicode[STIX]{x1D668}_{lub}^{SD}$
). Subsequently, the particles are pushed away from the centre due to the normal stress (
$\text{d}\overline{y}/\text{d}t>0$
), and a peak in
$\unicode[STIX]{x1D634}^{12}$
is soon achieved (the position of the peak is the same regardless of the calculation method). At
$t=0~(\overline{x}=0)$
, just after the direct inter-particle LJ interaction vanishes, hydrodynamic drag forces start to pull the particles closer to the centre line (
$\text{d}\overline{y}/\text{d}t<0$
). Note that the hydrodynamic stress
$\unicode[STIX]{x1D668}_{H}^{SPM}$
significantly underestimates the full stress during the particle collision (
$-1.3\lesssim \dot{{\it\gamma}}t<0$
), when the direct inter-particle LJ interactions are significant. Due to the nature of our method, it is clear that we could not expect to accurately take into account the lubrication effects when the surface-to-surface distance is less than the grid spacing. Including the direct inter-particle contribution to the stress gives much better agreement with SD, particularly when we compare against the results obtained using the effective radius
$a_{eff}$
(fitted to the stress outside of the collision region). However, there is still a relatively large error for the peak stress at
$\dot{{\it\gamma}}t\simeq -1$
(
${\approx}20{-}30\,\%$
). If required, one can consider more accurate representations for the lubrication forces, for example, by using a steeper inter-particle repulsion, or even by adding the lubrication corrections directly (Nguyen & Ladd Reference Nguyen and Ladd2002).
Finally, we show the dependence of the stress calculations on the particle resolution
$a/{\it\Delta}$
. Outside of the collision region (
$r>2a$
), we have seen that SD and Batchelor’s results coincide. However, full agreement with our SPM calculations was obtained only if we used an effective hard-sphere radius
$a_{eff}$
for both SD and Batchelor’s calculations. Using the same hard-sphere radius
$a$
, the SPM results slightly underestimate the stresses. This is reasonable given the diffuse nature of the particle interface within the SPM formulation. As further evidence for this, we show that by increasing the particle resolution
$a/{\it\Delta}$
, the difference between the stresses given by the SPM, and those obtained from Batchelor’s solution using the same hard-sphere radius
$a$
, is reduced (i.e.
$a/a_{eff}\rightarrow 1$
). In figure 8, the relative difference between these two stresses
${\rm\Delta}\unicode[STIX]{x1D634}^{12}=(\unicode[STIX]{x1D634}^{12,DNS}-\unicode[STIX]{x1D634}^{12,Batchelor})/\unicode[STIX]{x1D634}^{12,Batchelor}$
is plotted as a function of particle resolution
$a/{\it\Delta}$
, showing that the deviation becomes smaller at large enough particle radius. Thus, it is clear that the deviation originates from the use of the diffuse fluid–particle interface. It is expected to vanish if we employ a large enough radius.
3.4 Viscosity of concentrated colloidal dispersions
Having established the accuracy with which the SPM method is able to resolve the hydrodynamics of particles in shear flow, we now consider the viscosity of dense colloidal suspensions. It is well known that in the presence of shear, colloidal dispersions will exhibit several interesting non-Newtonian behaviours, such as shear thinning or shear thickening. While a complete understanding of such behaviour is still lacking, computer simulations have significantly helped to establish the mechanisms behind these phenomena. Of particular importance is the ability to establish a link between the microstructure of the dispersion and the macroscopic rheological properties. The present method is well suited to study such problems. Indeed, previous work by Iwashita & Yamamoto (Reference Iwashita and Yamamoto2009) using the SPM with a zig-zag velocity profile has successfully studied the shear-thinning behaviour of Brownian colloidal particles over a wide range of Reynolds (
$10^{-4}<Re<10^{0}$
) and Péclet numbers (
$10^{-2}<Pe<10^{1}$
), where the latter is defined as
$$\begin{eqnarray}\displaystyle Pe=6{\rm\pi}{\it\eta}_{f}a^{3}\dot{{\it\gamma}}/k_{\mathit{B}}T. & & \displaystyle\end{eqnarray}$$
Figure 9. The viscosity of a dispersion of spherical colloidal particles as a function of Péclet number
$Pe$
, for various particle concentrations
${\it\Phi}$
. The data obtained using the proposed Lees–Edwards SPM method (filled symbols), is compared with the data previously obtained data using an externally driven zig-zag shear flow (open symbol) by Iwashita & Yamamoto (Reference Iwashita and Yamamoto2009). The dashed lines show the fit to the empirical formula given in (3.14). Error bars, representing the standard deviation of the mean, have been drawn for a selection of the simulation results.
We have compared our results using the Lees–Edwards SPM, to those obtained using the zig-zag flow (Iwashita & Yamamoto Reference Iwashita and Yamamoto2009). To facilitate the comparison we have used the same parameters. The size of the system is
$L/{\it\Delta}=64$
, with spherical particles of radius
$a/{\it\Delta}=4$
, interface thickness
${\it\zeta}/{\it\Delta}=2$
(
${\it\epsilon}=0.0625$
). The temperature is
$k_{\mathit{B}}T=7$
, as determined by equilibrium calculations in the absence of shear. Simulations were performed for volume fractions of
${\it\Phi}=0.31,0.41,0.46,0.51$
and 0.56, corresponding to particle numbers of
$N=300$
,
$400$
,
$450$
,
$500$
and
$550$
, respectively, for shear rates in the range
$10^{-4}\leqslant \dot{{\it\gamma}}\leqslant 10^{-1}$
. The total time
$t_{max}$
, corresponding to
$2\times 10^{6}$
steps, was the same for all systems. Thus, the total strain
${\it\gamma}=\dot{{\it\gamma}}t_{max}$
varied from
$10^{1}$
for the lowest shear rate, to
$10^{4}$
for the largest shear rate. Figure 9 summarizes our results for the dispersion viscosity
${\it\eta}$
as a function of
$Pe$
, together with the previously reported zig-zag results, and the corresponding fits to the following empirical formula (Iwashita & Yamamoto Reference Iwashita and Yamamoto2009)
$$\begin{eqnarray}\displaystyle {\it\eta}(Pe,{\it\Phi})={\it\eta}_{\infty }+\left(\frac{{\it\eta}_{0}-{\it\eta}_{\infty }}{1+b^{-1}({\it\Phi})Pe}\right), & & \displaystyle\end{eqnarray}$$
where
${\it\eta}_{0}({\it\eta}_{\infty })$
is the low (high) shear rate limiting viscosity, and
$b({\it\Phi})$
is a fitting parameter. As expected, upon increasing the
$Pe$
, and thus the relative importance of the hydrodynamic interactions over the thermal fluctuations, the shear-thinning behaviour is recovered for particle concentrations
${\it\Phi}\gtrsim 0.46$
. Excellent agreement is obtained between both methods, at least for
$Pe\gtrsim 10^{-1}$
and
${\it\Phi}\lesssim 0.51$
, which indicates that the use of the zig-zag profile is not a significant problem for rheological measurements of moderately dense spherical particle dispersions. At smaller
$Pe$
the results show clear differences for
${\it\Phi}=0.56$
, which indicates that the presence of the kinks in the zig-zag flow profile becomes important, and their effects on the structure and dynamics of the dispersion needs to be considered.
4 Conclusion
We have successfully extended the SPM formalism to compute the rheological properties of rigid bodies under simple shear. We have also shown that the methodology is applicable to bodies of arbitrary shape, not just to spherical particles. To simulate bulk systems, we use the appropriate Lees–Edwards periodic boundary conditions, which in turn require that the fluid equations of motion be expressed in an oblique time-dependent reference frame. We have given a general derivation for the fluid equations of motion, which show how the current approach may be extended to more complex flows (e.g. oscillatory shear). Our method provides explicit expressions for the local viscous shear stress, as well as the particle contribution to the stress. We have validated our results by analysing three well-known flow problems for particles in simple-shear flow: (1) a single particle, (2) a single rigid chain, and (3) two colliding particles. In addition, we have examined the shear-thinning behaviour of concentrated colloidal dispersions. In the three simple test cases, the simulation results for the fluid velocity and shear stress are quantitatively compared with the analytical solutions and the SD results. We obtain very good agreement for all the cases considered. However, we note that certain differences exist due to the diffuse nature of the particle/fluid interface. Such differences are localized around the particle, and will increase as the flow around the surface becomes more important. For the case of a single particle, the local stress will deviate from the analytical solution, as the particle boundary is not uniquely defined in our simulations. Furthermore, this deviation can depend on the relative orientation of the particle (as seen when calculating the stress of the single sphere or the rigid chain in the Jeffery Orbits), and is found to be larger along the compression axis. As particles will tend to align in this direction (Foss & Brady Reference Foss and Brady2000; Kumar Reference Kumar2010), the extent to which this discrepancy in the stress will affect the structure or dynamics in concentrated dispersions is a point which still needs to be considered. In addition, when two (or more) particles approach each other, we are not able to accurately reproduce the lubrication effects if the surface-to-surface distance is less than the grid spacing. If necessary, one could also add the lubrication corrections directly (Nguyen & Ladd Reference Nguyen and Ladd2002), although this becomes quite involved when non-spherical particles are considered. Overall, we consider such errors, which are typically
${\sim}10\,\%$
, to be acceptable, considering the generality of the method and the relative ease with which it can be implemented. Our method has been developed to provide a general framework which is applicable to complex host solvent from low to moderate Reynolds numbers, with arbitrary particle geometries, and is consistent with SD and its variants.
Acknowledgements
This work was supported by a Grant-in-Aid for Scientific Research on Innovative Areas ‘Synergy of Fluctuation and Structure: Foundation of Universal Laws in Nonequilibrium Systems’ (grant nos 25103010 and 26103522) from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan and a KAKENHI grant (no. 26247069) from the Japan Society for the Promotion of Science (JSPS). The authors would also like to thank Dr R. Seto for fruitful discussions regarding the Stokesian dynamics calculations, which have been performed with the open source RYUON simulator developed by K. Ichiki. Some of the numerical calculations have been performed on an SGI Altix ICE 8400EX at the ISSP of the University of Tokyo.
Appendix A. Navier–Stokes equation in a moving coordinate system
The (oblique) basis vectors are defined such that they move with the average flow velocity, and satisfy the periodic boundary conditions at all times. Thus, under shear flow with mean velocity
$\boldsymbol{U}(t)=\dot{{\it\gamma}}(t)y\boldsymbol{e}_{x}$
(
$\dot{{\it\gamma}}(t)=\text{d}{\it\gamma}(t)\,/\text{d}t\,$
the imposed shear rate at time
$t$
,
${\it\gamma}(t)$
the shear strain, and
$\boldsymbol{e}_{x}$
the Cartesian
$x$
-axis), the velocity of the lattice points (coordinate-flow velocity) must be equal to the shear-flow velocity
$\boldsymbol{U}$
. Since the basis vectors are no longer ortho-normal, it becomes necessary to distinguish between contravariant and covariant tensor components. The covariant and contravariant transformation matrices,
${\it\bf\Lambda}$
and
${\it\bf\Lambda}^{\prime }$
, respectively, which define the transformation between the canonical Cartesian frame, and the time-dependent oblique frame, are given by Kobayashi & Yamamoto (Reference Kobayashi and Yamamoto2011)
$$\begin{eqnarray}\displaystyle {\it\bf\Lambda}=\left(\begin{array}{@{}ccc@{}}1 & {\it\gamma}(t) & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{array}\right),\quad {\it\bf\Lambda}^{\prime }=\left(\begin{array}{@{}ccc@{}}1 & -{\it\gamma}(t) & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{array}\right), & & \displaystyle\end{eqnarray}$$
with
${\it\bf\Lambda}\boldsymbol{\cdot }{\it\bf\Lambda}^{\prime }={\it\bf\Lambda}^{\prime }\boldsymbol{\cdot }{\it\bf\Lambda}=\unicode[STIX]{x1D644}$
. In what follows, we will distinguish between tensor components measured in the fixed Cartesian lab frame and the oblique frame, by using a caret or hat (
$\hat{\cdot }$
) for the latter. Thus, the oblique covariant (contravariant) basis sets
$\hat{\boldsymbol{E}}_{{\it\mu}}$
(
$\hat{\boldsymbol{E}}^{{\it\mu}}$
), and the corresponding components
$\hat{x}_{{\it\mu}}$
(
$\hat{x}^{{\it\mu}}$
) of a given vector
$\boldsymbol{x}=x^{{\it\mu}}\boldsymbol{e}_{{\it\mu}}=\hat{x}^{{\it\mu}}\hat{\boldsymbol{E}}_{{\it\mu}}=\hat{x}_{{\it\mu}}\hat{\boldsymbol{E}}^{{\it\mu}}$
, are defined as Carroll (Reference Carroll2004)
$$\begin{eqnarray}\displaystyle \hat{\boldsymbol{E}}_{{\it\mu}}={\it\Lambda}_{\phantom{{\it\nu}}{\it\mu}}^{{\it\nu}}\boldsymbol{e}_{{\it\nu}},\quad \hat{\boldsymbol{E}}^{{\it\mu}}={\it\Lambda}_{\phantom{\prime {\it\mu}}{\it\nu}}^{\prime {\it\mu}}\boldsymbol{e}^{{\it\nu}}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \hat{x}_{{\it\mu}}={\it\Lambda}_{\phantom{{\it\nu}}{\it\mu}}^{{\it\nu}}x_{{\it\nu}},\quad \hat{x}^{{\it\mu}}={\it\Lambda}_{\phantom{\prime {\it\mu}}{\it\nu}}^{\prime {\it\mu}}x^{{\it\nu}}, & & \displaystyle\end{eqnarray}$$
${\it\mu}=1,2,3$
),
$\boldsymbol{e}_{{\it\mu}}=\boldsymbol{e}^{{\it\mu}}$
by definition, and
${\it\Lambda}_{\phantom{{\it\nu}}{\it\mu}}^{{\it\nu}}=[{\it\bf\Lambda}]_{{\it\nu}{\it\mu}}$
and
${\it\Lambda}_{\phantom{\prime {\it\mu}}{\it\nu}}^{\prime {\it\mu}}=[{\it\bf\Lambda}^{\prime }]_{{\it\mu}{\it\nu}}$
denote the matrix components of
${\it\bf\Lambda}$
and
${\it\bf\Lambda}^{\prime }$
(A 1a,b
), respectively. The metric tensors for this oblique coordinate space are then 
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D60E}^{{\it\mu}{\it\nu}} & = & \displaystyle \hat{\boldsymbol{E}}^{{\it\mu}}\boldsymbol{\cdot }\hat{\boldsymbol{E}}^{{\it\nu}}=\left(\begin{array}{@{}ccc@{}}1+{\it\gamma}^{2}(t) & -{\it\gamma}(t) & 0\\ -{\it\gamma}(t) & 1 & 0\\ 0 & 0 & 1\end{array}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D60E}_{{\it\mu}{\it\nu}} & = & \displaystyle \hat{\boldsymbol{E}}_{{\it\mu}}\boldsymbol{\cdot }\hat{\boldsymbol{E}}_{{\it\nu}}=\left(\begin{array}{@{}ccc@{}}1 & {\it\gamma}(t) & 0\\ {\it\gamma}(t) & 1+{\it\gamma}^{2}(t) & 0\\ 0 & 0 & 1\end{array}\right).\end{eqnarray}$$
$\hat{{\rm\nabla}}_{{\it\mu}}=\partial /\partial \hat{x}^{{\it\mu}}$
and
$\hat{{\rm\nabla}}^{{\it\mu}}=\unicode[STIX]{x1D60E}^{{\it\mu}{\it\nu}}\partial /\partial \hat{x}^{{\it\nu}}$
).
When the coordinate system used to solve the Navier–Stokes equation is time-dependent, which is the case we are interested in here, the proper form for the equations of motion is given by Luo & Bewley (Reference Luo and Bewley2004), Venturi (Reference Venturi2009)
$$\begin{eqnarray}\displaystyle \frac{{\it\delta}u^{{\it\mu}}}{{\it\delta}t}={\it\rho}^{-1}{\rm\nabla}_{{\it\nu}}{\it\sigma}^{{\it\nu}{\it\mu}}+{\it\phi}f_{p}^{{\it\mu}}, & & \displaystyle\end{eqnarray}$$
where
${\it\delta}/{\it\delta}t$
denotes the intrinsic time derivative, defined for a contravariant quantity as
$$\begin{eqnarray}\displaystyle \frac{{\it\delta}A^{{\it\mu}}}{{\it\delta}t}=\frac{\partial A^{{\it\mu}}}{\partial t}+(u^{{\it\nu}}-U^{{\it\nu}}){\rm\nabla}_{{\it\nu}}A^{{\it\mu}}+A^{{\it\nu}}{\rm\nabla}_{{\it\nu}}U^{{\it\mu}} & & \displaystyle\end{eqnarray}$$
with
$U^{{\it\mu}}\equiv -\partial x^{{\it\mu}}/\partial t$
the velocity of the coordinate flow. For simple-shear flow, the coordinate flow is defined as the mean shear-flow velocity
$$\begin{eqnarray}\displaystyle \boldsymbol{U}=\dot{{\it\gamma}}(t)y\,\boldsymbol{e}_{1}=\dot{{\it\gamma}}(t)\hat{x}^{2}\hat{\boldsymbol{E}}_{1}. & & \displaystyle\end{eqnarray}$$
In terms of the peculiar velocity,
${\bf\xi}=\boldsymbol{u}-\boldsymbol{U}$
, (A 4) can then be written in the oblique frame as
$$\begin{eqnarray}\displaystyle (\hat{\partial }_{t}+\hat{{\it\xi}}^{{\it\nu}}\hat{{\rm\nabla}}_{{\it\nu}})\hat{{\it\xi}}^{{\it\mu}}={\it\rho}^{-1}\hat{{\rm\nabla}}_{{\it\nu}}\hat{{\it\sigma}}^{{\it\nu}{\it\mu}}+\hat{{\it\phi}}\hat{f}_{p}^{{\it\mu}}-\left(\frac{\text{d}\dot{{\it\gamma}}(t)\,}{\text{d}t\,}\hat{x}^{2}+2\dot{{\it\gamma}}(t)\hat{{\it\xi}}^{2}\right){\it\delta}^{{\it\mu},1}, & & \displaystyle\end{eqnarray}$$
where
$\hat{\partial }_{t}$
gives the partial time derivative at a fixed point in oblique space
$$\begin{eqnarray}\displaystyle \hat{\partial }_{t}\equiv \left.\frac{\partial }{\partial t}\right|_{\hat{x}^{{\it\mu}}}=\left.\frac{\partial }{\partial t}\right|_{x^{{\it\mu}}}+\dot{{\it\gamma}}(t)y\frac{\partial }{\partial x}. & & \displaystyle\end{eqnarray}$$
The last term on the right-hand side of (A 7) arises because of the coordinate flow. It is instructive to consider the meaning of each of these terms independently. One of the factors in
$\dot{{\it\gamma}}(t)\hat{{\it\xi}}^{2}$
comes from the explicit time dependence of the basis vectors, the other from the advection of the coordinate flow, and the term proportional to
$\partial \dot{{\it\gamma}}/\partial t$
from the time dependence of the imposed shear. In the case of oscillatory-shear flow, this last term is problematic, as it depends explicitly on
$\hat{x}^{2}$
, rendering the equation incompatible with the proposed periodic boundary conditions. However, under the assumption that the
$\boldsymbol{k}=0$
component of the fluid velocity instantaneously follows the applied shear, we can safely set this term equal to zero, i.e.
$\hat{\partial }_{t}\hat{U} ^{{\it\mu}}=0$
. The fact that this approximation is necessary is not surprising, as we are dealing with an unbounded system under (Lees–Edwards) periodic boundary conditions. As there is no boundary or external force driving the system, it is not possible to consider the response of the fluid to the imposed flow. Experimentally, this assumption corresponds to situations in which the oscillatory Reynolds number
$Re_{{\it\omega}}={\it\omega}L^{2}/{\it\nu}\ll 1$
(with
${\it\omega}$
the driving frequency and
$L$
the characteristic length over which the flow is varying). In such quasi-steady situations, the Stokes boundary layer
${\it\delta}\simeq \sqrt{{\it\nu}/{\it\omega}}$
is much larger than the system size, and the velocity and stress profiles are equivalent to the steady case
${\it\omega}=0$
(e.g. steady Couette flow). Thus, when considering oscillatory-shear flows,
${\it\delta}$
must be much larger any length scale in our system.
Finally, we arrive at
$$\begin{eqnarray}\displaystyle (\partial _{t}+\hat{{\it\xi}}^{{\it\nu}}\hat{{\rm\nabla}}_{{\it\nu}})\hat{{\it\xi}}^{{\it\mu}}={\it\rho}^{-1}\hat{{\rm\nabla}}_{{\it\nu}}\hat{{\it\sigma}}^{{\it\nu}{\it\mu}}+\hat{{\it\phi}}\hat{f}_{p}^{{\it\mu}}-2\dot{{\it\gamma}}(t)\hat{{\it\xi}}^{2}{\it\delta}^{{\it\mu},1}, & & \displaystyle\end{eqnarray}$$
which gives the time evolution of the (contravariant) components of the peculiar fluid velocity
$\hat{{\it\xi}}^{{\it\mu}}$
, with respect to a time-dependent coordinate system
$\hat{x}^{{\it\mu}}$
, which is being instantaneously advected by the linear flow field
$\boldsymbol{U}$
(A 6). Our final form for the Navier–Stokes equation, (A 9), is analogous to the SLLOD equations of motion for the peculiar momentum of a particle in an arbitrary homogeneous and divergence-free flow (Todd & Daivis Reference Todd and Daivis2007; Evans & Morriss Reference Evans and Morriss2008). In this case, the same approximation of an instantaneous response is assumed. In contrast with (A 7), (A 9) contains no explicit dependence on the coordinates and, as such, can be used together with the Lees–Edwards boundary conditions to describe periodic systems under both steady and oscillatory shear.
Appendix B. Computational Algorithm
In what follows we provided a detailed description of the time-stepping scheme we have used. First, we update the fluid velocity by solving for the advection and viscous stress, and update the particle configuration. We solve (A 9) in the absence of any particle constraint force
$$\begin{eqnarray}\displaystyle \hat{{\it\xi}}^{{\it\mu}\ast }=\hat{{\it\xi}}^{{\it\mu}}+\int _{t_{n}}^{t_{n+1}}\text{d}s\,[\hat{{\rm\nabla}}_{{\it\nu}}({\it\rho}^{-1}\hat{{\it\sigma}}^{{\it\nu}{\it\mu}}-\hat{{\it\xi}}^{{\it\nu}}\hat{{\it\xi}}^{{\it\mu}})-2\dot{{\it\gamma}}(t)\hat{{\it\xi}}^{2}{\it\delta}^{{\it\mu},1}] & & \displaystyle\end{eqnarray}$$
under the incompressibility condition (
$\hat{{\rm\nabla}}_{{\it\mu}}\hat{{\it\xi}}^{{\it\mu}}=0$
). Next, we solve the Newton–Euler equations to obtain the new particle positions and orientations
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{R}_{i}^{n+1}=\boldsymbol{R}_{i}^{n}+\int _{t_{n}}^{t_{n+1}}\text{d}s\,\boldsymbol{V}_{i}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D618}_{i}^{n+1}=\unicode[STIX]{x1D618}_{i}^{n}+\int _{t_{n}}^{t_{n+1}}\text{d}s\,\text{skew}({\it\bf\Omega}_{i})\boldsymbol{\cdot }\unicode[STIX]{x1D618}_{i}. & \displaystyle\end{eqnarray}$$
The fluid updated fluid velocity
$\hat{{\it\xi}}^{{\it\mu}\ast }$
is transformed from the oblique frame to the ortho-normal Cartesian frame. The updated fluid and particle velocity fields, are then
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{\ast }(\boldsymbol{x})=\big(\dot{{\it\gamma}}(t_{n+1})y{\it\delta}^{{\it\mu},1}+{\it\Lambda}_{\phantom{{\it\mu}}{\it\nu}}^{{\it\mu}}\hat{{\it\xi}}^{{\it\nu}\ast }(\boldsymbol{x})\big)\boldsymbol{e}_{{\it\mu}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\phi}\boldsymbol{u}_{p}^{\ast }(\boldsymbol{x})=\mathop{\sum }_{i}{\it\phi}_{i}^{n+1}[\boldsymbol{V}_{i}^{n}+{\it\bf\Omega}_{i}^{n}\times \boldsymbol{r}_{i}^{n+1}], & \displaystyle\end{eqnarray}$$
with
$\boldsymbol{r}_{i}^{n+1}=\boldsymbol{x}-\boldsymbol{R}_{i}^{n+1}$
. After solving the Navier–Stokes equation, the oblique velocity
$\hat{{\it\xi}}^{{\it\nu}\ast }$
is not defined on the Cartesian grid, but instead on the moving oblique grid. We need to interpolate the values between the two grids. We abandon the linear interpolation scheme used previously (Kobayashi & Yamamoto Reference Kobayashi and Yamamoto2011) to map between the oblique and rectangular grid points. Instead, we use a higher-order interpolation, which reduces the fluctuations in the local stress calculations. Thus, a field variable
$f(\boldsymbol{x}_{i,j,k})$
, evaluated at a fixed Cartesian grid point
$\boldsymbol{x}_{i,j,k}={\it\Delta}(i,j,k)$
(
$i,j,k$
integers and
${\it\Delta}$
the grid spacing), is interpolated using the values at all
$N_{x}$
oblique grid points
$\hat{\boldsymbol{x}}_{\hat{\imath },\hat{\jmath },\hat{k}}$
, along the flow direction
$\boldsymbol{e}_{x}$
$$\begin{eqnarray}\displaystyle f(\boldsymbol{x}_{i,j,k})=\mathop{\sum }_{\hat{\imath }=0}^{N_{x}}{\it\alpha}(i\leftarrow \hat{\imath })f(\hat{\boldsymbol{x}}_{\hat{\imath },\hat{\jmath },\hat{k}}), & & \displaystyle\end{eqnarray}$$
where, according to (A 2b
),
$\hat{\jmath }=j$
and
$\hat{k}=k$
, and the interpolation weights
${\it\alpha}(i\leftarrow \hat{\imath })$
are chosen to obtain a periodic cubic spline fit with continuous first and second derivatives (Hildebrand Reference Hildebrand1987).
Second, we compute the hydrodynamic force
$\boldsymbol{F}^{H}$
and torque
$\boldsymbol{N}^{H}$
exerted by the fluid on the particles. Assuming momentum conservation, the time-integrated hydrodynamic force and torque are given directly by the momentum exchange over the particle domain
$$\begin{eqnarray}\displaystyle & \displaystyle \left[\int _{t_{n}}^{t_{n+1}}\text{d}s\,\boldsymbol{F}_{i}^{H}\right]=\int \text{d}\boldsymbol{x}\,{\it\rho}{\it\phi}_{i}^{n+1}(\boldsymbol{u}^{\ast }-\boldsymbol{u}_{p}^{\ast }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left[\int _{t_{n}}^{t_{n+1}}\text{d}s\,\boldsymbol{N}_{i}^{H}\right]=\int \text{d}\boldsymbol{x}\,[\boldsymbol{r}_{i}^{n+1}\times {\it\rho}{\it\phi}_{i}^{n+1}(\boldsymbol{u}^{\ast }-\boldsymbol{u}_{p}^{\ast })]. & \displaystyle\end{eqnarray}$$
Third, we update the particle velocities and angular velocities
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{V}_{i}^{n+1}=\boldsymbol{V}_{i}^{n}+M_{i}^{-1}\int _{t_{n}}^{t_{n+1}}\text{d}s\,[\boldsymbol{F}_{i}^{H}+\boldsymbol{F}_{i}^{C}+\boldsymbol{F}_{i}^{Ext}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\bf\Omega}_{i}^{n+1}={\it\bf\Omega}_{i}^{n}+\unicode[STIX]{x1D610}_{i}^{-1}\boldsymbol{\cdot }\int _{t_{n}}^{t_{n+1}}\text{d}s\,[\boldsymbol{N}_{i}^{H}+\boldsymbol{N}_{i}^{\text{C}}+\boldsymbol{N}_{i}^{Ext}], & \displaystyle\end{eqnarray}$$
which gives the final particle velocity field
$\boldsymbol{u}_{p}^{n+1}$
, where
$\boldsymbol{F}^{C}$
(
$\boldsymbol{N}^{C}$
) and
$\boldsymbol{F}^{Ext}$
(
$\boldsymbol{N}^{Ext}$
) represent the inter-particle and external forces (torques). Without loss of generality, we will neglect all external fields, such that
$\boldsymbol{F}_{i}^{Ext}=\boldsymbol{N}_{i}^{Ext}=0$
. For what follows, we also define the hydrodynamic and colloidal velocity increments as
$$\begin{eqnarray}\displaystyle & \displaystyle {\rm\Delta}\boldsymbol{V}_{i}^{X}=M_{i}^{-1}\int _{t_{n}}^{t_{n+1}}\text{d}s\,\boldsymbol{F}_{i}^{X}\quad (X=\text{H, C}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\rm\Delta}{\it\bf\Omega}_{i}^{X}=\unicode[STIX]{x1D610}_{i}^{-1}\boldsymbol{\cdot }\int _{t_{n}}^{t_{n+1}}\text{d}s\,\boldsymbol{N}_{i}^{X} & \displaystyle\end{eqnarray}$$
such that
$\boldsymbol{V}_{i}^{n+1}=\boldsymbol{V}_{i}^{n}+{\rm\Delta}\boldsymbol{V}_{i}^{H}+{\rm\Delta}\boldsymbol{V}_{i}^{C}$
, and likewise for the angular velocity update. Finally, we impose the particle rigidity on the total fluid velocity through the body force
${\it\phi}\boldsymbol{f}_{p}$
in the Navier–Stokes equation
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{n+1}=\boldsymbol{u}^{\ast }+\left[\int _{t_{n}}^{t_{n+1}}\text{d}s\,{\it\phi}\boldsymbol{f}_{p}\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left[\int _{t_{n}}^{t_{n+1}}\text{d}s\,{\it\phi}\boldsymbol{f}_{p}\right]={\it\phi}^{n+1}(\boldsymbol{u}_{p}^{n+1}-\boldsymbol{u}^{\ast }) & \displaystyle\end{eqnarray}$$
and we again impose the divergence-free condition on the total velocity field
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{n+1}=0$
. The new velocity field is then transformed back to oblique space, and the procedure is repeated for the subsequent time steps. As pointed out by Onuki (Reference Onuki1997), whenever
$\dot{{\it\gamma}}t=1$
, we should reset the grid and the strain (
${\it\gamma}(t)\rightarrow \dot{{\it\gamma}}t-1$
) in order to maintain numerical stability. This procedure introduces additional problems of its own, such as the loss of kinetic energy and dissipation rate, but this is significant only at very high shear rates, and is not a problem for the systems considered in this work (Brucker et al.
Reference Brucker, Isaza, Vaithianathan and Collins2007).
Appendix C. Calculation of particle stress
Within the SPM, the particle contribution to the stress is given in terms of the body force
$\boldsymbol{f}_{p}$
as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D668}=\frac{1}{V}\int \text{d}\boldsymbol{r}\,[{\it\delta}{\bf\sigma}-\boldsymbol{r}{\it\rho}{\it\phi}\boldsymbol{f}_{p}+\boldsymbol{r}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}({\it\rho}\boldsymbol{u})]. & & \displaystyle\end{eqnarray}$$
It is convenient to separate the contributions arising from the hydrodynamic interactions, from those arising from the direct inter-particle colloidal forces, such that
${\it\phi}\boldsymbol{f}_{p}={\it\phi}\boldsymbol{f}_{p}^{H}+{\it\phi}\boldsymbol{f}_{p}^{C}$
. By definition, we have
$$\begin{eqnarray}\displaystyle \left[\int _{t_{n}}^{t_{n+1}}\text{d}s\,{\it\phi}\boldsymbol{f}_{p}\right] & = & \displaystyle {\it\phi}^{n+1}(\boldsymbol{u}_{p}^{n+1}-\boldsymbol{u}^{\ast })\nonumber\\ \displaystyle & = & \displaystyle {\it\phi}^{n+1}(\boldsymbol{u}_{p}^{n+1}-\boldsymbol{u}_{p}^{\ast })-{\it\phi}^{n+1}(\boldsymbol{u}^{\ast }-\boldsymbol{u}_{p}^{\ast }).\end{eqnarray}$$
The first term on the right-hand side can then be expressed in terms of the hydrodynamic and colloid (inter-particle) forces, as
$$\begin{eqnarray}\displaystyle {\it\phi}^{n+1}(\boldsymbol{u}_{p}^{n+1}-\boldsymbol{u}_{p}^{\ast }) & = & \displaystyle \mathop{\sum }_{i}{\it\phi}_{i}^{n+1}[{\rm\Delta}\boldsymbol{V}_{i}^{H}+{\rm\Delta}{\it\bf\Omega}_{i}^{H}\times \boldsymbol{r}_{i}^{n+1}]\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{i}{\it\phi}_{i}^{n+1}[{\rm\Delta}\boldsymbol{V}_{i}^{C}+{\rm\Delta}{\it\bf\Omega}_{i}^{C}\times \boldsymbol{r}_{i}^{n+1}].\end{eqnarray}$$
The individual contributions to the body force, to first order in time, are then given by
$$\begin{eqnarray}\displaystyle & \displaystyle h{\it\phi}\boldsymbol{f}_{p}^{H}=\mathop{\sum }_{i}[{\rm\Delta}\boldsymbol{V}_{i}^{H}+{\rm\Delta}{\it\bf\Omega}_{i}^{H}\times \boldsymbol{r}_{i}^{n+1}]-{\it\phi}^{n+1}(\boldsymbol{u}^{\ast }-\boldsymbol{u}_{p}^{\ast }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h{\it\phi}\boldsymbol{f}_{p}^{C}=\mathop{\sum }_{i}{\it\phi}_{i}^{n+1}[{\rm\Delta}\boldsymbol{V}_{i}^{C}+{\rm\Delta}{\it\bf\Omega}_{i}^{C}\times \boldsymbol{r}_{i}^{n+1}]. & \displaystyle\end{eqnarray}$$
Appendix D. Calculation of local stress field
When the colloidal dispersion is simulated using the SPM, the Navier–Stokes equation for the total fluid velocity is solved using a Fourier spectral scheme in oblique space. As the velocity field is given in this oblique (dual) wavevector space, it is computationally more efficient to directly compute the local stress (pressure) field in this space as well, and then perform an inverse Fourier transform to recover the real components of the stress over the fixed Cartesian grid. Based on the definition of the basis vectors and the tensorial quantities defined in (A 2a
) and (A 2b
), the contravariant components of the wavevectors
$\hat{k}^{{\it\mu}}$
are represented in terms of covariant wavevectors
$\hat{k}_{{\it\mu}}=2{\rm\pi}n/L_{{\it\mu}}$
, where
$n$
is an integer and we assume a rectangular simulation box of dimensions
$L_{{\it\mu}}$
, as
$\hat{k}^{{\it\mu}}={\hat{G}}^{{\it\mu}{\it\nu}}\hat{k}_{{\it\nu}}$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{k}^{1}=[1+{\it\gamma}^{2}(t)]\hat{k}_{1}-{\it\gamma}(t)\hat{k}_{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{k}^{2}=\hat{k}_{2}-{\it\gamma}(t)\hat{k}_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{k}^{3}=\hat{k}_{3}. & \displaystyle\end{eqnarray}$$
Application of the incompressibility condition in (2.5b ) to the Navier–Stokes equation (A 9) leads to
$$\begin{eqnarray}\displaystyle \hat{{\rm\nabla}}_{{\it\mu}}\left[\hat{{\rm\nabla}}_{{\it\nu}}(\hat{{\it\xi}}^{{\it\mu}}\hat{{\it\xi}}^{{\it\nu}})+\frac{{\it\eta}}{{\it\rho}}\hat{{\rm\nabla}}^{2}\hat{{\it\xi}}^{{\it\mu}}+2\dot{{\it\gamma}}\hat{{\it\xi}}^{2}{\it\delta}^{{\it\mu},1}\right]+\hat{{\rm\nabla}}^{2}\frac{\hat{p}}{{\it\rho}}=0, & & \displaystyle\end{eqnarray}$$
which allows us to calculate the pressure
$\hat{p}(\hat{k})$
as
$$\begin{eqnarray}\displaystyle \hat{p}(\hat{k})=\text{i}{\it\rho}\frac{\hat{k}_{{\it\mu}}}{\hat{k}^{{\it\lambda}}\hat{k}_{{\it\lambda}}}\left(\text{i}\hat{k}_{{\it\nu}}\hat{V}^{{\it\mu}{\it\nu}}+\frac{{\it\eta}}{{\it\rho}}\hat{k}^{{\it\nu}}\hat{k}_{{\it\nu}}\hat{{\it\xi}}^{{\it\mu}}+2\dot{{\it\gamma}}\hat{{\it\xi}}^{2}{\it\delta}^{1,{\it\mu}}\right), & & \displaystyle\end{eqnarray}$$
where
$\hat{V}^{{\it\mu}{\it\nu}}(\hat{\boldsymbol{k}})=\int \hat{{\it\xi}}^{{\it\mu}}\hat{{\it\xi}}^{{\it\nu}}\exp (-\text{i}\hat{k}_{{\it\alpha}}\hat{x}^{{\it\alpha}})\text{d}\hat{\boldsymbol{x}}$
is the Fourier transform of the dyadic product of the (peculiar) velocity field with itself. As usual for incompressible systems, the pressure
$\hat{p}$
itself is meaningless, it just serves to obtain a solenoidal velocity field.
