## 2 The Microstructure Model Considering Nonaffine Deformations
### Local and Global Constitutive Relations  
At particle level, for the elastic regime considered in this paper, the interaction between particles in contact can be regarded as linearly elastic. Hence, ignoring the effect of cross terms, the local constitutive relations for particle pairs (i.e., contact law) can be given as:  
\[\mathbf{\dot{f}}=\mathbf{K}^{\zeta}\cdot\mathbf{\Delta}^{\zeta\zeta},\mathbf{m}^{\zeta}= \mathbf{G}^{\zeta}\cdot\mathbf{\theta}^{\zeta},\ \zeta=\text{a, r, g}, \tag{12}\]  
where \(\mathbf{\dot{f}}\) is the contact force and \(\mathbf{m}^{\zeta}\) is the contact moment. The linear stiffness \(\mathbf{K}^{\zeta}\) and rotation stiffness \(\mathbf{G}^{\zeta}\) can be expressed as [19; 20]  
\[\mathbf{K}^{\zeta}=\Big{(}K_{\text{n}}^{\zeta}-K_{\text{s}}^{\zeta}\Big{)}\mathbf{n }\otimes\mathbf{n}+K_{\text{s}}^{\zeta}\mathbf{I}, \tag{13}\]  
\[\mathbf{G}^{\zeta}=\Big{(}G_{\text{n}}^{\zeta}-G_{\text{s}}^{\zeta}\Big{)}\mathbf{ n}\otimes\mathbf{n}+K_{\text{s}}^{\zeta}\mathbf{I}, \tag{14}\]  
where \(K_{\text{n}}^{\zeta}\) and \(K_{\text{s}}^{\zeta}\) are linear displacement stiffnesses at the contact in normal and tangential directions, respectively. \(G_{\text{n}}^{\zeta}\) and \(G_{\text{s}}^{\zeta}\) are rotation stiffnesses in normal and tangential directions, respectively. \(\mathbf{n}\) is the contact normal vector, which can be expressed as \(\mathbf{n}=\cos\gamma\mathbf{x}_{1}+\sin\gamma\cos\varphi\mathbf{x}_{2}+\sin\gamma\sin \varphi\mathbf{x}_{3}\) in global coordinate systems (see Figure 2), and \(\mathbf{I}\) is the second order identity tensor. Consequently, ignoring the effect of cross terms, the microscale deformation energy increment \(\text{d}w^{\epsilon}\) at a contact can be reasonably defined as follows:  
\[\text{d}w^{\epsilon}=\sum_{\zeta}\Big{(}\mathbf{K}^{\zeta}\cdot\mathbf{\Delta}^{ \epsilon\zeta}\cdot\text{d}\mathbf{\Delta}^{\epsilon\zeta}+\mathbf{G}^{\zeta} \cdot\mathbf{\theta}^{\epsilon\zeta}\cdot\text{d}\mathbf{\theta}^{\epsilon\zeta}\Big{)}. \tag{15}\]  
At continuum level, based on the macroscale deformations analysis previously, the macroscale deformation energy density increment \(\text{d}W\) can be expressed as a function of the macroscale continuum kinematic measures, i.e., \(\text{d}W=\text{d}W(\epsilon,\mathbf{R},\
abla\mathbf{A},\kappa,\mathbf{D}, \
abla\mathbf{B})\). The following form is adopted in this paper, ignoring the effect of cross terms:  
\[\text{d}W=\mathbf{\sigma}:\text{d}\epsilon+\mathbf{\tau}:\text{d}\mathbf{R}+\mathbb{T }^{3}\vdots\text{d}(\
abla\mathbf{A})+\mathbf{\mu}:\text{d}\kappa+\mathbf{v}:\text{d} \mathbf{D}+\mathbb{V}^{3}\vdots\text{d}(\
abla\mathbf{B}). \tag{16}\]  
Based on energy conservation, the macroscale deformation energy density should be equal to the volumetric average of the sum of the microscale deformation energy over all contacts for the elastic conditions where the energy dissipation is not considered. Meanwhile, for a material point with a very large number of particles, it is reasonable to use an integral form to replace the summation. Thus, we have  
\[\text{d}W=\frac{N}{V}\iint\limits_{V}\sum_{\zeta}\Big{(}\mathbf{K}^{\zeta} \cdot\mathbf{\Delta}^{\epsilon\zeta}\cdot\text{d}\mathbf{\Delta}^{\epsilon\zeta}+ \mathbf{G}^{\zeta}\cdot\mathbf{\theta}^{\epsilon\zeta}\cdot\text{d}\mathbf{\theta}^{ \epsilon\zeta}\Big{)}\ \xi(\gamma,\ \varphi)\sin\gamma\text{d}\gamma\text{d}\varphi, \tag{17}\]  
Figure 2: Relations of local and global coordinate systems. in which \(N\) is the total number of contacts and \(\xi(\gamma,\ \varphi)\) denotes the probability density function of contact distribution. For a transversely isotropic granular media, \(\xi(\gamma,\ \varphi)\) can be approximated by the second-order spherical harmonic function, given in the tensor form by  
\[\xi(\gamma,\ \varphi)=\frac{1}{4\pi}(\updelta+\mathbf{F}):(\mathbf{n}\otimes\mathbf{n}), \tag{18}\]  
where \(\updelta\) is the Kronecker tensor, and \(\mathbf{F}\) denotes the deviatoric tensor of contact distribution, given by  
\[\mathbf{F}=