Simulation of dense non-Brownian suspensions with the lattice Boltzmann method: Shear jammed and fragile states

TL;DR

This study uses the lattice Boltzmann method to simulate dense non-Brownian suspensions, revealing shear-jammed and fragile states, and providing insights into their mechanical responses under different shear protocols.

Contribution

First simulation of shear-jammed and fragile states in dense suspensions incorporating hydrodynamics and frictional contacts using the lattice Boltzmann method.

Findings

Reproduced discontinuous shear thickening in 3D systems.

Identified shear-jammed states with finite storage modulus.

Observed fragile states with mixed fluid-like and solid-like responses.

Abstract

Dense non-Brownian suspensions including both the hydrodynamic interactions and the frictional contacts between particles are numerically studied under simple and oscillatory shears in terms of the lattice Boltzmann method. We successfully reproduce the discontinuous shear thickening (DST) under a simple shear for bulk three-dimensional systems. For our simulation of an oscillatory shear in a quasi-two-dimensional system, we measure the mechanical response when we reduce the strain amplitude after the initial oscillations with a larger strain amplitude. Here, we find the existence of the shear-jammed state under this protocol in which the storage modulus $G^{\prime}$ is only finite for high initial strain amplitude $\gamma_0^{I}$. We also find the existence of the fragile state in which both fluid-like and solid-like responses can be detected for an identical area fraction and an initial strain amplitude $\gamma_0^{I}$ depending on the initial phase $\Theta$ (or the asymmetricity of the applied strain) of the oscillatory shear. We also observe the DST-like behavior under the oscillatory shear in the fragile state. Moreover, we find that the stress anisotropy becomes large in the fragile state. Finally, we confirm that the stress formula based on the angular distribution of the contact force recovers the contact contributions to the stress tensors for both simple and oscillatory shears with large strains.

Equations

$$\frac{d}{dt}\left(\begin{array}{c}m\mathbf{U}\\ I\mathbf{\Omega }\end{array}\right)=\left(\begin{array}{c}{\mathbf{F}}^h+{\mathbf{F}}^c+{\mathbf{F}}^r\\ {\mathbf{T}}^h+{\mathbf{T}}^c\end{array}\right),$$

$${\mathbf{F}}_i^c=\sum _{i\ne j}({\mathbf{F}}_{ij}^{\text{nor}}+{\mathbf{F}}_{ij}^{\text{tan}}),$$

$${\mathbf{T}}_i^c=\sum _{i\ne j}{a}_i{\mathbf{n}}_{ij}\times {\mathbf{F}}_{ij}^{\text{tan}},$$

$${\sigma }_{\alpha \beta }^c=-\frac{1}{2V}\sum _i\sum _{j\ne i}(r_{ij,\alpha }{F}_{ij,\beta }^c+r_{ij,\beta }{F}_{ij,\alpha }^c),$$

$${\mathbf{F}}_{ij}^r=F_0\exp (-h/\lambda ){\mathbf{n}}_{ij},$$

$${\sigma }_{\alpha \beta }^r=-\frac{1}{2V}\sum _i\sum _{j\ne i}(r_{ij,\alpha }{F}_{ij,\beta }^r+r_{ij,\beta }{F}_{ij,\alpha }^r).$$

$${\sigma }_{\alpha \beta }={\sigma }_{\alpha \beta }^h+{\sigma }_{\alpha \beta }^c+{\sigma }_{\alpha \beta }^r.$$

$${\zeta }_N(\theta )$$

$${\zeta }_T(\theta )$$

$${\xi }_N(\theta )$$

$${\xi }_T(\theta )$$

$${\int }_0^{\pi }\rho (\theta )d\theta$$

$${\hat{\sigma }}_{\alpha \beta }^{(3)}=\frac{3{\varphi }_{3}Z⟨F^{\text{nor}}⟩}{2\pi a^2}{\int }_0^{\pi }\left[{\zeta }_N(\theta )n_{\theta ,\alpha }-{\zeta }_T(\theta )t_{\theta ,\alpha }\right]n_{\theta ,\beta }d\theta ,$$

$$\gamma (t)={\gamma }_0[\cos \Theta -\cos (\omega t+\Theta )],$$

$$G^{\prime }$$

$$G^{\prime \prime }$$

$${\hat{\sigma }}_{\alpha \beta }=\frac{\varphi Z⟨F^{\text{nor}}⟩}{\pi a}{\int }_0^{\pi }\left[{\zeta }_N(\theta )n_{\theta ,\alpha }-{\zeta }_T(\theta )t_{\theta ,\alpha }\right]n_{\theta ,\beta }d\theta .$$

$$n(\mathbf{r}+{\mathbf{c}}_{\mathbf{q}}\Delta t,{\mathbf{c}}_{\mathbf{q}},t+\Delta t)=n(\mathbf{r},{\mathbf{c}}_{\mathbf{q}},t)+{\Delta }_{\mathbf{q}}[n(\mathbf{r},{\mathbf{c}}_{\mathbf{q}},t)],$$

$$\rho =\sum _{\mathbf{q}}n({\mathbf{c}}_{\mathbf{q}}),\mathbf{j}=\sum _{\mathbf{q}}n({\mathbf{c}}_{\mathbf{q}}){\mathbf{c}}_{\mathbf{q}},\mathbf{\Pi }=\sum _{\mathbf{q}}n({\mathbf{c}}_{\mathbf{q}}){\mathbf{c}}_{\mathbf{q}}{\text{c}}_{\mathbf{q}}.$$

$${\Delta }_{\mathbf{q}}(n)=\sum _p{\mathcal{L}}_{\mathbf{qp}}{n}^{\mathrm{neq}}({\mathbf{c}}_{\mathbf{p}}),$$

$$n(\mathbf{r},{\mathbf{c}}_{{\mathbf{q}}^{\prime }},t+\Delta t)=n^{\ast }(\mathbf{r},{\mathbf{c}}_{\mathbf{q}},t)-\frac{2C^{\mathbf{q}}{\rho }_0{\mathbf{u}}_b\cdot {\mathbf{c}}_{\mathbf{q}}}{c_s^2},$$

$${\mathbf{u}}_b=\mathbf{U}+\mathbf{\Omega }\times ({\mathbf{r}}_b-\mathbf{R}),$$

$$\bm{f}^{h}(\bm{r}_{b},\bm{c}_{\bm{q}},t)=\frac{2\Delta x^{3}}{\Delta t}\Bigg{[}2n^{*}(\bm{r},\bm{c}_{\bm{q}},t-\Delta t)-\frac{2C^{\bm{q}}\rho_{0}\bm{u}_{b}\cdot\bm{c}_{\bm{q}}}{c_{s}^{2}}\Bigg{]}\bm{c}_{\bm{q}}.$$

$${\mathbf{F}}^h$$

$${\mathbf{T}}^h$$

$$\left(\begin{array}{c}{\mathbf{F}}_i^h\\ {\mathbf{T}}_i^h\\ {\mathbf{T}}_j^h\\ {\mathbf{\sigma }}_i^h\\ {\mathbf{\sigma }}_j^h\end{array}\right)=-\left(\begin{array}{ccc}{\mathbf{A}}_{ii}& -{\mathbf{B}}_{ii}& {\mathbf{B}}_{jj}\\ {\mathbf{B}}_{ii}& {\mathbf{C}}_{ii}& {\mathbf{C}}_{ij}\\ -{\mathbf{B}}_{jj}& {\mathbf{C}}_{ij}& {\mathbf{C}}_{jj}\\ {\mathbf{G}}_{ii}& {\mathbf{H}}_{ii}& {\mathbf{H}}_{ij}\\ -{\mathbf{G}}_{jj}& {\mathbf{H}}_{ji}& {\mathbf{H}}_{jj}\end{array}\right)\left(\begin{array}{c}{\mathbf{U}}_{ij}\\ {\mathbf{\Omega }}_i\\ {\mathbf{\Omega }}_j\end{array}\right),$$

$$A_{\alpha \beta }^{(ij)}=A_{\beta \alpha }^{(ji)}.$$

$$A_{\alpha \beta }^{(ij)}$$

$$B_{\alpha \beta }^{(ij)}$$

$$C_{\alpha \beta }^{(ij)}$$

$$G_{\alpha \beta \gamma }^{(ij)}$$

$$+Y_{ij}^G(n_{\alpha }^{(ij)}{\delta }_{\beta \gamma }+n_{\beta }^{(ij)}{\delta }_{\alpha \gamma }-2n_{\alpha }^{(ij)}{n}_{\beta }^{(ij)}{n}_{\gamma }^{(ij)}),$$

$$H_{\alpha \beta \gamma }^{(ij)}$$

$$X_{ii}^A$$

$$X_{ii}^G$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: *Yukawa Institute for Theoretical Physics, Kyoto University, Kitashirakawaoiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan; E-mail: pradipto@yukawa.kyoto-u.ac.jp *

1 Introduction

The behavior of suspended particles in solvents (suspensions) has been studied since Einstein published his seminal work in 19051. However, our understanding of rheological properties of suspensions is still limited. One of the interesting phenomena in the rheology of dense suspensions is the discontinuous shear thickening (DST) which is an abrupt jump of the viscosity at a critical shear rate $\dot{\gamma}$ if the volume fraction $\phi$ is larger enough. The DST can be observed easily in a mixture of cornstarch and water. The DST has also industrial applications such as traction controls and protective vests. Since the first observation of it about 90 years ago2, the DST has been studied extensively3, 4, 5. Some experiments clarified the roles of shape of particles6, 7, and the boundary effects on the DST8, 9.

We have substantial progresses on simulations of dense suspensions. By using the Lubrication-Friction Discrete Element Method (LF-DEM) 10, 11, 12 which is a simplified version of the Stokesian dynamics 13, Seto and his coworkers confirmed an important role of frictional contact force between particles for the DSTs in non-Brownian10, 11 and Brownian suspensions12. They successfully reproduced the DST which quantitatively agrees with experimental results. A good review along this line14 has been published recently. Despite these efforts, the theoretical understanding is still limited. By using a phenomenology that taking into account frictional contacts between grains15, Thomas et al. suggested that the anisotropy of the radial distribution function plays an important role in the DST16. Furthermore, the DST is also found in dry frictional granular particles17.

The isotropic jammed state18, 19 in which the system cannot flow above a critical volume fraction is well-defined, at least, for dry granular particles in the zero shear limit. Bi et al. introduced the shear-jammed state below the critical fraction of the isotropic jamming and the fragile state in an experiment for two-dimensional dry granular materials under a pure shear20. However, the definitions of these states are not clear if we discuss different setups. Recently, Otsuki and Hayakawa21 have proposed one definition of the shear jamming in oscillatory shear as a memory effect of the initial shear after reducing the strain amplitude. They have also adopted a new definition of the fragile state in which the system has or loses the rigidity depending on the initial phase (or the type of oscillatory shear) of the oscillatory shear. Moreover, they demonstrated that a DST-like phenomenon as a discontinuous jump of the dynamic viscosity can be observed in the fragile state21.

Until recently, there are only a few studies in suspensions under oscillatory shear, numerically 22, 23, 24 and experimentally25. Ness et al. also studied suspensions under an oscillatory shear below the jamming point where they observed the strain hardening and its frequency dependence in 3D systems in their experiment and numerical simulation26. One drawback of the oscillatory shear is that the storage and loss moduli 27 beyond the linear response regime are not well-defined 28. There are several proposals to handle the nonlinear response using so called the Fourier transform rheology29, 23 or Chebyshev polynomials decomposition30, 23. Instead of taking into account the nonlinearities, Ref. 21 adopts a protocol in which they initially use a finite or a large strain amplitude and then reduce the amplitude to be in the linear response regime.

The situations are more controversial if we discuss the shear-jammed and the fragile states in suspensions. The definitions of the shear-jammed state, and the fragile state, as well as their relations with the DST in suspensions are still unclear. It should be noted that the DST takes place below the shear jamming density in suspensions11. The distinction between shear jamming and the DST also exists in the experimental results by Peters et al.31, which suggested that the DST might correspond to the fragile state31. Meanwhile, another experiment suggests that DST takes place at the lower onset of the shear-jammed states32. In this paper, we demonstrate that the definitions of the shear jammed and the fragile states proposed by Ref. 21 can be used even in suspensions.

The LF-DEM only considers the lubrication interaction with ignoring the long range hydrodynamic interaction between particles. Although many people believe that the LF-DEM is sufficient in describing the rheology of dense suspensions near the DST and the shear jamming, the existence of the long-range interaction might affect the critical behavior slightly lower densities or the critical densities of such rheological transitions. Therefore, we adopt the lattice Boltzmann method (LBM) for suspensions33, 34 with its lubrication correction35 in order to calculate the hydrodynamic interaction between the particles. There are also several studies on suspensions using the LBM36, 37, 38, 39, 22, where this method has been confirmed to be efficient and accurate for sedimentation simulations40. Although we can apply the LBM to suspensions embedded in a fluid with finite Reynolds number, we restrict our analysis to suspensions embedded in a fluid with low Reynolds number limit in this paper.

The outline of this paper is as follows. In the next section, we explain our simulation method briefly. In section 3, we present the results of our simulations under the simple shear to verify the relevance of our LBM as well as some new aspects of our analysis. Then the results for the oscillatory shear are presented as the main result in section 4. In the last section, we summarize our results and discuss some future perspective. In Appendices, we describe some details of our analysis and some supplemental information to the main text.

2 Simulation Method

Let us consider a suspension consisting of spherical particles. A set of equations of motion of the suspended particles is given by

[TABLE]

where $m$, $I$, $\bm{U}$, and $\bm{\Omega}$ are the diagonal matrix for mass of the particles, the diagonal matrix for the moment of inertia of particles, the translation velocity, and the angular velocity, respectively. The forces are the sum of the hydrodynamic force $\bm{F}^{h}$, the contact force $\bm{F}^{c}$, and the electrostatic repulsive force $\bm{F}^{r}$, while the torques are the sum of hydrodynamic torque $\bm{T}^{h}$ and contact torque $\bm{T}^{c}$.

In the LBM, the hydrodynamic force $\bm{F}^{h}$ is calculated by computing the discrete distribution function on lattices. Due to this lattice-based calculation, one needs to discretize the unit of length and time into the lattice units $\Delta x$. This discretization is related to the stability and accuracy of the scheme. This paper adopts $\Delta x=0.24a_{\text{min}}$, where $a_{\rm min}$ is the radius of the smallest suspended particle. Details of the LBM can be seen in Appendix A. On the other hand, the LF-DEM includes only two-body lubrication interactions between particles of the Stokes flow.

Perfect spherical hard particles in the Stokes flow do not allow any contact since the lubrication force diverges if the gap between particles becomes zero. Here, we introduce a cutoff length $\delta$ of the lubrication forces to allow the contact between rough particles. Our simulations adopts $\delta=0.01a_{\text{min}}$, which is used in Ref. 11 ***the viscosity for $\delta=0.01a_{\text{min}}$ is about $20\%$ smaller than that for $\delta=0.001a_{\text{min}}$ under stress control simulation 11.

We adopt the linear spring model for the contact interaction between particles commonly used in discrete element method (DEM) 41 which involves both the normal and the tangential contact forces. Note that we omit the dissipative contact force for the tangential part. For particle $i$, the contact force $\bm{F}^{c}_{i}$ and torque $\bm{T}^{c}_{i}$ are, respectively, written as

[TABLE]

[TABLE]

where $a_{i}$ is the radius of the particle $i$. The normal force is expressed as $\bm{F}^{\text{nor}}_{ij}=(k_{n}\delta_{ij}^{n}+\eta^{(n)}u_{ij}^{(n)})\bm{n}_{ij}$ where $k_{n}$ is the spring constant, $\delta_{ij}^{n}$ is the normal overlap, $\bm{n}_{ij}$ is the normal unit vector between particles, $u_{ij}^{(n)}$ is the normal velocity of the contact point, and $\eta^{(n)}=\sqrt{m_{0}k_{n}}$ is the damping constant, where $m_{0}$ is the average mass of the particles. The stick tangential counterpart is represented as $\tilde{\bm{F}}^{\text{tan}}_{ij}=k_{t}\delta_{ij}^{t}\bm{t}_{ij}$ where $k_{t}$ is tangential spring constant equals to $0.2k_{n}$, $\delta_{ij}^{t}$ is the tangential compression and $\bm{t}_{ij}$ is the tangential unit vector at the contact point between particles $i$ and $j$. We also adopt the Coulomb friction rules, where $|\bm{F}_{ij}^{\text{tan}}|$ is replaced by $\mu|\bm{F}_{ij}^{\text{nor}}|$ if $|\tilde{\bm{F}}^{\text{tan}}_{ij}|\geq\mu|\bm{F}^{\text{nor}}_{ij}|$ for slip contacts and use $|\tilde{\bm{F}}^{\text{tan}}_{ij}|$ if $|\tilde{\bm{F}}^{\text{tan}}_{ij}|\leq\mu|\bm{F}^{\text{nor}}_{ij}|$ for stick contacts, whereas $\delta_{ij}^{t}$ is updated each time with the relative tangential velocity41. Here, $\mu$ is the friction coefficient for the particle interaction. All results of our simulation presented in this paper are obtained for $\mu=1.0$. Note that the normal spring constant $k_{n}$ determines the time scale $\tau$ as $\tau=\sqrt{m_{0}/k_{n}}$. In our simulation, $\tau$ is related to the time step of the LBM $\Delta t$ as $\Delta t=0.6\tau$.

The stress contribution from the contact forces is given by

[TABLE]

where $r_{ij,\alpha}$ is the distance between the center of masses of two contacting particles in $\alpha$ direction. Although there are asymmetric contributions of the stress tensor for frictional particles 42, 43, we have ignored such contributions because the asymmetric stress tensor is much smaller than the symmetric part.

To stabilize the suspensions, we introduce a double-layer (repulsive) electrostatic force in the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory44, 45, 46

[TABLE]

where $F_{0}=k_{B}T\lambda_{B}\hat{Z}^{2}(e^{a_{\text{min}}/\lambda}/(1+a_{\text{min}}/\lambda))^{2}/h^{2}$ with the charge number $\hat{Z}$, the Bjerrum length $\lambda_{B}$ and the Debye-Hückel length $\lambda$. Note that $\lambda_{B}$ can be expressed as $\lambda_{B}=e^{2}/(4\pi\epsilon_{0}\epsilon_{r}k_{B}T)$ where $e$, $\epsilon_{0}$, $\epsilon_{r}$, and $k_{B}$ are the elementary charge, the vacuum permittivity, the dielectric constant, and the Boltzmann constant, respectively 46. Here, we adopt the Debye length $\lambda=0.02a_{\text{min}}$. This stabilizing force is relevant to get quantitative agreements with the experiments as shown in Ref.12. The contribution to the stress from this double-layer force is given by

[TABLE]

The total shear stress contains all contributions from the hydrodynamic stress, the contact stress, and the electrostatic stress as

[TABLE]

We use bi-disperse particles which includes equal number of particles with ratio $a_{\text{max}}=1.4a_{\text{min}}$ to prevent them from crystallization. In both the simple and the oscillatory cases, we apply shears by moving bumpy walls made of particles (see Figs. 1 and 2), whereas we use the periodic boundary conditions for the other directions.

3 Simple shear simulation

This section consists of four parts. In the first part, we explain the setup of our simulation. In the second part, we demonstrate that our LBM simulation successfully reproduces the DST. In the third part, we illustrate the role of anisotropy of the local contact force on the particles, which can capture the behavior of the macroscopic stress tensor including the DST. In the last part, we compare the results of the LBM with those of the LF-DEM.

3.1 Setup

Our simulation for simple shear flows contains $N=450$ particles in a three-dimensional box. We use the volume fraction $\phi_{3}=2N\pi(a_{\text{max}}^{3}+a_{\text{min}}^{3})/(3V)$ to characterize the density of suspensions in 3D systems. The viscosity $\eta=\sigma_{xy}/\dot{\gamma}$ with the shear rate $\dot{\gamma}$ is time averaged between the strains $\gamma:=\dot{\gamma}t=5$ and $\gamma=11$. Then, we plot the viscosity against the dimensionless shear rate $\dot{\gamma}^{*}=6\pi\eta_{0}a_{\text{min}}^{2}\dot{\gamma}/F_{0}$ (see Fig. 1), where $\eta_{0}$ is the viscosity of the solvent fluid 11, 10. Note that we do not have to specify the value of $F_{0}$ in the simulation of simple shears because we only use $\dot{\gamma}^{*}$ in our simulation. This dimensionless shear rate $\dot{\gamma}^{*}$ is analogous to the Stokes number $St=m\dot{\gamma}/6\pi\eta_{0}a_{\text{min}}$. However, since the particles’ inertia is irrelevant, using $\dot{\gamma}^{*}$ is more appropriate than $St$. Using $\dot{\gamma}^{*}$ has been adopted in the previous numerical simulations 10, 11, 12.

3.2 Discontinuous shear thickening

First, we present the results for the viscosity under a simple shear where we observe discontinuous jumps of the viscosity in dense situations (e.g., $\phi_{3}=0.57$) above critical shear rates as shown in Fig. 3. These results semi-quantitatively agree with the DST observed in the LF-DEM simulations 10, 11 and experiments 47, 48, 49, 50, 51, 52, 53, 5, 54. The onset of the DST corresponds to the beginning of the frictional contact between grains. At low shear rate regime, the particles are prevented from contact by the hydrodynamic lubrication and the electrostatic repulsive forces.

3.3 Anisotropies of contact forces and stress

We quantify the role of the contact forces by analyzing the angular distributions of the contact forces55 $\zeta_{N}$ and $\zeta_{T}$ for the normal and the tangential parts, respectively, defined as

[TABLE]

The contact angle $\theta$ satisfies $\theta={\rm tan}^{-1}(r_{ij,y}/r_{ij,x})$ which is ranged between [math] and $\pi$ and calculated counterclockwise from the $x$-axis. Here, $\rho(\theta)$ is the angular distribution of the contact orientations. $\xi_{N}(\theta)$ and $\xi_{T}(\theta)$ are the angular distributions of intensities of normal and tangential forces, respectively, defined as

[TABLE]

Here, $F^{\text{nor}}(\theta)$ and $F^{\text{tan}}(\theta)$ are the normal and tangential forces in the direction of $\theta$ and $\langle F^{\text{nor}}\rangle$ is the average normal forces $\langle F^{\rm nor}\rangle=\int_{0}^{\pi}d\theta\rho(\theta)F^{\rm nor}(\theta)$. The angular distribution $\rho(\theta)$ satisfies the following normalization:

[TABLE]

Figures 4(a) and 4(b) present the angular distributions of the normal and tangential contact forces, respectively. Our results show that the anisotropy exists in the normal part, where $\zeta_{N}(\theta)$ takes maxima in the directions of compression ($\theta\approx 3\pi/4,7\pi/4$) and becomes minima in the expansion directions ($\theta\approx\pi/4,5\pi/4$). The tangential counterpart also shows anisotropy by exhibiting an 8-shaped figure, where the positive $\zeta_{T}(\theta)$ (counterclockwise rotations) is most likely oriented at the directions of shear gradient $\theta=\pi/2$ while the negative $\zeta_{T}(\theta)$ (clockwise rotations) is at the directions of shear. We can express the stress tensor as a function of these angular distributions as 55, 56

[TABLE]

where $Z$ is the average contact number, $\bm{n}_{\theta}=(\cos\theta,\sin\theta)$ and $\bm{t}_{\theta}=(-\sin\theta,\cos\theta)$. The derivation of this formula can be seen in Appendix C. Note that three-dimensional effect only appears as the geometrical factor in Eq. (13) because of the azimuthal symmetry. This is acceptable since the dominant part of the shear stress comes from the $x$ (shear) and $y$ (shear gradient) directions. In other words, we neglect the motion of particles in $z$-direction due to tangential forces. This formulation (Eq. (13)), in addition to the calculation of the lubrication interactions in the Appendix A 57, can reproduce the viscosity above the critical shear rate of the DST under the simple shear as shown in the dashed line in Fig. 3. Note that the main contribution for the shear rate larger than the critical one comes from the contact stress in Eq. (13). However, Eq. (13) can only be used for large shear rates since almost no contact exists for low shear rates.

We further explore the role of anisotropy as suggested from the angular distribution of the contact forces in Fig. 4 by analyzing the anisotropy of the stress tensor with the aid of $\Delta\lambda/P$. Here we introduce the stress anisotropy $\Delta\lambda$ as $\Delta\lambda=\sigma_{1}-\sigma_{3}$, where $\sigma_{1}$ and $\sigma_{3}$ are the maximum and the minimum eigenvalues of the stress tensor in 3D system. On the other hand, the pressure $P$ is calculated as $P=-(\sigma_{1}+\sigma_{2}+\sigma_{3})/3$ where $\sigma_{2}$ is the second largest eigenvalue of the stress. We have confirmed that this stress anisotropy increases as $\dot{\gamma}^{*}$ increases as shown in Fig. 5. This behavior is similar to in a stress-controlled simulation under a simple shear 16.

3.4 LBM vs LF-DEM

In this part, let us compare the results of our LBM simulation with those of the LF-DEM to clarify the role of long-range hydrodynamic interaction between particles. As seen in Fig. 6, one can get qualitatively similar results without using the LBM for high shear regime, though the viscosity by the LBM in the low shear regime quite differs from that by the LF-DEM. We also find that the slope of the continuous shear thickening obtained from the LBM is smaller than that for the LF-DEM. This might suggest the shift of the critical density of the DST if we consider the full hydrodynamic interactions. Thus, the LF-DEM gives only qualitative results. The long-range hydrodynamic interactions between particles play some roles if we are interested in a wide range of parameters. This is one merit to introduce the LBM simulation even for the rheology of dense suspensions.

4 Oscillatory shear simulation

This section summarizes the results of our simulation under oscillatory shear flows, which consists of five parts. In the first part, we explain the setup and the protocol of our simulation. In the second part, we summarize the results of our simulation using the storage and the loss moduli to characterize the mechanical responses under our protocol. In the third part, similar to the simple shear case, we illustrate the role of anisotropy of the local contact force to the macroscopic stress anisotropy. In the fourth part, we discuss the existence of the fragile phase which depends on the initial phase of the applying the oscillatory shear. In the last part, we compare our results of the LBM with those of the LF-DEM.

4.1 Setup and Protocol

Our simulation for oscillatory shears contains $N=450$ particles confined in a quasi-two-dimensional box with $L_{z}=4a_{\rm max}$ with the radius $a_{\rm max}$ of the largest particle. We discuss size dependence for simulations by the LF-DEM in the Appendix D. The motion of particles is considered as complete two-dimensional one to keep the monolayer configuration. Therefore, instead of using $\phi_{3}$ we use the area fraction $\phi=N\pi(a_{\text{max}}^{2}+a_{\text{min}}^{2})/(2L_{x}L_{y})$ to characterize the density in our monolayer system. Note that we simulate 3D spheres instead of 2D disks since the hydrodynamic interactions are only well defined in three dimensional systems. One of the reasons to use the monolayer configurations is to save computational time, which is commonly used by some previous simulations 16, 21. The second reason to use the configurations is the easiness of the visualization of the microstructure of particles. Furthermore, monolayered suspensions have also been studied in some experiments on lipid domains 58 and on fluid surfaces 59, 60, 61, 62. The snapshot of the monolayer configuration is shown in Fig. 2.

We implement a time-dependent oscillatory shear strain as

[TABLE]

where $\gamma_{0}$ is the strain amplitude, $\omega$ denotes the angular frequency, and $\Theta$ is the initial phase which characterizes the asymmetricity of the strain. This is equivalent as implementing a strain rate $\dot{\gamma}(t)=\omega\gamma_{0}\sin(\omega t+\Theta)$. We measure the mechanical response in terms of 27

[TABLE]

where $G^{\prime}$ is the storage modulus corresponding to the rigidity and $G^{\prime\prime}$ is the loss modulus which is proportional to the dynamic viscosity $\eta:=G^{\prime\prime}/\omega$. For later discussion we have introduced the dimensionless storage and loss moduli $G^{\prime}_{*}$ and $G^{\prime\prime}_{*}$ as $G^{\prime}_{*}=G^{\prime}/k_{n}$ and $G^{\prime\prime}_{*}=G^{\prime\prime}/k_{n}$. As a result, the dimensionless viscosity $\eta^{*}$ is given by $\eta^{*}=G^{\prime\prime}_{*}/\omega\tau$. We have also introduced the dimensionless stress $\sigma^{*}_{\alpha\beta}$ as $\sigma^{*}_{\alpha\beta}=\sigma_{\alpha\beta}/k_{n}$. As mentioned in the introduction, this formulation is only reliable in the linear response regime. Some authors use different definitions of $G^{\prime}$ and $G^{\prime\prime}$ to handle the nonlinear responses23, 26, 29, 30. Note that the stress-strain curve in the linear response regime is not expressed as a simple straight line because of the existence of $G^{\prime\prime}$. The linear response means that $G^{\prime}$ and $G"$ are independent of the strain amplitude. To avoid uncertainty of the definition of the linear response functions $G^{\prime}$ and $G^{\prime\prime}$ in the nonlinear regime, we follow the protocol introduced by Otsuki and Hayakawa 21 in which we reduce the strain amplitude $\gamma_{0}$ to be in the linear response regime as $\gamma_{0}^{F}$ for observation after initial cycles with the large initial strain amplitude $\gamma_{0}^{I}$.

In other words, our protocol tries to extract the memory effect of the initial oscillation. When $\gamma_{0}^{F}\leq 10^{-3}$, the storage modulus $G^{\prime}$ is almost independent of $\gamma_{0}^{F}$ (Fig. 7), which confirms that our observation of the mechanical responses is in the linear response regime. Therefore, all oscillatory shear simulations adopt $\gamma_{0}^{F}=10^{-4}$.

We have simulated range of $\omega$ as shown in Fig. 8 to check its dependencies and confirmed that the response is almost independent of $\omega$ for $\omega\tau\leq 1.1\times 10^{-3}$. Therefore, we use $\omega\tau=1\times 10^{-3}$ in our simulations. For electrostatic repulsive force between particles in the oscillatory shear, we adopt $F_{0}=2\times 10^{-4}k_{n}$. The control parameters of our simulation are the initial strain amplitude $\gamma_{0}^{I}$, the area fraction $\phi$, and the initial phase $\Theta$. Our results below are obtained after ten initial cycles and averaged over nine final observation cycles after the reduction of the strain amplitude. We average the results for three ensembles of different initial configurations. The convergences for the number of initial cycles $N^{I}_{c}$ and observation cycles $N^{F}_{c}$ can be seen on the Appendix B.

4.2 Mechanical responses

Let us look at the storage and the loss moduli in the linear response regime under the oscillatory shear. Our obtained results of the storage modulus and the viscosity under oscillatory shear are shown in Figs. 9 and 10, respectively. One can observe the existence of finite storage modulus for all $\gamma_{0}^{I}$ on $\phi=0.84$. For $\phi=0.82$, the storage and the loss moduli have discontinuous jumps at critical value of $\gamma_{0}^{I}$. Note that the storage modulus below the critical $\gamma_{0}^{I}$ is almost zero while the loss modulus in this region is small but finite.

From these results, we can identify the jammed, unjammed, and shear-jammed states with the aid of $G^{\prime}$. First we introduce a threshold value of the storage modulus, $G^{\prime,\text{th}}=10^{-3}k_{n}$. If $G^{\prime}>G^{\prime,\text{th}}$, we regard the system as having finite storage modulus. The unjammed state is when $G^{\prime}<G^{\prime,\text{th}}$. The isotropic jammed state is when $G^{\prime}>G^{\prime,\text{th}}$ for all $\gamma_{0}^{I}$ at a given $\phi$ or $G^{\prime}>G^{\prime,\text{th}}$ for low $\gamma_{0}^{I}$. Meanwhile, we define the shear-jammed state as the state when $G^{\prime}>G^{\prime,\text{th}}$ for high $\gamma_{0}^{I}$ and unjammed for low $\gamma_{0}^{I}$. From the results of the viscosity in Fig. 10, we observe finite jumps for $\phi=0.81$ and $\phi=0.82$, while the viscosity is almost independent of $\gamma_{0}^{I}$ on $\phi=0.84$. The jumps of the viscosity observed here correspond to the DST in the simple shear2, 3, 4, 10.

Both of the abrupt increases on the storage modulus and the viscosity can be understood as the appearance of a percolating force chains after the reduction of the strain amplitude. This clearly can be seen in Figs. 11 and 12 where we visualize the force chains before and after the reduction for $\phi=0.82$, $\gamma_{0}^{I}=0.1$ in the shear-jammed state and $\gamma_{0}^{I}=0.08$ in the unjammed state in Figs. 11 and 12. It is easy to find that the percolated force chains survive after the reduction of the strain amplitude in the shear-jammed state, while they disappear after the reduction in the unjammed state.

4.3 Anisotropy of the contact forces

Here, we explore the anisotropy of the stress tensor and the contact forces in our simulation by analyzing the angular distributions of the contact forces (Fig. 13) as in the simple shear case (Eqs.(8) and (9)). The angular distribution of the contact force for large $\gamma_{0}^{I}$ are almost equivalent to those observed in the simple shear. Moreover, the angular distribution of the normal contact force is almost isotropic for isotropic jammed state, while the stress anisotropy is clearly visible in the shear jammed state. For two-dimensional cases, we replace Eq. (13) by 55

[TABLE]

Figure 14 compares the theoretical shear stress in Eq.(17) with the results of our simulation. This agreement confirms that the contact contributions are dominant. However, this analysis cannot be applied to the unjammed states since almost no contact exists in these regimes, whereas only the hydrodynamic contributions exist as shown in Fig. 15.

4.4 Initial phase dependence and fragile state

Now let us clarify properties of the fragile state within our simulations. Originally, the fragile state is defined as the state where the system can only sustain load in a particular direction 63, 20. This suggests that both solid-like and fluid-like responses can be observed depending on the the initial phase $\Theta$ or the asymmetricity of the strain introduced in Eq. (14) in the fragile state. Therefore we try to explore how this duality can happen for given control parameters $\phi$, $\gamma_{0}^{I}$, and the initial phase $\Theta$. In our protocol, the different initial phase $\Theta$ essentially corresponds to the different position of the observation in the stress-strain curve (see Fig. 18). If the state is fragile, we expect that the response depends on $\Theta$.

In Fig. 16, we plot the storage modulus against $\gamma_{0}^{I}$ for $\phi=0.82$ and $\Theta=0$ and $\Theta=\pi/2$. Here we can see that the points at $\gamma_{0}^{I}=0.1$ and $\gamma_{0}^{I}=0.15$ have $G^{\prime}>G^{\prime,\text{th}}$ for $\Theta=0$ while they have $G^{\prime}<G^{\prime,\text{th}}$ for $\Theta=\pi/2$. Hence, we confirm that this $\Theta$-dependence can be used for the definition of the fragile state as in the dry granular case 21. We also analyze the stress anisotropy as in the simple shear by using $\Delta\lambda=\sigma_{1}-\sigma_{2}$ divided by the pressure $P=-\frac{1}{2}(\sigma_{1}+\sigma_{2})$. Here, $\sigma_{1}$ and $\sigma_{2}$ are the maximum and the minimum eigenvalues of the stress tensor, respectively. In Fig. 16, we plot the stress anisotropy averaged in the last five cycles after the reduction. We confirm that the stress anisotropy in this fragile state is much larger than that in the shear jammed state as expected in Refs. 20, 21. This observation agrees with the experiment of granular materials 64, the simulations of frictionless dry granular particles65 and frictional dry granular disks 21. It is reasonable that anisotropy is large in the fragile state while it is small in the shear jammed state, because the shear jammed state has more rigid percolating force chains to sustain in all directions while the fragile state has weak force chains which are connected only in the compressive direction.

The $\Theta$-dependence of $G^{\prime}$ can be seen clearly in Fig. 17 for $\phi=0.82$ and $\gamma_{0}^{I}=0.1$. Note that Fig. 17 is obtained from the LF-DEM simulation. Furthermore, the change of the observation point depending on $\Theta$ can be seen from the stress-strain curve in Fig. 18 which is obtained from the LBM. This $\Theta-$dependent relationship between the initial stress-strain curve and the observation point is the origin of the $\Theta$-dependent response.

To summarize the results of our oscillatory shear simulation, we draw the phase diagram for the dense monolayer suspensions under oscillatory shears (with our control parameters, area fraction $\phi$ and initial strain amplitude $\gamma_{0}^{I}$) in Fig. 19. Note that this phase diagram is drawn by using the LF-DEM because of the limitation of our computational resources. Our phase diagram is quite similar to that for the dry granular materials in the corresponding protocol 21. Furthermore, our distinction within the shear-jammed and isotropic jammed states is also similar to the experimental results by controlling the shear rate for cornstarch 32 with additional fragile states on the onset of the shear jamming states.

4.5 LBM vs LF-DEM

In this part, we clarify the difference between the LF-DEM and the LBM simulation. From Fig. 20, we can observe a qualitative agreement between both approaches for high shear regime. It is natural that the hydrodynamic long range interactions play a minor role for the emergence of exotic states such as the shear jamming and the fragile states. Nevertheless, the LF-DEM cannot be used for low density regime and the low shear regime because it does not contain the full hydrodynamic interactions.

5 Discussions and Conclusions

In conclusion, we have numerically studied dense non-Brownian suspensions under the simple shear and the oscillatory shear. We have successfully reproduced the DST in our simple shear simulation. For oscillatory shear, we have successfully characterized the unjammed, the jammed, the shear-jammed, and the fragile states as a memory effect of the initial shear by the reduction of the strain amplitude. We have drawn the phase diagram for monolayered dense non-Brownian suspensions on the plane of the initial strain amplitude and the area fraction. We have found that the viscosity under oscillatory shear also has an abrupt increase which corresponds to the DST. Furthermore, we have defined the fragile state as a state which depends on the asymmetricity of the strain. Moreover, we have also analyzed the angular distributions of the contact forces characterizing the DST and the shear jamming under the simple shear and the oscillatory shear, respectively. We have expressed the shear stress, the storage modulus at high initial strain amplitude, and the viscosity for large shear rate as functions of these angular contact distributions. The agreement between this phenomenology and the results of our simulation clarifies that the contact contributions are much larger than those of the hydrodynamics in the shear-jammed regime under the oscillatory shear and for large shear rate under the simple shear. This agreement also indicates that one needs to construct a theory for the angular distribution66 of contact forces, i.e., derive the $\zeta_{N}$ and $\zeta_{T}$ theoretically to recover the mechanical responses at high strain regime and the viscosity for large shear rate. Finally, we have also illustrated that the stress anisotropy takes its maximum on the onset of the shear jamming, where it is corresponds to the fragile state and DST.

Let us discuss our future perspectives. There are several unsolved problems which have not been addressed in this paper. First, the effect of dilatancy is missing in our analysis because the system is incompressible. However, the effect of dilatancy becomes significant in sand beaches, which are mixtures of air, water, and grains. Recently, dilatancy and compaction of dry granular materials are studied under an oscillatory shear and in a pressure-control simulation67, where the anisotropy of the stress tensor also plays an important role. Therefore, we need to perform the simulation of a mixture of air, water, and grains. Second, one can also investigate the rheology of an intruder particle and measure the drag force and effective friction coefficient, which has been performed experimentally on the granular hydrogels immersed in waters68. Finally, we can also numerically study a tapping problem or impact-activated solidification which has been experimentally showed that it is related to the formation of the jamming front69, 70, 71. One may simulate this process by the coupled DEM-LBM for the free surface simulation72.

Appendix A Details of the simulation method

Lattice Boltzmann Method

To calculate the hydrodynamic forces acting on the particles in the LBM, we introduce solid nodes inside the particles, fluid nodes on the surrounding fluid, and boundary nodes on the surface of the particles. The hydrodynamic field is calculated from the time evolution of the discrete distribution function $n(\bm{r},\bm{c}_{\bm{q}},t)$ which has the dimension of mass density. At each fluid node $\bm{r}$, $n(\bm{r},\bm{c}_{\bm{q}},t)$ is updated as

[TABLE]

where $\bm{c}_{\bm{q}}$ is the lattice velocity with direction $\bm{q}$, $\Delta t$ is the time step, and $\Delta_{\bm{q}}$ is the collision operator that depends on all $n(\bm{r},\bm{c}_{\bm{q}},t)$. The left-hand side of Eq. (18) expresses the change of the distribution function by the streaming effect, whereas the right-hand side expresses the change by the collision. Our simulation uses nineteen directions of $\bm{c}_{\bm{q}}$ (D3Q19 lattice for 19 directions/quadratures in 3 dimensions) 73. Some moments of $n(\bm{c}_{\bm{q}})$ which is the abbreviation of $n(\bm{r},\bm{c}_{\bm{q}},t)$ are related to the hydrodynamic variables as

[TABLE]

In this paper, we use the linearized collision operator near the equilibrium distribution

[TABLE]

where $n^{\rm neq}(\bm{c}_{\bm{p}})=n(\bm{c}_{\bm{p}})-n^{\rm eq}(\bm{c}_{\bm{p}})$ is the deviation from the equilibrium distribution and $\Delta_{\bm{q}}(n)$ satisfies $\Delta_{\bm{q}}(n^{\rm eq})=0$. It is not necessary to specify an explicit collision operator $\mathcal{L}_{\bm{qp}}$ for the calculation of the hydrodynamic variables introduced in Eq. (19). Instead, it is sufficient to consider the eigenvalues and eigenvectors of $\mathcal{L}_{\bm{qp}}$ for the collision process. Details of this collision operator are not discussed in this paper and can be found in Refs. 33, 34, 35.

The discrete distribution function on the boundary node is updated with a bounce-back rule on the surface of each particle. At each boundary node, there exists virtual (post-collisional) distribution function $n^{*}(\bm{r},\bm{c}_{\bm{q}},t)=n(\bm{r},\bm{c}_{\bm{q}},t)+\Delta_{\bm{q}}(\bm{r},t)$ instead of using Eq. (18). Moreover, we need an additional term for the time updating of the boundary nodes

[TABLE]

where $\bm{c}_{\bm{q}^{\prime}}=-\bm{c}_{\bm{q}}$, i. e. $n(\bm{r},\bm{c}_{\bm{q}^{\prime}},t)=n_{\bm{q}}(\bm{r},-\bm{c}_{\bm{q}},t)$, and $c_{s}$ is the lattice sound speed, $c_{s}=\sqrt{c^{2}/3}$ with $c=\Delta x/\Delta t$. Here, we choose the coefficient $C^{\bm{q}}$ equals to $1/3$ if the vector $\bm{q}$ is at rest, $1/18$ at the nearest neighbor site and $1/36$ at the second nearest neighbor site, respectively 33, 34, 35. The velocity of the boundary node $\bm{u}_{b}$ satisfies the relation

[TABLE]

where $\bm{r}_{b}$ is the position of the boundary nodes $\bm{r}_{b}=\bm{r}+\frac{1}{2}\bm{c}_{\bm{q}}\Delta t$ (halfway between the fluid and the solid nodes), where $\bm{R}$ is the center of mass of the particle. As a result, the hydrodynamic force $\bm{f}^{h}(\bm{r}_{b},\bm{c}_{\bm{q}},t)$ exerted on the boundary node is calculated from the momentum transferred in Eq. (21) as

[TABLE]

For simplicity, we abbreviate the force at each time $\bm{f}^{h}(\bm{r}_{b},\bm{c}_{\bm{q}},t)$ to be $\bm{f}^{h}(\bm{r}_{b})$. Then, the hydrodynamic force at each time $\bm{F}^{h}$ and torque $\bm{T}^{h}$ exerted on the particle are just sum of the forces and torques on all boundary nodes on the particle as

[TABLE]

With the aid of the hydrodynamic force $\bm{f}^{h}(\bm{r}_{b})$, the $\alpha,\beta$ element of the hydrodynamic stress exerted on $i-$th particle $\sigma_{i,\alpha\beta}^{h}$ is given by $\sigma^{h}_{i,\alpha\beta}=-\sum_{b}r_{b,\alpha}f^{h}_{\beta}(\bm{r}_{b})/V$ where $V=L_{x}L_{y}L_{z}$ is the volume of the system with the linear dimension $L_{\alpha}$ of $\alpha$ direction. Then, we sum up all the contributions from particles to obtain the hydrodynamic stress $\sigma^{h}_{\alpha\beta}=\sum_{i}\sigma^{h}_{i,\alpha\beta}$.

Note that there are two technical difficulties to adopt this formulation with the aid of the boundary nodes in simulating dense suspensions. (i) The formulation based on boundary nodes crashes when the particles are overlapping. This forces us to introduce a contact radius larger than the mapped radius for the boundary nodes with ratio $a_{\text{map}}=0.95a_{\text{contact}}$. The mentioned radius in this paper is the contact radius. (ii) This scheme works well for the gap between particles $h>0.5\Delta x$ and is inaccurate for $h\leq 0.5\Delta x$ due to shared nodes between two particles in the opposite case 34. Therefore, we need to incorporate the lubrication correction when the gap between particles is small. This correction is calculated by the grand-resistance matrix formulation of pairwise lubrication interaction 57, 35, 74

[TABLE]

where $\bm{F}^{h}_{i}$ is the hydrodynamic (lubrication) force on particle $i$ and for particle $j$, $\bm{F}^{h}_{j}=-\bm{F}^{h}_{i}$. $\bm{U}_{ij}=\bm{U}_{i}-\bm{U}_{j}$ is the relative velocity. We adopt the notation in Ref. 57, 74 such as $\bm{A}_{ij}=(A^{(ij)}_{\alpha\beta})$ and due to the Lorentz reciprocal theorem, we have the symmetry relations such as

[TABLE]

For axisymmetric geometries, the coefficients can be expressed in terms of scalar functions as

[TABLE]

where $n^{(ij)}_{\alpha}$ is the normal unit vector between particles $i$ and $j$ in the $\alpha$ direction and $\epsilon_{\alpha\beta\gamma}$ is Levi-Citiva symbol, i. e. $\epsilon_{\alpha\beta\gamma}=1$ for an even permutation of $(\alpha,\beta.\gamma)$, $\epsilon_{\alpha\beta\gamma}=-1$ for an odd permutation of $(\alpha,\beta.\gamma)$, and $\epsilon_{\alpha\beta\gamma}=0$ otherwise. The scalar functions $X$ and $Y$ are functions of interparticle gap $h$. For two spheres of arbitrary size with the leading order only, the scalar functions are written as

[TABLE]

where $\beta$ is the ratio of the particle radius defined as $\beta=a_{i}/a_{j}$ and $\delta$ is the roughness length.

Appendix B Dependences on number of initial and observation cycles

In this section, we verify that the results for the mechanical responses in the main text is not merely transients. First, we plot the storage modulus for $\phi=0.82$ against initial strain amplitude $\gamma_{0}^{I}$ to see how the results depend on number of observation cycles $N_{c}^{F}$. Our result in Fig. 21 indicates that the results converges after 5 cycles. Meanwhile, for the number of initial cycles, we also get convergence after eight initial cycles at least, for shear jammed states as seen in Fig. 22.

Appendix C Brief derivation of the stress formulation

In this section we describe the derivation of the stress formulation by following Ref. 56. Let us consider a system consisting of $N$ particles with average radius $a$. The particle radius is assumed to be uniform. From the contact angle $\theta={\rm tan}^{-1}(r_{ij,y}/r_{ij,x})$ (Fig. 23) which is ranged between [math] and $2\pi$ and calculated counterclockwise from the $x$-axis, one can define $\rho_{\theta}$ as the angular distribution of the contact orientations. $\rho_{\theta}$ satisfies the following normalization

[TABLE]

Note that the azimuthal angle is ignored here. Then, $F^{N}n_{\alpha}$ is the $\alpha-$ component of the normal contact forces vector on a contact point. Suppose that the particles are in mechanical equilibrium, the force balance in $\theta$ directions

[TABLE]

is satisfied, where $Z$ is the coordination number or the average number of contact points on each particle. Here, $n_{\alpha}\rho_{\theta}$ represents the expectation value of the contact unit vector. The torque balance is written as

[TABLE]

Suppose that the material is subject to a macroscopically uniform stress. Then, one can assume the first order deformation of the material

[TABLE]

under the deformation the point $x_{\alpha}$ to $x^{\prime}_{\alpha}$. Thus, the displacement $u_{\alpha}=x^{\prime}_{\alpha}-x_{\alpha}$ is given by

[TABLE]

where $\delta_{\alpha\beta}$ is the Kronecker delta. The distortion tensor $\Gamma_{\alpha\beta}$ contains the symmetric and anti-symmetric part as

[TABLE]

where $E_{\alpha\beta}$ is the strain (symmetric) tensor and $R_{\alpha\beta}$ is the rotation (anti-symmetric) tensor.

Here, we assume that the particles also experienced a virtual deformation. This virtual deformation distorts a spherical particle into an ellipsoid due to the macroscopic deformation. Therefore, the small displacement $\epsilon_{\alpha}$ in the contact direction $n_{\alpha}$ is given by

[TABLE]

Another assumption is that the contact forces are not changed during this virtual deformation. Then, the virtual work done by the contact forces can be expressed in terms of angular integral (Eq. (43)) as

[TABLE]

By the virtue of Eq. (47) and the torque balance equation (Eq. (44)), one can rewrite the virtual work Eq. (49) as

[TABLE]

Now, the number of particles in unit volume is $\phi_{3}/(4/3)\pi a^{3}$, where $\phi_{3}$ is the solid volume fraction of the material in three dimensions. Hence, the virtual work done in unit volume $W$ is

[TABLE]

Meanwhile, the virtual work done per unit volume by the virtual strain $E_{\alpha\beta}$ under stress $\sigma_{\alpha\beta}$ is

[TABLE]

Therefore, the stress tensor can be expresssed as

[TABLE]

Then, if we define the angular distributions of intensities of the normal contact forces $\xi_{N}=F^{N}/\langle F^{N}\rangle$, and the angular distributions of the normal contact forces $\zeta_{N}=\rho_{\theta}\xi_{N}$, where $\langle F^{N}\rangle$ is the average normal forces, one could have

[TABLE]

Finally, one can also take into account the tangential forces by altering Eq. (54) as

[TABLE]

where $t_{\alpha}$ is the tangential unit vector. Note that for a system of disks, the number of particles in the unit area is $\phi/\pi a^{2}$, with $\phi$ is the area fraction. Therefore, for 2-dimensional disks, Eq. (55) is rewritten as

[TABLE]

Appendix D System size dependence

In this Appendix, we examine the system size dependence for the oscillatory shear simulations as shown in Fig. 24 by using LF-DEM. Here, one does not have to worry about serious finite size effects for $N\geq 300$. Because of the limitation of our computational resources we have not examined the finite size effects for our LBM simulation. Nevertheless, we believe that the finite size effect is not severe for $N\geq 300$ even for the LBM simulation because the results of the LBM simulation are not quite different from those obtained from the LF-DEM.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

One of the authors (P) expresses his gratitude to A. J. C. Ladd for providing the Lattice Boltzmann susp3d code and S. Takada for his 3D-DEM code. The authors thank M. Otsuki, R. Seto, A. Santos, and D. Ishima for essential comments and fruitful discussions. This work is partially supported by the Grant-in-Aid of MEXT for Scientific Research (Grant No. 16H04025) and the Programs YITP-T-18-03 and YITP-W-18-17. All numerical calculations were carried out at the Yukawa Institute for Theoretical Physics (YITP) Computer Facilities, Kyoto University, Japan.