1 Introduction

A granular system is an assembly of large size particles in which the thermal fluctuation is irrelevant. The physics of granular matter has attracted much attention due to many fascinating effects presented in granular matter and its ubiquity in nature and everyday life [1, 2]. From the theoretical point of view, a wide range of complex behaviors is far from fully understood for granular matter. Many unusual behaviors, such as the stress response of granular solids caused by a point load and the Janssen effect, are observed in static or quasi-static granular matter [3,4,5,6]. Many of these behaviors are related to the physical laws of stress distribution in granular medium. In the conventional elasticity theory, the stress distribution in an elasticity solid can be calculated by solving its governing equations with some proper boundary conditions [7]. However, the constitutive relation still hasn’t been obtained to complete the governing equations for granular solids except for some special cases. One way for the research of stress distribution is to consider a granular solid as an elastic medium. Based on some experimental facts and mechanical principles, it is possible to propose a proper strain–stress constitutive relation to close the governing equations [8].It is also meaningful to investigate the granular medium from different length scales. From the microscopic perspective, a static granular system is a particulate packing in mechanical equilibrium and it is straightforward to obtain the equations for intergranular forces according to classical mechanics. A granular packing can also be modeled as a continuous medium from the large length scale and the concept of stress tensor is more convenient for the description of mechanical behaviors of granular medium. The stress tensor on a grain can be defined as the density of force moment tensor on this grain and the stress tensor field in granular medium can be calculated by using micro-level variables such as contact forces and contact vectors [9, 10]. However, one problem which hasn’t been completely solved is how to derive the constitutive equation based on micro-mechanical analysis of the packings.

The stress field of a dry granular assembly due to external loads may be completely determined without knowing the details of the deformation of grains. A typical case is the isostatic granular system whose intergranular forces can be determined by some known boundary forces and the balance equations of intergranular forces and torque moments.To the isostatic granular packing, its average contact number has to be a particular value \(Z_{C}\). Many theoretical methods were proposed to derive the governing equations of stress fields in isostatic granular packings. An interesting method is based on Edwards statistical mechanics of jammed granular matter [10,11,12]. The results of this method show that the well-known equations of stress tensor \({\mathbf {Q}}\), i.e. \({\mathbf {Q}}={\mathbf {Q}}^{T}\) and \(\nabla \cdot {\mathbf {Q}}+\rho {\mathbf {g}}={\mathbf {0}}\) root in the balance equations of torque moments and intergranular forces respectively. By using the Edwards statistical method, the constitutive relations of stress tensor can also be derived for some simple cases. The most significant consequence is that these constitutive relations appear in the form of stress-geometry equations instead of stress-strain equations. Since three dimensional packings are complicated, many theoretical works focus on the governing equations of two dimensional granular packings [13,14,15,16]. For packings without body force, gauge variables such as loop forces were introduced to simplify the theoretical analysis on balance equations and the derivation for constitutive equations. Blumenfeld networks [13] and Voronoi networks [15] were used to describe the geometrical structure of granular packings and some fabric tensors were proposed to characterize the local geometry around a grain. The stress-geometry equations which were obtained in these works [13, 16] provide the missing link between the stress tensor and the local geometry.

After the governing equations were obtained, the stress field in granular medium can be solved when the proper boundary conditions are provided. The Blumenfeld stress-geometry equation was used to close the governing equations and the solutions were used to investigate stress chains in granular medium [17, 18]. The existed theoretical researches show that the packing fabric tensor plays an essential role to the stress transmission in isostatic granular packings. However, the values of geometrical tensors usually fluctuate around zero mean and this means the coarse-graining of stress-geometry equation far from straightforward. Granular packings are not necessarily isostatic, therefore the theories for obtaining governing equations should be extended to hyperstatic systems in which the equilibrium equations are insufficient for the determination of intergranular forces. Many of the theoretical results should also be verified by performing experiments or computer simulations. Recently, we proposed a multiscale method to derive the governing equations for rhombic lattice packings of disks with friction and obtained a constitutive equation with a coefficient called friction-geometry tensor [19]. The Janssen coefficient \(\lambda \) was also given by using the eigenvalues of the friction-geometry tensor. Based on these works, it is natural to use the governing equations to investigate the stress transmission in rhombic disk packings with some typical boundary conditions.

In this paper, the governing equations for the components of stress tensor \({\mathbf {Q}}\) are given and the hyperbolic nature of these equations shows that the stresses transmit through rhombic disk packings in the behavior of traveling wave propagation. We propose a method which is called simultaneous diagonalization to decouple the stress equations. Then, we use the method of characteristic curves to derive the exact solutions for stress transmission in rhombic disk packings with friction. The theoretical results show that the square root of Janssen coefficient \(\sqrt{\lambda }\) plays a role of wave speed for the stress transmission. To verify the theoretical results, DEM simulations are performed to investigate the Janssen effect and the stress response under a point load for the rhombic disk packings. We discuss the differences between the DEM simulations and the theoretical results in detail.

2 Theoretical model and solutions

Fig. 1
figure 1

A two dimensional granular medium formed by a rhombic lattice packing of disks. A local area of the rhombic packing is picked out and presented on the left side of this figure in which \({\mathbf {k}}_{1},{\mathbf {k}}_{2}\) are the unit lattice vectors of this packing. \(\theta \) is the angle between \({\mathbf {k}}_{1}\) and \({\mathbf {k}}_{2}\). \(\alpha \) is the angle between \({\mathbf {k}}_{1}\) and the unit vector \({\mathbf {e}}_{x}\). \({\mathbf {g}}\) is the gravitational acceleration

Full size image

A granular packing can be modeled as a continuous medium from the large length scale. The governing equations of stress distribution in granular media are essential for understanding the macroscopic behaviors of granular matter. In this paper, we study the stress distribution in rhombic lattice packings of disks with friction. As shown in Fig.1, a rhombic lattice packing is formed by the array of disks with radius \(r_{0}\) and mass \(m_{0}\). Each disk is only subject to the gravity \(m_{0}{\mathbf {g}}\) and contact forces. \(\{{\mathbf {e}}_{x},{\mathbf {e}}_{y}\}\) is the orthonormal basis of the Cartesian coordinate. \(\{{\mathbf {k}}_{1},{\mathbf {k}}_{2}\}\) are the unit lattice vectors. The angle \(\alpha \) indicates the orientation of the rhombic lattice with respect to the Cartesian coordinate and the angle \(\theta \) represents the structure of this lattice. Let \(\rho =\dfrac{m_{0}}{4r_{0}^{2}sin\theta }\) be the mass density of this packing and \(\mu \) be the friction coefficient of contacts between disks. For each contact, the friction angle \(\beta \) is defined by \(\mathbf {\tan \beta }=\dfrac{f_{\tau }}{f_{N}}=\mu \mathbf {\eta }\) where \(f_{\tau }\) and \(f_{N}\) are the tangential and normal component of the contact force respectively and \(\eta \) is called friction mobilization which satisfies \(\mathbf {\mathrm {-1}\le \eta }\le 1\). We suppose that the values of friction mobilization of all contacts between disks approximate a constant, so do the values of friction angle. The governing equations of stress tensor field \({\mathbf {Q}}\) in the granular medium in Fig.1 can be derived [19] and we have

$$\begin{aligned} {\left\{ \begin{array}{ll} \nabla \cdot {\mathbf {Q}}+\rho {\mathbf {g}}={\mathbf {0}}\\ {\mathbf {Q}}={\mathbf {Q}}^{T}\\ \mathbf {Q:G=0} \end{array}\right. } \end{aligned}$$
(1)

where \({\mathbf {G}}\) is called friction-geometry tensor and is defined as

$$\begin{aligned} {\mathbf {G}}={\mathbf {k}}^{1}{\mathbf {k}}^{2}+{\mathbf {k}}^{2}{\mathbf {k}}^{1}-\dfrac{\sin \beta }{\sin (\theta -\beta )}({\mathbf {k}}^{1}{\mathbf {k}}^{1}+{\mathbf {k}}^{2}{\mathbf {k}}^{2}) \end{aligned}$$
(2)

In this formula, \(\{{\mathbf {k}}^{1},{\mathbf {k}}^{2}\}\) are the dual unit lattice vectors of \(\{{\mathbf {k}}_{1},{\mathbf {k}}_{2}\}\); in other words, they satisfy \({\mathbf {k}}^{i}\cdot {\mathbf {k}}_{j}=\sin \theta \delta _{ij}\) where \(\delta _{ij}\) is the Kronecker delta. The third equation in Eq.(1) is the constitutive equation for rhombic disk packings and the detailed derivation process of this equation can be found in the article ( [19]). Since \({\mathbf {G}}\) is a symmetric tensor, if we choose the basis \(\{{\mathbf {e}}_{x},{\mathbf {e}}_{y}\}\) to make \(\alpha =\dfrac{\pi }{2}-\dfrac{\theta }{2}\), \({\mathbf {G}}\) can be written in its diagonal form as

$$\begin{aligned} {\mathbf {G}}=G_{xx}{\mathbf {e}}_{x}{\mathbf {e}}_{x}+G_{yy}{\mathbf {e}}_{y}{\mathbf {e}}_{y} \end{aligned}$$
(3)

where \(\{{\mathbf {e}}_{x},{\mathbf {e}}_{y}\}\) are the eigenvectors of the friction-geometry tensor \({\mathbf {G}}\) and \(\left\{ G_{xx},G_{yy}\right\} \) are its eigenvalues. Suppose that \({\mathbf {Q}}=Q_{xx}{\mathbf {e}}_{x}{\mathbf {e}}_{x}+Q_{xy}{\mathbf {e}}_{x}{\mathbf {e}}_{y}+Q_{yx}{\mathbf {e}}_{y}{\mathbf {e}}_{x}+Q_{yy}{\mathbf {e}}_{y}{\mathbf {e}}_{y}\). Then, the equation \({\mathbf {Q}}={\mathbf {Q}}^{T}\) requires that \(Q_{xy}=Q_{yx}\). Using these notations, the constitutive equation in Eq.(1) can be written as

$$\begin{aligned} Q_{xx}=\lambda Q_{yy} \end{aligned}$$
(4)

where \(\lambda \) is called Janssen coefficient and we have

$$\begin{aligned} \lambda =\dfrac{-G_{yy}}{G_{xx}}=\tan \dfrac{\theta }{2}\tan (\dfrac{\theta }{2}-\beta ) \end{aligned}$$
(5)

It will be shown below that \(\lambda \) is an important parameter for the macroscopic mechanical properties of rhombic disk packings. Based on the above discussion, the equations for stress components \(Q_{yy}\) and \(Q_{xy}\) can be derived as

$$\begin{aligned} {\left\{ \begin{array}{ll} \lambda \frac{\partial Q_{yy}}{\partial x}+\frac{\partial Q_{xy}}{\partial y}=-\rho g_{x}\\ \frac{\partial Q_{yy}}{\partial y}+\frac{\partial Q_{xy}}{\partial x}=-\rho g_{y} \end{array}\right. } \end{aligned}$$
(6)

where \(g_{x}={\mathbf {g}}\cdot {\mathbf {e}}_{x}\) and \(g_{y}={\mathbf {g}}\cdot {\mathbf {e}}_{y}\) are the components of the gravity acceleration.

For convenience, we introduce differential operators \(D_{x}=\frac{\partial }{\partial x}\) and \(D_{y}=\frac{\partial }{\partial y}\) . Equation (6) can be written as

$$\begin{aligned} \left( \begin{array}{cc} \lambda D_{x} &{} D_{y}\\ D_{y} &{} D_{x} \end{array}\right) \left( \begin{array}{c} Q_{yy}\\ Q_{xy} \end{array}\right) =-\rho \left( \begin{array}{c} g_{x}\\ g_{y} \end{array}\right) \end{aligned}$$
(7)

It’s easy to see

$$\begin{aligned} \left( \begin{array}{cc} \lambda D_{x} &{} D_{y}\\ D_{y} &{} D_{x} \end{array}\right) =\left( \begin{array}{cc} \lambda &{} 0\\ 0 &{} 1 \end{array}\right) D_{x}+\left( \begin{array}{cc} 0 &{} 1\\ 1 &{} 0 \end{array}\right) D_{y}=A_{x}D_{x}+A_{y}D_{y} \end{aligned}$$
(8)

where\(A_{x}\) and \(A_{y}\) denote respectively the two real matrices in this formula. The coupling of \(Q_{yy}\) and \(Q_{xy}\) makes the difficulty to solve Eq.(7). In order to decouple \(Q_{yy}\) and \(Q_{xy}\), we suppose that there exist two invertible 2-by-2 matrices B and C which could realize the simultaneous diagonalization for \(A_{x}\) and \(A_{y}\). That is to say

$$\begin{aligned} BA_{x}C=\left( \begin{array}{cc} \eta _{x1} &{} 0\\ 0 &{} \eta _{x2} \end{array}\right) ,BA_{y}C=\left( \begin{array}{cc} \eta _{y1} &{} 0\\ 0 &{} \eta _{y2} \end{array}\right) \end{aligned}$$
(9)

By eliminating matrix B , one finds

$$\begin{aligned} C^{-1}A_{x}^{-1}A_{y}C=\left( \begin{array}{cc} \eta _{y1}/\eta _{x1}\; &{} 0\\ 0 &{} \eta _{y2}/\eta _{x2}\; \end{array}\right) \end{aligned}$$
(10)

After the eigenvalues and eigenvectors of the matrix \(A_{x}^{-1}A_{y}=\left( \begin{array}{cc} 0 &{} 1/\lambda \;\\ 1 &{} 0 \end{array}\right) \) been calculated out, the solutions for Eq.(9) can be given. We have

$$\begin{aligned} \left\{ \begin{array}{c} \eta _{y1}=\frac{1}{\sqrt{\lambda }}\eta _{x1},\eta _{y2}=\frac{-1}{\sqrt{\lambda }}\eta _{x2}\\ C^{-1}=\frac{1}{2\sqrt{\lambda }}\left( \begin{array}{cc} \sqrt{\lambda } &{} 1\\ \sqrt{\lambda } &{} -1 \end{array}\right) ,B=\frac{1}{2\lambda }\left( \begin{array}{cc} \eta _{x1} &{} \sqrt{\lambda }\eta _{x1}\\ \eta _{x2} &{} -\sqrt{\lambda }\eta _{x2} \end{array}\right) \end{array}\right. \end{aligned}$$
(11)

where \(\eta _{x1}\) and \(\eta _{x2}\) are arbitrary nonzero real numbers. Since B and C are invertible, Eq.(7) is equivalent to the following equation

$$\begin{aligned} B(A_{x}D_{x}+A_{y}D_{y})C\cdot C^{-1}\left( \begin{array}{c} Q_{yy}\\ Q_{xy} \end{array}\right) =-\rho B\left( \begin{array}{c} g_{x}\\ g_{y} \end{array}\right) \end{aligned}$$
(12)

Let \(\left( \begin{array}{c} f_{1}\\ f_{2} \end{array}\right) =C^{-1}\left( \begin{array}{c} Q_{yy}\\ Q_{xy} \end{array}\right) \),\(B'=\frac{1}{2\lambda }\left( \begin{array}{cc} 1 &{} \sqrt{\lambda }\\ 1 &{} -\sqrt{\lambda } \end{array}\right) \),\(\left( \begin{array}{c} g_{1}\\ g_{2} \end{array}\right) =B'\left( \begin{array}{c} g_{x}\\ g_{y} \end{array}\right) \) . The above equation can be written as

$$\begin{aligned} \left\{ \begin{array}{c} D_{x}f_{1}+(\sqrt{\lambda })^{-1}D_{y}f_{1}=-\rho g_{1}\\ D_{x}f_{2}-(\sqrt{\lambda })^{-1}D_{y}f_{2}=-\rho g_{2} \end{array}\right. \end{aligned}$$
(13)

These two first order partial differential equations are unidirectional wave equations with sources and can be worked out separately by using the method of characteristic curves. As an example, we elaborate the solution for the equation of \(f_{1}(x,y)\). Considering a curve \(\Gamma _{1}\) which is on the x-y plane and has the parametric form (x(s), y(s)) , this curve satisfies

$$\begin{aligned} \left\{ \begin{array}{c} \frac{dx(s)}{ds}=1,\frac{dy(s)}{ds}=(\sqrt{\lambda })^{-1}\\ x(s_{0})=x_{0},y(s_{0})=y_{0} \end{array}\right. \end{aligned}$$
(14)

where \((x_{0},y_{0})\) is a point on a boundary of the solution region of Eq.(13). Apparently, the parametric form of \(\Gamma _{1}\) is

$$\begin{aligned} \left\{ \begin{array}{c} x(s)=x_{0}+(s-s_{0})\\ y(s)=y_{0}+(\sqrt{\lambda })^{-1}(s-s_{0}) \end{array}\right. \end{aligned}$$
(15)

On the curve \(\Gamma _{1}\) , the rate of change of \(f_{1}(x(s),y(s))\) satisfies

$$\begin{aligned} \frac{df_{1}}{ds}=\frac{dx}{ds}D_{x}f_{1}+\frac{dy}{ds}D_{y}f_{1}=-\rho g_{1} \end{aligned}$$
(16)

Hence, we have

$$\begin{aligned} f_{1}(x(s),y(s))=f_{1}(x(s_{0}),y(s_{0}))-\rho g_{1}(s-s_{0}) \end{aligned}$$
(17)

Using Eqs.(15) and (17) and eliminating s, the characteristic curve \(\Gamma _{1}\) and the value of \(f_{1}(x,y)\) on \(\Gamma _{1}\) can be derived as

$$\begin{aligned} \left\{ \begin{array}{c} x=x_{0}+\sqrt{\lambda }(y-y_{0})\\ f_{1}(x_{0}+\sqrt{\lambda }(y-y_{0}),y)=f_{1}(x_{0},y_{0})-\rho g_{1}\sqrt{\lambda }(y-y_{0}) \end{array}\right. \end{aligned}$$
(18)

Similarly, the characteristic curve and the value of \(f_{2}(x,y)\) on this curve can be given by

$$\begin{aligned} \left\{ \begin{array}{c} x=x_{0}-\sqrt{\lambda }(y-y_{0})\\ f_{2}(x_{0}-\sqrt{\lambda }(y-y_{0}),y)=f_{2}(x_{0},y_{0})+\rho g_{2}\sqrt{\lambda }(y-y_{0}) \end{array}\right. \end{aligned}$$
(19)

If we regard y as time “ t ”, \(f_{1}(x,y)\) and \(f_{2}(x,y)\) are unidirectional wave functions with velocity \(\sqrt{\lambda }\). These two unidirectional waves travel along “ \(+x\) ” direction and “ \(-x\) ” direction respectively.

After the values of \(f_{1}(x,y)\) and \(f_{2}(x,y)\) been calculated out by using the known boundary conditions and the above two formulae, one can obtain the components of the stress tensor \({\mathbf {Q}}\) according to

$$\begin{aligned} \left( \begin{array}{c} Q_{yy}\\ Q_{xy} \end{array}\right) =\left( \begin{array}{cc} 1 &{} 1\\ \sqrt{\lambda } &{} -\sqrt{\lambda } \end{array}\right) \left( \begin{array}{c} f_{1}\\ f_{2} \end{array}\right) ,Q_{xx}=\lambda Q_{yy} \end{aligned}$$
(20)

Apparently, each of the component of \({\mathbf {Q}}\) is a wave function since \(f_{1}\) and \(f_{2}\) are unidirectional wave functions. In fact, by eliminating \(Q_{xy}\) in Eq.(6), the partial differential equation of \(Q_{yy}\) can be derived as

$$\begin{aligned} \lambda D_{x}^{2}Q_{yy}-D_{y}^{2}Q_{yy}=0 \end{aligned}$$
(21)

This is the wave equation in its standard form and the similar equations are also hold for \(Q_{xy}\) and \(Q_{xx}\). The above theoretical results indicate that the stresses transmit in the rhombic disk packings in the behavior of traveling wave with sources.The sources are caused by the gravity force.

3 Comparisons of theoretical solutions and DEM simulations

Fig. 2
figure 2

The DEM model for stress transmission in rhombic disk packings. The rhombic packing is formed by disks with same size. The rectangular boundary of this packing with length \(L_{x}\) and \(L_{y}\) is confined by four fixed walls,the left wall \(W_{L}\), the bottom wall \(W_{B}\), the right wall \(W_{R}\) and the top wall \(W_{Ta}\). A movable wall \(W_{Tb}\) with length \(\delta \) and center position \((x_{b},0)\) is prepared to realize a point load on the top of the packing

Full size image

The stress response caused by a point load [3, 4] and the Janssen effect [5, 6, 20] are often used to demonstrate the unusual mechanical behaviors of granular solids. DEM simulations can be used to investigate these two phenomenons in rhombic disk packings. We elaborate the protocol of our DEM simulations as follows. (1)As shown in Fig.2, a rhombic lattice packing of disks with radii \(r_{0}\) is generated and it is shaped in a rectangular with length \(L_{x}\) and height \(L_{y}\).The position of the center of a disk is calculated according to

$$\begin{aligned} {\left\{ \begin{array}{ll} x_{i,j}=r_{0}+2r_{0}(i-1)\sin \dfrac{\theta }{2}\\ y_{i,j}=r_{0}+2r_{0}(j-1)\cos \dfrac{\theta }{2} \end{array}\right. } \end{aligned}$$
(22)

where i and j are the index of this disk. Note that i and j should be both odd or even numbers. (2) The left wall \(W_{L}\), the bottom wall \(W_{B}\) and the right wall \(W_{R}\) are then generated as a fixed container. On the top of the rhombic packing, a fixed wall \(W_{Ta}\) with length \(L_{x}\) is generated and it is used to restrict the top layer of disks so that they couldn’t move out of the packing. To realize a point load on the top of the packing, we generate a movable wall \(W_{Tb}\) with length \(\delta \) and the initial position of its center is at \((x_{b},0)\). (3) The linear contact force is used to model the normal and shear interaction of each contact. The shear force represents the friction of each contact and it satisfies the Coulomb slip condition. In this paper, \(\mu \) denotes the friction coefficient between two disks while \(\mu ^{*}\) denotes the one between a disk and the container. Damping forces are involved to dissipate the kinetic energy of the granular system to obtain a steady-state packing in the DEM simulation. The gravitational acceleration is set to be \(\mathbf {g={\mathrm {g}}e_{y}}\) where \(\mathbf {\mathrm {g}}\) is a variable with a range \(0\le g\le 9.8m/s^{2}\). (4)If we perform the investigation for Janssen effect, \(W_{Tb}\) is fixed and \(\mathbf {\mathrm {g}}\) is increased from 0 to \(9.8m/s^{2}\) with step by step manner. After each change of the value of g, a DEM dynamic process with energy dissipation is applied to obtain a static lattice packing with the boundary friction. (5)When we perform the investigation for the stress response due to a point load, \(\mathbf {\mathrm {g}}\) is set to be zero and \(W_{Tb}\) is moved downward with step by step manner while other walls are fixed. After each step of the movement of \(W_{Tb}\), a DEM dynamic process with energy dissipation is applied to obtain a static lattice packing with a top point load. (6)We note that the displacement of \(W_{Tb}\) and the increment of \(\mathbf {\mathrm {g}}\) are very small for each step of the DEM simulation. The lattice structure of the disk packing will not be destroyed during the whole DEM simulation.

Based on Eqs.(18), (19) and (20), it is not hard to derive the exact solutions of stress transmission in rhombic disk packings with proper boundary conditions. In the following, we use these formulae to investigate Janssen effect and the stress response due to a point load in rhombic disk packings. In theory, both of these problems can be related to the solutions of Eq.(1) with the following boundary conditions

$$\begin{aligned} \left\{ \begin{array}{c} Q_{yy}(x,0)=-p(x),Q_{xy}(x,0)=-q(x)\\ Q_{yx}(\frac{L}{2},y)=\mu ^{*}\eta ^{*}\lambda Q_{yy}(\frac{L}{2},y)\\ Q_{yx}(-\frac{L}{2},y)=-\mu ^{*}\eta ^{*}\lambda Q_{yy}(-\frac{L}{2},y) \end{array}\right. \end{aligned}$$
(23)

where \(\eta ^{*}\) is the friction mobilization for contacts between a disk and the container. Note that the last two boundary conditions in Eq.(23) are given by the consideration of Coulomb’s law of static friction. Here, we make a hypothesis that all the friction mobilizations have the same value, but this may not be the case in experiments or DEM simulations. We could carry out the following three steps to obtain the distribution for stress components \(\left\{ Q_{xx},Q_{yy},Q_{xy}\right\} \). (1) Using Eq.(20) and (23) to derive the boundary conditions for \(\left\{ f_{1},f_{2}\right\} \). (2) Using Eq.(18) and (19) to calculate the value of \(\left\{ f_{1},f_{2}\right\} \) at a specified position \(\left\{ x,y\right\} \). (3) Calculating the value of \(\left\{ Q_{xx},Q_{yy},Q_{xy}\right\} \) at \(\left\{ x,y\right\} \) by using Eq.(20). Performing the first step, the top boundary condition of \(\left\{ f_{1},f_{2}\right\} \) can be derived as

$$\begin{aligned}&\left( \begin{array}{c} f_{1}(x,0)\\ f_{2}(x,0) \end{array}\right) =C^{-1}\left( \begin{array}{c} Q_{yy}(x,0)\\ Q_{xy}(x,0) \end{array}\right) \nonumber \\&\quad =\frac{-1}{2\sqrt{\lambda }}\left( \begin{array}{c} \sqrt{\lambda }p(x)+q(x)\\ \sqrt{\lambda }p(x)-q(x) \end{array}\right) \end{aligned}$$
(24)

Similarly, the left and the right boundary conditions of \(\left\{ f_{1},f_{2}\right\} \) can be obtained and be written as

$$\begin{aligned} \left\{ \begin{array}{c} f_{1}(-\frac{L}{2},y)=\gamma f_{2}(-\frac{L}{2},y)\\ f_{2}(\frac{L}{2},y)=\gamma f_{1}(\frac{L}{2},y) \end{array}\right. \end{aligned}$$
(25)

where \(\gamma \) is defined as

$$\begin{aligned} \gamma =\frac{1-\mu ^{*}\eta ^{*}\sqrt{\lambda }}{1+\mu ^{*}\eta ^{*}\sqrt{\lambda }} \end{aligned}$$
(26)

We discuss the physical meaning of Eq.(25) as following. If we regard y as time “ t ”, the unidirectional wave \(f_{2}\) travel along “ \(-x\) ” direction and a reflected wave \(f_{1}\) is produced when \(f_{2}\) meets the left boundary. The first equation in Eq.(25) shows that only part of the strength of \(f_{2}\) is reflected due to the friction between the left wall and the disk packing. We know that \(\sqrt{\lambda }\) acts like wave velocity. Apparently, \(\gamma \) has the meaning like reflectivity and it decreases as the effective friction coefficient \(\mu ^{*}\eta ^{*}\) increases. The second equation in Eq.(25) shows that the unidirectional wave \(f_{1}\) travel along “ \(+x\) ” direction and a reflected wave \(f_{2}\) is produced when \(f_{1}\) meets the right boundary. Based on the above boundary conditions of \(\left\{ f_{1},f_{2}\right\} \), the value of \(\left\{ Q_{xx},Q_{yy},Q_{xy}\right\} \) at any specified position \(\left\{ x,y\right\} \) could be calculated out by performing the second and the third steps. However, the analytical expression of \(\left\{ f_{1}(x,y),f_{2}(x,y)\right\} \) is not the composition of elementary functions if the wave \(\left\{ f_{1},f_{2}\right\} \) is reflected many times. Therefore, we calculate the numeric value of \(\left\{ Q_{xx},Q_{yy},Q_{xy}\right\} \) according to the above three steps method rather than write down the analytical expression in the following discussions. If g is not zero,we need to calculate \(\left\{ g_{1},g_{2}\right\} \) in Eqs.(18) and (19) by using

$$\begin{aligned} \left( \begin{array}{c} g_{1}\\ g_{2} \end{array}\right) =B'\left( \begin{array}{c} g_{x}\\ g_{y} \end{array}\right) =\left( \begin{array}{c} \frac{g}{2\sqrt{\lambda }}\\ \frac{-g}{2\sqrt{\lambda }} \end{array}\right) \end{aligned}$$
(27)

When grains are added to a container, the weight applied to the bottom is partially supported by the lateral walls due to the frictional interactions between the grains and lateral walls, therefore it is smaller than the total weight of the granular column. This is called Janssen effect [5, 6, 20]. In the following, we investigate the Janssen effect through the DEM simulations of rhombic disk packings and compare the DEM results with the theoretical solutions. The stress field \({\mathbf {Q}}(x,y)\) is calculated and shown in Fig.3 according to DEM data and theoretical solutions respectively. The normal pressure P on the bottom boundary is an important quantity for the investigation of Janssen effect and it varies with the depth of the bottom boundary y which is defined as the distance between the top wall and the bottom wall. The function P(y) can be calculated by using the following formula

$$\begin{aligned} P(y)=\frac{1}{L_{x}}\int _{-L_{x}/2}^{L_{x}/2}Q_{yy}(x,y)dx \end{aligned}$$
(28)

The DEM and theoretical result for P(y) are shown in Fig.4. Theoretical results shown in Figs.3 and 4 indicate that \({\mathbf {Q}}(x,y)\) and P(y) vary smoothly with spatial position and the saturation of P(y) can be given by our theory. However, the theoretical solutions are not completely consistent with the DEM results. The DEM simulations show that the stress field \({\mathbf {Q}}(x,y)\) doesn’t varies very smoothly and P(y) is not a monotone increasing function. The causes for the differences between the DEM simulations and the theoretical results will be discussed below.

Fig. 3
figure 3

The stress field in the rhombic disk packing for Janssen effect. \({\mathbf {g}}\) is the gravitational acceleration. a The DEM result with \(\mu =\mu ^{*}=0.5\). b The theoretical result with \(\mu \eta =0,\mu ^{*}\eta ^{*}=0.5\)

Full size image
Fig. 4
figure 4

The bottom pressure P varies with the depth y. Each curve is obtained according to the stress field \({\mathbf {Q}}(x,y)\) shown in Fig.3 respectively

Full size image

Next, we perform DEM simulations to investigate the stress response of the rhombic packing due to a point load and compare the DEM results with the theoretical results. In Fig.5 , we show the DEM results for force chains in rhombic disk packings each under a point load with different center position. The stress fields are shown in Fig.6 according to the DEM results and the theoretical solutions respectively. These figures indicate that the theoretical result and the DEM simulation are similar but not completely consistent. As shown in the DEM results (see Figs.5 and 6a), the force chain and the stress propagate mainly along two directions in the rhombic lattice packing and they will be partially reflected when they meet the boundaries. These phenomenons could be reproduced by our theoretical results (see Fig.6b–c). However, there appear many weak force chains around the strong force chains in Fig.5. This means that the stress will be scattered as it propagate through the rhombic packing. The DEM results also show the stress response bouncing off the bottom boundary and then propagating back up the packing. These phenomenons haven’t been reproduced in the current theoretical results.

Fig. 5
figure 5

DEM results for force chains in the rhombic packings each under a point load( \(\mu =\mu ^{*}=0.5\)). The magnitudes of the point loads \(F_{top}\) are the same for these three cases and we set \(F_{top}\approx 1N\)

Full size image
Fig. 6
figure 6

Stress fields in rhombic disk packings each under a point load. \(F_{top}\approx 12N\) is the magnitude of the point load . a The DEM result with \(\mu =\mu ^{*}=0.5\). b The theoretical result with \(\mu \eta =\mu ^{*}\eta ^{*}=0\).(c) The theoretical result with \(\mu \eta =\mu ^{*}\eta ^{*}=0.5\)

Full size image

The differences between the DEM simulations and the current theoretical results could be understood by considering the failures of the basic assumptions in our theoretical model. In our theoretical calculations, the boundary conditions (see Eq.(23)) could make Eq.(6) be the well-posed problem, but the bottom boundary is not involved in the boundary conditions. Hence, it doesn’t show the stress chain bouncing off the bottom boundary in the theoretical results (see Fig.6b–c). To reproduce this phenomenon, Eq.(23) should be modified to involve the bottom boundary condition. However, Eq.(6) are hyperbolic equations, the stress boundary conditions are not independent for the four boundaries of the rhombic packing. Therefore, it is not straightforward to set the boundary conditions with the bottom boundary. In this paper, we also adopt the theoretical hypothesis that the effective friction coefficients \(\mu \eta \) and \(\mu ^{*}\eta ^{*}\) are constants with the variation of spatial position. Therefore, in our theory, the Janssen coefficient \(\lambda \) is a constant and Eq.(6) are constant coefficient differential equations. This makes that the stress field \({\mathbf {Q}}(x,y)\) varies smoothly with spatial position (see Fig.3b) and the stress propagates in the rhombic packing without scattering (see Fig.6b–c).

In practical situations, the friction of each contact is usually mobilized by random factors and this causes that the friction mobilization \(\eta \) and \(\eta ^{*}\) vary with spatial position. In Fig.7, we show the probability distribution of the friction mobilization \(\hbox {P}_{D}(\eta )\) for the Janssen effect and the stress response under a point load respectively. The local maximum value points A1, A2, B1, B2, B3 provide the key information of the two curves in Fig.7. A2 is the maximum value point of \(\hbox {P}_{D}(\eta )\) of the Janssen effect. Let \(\eta (A2)\) be the value of the friction mobilization \(\eta \) at point A2. \(\eta (A2)\approx 1\) means that the frictional forces of most contacts are almost full mobilized in the DEM simulation of the Janssen effect. B2 is the maximum value point of \(\hbox {P}_{D}(\eta )\) of the point load response. Unlike the Janssen effect, \(\eta (B2)\approx 0\) and this means that most frictional forces are almost never mobilized in the DEM simulation of the point load response. \(\eta (B2)\approx 0\) also explains why the theoretical result with \(\mu \eta =\mu ^{*}\eta ^{*}=0\) is closer to the DEM result than the theoretical result with \(\mu \eta =\mu ^{*}\eta ^{*}=0.5\)(see Fig.6a–c). The existence of other local maximum value points A1, B1, B3 confirms that the friction mobilization \(\eta \) (or the friction angle \(\beta \)) is not a constant and it varies with spatial position. In our DEM simulations of the rhombic packings without friction, the weak force chains can also be found around the strong force chains. This means that the rhombic disk packing will be distorted slightly when the external stress acts on the packing. Hence, the angle \(\theta \) varies slightly with spatial position. According to Eq.(5), we conclude that the Janssen coefficient \(\lambda \) usually varies randomly with spatial position in the rhombic packing which is generated by the DEM simulation. Due to this effect, Eq.(6) become the differential equations with random coefficient. Therefore, as shown in the DEM results, the stress will be scattered as it propagates through the rhombic packing.

Fig. 7
figure 7

The probability distribution of the friction mobilization \(\hbox {P}_{D}(\eta )\) and (inset) the assumed direction of contact forces on each internal disk in theory. The dashed line with circles is obtained from the DEM results for Janssen effect (see Fig.3a). The dotted line with stars is obtained from the DEM results for point load response (see Fig.6a). A1, A2, B1, B2, B3 are the local maximum value points of the two curves respectively and their physical implications are discussed in the main text. \({\mathbf {k}}_{1},{\mathbf {k}}_{2}\) are the unit lattice vectors of a rhombic packing. \({\mathbf {f}}_{\tau }\) and \({\mathbf {f}}_{N}\) represent the tangential and normal component of contact forces respectively. The negative value of \(\eta \) means that the real direction of the tangential component of a contact force is opposite to the assumed direction of \({\mathbf {f}}_{\tau }\) shown in the inset

Full size image

The mechanical properties of a packing mainly depend on those factors such as the contact network of the packing, the distribution of friction mobilization, and so on. In this work, we focus on the stress transmission in rhombic disk packings which have simple contact structures. However, the contact network of a real granular packing usually appears as other kinds of crystal lattice or amorphous structures. If the grains are on other kinds of crystal lattice and the packing is in isostatic state, the constitutive equation of such granular packing is usually an algebraic equation like \(\mathbf {Q:G=0}\) [21]. In this equation, the friction-geometry tensor \({\mathbf {G}}\) should represent the contact structure of this lattice packing and it doesn’t vary with spatial position. Therefore, our theoretical approach could be applied to investigate the stress transmission in other isostatic lattice packings. Even for the packing with periodic contact network, if it is in hyperstatic state, the uncertainty of contact force will cause complex behaviors of stress transmission. In amorphous granular packings, due to the randomness of their contact structures, the stress transmit in a very different pattern with respect to the lattice packings [4]. For these complex situations, \({\mathbf {G}}\) may varies with spatial position, or the constitutive equation may become other forms of equation(for example, the differential equation [16, 21]), or the proper constitutive equation may be absent altogether. Hence, our theoretical approach couldn’t be directly applied to these complex situations.

4 Discussion and conclusion

The rhombic disk packing is the simplest structure for the investigation of stress transmission in granular matter. Based on the hypothesis that the friction mobilization \(\eta \) has the same value for different contacts, we derived the governing equations of stress transmission in rhombic lattice packings of disks with friction [19]. In the current article, the exact solutions of the governing equations for stress tensor field \(\mathbf {Q{\text {(x,y)}}}\) are derived and they are used to calculate the stress distribution in the rhombic disk packing with lateral and top boundaries. Our theoretical results show that the stress transmit in lattice packings with the behavior of traveling wave with sources. The square root of Janssen coefficient \(\sqrt{\lambda }\) plays a role of wave speed for the stress transmission.The traveling stress wave will be reflected when it meets the boundaries and it will be partially absorbed due to the friction between the boundaries and the rhombic packing. DEM simulations and theoretical solutions are also used to investigate the Janssen effect and the stress response under a point load for the rhombic disk packings. The comparisons of DEM results and theoretical results indicate that our theoretical solutions could describe the main features of stress distribution in the rhombic disk packings, but they also have some obvious differences. The bottom boundary is not involved in our boundary condition settings. A rhombic disk packing with friction is usually in hyperstatic state in which the friction mobilization \(\eta \) and \(\eta ^{*}\) vary with spatial position. A rhombic disk packing may also be distorted slightly when the external stress acts on the packing. These factors cause the deviation between the DEM simulations and our theoretical results. Therefore, these complex factors should be considered in the future theoretical works to obtain more reasonable results.