DEM GPU studies of industrial scale particle simulations for granular flow civil engineering applications

The use of the Discrete Element Method (DEM) for industrial civil engineering industrial applications is currently limited due to the computational demands when large numbers of particles are considered. The graphics processing unit (GPU) with its highly parallelized hardware architecture shows potential to enable solution of civil engineering problems using discrete granular approaches. We demonstrate in this study the pratical utility of a validated GPU-enabled DEM modeling environment to simulate industrial scale granular problems. As illustration, the flow discharge of storage silos using 8 and 17 million particles is considered. DEM simulations have been performed to investigate the influence of particle size (equivalent size for the 20/40-mesh gravel) and induced shear stress for two hopper shapes. The preliminary results indicate that the shape of the hopper significantly influences the discharge rates for the same material. Specifically, this work shows that GPU-enabled DEM modeling environments can model industrial scale problems on a single portable computer within a day for 30 seconds of process time.


Introduction
Granular materials are one of the most widely used materials in civil engineering.The geometrical properties that include particle shape, particle size distribution, surface roughness (texture) have a significant effect on the packing density of the particles.More importantly these properties significantly affect the macroscopic flow behavior i.e. they influence the mechanical properties that include dissipative and restoring forces between contact pairs.The Discrete Element Method (DEM) proposed by Cundall and Stack [1] for geotechnical applications, is starting to mature to the degree that numerous and diverse industrial civil engineering applications are being considered [11].The limiting factor has been the limited number of particles that can be viably simulated for industrial applications.Even though numerous assumptions such as limiting particle shapes to spheres have been employed, the current processing ability of the CPU has limited the advances for industrially relevant simulations.
Recent developments in simulating Discrete Element Models on Graphical Processor Unit (GPU) hardware bring new opportunities to tackle industrial granular applications.In particular, Govender et al [2,3] have simulated industrial processes that include ball mill simulations and hopper discharge using 16 million particles on a single portable computer allowing new computational performance levels to be achieved.In recent collaborative stude-mail: patrick.pizette@imt-lille-douai.fr ies between Mines-Douai (France), CSIR (South Africa) and University of Pretoria (South Africa) the GPU based DEM simulations have been extensively validated for silo discharge and milling simulations using both spherical and polyhedral shaped particles [4,5,7].These studies clearly demonstrated that GPU based DEM can model complex particle shapes in a reasonable computational time when the number of particles are sufficiently large to be practically valid.These studies in particular proposed a systematic comparison between DEM-GPU simulations and experiments for idealized particles (spherical and convex polyhedra) [5,9].This allows for systematic validation of developed GPU codes as the particle complexity increases.
As an extension to the previous studies, this study considers the discharge of an industrial storage silo in an industrial concrete facility in Paris, France.This paper demonstrates the pratical utility of the GPU using a validated GPU-enabled DEM simulation environment to solve an industrial scale granular problem within a day for which purposes we selected 8 million particles.

DEM GPU Approach: BlazeDEM3D-GPU Framework
The computational time requirement is the primary limitation on the application of DEM to industrial problems in particular in civil engineering [11].The highly parallelized hardware architecture of GPU (Graphic Processor Unit) can offer new opportunities for industrial applications.In this study, we use the BlazeDEM3D-GPU framework developed by Govender et al. [2,3,6].This numerical framework has been used and experimentally validated with success to simulate with sphere or more complex shape as for example regular polyhedron silo discharge [4,5,7] and milling operations [5,6].In this paper a laptop with a GTX 980 M GPU and 8GB of DDR5 RAM is used to conduct DEM simulations within a day for 30 seconds of process time using 8 millions of particles.

Contact Law Model
In this study, a linear spring dashpot model is used.The normal force between two particles in contact is calculated by : where δ , δ and n ij are respectively the interpenetration distance, the normal relative velocity and the normal contact vector between two particles in contact.K n is the contact stiffness and C n is the dissipative term.These two coefficients are determined from the coefficient of restitution (COR) and time duration [6,10].The tangential force is computed by a viscous damping model : where C t is the tangenial viscous coefficient and V T the relative tangential velocity computed at the contact point.This tangential force is limited by the Coulomb criterion.Thus, effective friction coefficients µ pp and µ pw are defined respectively for particle-particle and particle-wall.Note that a generated rolling resistance torque is added to counteract the relative rotational motion between two contacting particles involving a rolling resistance coefficient µ roll .

Industrial problematics
To demonstrate the power of a GPU based simulation, we simulate a real-life concrete hopper using around 17 mil-  In particular, we focused on a DEM simulation relating to the new concrete facility of EQIOM located in Paris intra muros in the ecodistrict of Clichy-Batignolles.As shown in Figure 1(a), this concrete facility is composed of several elevated storage silos used to store sand or coarse aggregates (gravel, crushed stone, etc).Concrete is manufactured by first feeding the constituents into a central mixing drum.The required amount by mass of the constituents are fed by gravity.The hoppers of the silos are composed of two coaxial conical shapes with two angles and a cylindrical extension between the two cones.The cylindrical silos have the same design with a diameter of 8 m and a height of 17 m.This study specifically focusses on silo 2 and 6.The two conical hoppers differ though, as the hopper for silo 6 has a cylindrical extension as depicted in Figure 1(b).

Industrial storage silo modeling
The granular material was modeled as 20 mm or 40 mm diameter spheres for the lower or upper limits of an equivalent size for the 20/40-mesh gravel.As illustrated in Figure 2(a), the initial packing was prepared by two stages: firstly, the particles were generated in square box inside the silo.Secondly the particles fell freely in order to obtain a dense packing of particles.The offset between the box of particles with respect to the silo axis controls the asymmetry of the packing.It was ensured that the kinetic energy of the system was close to zero before extracting the final packing.As shown in the Figure 2(b), due to a slight offset in the initial particle packing position with respect to the silo axis, we obtained a non-symmetrical dense packing.
In this study, two packings A and B were generated with just over 8 million 40 mm diameter particles or just over 17 million 20 mm diameter particles (as tabulated in table 1).The discharge is simulated by opening the hatch located at the bottom of last conical hopper.The number of particles amount passing through the orifice may be quantified in terms of number, mass or volume.The particles are systematically removed from the simulation after they have fallen a specified distance of 10 cm.The simulation model parameters are : COR = 0.45, µ pp = µ pw = 0.3, µ roll = 0.01, density = 2.4g/cm 3 for a time-step of 1.10 −4 s.
The linear and angular magnitude velocity fields after 15s depicted in Figure 3.The maximum velocity and velocity gradients correspond to the zone near the orifice as depicted in Figure 3(a).Above the conical hopper in the silo the flow is slower and more homogeneous than the hopper zone where the flow is dominated by shear behaviour.In particular the cross-sectional linear velocity in the hopper along the radius indicates velocity gradients that result in shear dominated flow between the granular particles.Higher in the hopper form blocks of uniform velocity that indicate domains of uniform flow in the granular material as paths of force chains develop in the material.These uniform domains are also highlighted by the angular velocity field in Figure 3(b)) where domains with higher rotation neighbour domains with lower rotation to develop local shear layers between these domains.

Hopper flow comparisons between silo 2 and 6
The height of the cylindrical extension is the main difference between the two storage silos considered.The resulting flow patterns is important as the silos are used to store several types of material with different particles sizes.Figure 5 shows the influence of this extension on the flow profile in the hopper for the 40 mm diameter spherical particles for the two silos for the first 30 seconds.
In particular, the flow pattern was followed by discriminating the particle color as a function of the position in sucessive layers (each layer is composed of around 500000 particles).The observed flow pattern in silo 6 highlights mass flow that are typical of hoppers with The mass flow rate for the two hoppers and for the two particle sizes are shown in Figure 6.This comparison confirms that the cylinder extension enhance the granular flow with a flow rate for silo number 6.This is in accordance with the flow patterns at 20s and 30s.

Applied mechanical stress on the silo
The understanding of the mechanical loading and the flow behaviour is an important key in the silo designs [12].Practioners usually used empirical and continuum models to design hoppers and silos.Applying DEM at an industrial scale may allow for better designs with more fundamental understanding of the structural and particle loading in the various flow regimes.Towards this aim a mechanical loading study is conducted in which we compute the stress applied by the particles on the silos.
The resulting stress tensor of the particles on the silo boundaries is calculated from Love's formulation [8,11] : where the summation is made on all contacts N BF act along the boundaries.V is the particle volume, F i is the i th component of the total contact force and XC j is the j th component of the contact location.The Von Mises stress is computed from the computed stress tensor and plotted in Figure 7.The results show that the maximal shear zones are located near the junction surfaces of the silo and the hopper.The cylindrical extension reduces the shear stress near the orifice.Thus, the stress field exhibits some localizations of shear stress around the cohesive block zones that it is in agreement with the linear and angular velocity observations of the previous section.

Conclusion
A DEM GPU approach based with the recently developed BlazeDEM3D-GPU framework has been used to study the discharging flow of an industrial storage silos.The DEM simulation abilities are considerably increased by the GPU architecture allowing for typical industrial facility.In particular, the flow rate and loading have been studied for an industrial scale problem.

Figure 1 .
Figure 1.(a) Overview of the concrete central and the schematics of design of two of the storage silos (size in mm)

Figure 2 .
Figure 2. (a) Illustration of filling (color by ID) and (b) silo 2 at the initial packing (Range by ID)

Figure 3 .
Figure 3. Silo 2 after 15s of discharge time for the particle case B depicting a cross-sectional view of (a) the magnitude of linear velocity (b) the magnitude of angular (magnitude)

Figure 4 .Figure 5 .
Figure 4. Cross-sectional views of linear velocity fields between the two shaped hoppers for silos 2 and 6 for case B

Figure 6 .Figure 7 .
Figure 6.Flow rates (%) as a function of time (second) for the silo 2 and for the silo 6 and for the two case of packing

Table 1 .
Summary of DEM packing characteristics