Simple Particle Model for Low-Density Granular Flow Interacting with Ambient Fluid

To understand the time evolutions of frontal speed and shape in a low-density granular flow, we propose a simple particle model. This model solves the equation of motion for each particle and simulates the time evolution of low-density granular flow. Spherical particles constituting a low-density granular flow slide on a slope at a steeper angle than the angle of repose. The particle motion is determined based on three forces: gravity as the driving force, repulsive force due to particle collision, and drag force due to the particle interaction through the ambient fluid. Two-dimensional numerical simulations of this model are conducted on the slope: the x–y plane parallel to the slope and the x–z plane perpendicular to the slope. In the x–y plane, particles aggregate at the moving front of the granular flow, and subsequently, flow instability occurs as a wavy pattern. This flow pattern is caused by the interparticle interaction arising from the drag force. Additionally, a vortex convection of particles is formed inside the aggregations. Simultaneously, particle aggregation is also found at the moving front of the granular flow in the x–z plane. The aggregation resembles a head–tail structure, where the frontal angle against the slope approaches 60 ∘ from a larger angle as time progresses. Comparing the numerical result by varying the particle size reveals that the qualitative dynamics of the granular flow are independent of particle size. Although the model is not realistic, our study presents a new particle-based approach that elucidates the dynamics of low-density granular flow.


Introduction
Snow or rock avalanches in nature are generally regarded as a class of massive slide phenomena involving gravity and density currents. These flows slide down a slope as a mixture of solid and fluid, and exhibit various patterns and complicated inner structures. For instance, powder snow avalanches are roughly divided into two structural regions according to field observations [1,2]: dense and dilute regions. The former is formed in the vicinity of a slope, whereas the latter develops above the dense region. In general, the dense region consists of coarse or high-density particles that are packed densely in the flow. The dilute region (the powder cloud) comprises fine particles or low-density particles that are suspended in the air. The driving force of both flows is gravity. However, the resisting force is different: inner and basal frictions govern the dense flow, whereas drag force governs the powder cloud. Therefore, these resisting forces lead to a major difference in the structure and dynamics.
A number of laboratory experiments on dense flow have been performed [3][4][5][6][7][8][9]. Dense granular flows consisting of monodisperse glass beads form the fingering pattern at the front and the streaky

Concept of the Proposed Model
We propose a simple particle model for low-density granular flow based on the equation of motion for particles. Realistically, the granular flows of natural phenomena are roughly affected by physical factors, namely, gravity and basal and inner frictions. The motion of each particle in the flow is also affected by the particle size and turbulence intensity. Thus, the dynamics of granular flow are complex, and understanding the mechanism becomes difficult. Our model aims to significantly simplify the complicated behavior of low-density granular flows interacting with the ambient fluid. Although our model is not realistic at this time, we aim to reproduce the flow patterns found in previous experiments with ping pong balls and polystyrene particles [20,23].
As a simple case, we consider that a low-density granular flow slides on a slope at an angle steeper than the angle of repose ( Figure 1). This model is based on three basic assumptions for the simplification of the modeling: (i) The granular flow consists of hard spherical particles composed of a material of low specific density. (ii) Only the translational motion of the particles is considered (i.e., the rotational motion is ignored). (iii) Three types of forces act on the particles: (a) gravity as the driving force for granular flow, (b) the repulsive force between the particles due to particle collision, and (c) the drag force due to the ambient fluid. Figure 1 shows the schematic of this model, where θ is the incline angle, the x-axis is the inclination direction, the y-axis is the lateral direction, and the z-axis is perpendicular to the slope. Figure 1. Schematic of a simple particle model. The x-y plane is parallel to the slope with the incline angle of θ (x: inclination direction; y: lateral direction). The z-axis is perpendicular to the slope.

Particle Model
To derive the governing equation for the particles, we formulate the force acting on the ith particle is the repulsive force between the particles, and F d i is the drag force due to the ambient fluid. Hereinafter, the radius, coordinate, and velocity of the ith particle are denoted by a i , r i , and v i , respectively, and the true density of the particles is denoted by ρ p .
Gravity F g : The gravity in the vertical direction of the ith particle F g i is calculated by considering the density difference between the particle and fluid. Subsequently, it is divided into the x and z components depending on the incline angle θ, as follows: where e θ is the transformation vector from the standard coordinate system into the incline coordinate system, as shown in Figure 1, V i is the volume of the ith particle, g is the gravitational acceleration, and ρ f is the density of the fluid. Repulsive force F r : Because the rotational motion of the particles is ignored, only the normal force of collision between the ith and jth particles is calculated as the repulsive force of the ith particle F r i (and F r j ), which is represented by the elastic spring force F r i as follows, Here, r ij = |r j − r i | is the interparticle distance and n ij is the unit vector joining the centers of two particles. Further, k n is the linear spring constant and δ ij is the overlap between two particles. For the collision between a particle and the slope, we calculate the overlap distance between them. Then, the normal repulsive force acting on the particle is determined using Equation (2b). Drag force F d : In this model, we estimate the drag force from the fluid velocity field generated by the particle motion. Concretely, the ambient fluid is assumed to be Stokes flow (i.e., we assume a low Reynolds number; Re 1), although the particle inertia is not negligible and the turbulent eddy occurs in the actual granular flow. The Stokes flow containing two spherical particles can be solved theoretically when a sufficient distance exists between two particles. The exact solution is known as the Rotne-Prager tensor J(r) [35]. This tensor can solve the fluid velocity field based on the velocity of the jth particle. The fluid velocity at the coordinate of the ith particle is added to the velocity of the ith particle as the induced velocity, because the particle completely follows the fluid flow in Stokes flow. Using Stokes' drag (Equation (3a)) and the Rotne-Prager tensor (Equation (3c)), we express the interparticle interaction through the ambient fluid, that is, the drag force of the ith particle F d ij is generated by the force acting on the jth particle F g+r j = F g j + F r j , as follows: In these equations, µ is the viscosity coefficient of the fluid, u i (j) is the velocity induced by the jth particle at the coordinate of the ith particle (Figure 2), and I is the unit tensor. Particularly, rr in Equation (3c) denotes the tensor corresponding to I. For instance, rr xy = (x i − x j )(y i − y j ). Induced velocity field u(j) generated by gravity and the repulsive force of the jth particle, F g+r j , using Equation (3b). According to the force direction, the ambient fluid is pushed and pulled as shown by the arrows. The induced velocity at the ith particle is expressed as u i (j).
According to the assumption for the drag force (i.e., Stokes flow), we derive the governing equation for the ith particle from Stokes' drag as follows, where N is the number of particles. In Equation (4a), the second term indicates the long-range interaction between particles through the ambient fluid, which resembles the asymmetric potential ( Figure 2). A limitation of this model is that the Rotne-Prager tensor overestimates the interaction between particles when they are close to each other (its condition is out of the theory [35]). The simulation of granular flow with this tensor is not realistic, but it could be helpful to understand the time evolution of granular flow strongly affected by the drag force.

Simulation Setup
The numerical modeling for the low-density granular flow is performed in a three-dimensional space, as shown in Figure 1. However, we use the time complexity of O(N 2 ) in our model to calculate the drag force F d ij in Equation (3). The repulsive force F r i in Equation (2) additionally requires the small timestep dt used in the numerical simulation. Therefore, this model is unsuitable for calculations involving a large number of particles. In this study, the number of particles is set as N = 2000. Subsequently, numerical simulations are conducted in two-dimensional planes: the x-y plane parallel to the slope and the x-z plane perpendicular to the slope. Figure 3 shows three different types of initial conditions for the particles in the x-y and x-z planes: (a) circular shape, (b) rectangular shape, and (c) triangular shape. Here, the incline angle is fixed as θ = 45 • . These shapes are determined by reference to previous experiments with polystyrene particles [23] although the size is small in comparison to the work in [23] because of the small N. The circle radius R, rectangle width W, and triangle depth D are set as constant values depending on the particle size used in the simulations, as shown later. Particularly, in the case of Figure 3b, the length of the rectangle in the x direction is fixed as 10 particles. In each shape, 2000 particles are packed randomly without overlap. The x-y plane is parallel to the slope, whereas the x-z plane is the vertical cross section on the slope, where θ = 45 • is the incline angle. In each initial condition, the circle radius R, rectangle width W, and triangle depth D are constant depending on the particle size, such that the particle volume fraction Φ p of 2000 particles is approximately 0.5.
Regarding the specific materials of the particle and fluid, we assume that the granular flow consists of low-density particles, such as polystyrene in air, to allow a qualitative comparison with the previous experiments [23]. The particle density and linear spring constant for the repulsive force are given as ρ p = 20 kg m −3 and k n = 10 N m −1 , respectively, whereas the density and viscosity coefficient of the fluid are ρ f = 1.2 kg m −3 and µ = 1.82 × 10 −5 Pa s, respectively. We determine the value of k n so that the overlap distance between particles can be small. As the control parameter, we vary only the particle radius, such that a p = 1, 2.5, and 5 mm (we consider three values). To avoid the crystallization of particles, particles with a polydispersity of ±5% are used in simulations.
The sizes of the initial shapes R, W, and D, as shown in Figure 3, were determined depending on a p = 1, 2.5, and 5 mm: R = 0.065, 0.17, and 0.34 m; W = 0.35, 0.85, and 1.7 m; and D = 0.115, 0.3, and 0.6 m. These sizes normalized by a p areR = 65-68,W = 340-350, andD = 115-120, respectively. Further, the particle volume fraction is approximately 0.5. We perform 20 simulations by changing the initial position of the particles foreach initial shape and each particle radius.

Nondimensionalization of the Proposed Model
To provide consistent interpretation of numerical simulations, we conduct the nondimensionalization of our model before the results. First, we derived the speed of the single particle v s in the case of N = 1; that is, this is obtained by referring to gravity alone. The particle coordinates r, particle velocity v, and time t are normalized using the particle radius a p and v s as follows,r Next, the governing equations for the ith particle (Equation (4)) are rewritten as follows, whereF g+r i =F g i +F r i is the sum of the dimensionless gravity and dimensionless repulsive force of the ith particle, andũ i (j) is the dimensionless induced velocity generated by the jth particle at the coordinate of the ith particle. Furthermore, the dimensionless drag forceF d In this equation, because the dimension of the Rotne-Prager tensor J(r) is the inverse of r, we obtain the equality J(r) = J(r)/a p . The dimensionless parameter is include inF g+r i in Equations (6a) and (7b). Finally, the dimensionless gravityF g i is expressed using Equation (1) as follows, where θ is the incline angle of the slope, and e θ is the vector defined in Equation (1a). The dimensionless repulsive forceF r i is also expressed using Equation (2) as follows: where γ is the dimensionless parameter in our model. The physical meaning of γ is the ratio of the excluded volume effect between particles divided by gravity. We substitute the values used for the simulations in γ and assume δ ij = 0.01a p . Then, γ is calculated as 130, 20, and 5 (a p = 1, 2.5, and 5 mm). It means that the repulsive force is larger than gravity in our simulation.

Circular Shape in the x-y Plane
The typical time evolutions of the granular flow in the circular shape are shown in Figure 4. Here, the particle radius is a p = 5 mm. In the early stage of the simulations, the particle located at the central part of the circle migrates rapidly to the moving front of the granular flow as time progresses. Consequently, the initial uniform distribution of particles shifts to the local distribution at the moving front of the granular flow (Figure 4a,b). This pattern of particle aggregation appears to be a single head resembling a crescentic form, similar to the pattern of a low-density granular flow [20,23]. After the formation of the head, the particle moves away gradually from the aggregated part. Subsequently, the head shrinks with time (Figure 4c,d). This head formation process is also confirmed for the other particle radiuses although the time required for the head formation is different. To verify the key factor for the head formation, we perform the simulation without the drag force. In this case, the initial pattern is maintained because the basal friction and particle inertia are ignored in our model. Therefore, the drag force generates a long-range interaction between particles through the velocity induced by the ambient fluid u i (j) in Equation (3), which is essential to form the head. The reason for the head formation, namely owing to the induced velocity, is clear. In the initial stage of the simulations (Figure 4a), the induced velocity is high at the central part because it is inversely proportional to the interparticle distance. Therefore, particle aggregation occurs in the moving front of the granular flow, and the aggregation migrates at a faster speed.
For the migration speed of the granular flow, we measure the frontal speed v f from the displacement in the x direction of the forefront particle position. In our model, v f is defined as follows, where dt is the time step. Consequently, the frontal speed is normalized by the speed of the single particle (Equation (5b)) asṽ f = v f /v s . Because v s increases with a p , the smaller time step dt is required for the calculation of the repulsive force in the larger a p . The numerical simulations continue until the granular flow reaches the state shown in Figure 4d: t end = 20, 10, and 5 ms (a p = 1, 2.5, and 5 mm). Figure 5a shows the temporal variations in the normalized frontal speedsṽ f (t) measured at different particle radii a p . The solid curve and thin band denote the mean and standard deviation of 20 simulations, respectively. The standard deviation is small compared to the mean even though the size distribution and initial position of the particles are changed. In all cases,ṽ f increases significantly in the early stage and subsequently decreases gradually with time because of the head shrinkage.
Additionally,ṽ f ranges from 30 to 70, independent of a p . Therefore, the head size (i.e., number of particles) causes the change inṽ f (Figures 4 and 5a). The head size is discussed in detail in Section 4.2. The movements of the particles in the head can be visualized by their respective velocity vectors (Figure 5b). This particle location corresponds to that in Figure 4c; that is,t = 19.5 and a p = 5 mm. To study the particle velocity relative to the frontal speed of the granular flow v f , the x component of the vectors is denoted by v x − v f , where v x is the x component of the particle velocity. Consequently, a granular vortex convection is observed in the head as well as the particles away it. In particular, the particle at the rear of the head moves inward to the forefront of the head, whereas the particle at the front of the head is pushed outward. Some particles return to the head again. However, most of the particles are left behind. Thus, the head shrinks with time.

Rectangular Shape in the x-y Plane
The typical time evolutions of granular flow in the rectangular shape are shown in Figure 6a-c. Here, the particle radius is a p = 5 mm. The simulation time at each a p is the same as that for the circular shape: t end = 20, 10, and 5 ms (a p = 1, 2.5, and 5 mm). In the early stage of the simulations, the initial straight shape at the moving front of the granular flow is maintained although the particle moves forward slowly to the front due to the long-range interaction between the particles via the drag force. Subsequently, the straight shape with particle aggregation deforms into a wavy pattern (Figure 6b); that is, flow instability occurs. This destabilization of the moving front is caused by the inhomogeneous distribution of the size and coordinates in the initial condition of the particles (Figure 6a). Because the aggregation is faster than the other parts, the wavy pattern elongates in a downward direction. Eventually, the multihead structure is formed from the straight shape (Figure 6c).  Figure 6d shows the front speed in the x direction normalized by the speed of the single particle, v f (t), at each particle radius a p . The mean and standard deviation of 20 simulations are denoted by the solid curve and thin band, respectively. In all cases,ṽ f decreases gradually in the early stage. Subsequently, the decreasing rate ofṽ f tends to be relatively high. Additionally,ṽ f ranges from 20 to 28 independent of a p , and the formation of the multihead structure is confirmed at the end of each simulation (Figure 6c). The formed pattern (number of heads and head size) differs depending on the initial condition. However, the wavelength reproducibility is good. This point is discussed in detail in Section 4.2.
In Figure 6d, the rate of decrease inṽ f is linked with the flow pattern, as shown in Figure 6a-c. For example, during the deformation from the straight shape (Figure 6a) to the wavy pattern (Figure 6b), the number of particles constituting the aggregation remains nearly unchanged. Thus,ṽ f is maintained in the early stage of the simulations. Simultaneously, the number of particles decreases locally with the elongation of the heads from the wavy pattern. Consequently,ṽ f decreases at a relatively high rate after the middle stage of the simulations.

Triangular Shape in the x-z Plane
The typical time evolutions of the granular flow in the triangular shape are shown in Figure 7a,b. Here, the incline angle and particle radius are θ = 45 • and a p = 5 mm, respectively. The simulation time at each a p is set as t end = 20, 10, and 5 ms (a p = 1, 2.5, and 5 mm). Because the particles are up in the air initially (Figure 3c), they immediately accumulate near the slope according to the gravity. Subsequently, a large-scale cluster consisting of particles emerges from the accumulation, and moves down a slope at a high speed (Figure 7a). This pattern appears to be the head-tail structure, where the thickness in the z direction of the front part is more than twice that of the rear part. This profile is similar to the flow height of the experiments with ping pong balls [19]. After the head-tail structure is formed, the head shrinks with time ( Figure 7b). In the x-z plane, the front speed in the x direction v f is measured according to the displacement of the forefront particle on the slope (i.e., the first layer). Figure 7c shows the normalized front speedṽ f (t) at different particle radiuses a p . The solid curve and thin band denote the mean and standard deviation of 20 simulations, respectively. In each case,ṽ f increases significantly with the cluster formation in the early stage. Subsequently,ṽ f decreases with time because of the head shrinkage. The standard deviation is negligibly small after the head-tail structure is formed. This implies that the identical structure is formed independently of the initial condition.
To verify the particle movement in the head-tail structure, we focus on the x component of the particle velocity. In Figure 7a, the speed of the white particle is slower compared to the front speed v f , whereas that of the gray particle is faster compared to v f . The particle speed in the front part of the granular flow (except the upper part of the head) exceeds v f , whereas that at the rear part (i.e., the tail) is smaller than v f . Therefore, the particle moves from the upper part of the head to the tail. These particle movements lead to the head shrinkage.

Definition of the Head in the x-y Plane
As mentioned in Sections 3.1 and 3.2, the front speed of the granular flow relates to the head size. The multihead structure is formed from the straight shape in the x-y plane, which is parallel to the slope. Hereinafter, we define the head based on the particle number density to discuss the temporal variation of the head size and the size distribution of the head.
Our idea is to detect the particle with the relatively high particle number density from the flow pattern, as shown in Figure 8. As a simple index, we count two different types of numbers around the ith particle, n i and n´i as follows, where r ij is the interparticle distance between the ith and jth particles. Here, n in Equation (11a) is the number of particles around the ith particle, whereas n´indicates the number of indirect neighbors of n.
In other words, n´yields the quasi-local particle number density. According to the value of n or n´, the particle constituting the head is detected, as shown in Figure 8a,b. The procedure for calculating n (or n´) is described below. The value of n depends on the aggregation size (i.e., the head). Therefore, we set the criterion for judging the constituent particle of the head using the maximum value n max calculated for each aggregated part. If the ith particle satisfies n i /n max > 1/4, this particle is assumed to be the constituent of the head. We apply the same procedure to n´. Figure 8a,b shows the head defined according to n and n´, respectively. The blue particle satisfies the criterion above; that is, it is a constituent of the head. Compared to the criteria for n and n´, the difference is obvious. The head extends backward along the outside edge of the granular flow in the case of n (Figure 8a), whereas the head is located at the front of the granular flow in the case of n (  Figure 8b). In this study, we adopt the criterion with n´as the definition of the head.
To measure the head size, we select four targets from the constituent particles of the head, as shown in Figure 8c. Two green particles are selected as the outmost positions in the y direction, and the others are selected as the forefronts of the outer and inner positions in the x direction. Using these particles, the width and layer thickness of the head, w h and l h , respectively, are measured (Figure 8c). The calculations for the temporal variation in the head size and the head size distribution are explained in Section 4.2.

Characteristics of Head Size in the x-y Plane
Based on the quasi-local particle number density n´in Equation (11b), we study the temporal variation of the head size from the granular flow in the circular shape (see Section 3.1). Because the simulation time is different for each particle radius a p , the width and layer thickness of the head, w h and l h , respectively, are measured at different time intervals depending on a p : ∆t = 1, 0.5, and 0.25 ms (a p = 1, 2.5, and 5 mm). Therefore, the number of data is fixed as 20 in each simulation. Figure 9a shows the time evolutions of the number of particles constituting the head defined by the criterion n in Equation (11b). The point and error bar denote the mean and standard deviation of 20 simulations, respectively. In the early stage, the head consists of almost half of the total number of particles (N = 2000), but the number of constituent particles rapidly decreases with time. Figure 9b shows the time evolutions of the normalized width and thickness of the head,w h (t) andl h (t), which are denoted by the open and filled points, respectively. The normalized width and thickness decrease with time, leading to the shrinking of the head. Also, these decreasing trends ofw h andl h are very similar to that of the number of particles constituting the head. To evaluate the relationship between the head size and the frontal speed, as shown in Figures 4 and 5a, we compare the number of particles constituting the head with the normalized frontal speedṽ f (Figure 9c). The frontal speed of the granular flow becomes slower with the decrease in the number of particles constituting the head. Therefore, we find a good positive relation between the head size and frontal speed. Next, we focus on the characteristics of the head shape in the circular and rectangular shapes. Figure 10a shows the relationship betweenw h andl h . The point and error bar denote the mean and standard deviation of 20 simulations, respectively. During the head shrinkage, the ratio between width and layer thickness is maintained (w h /l h = 4), except during the early stage. We find that the similarity law for the head shape holds in the low-density granular flow. Moreover, we verify the head size distribution on the multihead structure of the rectangular shape (see Section 3.2). In this analysis, the head size is measured at the end of the simulations, where the head elongates downward, as shown in Figure 6c. In Figure 6c, the blue particles denote the constituent of the head, and multiple heads are detected. The measurement of the layer thickness l h is difficult in this situation. Therefore, only the width w h is measured as the head size. For the number of heads, the mean and standard deviation of 20 simulations at a p = 1, 2.5, and 5 mm are (17.20, 1.69), (13.95, 1.94), and (14.85, 1.35), respectively. Figure 10b shows the probability density of the normalized widthw h . Here,w h is distributed symmetrically and the distribution shape is independent of a p . This distribution is fitted well by the Gaussian distribution indicated by the dashed black line in Figure 10b. Although the fluctuation inw h might be caused by the inhomogeneity of the initial condition, the flow instability at the moving front of the granular flow exhibits the characteristic wavelength. For the width of the head w h , we compare our result qualitatively with the previous experimental findings using polystyrene particles [23]. In the paper, flow instability resembling a wavy pattern was observed, and the frontal radius was estimated by fitting the circle against each head (see Table 1 in [23]). In this study, we assume that w h corresponds to the diameter of the frontal shape (i.e., 2×radius) in the previous experiments. Figure 10c shows the relationship between particle radius a p and w h in this study and that in the previous experiments. The point and error bar denote the mean and standard deviation, respectively. In particular, the data in this study are calculated from the distribution shown in Figure 10b. Although the measured w h in this study is much smaller than that in the experiments because of the small number of particles in the two-dimensional planes, each set of data fits the linear function well (dashed line in Figure 10c). We find that w h is proportional to a p in both cases.

Frontal Angle in the x-z Plane
As mentioned in Section 3.3, the low-density granular flow is simulated in the x-z plane and the head-tail structure is formed. The shape of the head was observed from the side view of the granular flow in previous experiments with ping pong balls or polystyrene particles [20,25,32]. Notably, the frontal angle of the head generally exhibits a high elevation angle to the slope and then becomes a constant angle (60 • ) according to the kinematic theory (see in [32]). Therefore, we measure the frontal angle of the head from the flow pattern, as shown in Figure 7b,c, to check if the result agrees with this fact.
The definition of frontal angle φ f is based on the coordinate of the particle located on the surface of the head. Figure 11a shows the typical shape of the head and the corresponding φ f . To measure φ f , we first select the forefront particle of the granular flow in the x direction, namely, the green particle in Figure 11a. Next, we search the constituent of the head surface around the green particle, as shown in Figure 11b, where the search range spans the three particle radii and the search angle against the slope is ϕ = 0 • -180 • . By repeating this search, the red particles in Figure 11a are detected. A fitting line is obtained from the coordinates of the green and red particles using the least-squares method.
Subsequently, φ f is estimated. It is noteworthy that the blue particles below the green particle are excluded from the estimation. (a,b) The green particle is at the forefront of the granular flow, and the neighbor particle is selected within 3a p and ϕ = 0 • -180 • . The red particle is selected in a similar way, and the green particle contributes to φ f , whereas the blue particle is excludable. The dashed black line is drawn using the least-squares method against the coordinates of the red and green particles. (c) The mean and standard deviation of 20 simulations are denoted by a point and error bar, respectively. Figure 11c shows the temporal variation of the frontal angle φ f at each particle radius a p in the triangular shape. The mean and standard deviation of 20 simulations are denoted by the point and error bar, respectively. In the initial condition, as shown in Figure 7a, φ f is estimated as approximately 90 • . Subsequently, φ f approaches 60 • from the larger angle as time progresses, and the standard deviation of φ p decreases with time in all cases. This convergence angle of 60 • is consistent with the frontal angle reported in previous studies [32]. However, φ f = 60 • in this study results in the crystallization of particles in the two-dimensional simulation (Figure 7b,c), although we use particles with a polydispersity of ±5% to avoid this phenomenon. Even though a three-dimensional simulation is performed, such a crystallization is to be expected.

Conclusions
In this study, three types of pattern formations for a low-density granular flow were reproduced qualitatively using the single particle model: a crescent-shaped head, a wavy pattern due to flow instability, and a head-tail structure. This model considered the gravity, repulsive force, and drag force due to the fluid as forces acting on the particle, thus providing the minimum number of factors to simplify this natural but complicated phenomenon. The frontal speed of the head decreased with time and was related positively to the head shrinkage (or number of constituent particles). Moreover, granular vortex convection was formed at the head in the x-y plane parallel to the slope, whereas the particle moved from the upper and rear parts of the head to the tail in the x-z plane perpendicular to the slope. According to the quasi-local particle number density, the width and layer thickness of the head were defined as the characteristic head size. The ratio (width/layer thickness) was maintained (≈4) while the head shrunk, and the size distribution of the width obeyed the Gaussian distribution. For a comparison with previous experiments regarding the low-density granular flow, we found a linear relationship between the head width and particle radius. Additionally, the frontal angle of the head converged on 60 • . This model is not realistic yet. However, we believe we pioneered a particle-based approach to elucidate the dynamics of low-density granular flow in this study.