Impact on granular bed: validation of discrete element modeling results by means of two-dimensional finite element analysis

The study of processes characterized by impulsive nature (i.e. impacts) can be considered of great interest in both physics and engineering disciplines: in the geotechnical field, for instance, their effect on the interaction between soil and structures need to be investigated. The present work aims at the validation, by means of two-dimensional finite element simulations, of a methodology of force calibration which uses results obtained from three-dimensional discrete element analysis for predicting the stress at the base of a granular bed, retained by a movable wall, arising when the system is hit by a projectile. To approach this problem, the low-velocity impact has been modeled as a punctual impulsive force on a granular packing.


Introduction
When modelling the interaction between granular soil and a retaining wall, under both static and dynamic conditions, Finite Element Method (FEM) [1] analysis is often the most practical way to follow, since alternative methods may require great computational efforts to solve boundary value problems. In the last decades many efforts have been made to investigate the soil response and the soil-structure interaction induced by periodic stresses (i.e. earthquakes) [2,3], using FEM analysis [4,5]. Nevertheless, predicting this interaction becomes more complex when an impulsive force acts on the system, as for the case of impacts. Calculating the effect of impacts on the stress state of a soil layer requires both theoretical and experimental approaches. Recent works have concerned both laboratory and numerical experiments [6][7][8][9][10][11][12][13][14]. Impact tests have been made in a centrifuge by Gibbings and Vasile [15] and low velocity impacts into simulated regolith have been performed under microgravity conditions by Colwell [16]. Ho and Masuya [17] approached the calculation of impact load by rockfall on energy absorbing sand-cushioning layer using Finite Element Method. They used a constitutive model based on the Drucker-Prager yield surface, which is a smooth version of the Mohr-Coulomb yield surface, and selected suitable values for the elastic constants and the Mohr-Coulomb constants of the sand. They also found good agreement between the simulations and experimental results for geocells filled with sand obtained by Lambert et al. [18], for example in terms of maximum values of the impact forces, transmitted forces and impact durations.
Although such dynamic process has been extensively studied and hence widely discussed in the literature, attention must be paid to the influence of the kinematic constraints on the stress distribution resulting from an impact, like for example in the presence of a moving retaining wall. In the last years, with the appearance of the Discrete Element Method (DEM) for simulating the behaviour of granulated materials [19], this estimate has benefited from new computational tools. In this case, the uncertainty of the results of simulations can be related, among other factors, to the interaction law governing the collision between particles [20] and to the theoretical assumptions on the duration of a contact [21]. Discrete Element Method was employed to model impact problems for orthotropic media by Liu et al. [22], who used the Mohr-Coulomb type criterion to judge the failure of concrete and added a contact discrete model in the algorithm, capable of calculating the transient process from continuum to non-continuum. DEM simulations of impacts on granular bed have been carried out by Ciamarra et al. [23] and Nishida and Tanaka [24].
To get more confident predictions using DEM simulations, combined approaches making use of FEM simulations can be followed. Although both calculation methods are characterized by inherent limitations, their combined use can corroborate each other. In the last years, several methods for coupling DEM and FEM have been proposed [25], for example by combining them within the same algorithm of calculation [26]. The two methods can also pertain different scales. For example, following a multi-scale approach, Desrues et al. [27] propose a coupling method, in which the DEM calculation refers to the microscale and is aimed at deducing numerically the stress-strain response of the Representative Elementary Volume (REV), and the FEM is used at the macroscale implementing in every Gauss point of all the elements of the mesh a specific REV, providing the stress evolution in the Gauss point during computation.
Alternatively, a DEM-FEM combined approach could also consist of a FE computation for which the macroscopic parameters of the porous medium must be calculated, based to micro-macro relationships [28], from mechanical parameters (e.g.: normal and tangential stiffness) deduced at a particle level. More recently, Nishiura et al. [29] pursued an approach for investigating the dynamic response of ballasted railway track coupling FEM and the quadruple discrete element method (QDEM, [30]), which is based on a constitutive law of a four-particle element. They used QDEM simulations to model the sleeper and the ballast particles and FEM analysis to model the wheel and the rail: as an input for the QDEM simulation, they used the FEM time-series data of the displacement at the nodes connecting the sleeper and the rail pad.
In general, the impact force assumed as an input in the FE analysis cannot be known a priori, unless it is measured during a real impact event and if we exclude the theoretical way for its calculation, generally based on oversimplified assumptions. Conversely, deducing the impact force value from the effects produced could be a promising approach. In this work, a coupled DEM-FEM approach is pursued for deducing the force related to an impact at the boundary of a layer made of spheres and retained by a rigid wall, which is free to move, and consequently for corroborating the DEM simulation of the impact effects by means of FEM results. In this work, DEM and FEM are not jointed within unique algorithm but they are used at different steps, as explained in the following. Since the focus of this work is mainly on the methodological path proposed to deduce the impact force, it has been chosen to implement in the FEM simulations a simple constitutive model (i.e. elastic-perfectly plastic Mohr-Coulomb model), which cannot take into account stress and strain rate effects. From a kinematic point of view, it has been considered of interest the case in which the wall is stable with respect to the rotational collapse and it can only move horizontally ( Fig. 1) until the backfill thrust, which decreases as the wall moves, becomes equal to the interface friction resistance at the wall base.
The system is then out of mechanical equilibrium, since the motion of the wall (mechanical instability), which we hypothesized to start only after the filling stage (i.e. constrained condition during filling), allows the stress regime within the backfill to change with respect to the initial state. The strategy of calculation can be summarized as follows: a. A FE setup of the same size as that employed for DEM is created, with suitable material and interface properties, in order to satisfy both stress distribution and wall Fig. 1 Sketch of the simulated system displacements, with reference to the stage characterized by static actions (filling of the container and reduction of horizontal stresses), as calculated by DEM. b. A FE dynamic analysis is carried out to back-analyze the impulsive force acting at the top of the layer on the basis of the wall displacements induced by the impact, that is forcing the correspondence between displacements from DEM and displacements from FEM. c. Using the impact force deduced at the previous step and implemented in FEM as impulsive load, the maximum force exerted by the granular material at the base of the container during the impact is calculated and compared with that coming from DEM, in order to assess the performance of the proposed calibration strategy and, therefore, to corroborate the DEM results.
Furthermore, a comparison between the impact force numerically deduced by FEM and that calculated by the impulse-momentum theorem, based on the time of collision implemented in DEM at the particle scale, is made.

Simulation model
Discrete Particle Method (DPM) simulations have been performed using Mercury DPM software [31,32]. To compute the dissipative collisions between grains, many models have been created over time. The linear spring-dashpot model, reported in Eq. (1), refers to the normal damped harmonic oscillator force: where γ n is a damping constant and k n is related to the stiffness of a spring whose elongation is ξ, the deformation of the grain [33]. The coefficient of restitution and the time of collision can be expressed by Eqs. (2) and (3).
To derive the interaction force of elastic spheres, the Hertz theory of elastic contact [34] can be used (4) (1) F n = −k n − ṅ.
(2) e n = exp − n 2m eff t n . (3) where ξ is the deformation of the material, Y is the Young Modulus, ν is the Poisson ratio and R eff = (1/R i + 1/R j ) −1 is the effective radius of the colliding spheres [19]. Moreover, to derive the interaction force of viscoelastic particles, the Hertz law has been generalized by Brilliantov et al. [35] (5).
where A is the dissipative constant and it is a function of the material viscosity [19,35].
In the performed simulations, the linear spring-dashpot model has been used, while the Coulomb's law of friction has been implemented for the sliding case. In particular, the coefficient of restitution and the duration of collision can be directly defined by the user and, thus, the spring constant k n and the damping coefficient γ n of the material can be derived by Eqs. (2) and (3). Once the normal coefficients k n and γ n are known, the spring-dashpot normal force can be calculated by Eq. (1). All particles within the system are subject to gravity, acting downwards in the negative z direction.

Simulated system
The system modelled is a 10 cm × 10 cm × 16 cm container made of four lateral walls and open from the top, filled by a polydisperse sample of spheres, where a lateral wall is allowed to move after filling (Fig. 2).
The total number of spherical particles, whose parameters are summarized in Table 1, is equal to 12,568. The time step implemented to solve the equations of motion is set to be equal to 0.02 × t c . The mass of the movable wall is equal to 0.2606 kg and the coefficient of sliding friction between the wall and the base μ wall is equal to 0.58. At first, all particles are generated randomly without overlap above the base of the container. Consequently, they become subject to the gravitational field (g = 9.81 m/s 2 ) and they fall vertically. During the filling of the container all the lateral walls are fixed and only at the end of this process the movable wall is released. The wall motion is caused by the pressure exerted by the granular material against the wall. Although the earth pressure distribution in presence of retaining structures can be, in general, calculated both theoretically and numerically [36], under these kinematic conditions it can be convenient to use a numerical approach [37]. The wall mass and the friction between the wall and the base act together like a regulator, therefore when the force exerted by the granular material on the wall is balanced by the frictional resistance, the wall stops. In other words, the wall moves until the mechanical equilibrium between the granular material thrust and the friction resistance at the wall-base interface is reached. Once the equilibrium is achieved, the impact process can begin. To reproduce the impact, a spherical projectile is dropped from a fixed height of 10 cm above the top level of the granular material with a given initial velocity v 0 of 10 m/s. The projectile is characterized by a radius r p = 0.5 cm and by a density ρ p = 7850 kg/m 3 . The entire process, i.e. filling stage and impact, is set to simulate 5 s of real time. The stress variation generated, during an impact, at the base of the container, characterized by the above-mentioned mechanical instability (i.e. wall motion) has been investigated. Figure 3 shows DEM results in terms of both maximum force exerted by the material at the base of the granular container and wall horizontal displacements, during both the filling of the container and the impact of the projectile.
During the filling process, the maximum force exerted by the granular material on the base of the container starts to increase up to a peak of about 12 N, due to the fall of the particles (Fig. 3a). As previously mentioned, at this first stage the walls are fixed and the system reaches a first condition of equilibrium. After 0.6 s the lateral wall is released and it starts to move due to the force exerted by the material on it (Fig. 3b). When the system reaches the equilibrium the wall stops and, right after, the impact process begins. When the projectile hits the granular bed, a sharp increase of the force experienced at the base can be detected (Fig. 3a). During the penetration of the projectile into the granular bed, the force starts to decrease almost up to zero. It happens because the particles at the base of the container move upwards and loose the contact with the container base instantaneously, due to their reaction to the momentum transferred by the impact process. The gravity pulls down again the particles and it can be detected in the last plateau of the function (Fig. 3a). In parallel, the lateral wall starts to move again due to the perturbation of the system, from its equilibrium state, induced by the impact process. The wall moves until it reaches the new condition of equilibrium (Fig. 3b).

Calibration of the FE model under static conditions
FE simulations of the soil-wall interaction are performed using the finite element software PLAXIS 2D [38], assuming that the calculation is not greatly affected by 3D effects: the quantification of these effects was in any case out of the scope of this investigation. The numerical model consists of a retaining wall founded on a rigid soil foundation and supporting the dry backfill material. The 2D FE domain is discretized by 2103 15-noded triangular elements under plane strain conditions. The plane-strain boundary condition corresponds to a negligible displacement of the system in the z-direction, therefore, along the longitudinal axis, the strain is assumed to be zero ε z = 0. Although PLAXIS 2D allows two-dimensional analysis, stresses are based on a three-dimensional cartesian coordinate system [38]. In these simulations no viscous or absorbing boundary conditions have been used, for consistency with the DEM calculation. The geometrical model is characterized by a solid foundation, supporting the 10 cm high and 1 cm wide wall, retaining the granular material (Fig. 4).
The static stage is supposed to simulate the filling phase of the box, during which the wall can move and, therefore, allows the reduction of the horizontal stresses [39]. The interaction of the granular material with both vertical sides (left side and movable wall) is modelled by means of interface elements (Fig. 4). The interface element allows the separation and the slippage of the soil interacting with the wall in both normal and tangential direction. The shear response of the interface is governed by the elastic-perfectly plastic Mohr-Coulomb failure criterion. Moreover, full fixities are applied at the base of the soil foundation. The wall behavior is modelled as stiff non-porous linear elastic model, assuming Young modulus E′ wall equal to 10 9 kPa and Poisson's ratio ν = 0.25. The wall stiffness has been set very high to avoid wall straining. To reproduce the same friction resistance between the wall and the foundation modelled in DEM simulations, a friction angle δ wall-foundation of 3.21° and a wall mass equal to 2.7 kg have been imposed. The soil foundation stiffness E′ foundation is equal to 10 9 kPa, according to the assumption of rigid base of the box. The backfill material is modelled as elastic-perfectly plastic Mohr-Coulomb model, whose parameters (Table 2) are back-analyzed based on the horizontal displacement of the wall during the static stage computed in DEM (with tolerance lower than 10%), except on the unit weight γ backfill . It has been evaluated as a function of the geometry of the box containing 12568 spheres, characterized by density ρ equal to 1.13 g/cm 3 . Although the DEM model is a polydisperse system in which grains have radii r = 0.25 ± 0.05 cm, the weight of the continuum medium γ backfill has been calculated considering a constant radius of 0.25 cm and γ backfill = 5.7 kN/m 3 . Consistently, the initial porosity of FEM model (n = 0.49) is slightly higher than the initial porosity of the DEM (n = 0.45). The calculation has been averaged on the entire 3D volume of the box and implemented in the FE calculation, assuming the 2D section as representative of the granular system in terms of three-dimensional porosity.
More specifically, the shear resistance parameters and the elastic properties are selected in order to allow the generation of an active mechanism behind the moving wall, without driving to the collapse of the system, as it can occur for badly designed retaining walls (e.g.: sliding and rotation of the wall). Thus, the identified friction angle is 27°, which is larger than that corresponding to the friction coefficient between spheres assumed in DEM simulations   (Table 1), consistently to relationships between macroscopic and microscopic friction angles [28]. The backfill stiffness is determined after parametric analyses, providing a Young modulus E backfill of 1.8 × 10 3 kPa. The Poisson's ratio ν is assumed equal to 0.353, according to the elastic relation (6) given by Tschebotarioff [40]: where K 0 is the coefficient of earth pressure at rest, which refers to the condition where no lateral deformation within the soil mass can occur and is the ratio of the horizontal (σ h ) to vertical (σ v ) effective stress [41,42]. In granular material the magnitude of K o depends on the amount of frictional resistance mobilized at contact points between particles [43]. When a granular soil is loaded for the first time, the frictional forces at the contacts act in such a direction that σ h is less than σ v (K o < 1). Jaky, on experimental basis, suggests that in granular soils the coefficient of lateral earth pressure at rest can be estimated by Eq. (7) [44].
where ′ is the angle of internal friction.
The wall-backfill interfaces have been assumed rigid, i.e. a friction angle δ backfill equal to the backfill friction angle ′ backfill is considered (strength reduction factor R inter = tanδ/ tanϕ′ = 1, where δ backfill is assumed equal to ′ backfill and represents the wall-soil friction angle), according to DEM simulations. This assumption is justified, provided the wall is much harder than the retained granular material.
Numerical results obtained after the static phase are shown in Fig. 5 in terms of the deformed mesh at the end of the filling of the box (Fig. 5a). The filling process gives rise to the pressure behind the retaining wall: the resulting horizontal displacement, the only movement allowed to the wall, mobilizes the earth active pressure [38]. From the distribution of the plastic points (Fig. 5b) a slip zone can be recognized, indicating the onset of the backfill collapse, i.e. the active thrust wedge behind the retaining wall predicted by the Coulomb's theory [45]. This global mechanism can be invoked even if not all the stress points in soil elements, inside the thrust wedge, are in plastic state and the distribution of plastic points is not uniform, as can be observed in Fig. 5b.
In Fig. 6 the distribution of the vertical stress across the width of the backfill at the base is illustrated. The static force F base at the base of the system, due to the backfill weight, has been determined by integrating the vertical stress distribution across the horizontal section passing through the base of the backfill.
It is worth highlighting that FEM simulation provides results in agreement with DEM approach both in terms of force transferred at the base and in terms of horizontal displacement of the gravity wall. In fact, it might be observed a slight overestimation of the static force F base , equal to 7.28 N, which is higher than the one determined with DEM approach, i.e. 4.95 N, but still capable to catch the arching effect [46,47]. This discrepancy of about 2 N, however small  Fig. 6 Vertical stress across the width of the backfill at the base with respect to the peak assessed by DEM during the impact, can be explained considering the differences in the computational methods used. For instance, the FE modeling is not able to describe the inter-particle nature of the material and, therefore, the force-chain network characterizing the system.

FE back-analysis of the impact force
The dynamic phase is set to simulate the impact of a projectile on the backfill surface, still allowing the movement of the wall and monitoring the stress distribution at the base of the container filled by granular material. During the dynamic stage, no numerical damping is considered, thus Newmark parameters α N and b N are set equal to 0.25 and 0.5, respectively. When dealing with finite element modelling of impulsive phenomena, e.g. blasting or rapid impact compaction in soil, some issues arise since they can be governed by frequencies higher than those usually experienced during earthquakes. First of all, it is important to ensure that the waves propagating through the backfill are not filtered at high frequencies. To this purpose, Kuhlemeyer and Lysmer [48] suggested a maximum size of the element equal to oneeighth of the wavelength associated with the highest frequency of the input signal, herein assumed equal to 100 Hz [49]. Furthermore, the time step should be smaller than the ratio of the distance between two nodes and the wave propagation velocity of the soil [38]. In the present study, the average element size is about 3 × 10 −3 m, and it satisfies the above-mentioned criterions. Under dynamic conditions, the dissipative capacity of the soil should be taken into account. Wide geotechnical earthquake engineering literature, i.e. [2,50,51], highlighted the possibility to describe in numerical simulations the cyclic soil behaviour through linear viscoelastic model. In the current work we used this model due to its successful application to related problems, evidenced in, for example, [48,[52][53][54][55]. To consider the viscous damping, accordingly with the Eq. 8, Rayleigh parameters (α, β) have been implemented.

where [C] represents the damping matrix, [M] is the mass matrix and [K] is the stiffness matrix.
The coefficients α and β can be obtained through Eq. 9 [56]: where ω m and ω n are the angular frequencies related to the frequency interval f m -f n over which the viscous damping is equal to or lower than D [57,58]. Given the elastic properties of the modeled granular material, only a small amount of damping ratio has been considered, i.e. D = 10 −4 %. In  this case, the viscous damping is implemented according to the simplified Rayleigh formulation (Eq. 10), which considers the damping matrix as a function of the stiffness matrix only [59].
The Rayleigh coefficient β can be found accordingly to Eq. 11: where D is the target damping ratio (D = 10 −4 %) and ω* is the control angular frequency.
The control frequency f* is assumed equal to 100 Hz and it follows that the Rayleigh coefficient β is 3.18 × 10 −9 . The incidence of this assumption on the output of the simulations has been investigated. In particular, the control frequency of 10 Hz and the frequency range of 100-300 Hz have been imposed through the complete and simplified Rayleigh formulation, respectively. The discrepancy in terms of output calculation when different frequency values are implemented can be considered negligible with respect to the maximum force at the base of the granular bed. Rate effects have not been investigated, for the reasons already discussed in the Introduction.
The impact of the projectile is modelled as a punctual force variable with time according to a triangular function, applied to the middle of the backfill surface (Fig. 7a).
The assumed duration of the dynamic action has been estimated, looking at the DEM results reported in Fig. 3a, to be equal to 0.5 ms. More specifically, it has been put equal to the time taken to increase the initial value of the force exerted by the granular material at the base of the container (4.95 N) to the peak (52 N) and to decrease it again (Fig. 3).
In order to determine the peak value of the equivalent impact force, a FE back-analysis has been performed imposing the horizontal displacement of the wall (Fig. 7b) to be equal to that predicted by DEM (4.2 mm). The FE simulation has been considered to satisfy the required prediction when, at the end of the dynamic process, the wall cumulated an additional displacement, with respect to the static stage, of 4.17 mm.
Since it is not expected that a numerical model can simultaneously satisfy displacements field and stresses field predictions with the same accuracy, three stress points (one in proximity of the left side, one in the middle and one in proximity of the gravity wall) have been selected, with the aim of monitoring the stress distribution at the base of the system. FEM numerical results of the dynamic stage are illustrated in terms of variation with time of the vertical stress for the three points (Fig. 8).

Calculation of the force transmitted at the base and discussion
The force F base transmitted at the base of the backfill has been determined by integrating over the horizontal cross section of the backfill (L box = 10 cm x B box = 16 cm) the vertical stress time history averaged for the three above-mentioned stress points (Fig. 8). This is justified by the evidence that the vertical stress distribution at the base of the backfill can be considered quite uniform, due to the similarity of the curves in Fig. 9. The initial value of F base , as shown in the plot, is that arising from the static stage. After about 2.5 ms the wave propagation process starts to affect the vertical stress, with a sharp increase of the same, while the wall starts moving. In terms of prediction of the peak value of the force transmitted at base, as shown in Table 3, the discrepancy between FEM and DEM is about 5% (Fig. 9). Ho and Masuya [17], in their study already discussed in the Introduction, found a discrepancy between the numerical predictions of the maximum force transmitted through the sand cell (and the underlying concrete basement) and the experimental data from Lambert et al. [18] of 10-11%.
On the basis of the discrepancy found in the present work in terms of prediction of the force at the base, which can be considered to be quite low given the different computational approaches, it seems that the performance of the proposed DEM-FEM calibration strategy is very good. It should provide reliable prediction of the stress transmitted to the roof of a buried structure, covered by an overlying soil layer, when a rigid body hits its surface (i.e. effects of micrometeoroid impacts on a moonbase covered by a regolith layer).
Finally, the computational found value of the force of impact, arising from the back-analysis of the wall displacements, has been also compared to that derived by applying the impulse-momentum theorem, which gives a maximum force of impact equal to about 755 N, if the time interval is taken equal to the time of collision adopted in the DEM simulation (Table 1). This value differs considerably (32%) from the force of impact required to generate same wall displacements and similar force at the base of the container in both FE and DEM simulations, probably due to the simplicity of assumptions on which the theoretical calculation is based.

Conclusions
A DEM-FEM combined approach to simulate an impact process and its effect on a granular bed retained by wall has been proposed. First, the FE model has been calibrated by imposing a wall horizontal displacement in static conditions, as resulting from the thrusting granular material behind the wall after filling, almost equal to that calculated by DEM. Then, using FE calculation, a back-analysis of the impact force has been carried out by imposing the amount of wall displacement induced by the impact to be equal to that obtained by DEM. Finally, the so calculated impact force has been employed as an input for the FE dynamic analysis. The maximum force, calculated by FE analysis, at the base of the granular bed during impact is very similar to that calculated by DEM, with discrepancy of only about 5% in terms of force variation at the base. This discrepancy is even smaller than the one detected between calculated and measured values by other authors when impact tests are performed into granular media [17]. Therefore, the performance of the proposed novel approach with assessing the stress variation at the base of a granular bed, caused by the application of an impulsive load on its surface, seems to be very good.
It is worth emphasising that, although more complex constitutive models for the material could be implemented in the FE calculation, this study achieved its main goal, consisting in pointing out a new method to deduce the impact force generated by an impulsive load from the effects produced by it. When displacement-based, this calibration strategy seems to offer reliable prediction of the force values too. Further numerical investigation could provide information on the importance of 3D effects for this prediction. Furthermore, a full parametric study, implementing different impact energies by changing velocity and mass of the projectile in DEM simulations, could be developed in future work in order to point out the incidence of different kinetic scenarios on the accuracy of force prediction by FE analysis.
Besides the subjects studied in the present paper, there are many more factors which affect the stress at the base of a granular bed, such as, e.g. the specific particle-particle contact force which may include cohesive interaction, brittle, and other. More advanced studies should also consider these factors to obtain more realistic descriptions.
Funding This study was not funded by grants.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.