⇦Chinese Physics Letters, 2019, Vol. 36, No. 2, Article code 024501 Dynamical Study of Granular Flow through a Two-Dimensional Hopper * Lina Yang (杨丽娜)1,2, Yu-Qi Chen (陈裕启)1,2** Affiliations 1Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190 2School of Physical Sciences, University of Chinese Academy of Science, Beijing 100049 Received 30 October 2018, online 22 January 2019 *Supported by the Key Research Program of Frontier Science of the Chinese Academy of Sciences under Grant No QYZDY-SSW-SYS006.
**Corresponding author. Email: ychen@itp.ac.cn Citation Text: Yang L N and Chen Y Q 2019 Chin. Phys. Lett. 36 024501    Abstract The mass flow rate of a granular flow through an aperture under gravity is a long-standing challenge issue in physical science. We show that for steady flow field close to laminar flow, the dynamical equations together with the continue equation and Mohr-circle description of the stress are closed, and hence solvable. In a case of streamline guided by the two-dimensional hopper, we obtain a consistent condition and use it to determine the stress and the velocity distribution. Our result indicates that 3/2 power scaling behavior is recovered with a coefficient $C(\mu,\alpha)$ being a function of frictional coefficient and the hopper angle. It is found that the predicted coefficient $C(\mu,\alpha)$ is compatible with previous studies. DOI:10.1088/0256-307X/36/2/024501 PACS:45.70.Mg, 47.57.Gc, 45.70.-n © 2019 Chinese Physics Society Article Text The mass flow rate of a granular flow through an aperture under gravity is a long standing puzzle in physical science. So far, numerous works have been carried out and considerable progress has been made. Beverloo[1] proposed an empirical formula for the mass flow rate based on a throughout experimental observation. Digital particle image velocimetry measurements (DPIV)[2] were developed to allow the precise measurement of the flow field by virtue of computer technology.[3] In addition, Choi et al.[4] proposed a Gaussian solution for velocity distribution in a silo to describe the velocity distribution. The other method is based on the numerical simulation of granular flow, which is largely applied in industry. By analyzing the data, one can figure out the approximate solution of the velocity distribution,[5,6] pressure distribution, etc. An alternative way to study the granular flow is to describe it as a continuous flow, which satisfies the dynamical flow equations.[7] A typical work along this line was carried out by Savage et al.[8,9] Starting from the partial differential equations describing the flow, by averaging contributions from one of two degrees of freedom they obtained an ordinary differential equation and then gained solution. However, this kind of approximation was oblique from the theoretical point of view. In this Letter, we study a granular flow in the framework of continuous flow based on the following considerations. Notice that the researches in both experiments and numerical simulations support a steady flow field with streamlines which is close to a laminar flow. This means that there is a steady velocity distribution in the flow field. The granular flow along the streamline is under both the gravitational force and the inner stress force to be determined. A velocity difference between two different streamlines implies that the granular flow is under the most sliding frictional force, which makes the inner pressure force predictable by the Mohr–Coulomb yield condition. Once the stress force is determined, one can predict the granular flow using dynamical methods. We then show that those dynamical equations are self-contained and closed, and hence are solvable under certain boundary conditions. To simplify the problem, we consider the granular flow through a hopper at the bottom with a silo which plays a guider for the streamline. Furthermore, as in Ref. [9] we take a coarse hopper so that the velocity on the hopper surface vanishes. This boundary condition is supported by the experimental observation,[10] where the velocity drops very quickly at the side edge of the outlet. We also assume that the size of the granular is much smaller than the size of the silo. Giving the streamline, we obtain a consistent condition which is used to obtain an approximate solution for the dynamical equations of the granular flow. We now write down all equations satisfying the flow field. For a steady flow with the same density $\rho$, the first equation is the equation of continuity; i.e., $$ \frac{\partial V_x}{\partial x}+\frac{\partial V_y}{\partial y}=0,~~ \tag {1} $$ where $V_x$ and $V_y$ are the components of the velocity in $x$ and $y$ directions, respectively. For steady streamlines with a streamline slope angle $\theta$, $$ V_x=V_y\cot\theta.~~ \tag {2} $$ Thus Eq. (1) can be rewritten as $$ \frac{\partial V_y}{\partial y}+\frac{\partial V_y}{\partial x}\cot\theta+V_y\frac{\partial\cot\theta}{\partial x}=0.~~ \tag {3} $$ The dynamical equation can be obtained by applying Newton's second law of motion to the flow. Given $\rho{\boldsymbol F}$ being the body force per unite volume and $P$ being the surface force per unit area which has a normal direction ${\boldsymbol n}$, the dynamical equation is $$ \rho({\boldsymbol V}\cdot\nabla){\boldsymbol V}=\rho{\boldsymbol F}+\nabla\cdot P.~~ \tag {4} $$ In our case, the gravity is the only source of the body force. The surface force is a second-rank stress tensor. In two dimensions, it can be expressed by a $2\times 2$ symmetrical matrix $$ P=\left[ \begin{matrix} \sigma_{xx} & \tau \\ \tau & \sigma_{yy} \end{matrix}\right],~~ \tag {5} $$ where $\sigma$ and $\tau$ are normal and shear stresses, respectively. The dynamical Eq. (4) can then be written in a component form $$\begin{alignat}{1} \rho\Big(V_x \frac{\partial V_x}{\partial x}+V_y \frac{\partial V_x}{\partial y}\Big)=\,&\frac{\partial\sigma_{xx}}{\partial x}+\frac{\partial\tau}{\partial y},~~ \tag {6} \end{alignat} $$ $$\begin{alignat}{1} \rho\Big(V_x \frac{\partial V_y}{\partial x}+V_y \frac{\partial V_y}{\partial y}\Big)=\,&-\rho g+\frac{\partial\tau}{\partial x}+\frac{\partial\sigma_{yy}}{\partial y}.~~ \tag {7} \end{alignat} $$ The basic equations that the steady granular flow satisfies are dynamical Eqs. (6) and (7), together with continuity Eq. (1). We can see that there are five physical variables $V_x$, $V_y$, $\sigma_{xx}$, $\sigma_{yy}$ and $\tau$ to be determined. However, there are only three differential equations. Hence, these equations alone are not closed. In the case of steady flow that we consider, there are two additional conditions to constrain the stress tensor, which means that these equations are closed. We notice that for the steady flow, there is a velocity difference along the interface which is perpendicular to the streamline. This means that there exists a relative slide between different streamlines. From the Mohr–Coulomb yield condition, we know that $|\tau|\le\mu\sigma+c$, where $\mu$ stands for the tangent value of the internal friction angle $\beta$, and $c$ is the cohesion. This states that a slide may occur when the shear stress $\tau$ takes a maximal values under the normal stress, i.e., $|\tau|=\mu\sigma$. Here, for simplicity only, we take the granular material as cohesionless; that is, $c=0$. The stress tensor in the coordinate of the streamline can be written as a matrix $$ P=\left[\begin{matrix} \sigma & \lambda\mu\sigma \\ \lambda\mu\sigma & (2\mu^2+1)\sigma \\ \end{matrix} \right],~~ \tag {8} $$ where $\lambda=\pm1$, with the sign determined by the velocity gradient along the direction perpendicular to the streamline. For $x>0$, the velocity gradient is negative. Hence, we have $\lambda=-1$. Similarly, for $x < 0$, $\lambda=+1$. From Eq. (8) we can see that there is only one independent component $\sigma$ left in the stress tensor. To obtain the stress tensor in the $x$–$y$ direction, we make a coordinate rotation. We then have the stresses in the conversional $x$ and $y$ directions with $x>0$, $$\begin{alignat}{1} \sigma_{xx}=\,&(1+\mu^2)\sigma+\mu^2\sigma\cos(2\theta)-\mu\sigma\sin(2\theta),~~ \tag {9} \end{alignat} $$ $$\begin{alignat}{1} \sigma_{yy}=\,&(1+\mu^2)\sigma-\mu^2\sigma\cos(2\theta)+\mu\sigma\sin(2\theta),~~ \tag {10} \end{alignat} $$ $$\begin{alignat}{1} \tau=\,&\mu^2\sigma\sin(2\theta)+\mu\sigma\cos(2\theta).~~ \tag {11} \end{alignat} $$ The left hand side of Eqs. (6) and (7) can be simplified by Eqs. (2) and (3) as follows: $$\begin{alignat}{1} \rho\Big(V_x \frac{\partial V_x}{\partial x}+V_y \frac{\partial V_x}{\partial y}\Big)=\,&\rho V_y^2\frac{\partial\cot\theta}{\partial y},~~ \tag {12} \end{alignat} $$ $$\begin{alignat}{1} \rho\Big(V_x \frac{\partial V_y}{\partial x}+V_y \frac{\partial V_y}{\partial y}\Big)=\,&-\rho V_y^2\frac{\partial\cot\theta}{\partial x}.~~ \tag {13} \end{alignat} $$ With Eqs. (9)-(11), finally Eqs. (6) and (7) are simplified to $$\begin{alignat}{1} \rho V_y^2\frac{\partial\cot\theta}{\partial y}=\,&k_{11}\frac{\partial\sigma}{\partial x}+k_{12}\frac{\partial\sigma}{\partial y}+k_{13}\sigma,~~ \tag {14} \end{alignat} $$ $$\begin{alignat}{1} -\rho V_y^2\frac{\partial\cot\theta}{\partial x}+\rho g=\,&k_{21}\frac{\partial\sigma}{\partial x}+k_{22}\frac{\partial\sigma}{\partial y}+k_{23}\sigma ,~~ \tag {15} \end{alignat} $$ where $$\begin{align} k_{11}=\,&1+\mu^2+\mu^2\cos(2\theta)-\mu\sin(2\theta),~~ \tag {16} \end{align} $$ $$\begin{align} k_{12}=\,&k_{21}=\mu^2\sin(2\theta)+\mu\cos(2\theta),~~ \tag {17} \end{align} $$ $$\begin{align} k_{22}=\,&1+\mu^2-\mu^2\cos(2\theta)+\mu\sin(2\theta),~~ \tag {18} \end{align} $$ $$\begin{align} k_{13}=\,&\frac{\partial k_{11}}{\partial x}+\frac{\partial k_{12}}{\partial y},~~ \tag {19} \end{align} $$ $$\begin{align} k_{23}=\,&\frac{\partial k_{21}}{\partial x}+\frac{\partial k_{22}}{\partial y}.~~ \tag {20} \end{align} $$ The basic equations to solve our problem are Equations (14) and (15), together with Eq. (3). These three equations must be satisfied at every point in the field. We now see that there are only three unknown variables in the three basic equations: $V_y$, $\theta$ and $\sigma$. Thus the equations are closed and can be solved under proper boundary conditions. Savage et al. derived the equations to describe granular flow along radial direction through an aperture under gravity under the Mohr–Coulomb yield condition.[8,9] The equations obtained here generalize the equations of Savage et al. to be suitable to describe an arbitrary laminar flow under the Mohr–Coulomb yield condition in two dimensions. The basic Eqs. (14) and (15) together with Eq. (3) are nonlinear partial differential equations. We combine them to simplify those equations to make them easier to be solved. From Eqs. (14) and (15), we can rewrite the expressions of $\frac{\partial\sigma}{\partial x}$ and $\frac{\partial\sigma}{\partial y}$ as follows: $$\begin{align} \frac{\partial\sigma}{\partial x}=\,&c_{11}\sigma+c_{12},~~ \tag {21} \end{align} $$ $$\begin{align} \frac{\partial\sigma}{\partial y}=\,&c_{21}\sigma+c_{22},~~ \tag {22} \end{align} $$ where $$\begin{align} c_{11}=\,&k(k_{12}k_{23}-k_{13}k_{22}),~~ \tag {23} \end{align} $$ $$\begin{align} c_{12}=\,&k k_{22}\rho V_y^2\frac{\partial\cot\theta}{\partial y}+ k k_{12}\rho\Big(V_y^2\frac{\partial\cot\theta}{\partial x}-g\Big),~~ \tag {24} \end{align} $$ $$\begin{align} c_{21}=\,&k(k_{13}k_{21}-k_{11}k_{23}),~~ \tag {25} \end{align} $$ $$\begin{align} c_{22}=\,&-k k_{21}\rho V_y^2\frac{\partial\cot\theta}{\partial y}-k k_{11} \rho\Big(V_y^2\frac{\partial\cot\theta}{\partial x}-g\Big),~~ \tag {26} \end{align} $$ $$\begin{align} k=\,&k_{11}k_{22}-k_{12}k_{21}=(\mu^2+1)^{-1}.~~ \tag {27} \end{align} $$ Furthermore, the second partial derivative of $\sigma$, that is, $\frac{\partial^2 \sigma}{\partial x\partial y}$, can also be derived from Eqs. (21) and (22), $$\begin{align} \frac{\partial^2\sigma}{\partial x\partial y}=\,&c_{11}\frac{\partial\sigma}{\partial y}+\frac{\partial c_{11}}{\partial y}\sigma+\frac{\partial c_{12}}{\partial y},~~ \tag {28} \end{align} $$ $$\begin{align} \frac{\partial^2\sigma}{\partial x\partial y}=\,&c_{21}\frac{\partial\sigma}{\partial x}+\frac{\partial c_{21}}{\partial x}\sigma+\frac{\partial c_{22}}{\partial x}.~~ \tag {29} \end{align} $$ Substituting $\frac{\partial\sigma}{\partial x}$ and $\frac{\partial\sigma}{\partial y}$ with Eqs. (21) and (22) we have $$\begin{alignat}{1} \frac{\partial^2\sigma}{\partial x\partial y}=\,&\Big(c_{11}c_{21}+\frac{\partial c_{11}}{\partial y}\Big)\sigma+c_{11}c_{22}+\frac{\partial c_{12}}{\partial y},~~ \tag {30} \end{alignat} $$ $$\begin{alignat}{1} \frac{\partial^2\sigma}{\partial x\partial y}=\,&\Big(c_{11}c_{21}+\frac{\partial c_{21}}{\partial x}\Big)\sigma+c_{12}c_{21}+\frac{\partial c_{22}}{\partial x},~~ \tag {31} \end{alignat} $$ from which we then obtain the expression of $\sigma$ as follows: $$\begin{align} \sigma=-\frac{c_{11}c_{22}-c_{12}c_{21}+\frac{\partial c_{12}}{\partial y}-\frac{\partial c_{22}}{\partial x}}{\frac{\partial c_{11}}{\partial y}-\frac{\partial c_{21}}{\partial x}}.~~ \tag {32} \end{align} $$ With Eq. (32), we can obtain the relation between $V_y$ and $\cot\theta$ from Eqs. (14), (15), (30) or (31). Generally, the velocity distributions are very complicated in real flow field. To make these equations be solved easily, we choose a process with a relatively simpler steady flow.
cpl-36-2-024501-fig1.png
Fig. 1. Illustration for the hopper and the granular flow.
To this end, we choose a hopper, as shown in Fig. 1. The hopper is divided into two parts: the upper region and the lower one. The streamline is guided by the shape of the hopper as shown by the arrow line in Fig. 1. This can be thought of as a leading order approximation of streamline. To make this approximation valid, the angle $\pi/2-\alpha$ should not be too large. Furthermore, we take a coarse hopper as in Ref. [9], so that the velocity on the hopper surface vanishes. Under this boundary condition, the solution of the equations for granular flow can be considerably simplified. We first solve the equation for the upper region of the hopper. Since it is a vertical flow, we have $\theta=\pi/2$ and $V_x=0$. It follows from Eq. (3) that $\frac{\partial V_y}{\partial y}=0$. Thus, $V_y$ is a function of $x$ only; that is, $V_y=v(x)$. We can solve the distribution of stress $$ \sigma=\frac{\mu\rho g |x|}{\mu^2+1}+\frac{\rho g y}{\mu^2+1}+c_0,~~ \tag {33} $$ where $c_0$ is a constant which can be determined by the boundary condition of the stress. This shows that the stress at the center is smaller than that at other places and the stress at higher place is larger. However, there is no constraint on the value of $v(x)$. We now turn to analyze the lower region of the hopper. In the lower region of the hopper, we assume that the streamline is a set of straight lines converging to the center as a radial flow. As pointed out previously, at the leading order approximate streamline, the function of $\cot\theta$ is $$ \cot\theta=\frac{x\cot\alpha}{y\cot\alpha+r},~~ \tag {34} $$ where $r$ is half the width of the aperture, and $\alpha$ is the slope angel of the wall. Equation (3) becomes a first-order partial differential equation of $V_y$, $$ \frac{\partial V_y}{\partial y}+\frac{\partial V_y}{\partial x}\frac{x\cot\alpha}{y\cot\alpha+r}+V_y\frac{\cot\alpha}{y\cot\alpha+r}=0.~~ \tag {35} $$ This means that we can obtain the velocity distribution among the lower region by solving Eq. (35) under the boundary condition $V_y=0$ at the surface. On the interface between the upper and the lower regions, the continuity equation implies that the direction of the velocity changes while its absolute value remains unchanged. The velocities in two regions are related at the boundary as follows: $$ V_y^2=v^2(x)=V_y(x,y)^2(1+\cot^2\theta)|_{y=f(x)},~~ \tag {36} $$ where $y=f(x)$ stands for the interface and may have many possible shapes, such as the convex and concave arcs shown in Fig. 1. In the real case, there must be a transition region so that the streamline is continuous and smooth through the whole field. The shape of the interface can be parameterized by a simple function $f(x)=q(x^2-a^2)+\frac{a-r}{\cot\alpha}$ with $a$ being half the width of the upper region of hopper and $q$ being a free parameter to reflect the variation of the shape of the interface. We now start to solve the three basic equations under the coarse boundary condition. Since the continuity equation becomes a first-order linear partial different equation, it can be solved using its characteristic equation. Solving Eq. (35) under the boundary condition, we obtain the velocity distribution in the lower region as follows: $$\begin{align} V_y=\frac{(f(x_{\rm b})+\frac{r}{\cot\alpha})G(x_{\rm b})}{r+y\cot\alpha},~~ \tag {37} \end{align} $$ where $G(x)=-\sqrt{\frac{v^2(x)}{1+\frac{x}{f(x)+\frac{r}{\cot\alpha}}}}$, and $x_{\rm b}$ is the value of $x$ at the interface. Comparing Eqs. (30) and (31), we find that $\frac{\partial c_{11}}{\partial y}=\frac{\partial c_{21}}{\partial x}$, which arises from the radial streamline. This implies that $c_{11}c_{22}+\frac{\partial c_{12}}{\partial y}$ in Eq. (30) and $c_{12}c_{21}+\frac{\partial c_{22}}{\partial x}$ in Eq. (31) must be equal. This is the self-consistent condition of the stress. It conducts $$ \frac{2\cot\alpha\rho}{\mu^2+1}\Big(S_1+S_2V_y^2+S_3\frac{\partial V_y^2}{\partial x}\Big)=0,~~ \tag {38} $$ where $$\begin{align} S_1=\,&\frac{2\mu g(r+y\cot\alpha-x\mu\cot\alpha)}{((r+y\cot\alpha)^2+(x\cot\alpha)^2},~~ \tag {39} \end{align} $$ $$\begin{align} S_2=\,&\frac{2x\cot^2\alpha}{(r+y\cot\alpha)^3},~~ \tag {40} \end{align} $$ $$\begin{align} S_3=\,&\frac{(r+y\cot\alpha)^2+(x\cot\alpha)^2}{2(r+y\cot\alpha)^3}.~~ \tag {41} \end{align} $$ By substituting Eq. (37) into Eq. (38), we have a first-order differential equation of $v(x)$. We find that the solution of $v(x)$ depends on values of both $x$ and $y$. We may set $y$ at a fixed value $y_0$. Under the boundary condition $v(\pm a)=0$, and we eventually obtain the distribution of velocity in the lower region $$\begin{align} V_y=\,&-[2\mu g(r+y_0\cot\alpha)^3(r+y\cot\alpha-x)(2(r\\ &+y\cot\alpha)-\mu\cot\alpha(r+y\cot\alpha+x))]^{1/2}/[(r\\ &+y\cot\alpha)^2+(x\cot\alpha)^2].~~ \tag {42} \end{align} $$ We find that the dependences of parameters $q$ and $a$ in the parameterized interface vanish. This means that our results are insensitive to the shape of the interface. The velocity distribution is shown in Fig. 2(a). Here we take $\mu=\tan{35^{\circ}}$, $y_0=0.01$ m, $r=0.01$ m, and $\cot\alpha=1$ to carry out the actual calculation.
cpl-36-2-024501-fig2.png
Fig. 2. (a) The velocity distribution $V_y$. (b) The acceleration distribution $a_y$. (c) The velocity distribution at outlet. (d) The velocity distribution $|V|$ at the interface.
With radial stream-tube flow assumption, Brown[11] predicted that a free-fall arch is formed by considering the energy decrease during granular flow along the tube and reaching minimal value at the arch. With a further assumption that the arch is a round one, Brown predicted the flow rate. His results seem to fit experimental measurements well.[11] By examining the acceleration to reach gravitational acceleration $g$ from the velocity distribution, we also find that there exists such an arch but not a round one. The velocity of granular flow above the free-fall arch in the lower region is described by Eq. (42). Below the arch, the granular flow undergoes the free-fall motion. The lower region is further divided into two parts by the arch. We use several methods to obtain the velocity at the outlet of the hopper. For the part below the arch, we take velocity at the arch as the initial condition and use the free-fall motion to predict the velocity at $y=0$, while for the part above the arch, we just use Eqs. (2) and (42) to calculate the velocity. The velocity distribution at the outlet where $y=0$ is shown in Fig. 2(c). Notice that the vertical velocity from the work of Janda et al.[10] is slightly flatter than ours. The difference may be caused by different shapes of the aperture and property of granular material. These may lead to the different boundary conditions and stress conditions that can affect the results. In addition, the assumption of the streamline indicates that the velocity of the flow through a coarse hopper is converging. This will cause a small increase near the center. As pointed out previously, this result shows that the velocity distribution in the lower region is independent of the shape of the interface. We show the velocity distribution of $V$ at a parameterized interface with $q=0$ and $a=0.06$ m in Fig. 2(d). From the data in Fig. 2(b) we can obtain the free-fall arch where $a_y=-9.8$ m/s$^2$. Below the free-fall arch, the granular flow undergoes the free-fall motion as we refer above. Hence the acceleration distribution at the outlet of the hopper is constant, that is, $a_y=-9.8$ m/s$^2$, as shown in Fig. 3(c).
cpl-36-2-024501-fig3.png
Fig. 3. (a) The relation between mass flow rate $W$ and the half width $r$ in double logarithmic coordinates. (b) The numerical result of free-fall arch determined by gravity acceleration. (c) The coefficient of mass flow rate $C(\mu,\alpha)$ with different cone half angles $\delta$.
The mass flow rate given by Beverloo et al.[1] in three dimensions is reduced in a form of $W=C\rho_{\rm B}\sqrt{g}r^{3/2}$ in two dimensions,[12] where $\rho_{\rm B}$ is the bulk density of the particle. In Fig. 3(a), the relation between mass flow rate and the half width is shown in double logarithmic coordinates. From the figure we see that it is a straight line with slope 3/2. It is in agreement with the dimensional analysis. Here, we calculate the mass flow rate at $y_0=1.5r$. The mass flow rate $W$ mainly depends on the density $\rho_{\rm B}$, the gravitational acceleration $g$ and half width $r$, and $W$ also depends on the hopper's parameter $\cot\alpha$ and the frictional coefficient $\mu$. A dimensionless quantity $C(\mu,\alpha)$ is defined by $$ C(\mu,\alpha)=\frac{W}{\rho_{\rm B}\sqrt{g}r^{3/2}},~~ \tag {43} $$ as functions of $\mu$ and $\alpha$. Here, we focus on a situation of $\mu\cot\alpha=1$. The cone half angel $\delta$ of the hopper is related to the slope angel of wall $\alpha$ by $\delta=\pi/2-\alpha$. The relationship between the coefficient of mass flow rate $C(\mu,\alpha)$ and different cone half angles $\delta$ is shown in Fig. 3(c). It is interesting to compare our result with Brown's conclusion.[11] First, Brown predicted an free-fall arch, while we also find that there exists such an arch. Thus, our result supports Brown's free-fall arch result. However, Brown assumed that the arch is round, while the arch obtained in our work is not exactly round. In principle, in our approach, this arch can be obtained by solving Eqs. (14) and (15), together with Eq. (3) rigorously by setting $\sigma=0$. However, this is hard to carry out. Alternatively, in this study, we obtain the arch by examining the acceleration from the obtained velocity distribution under radical flow assumption. We then find that the obtained arch is not exactly round. Actually, in our approach the shape of the arch is dynamical dependent while the round arch in Brown's work is a conjecture. We also find that our predicted flow rate is in agreement with Brown's when $\delta$ is small. However, when the cone half angle of hopper is increasing, the difference between these two solutions becomes much larger. This is not surprising. Actually, as we discussed previously, the radical flow assumption is not so good for a large angle hopper where we may not expect that the obtained result is very reliable. Further work needs to be carried out to gain a more reliable result for this case with our approach. The difference between our work and Savage[8] is mostly caused by taking different means of approximation. The dynamical equation of Savage for the two-dimensional flow is an ordinary differential equation. In summary, the steady flow in the two-dimensional hopper has been solved by setting the straight streamline in the leading-order approximation with streamline fixed by the shape of the hopper. The velocity distribution roughly is in agreement with the experimental measurements at a smaller cone half angle $\delta$. Our results show that the velocity distribution in the lower region is insensitive to the interface at the head of hopper. This result is reliable for smaller $\delta$. However, in a general case, the streamline should be determined by solving the basic equations rigorously, which is hard due to the nonlinear partial differential equations. To achieve this goal, further work is under investigation. References Optical technique DPIV in measurements of granular material flows, Part 1 of 3—plane hoppersJamming during the discharge of granular matter from a siloVelocity profile of granular flows inside silos and hoppersSize Scaling of Velocity Field in Granular Flows through AperturesInfluence of granule velocity on gravity-driven granular flowThe mass flow of granular materials derived from coupled velocity-stress fieldsGravity flow of coarse cohesionless granular materials in conical hoppersFlow Rate of Particles through Apertures Obtained from Self-Similar Density and Velocity ProfilesMinimum Energy Theorem for Flow of Dry Granules through AperturesGravity discharge rate of fine particles from hoppers
[1]Beverloo W A, Leniger H A and Velde J V D 1961 Chem. Eng. Sci. 15(3) 260
[2] Sielamowicz I, Blonski S and Kowalewski T A 2005 Chem. Eng. Sci. 60 589
[3] Zuriguel I, Garcimartín A, Maza D, Pugnaloni L A and Pastor J M 2005 Phys. Rev. E 71 051303
[4] Choi J, Kudrolli A and Bazant M Z 2005 J. Phys.: Condens. Matter 17 S2533
[5] Hu G, Lin P, Zhang Y, Li L, Yang L and Chen X 2017 arXiv:1708.09647[cond-mat.soft]
[6] Huang D, Sun G and Lu K 2011 Phys. Lett. A 375 3375
[7]He R, Liu R, Lu K Q, Hou M Y, Zhang Q Y and Peng Z 2007 Acta Phys. Sin. 56 4708 (in Chinese)
[8] Savage S B 1965 Br. J. Appl. Phys. 16 1885
[9] Savage S B and Sayed M 1981 Z. Angew. Math. Phys. ZAMP 32 125
[10] Janda A, Zuriguel I and Maza D 2012 Phys. Rev. Lett. 108 248001
[11] Brown R L 1961 Nature 191 458
[12] Spink C D and Nedderman R M 1978 Powder Technol. 21 245