Microstructure and thickening of dense suspensions under extensional and shear flows

TL;DR

This study uses numerical simulations to explore how microstructure influences the non-Newtonian rheological behavior of dense suspensions under extensional and shear flows, revealing flow-type dependence and normal stress differences.

Contribution

It provides new insights into the microstructure-flow relationship and flow-type dependence in dense suspensions under different flow conditions.

Findings

Flow-type dependence varies below and above thickening regimes.

Normal stress differences are present regardless of flow-type dependence.

Microstructure influences rheological properties significantly.

Abstract

Dense suspensions are non-Newtonian fluids which exhibit strong shear thickening and normal stress differences. Using numerical simulation of extensional and shear flows, we investigate how rheological properties are determined by the microstructure which is built under flows and by the interactions between particles. By imposing extensional and shear flows, we can assess the degree of flow-type dependence in regimes below and above thickening. Even when the flow-type dependence is hindered, nondissipative responses, such as normal stress differences, are present and characterise the non-Newtonian behaviour of dense suspensions.

Equations

$$\rho \left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}\right]=\nabla \cdot \mathbf{\sigma }\text{ with }\nabla \cdot \mathbf{u}=0,$$

$${\mathbf{F}}_{\mathrm{H}}+{\mathbf{F}}_{\mathrm{C}}+{\mathbf{F}}_{\mathrm{R}}=0\text{and}{\mathbf{T}}_{\mathrm{H}}+{\mathbf{T}}_{\mathrm{C}}=0.$$

$$\left(\begin{array}{c}{\mathbf{F}}_{\mathrm{H}}\\ {\mathbf{T}}_{\mathrm{H}}\end{array}\right)=-\text{mathsfbi}R\cdot \left(\begin{array}{c}\Delta \mathbf{U}\\ \Delta \Omega \end{array}\right)+\text{mathsfbi}{R}^{\prime }:\text{mathsfbi}{D}_N,$$

$$\bm{\sigma}=-p_{0}\mathsfbi{I}+2\eta_{0}\mathsfbi{D}+V^{-1}\biggl{(}\sum_{i}\mathsfbi{S}_{\mathrm{D}}^{(i)}+\sum_{i>j}\mathsfbi{S}_{\mathrm{P}}^{(i,j)}\biggr{)},$$

$$\text{mathsfbi}\hat{D}={\hat{\mathbf{d}}}_1{\hat{\mathbf{d}}}_1-{\hat{\mathbf{d}}}_2{\hat{\mathbf{d}}}_2,\text{mathsfbi}\hat{W}={\beta }_3({\hat{\mathbf{d}}}_2{\hat{\mathbf{d}}}_1-{\hat{\mathbf{d}}}_1{\hat{\mathbf{d}}}_2).$$

$$\bm{\sigma}(\dot{\varepsilon},\beta_{3})=-p(\dot{\varepsilon},\beta_{3})\mathsfbi{I}+\dot{\varepsilon}\bigl{[}\kappa(\dot{\varepsilon},\beta_{3})\mathsfbi{\hat{D}}+\lambda_{0}(\dot{\varepsilon},\beta_{3})\mathsfbi{\hat{E}}+\lambda_{3}(\dot{\varepsilon},\beta_{3})\mathsfbi{\hat{G}}_{3}\bigr{]},$$

$$\varphi\equiv\arctan\Bigl{[}\lambda_{3}\Big{/}\Bigl{(}\kappa+\sqrt{\kappa^{2}+\lambda_{3}^{2}}\Bigr{)}\Bigr{]}.$$

$$\eta (2\dot{\epsilon })=\kappa (\dot{\epsilon },1)/2,N_1(2\dot{\epsilon })=-2\dot{\epsilon }{\lambda }_3(\dot{\epsilon },1),\text{and}{N}_2(2\dot{\epsilon })=\dot{\epsilon }[{\lambda }_3(\dot{\epsilon },1)-(3/2){\lambda }_0(\dot{\epsilon },1)].$$

$$\text{mbox}St\equiv \frac{2{\rho }_{\mathrm{p}}{a}^2\dot{\epsilon }}{{\eta }_0}=\frac{{\rho }_{\mathrm{p}}{F}_{\mathrm{R}}(0)}{6\pi {\eta }_0^2}\frac{\dot{\epsilon }}{{\dot{\epsilon }}_0},$$

Full text

1 Introduction

Suspensions, namely mixtures of solid particles and a viscous liquid, can be considered as an incompressible fluid as long as the volume fraction $\phi$ of solid particles is less than a certain value, the jamming point, above which a solid-like behaviour is observed. The behaviour of suspensions is not usually captured by simple Newtonian models. As primary example of non-Newtonian effect, the viscosity can vary with the shear rate, exhibiting shear thinning and shear thickening (Laun, 1984; Barnes, 1989; Bender & Wagner, 1996; Guy et al., 2015). Moreover, nonvanishing normal stress differences $N_{1}$ and $N_{2}$, another hallmark of non-Newtonian behaviour, are often observed (Laun, 1994; Lootens et al., 2005; Lee et al., 2006; Couturier et al., 2011; Dbouk et al., 2013; Cwalina & Wagner, 2014). Discontinuous shear thickening is a particularly intriguing phenomenon of dense suspensions and the underlying mechanism raised a significant debate (Brady & Bossis, 1985; Hoffman, 1998; Melrose & Ball, 2004; Fall et al., 2008; Brown & Jaeger, 2009). Analysing the rheology of suspensions is a difficult task since forces of various nature act among particles and the system lives mostly far from thermodynamic equilibrium. Particle simulations have been used to explore the microstructure emerging among particles in various flows and to estimate the importance of different interactions. Several particle simulations recently succeeded in reproducing shear thickening by taking into account direct contact forces (Fernandez et al., 2013; Heussinger, 2013; Seto et al., 2013). These works support the “stress-induced friction” scenario (Wyart & Cates, 2014; Mari et al., 2014) and the contribution of contact forces was also confirmed in experiments (Lin et al., 2015; Clavaud et al., 2017). Thus, the particle-scale mechanism of shear thickening is, to a great extent, understood.

However, particle-scale simulations are not capable of reproducing engineering-scale flows of dense suspensions due to the practical limits imposed on the system size by computational tractability. For this reason, it is important to develop effective continuum models through the design of suitable non-Newtonian constitutive relations. Besides laboratory experiments, particle-scale simulations are an important source of indications for the development of such models. A complete model should describe the fluid response under any flow condition (Miller et al., 2009), not only in the simple shear flows in which most experimental and computational data are retrieved. Indeed, the response of non-Newtonian fluids can depend on the type of flow, as exemplified by the observations of shear thinning and extensional thickening in some viscoelastic fluid. Particularly important is the class of extensional flows of suspensions, for which few rheological characterisations are available (Dai & Tanner, 2017) and the sole computational investigation of which the authors are aware was performed by Sami (1996), who studied semidilute Brownian suspensions. (We note that in his analysis flow-type dependence was not evidenced.) A related computational method to treat hydrodynamic interactions in diluted suspensions was introduced by Ahamadi & Harlen (2008). For important developments regarding emulsions of deformable droplets, we refer the reader to the work of Zinchenko & Davis (2015).

To study the material response, we simulate motions of particles in the bulk region under prescribed flow conditions. As usual, periodic boundary conditions are employed to minimise finite-size effects. The Lees–Edwards boundary conditions (Lees & Edwards, 1972) are commonly used to impose simple shear flows in many contexts, including suspension rheology (Bossis & Brady, 1984; Mari et al., 2014). In this work, we also apply the Kraynik–Reinelt boundary conditions (Kraynik & Reinelt, 1992; Todd & Daivis, 1998), originally devised to impose planar extensional flows in nonequilibrium molecular dynamics simulations. With these we can provide a first assessment of the flow-type dependence of the response in dense suspensions.

In §2.1 and §2.2 we describe our simulation technique which operates in the inertialess approximation. To compare consistently the results under different flow conditions, we employ the rheometric framework introduced by Giusteri & Seto (2017) (summarised in §2.3) which defines, for the case of planar flows, a dissipative response function, $\kappa$, and two non-dissipative response functions, $\lambda_{0}$ and $\lambda_{3}$. Those are defined for any flow type (simple shear, extensional, and mixed flows) and offer a unified description of the material response. The results of our analysis, discussed in §3, highlight the presence of flow-type dependence in the microstructure and in the non-Newtonian effects observed for dense suspensions.

2 Methods

2.1 Bulk rheology with periodic boundary conditions

Non-Newtonian incompressible fluids obey the differential equations

[TABLE]

where $\bm{u}$ is the velocity field and $\rho$ the density. To close the system of equations, the stress tensor $\bm{\sigma}$ must be given in terms of the velocity gradient through a constitutive prescription. The local value of the stress tensor describes the material response and is determined by the local history of deformation. To investigate such response, we consider small volume elements in which the velocity gradient $\nabla\bm{u}$ is approximately uniform. By simulating motions of particles in the volume element with fixed $\nabla\bm{u}$, we can find the typical stress for a certain deformation history.

Time-dependent periodic boundary conditions allow to impose $\nabla\bm{u}$ and effectively simulate the bulk behaviour. Since we will consider planar flows in a 3D geometry, we can describe our methods considering the 2D projections of the computational cells. The cell frame vectors $\bm{l}_{1}(t)$ and $\bm{l}_{2}(t)$ (see Figure 1) are prescribed to follow the velocity field $\bm{u}=\nabla\bm{u}\cdot\bm{r}$ and periodic images of a particle at $\bm{r}$ are given by $\bm{r}^{\prime}=\bm{r}+i\bm{l}_{1}(t)+j\bm{l}_{2}(t)$ with ($i,j=\pm 1,\pm 2,\dotsc$). For simple shear flows ($\nabla\bm{u}=\dot{\gamma}\bm{e}_{y}\bm{e}_{x}$), this is equivalent to the Lees–Edwards boundary conditions. The initial periodic cells are rectangles in the flow plane (blue in Figure 1 (a)). A simple shear flow deforms the cells to parallelogram shapes (red in Figure 1 (a)). To avoid significantly deformed periodic cells, the initial rectangular cells can be recovered as shown in Figure 1 (a). To impose planar extensional flows ($\nabla\bm{u}=\dot{\varepsilon}\bm{e}_{x}\bm{e}_{x}-\dot{\varepsilon}\bm{e}_{y}\bm{e}_{y}$) for long times, we employ the Kraynik–Reinelt periodic boundary conditions (Kraynik & Reinelt, 1992; Todd & Daivis, 1998). If the initial master cell is a regular square oriented at a certain angle $\theta^{\ast}$ from the extension axis ($x$-axis), the deformed parallelogram cell after a certain strain $\varepsilon_{\mathrm{p}}$ can be remapped to the initial regular shape as shown in Figure 1 (b).

2.2 Inertialess particle dynamics for suspensions

We numerically evaluate the stress tensor $\bm{\sigma}$ by using particle simulations with deforming periodic cells. Our simulation is analogous to rate-controlled rheological measurements in the sense that time-averaged stress responses $\langle\bm{\sigma}\rangle$ are evaluated for imposed velocity gradients $\nabla\bm{u}$.

We consider non-Brownian, density matched, and dense suspensions. Suspended particles interact with each other in several ways. As discussed in Mari et al. (2014), we take into account contact forces $\bm{F}_{\mathrm{C}}$ (and torques $\bm{T}_{\mathrm{C}}$) and stabilising repulsive forces $\bm{F}_{\mathrm{R}}$, besides hydrodynamic interactions $\bm{F}_{\mathrm{H}}$ and $\bm{T}_{\mathrm{H}}$. Since the inertia of sufficiently small particles is negligible in comparison to the hydrodynamic drag forces, the particles obey the quasi-static equations of motion

[TABLE]

Here, forces $\bm{F}$ and torques $\bm{T}$ represent the set of forces and torques for $N$ particles. Flows around microscale particles are dominated by viscous dissipation and the inertia of the fluid is negligible, so they are described by the Stokes equations. The imposed velocity gradient $\nabla\bm{u}$ gives the background flow via the velocity $\bm{u}(\bm{r})$, vorticity $\bm{\omega}\equiv\nabla\times\bm{u}$, and rate of deformation tensor $\mathsfbi{D}$ such that $\nabla\bm{u}\cdot\bm{r}=\bm{u}(\bm{r})=\mathsfbi{D}\cdot\bm{r}+(\bm{\omega}/2)\times\bm{r}.$ In this case, the hydrodynamic interactions can be expressed as the sum of linear resistances to the relative velocities $\varDelta\bm{U}^{(i)}\equiv\bm{U}^{(i)}-\bm{u}(\bm{r}^{(i)})$, angular velocities $\varDelta\bm{\varOmega}^{(i)}\equiv\bm{\varOmega}^{(i)}-\bm{\omega}/2$, for $i=1,\dotsc,N$, and imposed deformation $\mathsfbi{D}$ via

[TABLE]

where $\varDelta\bm{U}$ and $\varDelta\bm{\varOmega}$ represent the set of relative velocities for $N$ particles, $\mathsfbi{D}_{N}$ is block-diagonal with $N$ copies of $\mathsfbi{D}$, and $\mathsfbi{R}$ and $\mathsfbi{R}^{\prime}$ are resistance matrices which can, in principle, be derived from the Stokes equations once the particle configurations are given. In dense suspensions, the long-range hydrodynamic interactions are screened by crowds of particles. Therefore, we may approximately construct the resistance matrices by including only the contributions of Stokes drag and lubrication forces.

In real suspensions, the lubrication singularity in $\bm{F}_{\mathrm{H}}$ is absent due to factors such as the surface roughness of particles—direct contacts are not forbidden. Hence, we include the contact interactions $\bm{F}_{\mathrm{C}}$ and $\bm{T}_{\mathrm{C}}$ in (2). The contact forces between solid particles depend on the nature of the particle surfaces. This is effectively encoded in the friction coefficient $\mu$ that enters a simple friction model. By denoting with $F_{\mathrm{C}}^{n}$ and $F_{\mathrm{C}}^{t}$ normal and tangential forces, respectively, we prevent sliding if $F_{\mathrm{C}}^{t}\leq\mu F_{\mathrm{C}}^{n}$. The normal force depends on the overlap between particles through an effective elastic constant and the tangential force depends on the sliding displacement in a similar way. The details of the employed model are given in Mari et al. (2014).

The presence of the stabilising repulsive force $\bm{F}_{\mathrm{R}}$ in (2) generates the rate dependence of rheological properties in such suspensions. Indeed, while reaching the same strain, $\bm{F}_{\mathrm{R}}$ can work more to prevent particle contacts under lower deformation rates, but less under higher rates. As a result, the number of contacts depends on the rate of the imposed flow. In colloidal suspensions, Brownian forces may play a similar role, as discussed by Mari et al. (2015a).

The bulk stress tensor is obtained as

[TABLE]

where $\eta_{0}$ is the viscosity of the solvent, $\mathsfbi{S}_{\mathrm{D}}^{(i)}$ is the stresslet on particle $i$ due to $\mathsfbi{D}$, and $\mathsfbi{S}_{\mathrm{P}}^{(i,j)}$ is the stresslet due to non-hydrodynamic interparticle forces between particle $i$ and $j$ (Mari et al., 2015b). Note that, since the hydrostatic pressure $p_{0}$ is arbitrary, we set $p_{0}=0$. However, the last term in (4) is not traceless and thus contributes to the total pressure $p\equiv-(1/3)\operatorname{Tr}{\bm{\sigma}}$.

2.3 General response functions for steady flows of non-Newtonian fluids

The stress $\bm{\sigma}$ is a tensorial quantity and we need a procedure to extract from it the relevant information in terms of scalar quantities. We are interested in comparing the material response under different types of imposed flow conditions. For this reason, we need a framework in which it is possible to identify the dependence on the flow type of each independent non-Newtonian effect. To this end, we use the framework introduced by Giusteri & Seto (2017), in which the characteristic rate of the imposed flow is defined independently of the flow type and a complete set of response functions is given. These functions generalise to any flow type standard quantities such as viscosity and normal stress differences.

The velocity gradient $\nabla\bm{u}$ is decomposed into symmetric and antisymmetric parts as $\nabla\bm{u}=\mathsfbi{D}+\mathsfbi{W}$. In the planar case, with $\mathsfbi{D}\neq\mathsfbi{0}$, we denote by $\dot{\varepsilon}$ the largest eigenvalue of $\mathsfbi{D}$ and express $\mathsfbi{\hat{D}}\equiv\mathsfbi{D}/\dot{\varepsilon}$ and $\mathsfbi{\hat{W}}\equiv\mathsfbi{W}/\dot{\varepsilon}$ on the basis of the eigenvectors $\hat{\bm{d}}_{1}$ and $\hat{\bm{d}}_{2}$ of $\mathsfbi{D}$ (corresponding to the eigenvalues $\dot{\varepsilon}$ and $-\dot{\varepsilon}$) as follows:

[TABLE]

The non-vanishing and positive rate $\dot{\varepsilon}>0$ is used to set the time scale of deformation in any flow type. With this definition, the standard rate $\dot{\gamma}$ for simple shear corresponds to the value $2\dot{\varepsilon}$. The vorticity $\omega_{z}$ is represented by the dimensionless parameter $\beta_{3}$ through $\omega_{z}=2\dot{\varepsilon}\beta_{3}$. Note that planar extensional flows are characterised by $\beta_{3}=0$, and simple shear flows by $\beta_{3}=1$.

A general representation of the stress tensor in planar flows is then given by

[TABLE]

where $\mathsfbi{\hat{E}}\equiv-(1/2)(\hat{\bm{d}}_{1}\hat{\bm{d}}_{1}+\hat{\bm{d}}_{2}\hat{\bm{d}}_{2})+\hat{\bm{d}}_{3}\hat{\bm{d}}_{3}$, $\hat{\bm{d}}_{3}$ is the eigenvector of $\mathsfbi{D}$ orthogonal to the flow plane, and $\mathsfbi{\hat{G}}_{3}\equiv\hat{\bm{d}}_{1}\hat{\bm{d}}_{2}+\hat{\bm{d}}_{2}\hat{\bm{d}}_{1}$ is introduced to complete an orthogonal basis for the space of symmetric tensors for planar flows. The functional dependence of $\kappa$, $\lambda_{0}$, and $\lambda_{3}$ on the two kinematical parameters $\dot{\varepsilon}$ and $\beta_{3}$ needs to be determined to characterise the response in generic flows. We remark that the response functions $\kappa$, $\lambda_{0}$, and $\lambda_{3}$ can depend on any other quantity that characterise the system. For instance, in §3 we will also show their dependence on the volume fraction $\phi$. The function $\kappa$ is the only one to carry information about dissipation. We therefore refer to it as the dissipative response function, a generalised viscosity. The functions $\lambda_{0}$ and $\lambda_{3}$ carry information about non-dissipative responses and we call them non-dissipative response functions. The presence of a nonvanishing $\lambda_{0}$ leaves the eigenvectors of the stress $\bm{\sigma}$ aligned with those of $\mathsfbi{D}$, as happens in Newtonian fluids, but gives a contribution to the stress in the form of a modified pressure which is isotropic in the flow plane but different in the direction normal to the flow plane. On the other hand, a nonvanishing $\lambda_{3}$ corresponds to a rotation of the eigenvectors of $\bm{\sigma}$ in the flow plane with respect to those of $\mathsfbi{D}$, determining the reorientation angle

[TABLE]

For the sake of comparison, the shear viscosity $\eta$ and normal stress differences $N_{1}$ and $N_{2}$, defined for simple shear flows with $\beta_{3}=1$ as functions of $\dot{\gamma}=2\dot{\varepsilon}$, are given by

[TABLE]

Moreover, the extensional viscosity, defined for extensional flows with $\beta_{3}=0$, is given by $\eta_{\mathrm{E}}(\dot{\varepsilon})=2\kappa(\dot{\varepsilon},0)$. So that the Trouton ratio $\eta_{\mathrm{E}}(\dot{\varepsilon})/\eta(2\dot{\varepsilon})$ equals $4$ only if $\kappa(\dot{\varepsilon},0)=\kappa(\dot{\varepsilon},1)$.

We want to stress that, to arrive at (6), no a priori assumption is made on the list of quantities on which the material functions can depend. Hence, the stress tensor for any non-Newtonian fluid model under steady flow conditions can be expressed in the form (6). For example, since $\mathsfbi{D}^{2}$ in planar flows is a linear combination of $\mathsfbi{I}$ and $\mathsfbi{\hat{E}}$, the class of Reiner–Rivlin fluids corresponds to choosing $\lambda_{3}=0$, and assuming $\kappa$ and $\lambda_{0}$ independent of $\beta_{3}$. Similarly, second-order fluids under steady shear flows would produce constant values of $\kappa$ and $\lambda_{0}$, and entail $\lambda_{3}\propto\beta_{3}\dot{\varepsilon}$, since, under such flows, $\mathsfbi{W}\cdot\mathsfbi{D}-\mathsfbi{D}\cdot\mathsfbi{W}=\beta_{3}\dot{\varepsilon}^{2}\mathsfbi{\hat{G}}_{3}$. A detailed discussion of the representation (6) and its relation to fluid models are given in Giusteri & Seto (2017).

3 Results

To obtain the numerical results, we mainly performed 50 independent 3D simulations with $2000$ particles. The periodic cells are initially cuboids with ratio $5:5:1$. We also performed some simulations with 4000 particles using double-sized cells ($5:5:2$) to confirm the absence of significant finite-size effects (data are not shown). Regarding the friction coefficient, we set $\mu=1$ since it is the value that, in a previous paper (Mari et al., 2015a), was found to give good agreement with the experimental data by Cwalina & Wagner (2014). It is worth mentioning that Tanner & Dai (2016) showed that $\mu=0.5$ gives a better quantitative agreement with different experimental data. Nevertheless, such a fine tuning of $\mu$ is not necessary for our qualitative analysis.

The short-range repulsive force is given by $|\bm{F}_{\mathrm{R}}|=F_{\mathrm{R}}(0)\exp[-(r-2a)/(0.02a)]$, with $a$ the particle radius. A reference rate is set as $\dot{\varepsilon}_{0}\equiv F_{\mathrm{R}}(0)/(12\pi\eta_{0}a^{2})$. To estimate the importance of the inertial effects, we can use the Stokes number given by

[TABLE]

with $\rho_{\mathrm{p}}$ the particle density. Hence, inertial effects can be neglected if $\mathrm{St}\ll 1$, that is when the ratio $\dot{\varepsilon}/\dot{\varepsilon}_{0}$ is much smaller than $6\pi\eta_{0}^{2}/\rho_{\mathrm{p}}F_{\mathrm{R}}(0)$. The preceding threshold determines, for each specific system, the region in which the rheology curves obtained with our simulations can be expected to be in agreement with real data.

3.1 Dissipative response function $\kappa$

For the case of monodisperse suspensions, the dissipative response function $\kappa$ significantly increases with the rate $\dot{\varepsilon}$ both in simple shear and extensional flows (Figure 2). Not only shear thickening but also extensional thickening occurs. However, below thickening, there is a clear flow-type dependence. The value of $\kappa$ in extensional flow is much higher than the one obtained in simple shear flow (the Trouton ratio is much larger than $4$). On the other hand, above thickening, the values of $\kappa$ in extensional and shear flows are almost indistinguishable (the Trouton ratio is very close to $4$) and the flow-type dependence is hindered.

The significant discrepancy observed below thickening is due to shear-induced ordering, which can occur only in simple shear flows, as we confirm by analysing the pair distribution functions in §3.4. Since streamlines of a simple shear flow are straight and parallel to each other, particles tend to be arranged in chain-like structures along the flow direction. We observe a gradual decrease of $\kappa$ over time (strain thinning) in simple shear flows, which indicates the growth of the ordered structure. It should be noted that the shear-induced ordering is enhanced by the periodic boundary conditions, since linear chains may connect with their own periodic images. By contrast, the streamlines of extensional flows are never parallel to each other. Therefore, there is no obvious ordered structure compatible with extensional flows. Indeed, we neither observe strain thinning nor any ordered microstructure in the extensional flow simulation.

In the thickened regime, frictional contact forces are constantly activated. Such contact forces are so strong that particles are easily prevented from following the background flow, thus ordered structures cannot be developed. As long as the disordered structure is maintained under simple shear flows, the value of $\kappa$ remains very close to that observed in extensional flows.

The shear-induced ordering can be hindered by mixing particles with different sizes. To see this effect, we consider two types of bidisperse suspensions with different size ratios: $a_{2}/a_{1}=1.2$ and $a_{2}/a_{1}=1.4$ (named “weak” and “strong”, respectively). Two populations occupy the same volume fractions, i.e., $\phi_{1}=\phi_{2}=\phi/2$. In the weakly bidisperse suspensions (Figure 3 (a)), although the differences clearly become smaller, some flow-type dependence can still be seen, especially for $\phi=0.5$. In the strongly bidisperse suspensions (Figure 3 (b)), we no longer see a noticeable flow-type dependence—Trouton ratios are always close to $4$.

3.2 Pressure and anisotropic response

The total stress tensor $\bm{\sigma}$ is usually split into two parts: isotropic pressure term and traceless extra-stress term. Though only the extra-stress term determines the flows of incompressible fluids, the pressure $p$ is also a part of the material response. As seen in Figure 4 (a) for monodisperse suspensions in extensional flows, the pressure term $p$ varies in a similar way as $\dot{\varepsilon}\kappa$; the ratio $\dot{\varepsilon}\kappa/p$ remains of the order of unity even when $\kappa$ significantly increases by thickening. In our simulation, the volume of the periodic cells is fixed, therefore the system can never dilate. However, such increase of $p$ with $\dot{\varepsilon}\kappa$ suggests that extensional thickening (and shear thickening) of suspensions is a phenomenon related to that of dilatancy in granular materials.

The pressure term $p$ contributes isotropically to $\bm{\sigma}$ by definition. However, there is another contribution to the stress $\bm{\sigma}$ sharing the same origin. The non-dissipative response associated with dilatancy can be anisotropic and activate the response function $\lambda_{0}$. The dimensionless ratio $\dot{\varepsilon}\lambda_{0}/p$ represents such anisotropy. Its positive values reported in Figure 4 (b) indicate that the in-plane pressure is higher than the out-of-plane pressure. However, this anisotropy is not very strong, as it would be if the pressure dilatancy were only present in the flow plane.

3.3 Reorientation angle of stress eigenvectors

Besides the ordering in simple shear flow, we can see some flow-type dependence in the reorientation angle $\varphi$, defined in (7). In extensional flows, the principal axes of the stress tensor $\bm{\sigma}$ must be parallel to the eigenvectors of $\mathsfbi{D}$ due to symmetry considerations. Indeed, the reorientation angle $\varphi$ fluctuates around zero in those simulations. In simple shear flows, the shear-induced ordering is accompanied by large negative values of $\varphi$ (Figure 5 (a)). On the other hand, in the disordered states above thickening (Figure 5 (b)) and with strong bidispersity (Figure 5 (c)), $\varphi$ is always rather small but non-zero. In our inertialess simulation, this finite flow-type dependence indicates some characteristic microstructure (see §3.4) due to the presence of vorticity in simple shear, which is absent in extensional flows.

It is worth commenting on the dependence of the angle $\varphi$ on the volume fraction $\phi$ (Figure 5). In the thickened regime, corresponding to higher values of $\dot{\varepsilon}$, the angle $\varphi$ is always positive for $\phi=0.5$. The values of $\varphi$ become smaller and can take slightly negative values as $\phi$ increases. This behaviour is consistent with some experimental measurements of $N_{1}$, which is proportional to $-\lambda_{3}$. When the volume fraction is not very high, negative values have been observed for $N_{1}$ (Lee et al., 2006; Cwalina & Wagner, 2014), corresponding to positive $\varphi$, while the sign of $N_{1}$ turns positive (negative $\varphi$) at higher volume fractions (Lootens et al., 2005; Dbouk et al., 2013).

3.4 Microstructure

As discussed in the modelling section §2.2, it is reasonable to neglect particle and fluid inertia in the particle-scale dynamics. The response of such inertialess material elements to an imposed flow essentially depends on the microstructure built by the particles during the flow (Morris, 2009). To measure the correlation of particle positions, we evaluate the pair distribution function $g(\bm{r})\equiv P_{1|1}(\bm{r}|\bm{0})/n$, where $n$ is the average number density of particles and $P_{1|1}(\bm{r}|\bm{0})$ is the conditional probability of finding a particle at $\bm{r}$ with the condition that another particle is at the origin $\bm{0}$. Figure 6 (a) and (b) show $g(\bm{r})$ in the flow-plane slice $|z|<0.1a$. We also consider angular distributions $g_{\mathrm{c}}(\theta)$ for contacting (and nearly contacting) particles such that $|\bm{r}|<2.02a$. The angle $\theta$ is measured from the $\hat{\bm{d}}_{1}$ axis.

As seen in Figure 6 (a), a stripe-patterned correlation $g(\bm{r})$ appears for the monodisperse suspensions in simple shear flow below thickening ($\dot{\varepsilon}/\dot{\varepsilon}_{0}=0.01$). The periodic peaks and striped correlation indicate the formation of chain-like structures by the particles. Once such chain-like structure is formed, particle interactions are rather weak, which leads to significantly low values of $\kappa$ as seen in Figure 2. The microstructure is totally different above thickening ($\dot{\varepsilon}/\dot{\varepsilon}_{0}=0.02$). The long-range correlation is no longer seen. The correlation pattern indicates some disordered anisotropic microstructure. In Figure 6 (c), the angular contact distribution $g_{\mathrm{c}}(\theta)$ clearly shows that the number of contacting particles remarkably increases for $\dot{\varepsilon}/\dot{\varepsilon}_{0}\geq 0.02$. This observation is consistent with the idea that shear thickening is caused by the development of the contact network (Seto et al., 2013).

These results can be directly compared with those for the extensional flow simulation. As seen in Figure 6 (b), even below thickening ($\dot{\varepsilon}/\dot{\varepsilon}_{0}=0.01$), there is no long-range correlation in $g(\bm{r})$. The distribution pattern has horizontal and vertical mirror symmetries and no vorticity skews the correlation in the extensional flow. The distribution pattern does not change much above thickening ($\dot{\varepsilon}/\dot{\varepsilon}_{0}=0.02$). But a clear difference is present in the angular contact distribution $g_{\mathrm{c}}(\theta)$ (Figure 6 (d)). A flame-shaped distribution transforms into a fan-shaped distribution at the extensional thickening transition. Thus, just below the transition, we can find contacting particles only around the directions of the compression axis; nevertheless, the width of the flame shape indicates that, differently from what we observed under shear, the contact chains do not correspond to stable ordered chains of particles. Rather, they are constantly rebuilt among new neighbouring particles. Such contact chains which are roughly parallel and oriented along the compression axis do not contribute to the viscosity significantly (Figure 2). By contrast, above the thickening transition, contacting particles can be found in all directions, even in the directions of the extension axis ($\theta=0$ and $\pi$); such distribution suggests an anisotropic network structure for the pattern of contacts, which enhances the viscosity. Thus, we can describe the essence of extensional thickening as a contact-chain to contact-network transition.

The fact that such a transition occurs in extensional flows without significant change of the long-range correlation (always absent) indicates that also in simple shear flows the main responsible for thickening is the contact-chain to contact-network transition. Indeed, while we observe a concurrent order-disorder structural transition in monodisperse suspensions under shear, this is not present in strongly bidisperse suspensions, which nevertheless display a strong thickening behaviour.

4 Conclusions

We numerically explored the non-Newtonian character of dense suspensions, which has a different origin from that of viscoelastic fluids. This character is manifested in three main aspects: rate dependence, non-dissipative responses, and flow-type dependence. Analysing thickening in both extensional and simple shear flows, we were able to confirm that the contact-chain to contact-network transition is its main cause. Non-dissipative responses, such as normal stress differences, are present in any flow regime. Flow-type dependence is evident in monodisperse suspensions below thickening, where ordering occurs under simple shear. Preventing ordering through thickening or polydispersity hinders (but does not cancel) the flow-type dependence.

Acknowledgements

The authors acknowledge support from the Okinawa Institute of Science and Technology Graduate University. The research of Ryohei Seto is partially supported by JSPS KAKENHI Grant Number JP17K05618.