Simulation of Solid Particle Interactions Including Segregated Lamination by Using MPS Method

A new MPS (Moving Particle Semi-implicit) method is developed to simulate the behaviors and interactions of multiple fine solid particles as a continuum. As fluid particles are affected by viscosity, so solid particles are affected by friction. The solid particle dynamics for landslides, dumping, and gravity sorting etc. which can be difficult to simulate using conventional MPS methods, are modeled in this paper using the developed multi-solid-particle MPS method that benefits from drawing comparisons with the corresponding fluid particle behaviors. The present MPS results for dumping solid particles are verified against the corresponding DEM (Discrete Element Method) results. The shape and angle of repose for solid particles are shown to be highly dependent on the friction coefficient between grains. The peculiar phenomenon of segregated lamination (gravity sorting) among grains of different densities has been successfully reproduced using the multi-solid-particle MPS method. Lamination quality is found to be dependent on the densities and frictional coefficients of the constituent particles. The behavior of heterogeneous mixtures of multiple solid and liquid particles are also compared and discussed. This newly developed tool offers a window into the physical dynamics of sedimentology that the broader geoscience community might benefit from.

1 Introduction CFD (Computational fluid dynamics) has benefited greatly over the past decades from improved algorithms and hardware capabilities. Three approaches to the field of CFD have emerged. The conventional Eulerian approach employs a fixed-grid system as its computational domain. The Lagrangian approach is grid-less and assigns physical properties to the representative particles. The third approach is a hybrid method that employs both Eulerian and Lagrangian system. There are advantages and disadvantages to each.
problems. The resulting stack shapes, which are dependent on friction coefficients and initial conditions, are observed. For the verification of the developed MPS-based method, another numerical simulation by a widely used DEM (discrete element method) computer program was carried out under the same conditions of the dumping problem. The DEM results agree well with our MPS results. DEM is the well-established method for analyzing microscopic behaviors of discontinuous media. Zhou et al. [Zhou, Xu, Yu et al. (2002)] and Chen et al. [Chen, Liu, Zhao et al. (2015)], for example, performed extensive parametric studies for the stacking shape of similarly dumped mono-sized coarse spheres by DEM and the DEM results were validated by their experiments. Finally, stratification by gravity sorting of multiple solid grains of different density is simulated using the developed MPS with various friction coefficients. The above results are compared with the corresponding multi-liquid simulations. For the purposes of this study, which concerns itself with the behavior of fine, loose grains, the consideration of other soil parameters, such as elastic modulus, Poisson's ratio, and undrained shear strength, is not significant and thus disregarded.

MPS for solid particles
The MPS method is based on a Lagrangian approach. It was designed originally to address the kinds of fluid dynamics problems such as presented in "broken dam," sloshing, and similar scenarios. It employs the (1) where r is the density, µ is the dynamic viscosity, t is time, u  is the fluid velocity, P is pressure, σ is the surface tension coefficient, κ is the curvature for surface tension, n  is the unit normal vector of interface, F  is the external force per unit volume including gravitational force, ∇ is the gradient, and 2 ∇ is Laplacian operator. Though the MPS approach mainly has been used in fluid dynamics, it can be modified for applications in solid particle interaction such as landslides and dumping problems. Compared to liquid particles, which uses the terms of viscous diffusion and surface tension on the right-hand side, solid particles as a continuum use terms of frictional force and drag force.
where f F  is the frictional force and D F  is the drag force. The drag force can be applied to particles exposed to ambient flow; otherwise D F  is zero. The forces in the right hand side of (3) are for unit volume. In the present study of solid particle interactions in air, the drag force can be disregarded. The normal forces appearing in the DEM formulation are indirectly accounted for by employing a proper collision model. While DEM actually models all the particles, the present MPS method assigns the representative physical properties of medium at dense collocation points in a weighted manner. In order to account for interaction among particles, a special numerical treatment, which weighs the influence of a central particle on neighboring ones through kernel function up to re, is required. In this way, a kernel function is employed to measure the influence on neighboring particles.
e e e r r r r w r r r r r where subscription ij r denotes ( ) j i r r − , i r is the position of central particle i , and e r is the effective range. In this study, e r is set to 0 2.1 l × , where 0 l is the initial particle distance [Lee, Park, Kim et al. (2011)]). According to Eq. (4), by definition when the distance between the center particle and its neighbors exceeds its effective range, the neighboring particles have no influence on the central one because the magnitude of the kernel function is reduced to zero. Another special way of treating particles is by particle number density. The summation of kernel functions within an effective range can be regarded as fluid density. The combination of the kernel function and the concept of particle number density allows the MPS system to be treated as a continuum. The particle number density can be expressed by: where ij r  is the particle distance vector between center and neighboring particles. The first term on the right-hand side of Eq. (3) denotes the pressure gradient which represents a local weighted average of the gradient vectors acting between the central particle i and neighboring particle j . The gradient model can be expressed as follows: where denotes the particle interaction model, φ is an arbitrary physical quantity, d is dimension number, r  is the particle's position vector, and 0 n is the initial particle number density. ], for example, gave further details that can be applied to general MPS formulations including pressure calculations. One of the issues in using the Lagragian MPS approach is the estimation of force when particles collide. As particles get close, the overlapping range of each particle can introduce numerical errors. Therefore, a proper collision model is used to adjust the repulsive force and to avoid a sudden increase in particle number density. When the particles are initially arranged, they are uniformly distant. As they move with a velocity, the distance between them can change and a collision model is activated once their distance goes below a threshold value of 0 a l × where a is a constant. After applying the conservation of momentum, the repulsive velocity can be estimated using restitution coefficient b . The parameters a and b can be assigned using numerical tests, and in this study they are set at 0.97 and 0.2, respectively. The third-term of the right-hand-side of Eq. (3) denotes the frictional force on a central particle. The frictional force can be measured using the normal force and frictional coefficient. One of the normal force elements is inertial force due to acceleration. For that, it is necessary to calculate the relative acceleration between center and neighboring particles in the normal direction. Then, the corresponding frictional force can be determined as shown in Fig. 1(a). The frictional force due to particle acceleration can be calculated as follows: where 0 l is the initial (or starting) particle distance, µ denotes the frictional coefficient, n a is the normal component of acceleration for each particle, and ij t  denotes the unit tangential vector of the frictional force. An additional component of the frictional force arises from the accumulated weight of particles above the central particle. Since the particles are stacked, the effect of the accumulated weight can be compared to hydrostatic pressure, and should be considered. As shown in Fig. 1(b), particles above a central particle are accounted for in the calculation of accumulated weight. To avoid overestimation, a term for the specific range, or x-interval, is introduced. When the neighboring particles above a central particle are out of the x-interval range, it is excluded from the accumulated effects. The accumulated weight can be calculated as follows: where , h i p is the accumulated weight of stacked particles, , s ij p is the static weight averaged from the accumulated weight acting on a central particle from its neighbors, f r is the effective range of the central particle in the frictional force calculation, and x z is the particle's position on x-and z-axes, respectively. In order to avoid overestimation of the stacked effects, the f r is set to 0.65, which can be adjusted as needed. With a calculated accumulated weight, the frictional force due to static weight can be obtained from: where ij θ is given in Fig. 1.
By combining Eqs. (7) and (10) the total frictional force on a central particle can be calculated as follows: Frictional force resists forward movement and can only reduce or stop particle motion. Particle motion begins when the inertial or gravitational forces exceed the force of friction. In this case, static coefficients of friction are used. Once particles move, dynamic friction coefficients may be smaller than static ones. However, for simplicity in this paper, the differences between static and dynamic frictional coefficients are assumed to be minimal, and are disregarded.

Applications and numerical examples
The MPS method developed here is next applied in three separate scenarios. The brokendam, cargo-dump, and lamination by gravity sorting. In each application, the method lends insights into the physical behavior of particles. Generally, a liquid will readily spill out onto a surface due to the absence of significant frictional forces acting between particles. This is not the case for solid particles. Because of this, stacked solid particles tend to assume various forms such as that of a mountain, dune, etc. The angle of repose for the stacked particles will have a range depending on the particle's properties. In geoscience applications, the coefficient of friction between particles is regarded as foremost among these properties [Clover (1995)].
The following "dam break" scenario compares the behavior of water with that of solid particles (Fig. 2). The container is 0.6 m high, 0.6 m long, lacks a lid, and the rectangular volume is restrained on the right-hand side by a an infinitely thin gate. The materials, whether water or solid particles, have a height of 0.3 m and length of 0.15 m. When the gate is opened at infinite speed, the force of gravity acts on the material. Fig. 3 depicts snapshots taken at various times measured in seconds following the release of the contents. In the first case, the column of water collapses rapidly and collides forcefully against the far wall. The momentum of the moving water generates an upward-directed jet that shoots up the wall. Water behavior simulated in this case (Fig. 3a) was validated against experimental, and other CFD results in authors' previous papers e.g. Lee et al. [Lee, Park, Kim et al. (2011)], Bakti et al. [Bakti, Kim, Kim et al. (2016)]. The fastmoving water took less than 0.4 s to traverse the tank and slam against the far wall.
Eventually the water surface settled to a planar surface parallel to the bottom of the tank. Solid particles under the hindered influence of friction behave quite differently (Fig. 3b), given the same starting conditions. When the barrier is removed, the grains roll down and pile up against one another to form a slope under the influence of friction forces. The development and movement of the free-surface are much slower compared to the water case. The final dune-shape is reached only after the gravitational force acting on each particle no longer exceeds the frictional force. The dune-face is characterized by an angle of repose, which is not quite planar but steepened upward. The "dam break" simulation above (shown in Fig. 2) is a classical problem for the particle-based Lagrangian CFD solution. The total number of particles used is 4,300 and 1,800 for fluid and solid particles, while others are used for the wall components. The properties of fluid and solid particles are given in Tab. 1.  The "dam-break" example shows that both simulations for water and fine solid grain reasonably correspond to the expectation from physics and typical observation in nature.
Next, let us move on to "dumping" problem. The dumping problem shown in Fig. 4 gives a second comparison of fluid and solidparticle behavior. In this simulation, the bottom of a container suspended 17 cm above the floor of a rectangular tank is opened suddenly and its contents spilled out under the force of gravity. For this simulation, 2,938 particles were used for the dumped material and 3,652 particles were used for wall and dummy particles. Water properties and solidparticle properties are the same as those used in the previous example. The total time of the simulation was 5 s and calculations were performed at 0.001 s interval. Results from the dumping-problem are summarized in Fig. 5. In the case of water, exit from the container is rapid and the vertical momentum upon striking the floor of the tank is converted into horizontal flow in both directions. Backwash from waters colliding with tank walls generates overturning breakers that eventually lose momentum and settle to an equilibrium state. In contrast, solid particles dropped from the same simulated apparatus do not travel far from their first point of impact due to the force of friction. Subsequent particles falling on the loose pile possess diminished kinetic energy as the simulation progresses due to the growing height of the pile and reduction in fall distance. Yet this diminished kinetic energy is offset in part by the enhanced ability of particles to roll, slide and bounce down the pile flanks as the stack slope increases. In the end, the particle stack resembles a dune of uniform slope. Equilibrium is reached once the dune height reaches the mouth of the opening.

Figure 4: Schematic view for cargo dump problem
Runs were made using particles of various coefficients of friction, as seen in Fig. 5. Note how the angle of repose is dependent on the coefficient of friction. At lower coefficients such as 0.1, the dumped material spreads more widely, allowing for the chamber to empty all its contents without clogging. The relationship can be stated as follows: where θ is the angle of repose, and i µ is the static frictional coefficient for particle i .
According to Clover [Clover (1995)], the range of repose angle of the same soil can vary due to initial conditions. For example, comparing Figs. 6(a) and 6(b), which use particles identical in density and coefficient of friction, yet because the simulations differ in initial condition, the resulting angle of repose becomes different. Fig. 6(b) shows measured repose angles with various friction coefficients.  Fig.7 is compared with Fig. 5, the result in the case of the raised container corresponds to a stack that no longer reached the height of the release point; yet a small portion of grains of sufficiently high coefficient of friction remained inside the container. Again, the higher frictional coefficients result in steeper angles of repose, like in the previous case. Notice also, that at the lowest frictional coefficient of 0.1, the cusp of the stack is more rounded and the stack more spread out due to the higher kinetic energy of the falling particles. Fig. 8a compares the two stack shapes. Again, the "dumping" example shows that both simulations for water and fine solid grain reasonably correspond to the expectation from physics and typical observation in nature although the current MPS method does not model all the details of micro-scale particle shapes, sizes, and properties.  [Zhou, Xu, Yu et al. (2002)] and Chen et al. [Chen, Liu, Zhao et al. (2015)], their DEM simulations were compared with their experiments and they concluded that the DEM simulation methods can reasonably reproduce the experimental results and thus are adequate in producing the relevant physics. So, in order to validate the developed methodology, we compare the sample case of the present MPS simulation with the DEM simulation using the same sets of parameters.
For the DEM simulation, the open-source software, LIGGGHTS-PUBLIC was used. According to this program, the particle's motions are updated in time domain by calculating momentum equations including inertia and interaction forces with neighboring particles or walls, which is very similar to the present case. The damped Hertzian contact model was employed to estimate interaction forces between particles which can be described by the following equation: where ij n δ is the overlap distance of two particles, i , j , ij t δ is tangential displacement vectors, v ij n and v ij t are normal and tangential components of the relative velocity of two particles, n k and t k are normal and tangential stiffness, and n γ and t γ are viscoelastic damping constants for normal and tangential contacts, respectively. Therefore, the first and second brackets represent the normal and tangential forces respectively. The parameters used to find damping and stiffness coefficients are tabulated as shown in the following table. Fig. 8b shows the DEM results for the same cases of Fig. 8a.   When the individual modeling of DEM is used for very fine particles, the CPU time will greatly be increased compared to SPM. When fine particles move with ambient fluid motions, the present MPS can straightforwardly model the solidfluid combined case while DEM can hardly be employed for this kind of application.

Segregated-lamination simulation by multi-particle MPS
The phenomenon of segregated lamination by gravity sorting is next simulated using the same MPS method, with the added complication of a heterogeneous mixture of particles. Previous simulations used only homogenous particles. Interesting experiments of segregated lamination using mixed heterogenous sands dropped in air from a container were reported in Julien et al. [Julien, Lan, Berthault et al. (1993, 1994]. Fig. 9 is a photo of results from one of these lab experiments using sands that were bimodal in density distribution. Julien et al. [Julien, Lan, Berthault et al. (1993, 1994] reported clearer lamination when using a mix of two sands that differed in both density and grain diameter. Bouncing and rolling distances of the two grain populations were not the same, resulting in well-defined lamination, which improves with greater rolling distances. The gravity sorting of heterogeneous grains into laminae has also been documented in a water flow in a laboratory water flume. This finding could have important implications for understanding sediment behavior, including the rapid accumulation of vertically stacked laminae under flood-conditions. Figure 9: Experimental lamination of a heterogeneous mixture of two sands that differ in density and grain-diameter In this regard, the solid MPS method is next applied to investigate segregated lamination in air, using the same numerical apparatus as those from previous "cargo dump" runs. Three particle populations, each different in density yet identical in size, were selected. Initially, three different particles are arranged by 6 layers. When the contained moth is suddenly opened, the three particles are mixed and fall simultaneously as a heterogeneous mixture. The mixed particles fall and accumulate as a stack, gravity sorting (segregated lamination) of grains takes place spontaneously. The simulated results strongly resembled Julien et al.'s [Julien, Lan, Berthault et al. (1994)] lab results. Fig. 10 shows starting conditions for the simulation of three populations of particles arranged in six layers. The particles differ in density and frictional coefficient as listed in Tab. 3 but are identical in size. Fig. 11 depicts two snapshots in the progress of the run, and Fig. 12 shows the two snapshots for sorted individual grain population. In Fig. 12 for 5 s snapshot, we can clearly see that all 3 particles are falling simultaneously as a heterogeneous mixture. If it were three specific liquids with the three different densities, the eventual separation by 3 layers is caused by buoyancy forces. Solid particles behave much differently in this respect, that they sort into a continuous series of discrete, density specific laminae.   13 shows a repeat of the run shown in Fig. 11 with the change in frictional coefficients. Particle densities remained the same and the new frictional coefficients are applied uniformly to all particles. This resulted in a new shape for the stacked particles, but the tendency to segregate into laminae remains true regardless of the change in frictional coefficients. Clearly, density difference plays a far greater role than friction in the phenomenon of gravity sorting. Finally, starting from the same initial arrangement of Fig. 10 with the same density combination of fluids as the solid-particle case, the behaviors of low-viscosity and highviscosity heterogeneous fluids are simulated and compared against the solid-particle cases. The properties of the constituent fluid particles are listed in Tab. 4. Fig. 14 shows the results. In the case of fluids, viscosity functions in a way that is comparable to the role that friction plays in solid particles. Even highly viscous fluids like muds do not result in sloped surfaces characterized by an angle of repose, as solid particles do. Instead they segregate via buoyancy forces into respective layers, with flat and horizontal bounding surfaces. Low viscosity liquids reach equilibrium more quickly than highviscosity ones, but both result in horizontal and flat bounding surfaces.  Figure 14: Results of multiple liquid particle drop 6 Concluding remarks A new Moving Particle Semi-implicit method has been developed by the authors to simulate the behavior and interaction of multiple solid particles as a continuum. Compared to fluids, solid particles are affected by frictional forces rather than by viscosity. The new method represents a modification of existing MPS equations for application in solid particle problems and successfully simulates particle behavior in the broken-dam and cargo-dump types of problems. The present MPS results for dumping solid particles were validated against the corresponding DEM results, which were previously verified against experiments. The developed MPS method also reproduces via simulation the peculiar behavior of lamination-by-density-sorting from a heterogeneous mixture, a little-understood and yet experimentally-proven real phenomenon. Sedimentologists involved in numerical modeling can benefit from the new model, which can potentially determine upper limits for rates at which laminated sediments can accumulate, for example during large-scale flood situations. By combining with fluid MPS method, the developed program can also be applied to tsunami generation by land-slides.