Since the formulation of the discrete element method (DEM) by Cundall and Strack (1979), DEM-based simulations have gained significant popularity in recent years. This is largely due to their ability to provide a direct understanding of the granular material interactions and behaviors at microstructural level (e.g., (Guo and Curtis, 2015; Kieckhefen et al., 2020; Kamrin et al., 2024)).

Moreover, coupling the computational fluid dynamics (CFD) calculations of the fluid phase with the DEM computations of the granular phase has provided an enhanced platform to understand mixed solid-liquid-gas behavior under more realistic conditions such as granular materials immersed in an ambient fluid e.g., (Boyer et al., 2011; Amarsid et al., 2017; Cui et al., 2020). For the collapse of a granular column, for example, its run-out may behave very differently, depending on whether the use of coarse or fine grains are immersed partially in water (Si et al., 2018; He et al., 2022). On the other hand, the coupled CFD-DEM simulation has the adoption of much larger time steps in simulation by incorporating the concept of a “parcel” to represent a cloud of particles with similar properties (Bérard et al., 2020). In this approach, the discrete element method models the granular particle system, while the volume-averaged Navier-Stokes equations solve the fluid flow through a multiphase-Lagrangian framework (Crowe et al., 2012; Vijayan et al., 2021; Li et al., 2020; Cui et al., 2022).

The continuum simulation of chute flow is based on solving a depth-averaged granular flow model of the shallow water type (Stoker, 1957; Eglit and Shahinpoor, 1983; Savage and Hutter, 1989). By averaging the properties of the granular flow across the depth of the flow, this model reduces the computational complexity associated with simulating three-dimensional simulations. With the normal velocity component w in the z direction omitted, the mass and momentum equations of dimensionless form can be expressed as follows (Gray et al., 1999; Cui and Gray, 2013; Tregaskis et al., 2022). ∂ h ∂ t + ∇ ⋅ h u ̄ = 0 , (5) ∂ ∂ t h u ̄ + ∇ ⋅ h u ̄ ⊗ u ̄ + ∇ 1 2 g h 2 ⁡ cos ⁡ ζ = S , (6) where h is the flow thickness, u ̄ = ( u ̄ , v ̄ ) represents the velocity and its components in the x and y directions, respectively, and the gradient operator is ∇ = ( ∂ / ∂ x , ∂ / ∂ y ) . The source term S = ( S x , S y ) on the right-hand side accounts for the effects of gravitational force and frictional resistance exerted on the basal surface, and is given in Equation 7 by S = g h i ⁡ sin ⁡ ζ − μ u ̄ / | u ̄ | cos ⁡ ζ , (7) where μ = tan ⁡ δ is the coefficient of friction with δ being the basal friction angle, | u ̄ | = u ̄ 2 + v ̄ 2 1 2 , and i is the unit vector in the x − direction.

Using the conservative variables h , m = h u ̄ and n = h v ̄ , we can rewrite the system of Equations 5 and 6 into a non-strict hyperbolic form in Equations 8 ∂ U ∂ t + ∂ E ∂ x + ∂ F ∂ y = S , (8) where U = ( h , m , n ) T with the superscript T denoting the transpose to a row vector. The respective fluxes and source term vector are given in Equation 9 below E = m , m 2 / h + 1 2 g h 2 ⁡ cos ⁡ ζ , m n / h T , F = n , m n / h , n 2 / h + 1 2 g h 2 ⁡ cos ⁡ ζ T , S = 0 , S x , S y T . (9)

To align with the dimensional DEM simulations, the following scaling is used: the length scale is set to l = 0.03 meters, velocities are scaled by g l = 0.54   m s − 1 , and time is scaled by l / g = 0.055 seconds. Additionally, the Froude number Fr is defined as the ratio of the flow speed | u ̄ | to the wave speed c = h ⁡ cos ⁡ ζ , such that Fr = | u ̄ | / c .

Without considering the complexity of the flow released from the hopper and gate, we focus our simulation on a region of the flow that is sufficiently developed along the slope of the chute. As shown in Figure 1, a typical chute flow set-up consists of a hopper, a release gate, and a chute with a bed and sidewalls. The chute has a downslope length L (e.g., 1,000 mm), a cross-slope width W (e.g., 300 mm), and an adjustable inclination ζ relative to the horizontal. The sidewalls of the chute are usually made of transparent material and have a height H (e.g., 100 mm).Although our current focus is on numerical simulations, we follow the experimental setup for chute flows, such as described in Cui and Gray (2013); Cui et al. (2022), to define our computational domain. This alignment ensures that the initial flow thickness h 0 and velocity u 0 (e.g., obtained through high speed camera measurements using the particle image velocimetry (PIV) method) are well defined at a specific location in the chute. In Figure 1, this location is set at x = X 0 with the Cartesian coordinate ( x , y , z ) aligned with the release gate. Such initial conditions are crucial for both DEM and continuum simulations. Consequently, the continuum simulation domain can be defined within a basal surface region with a downslope length of X 1 and a cross-slope width of Y 1 , as illustrated in the top-right inset in Figure 1. In our case, we select X 1 = 150 mm and Y 1 = 20 mm, as no obstacles are used to divert the flow.

Following the selection of the simulation domain, the DEM simulation is conducted in a three-dimensional region where x ∈ [ 0,150 ] mm, y ∈ [ 0,20 ] mm, and z ∈ [ 0 , h 0 ] mm within the chute. The initial conditions are set h 0 = 22 mm and u 0 = 0.5492 m/s with an inclination angle ζ = 38 ° . This setup results in an initial mass flow rate of 0.1961 kg/s at the inlet. In the continuum simulation, the same initial conditions of h 0 and u 0 are applied. Key physical and control parameters used for the DEM simulation are summarized in Table 1.

The governing Equation 8 is solved numerically using a finite-difference method (Cui, 2014; Cui, 2021) on a uniform H-grid covering the domain defined by x ∈ [ 0,150 ] mm and y ∈ [ 0,20 ] mm. With a grid interval of 1 mm, the computational grid consists of 151 × 21 points for the simulation. Due to the simplicity of the domain without obstacles, only the initial conditions of h 0 = 22 mm, u ̄ 0 = 0.5492 m/s, and v ̄ 0 = 0 m/s are needed at the upstream boundary. The other boundaries develop naturally as the flow propagates downstream. Since the numerical method is explicit, the time step is automatically controlled according to the CFL (Courant-Friedrichs-Lewy) stability condition (Anderson, 1995).Finally, the coupled CFD-DEM simulations were conducted using the Simcenter Star-CCM + package, and our continuum simulations utilized an in-house code developed for the finite difference method based on an explicit non-oscillatory central (NOC) scheme (Nessyahu and Tadmor, 1990).

Since DEM and continuum simulations are based on fundamentally different modeling frameworks in both spatial and temporal domains, it is particularly effective to first investigate the time-dependent development of the simulations. This approach provides an overall perspective, as well as a qualitative comparison, between the two methods. To achieve this, we simulate spherical particles with a diameter of 1 mm flowing down a 150 mm slope inclined at ζ = 38 ° .

As shown in Figure 2, the DEM solutions on the left-hand side generally agree well with the continuum solutions on the right-hand side, particularly at later stages of development, say, from t = 0.152 seconds onward. Even at earlier times, such as t = 0.070 and 0.1 s, the flow evolution is largely consistent between the DEM and continuum simulations. However, the results at t = 0.035 seconds show an interesting difference, especially in the moving fronts. This discrepancy arises because, in the DEM simulation, a sheet of particles is injected at the inlet surface, and the particles then move down the slope in a solid-like state. In contrast, the continuum simulation produces a fluid-phase solution from the inlet, causing its moving front to become infinitely thin. Overall, due to mass conservation, the flow thickness decreases as its velocity increases while propagating further downslope.

Figure 3 provides a more quantitative comparison for the spatial-temporal averaged downslope velocity <u>. This velocity is calculated during the simulation in a form of ∑ j = 1 n i u i / N , where u i is the downslope velocity of a particle, n i is the temporal number of particles in the simulation, and N represents the total number of particles in the steady-state solution. Normalizing <u> by N allows its DEM solution to be directly comparable with the continuum solution, which is averaged over the entire computational domain–either the number of particles in the DEM steady-state solution (e.g., N dem = 34756 ) or grid points in the continuum simulation (e.g., N con = 151 × 21 = 3171 ). In this figure, both the DEM and continuum solutions converge to a steady state around the same time scale, approximately 0.25 s, though they exhibit notably different build-up histories. The figure also shows the steady-state downslope velocities obtained from the DEM solution at the free surface, mid-depth, and basal surface. The free-surface velocity aligns closely with the continuum solution, consistent with the granular flow being treated as a plug-flow type (Cui and Gray, 2013). On the other hand, the DEM solution at mid-depth closely matches the overall velocity profile, indicating a macroscopic representation of the particle velocities. Given that N dem and N con are completely independent measures of particle/grid details, the result in Figure 3 is particularly promising for assessing the microscopic and macroscopic behaviors.

A key advantage of continuum simulations over DEM simulations is their significantly lower computational cost. In our case, under similar computing conditions, the continuum simulation required only 42 s on a 151 × 21 grid with a time step of 0.00055 s, whereas the DEM simulation took 8,750 s for up to 34,756 particles with an implicit time step of 0.001 s, both simulating up to 1 s of physical time. This demonstrates that the continuum simulation is over 200 times faster than the DEM. For a three-dimensional simulation, the continuum approach could still maintain a computational advantage, being approximately 10 times faster on a 151 × 21 × 23 grid, with 23 grid points extended in the z − direction to cover a range of 22 mm ( = h 0 ) .

Further details of velocities obtained from the DEM steady-state solution are shown in Figure 4. In this figure, the velocities have been averaged in the cross-slope y direction to produce a velocity field, < u ′ > and < w ′ > , in the downslope x and normal z directions, respectively. Since there is no obstacle in the flow field, the cross-slope velocity v remains close to zero. In Figure 4A, the velocity field is overlaid with the magnitude of the velocity, showing a gradual increase from the basal surface at z = 0 to the free surface. Also shown in Figure 4A is the steady-state free-surface height of the flow obtained from the continuum simulation, denoted in green circular symbols. The agreement of the free-surface heights between the DEM and continuum simulations is overly good, given the fact that these results are obtained independently, using only the same initial inlet conditions defined by h 0 and u 0 fogure 4b presents the velocity profiles along the flow thickness direction at various downslope locations x = 30,50 , … , 150 mm, obtained from the DEM simulation. Due to the “no-slip” viscous condition applied at the base of the slope, a significant velocity change is observed near the basal boundary, resembling a “boundary layer”. Beyond this layer, the velocity increases roughly linearly toward the free-surface. The corresponding depth-averaged velocity, u ̄ , obtained from the continuum simulation is also shown in green circular symbols in Figure 4B. Since these velocities are based on the free-surface flow for a basal friction angle δ = 23 ° where μ = tan ⁡ δ (see, (7)), they closely match but are slightly higher than the free-surface solution of DEM. This may be generally expected in continuum simulations, as granular materials are assumed dry, cohesionless, and isotropic. By increasing the value of δ to 30 ° , the continuum’s free-surface velocities align well with the DEM velocity profiles, as illustrated by white-filled circular symbols labelled δ = 30 ° that cross the DEM profiles in Figure 4B (note, the vertical sense of h has been discarded for this δ = 30 ° case). In other words, analyzing velocity profiles of DEM can help better calibrate the basal friction coefficient, μ , in continuum simulations, leading to a more consistent solution for granular flows.

Accurately capturing the velocity profile is crucial for modelling the rheological behaviour of granular flows, particularly in continuum simulations (Pouliquen, 1999; GDR, 2004; Franci and Cremonesi, 2019; Tregaskis et al., 2022), while also serving as a metric for assessing the robustness and effectiveness of DEM simulations. As particle characteristics become more complex–such as in shape, size, moisture content, attrition, breakage, and mixing–modelling inter-particle contact and interaction becomes increasingly challenging (Guo and Curtis, 2015; Zhong et al., 2016; Vakis et al., 2018). However, detailed insights into particle velocity profiles are vital for uncovering critical physical insights into particle dynamics and behaviors. Since spherical particles represent the most typical shape, our results could offer a benchmark solution for further studies, including the application of machine learning.

Using the Random Forest method as an example, we briefly explain how machine learning can be applied to granular flow data analysis (Breiman, 2001; Zaki and Meira, 2014; Clarke et al., 2009) from training, optimization, to prediction.

Let X = { x 1 , x 2 , … , x n } be the input dataset, with each feature vector x i = [ x i 1 , x i 2 , … , x i m ] T in R m . In our case, the feature vector x i can be defined to include particle-related information such as position, velocity, particle ID, forces, stresses, and volume fraction at each output time step from the DEM simulation. For example, if we use a time interval of 0.02 s over a total simulation time of 1 s, we would generate N = 100 feature vectors, forming the input (or training) dataset X , stored in the space R m .

Let Y = { y 1 , y 2 , … , y n } represent the corresponding set of labels or target values, where y i ∈ R for regression tasks. In this context, we test RandomForest, LinearRegression, and RidgeRegression models. The RandomForest model operates by building multiple decision trees, denoted by T j , during training process and then aggregating their results to improve prediction accuracy. This can be expressed in Equation 10 as: y ̂ = F x = 1 B ∑ j = 1 B T j x (10) where B is the total number of trees, T j ( x ) is the prediction from the j -th decision tree, and F ( x ) can be a generic machine learning model.

Consequently, a loss function, e.g., the Mean Squared Error (MSE), as given in Equation 11 as follows L Y , Y ̂ = 1 n ∑ i = 1 n y i − y ̂ i 2 (11)

can be used to measure the discrepancy between the true labels Y and the predicted labels Y ̂ .

Without repeating further details of the training process (Breiman, 2001), the model is evaluated using a test dataset ( X test , Y test ) . Metrics such as the R 2 -score can be used to assess the model’s effectiveness, as demonstrated in Figure 5.

The trained RandomForest model is then applied to predict outputs for new, unseen data points x new , using Equation 12: y ̂ new = 1 B ∑ j = 1 B T j x new . (12)

When training a machine learning model, it is essential to split the dataset ( X , Y ) into a training set ( X train , Y train ) and a test set ( X test , Y test ) . Typically, around 80 % of the data is allocated to training, while the remaining 20 % is reserved for testing. This split allows the model to learn patterns from the training data and then assess its performance on the test data to evaluate how well it has generalized those patterns, which can be measured by the R 2 -score for the loss function (11).

In Figure 5, three machine learning models–RandomForest, LinearRegression, and RidgeRegression–are used to train the granular flow dataset with different numbers of feature vectors at N = 30,50,75,100,150 . Notably, the RandomForest model consistently outperforms both the Linear and Ridge regression models, which rely on fitting a linear relation to the observed model. Additionally, increasing the number of feature vectors N does not appear to significantly improve the model scores.

Figure 6 shows the predicted output, y ̂ new , generated by the Random Forest, Ridge Regression, and Linear Regression models for the granular flow dataset. In this case, the models are used to predict particle positions. The images in the top panel, (a) and (b), are drawn from the test data and represent the true steady-state solution of DEM, as shown in Figure 2 e.g., at t = 0.5 seconds. The corresponding predictions by the Random Forest model are shown in Figures 6C, D, by the Ridge Regression model in (e) and (f), and by the Linear Regression in (g) and (h), respectively. The images on the left (a,c,e,g) and right (b,d,f,h) panels show the same result from different viewing angles.

Even with the RandomForest model, the predicted result shows a significant discrepancy from the true solution, highlighting the challenges of applying machine learning to accurately predict even relatively simple granular flow patterns. One possible reason for this disagreement is that the original training dataset was selected from the initial injection phase of the flow through to its steady state. This may have complicated (deliberately) the training process due to the continuously evolving flow patterns, as seen in Figure 2. Nevertheless, there is ongoing development in machine learning techniques to address such complexities, with methods like physics-informed neural networks (PINNs) (Karniadakis et al., 2021; Wang et al., 2021; Mao et al., 2020) and graph neural networks (GNNs) (Choi and Kumar, 2024) being tested for tackling complex fluid mechanics problems.