Scaling laws for frictional granular materials confined by constant pressure under oscillatory shear

TL;DR

This study numerically investigates the rheological behavior of 2D frictional granular materials under oscillatory shear, revealing scaling laws for moduli and linking them to contact force distributions and physical mechanisms.

Contribution

It introduces new scaling laws for the rheology of frictional granular materials under oscillatory shear and explains their physical origins.

Findings

Identified scaling laws for storage and loss moduli

Linked scaling exponents to contact force distributions

Explained the physical mechanisms behind the scaling laws

Abstract

Herein, we numerically study the rheology of a two-dimensional frictional granular system confined by constant pressure under oscillatory shear. Several scaling laws for the storage and loss moduli against the scaled strain amplitude have been found. The scaling laws in plastic regime for large strain amplitude can be understood by the angular distributions of the contact force. The scaling exponents are estimated by considering the physical mechanism.

Equations

$${\gamma }_{0,\text{eff}}:=\frac{2A}{L_0}.$$

$$m_w\frac{d{v}_{w,y}^{\pm }}{dt}=\pm (P_w^{\pm }-P)L_x-{\xi }_{\text{d}}{v}_{w,y}^{\pm },$$

$$m_i\frac{d^2{\mathbf{r}}_i}{d{t}^2}=\sum _{j\ne i}{\mathbf{f}}_{ij},$$

$$I_i\frac{d{\omega }_i}{dt}=\sum _{j\ne i}{T}_{ij},$$

$$T_{ij}=-\frac{d_i}{2}{\mathbf{f}}_{ij}\cdot {\mathbf{t}}_{ij}$$

$${\mathbf{f}}_{ij}=\left({\mathbf{f}}_{ij,n}+{\mathbf{f}}_{ij,t}\right)\Theta (d_{ij}-r_{ij}),$$

$${\mathbf{f}}_{ij,n}$$

$${\mathbf{f}}_{ij,t}$$

$${\sigma }_{\alpha \beta }^{\text{sym}}:=-\frac{1}{2L_x{L}_y}\sum _i\sum _{j>i}(x_{ij,\alpha }{f}_{ij,\beta }+x_{ij,\beta }{f}_{ij,\alpha }).$$

$$\sigma (t)$$

$$G^{\prime }({\gamma }_{0,\text{eff}},P)$$

$$G^{\prime \prime }({\gamma }_{0,\text{eff}},P)$$

$$\tilde{\sigma }(t):={\sigma }_{xy}^{\text{sym}}(t)-\overline{{\sigma }_{xy}^{\text{sym}}}.$$

$$G^{\prime }$$

$$G_{\text{res}}^{\prime }(\hat{P}):$$

$$\underset{x\to 0}{\lim }{\mathcal{G}}^{\prime }(x)$$

$${\beta }_1^{\prime }=1.00\pm \mathrm{0.03.}$$

$$G_{\text{res}}^{\prime }/k_n=G_1^{\prime }+a{\hat{P}}^{1/2},$$

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

$$\underset{x\to 0}{\lim }{\mathcal{G}}^{\prime \prime }(x)$$

$${\beta }_1^{\prime \prime }=1.05\pm \mathrm{0.01.}$$

$$G^{\prime \prime }/k_n\sim \hat{\Omega }{\hat{\xi }}_n({\gamma }_{0,\text{eff}}/\hat{P}<x_c),$$

$$G^{\prime \prime }/k_n\sim \frac{{\gamma }_{0,\text{eff}}}{\hat{P}}(x_c<{\gamma }_{0,\text{eff}}/\hat{P}<1).$$

$${\sigma }_{\alpha \beta }\approx -2P{\int }_0^{2\pi }d\theta \left[{\zeta }_N(\theta )n_{\theta ,\alpha }-{\zeta }_T(\theta )t_{\theta ,\alpha }\right]n_{\theta ,\beta }$$

$$\sigma$$

$$=-\pi P(b-c).$$

$$G^{\prime }$$

$$G_{\text{res}}^{\prime }(\hat{P}):$$

$$\underset{x\to 0}{\lim }{\mathcal{G}}^{\prime }(x)$$

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

$$\underset{x\to 0}{\lim }{\mathcal{G}}^{\prime \prime }(x)$$

$$R:=$$

$$\hat{R}:=\frac{\tilde{R}}{\tilde{\sigma }(\Omega t=\pi /4)},$$

$$0.4\le \frac{d\log G^{\prime \prime }}{d\log {\gamma }_{0,\text{eff}}}\le \mathrm{0.8.}$$

$$G_0^{\prime \prime }:=\underset{{\gamma }_{0,\text{eff}}\to 0}{\lim }{G}^{\prime \prime }.$$

Full text

I Introduction

Amorphous materials consisting of repulsive and dissipative grains such as granular materials Gra , colloidal suspensions Col , bubbles For , and emulsions Zhang05 exhibit characteristic viscoelastic properties. Such materials behave as solids if the grains are densely packed, while they behave as liquids if the grains are not sufficiently dense. There is a jamming transition at a critical density where the viscosity diverges and the rigidity emerges above the critical density at zero load limit. Since the pioneer work by Liu and Nagel LNjam , studies on the jamming transition have attracted considerable attention among physicists. The jamming transition for frictionless grains is known as a mixed transition in which the coordination number changes discontinuously while the pressure and the shear stress change continuously at the jamming point Djam ; Mjam1 ; Tjam ; Ojam ; Mjam3 ; Mjam . Such a continuous change of the stress near the jamming point satisfies a number of scaling laws Esca1 ; Osca ; Hsca ; Wsca2 ; Esca2 ; OHsca .

It is remarkable that the introduction of mutual frictions between the grains drastically changes the behavior of the jamming transition including the existence of a discontinuous shear thickening (DST)  SHsca ; Ssca2 ; dst1 ; Sisca ; OHdst1 ; SMdst ; dst2 ; OHdst2 ; Kdst . Moreover, the mutual friction cannot be ignored for rigid grains to be consistent with Newton’s equation of motion. We have recognized that shear jammed states, commonly observed in frictional systems, differ from the conventional isotropic jammed states Zsj ; Exsj1 ; Exsj2 ; Exsj3 ; Nusj1 ; Nusj2 ; OHsj ; Pradipto . Several reserchers are interested in the mutual relationship between shear jamming and DST OHsj ; Pradipto ; Fall15 ; Peters16 . We realize that some shear jamming has originated from a memory effect of mechanical training OHsj ; Pradipto ; Kumar ; Jin17 ; Urbani ; Jin18 .

If an oscillatory shear is applied to such a material, the rigidity and viscosity can be measured simultaneously LAOS1 ; LAOS2 ; HyunRev ; LAOS3 . Recently, scaling laws on the rigidity and viscosity for frictionless grains under oscillatory shears have been studied Tosc ; Otsuki14 ; BBosc ; BSosc ; BTosc . In particular, it is interesting that two distinct regions exist, softening and yielding, when considering large strain amplitude BSosc . Previous studies have discussed scaling laws for frictional grains under simple shear OHdst1 and oscillatory shear OHdst2 ; however, thus far, most of these studies are based on situations under constant volume conditions.

Dilatancy is widely known as one of the typical characteristics of granular systems in which the density decreases as the shear increases under constant pressure conditions Rdil ; Tdil ; Sdil ; MiDi04 ; Fdil ; Bdil . Because the dilatancy cannot be observed in a constant volume system, it is important to understand the rheological properties of granular systems confined by the constant pressure under oscillatory shear.

The purpose of this study is to determine the viscoelastic properties of an oscillatory sheared granular system. For this purpose, we numerically study a granular system confined by constant pressure under oscillatory shear to extract scaling laws for the storage and loss moduli. The results suggest the absence of the distinction between softening and yielding in frictional systems, which contrasts with that of frictionless cases BSosc .

The contents of this paper are as follows. In the next section, we explain the setup of our numerical simulation. Section III illustrates the relevance of several scaling laws for the storage and loss moduli with a fixed friction constant for $\mu=1$, where $\mu$ is Coulomb’s friction constant. In the plastic regime under large strain amplitude, the semi-quantitative behavior of the stress can be understood by the consideration of the angular distributions of contacts. In Sec. IV we demonstrate that the scaling laws observed for $\mu=1$ in Sec. III are essentially independent of $\mu$. In the last section, we discuss and summarize our results. In the appendices, we explain the technical details and present some supplemental information which has not been described in the main text.

II simulated system

Our system is a two-dimensional one containing $N$ frictional granular disks. To avoid crystallization, we adopt a tetradisperse system in which the numbers of grains for the diameters $d_{0},\ 0.9d_{0},\ 0.8d_{0},$ and $0.7d_{0}$ are $N/4$ Crys . We assume that the density of each grain is identical and, thus, the mass is proportional to the square of the diameter of the grain. The mass corresponding to the largest diameter $d_{0}$ is denoted by $m_{0}$ in this study.

We assume that gravity is negligible in our system corresponding to grains on a smooth horizontal plane. We adopt the periodic boundary condition in the shear ($x-$)direction, while we apply the external pressure through the wall in the $y-$direction. The grains are confined in a square box, i.e. $L_{0}:=L_{y}=L_{x}$ in the absence of any external pressure, where $L_{x}$ and $L_{y}$ are the linear dimensions in the $x-$ and $y-$directions, respectively. When the external pressure is applied, $L_{y}(t)$ is no longer a constant but depends on time. The position of the center of mass of the wall at $y=\pm L_{y}(t)/2$ in the $x-$direction $x_{G}^{+}(t)$ obeys $x_{G}^{\pm}(t)=\pm A\sin(\Omega t)$, where $A$ and $\Omega$ are the strain amplitude and the angular frequency, respectively (see Fig. 1). Because the linear dimension of the system size in the $y-$direction is not constant, we introduce the effective strain amplitude

[TABLE]

To describe the motion of the walls under the external pressure $P$, we adopt the following equation:

[TABLE]

where $\pm$, $m_{w}$, $v_{w,y}^{\pm}$, $P_{w}^{\pm}$, and $\zeta_{\textrm{d}}$ are the superscripts for the wall at $y=\pm L_{y}(t)/2$, the mass of the wall satisfying $m_{w}=N_{w}m_{0}$, $y-$component velocity of the wall at $y=\pm L_{y}(t)/2$, the inner pressure acting on the walls from the grains, and the damping constant, respectively.

The translational and rotational equations of the motion of the $i-$th grain whose mass, position, and moment of inertia are, respectively, $m_{i}$, $\bm{r}_{i}$, and $I_{i}=m_{i}d_{i}^{2}/8$, respectively, with the diameter $d_{i}$ are given by

[TABLE]

where we have introduced the contact force $\bm{f}_{ij}$ and the torque $T_{ij}$ from the $j-$th grain acting on the $i-$th grain. Note that $\omega_{i}$ is the $z-$component of the angular velocity of the $i-$th grain. The torque $T_{ij}$ introduced in Eq. (4) satisfies

[TABLE]

with the tangential unit vector $\bm{t}_{ij}$ at the contact point between the $i-$th and $j-$th grains. The contact force of each grain is described by the Cundall-Stack model dem1 ; dem2 . Considering the contact force $\bm{f}_{ij}$ between the $i-$th and $j-$th grains, we inntroduce $d_{ij}=(d_{i}+d_{j})/2$ and $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$ with $r_{ij}:=|\bm{r}_{ij}|$, the contact force $\bm{f}_{ij}$ is expressed as

[TABLE]

where we have introduced the normal contact force $\bm{f}_{ij,n}$, the tangential contact force $\bm{f}_{ij,t}$ and Heaviside’s step function $\Theta(x)$ satisfying $\Theta(x)=1$ for $x\geq 0$ and $\Theta(x)=0$ otherwise. The normal force consists of the elastic repulsive force represented by a linear spring (the spring constant $k_{n}$) and the dissipative force represented by a dashpot (the damping constant $\xi_{n}$) as

[TABLE]

where we have introduced $\bm{n}_{ij}=\bm{r}_{ij}/r_{ij}$, $\bm{v}_{ij,n}=(\bm{v}_{ij}\cdot\bm{n}_{ij})\bm{n}_{ij}$ with $\bm{v}_{ij}=d\bm{r}_{ij}/dt$ and $\delta_{ij,n}=d_{ij}-r_{ij}$. On the other hand, the tangential force $\bm{f}_{ij,t}$ contains both a sticking state and a slip state, where the switching from the stick to the slip takes place if the magnitude of the tangential force $f_{ij,t}=|\bm{f}_{ij,t}|$ exceeds the critical condition as

[TABLE]

where $\delta_{ij,t}$ is the tangential displacement during a contact time. We have also introduced Coulomb’s friction coefficient $\mu$, the tangential spring constant $k_{t}$, the tangential damping constant $\xi_{t}$, $f_{ij,n}=|\bm{f}_{ij,n}|$ and the relative tangential velocity at the contact point $\bm{c}_{ij}:=\bm{v}_{ij}-\bm{v}_{ij,n}+\bm{t}_{ij}(d_{i}\omega_{i}+d_{j}\omega_{j})/2$.

For our simulation, we adopt $k_{t}=k_{n}/2$, $\xi_{n}=(m_{0}k_{n})^{-1/2}$, $\xi_{t}=\xi_{n}$, and $\xi_{\textrm{d}}=\xi_{n}$ which follow Refs. OHsj ; Cruz . The control parameters of our simulation are $\hat{P}:=P/k_{n}$, $\gamma_{0,\rm{eff}}$ and $\Omega\sqrt{m_{0}/k_{n}}$. We cover the range of $\hat{P}$ from $2.0\times 10^{-5}$ to $6.0\times 10^{-3}$. In most of our simulations we adopt the number of grains $N=4,000$, and the frequency $\Omega/(2\pi)=1.0\times 10^{-4}\sqrt{k_{n}/m_{0}}$. We have confirmed that the rigidity and the viscosity are independent of $N$ in the range of $4,000\leq N\leq 20,000$. We have checked the dependence of $\Omega$ in the range of $4.0\times 10^{-5}\sqrt{k_{n}/m_{0}}\leq\Omega/(2\pi)\leq 1.0\times 10^{-3}\sqrt{k_{n}/m_{0}}$ (see Appendix B) and confirmed that the storage modulus is independent of $\Omega$ and the loss modulus depends on $\Omega$ in the small strain regime but does not depend on $\Omega$ for large strain.

To make the initial configuration isotropic, we randomly place all grains whose diameters are $0.6d_{0},\ 0.5d_{0},\ 0.4d_{0},\textrm{and}\ 0.3d_{0}$ corresponding to tetradispersity of grains confined in a square box with the initial density $\phi_{\textrm{ini}}=0.24$ without any overlap and shear. After we increase the diameter of each grain by $0.02d_{0}$, we wait for the system to relax to a steady state in which the potential energy $V_{n}:=k_{n}\sum_{i\neq j}\delta_{ij,n}^{2}/2$ for the normal compression satisfies $V_{n}/(Nk_{n}d_{0}^{2})<2.0\times 10^{-7}$. Subsequently, we increase the diameter of each grain by a further $0.02d_{0}$ again. After repeating this process, we reach a desired state in which the area fraction is $\phi=0.82$. We compress both walls by pressure $P$ to compactify the system above the jamming density to achieve an isotropic structure as the initial condition. After $V_{n}$ converges to the steady value, we apply oscillatory shear. We adopt the symplectic Euler method with the time step $\Delta t=0.05\sqrt{m_{0}/k_{n}}$.

We average the data over $N_{T}$ cycles after abandoning the data in the initial $N_{\textrm{ini}}$ cycles. Because the obtained results depend on $N_{\textrm{ini}}$ for $N_{\textrm{ini}}\leq 5$, we adopt $N_{\textrm{ini}}=10$ and $N_{T}=10$. Furthemore, we have confirmed that the storage and loss moduli are independent of $N_{\rm ini}\leq 190$ except for the low pressure cases and the critical region of the yielding transition. We set $t=0$ at the instant when $N_{\textrm{ini}}$ oscillatory shear ends.

Figure 2 shows the time evolution of the force chains during one cycle in which the force chains are unchanged if the strain amplitude is small; however, they are changed if the strain amplitude is large. The former regime, in which the granular materials behave as an elastic solid, is elastic, while the latter one, in which the granular materials behave as a plastic material with configuration changes of the grains, is plastic. The yielding transition is the transition between the two states.

Because the contact stress is dominant for the stress in dense granular systems Cruz ; rotTen , we obtain the symmetric contact stress $\sigma_{\alpha\beta}^{\textrm{sym}}$ as

[TABLE]

The asymmetric part in the stress tensor is associated with the coupled stress in frictional systems  Mitarai02 ; Gol . Nevertheless, we have confirmed that the asymmetricity of the stress tensor is so small that we can ignore the asymmetric stress tensor or the coupled stress (see Appendix A). We note that the effect of the shear band near the boundaries is not negligible for low pressure and large strain amplitude. In this case, the measured stress from Eq. (9) might be different from the local stress.

In systems under oscillatory shear, the storage modulus $G^{\prime}$ and the loss modulus $G^{\prime\prime}$ represent the rigidity and dynamic viscosity multiplied by $\Omega$, respectively. Here, the shear stress $\sigma(t):=\sigma_{xy}^{\textrm{sym}}$ can be decomposed into the elastic and viscous parts as

[TABLE]

where $\sigma^{\textrm{(ela)}}(t)$ and $\sigma^{\textrm{(vis)}}(t)$ are expressed as $\sigma^{\textrm{(ela)}}(t)=G^{\prime}\gamma(t)$ and $\sigma^{\textrm{(vis)}}(t)=\eta_{d}\dot{\gamma}(t)$ with the strain $\gamma(t)$ and the strain rate $\dot{\gamma}(t)$, respectively. It is known that the dynamic viscosity can be expressed as the loss modulus $G^{\prime\prime}$ as $\eta_{d}:=G^{\prime\prime}/\Omega$. It should be noted that their definitions are not unique in the nonlinear response regime LAOS1 ; LAOS2 ; HyunRev ; LAOS3 . In this study, we adopt the following expressions LAOS3 :

[TABLE]

where we have introduced $\gamma(t):=\gamma_{0,\textrm{eff}}\sin(\Omega t)$ and $\tilde{\sigma}$ is defined as

[TABLE]

Here $\overline{\sigma_{xy}^{\textrm{sym}}}:=\int_{0}^{\tau_{p}}dt\sigma_{xy}^{\textrm{sym}}(t)/\tau_{p}$ is the time average over the oscillatory period $\tau_{p}:=2\pi/\Omega$. 111 We have examined another protocol in which the oscillatory shear is given by $\gamma(t)=\gamma_{0,\textrm{eff}}(1-\cos(\Omega t))$, and found that the results of $G^{\prime}$ and $G^{\prime\prime}$ are unchanged for this oscillation except for the region near the bending point of $G^{\prime}$. The choice of $\gamma(t)$ is important to discuss the fragile phase OHsj ; Pradipto ; however, we are not interested in the fragile phase.

III Scaling laws for storage and loss moduli for $\mu=1$

In this section we illustrate the existence of several scaling laws for the storage and loss moduli at $\mu=1.0$. The results are basically independent of $\mu$ for $0.01\leq\mu\leq 1$ as will be shown in Sec. IV. We also demonstrate that the stress-strain curve observed in the plastic regime can be reproduced by a phenomenology.

III.1 Scaling law for the storage modulus

From our simulation, we find the following scaling form for the storage modulus (see Fig. 3):

[TABLE]

The scaling exponent $\beta_{1}^{\prime}$ which is a special case of the $\mu-$dependent exponent $\beta_{\mu}^{\prime}$ is evaluated as

[TABLE]

Note that $\beta_{\mu}^{\prime}$ is basically independent of $\mu$ (see Sec. IV).

Further, the residual storage modulus $G^{\prime}_{\textrm{res}}(\hat{P})$ in the low strain limit depends on the pressure as

[TABLE]

where $G^{\prime}_{1}=0.26\pm 0.01$ and $a=1.49\pm 0.12$ (see Fig. 4). This behavior can be interpreted as follows. It is widely known that the excess coordination number $\Delta z:=z-z_{\textrm{iso}}$ is proportional to $\sqrt{\hat{P}}$, i. e. $\Delta z\propto\hat{P}^{1/2}$, where $z$ and $z_{\textrm{iso}}$ are the coordination and the isostatic coordination numbers, respectively Ojam ; Mjam ; Djam ; Hrev . Because the shear stress is expected to be proportional to the excess coordination number in the elastic regime, we obtain $G^{\prime}_{\textrm{res}}-G^{\prime}_{\mu}\propto\hat{P}^{1/2}$, where $G_{\mu}^{\prime}$ is a constant, irrespective of the pressure.

We find that $G^{\prime}$ is independent of $\gamma_{0,\textrm{eff}}$ for small $\gamma_{0,\textrm{eff}}$, while $G^{\prime}$ is proportional to $\gamma_{0,\textrm{eff}}^{-1}$ for large $\gamma_{0,\textrm{eff}}$. These asymptotic forms in Eqs. (14)-(17) can be interpreted as follows. If the strain is small, the material behaves as the Hookean regime in which $G^{\prime}$ must be independent of $\gamma_{0,\textrm{eff}}$. On the other hand, the stress ratio $\sigma_{xy}/P$ approaches a constant $\mu_{\textrm{macro}}$ for the large strain (plastic) regime because the macroscopic dynamical friction constant $\mu_{\rm macro}$ is well-defined in granular materials, which leads to $\beta_{\mu}^{\prime}=1$. This simple description explains the crossover from the Hookean regime for small strains to the plastic regime of $G^{\prime}$ for large strains.

We also note that this argument is independent of $\mu$. Therefore, we expect that the scaling law characterized by Eqs. (14)-(16) is held for arbitrary $\mu$. We can also evaluate the bending point of $G^{\prime}$, which is the turning point from the Hookean regime to the plastic regime, as $\hat{P}/\gamma_{0,{\rm eff}}\sim 1$ by setting the macroscopic friction constant $\mu_{\rm macro}$ to be unity. Suprisingly, this crude estimation of the location of the bending point is reasonable. Indeed, the bending point appears to be located in the range $1<\gamma_{0,{\rm eff}}/\hat{P}<2$ which provides a good estimation of the bending point.

III.2 Scaling law for the loss modulus

The loss modulus also satisfies the scaling form (see Fig. 5):

[TABLE]

where the scaling exponents $\beta_{1}^{\prime\prime}$ which is a special cases of $\beta_{\mu}^{\prime\prime}$ and is evaluated as

[TABLE]

We have confirmed that the scaling exponent $\beta_{\mu}^{\prime\prime}$ is approaching independence of $\mu$ (see Sec. IV).

As with $G^{\prime}$, the loss modulus in the plastic regime satisfies $G^{\prime\prime}\sim P/\gamma_{0,\textrm{eff}}$ for large $\gamma_{0,\textrm{eff}}$ (i. e. $\beta_{\mu}^{\prime\prime}=1$ ) because the stress ratio approaches a constant in this regime. Figsures 3 and 5 suggest that the bending point of $G^{\prime}$ and the peak of $G^{\prime\prime}$ are located at $\gamma_{0,\textrm{eff}}/\hat{P}\simeq 1$ and $\gamma_{0,\textrm{eff}}/\hat{P}\simeq 4$, respectively. Although there is no reason for the location of the peak of $G^{\prime\prime}$ to be close to that of the bending point of $G^{\prime}$, it is interesting that both points take place at the points $\gamma_{0,\textrm{eff}}/\hat{P}\sim O(1)$. It is remarkable that $G^{\prime\prime}$ has three regimes: (i) $G^{\prime\prime}$ is independent of $\gamma_{0,\textrm{eff}}$ for small strain amplitudes, (ii) $G^{\prime\prime}$ is linear for $\gamma_{0,\textrm{eff}}$ for the intermediate strain amplitudes, and (iii) $G^{\prime\prime}\sim\hat{P}/\gamma_{0,\textrm{eff}}$ in the plastic regime for large $\gamma_{0,\textrm{eff}}$. The existence of the linear regime (ii) suggests that the relaxation process for elastic vibrations plays an important role in this regime (see Appendix C).

To explain the crossover of $G^{\prime\prime}$ from regime (i) for small strain ($\gamma_{0,\textrm{eff}}/\hat{P}<x_{c}$) to regime (ii) for intermediate strain $x_{c}<\gamma_{0,\textrm{eff}}/\hat{P}<1$, where we have introduced the crossover point $x_{c}$. Because the contact model of our simulation model is the Kelvin-Voight model, $G^{\prime\prime}$ in regime (i) satisfies (see Fig. 13 in Appendix B)

[TABLE]

where $\hat{\Omega}=\Omega\sqrt{m_{0}/k_{n}}$ and $\hat{\xi}_{n}=\xi_{n}/\sqrt{m_{0}k_{n}}$. On the other hand, $G^{\prime\prime}$ in the regime (ii) , as shown in Appendix C, is independent of $\Omega$ as

[TABLE]

From the balance Eqs. (22) and (23), at $x_{c}$ we estimate $x_{c}\simeq\hat{\Omega}\hat{\xi}_{n}$, which is verified from our simulation (see Fig. 14 in Appendix B).

III.3 The role of the angular distribution functions of contact forces

We have studied the behaviors of $G^{\prime}$ and $G^{\prime\prime}$ as well as the stress-strain curve in the plastic regime quantitatively. It is known that the approximate stress tensor can be expressed by the angular distributions of the normal contact force density $\zeta_{N}(\theta)$ and tangential contact force density $\zeta_{T}(\theta)$ in simple shear flows Cruz ; calten1 ; Kanatani ; calten3 . We try to apply this theory to grains under oscillatory shear. Here $\zeta_{N}(\theta)$ and $\zeta_{T}(\theta)$ are expressed as $\zeta_{N}(\theta)=\rho(\theta)F_{N}(\theta)/\langle F_{N}\rangle$ and $\zeta_{T}(\theta)=\rho(\theta)F_{T}(\theta)/\langle F_{N}\rangle$ with the angular contact distribution $\rho(\theta)$, the normal force $F_{N}(\theta)$, the tangential force $F_{T}(\theta)$ at the contact angle $\theta$, and the averaged contact force over the angle $\langle F_{\alpha}\rangle=\int_{0}^{2\pi}d\theta\rho(\theta)F_{\alpha}(\theta)$ for $\alpha=N$ or $T$. There are the normalization conditions for these functions as $\int_{0}^{2\pi}d\theta\rho(\theta)=1$, $\int_{0}^{2\pi}d\theta\zeta_{N}(\theta)=1$ and $\int_{0}^{2\pi}d\theta\zeta_{T}(\theta)=0$, respectively. Figure 6 shows the plot of $\zeta_{N}(\theta)$ and $\zeta_{T}(\theta)$ in our system for various $\gamma_{0,\textrm{eff}}$ at $\hat{P}=2.0\times 10^{-3}$. In our system, the behavior of $\zeta_{N}(\theta)$ or $\zeta_{T}(\theta)$ for large strain is similar to that observed in simple shear Cruz ; Kanatani , while $\zeta_{N}(\theta)$ is isotropic for small strain.

In our system, even under oscillatory shear, the stress tensor $\sigma_{\alpha\beta}$ for large strain can be approximately expressed as

[TABLE]

where we have used $P=-(\sigma_{xx}+\sigma_{yy})/2$. Here, $(n_{\theta,\ x},\ n_{\theta,\ y})=(\cos\theta,\ \sin\theta)$ and $(t_{\theta,\ x},\ t_{\theta,\ y})=(-\sin\theta,\ \cos\theta)$ are the normal and tangential unit vectors between contacting grains, respectively. The phenomenological stress is only applicable to the plastic regime because the phenomenological stress is proportional to the pressure.

Using Eq. (24) our remaining task is to estimate $\zeta_{N}(\theta)$ and $\zeta_{T}(\theta)$ theoretically. The shear stress that originated from the normal contacts must be proportional to $\sin(2\theta)$ because the viscosity is proportional to $\int d\bm{r}\frac{xy}{r}V_{n}^{\prime}(r)g_{\rm eq}(r)h(\bm{r})$ where $g_{\rm eq}(r)$ is the equilibrium radial distribution function and $h(\bm{r})$ is the shear contribution of the radial distribution function $g(\bm{r})=g_{\rm eq}(r)(1+\dot{\gamma}h(\bm{r}))$. Because of the orthogonality of the trigonometric functions, the angular distribution that contributes to the shear stress is proportional to $\sin(2\theta)$, while the equilibrium contribution should be independent of $\theta$. Similarly, the viscosity that originated from the tangential contacts contains $\cos^{2}\theta$ for the positive tangential contacts and $\sin^{2}\theta$ for the negative tangential contacts. As shown in Fig. 6, our simple description can recover the angular distributions of contacts with the introduction of fitting parameters $b$ and $c$. Further, $b$ and $c$ depend on $\gamma_{0,\textrm{eff}}$ and $\hat{P}$ ( see Appendix D). We then obtain the theoretical shear stress:

[TABLE]

The stress-strain (Lissajous) curve for large strain amplitudes is plotted in Fig. 7, where the stress from our simulation and the phenomenology from Eq. (25) are expressed by the symbols and the lines, respectively. It is remarkable that the phenomenological evaluation of the stress $\sigma$ from Eqs. (24) and (25) recovers the stress-strain curve for $\gamma_{0,\textrm{eff}}/\hat{P}\geq 1$, at least.

Because we can reproduce the stress-strain curve via phenomenology with the aid of $\zeta_{N}(\theta)$ and $\zeta_{T}(\theta)$, we can also reproduce $G^{\prime}$ and $G^{\prime\prime}$ via the phenomenological storage and loss moduli with the aid of Eqs. (11) and (12). The lines in Figs. 3(b) and 5(b) show the phenomenological $G^{\prime}$ and $G^{\prime\prime}$, respectively. Although the phenomenology cannot be used for small strain amplitude, we reproduce the quantitative behavior of $G^{\prime}$ and $G^{\prime\prime}$ for $\gamma_{0,\textrm{eff}}\gtrsim\hat{P}$.

IV Scaling laws for various friction coefficients

We have discussed the scaling laws only for $\mu=1.0$ in the previous section. In this section we discuss how the results of Eqs. (14)-(20) depend on $\mu$. In short, (i) the scaling laws in Eqs. (14)-(20) still hold for arbitrary $\mu$, and (ii) the exponents $\beta_{\mu}^{\prime}$ and $\beta_{\mu}^{\prime\prime}$ are independent of $\mu$.

For arbitrary $\mu$ the storage modulus satisfies the relation:

[TABLE]

and the loss modulus satisfies the relation:

[TABLE]

The estimated scaling exponents from our simulation are plotted in Fig. 8. The exponents $\beta^{\prime}_{\mu}$ and $\beta^{\prime\prime}_{\mu}$ are almost independent of $\mu$, which are approximately equal to unity.

The scaling plots of $G^{\prime}$ and $G^{\prime\prime}$ against the scaled the strain amplitude at $\mu=0.1$ are presented in Fig. 9. These scaling plots suggest that the scaling discussed in Sec. III can be used for arbitrary $\mu$.

To confirm our conjecture that the scaling is independent of $\mu$ we examine the case of $\mu=0.01$. It is interesting that the scaling laws seem to be valid even for $\mu=0.01$ without any changes of scaling exponents, though the distinction between the plastic and elastic regimes is not sharp in this case (see Fig. 10). Nevertheless, we should emphasize that the results of the small friction coefficient limit for frictional grains is completely different from those in frictionless cases (see Appendix E) in which the scaling laws do not exist as presented in Ref. BSosc . This singularity in the zero-friction limit is consistent with that observed in Ref. OHdst2 .

V discussion and conclusion

First, we can observe the dilatancy during the process because the volume can be changed. It is known that the packing fraction $\delta\phi:=\phi-\phi_{0}$ from the unbiased fraction $\phi_{0}$ satisfies $\delta\phi\propto-\dot{\gamma}\sqrt{m/P}$ in a two-dimensional pressure control system under steady simple shearCruz ; Fdil . It is also known that the volume is compactified if the strain is small Nowak98 . The detailed behavior of the dilatancy will be reported elsewhere.

Second, the origin of the macroscopic friction law, i. e. $\sigma/P$ tends to be independent of the strain for large strain limit should be clarified. This is equivalent to the microscopic determination of $b$ and $c$ in Eq. (25). Quantitative analysis to determine $b$ and $c$ will be considered elsewhere.

In summary, we have numerically studied a frictional granular system confined by a constant pressure under oscillatory shear. We confirmed the existence of scaling laws for the storage and loss moduli. We also verified the absence of the distinction between the softening and yielding regimes for the frictional system. These results are completely different from those for frictionless systems; however, the friction coefficient can be disregarded, particularly for $0.01\leq\mu\leq 1$. We found the crossover of the storage modulus from $G^{\prime}/G^{\prime}_{\textrm{res}}=1$ for small strain amplitude to $G^{\prime}/G^{\prime}_{\textrm{res}}\propto\hat{P}/\gamma_{0,\textrm{eff}}$ for large strain amplitude. We also found that the scaling of the loss modulus which has three regimes, in which (i) $G^{\prime\prime}$ is independent of the strain amplitudes for small strain amplitudes, (ii) $G^{\prime\prime}$ is a linear function of $\gamma_{0,\textrm{eff}}$ for intermediate strain amplitudes, and (iii) $G^{\prime\prime}\sim\hat{P}/\gamma_{0,\textrm{eff}}$ for large strain amplitudes. The phenomenological theory with the aid of the angular distributions of the contacts works in the plastic regime

Acknowledgement

One of the authors (DI) expresses his sincere gratitude to T. G. Sano, S. Takada, and K. Saitoh for their assistance with the granular simulation. The authors thank M. Otsuki, T. Kawasaki, R. Seto, and N. Oyama for fruitful discussions and usefull comments. 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.

Appendix A Coupled stress and asymmetric stress

In this appendix, we discuss the coupled stress $R$ which is the asymmetric part of the stress tensor Gol defined as

[TABLE]

where $T_{ij}:=x_{ij}f_{ij,y}-y_{ij}f_{ij,x}$. Here, $f_{ij,k}$ with $k=x$ or $k=y$, $x_{ij}$, and $y_{ij}$ are the $k-$component of $\bm{f}_{ij}$, $x-$component of $\bm{r}_{ij}$, and $y-$component of $\bm{r}_{ij}$, respectively. Let us introduce a normalized coupled stress $\hat{R}$:

[TABLE]

where $\tilde{R}:=R-\int_{0}^{\tau_{p}}dtR/\tau_{p}$. We have confirmed $\hat{R}\ll 1$ for high pressure as shown in Fig. 11, though $\hat{R}(t)$ exhibits a stress-strain curve similar to the symmetric stress in the main text. Although $\hat{R}(t)$ for large strain amplitudes is visible, the largest value of $\hat{R}$ is 0.07 in Fig. 12 (a). Therefore, we can safely ignore the contribution of the asymmetric stress tensor or the couple stress.

Appendix B Angular frequency dependence of $G^{\prime}$ and $G^{\prime\prime}$

In this appendix, we investigate how the storage and loss moduli depend on $\Omega$ by fixing $\hat{P}=2.0\times 10^{-3}$. We have confirmed that $G^{\prime}$ is almost independent of $\Omega$ for $\Omega\sqrt{m_{0}/k_{n}}/(2\pi)\leq 10^{-3}$ in Fig. 13 (a), while $G^{\prime\prime}$ strongly depends on $\Omega$ for small $\gamma_{0,\textrm{eff}}/\hat{P}$ in Fig. 13 (b).

Here, we clarify how $G^{\prime\prime}$ depends on $\Omega$ for the small strain regime ($\gamma_{0,\rm{eff}}/\hat{P}<x_{c}$). We have introduced a critical parameter $x_{c}$ as the middle point of $\gamma_{0,\textrm{eff}}/\hat{P}$ which satisfies

[TABLE]

Figure 14 (a) shows that $x_{c}$ is also proportional to $\Omega$.

We have also introduced $G^{\prime\prime}_{0}$ by

[TABLE]

From Fig. 14(b), we have confirmed that $G^{\prime\prime}_{0}\propto\Omega$, as mentioned in Sec.III.2.

Appendix C The loss modulus in the linear regime to $\gamma_{0,\rm{eff}}$

In this appendix, we explain the origin of the regime (ii) in which $G^{\prime\prime}$ is proportional to $\gamma_{0,\rm{eff}}/\hat{P}$. We also evaluate the value of $G^{\prime\prime}$ at the peak.

The behavior of $G^{\prime\prime}$ for $\gamma_{0,\textrm{eff}}/\hat{P}\lesssim 1$ can be interpreted as follows. We will focus on the viscous shear stress $\sigma^{\textrm{(vis)}}$ for this purpose. For simplicity, we ignore the difference between the peak point of $G^{\prime\prime}$ and the bending point of $G^{\prime}$, although these two points are slightly different in reality.

Here the shear stress $\sigma^{\textrm{(vis)}}_{<}$ in the regime (ii) satisfying $\gamma_{0,\rm{eff}}\lesssim\hat{P}$ might obey

[TABLE]

where $\tau$ is a typical time scale for elastic energy relaxation, and the subscript $<$ represents the regime $\gamma_{0,\textrm{eff}}/\hat{P}\lesssim 1$. We can assume that the function $f(\hat{P})$ obeys a power law as $f(\hat{P})\sim\hat{P}^{\Delta}$.

The estimation of the time scale $\tau$ is as follows. The strain energy $E$ for a grain with the stiffness $k_{n}$ under the strain $\gamma$ might be expressed as $E\sim k_{n}(\gamma\ell)^{2}$, where $\ell$ stands for a typical length scale. This restoring energy should be balanced with the kinetic energy $E\sim m(\ell/\tau)^{2}$ where $m$ and $\tau$ are the mass of the grain and the typical time scale to characterize the time differentiation. Therefore, we obtain $\tau\sim\sqrt{m_{0}/k_{n}}/\gamma$. In our case, we can replace $\gamma$ by $\gamma_{0,\rm{eff}}$. We then obtain the form

[TABLE]

For the plastic regime $\gamma_{0,\textrm{eff}}/\hat{P}\gtrsim 1$, on the other hand, the stress $\sigma^{\textrm{(vis)}}_{>}$ ($>$ stands for the regime $\gamma_{0,\textrm{eff}}/\hat{P}\gtrsim 1$) satisfies

[TABLE]

The stresses in Eqs. (36) and (37) must take the identical value as

[TABLE]

at the peak of $G^{\prime\prime}$ under the condition $\gamma_{0,\rm{eff}}\sim\hat{P}$. Therefore, we obtain the relation $\hat{P}^{\Delta+1}\sim 1$. Because the right hand side of this equation is independent of $\hat{P}$, we get $\Delta=-1$, and $\sigma^{\rm{(vis)}}\sim k_{n}\gamma_{0,\rm{eff}}^{2}/\hat{P}\sim P$ at the peak. Finally, we obtain $G^{\prime\prime}_{\textrm{max}}=k_{n}$. This indicates that the peak value of the loss modulus should be independent of $P$.

Appendix D Fitting parameters for $\zeta_{N}$ and $\zeta_{T}$

In this appendix, we present how the fitting parameters $b$ and $c$ introduced in Fig. 6 and Eq. (25) depend on $\hat{P}$ and $\gamma_{0,\textrm{eff}}$. Here, we numerically evaluate the strain dependence of $b$ and $c$ at $\hat{P}=2.0\times 10^{-3}$ in Fig. 15. We also plot the pressure dependence of $b$ and $c$ at $\gamma_{0,\textrm{eff}}=1.0$ in Fig. 16.

Appendix E Results for frictionless grains

In this appendix, we present the results of our simulation for the storage and loss moduli for a frictionless system. It is known that the rigidity of two-dimensional frictionless systems near the jamming point is proportional to $\hat{P}^{1/2}$ Djam ; Ojam ; Otsuki14 . Because of $P\propto\phi-\phi_{J}$ near the jamming point, $G^{\prime}/(\phi-\phi_{J})^{1/2}$ corresponds to $G^{\prime}/\hat{P}^{1/2}$, where $\phi$ and $\phi_{J}$ are the volume fraction and the jamming fraction, respectively. We expect that the scaling between $G^{\prime}/\hat{P}^{1/2}$ and $\gamma_{0.\rm{eff}}$ is held corresponding to Ref. Otsuki14 . Here, we plot $G^{\prime}/\hat{P}^{1/2}$ against $\gamma_{0,\textrm{eff}}/\hat{P}$ in Fig. 17. However, there is no universal scaling law for frictionless systems; this is significantly different from the results of finite $\mu$ (even for $\mu=0.01$) which exhibit definite scaling laws. The scaling corresponding to Ref. Otsuki14 is only visible for a weakly plastic regime under low pressure. Figure 17 suggests the existence of the softening region when $\gamma_{0,\textrm{eff}}/\hat{P}\simeq 1$ for low pressure cases. On the other hand, we find the yielding for $\gamma_{0,\textrm{eff}}/\hat{P}>10^{3}$. Notably, the existence of two distinct regions corresponds to the softening and yielding described by Ref. BSosc , which contains the background friction.

We plot $G^{\prime\prime}$ against $\gamma_{0,\textrm{eff}}/\hat{P}$ in Fig. 18. The scaling for $G^{\prime\prime}$ is also absent, which is completely different from that of $\mu=0.01$.