Dry granular masses impacting on rigid obstacles: numerical analysis and theoretical modelling

The assessment of the time evolution of the impact force exerted by dry flowing masses on rigid obstacles is mandatory for the dynamic design of sheltering structures and the evaluation of the vulnerability of existing structures. In this paper, the results of an extensive numerical campaign performed by employing a discrete element method (DEM) code are presented and the role of different geometrical factors (flow length, height and front inclination) and state parameters (porosity and velocity) on the impact force–time evolution is investigated. The impact process is studied to correlate local information with the macroscopic response and a physically based force–time function, generalising the formula already introduced by the authors for the assessment of maximum impact force, in which each parameter is correlated with the previously mentioned factors, is proposed.


Introduction
Flow-like landslides represent a dangerous geological hazard occurring in mountain territories all around the world [32]. To this category, ranging from dry debris flows to rock avalanches, rapid movements of granular materials containing a large amount of gravel and rock fragments belong. Due to the huge volumes and high velocities, dry flows can pose significant hazard to human lives and cause significant damages to buildings and infrastructures [41].
To mitigate the risk associated with this kind of hazardous phenomena, structural protection measures are usually designed in landslide prone areas, to stop, intercept or slow down the flow [31].
In most national standards, [42,65,66] technical guidelines [21,24,25,42] and scientific reports [40], the design of sheltering structures is mainly based on the pseudo-static approach where the maximum value of the impact force is statically applied to the barrier. In this case, the maximum impact force represents the most relevant input variable and is assumed to depend on flow height, front inclination and impact velocity [12]. Its evaluation is derived from empirical formulae inspired to hydrostatic [4,31,34,66], hydrodynamic [3,4,21,31,42,66], hybrid [3,20] or boulder impact theories [21,30,40,42], based on oversimplifying hypotheses about the behaviour of the granular material and on the use of highly dispersed empirical coefficients calibrated on real cases data [11].
As is well known, the use of pseudo-static approaches is too conservative and, nowadays, displacement-based design is preferred. This makes necessary, under the assumption of uncoupling the analysis of the sheltering structure from the simulation of the flowing mass propagation, the definition of the impact force temporal evolution and this is the final goal of this paper.
In the literature, the methods employed to evaluate the impact force and study the dynamic interaction between the granular mass and the obstacle are based on (i) historical cases, (ii) real and small-scale flume tests and (iii) numerical analyses.
(i) Impact forces recorded after real debris flow events are reported, for example, in Zhang [71] and Hu et al. [29], whereas real-scale experiments have been carried out by Bugnion et al. [9]. These tests are characterised by huge costs and efforts in the device installation and operations. Nevertheless, the impact flow geometry and kinematics cannot be precisely defined.
(ii) Small-scale flume tests have been performed all around the world in the last forty years and represent a very diffuse alternative for studying flow-barrier interaction [1, 14, 16, 17, 38, 39, 50-52, 55, 62, 64, 68]. These experimental test results provide information on some factors affecting the temporal evolution of the impact force, such as the mean value of the front velocity and flow height. These studies have shown that the impact is essentially dominated by the flow regime, in particular by the ratio of inertial to gravitational forces, which is conveniently described by Froude number, F r [5,16,31,37].
What is unfortunately absent in small-scale studies is a detailed description, useful for validating numerical simulations, of both flow kinematics (velocity profile), mass density and front geometry.
Moreover, due to the very small soil volumes employed, while the maximum impact force can be easily measured, obtaining realistic values of the impact force-time evolution it is very difficult, since the impacting mass is extremely short and the Froude numbers are above the typical real case values [62].
(iii) In recent years, numerical methods, capable of dealing with large displacements and simulating soilstructure interactions, have been developed. Some methods are based on continuum mechanics, adopting a Eulerian [50] integration strategy, hybrid Eulerian/Lagrangian techniques [15,[45][46][47] or meshless methods [28]. Nevertheless, these methods appear to be not sufficient to accurately reproduce the obstacle-flow interaction, where sudden changes in field variables, large deformation, intensive compaction and high strain rates occur. Continuum methods should require advanced constitutive models accounting for the dependence of the material response on grain packing and strain rate [15,56,57,69].
Discrete methods, (DEM) [18], have been widely used to model and quantify the dynamic interaction between dry masses and obstacles [2,12,13,17,43,44,49,63,67,70]. DEM allows to take natively into account large displacements, energy transformation, local variation of porosity and internal microstructure evolution, as well as to describe the dissipation mechanisms characterising the granular material in both frictional and collisional regimes [59].
To better understand the impact mechanics, very recently, the authors have published different works presenting the results of an extensive numerical campaign performed by using a discrete element 3D model [11-13, 15, 58]. The approach consists in generating the granular mass just in front of the obstacle, and the velocity is assigned as initial condition. The role of granular temperature/velocity fluctuations, although the authors are perfectly aware of its importance [56,59], for the sake of simplicity, is disregarded. Such an approach allows to investigate separately the role of each geometrical factor (length, flow height, front inclination) and state parameters (porosity and velocity) [12,13] and to better interpret the numerical results. On the negative side, the impact initial conditions may seem somehow independent and artificial, since they are not the outcome of a run-out simulation performed on the same material. In Calvetti et al. [12] the results have been interpreted from a macroscopic point of view in terms of maximum impact force and a new design formula has been introduced. Several design factors determining the impact force, namely the geometry of the sliding mass (length, width and flow height, inclination of the front) and the impact velocity, have been considered. In Calvetti et al. [11], the previous numerical results, have been compared with the most common approaches and empirical formulas available in the literature, as inferred either from small-scale flume tests or field measurements. Calvetti et al. [13] investigated the impact phenomenon from a micromechanical point of view, by correlating both particles kinematics and contact forces distribution as well as the evolution of the energy contents, with the observed macroscopic response, in particular during the first-time instants of the impact.
By employing the same numerical model, this paper focuses on an unexplored aspect that is the time evolution of the impact force.
The interaction mechanisms governing the phenomenon, from the impact inception until the complete arrest of the mass, are investigated. The DEM model allows (i) to obtain a data set, previously not available in the literature, highlighting the effect of the initial conditions on the shape of the force-time curve and (ii) to provide a more Normal contact stiffness k n =D (MPa) 300 Shear contact stiffness k s =D (MPa) 75 Interparticle friction coefficient l s 0.3 Obstacle friction coefficient l w 0.6 comprehensive framework, useful for critically interpreting experimental data. In this paper, the investigated range of parameters is based on data reported in the literature for real debris flow events [9,19,36,54,60] in terms of average impact velocity, mass porosity, front inclination and flow height. All the considered impacts are characterised by Froude numbers ranging from 0 to 3, consistent with those estimated for most common real-scale observations [16,31,62].
Following the same line proposed in Calvetti et al. [12] for the maximum impact force, a parametric formula describing the time evolution of the impact force as a function of input data is also introduced.
The paper is organised as follows: in Sect. 2, the DEM numerical model, the micromechanical parameters and the imposed initial conditions are detailed. In Sect. 3, the DEM numerical results are analysed and critically discussed to highlight the physical and mechanical processes responsible for the macromechanical response. In Sect. 4 a comprehensive formula for simulating the force-time evolution, based on the previously described mechanical interpretation of the impact process, is proposed.    [38], is analogous to that described in [11-13, 15, 58].
The impacting dry soil mass is simulated as an assembly of polydispersed spherical particles of density q s with a uniform diameter distribution ranging between D max and D min and average grain diameter D. The boundary and the obstacle are described as rigid wall elements.
Grain-grain and grain-wall local interactions are governed by normal and tangential contact stiffnesses (k n and k s ) and contact friction (l s ). As in Calvetti et al. [12,13] and Ceccato et al. [15], the mechanical response of a real granular material, with rough and angular grains, is reproduced by inhibiting particle rotations [10,67]. The authors have previously performed some analyses not inhibiting grain rotations [15] and did not observed significant quantitative differences in the dynamic contribution of the impact force. In fact, since the material is loose and flowing under very large velocities, shear localisation phenomena do not occur. Similar results could be obtained using irregular clumps instead of spheres, as was shown in Albaba et al. [2].
The choice of micromechanical parameters (Table 1) is based on the calibration procedure reported in Calvetti [10], where the ranges of model parameters to be used for simulating the response of a sand are provided, based on the reproduction of triaxial compression experiments. The authors have performed triaxial tests by employing the parameters reported in Table 1 (not reported for the sake of brevity) and obtained a response analogous to that obtained by Calvetti [10]. In particular, the residual friction angle is equal to / 0 = 36°and the critical state porosity n cr = 0.416.
No contact damping is used. The model geometry and the initial conditions are illustrated in Fig. 1.
Simulations are performed under plane-strain conditions: the dry granular mass is confined between two lateral smooth and rigid wall elements. Basing on preliminary sensitivity analyses, in order to prevent boundary effects, [12,13,15], the channel width b is imposed equal to 8 D.
The base wall is characterised by a friction coefficient l b . The obstacle is modelled as a vertical rigid wall of height H and friction coefficient l w . As shown in Ceccato et al. [15], the value of wall friction does not affect the impact force value.  The geometry of the impacting mass is characterised by a length L, height h 0 , front inclination a 0 , initial velocity u 0 and porosity n 0 .
The results of the sensitivity analyses aimed at investigating the influence of D/h 0 and D max =D min on numerical results are reported in Appendix D. The authors have shown that, in agreement also with the results of Ceccato et al. [15] and Cui et al. [17], by imposing D/h 0 B 0.067, the number of particles is large enough to obtain both sufficiently accurate description of the problem and a reasonable computational time (Fig. 27). Moreover, the effect of D max =D min on the temporal evolution of impact force is negligible (Fig. 28).
The influence of the input parameters was analysed by performing several tests for different combinations of initial conditions ( Table 2). The range of investigated impact conditions is consistent with data concerning real-scale phenomena [9,16,19,36,53,54,60]. Among the 80 numerical tests performed, only 19 tests are described in detail in the following sections. The relative geometrical parameters and initial conditions are reported in Table 3.
The same simulation procedure employed by Calvetti et al. [12,13] has been adopted: (i) initially, the dry granular mass is generated under zero gravity conditions. In this stage, particles are frictionless and the mass is free from contact forces; (ii) interparticle friction l s , gravity and initial impact horizontal velocity u 0 are assigned to all particles; (iii) the impact is simulated and the evolution of the horizontal component of the impact force is monitored over time.

DEM results
In this section, the DEM numerical results are discussed to correlate the phenomena occurring within the granular mass during the impact process with the time evolution of the impact force.
In Sect. 3.1 a reference impact is studied, whereas in Sect. 3.2, the role of each initial condition and geometrical parameter in affecting the dynamic phenomenon is analysed.

The reference impact
In this paragraph, the reference impact test (Test 1 of Table 3), characterised by a 0 = 60°, n 0 = 0.45, L = 60 m and u 0 = 8 m/s, is considered.
The force-time curve plotted in Fig. 2 is characterised by an initial pseudo-linear increase until a first peak force, f p , is attained at t = t p . Then, the response presents a quite irregular trend, characterised by numerous sub-peaks with oscillations of rather constant period and progressively decreasing amplitude. When oscillations nullify, the curve converges to a final constant residual value f res . As is shown here in the following, the shape of the force-time curve is mainly affected by three distinct phenomena: I. the evolution of the contact height; II. the compression and rarefaction waves propagation; III. the progressive arrest of the mass in front of the obstacle.

Evolution of contact height
As is evident in Fig. 3, the contact height h between the mass and the obstacle evolves with time t during the impact process. h initially increases about linearly, while the geometry of the front progressively flattens and part of the mass piles up behind the obstacle forming a triangular dead zone already observed by Calvetti et al. [13]. At t ¼ t (0.21 s in this case) h reaches the flow height h 0 . The linear increase in h for t\t is related with the planar shape of the flow front ( Fig. 4): from which: For t [ t, the incoming material flows along the surface of the dead zone, particles are deflected upwards [13] and h continues to increase, although more slowly, from h 0 to the final value h res (4.95 m). This process was also observed by Jiang et al. [39], and Ng et al. [51] in experimental flume

Compression and rarefaction waves propagation
As was observed by Calvetti et al. [13], Ceccato et al. [15] and Federico and Amoruso [22] using DEM, material point method (MPM) and finite element method (FEM) respectively, during the impact process when the mass hits the obstacle, a compression wave, generated at the mass front, travels backwards within the granular material with velocity u m , creating highly loaded horizontal force chains. At the first peak time t p , the dynamically created compressive force chains become unstable, this resulting in force chains disruption. This process, which is analogous to the typical buckling phenomenon [13,15], produces the onset of a sort of rarefaction wave propagating from the surface toward the bed with velocity u f , and the impact force drops (Fig. 2).
According to FEM and MPM analyses, based on a continuum approach, the rarefaction wave velocity depends uniquely on density and coincides with u m , while DEM simulations show that u r decreases with both Froude number (Fig. 14) and porosity [58].
To better interpret the trend characterising the forcetime curve, in Fig. 5, four numerical test results are compared (a 0 = 90°and 60°; confined and unconfined flow). The confined flow represents an ideal condition which is obtained by adding a horizontal wall at the top of the mass. For both front inclinations, the initial behaviour (i.e., t\t p ) is unaffected by the confinement condition.
As observed by Calvetti et al. [13] in case of vertical front, if the flow is confined, the impact process resembles that of a solid body impacting on a rigid wall and the force value, that is almost constant, is approximately equal to that calculated by employing the linear impulse theory f ¼ q Á u 0 Á u m Á h 0 . The compression wave velocity u m can be easily numerically assessed under confined conditions as described in Calvetti et al. [12,13] and Redaelli et al. [58] and, in case of n 0 = 0,45, u m = 160 m/s. In contrast, under unconfined conditions, the impact force drops at t p due to the buckling process.
In case of inclined fronts, if the flow is confined, the force continues to increase about linearly until t ¼ t (i.e. when the contact height attains h 0 ) getting a force value practically coincident with that obtained for a 0 = 90°. In the unconfined case (reference impact), the impact force starts decreasing when t [ t p , which is significantly smaller than t. In other terms, as is shown in Fig. 3 and in agreement with Calvetti et al. [12] and Ceccato et al. [15], the peak occurs when the contact height is smaller than the flow height (h t p À Á h 0 Þ. After the buckling phase contact forces can regenerate due to the pressure exerted by the material flowing behind the dead zone. The successive oscillations characterising the force-time curve (Fig. 2) are thus associated with a continuous regeneration (compression wave) and disruption (due to total or partial buckling) of the force chains network in the granular mass and are signature of the hybrid solid/fluid nature of the granular material [13,51,59].

Progressive arrest of the mass in front of the obstacle
During the impact process, the dead zone progressively extends backwards. As was observed by Shen et al. [63], Luo and Zhang [48], Xiao et al. [70], the increase in the dead zone thickness justifies the attenuation of subsequent impacts and the reduction of the force oscillation amplitude. When the final geometry of the deposit is attained, the residual force f res is obtained (425 kN/m for Test 1) (Fig. 2).
The trust coefficient k Ã res is defined as the ratio of f res to the hydrostatic pressure being q ¼ q s 1 À n 0 ð Þ the mass density and g the gravity acceleration. For Test 1, k Ã res ¼ 2:47, which is between the value of the active and passive trust coefficients (k a = 0.26 and k p = 3.85), calculated by assuming a critical state friction angle of 36°for the granular mass (Sect. 2).

Parametric analyses
In Fig. 6, the effect of initial velocity u 0 , initial porosity n 0 , flow height h 0 , front inclination a 0 and flow length L on the impact force evolution is summarised.
The role of impact velocity is particularly relevant both quantitatively and qualitatively (Fig. 6a): for decreasing u 0 , the maximum impact force strongly reduces, the peak time instant increases and both impact duration and residual force decrease. For low values of u 0 , the oscillations progressively disappear and the difference between the maximum impact force and the residual force becomes negligible.
Maximum impact force, peak time and the residual force slightly decrease with n 0 (Fig. 6b); on the contrary, the amplitude of oscillations increases with n 0 .
In Fig. 6b, both porosity and state parameter W 0 ¼ n 0 À n cr [7] are given, where n cr is the porosity at critical state corresponding to the atmospheric pressure. W 0 is the most important parameter governing the material response [69]. In a flowing material, W 0 is usually positive and tends to increase with flowing velocity [56,59]. For this reason, in the following, the authors decided to impose as reference porosity n 0 = 0.45, which is above the critical value for the considered grain size distribution (n cr = 0.416). Larger values are not considered since the response of the materials would be more similar to that of a granular gas, without contacts among particles [58], conditions typical of granular materials with not realistic flowing velocities. An extensive study on the role of initial porosity is in Redaelli et al. [58].
The flow height has a strong influence on the overall evolution of the impact force (Fig. 6c). In particular, as h 0 increases, maximum impact force, peak time, oscillation amplitude period and residual force markedly increase, whereas the initial slope of the force-time curve does not change.
The front inclination has a large influence on the initial evolution of the impact force (Fig. 6d). In fact, both maximum impact force and initial curve slope increase with a 0 . For large front inclinations, the first peak is attained very rapidly (virtually immediately for a 0 = 90°) and is much higher than the following peaks. However, all curves converge with time to a common trend, characterised by the same oscillation amplitude and residual force.
As is evident in Fig. 6e, both maximum impact force and peak time are not influenced by the flow length. On the contrary, the evolution with time of the impact force is severely influenced by L: (i) in case of sufficiently short Here below, the role of initial conditions on the three processes already discussed in Sect. 3.1 is studied.

Evolution of contact height
In Fig. 7, the evolution with time of the normalised contact height h Ã ¼ h=h 0 is plotted for impacts characterised by different velocities (Fig. 7a), heights (Fig. 7b) and front inclinations (Fig. 7c).
From a qualitative point of view, all graphs share a similar trend. Both initial slope and final value h Ã res ¼ h res =h 0 increase with u 0 (Fig. 7a) and decrease with h 0 (Fig. 7b). For low impact velocities (u 0 = 2 m/s in this case), h Ã res does not reach 1: the mass stops before the complete flattening of the front. The initial slope increases also with front inclination (Fig. 7c). The final value h Ã res is barely affected by a 0 . Note that for a 0 = 90°, h Ã ¼ 1 at t ¼ 0.
The geometry of the deposit (residual state) is plotted in Fig. 8 for different Froude numbers F r0 , being   [5,16,31,37], and two different front inclinations. In Fig. 9, the normalised residual height h Ã res , is plotted as a function of F r0 for different front inclinations. From these results it is possible to conclude that: • h Ã res increases with F r0 ; • when F r0 \ 1, h Ã res depends on a 0 , in particular, for F r0 ! 0 and a 0 = 90°, h Ã res ! 1 whereas when a 0 \ 40°, h Ã res ! 0; • for F r0 sufficiently high (F r0 [ 1Þ, h Ã res barely depends on a 0 .

Compression and rarefaction waves propagation
In this section, the effect of initial conditions on the oscillatory trend of the force-time curve is investigated.
In Fig. 10, the first peak time t p is plotted as a function of h 0 , u 0 and a 0 : (i) t p linearly increases with h 0 (Fig. 10a), (ii) t p decreases with u 0 approaching the limit t p ¼ 0 for sufficiently large u 0 values (Fig. 10b), (iii) t p decreases with a 0 and t p ¼ 0 in case a 0 ¼ 90° (Fig. 10c). This implies that the dependence of t p on h 0 , u 0 and a 0 is analogous to that one given in Eq. 2 for t.
The analysis of the numerical data allows the authors to state that the ratio t p =t is not affected by F r0 but depends on the material porosity/deformability. To emphasise it, in Fig. 11 t p =t is plotted for two different values of u 0 , h 0 , and a 0 = 40-60°and two values of state parameter W 0 , one positive and one negative: although a certain dispersion of data is observed, t p =t tends to increase with W 0 .
The effect of initial conditions on buckling is illustrated in Fig. 12, where the evolution of the normalised impact In case of a 0 ! 80°and larger impact velocities (u 0 ! 8 m/s), after the peak, both the normalised impact force and coordination number decrease very rapidly to zero (Fig. 12a) and remain zero for a long time, signature that the contact network is completely disrupted [59].
As was observed by Calvetti et al. [13], in this phase rarefaction has involved the entire thickness of the stratum. The material has virtually liquefied: vertical velocities of particles are predominant with respect to horizontal and a jet-like phenomenon tends to occur.
The buckling phenomenon tends to reduce for decreasing impact velocities and front inclinations (Fig. 12c, d). For u 0 \ 8 m/s and a 0 60°neither the impact force nor the coordination number nullify signature that the disruption of force chains is partial.
Buckling tends to reduce in the secondary oscillations, as is evident in Fig. 13, where the trend of the coordination number is plotted for impacts characterised by different u 0 values. Both Z, at t ¼ t p , and the amplitude of its oscillations increases with u 0 . If u 0 is sufficiently low (u 0 = 2 m/ s), the fluctuations are negligible.
In all impacts, after the first peak, impact force oscillations reach a regime condition characterised by a period T (Fig. 6), related to the time employed to generate (compression wave) and erase (rarefaction wave) the force chains network. As was already observed, T increases with both u 0 (Fig. 6a) and h 0 (Fig. 6c). The counter-intuitive influence of impact velocity is justified by the disruptive effect of buckling that increases with u 0 . To quantitatively interpret this evidence, the quantity u r ¼ h 0 = T 2 À Á is introduced. u r can be interpreted as the velocity of the rarefaction wave, travelling the distance h 0 to fully nullify the contact forces, which occurs in half of the oscillation period. As is shown in Fig. 14, u r is a decreasing function of F r0 and is not affected by front inclination.

Progressive arrest of the mass in front of the obstacle
Due to the progressive dissipation of kinetic energy, the impact force oscillations gradually disappear and the impact force gets its residual value f res . As was already discussed, f res mainly depends on u 0 , h 0 and L.
Several literature contributions discuss whether the residual earth pressure is close to either active or passive conditions [48]. Some authors suggest that the residual force is equal to passive thrust, since the material is in a compressive state [19,39,50]. However, other authors observe that the residual earth pressure ranges between the active and passive one [51,52].
To discuss this issue, in Fig. 15a, the thrust coefficient k Ã res (Eq. 4) is plotted as a function of F r0 . The dashed lines in Fig. 15a represent the passive and active thrust coefficients, k p and k a , respectively. It is evident that: • for nullifying values of F r0 (quasi-static flow), k Ã res tends to k a ; • for 0:4\F r0 \1, k Ã res remains approximately constant and close to the passive thrust coefficient; • for F r0 [ 1, k Ã res decreases and approaches the unit value for sufficiently large F r0 values. In this limit case, the response of the material is more similar to that of a fluid at rest in which the state of stress is approximately isotropic [59].
The influence of F r0 and ratio L=h 0 , on thrust coefficient k Ã res is illustrated if Fig. 15b

Impact loading function
In this section, the authors provide an impact force-time function to be employed for design purposes and they discuss the validity domain of the formula. Since the obstacle is assumed to be rigid, the formula provides an upper bound for impact force.

Analytical expression
The design function suggested here below stems from the physical and mechanical interpretation of the numerical DEM results of Sect. 3 and represents an improvement of the formula previously defined by the authors [12] for the maximum impact force. For the sake of generality, the impact force function is expressed in non-dimensional terms by introducing the following quantities: being T ¼ 2h 0 =u r the previously defined period (Fig. 14). u r can be calculated by using the interpolation function u r ¼ u r F r0 ð Þ, given in Appendix B (Eq. 20). The analysis of the numerical simulations suggests f Ã to be expressed as a linear combination of an inertial force f Ã i , associated with the integral of a uniform horizontal pressure distribution, and a residual force f Ã res , corresponding to the integral of a linear pressure distribution: where n t Ã ð Þ is a decay function describing the decay of kinetic energy: n t Ã ¼ 0 ð Þ¼1 and nðt Ã ! 1) = 0. The role of n t Ã ð Þ is to progressively reduce the inertial force contribution in favour of the residual one.

I. Residual force contribution
The residual force, that is the integral of the linear pressure distribution exerted by the mass at rest in front of the obstacle, can be written as: The values of h Ã res and k Ã res are calculated by using the interpolation functions h Ã res ¼ h Ã res F r0 ; a 0 ð Þ and k Ã res ¼ k Ã res F r0 ð Þ; given in Appendix A (Eqs. 13, 15, respectively), derived from Figs. 9 and 15a, respectively. The reliability of Eq. 8 is confirmed by Fig. 16 where DEM data are compared with interpolating results in terms of residual force.

II. Inertial force contribution
The simultaneous propagation of compressive and rarefaction waves during the impact process suggests the use of the following expression for the inertial force contribution: where.
• F r;u m is the intrinsic Froude number [12], defined as: being u m the compression wave velocity within the frictional mass (Sect. 3.1.2) [13]. • t Ã p is the non-dimensional peak time, introducing a phase shift in the oscillatory response (Sect. 3.1): where t is calculated according to Eq. 2 and c n is calculated by using the interpolation functions c n ¼ c n W 0 ð Þ given in Appendix B, deriving from Fig. 11  (Eq. 16). For a 0 = 90°, both t and t Ã p tend to zero. The accuracy of the analytical expressions for the peak time (Eq. 11) is testified by Fig. 17 where the dimensional numerical and calculated data are compared.
• h Ã t Ã ð Þ provides the evolution of the contact height with time, according to the trends reported in Fig. 7. Its expression, depending on both h Ã res and t Ã , is reported in Appendix B (Eq. 17).
The first term in Eq. 9 refers to the deformable solid body and derives from the linear impulse theory, whereas the second term is associated with the propagation of the rarefaction wave and defines the oscillating trend over time.
As F r0 decreases, the oscillation period tends to zero (Fig. 14) and the second term in Eq. 9 vanishes. In this limit case, Eq. 9 reduces to the formula of the elastic body model (Sect. 3.1.2) [10].

III. Decay function
The following expression has been assumed for n t Ã ð Þ : where c 1 and c 2 are non-dimensional coefficients which have been determined by interpolating the DEM impact force-time curve with Eq. 7 (Appendix C). c 1 (Eq. 21) regulates the shape of the decay curve, decreasing with both F r0 and a 0 , whereas c 2 (Eq. 23) governs the speed of decay and increases with both F r0 and a 0 .
In Fig. 18, the DEM force-time curve relative to the reference impact (Test 1) is compared with Eq. 7.
The dotted line, coinciding with the upper envelope of the force-time trend, corresponds to Eq. 7 without oscillations, whereas the dashed line, coinciding with the lower envelope, corresponds to the quasi-static contribution of Eq. 7. The good agreement between the numerical DEM data and the proposed impact loading function is evident. Differently to DEM simulations, the force computed by using the formula presents a second peak that is larger than the first peak. To explain this discrepancy, it is worth observing that, in DEM simulations, when the rarefaction wave propagates downward, the contact force chains network is completely disrupted and the re-solidification time 2). In case of large initial kinetic energies, this effect is so destructive that the material struggles in dissipating the kinetic energy and regenerating the force chains network. This phenomenon, that is even amplified by the absence of contact viscous damping, is not taken into account by the analytical formula. In Fig. 19 the DEM data discussed in Sect. 3.2 are compared with the impact loading function (Eq. 7) for different values of flow height (Fig. 19a, b) front inclination (Fig. 19c, d) and initial velocity (Fig. 19d, e). It is evident that the proposed formula, is capable of capturing both quantitatively and qualitatively the force-time evolutions, and of taking into account the dependence of the impact force on impact velocity, flow height and front inclination. The previously mentioned discrepancy tends to disappear when the impacting mass is characterised by a sufficiently small kinetic energy (small initial velocity or flow thickness) and the accordance between DEM results and analytical formula is excellent.

Validity domain
The impact loading function has been inferred from DEM numerical simulations and, for this reason, can be considered applicable within the validity domain defined by the ranges of initial conditions considered by the authors. The assumptions and the choices of input parameters are summarised here below.
This constraint is reasonable for decelerating flowing masses as it happens when the mass is approaching the sheltering structure placed downstream. This also suggests that in small-scale flume tests, where smaller values of front inclinations are observed [1,2,50], the length of the channel is in general too short to achieve steady-state conditions and the mass is still accelerating.
The authors have observed that, under this hypothesis, the evolution of the impact force is not affected by the length of the impacting mass (Fig. 15b). This information is of great importance since, in the literature, this parameter is usually disregarded in the impact force evaluation. In fact, in small-scale flume tests the ratio L h 0 is in general too short [1,2,39,50,62], and, as a consequence, (i) the residual force is underestimated and (ii) force fluctuations are not observed.
• A uniform velocity and porosity profile.
The velocity profile of the mass approaching the obstacle depends on the propagation phase and on boundary conditions. Experimental tests [6,33,39,47,61] have shown that this hypothesis is reasonable in case of masses flowing on a flat-frictional base, where the velocity profile tends to plug in the upper part, whereas a thin zone of higher shear is present close to the stationary surface. In case of extremely rough sliding surfaces, the velocity profile decreases gradually with depth and the hypothesis of uniform velocity profile is not realistic [8,39]. In Appendix D, the authors have studied different impacts by imposing different velocity profiles (Fig. 29a) and observed that maximum impact force [12], first buckling, residual force and the impulse value are negligibly affected by the shape of the velocity profile (Fig. 29b,c).
• Dry masses with D=h 0 0.067. Under this hypothesis, as shown by Cui et al. [17], particle size effects can be disregarded in the impact force assessment (Appendix D). But consequently, the Within the validity domain, the formula can be easily employed by the designer, since all the coefficients are only function of the initial conditions at the pre-impact instant of time (namely initial Froude number and front inclination). The only parameter that is not directly provided is the compression wave velocity. Nevertheless, several empirical relationships are available in the literature to correlate relative density and the velocity of propagation of compression wave in porous media [26,27].

Conclusions
In this paper, a discrete element method numerical code has been employed to study the dynamic interaction between a dry granular mass and a rigid obstacle. The impacting mass has been generated just in front of the obstacle at the pre-impact time instant. The DEM model proposed by the authors has allowed to obtain a data set, previously not available in the literature, highlighting the effect of the initial conditions (flow length, width and height, front inclination, impact velocity) on the temporal evolution of the impact force. Independently of initial conditions, the force-time curve exhibits a trend characterised by an initial increase in the impact force, until a first peak is attained, followed by different sub-peaks occurring with a precise oscillatory period and with a progressively decreasing amplitude, until a final constant residual value is got.
The mechanical and physical interpretation of the numerical DEM results put in evidence that (i) the granular mass involved in the impact process behaves like both a solid and a fluid material and (ii) the macroscopic response is governed by different concurring mechanisms: (a) the progressive increase in the contact height between the mass   The numerical results obtained have allowed also to critically interpret the experimental data available in the literature. In fact, experimental tests are characterised by the lack a precise description of flow conditions, flow kinematics and front geometry. For example, the authors have discovered that the flow length L, that is in general neglected in most of the empirical formulae, is an important parameter influencing the shape of the force-time curve and the residual force value. For this reason, smallscale tests are useful to provide qualitative information about the most important parameters affecting the impact force but fails in providing a reliable input for a quantitative assessment of the impact force.
Starting from the physical-mechanical parametric analysis of the impact process, the authors have proposed an impact loading function that can be employed by designers since each parameter has been correlated with either the imposed initial conditions or the initial Froude number, quantifying the dynamicity of the impact process.
These initial conditions are assumed to be obtained as output of the propagation analysis that can be performed by employing numerical simulations either with depth averaged models or fully 3D analyses.
The formula is capable of capturing both qualitatively and quantitatively the numerical force-time evolution for many real cases situations and can be employed as an input loading function for the dynamic design of sheltering structures.
As a consequence, the results obtained in the paper would have a paramount impact in the design of rigid protection structures, overpassing the limitations of pseudo-static approaches.

Appendix
Appendix A: Residual force contribution In Eq. 8, f Ã res is assumed to be a function of both h Ã res and k Ã res .
Here below, the expression for h Ã res ¼ h Ã res a 0 ; F r0 ð Þ and k Ã res ¼ k Ã res F r0 ð Þ are given. The DEM numerical data of Fig. 9 have been fitted (Fig. 20) by employing the following analytical expression: where parameter a 1 ¼ a 1 a 0 ð Þ is given by with a 0 expressed in degrees. Analogously, in Fig. 21 the comparison between DEM data of Fig. 15a and the following interpolating function is illustrated: Appendix B: Inertial force contribution The expression for c n ¼ c n W 0 ð Þ, employed in Eq. 11 to assess the peak time, has been obtained by fitting the DEM numerical data of Fig. 11 with the following straight line (Fig. 22): The evolution with time of h Ã is described by using the following equation: where hi = 1, when hi ! 0 and hi = 0, when hi\0. Fig. 23), governs the grow rate of h Ã for t Ã [ t and can be calculated by employing the following equation: where a 2 ¼ a 2 a 0 ð Þ is given by: Equation 17 can describe the trends of h Ã t ð Þ already discussed in Sect. 3.2 (Fig. 24).
To calculate u r , necessary for the estimation of period T, the DEM numerical data of Fig. 14 have been fitted by employing the following analytical expression (Fig. 25): Appendix C: Decay function The decay function n t Ã ð Þ (Eq. 12) is assumed to depend on two non-dimensional coefficients: c 1 regulating the shape of the decay curve, decreasing with both F r0 and a 0 , and c 2 , The values obtained for the coefficients c 1 and c 2 are reported in Fig. 26a and b, respectively, together with the proposed fitting equations: with a 3 ¼ a 3 a 0 ð Þ a 3 ¼ À0:0026 Á a 0 þ 0:55 ð22Þ with a 4 ¼ a 4 a 0 ð Þ a 4 ¼ 0:009 Á a 0 À 0:257 ð24Þ Appendix D: Parametric analyses Influence of D=h 0 In Fig. 27, the results of the sensitivity analysis aimed at investigating the influence of D=h 0 on the time evolution of the impact force, is shown. For a 0 = 60° (Fig. 27a) and in case of D=h 0 [ 0.067, the response is affected by evident numerical oscillations. In case of a 0 = 90°, a discrepancy is observed in the initial value of the impact force related to the interaction with the first line of grains in contact with the wall [1]. In both cases, it is possible to observe that, for D=h 0 0.067, the effect of D=h 0 on impact force is negligible in terms of maximum and residual impact forces, oscillation period and impulse. For this reason, in this paper, in order to obtain accurate results in a reasonable computational time, all simulations have been performed by employing an average particle diameter D/h 0 \ 0.067.

Influence of D max =D min
The results of the sensitivity analysis performed by using the same value of D and by changing D max =D min are shown in Fig. 28. It is evident that the effect of D max =D min is negligible in terms of maximum and residual impact forces, oscillation period and impulse, except for the case of the ideal unrealistic granular monodisperse medium D max =D min ¼ 1 ð Þ , where the residual force is slightly larger and the oscillations less marked.

Influence of velocity profile
In this Appendix, the influence of velocity profile on the temporal evolution of the impact force is studied. For the sake of simplicity and under the assumption that the mass is decelerating, flows with vertical fronts are considered. At any rate considering different inclinations would make arbitrary the definition of the velocity profile the front represents a process zone where unsteady conditions occur. Four different velocity profiles, shown in Fig. 29a, characterised by the same mean velocity value, are taken into consideration. Note that, at the inception of the impact, the overall momentum is the same, but kinetic energy is different for the three profiles. In Fig. 29b, c impacts characterised by u 0 = 2 and 8 m/ s, respectively, are shown. In both cases is evident that the maximum impact force (as was also observed by Calvetti et al. [12]), the first buckling and the residual force are not affected by the velocity profile. In contrast, secondary oscillations tend to be less marked in case the velocity decreases with depth. Nevertheless, the overall impulse value, that is the most important design factor, is negligibly affected by the shape of velocity profile.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement.
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 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.