Jamming in hard sphere suspensions under shear

TL;DR

This paper investigates how shear affects the structure and jamming transition of monodisperse hard sphere suspensions, revealing that shear-induced distortions influence the transition to a jammed state similar to random close packing.

Contribution

It introduces a theoretical framework to analyze shear-induced structural distortions and their impact on jamming in hard sphere suspensions.

Findings

Shear distorts the structure factor of hard spheres at low shear rates.

Distortions influence the transition to a jammed state.

Jamming under shear resembles the random close packing transition.

Abstract

We consider a system of monodisperse hard spheres immersed in a sheared fluid. We obtain the distortion of the structure factor of the hard spheres at low shear rates, within a Percus-Yevick like framework. The consequent distortion of the pair distribution function is shown to affect the transition of the hard sphere fluid into a jammed state, which is similar to the transition to the state of random close packing in the absence of shear.

Equations

$$\gamma \left[{\dot{\mathbf{r}}}_i-\mathbf{V}({\mathbf{r}}_i)\right]=-\nabla \sum _{j=1}^{N}w({\mathbf{r}}_i-{\mathbf{r}}_j)+{\mathbf{\eta }}_i(t),$$

$$\frac{\partial P}{\partial t}=OP,$$

$$O=\frac{1}{\gamma }\left\{\sum _{i,\text{k}}\frac{\partial }{\partial r_{i\text{k}}}\left[kT\frac{\partial }{\partial r_{i\text{k}}}+\frac{\partial U}{\partial r_{i\text{k}}}-\gamma V_{\text{k}}({\mathbf{r}}_i)\right]\right\},$$

$$U=1/2\sum _{i\ne j}w(|{\mathbf{r}}_i-{\mathbf{r}}_j|).$$

$${\rho }_{\mathbf{q}}=\frac{1}{\sqrt{N}}\sum _i\exp (-i\mathbf{q}\cdot {\mathbf{r}}_i),$$

$$\begin{array}{rl}\frac{kT}{\gamma }\sum _{\mathbf{k}}& \{{\text{k}}^2\left[1+\frac{\overline{\rho }}{kT}w(\mathbf{k})+{\rho }_{\mathbf{k}}\frac{\partial }{\partial {\rho }_{\mathbf{k}}}\right]-\nu {\text{k}}_z\frac{\partial {\rho }_{\mathbf{k}}}{\partial {\text{k}}_x}\frac{\partial }{\partial {\rho }_{\mathbf{k}}}+\\ & +N^{-1/2}\sum _{\mathbf{l}}\mathbf{k}\cdot \mathbf{l}[-{\rho }_{\mathbf{k}+\mathbf{l}}\frac{{\partial }^2}{\partial {\rho }_{\mathbf{k}}\partial {\rho }_{-\mathbf{k}}}+\\ & +\frac{\overline{\rho }}{kT}w(\mathbf{l}){\rho }_{\mathbf{k}-\mathbf{l}}{\rho }_{\mathbf{l}}\frac{\partial }{\partial {\rho }_{\mathbf{k}}}]\}\hat{P}=0,\end{array}$$

$$\begin{array}{rl}& \{\sum _{\mathbf{k}}{\text{k}}^2[1+\frac{\overline{\rho }}{kT}w(\mathbf{k})+\left(1+\frac{\overline{\rho }}{kT}w(\mathbf{k})\right){\rho }_{\mathbf{k}}\frac{\partial }{\partial {\rho }_{\mathbf{k}}}+\\ & +\frac{{\partial }^2}{\partial {\rho }_{\mathbf{k}}\partial {\rho }_{-\mathbf{k}}}]-\nu {\text{k}}_z\frac{\partial {\rho }_{\mathbf{k}}}{\partial {\text{k}}_x}\frac{\partial }{\partial {\rho }_{\mathbf{k}}}\}\hat{P}=0,\end{array}$$

$$1-\left(1+\frac{\overline{\rho }}{kT}w(\mathbf{q})\right)S(\mathbf{q})+\frac{\nu }{2}\frac{q_z}{q^2}\frac{\partial S(\mathbf{q})}{\partial q_x}=0.$$

$$g(\mathbf{r})=1+\frac{1}{(2\pi {)}^3\overline{\rho }}{\int }_{\mathbf{q}\ne 0}d\mathbf{q}(S(\mathbf{q})-1)e^{i\mathbf{q}\cdot \mathbf{r}}$$

$$1-(1+\lambda (\mathbf{q}))S(\mathbf{q})+\frac{\nu }{2}\frac{q_z}{q^2}\frac{\partial S(\mathbf{q})}{\partial q_x}=0.$$

$$g(\mathbf{r})=0\text{for }r<R.$$

$$\begin{array}{rl}& \lambda (\mathbf{q})=\sum _i{\lambda }_i(\mathbf{q}){\nu }^i,\\ & S(\mathbf{q})=\sum _i{S}_i(\mathbf{q}){\nu }^i,\end{array}$$

$$-(1+{\lambda }_{PY}(q))S_1(\mathbf{q})-S_{PY}(q){\lambda }_1(\mathbf{q})+\frac{1}{2}\frac{q_z{q}_x}{q^3}\frac{d{S}_{PY}(q)}{dq}=0.$$

$$\begin{array}{rl}& {\lambda }_1(\mathbf{q})=\frac{q_x{q}_z}{q^2}{f}_1(q),\\ & S_1(\mathbf{q})=\frac{q_x{q}_z}{q^2}{h}_1(q).\end{array}$$

$$h_1(q)=-S_{PY}^2(q)f_1(q)+\frac{1}{4q}\frac{d{S}_{PY}^2(q)}{dq}.$$

$${\int }_0^{+\infty }{q}^2\left[\frac{1}{4q}\frac{d{S}_{PY}^2(q)}{dq}-S_{PY}^2(q)f_1(q)\right]j_2(qr)dq=0$$

$${\nabla }^2{\hat{\lambda }}_1(\mathbf{r})=\frac{{\partial }^2}{\partial x\partial z}{\hat{f}}_1(r).$$

$${\hat{\lambda }}_1(\mathbf{r})=-\frac{{\partial }^2}{\partial x\partial z}\varphi (r),$$

$$\varphi (r)=-\frac{1}{4\pi }\int d{\mathbf{r}}^{\prime }\frac{{\hat{f}}_1(r^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}=\frac{Q}{r}.$$

$${\int }_0^R{r}^2{\hat{f}}_1(r)dr=0.$$

$$\hat{f}_{1}(r)=\left\{\begin{array}[]{ll}\sum_{j=0}^{n}a_{j}r^{j}&\textrm{for }r\leq R\vspace{0.2cm}\\ 0&\textrm{elsewhere},\end{array}\right.$$

Full text

Systems of interacting extended objects subject to an imposed macroscopic flow have been attracting an ever increasing amount of interest over the last four decades. Examples of such systems are emulsions, suspensions, colloids and foams. The imposed flow can be achieved by immersing the system in a sheared fluid or by shearing it directly, for example between moving plates. This kind of system has diverse applications in material fabrication, biological systems, paints, etc. and has been thoroughly investigated in several ways 05Sti .

Most of the extended objects mentioned above are deformable but some of the experiments and numerical simulations have been done on model hard sphere systems 93Koel ; 80ClAc . The present article deals with the effect of imposed flow on the correlations in a system of monodisperse hard spheres immersed in it. We are interested, in particular, in the distortion of the correlations and in shear-induced jamming due to macroscopic shear flow of the host fluid.

The distortion of the structure factor of the immersed objects in the presence of shear has been observed experimentally 80ClAc and explained theoretically 84Ron . A more interesting observation seen in simulations and experiments is shear-induced ordering that is accompanied, in some cases, by jamming of the particle system (88Ack -97Farr ).

In this article, we are going to address the problem of jamming due to shear by relating it to random close packing in a system of hard spheres in equilibrium and to the distortion of the structure factor due to the shear of the host fluid.

Let us start with an overdamped Langevin description of a system of $N$ identical particles contained in a cubic box of volume $V$ and periodic boundary conditions. The particles are immersed in a host liquid flowing with a prescribed velocity field $\mathbf{V}(\mathbf{r})$,

[TABLE]

where $\mathbf{r}_{i}$ is the radius vector of the $i$-th particle, $w(\mathbf{r})$ is the two-body potential and the Cartesian components of the noise, $\eta_{i}^{k}$, obey $\langle\eta_{i}^{k}(t)\rangle=0$ and $\langle\eta_{i}^{k}(t)\eta_{j}^{l}(t^{\prime})\rangle=2\sigma\delta_{ij}\delta_{kl}\delta(t-t^{\prime})$. The Langevin equations above can be directly turned in the Fokker-Planck equation for the distribution of all the particles in the configuration space,

[TABLE]

where the operator $O$ has the following form

[TABLE]

with $kT=\sigma/\gamma$ and

[TABLE]

The Fokker Planck equation above can be transformed into an equation governing the evolution of the distribution $\hat{P}(\{\rho_{\mathbf{q}}\})$ of the collective degrees of freedom, which are the Fourier coefficients of the density $\rho(\mathbf{r})=\sum_{i=1}^{N}\delta(\mathbf{r}-\mathbf{r}_{i})$ 02Edw .

Note that the Fourier coefficients are given in terms of the coordinates of the particles by

[TABLE]

so that $\rho_{\bm{0}}=\sqrt{N}$ is not a degree of freedom.

Since we are interested only in the steady state properties, we write down just the steady state equation $(\partial\hat{P}/\partial t=0)$,

[TABLE]

where $\overline{\rho}=N/V$, $w(\mathbf{k})$ is the Fourier transform of the two-body potential $w(\mathbf{r})$ and $\nu=\gamma C/kT$.

The equation above can be simplified by the random phase approximation (RPA), which consists of dropping the trilinear terms in the operator applied on the probability distribution in the equation above vinsc2003 . Thus,

[TABLE]

Multiplying Eq.(6) by $1/2\rho_{\mathbf{q}}\rho_{-\mathbf{q}}$ and integrating by parts over all the collective variables lead directly to a closed form equation for the structure factor $S(\mathbf{q})=\langle\rho_{\mathbf{q}}\rho_{-\mathbf{q}}\rangle$,

[TABLE]

At first sight, the above discussion seems totally irrelevant to the hard sphere system because the hard sphere potential does not have a Fourier transform. However, it has been shown SchVin02 that it is useful to view the hard sphere system as an ideal gas with the constraint that the pair distribution function

[TABLE]

vanishes for distances $r=|\mathbf{r}|$ smaller than the range of the hard sphere interaction, that is the sphere diameter $R$. The constraint is imposed by introducing a Lagrange multiplier $\lambda(\mathbf{k})$ instead of $w(\mathbf{k})\overline{\rho}/kT$ in Eq.(5) and (6).

Since the constraint holds only for $r<R$, $\lambda(\mathbf{k})$ is a Fourier transform of a function $\hat{\lambda}(\mathbf{r})$ that vanishes for $r\geq R$.

Thus, the full solution involves obtaining $S(\mathbf{q})$ as a functional of such a $\lambda(\mathbf{q})$ from the equation

[TABLE]

Then $\lambda(\mathbf{q})$ has to be determined from the condition

[TABLE]

It is easy to see that in the absence of shear the conditions described by Eq.(8), (9) (with $\nu=0$) and (10) and the condition on $\hat{\lambda}(\mathbf{r})$ are equivalent to the hard sphere Percus-Yevick equation PY .

The Percus-Yevick approximation leads to an analytic solution for $\lambda(q)$ which yields very good results for liquid densities of the system. Above those densities, though, the pressure predicted by that approximation shows no sign of termination of the liquid phase.

Fortunately, however, the pair distribution function, which is a probability distribution and therefore non-negative, vanishes first at some distance for a threshold value of the volume fraction, $\eta_{c}$. For higher volume fractions, the pair distribution function becomes negative for a range of distances, indicating the failure of the Percus-Yevick approximation at volume fractions higher than the threshold value schwartz .

Since the vanishing of the pair distribution function at some distance indicates caging, it has been suggested in the above reference that the threshold volume fraction, although approximate, has a physical meaning of jamming or random close packing. This suggestion is supported by the comparison of the threshold volume fraction for dimensions between 3 to 9 obtained from the Percus-Yevick equation schwartz and the corresponding values for random close packing (RCP) RCP or maximal random jammed (MRJ) MRJ systems.

In fact, this is one of the main motivations for the present work. It is well known that a jamming transition may occur in systems of objects subject to shear 08Lu . This jamming transition may be related, like in the absence of shear, to the vanishing of the pair distribution function. So, the study of the effect of shear on the pair distribution function is one of the main goals of the present article.

Therefore, we turn our attention back to the general case of the finite shear, described by Eq.(8)-(10) and $\hat{\lambda}(\mathbf{r})=0$ for $r>R$.

This is a formidable set of equations, much more complicated than the original hard sphere Percus-Yevick equation. We may, however, try to construct a solution as follows. Since the result for the hard sphere system without shear is known, we expect to be able to write our solution as an expansion in $\nu$ around the (known and exact) Percus-Yevick solution $(\lambda_{PY}(q),S_{PY}(q))$:

[TABLE]

where $(\lambda_{0}(\mathbf{q}),S_{0}(\mathbf{q}))\equiv(\lambda_{PY}(q),S_{PY}(q))$ with $q=|\mathbf{q}|$. This will result in $g(\mathbf{r})=\sum_{i}g_{i}(\mathbf{r})\nu^{i}$ with $g_{0}(\mathbf{r})=g_{PY}(r)$.

In the present article, we will just concentrate on the first order expansion of $\lambda(\mathbf{q})$ and $S(\mathbf{q})$. Higher orders follow the same strategy that will be presented in the following but, although conceptually straightforward, are considerably more complicated.

The zero order of the system is obviously satisfied by $(\lambda_{PY}(q),S_{PY}(q))$. From the first order in $\nu$ we obtain the following equation relating the couple $(S_{1}(\mathbf{q}),\lambda_{1}(\mathbf{q}))$,

[TABLE]

The form of the last term of the equation conveys the kind of dependence on the direction of $\mathbf{q}$, so that we can write

[TABLE]

Eq.(12) is then written for $f_{1}(q)$ and $h_{1}(q)$ and eventually provides the relation between these two quantities, which depend only on the absolute value of $\mathbf{q}$:

[TABLE]

The relation above gives $S_{1}(\mathbf{q})$ in terms of $\lambda_{1}(\mathbf{q})$, in analogy to the simple relation between $S_{PY}(q)$ and $\lambda_{PY}(q)$ in the absence of shear.

Now that we have reduced the problem to finding a function of a scalar value, we can turn to the evaluation of the pair distribution function as a functional of the unknown $f_{1}(q)$. This function will then be determined from condition (10) above.

It proves useful to express the directional dependence of the integrand on the right hand side of Eq.(8) in terms of spherical harmonics. We use the expansion of the plane wave in terms of spherical harmonics and the relation $q_{x}q_{z}/q^{2}=\sqrt{2\pi/15}\left[Y_{2}^{-1}(\theta_{\mathbf{q}},\varphi_{\mathbf{q}})-1/6Y_{2}^{1}(\theta_{\mathbf{q}},\varphi_{\mathbf{q}})\right]$. The angles $(\theta_{\mathbf{q}},\varphi_{\mathbf{q}})$ determine the direction of $\mathbf{q}$.

Eq.(8) results now in a single equation for the first order correction in $\nu$,

[TABLE]

for $r<R$, where $j_{2}(qr)$ is the spherical Bessel function of order two.

Let us denote by $\hat{f}_{1}(r)$ and $\hat{\lambda}_{1}(\mathbf{r})$ the functions whose Fourier transforms are $f_{1}(q)$ and $\lambda_{1}(\mathbf{q})$ respectively. We use a parametrization for $\hat{f}_{1}(r)$ that ensures that $\hat{\lambda}_{1}(\mathbf{r})$ vanishes for $r>R$.

The relation between $\lambda_{1}(\mathbf{q})$ and $f_{1}(q)$, following Eq. (13), implies

[TABLE]

By consequence, we obviously need to require that $\hat{f}_{1}(r)=0$ for $r>R$. But this is not the whole story: Eq.(16) can be solved to yield

[TABLE]

where the spherically symmetric potential $\phi(r)$ is written in terms of the spherically symmetric charge density $-\hat{f}_{1}(r)/4\pi$,

[TABLE]

The last equality on the right hand side of the above equation holds for $r>R$ and $Q$ is the total charge contained in a sphere of radius $R$. Thus, it follows that not only the charge density has to vanish outside the sphere, but also the total charge within the sphere has to be zero. This implies

[TABLE]

We use a parametrization for $\hat{f}_{1}(r)$ that ensures that $\lambda_{1}(\mathbf{r})$ vanishes for $r>R$,

[TABLE]

in terms of the coefficients $\{a_{j}\}$.

Eq.(15) is used now to find the coefficients $\{a_{j}\}$. When doing that, we must recall that those coefficients are constrained by the requirement $Q=0$ (Eq.(19)).

Once $f_{1}(r)$ is known, we can easily set $S_{1}(\mathbf{q})$ and $g_{1}(\mathbf{r})$ (practically, we choose $n=10$ independent unknowns $a_{j}$ and this ensures us that the root mean square of $g_{1}(\mathbf{r})$ in the interval $(0,R)$ is of the order of $10^{-4}$).

Now, we can use the solution $(\lambda_{1}(\mathbf{q}),S_{1}(\mathbf{q}))$ to solve the equation for the second order term in $\nu$ but we postpone that to future study.

In order to work with dimensionless quantities, we use the volume fraction $\eta=\pi R^{3}\overline{\rho}/6$ and we define the Péclet number $P=\nu R^{2}$. The pair distribution function is then expressed as $g(\mathbf{r})=g_{PY}(r)+P/R^{2}\,g_{1}(\mathbf{r})$.

Fig. 1(a) shows $g_{PY}(r)$ in the $(r_{x},r_{z})$ plane. Obviously, $g_{PY}(r)$ is just a function of $r_{x}^{2}+r_{z}^{2}$. Fig. 1(b) shows $g_{1}(\mathbf{r})$ in the same plane: it has, in addition, also an angular dependence. In both plots, a dashed line indicates the direction of $\mathbf{r}$ determined by the angles $(\theta_{\mathbf{r}},\varphi_{\mathbf{r}})=\left(3\pi/4,0\right)$.

Fig. 2 shows the behavior of $g\left(r,3\pi/4,0\right)$, that is the pair correlation function exactly in this direction, for $\eta=0.5$ and $P=1$. The choice of the spatial angle is not accidental: it corresponds to the direction where $g_{1}(\mathbf{r})$ attains its lowest value for $r$ close to the minimum point of $g_{PY}(r)$. The interest in this direction will be better explained in the following.

It is particularly interesting to notice that the positivity condition for $g(\mathbf{r})$ is satisfied for the value of the couple $(\eta,P)$ in Fig. 2.

As already said, this condition fails to be satisfied already in the case without shear for values of $\eta$ beyond a certain threshold and this event can be interpreted as occurrence of jamming in the suspension. By adding external shear, we simply expect to facilitate the phenomenon of jamming for smaller $\eta$. Practically, we want to see when $g(\mathbf{r})$ vanishes for the first time. This means, first of all, that the contribution of $g_{1}(\mathbf{r})$ must be negative and as significant as possible. Second, this has to happen close to the minimum point for $g_{PY}(r)$. As Fig. 1 shows, this happens for $(\theta_{\mathbf{r}},\varphi_{\mathbf{r}})=(3\pi/4,0)$.

Let us choose a fixed volume fraction and vary the shear. Fig. 3 shows the behavior of $g(r,3\pi/4,0)$ around its minimum point for $\eta=0.6$ and different values of the Péclet number. For small $P$ the curve is always positive but it becomes negative for increasing values of $P$, as expected.

This leads to our final result, that is a phase diagram in the plane $(P,\eta)$ as shown in Fig. 4.

The curve can be calculated for bigger values of $P$ but it shows a slight deviation from the linear behavior that has to be expected within the first order approximation.

ne In conclusion, by combining previous analytic results for two primary problems, it is possible to write a system of equations for hard sphere suspensions subject to external shear.

We propose a suitable form of the solution, as an expansion in the shear rate around the Percus-Yevick solution for non-sheared hard sphere suspension. This allows us to capture the distortion of the structure factor and its angular dependence for small shear. From this distortion we obtain the (deformed) pair distribution function and the volume fraction for which it vanishes for the first time.

Consequently, we derive the border between the fluid phase and the jammed phase in the suspension, as a function of volume fraction and shear.

The problem has still many aspects to investigate. A natural next step is the comparison with numerical and experimental results.

Moreover, a better approximation for the solution, hopefully valid for a wider range of shear rates, will be the object of future research.