Predictions of microstructure and stress in planar extensional flows of a dense viscous suspension

In a recent study, Jenkins and La Ragione (2015) determine the typical trajectory of a force-equilibrated pair of particles of a dense, two-dimensional suspensions of spheres subjected to a simple shearing flow. They evaluate the distribution function of such near-contacting neighbours along the trajectory and, using this distribution function and the expression for the force between the pair, they predict the particle pressure, the difference in normal stresses and the difference between the average rotation of the spheres and half the vorticity of the average velocity.

Here, we focus on extensional flows, also called pure shearing, of a dense layer of spheres and, as an extension of the previous work, also introduce the moment equilibrium. We employ a simplified Stokesian dynamics numerical simulation, perhaps more properly called lubrication dynamics (e.g., Ball and Melrose 1995), to guide the choice of the near-contacting pairs on a representative trajectory that contributes most to the inter-particle stress. We specify the boundary conditions on the representative trajectory, and determine the distribution of particles along it and the relationship between stress, microstructure and strain rate. We test these predictions against the results of the numerical simulations. We show that the relation between the stress and strain rate tensors involves the second and fourth moments of the particle distribution, and place this and other aspects of our approach in the context of the earlier models of Phan-Thien (1995), Stickel, et al. (2006), Goddard (2006), Gillissen and Wilson (2018, 2019) and Gilissen et al. (2019) that focus on the second moment and that of Chacko et al. (2018), who introduce a fourth-rank tensor to describe flow reversal.

The approximate satisfaction of force and torque balances for particles in the flow plays an important role in what we do. In that regard, we operate in the spirit of Nazockdast and Morris (2012a, 2012b, 2013) or that of the statistical characterization by Thomas et. al. (2018) of equilibrated particles sheared in two dimensions, but in the limit of dense flows of the planar extensional flow. The analysis must be extended to three-dimensional simple shear flows before it can be placed in relation to phenomenological relations that have resulted from experiments on dense three-dimensional shearing flows (Boyer, et al. 2011; Guazzelli and Pouliquen 2018).

A steady, planar extensional flow of a dense suspension of identical spheres with radius $a$ is characterized by an average rate of deformation tensor $\mathbfsf{D}$ with non-zero components $\mathsfit{D}_{11}=-\mathsfit{D}_{22}=\dot{\gamma}$, where $x_{1}$ and $x_{2}$ are the axes in the directions of greatest extension and compression, respectively, and $\dot{\gamma}$ is the constant shear rate. We focus on a typical pair of spheres and their near-contacting neighbours, and take $\boldsymbol{\hat{d}}^{(BA)}$ to be the unit vector directed from the centre of sphere $A$ to that of sphere $B$, with $\boldsymbol{\hat{d}}^{(AB)}=-\boldsymbol{\hat{d}}^{(BA)}$ (see Fig. 1). Then, with $\theta^{(BA)}$ the time-dependent angle between $\boldsymbol{\hat{d}}^{(BA)}$ and the $x_{2}$ axis,

and the components of the unit tangent vector, $\boldsymbol{\hat{t}}^{(BA)}=-\boldsymbol{\hat{t}}^{(AB)}$, perpendicular to it, are

or $t_{\alpha}^{(BA)}=\varepsilon_{\alpha\beta}d_{\beta}^{(BA)}$, where $\varepsilon_{12}=-\varepsilon_{21}=1$ and $\varepsilon_{11}=\varepsilon_{22}=0$. The unit vectors $\boldsymbol{\hat{d}}^{(BA)}$ and $\boldsymbol{\hat{t}}^{(BA)}$ are indicated in Fig. 1.

Knowledge of the distribution of near-contacting neighbours $A(\theta)$ and the contact forces along the trajectory permits the calculation of the macroscopic particle stress in the suspension. The stress tensor is, according to Cauchy (Love 1944, Appendix, Note B),

where $n$ is the number of particles per unit area and $F_{\alpha}$ is given by Eq. (30). Then, the two-dimensional viscosity is $2a\mu$. The dimensionless form, $\mathsfit{t}_{\alpha\beta}=\mathsfit{T}_{\alpha\beta}/(2a\mu\dot{\gamma})$, with $n=c/(\pi a^{2})$, is

The particle shear stress,

is

where $b$, $c_{1}$ and $c_{2}$ are given in terms of $k$ in Eqs. (24) and (28), respectively. The shear stress depends on the separation, $\bar{s}$, of near-contacting neighbours other than $B$, and on the area fraction, explicitly and through the coordination number, $k$. Because the direct contribution of the repulsive force to the integral is very small and the trigonometric factors associated with the other contributions are even about $\pi/4$, the shear stress is independent of the asymmetry of the particle distribution about $\pi/4$. In contrast, this asymmetry is crucial to the determination of the particle pressure.

The particle pressure,

is

This pressure also depends on $\bar{s}$ and $c$ and its existence is due to the asymmetry of $A$ about $\pi/4$. This asymmetry is due to that of the separation along the representative trajectory created by $\hat{F}$ and the influence of the asymmetry of the separation on the angular velocity, $\dot{\theta}$. The particle pressure and the mechanisms responsible for it are a focus of this paper; a particle shear stress may be calculated based on the average flow, although that determined here is several times less than this, because of the approximate satisfaction of equilibrium.

Particle stresses associated with motion along the representative trajectories are compared with those measured in Stokesian dynamics simulations after a discussion of the simulations.

We determine the trajectories of spherical particles in the flows by performing simulation with the same conditions as the theory (a monolayer with no inertia). We impose a planar extensional flow with shear rate $\dot{\gamma}$, ¹¹1Note: This is equivalent to the extensional rate $\dot{\varepsilon}$ in Seto et al., 2017.

A simulation box with periodic boundary conditions constantly deforms according to this velocity gradient $\mathbfsf{D}$. Significant deformations of the simulation box can be avoided by using the Kraynik–Reinelt periodic boundary conditions, which rearrange the deformed box to the original square box after a constant strain interval (Kraynik and Reinelt 1992, Todd and Daivis, 1998, Seto et al., 2017). Thus, the flow can be applied for a sufficiently long time to evaluate its steady states.

Due to the negligible inertia of the particles, translational and angular velocities can be determined by solving the force and torque balance equations for the respective particles ($i=1,\dotsc,N$):

Here, a vector, such as $\boldsymbol{F}_{\mathrm{H}}$, represents all $N$ particles, $\boldsymbol{F}_{\mathrm{H}}\equiv(\boldsymbol{F}_{\mathrm{H}}^{(1)},\dotsc,\boldsymbol{F}_{\mathrm{H}}^{(N)})$.

The hydrodynamic interactions in the Stokes, zero Reynolds number, regime are linear in the velocities

where $\mathbfsf{D}_{\mathrm{N}}$ is block diagonal of $N$ copies of $\mathbfsf{D}$. There exist several levels of approximations to construct the resistance matrices $\mathbfsf{R}_{\mathrm{FU}}$ and $\mathbfsf{R}_{\mathrm{FD}}$. Brady and Bossis (1988) constructed them using truncated multipole expansions for the far-field interactions and a pairwise solution for lubrication interactions. In this work, we focus on a special situation in which repulsive forces are very weak in comparison with viscous drag forces. Under such conditions, particles tend to approach their neighbours very closely.

Because the resistance coefficients diverge at contacts ($s=0$), the nearly touching hydrodynamic interactions dominate the dynamics. This is why we construct the approximate resistance matrices with the leading $1/s$ term in the normal component and the logarithmic term $\log(1/s)$ and following constants in the tangential component, using the solution for two nearly touching rigid spheres (Jeffrey and Onishi, 1984, Jeffrey, 1992). (A detailed description can be seen elsewhere, c.f. Mari el al., 2014.) The hydrodynamic interaction is effective only when $s<0.10$; thus, the resistance coefficients remain positive in this range.

The repulsive force employed in this work is the same as that used in Nott and Brady (1994):

where the range of repulsive force is set by a parameter $\lambda$. Because the repulsive force diverges as $F_{0}/s$ in the limit of contact, $s\to 0$, some force balance can occur at a finite gap. Because the repulsive force diverges as $F_{0}/s$ in the limit of contact, $s\to 0$, some force balance can occur at a finite gap. Thus, the gap $s$ remains positive, and contact forces do not appear in the current system. Note that the divergence in the lubrication coefficient does not guarantee the presence of a minimum $s>0$, thus it leads to a pathologic singularity in theoretical models (Ball and Melrose 1995).

By solving the force and torque balance equations (46) with the hydrodynamic interaction (47) and repulsive force (48), the linear and angular velocities $(\boldsymbol{U},\boldsymbol{\Omega})$ can be determined at each time step. Integrating these velocities $\boldsymbol{U}$ with a discretized time step, we obtain trajectories of particles. The particle stress tensor $\hat{\mathbfsf{T}}$ is given by the symmetrized first moment

with the pairwise forces $\boldsymbol{F}^{ij}\equiv\boldsymbol{F}^{ij}_{\mathrm{Lub}}+\boldsymbol{F}^{ij}_{\mathrm{R}}$ and relative positions $\boldsymbol{r}^{ij}\equiv\boldsymbol{r}^{i}-\boldsymbol{r}^{j}$ of all interacting particle pairs. Here, $V=2aL^{2}$ is the volume of the monolayer system. Normalizing the symmetrized first moment with the shear stress of the suspending fluid $2\mu\dot{\gamma}$, gives the dimensionless stress $\mathsfit{t}_{\alpha\beta}\equiv\hat{\mathsfit{T}}_{\alpha\beta}/2\mu\dot{\gamma}$. Thus, we have the dimensionless particle pressure $p\equiv-(\mathsfit{t}_{11}+\mathsfit{t}_{22})/2$ and the dimensionless particle shear stress $\tau\equiv(\mathsfit{t}_{11}-\mathsfit{t}_{22})/2$, respectively.

We simulate monolayer systems with 1000 spheres of radius $a$ at area fractions, $c$, between $0.52$ and $0.64$. We generate initial configurations with a simple algorithm using random numbers. To reduce effects of such artificially generated initial configurations, the post-processing analyses use steady state data from 10 to 50 strain units. The repulsive force is set to be very weak $F_{0}/\pi a^{3}\mu\dot{\gamma}=10^{-4}$ and short-ranged $\lambda/a=10^{-2}$.

In the implementation of the model, we take $\bar{s}=0.02$, $\theta_{0}=10^{-6}$, $s_{0}=s_{1}=0.10$ and assume that $k$ varies linearly with $c$ from 2.0, at $c=0.52,$ to 2.5, at $c=0.64$. These values and the relation for $k$ are plausible and they are influenced by those measured in the simulations. The value $\bar{s}$ gives values of the shear stresses that are close to those measured. The value of $\theta_{0}$ is the default absolute tolerance of the solver; smaller values of $\theta_{0}$ have little influence on the shear stress, but do slightly improve the prediction of the pressure. The initial and final values of the separation were those employed in the simulations, and the variation of near-contacting neighbours with concentration, $k(c)$, was that measured in the simulation.

Fig. 2 shows plots of the shear stress $\tau$ and pressure $p$ measured in the simulation and predicted by the model, over a range of area fraction $c.$ The stresses in the simulation increase in a similar manner with $c$, but the ratio $\tau/p$ decreases gradually. The predicted particle pressure is somewhat less than that measured in the simulations and the predicted shear stress is somewhat greater. The ratio of shear stress to pressure decreases with area fraction, as in the numerical simulations; but, because of under- and over-predicting, we have a greater value for the ratio.

Most of stress measured in the simulations is generated by closely approaching particles - defined as those with a separation less than one per cent of the particle radius. As seen in Fig. 3, more than 90% of shear stress is generated from particle pairs with $s<0.01$. Moreover, such near-contacting particles generates almost 100% of the pressure $p$. Finally, approximately 80 % of the shear stress $\tau$ comes from the normal force.

We can check the concentration of stress contribution in the very narrow range of $s$ using distribution maps. We calculate the spatial distribution in $\xi\equiv\log s$. The statistics are calculated with discretized bins $\xi_{k}\equiv\xi_{1}+(k-1)\Delta\xi$, $k=1,\dotsc,k_{\mathrm{max}}$. $\xi_{1}=\log{10^{-7}}$ and $\xi_{k_{\mathrm{max}}}=\log{10^{-1}}$. The results are plotted with $s$ in a logarithmic scale. Fig. 4(a) displays the distribution of shear stress, indicating that the stress tends to be concentrated near the stagnation point $(\theta,s)\approx(0,10^{-6})$. The region of concentration spreads until $\theta\sim\pi/4$. We also separately calculate the stress of (49) constructed with normal forces $\boldsymbol{F}^{ij}_{\mathrm{n}}\equiv\boldsymbol{F}^{ij}\cdot\boldsymbol{n}^{ij}\boldsymbol{n}^{ij}$ and tangential forces $\boldsymbol{F}^{ij}_{\mathrm{t}}\equiv\boldsymbol{F}^{ij}-\boldsymbol{F}^{ij}\cdot\boldsymbol{n}^{ij}\boldsymbol{n}^{ij}$. As shown in Fig. 3(b), 80 % of the shear stress $\tau$ indeed comes from the normal forces. Besides systematic motions due to the shearing deformation, particle motions fluctuate due to occasional configurations of surrounding particles.

Therefore, it is necessary to reconstruct averaged trajectories to compare with theoretical ones. To this end, we first calculate the averaged relative-velocity field $\langle\boldsymbol{U}^{(j)}-\boldsymbol{U}^{(i)}\rangle$ over all interacting pairs $i$ and $j$ in terms of the relative position coordinate $\Delta\boldsymbol{r}^{ij}=(2a+s,\theta)$. Owing to the symmetry of planar extension, the statistics are taken on a quadrant: $0<\theta<\pi/2$. Because we consider a situation that is very close to the singularity (Ball and Melrose 1995), the particles tend to approach very close to contact, i.e., a bundle of trajectories is compressed into an extremely narrow range of $s$. To avoid a loss of precision due to averaging, we carry out the statistical data binning with $\xi$ instead of $s$.

Once we evaluate the velocity field in the $\xi$–$\theta$ space, i.e., $(\langle\dot{\xi}\rangle=\langle\dot{s}/s\rangle,\langle\dot{\theta}\rangle)$, we can obtain trajectories as streamlines of the velocity field. In Fig. 4(b), trajectories of the system with $c=0.52$ and various initial positions are plotted. The trajectories are colored from blue to red, according to their contribution to the stress. We identify the band of significant trajectories as those with a separation of less than $10^{-2}$.

In Fig. 5, we show the particle number density, measured in the simulations for $c=0.64,$ with particle trajectories superposed. The near-contact distribution $A(\theta),$ normalized by the total number of particles, is obtained by integrating across the trajectories in the region below $s=10^{-2}$. The product of the integral of the particle probability distribution over the area of Fig. 5 below $s=10^{-2}$ and the total number of particles is the coordination number.

In Fig. 6(a), we plot the trajectories from the simulation, for the smallest separation in Fig. 4(b), and the predicted representative trajectories of the model for concentrations of 0.52 and 0.64. The representative trajectories are located within the band and have a shape similar to the individual trajectories at smaller separations. In Fig. 6(b), we show distributions, $A(\theta)$, measured in the simulation and predicted along the representative trajectories for two values of the concentration. These share the same features and have a similar agreement as the trajectories. The asymmetry of the distributions result from the influence of the repulsive force that breaks the symmetry of approach and departure.

We have employed information from the simulation on the variation in the coordination number as a function of concentration and the value of separation of near-contacting neighbours necessary to reproduce the measured particle shear stress. These, used in the model, gives it the capacity to generate particle trajectories that are representative of the stress-producing trajectories of the simulation and particle distributions along the representative trajectories with the appropriate asymmetry about $\pi/4$ to predict a reasonable variation of a particle pressure. .

We next indicate how the structure of the model can be used as the basis for a continuum theory of dense suspensions and to provide a context for existing theories.

As elaborated upon by Onat and Leckie (1988) and Advani and Tucker (1987,1990), the distribution of near-contacting neighbours can be represented by an infinite series with respect to basis functions, such as

and

The coefficients $\mathbfcal{Z}$ and $\mathbfcal{B}$ are completely traceless and completely symmetric tensors, related to the distribution through

and

These are the second and fourth moments of the distribution, respectively.

The stress of Eq.39) may be written in terms of these as

where

or, more compactly, as

where

The stress depends only on the second and fourth moments of the distribution, although an approximation of the distribution in terms of these does not provide a good representation of it. Because we predict the distribution function, we are able to capture the essential role played by the fourth moment. This places our theory in the context of the work of Chacko, et al. (2018), in which numerical simulations confirm the need of the fourth moment, in addition to the second moment, to describe stress reversal.

With knowledge of the distribution function, it is possible to evaluate the components of the second and fourth moments. In particular, when $F=0$, the only non-zero components are $\mathcal{B}_{1111}=\mathcal{B}_{2222}=-\mathcal{B}_{1122}$; when $F\neq 0$, then $\mathcal{Z}_{11}=-\mathcal{Z}_{22}$ are also different from zero. For example, when $c=0.64$, $k=2.50$ and $\hat{F}=10^{-4}$, their numerical values are $\mathcal{B}_{1111}=0.25$ and $\mathcal{Z}_{11}=-0.57$. Then, with Eq. (56), the dimensionless particle pressure is

Given the force and torque balances, it is possible to characterize the role played by the torque balance in determining features of the trajectory and the distribution of particles along it. Ignoring the torque balance is equivalent to taking $S=0$, or $c_{2}=0$ in Eq. (27). This has an important influence on $\dot{\theta}$. Then, because both the distribution of near-contacting neighbours and the interval over which it is defined depend upon $\dot{\theta}$, there is a dependence of the stress upon it. For example, when $c=0.64$ and the torque balance is ignored, $p$ is $14.81,$ rather than $8.51$. That is, the value of $p$ is affected by the balance through the distribution. In contrast, $\tau$ is 18.60, rather than 18.94 and changes little, because it is independent of the shape of the distribution.

In the absence of the knowledge of the distribution function, it is possible to develop evolution equations for the approximate determination of its moments (e.g., Prantil et al., 1993). Phan-Thien (1995) and Stickel et al. (2006) employ such an equation for the second moment, and break the symmetry of approach and departure by including a term in it that is proportional to $[\text{tr}(\mathbfsf{D}^{2})/2]^{1/2}$. Goddard (2006) introduces a memory integral for the second moment – a representation for the solution to its evolution equation – that breaks this symmetry by incorporating a term proportional to the $\mathbfsf{D}^{2}$. Gillissen and Wilson (2018, 2019) and Gilissen et al. (2019) distinguish between an anisotropic distribution of near contacts and an isotropic distribution of outer contacts, introduce a difference in association and disassociation of these contacts in the directions of compression and extension. This difference appears in the evolution equation of the second moment and provides the asymmetry necessary for normal stress differences. Theories of this type produce stress relations that are linear in the strain rate; that is, rate dependent. In contrast, we employ only a short-range repulsive force that is independent of the shear rate. Consequently, our stress relation contains contributions that are independent of rate.

We have considered a planar extensional flow of a dense layer of spheres in a viscous fluid. In addition to the viscous forces associated with the flow, we assumed that there was a short-range repulsive force between the spheres. We focused on pairs of spheres, assumed that their neighbours translate with the average flow, and required that they be in force and moment equilibrium with each other and their neighbours. We then assumed that the neighbourhoods of pairs with the same orientation were equal to the average over that orientation; this permitted us to write equations for the radial and angular velocity of the relative motion of a single pair as they began and ended an interaction, and the orientation of their line of centres with the axis of greatest compression of the flow.

The possible determination of the distribution of near-contacting particles along the trajectory then leads to expressions for the particle shear stress and pressure. Stokesian dynamics simulations provided a value of particle shear stress that permitted the determination of the angle of departure for the trajectories and the separation of near-contacting neighbours. The variation of the shear stress with area fraction suggested the variation of the number of near-contacting neighbours per particle with area fraction. With this information, numerical values of the particle distribution and the particle pressure could be calculated.

The simplicity of the micro-mechanical model makes it possible to understand the influence of the normal and tangential components of the viscous force and the short-range repulsive force on the shape of the trajectory and the distribution of particles along it. It also permits the derivation of a continuum theory for dense suspensions based on the micro-mechanics of equilibrated pairs of particles and a better understanding of those based on phenomenology. However, it is important to note that correlated motions of the particles along the axis of compression in the simulations have no counterpart in the model; these are likely to force approaching particles into trajectories with significantly smaller separations. The small value of $\theta_{0}$ necessary in the model may be a way for it to incorporate the influence of such particles. Other modelling issues that remain involve the parameters $k$ and ${\bar{s}}$. We believe that it is likely that measurements of the relationship between $k$ and $c$ in three dimensions will be robust and apply to different sharing flows; the determination of ${\bar{s}}$ is less certain, although methods employed by Nazockdast and Morris (2013) may permit its prediction.

The two-dimensional, simple shear flow that was considered earlier by Jenkins and La Ragione (2015) is much more complicated. In it, there are upstream and downstream trajectories, roughly, on either side of the direction of greatest rate of compression. The exact location of the limiting trajectory is influenced by the rotation of the average flow and, in the calculation, its location is difficult to determine. In the earlier study, the separation, was taken to be the average separation calculated by Torquato for a dense, equilibrated system of colliding, elastic spheres in two dimensions. This varies between 0.13 and 0.06 as c varies from 0.52 and 0.64. These values are much greater than the value 0.02 employed here. However, in the calculation for simple shear we obtain an excellent agreement between theory and simulation for the shear stress because the length of the particle interaction was taken to be much longer, and probably compensated for the greater separation. Given our experience with the simpler planar extension, we should return to the simple shear flow.